# When Conventional machine learning meets neuromorphic engineering: Deep
Temporal Networks (DTNets) a machine learning frawmework allowing to operate
on Events and Frames and implantable on Tensor Flow Like Hardware

Marco Macanovic2, Fabian Chersi2, Felix Rutard2, ,
, Sio-Hoi Ieng 2, Ryad Benosman123
———————————————————
1 University of Pittsburgh/Medical Center, 2 Sorbonne Universitas,
,
3 Carnegie Mellon University

###### Abstract

We introduce in this paper the principle of Deep Temporal Networks that allow
to add time to convolutional networks by allowing deep integration principles
not only using spatial information but also increasingly large temporal
window. The concept can be used for conventional image inputs but also event
based data. Although inspired by the architecture of brain that inegrates
information over increasingly larger spatial but also temporal scales it can
operate on conventional hardware using existing architectures. We introduce
preliminary results to show the efficiency of the method. More in-depth
results and analysis will be reported soon!

## Introduction

Temporal surfaces introduced in [1] allowed to learn temporal sequences of
events on different temporal and spatial scales called HOTS (for Hierarchy of
Temporal Surfaces). This learning method allows classification from the ouput
of event based sensor, but can also be extended to any signal including
conventional video sequences. Although promissing Neurmorphic event based
engineering is still lacking a convincing computational paltform that allows
to make full use of its potiential and mostly to take advantage of all the
properties of event based cameras. HOTS is also very much linked to event
based data and is not able to operate on frames and more important make use of
existing GPU hardware. Event-based cameras are offering many advantages that
frame-based cameras are not able to provide without an unreasonable increase
in computational resources. Low computation needs is especially achieved by a
lower redundant data acquisition, while presenting lower latency than standard
cameras, achieved via a highly precise temporal and asynchronous level
crossing sampling. The main idea behind HOTS is to introduce a temporal
context around incoming events to create a temporal context that we called a
time surface. This concept can be extended to multiple time scale binding and
integrating information on also larger sptial scales.

The first layers of the architecture being fed with the output of the camera
need to operate at time scales of a around hundred of microsecond to few
millisecond that is incompatible with existing hardware.In fact altough
processors can operate fast, random memory access as needed by sparse data
output from event based camera is a major limitation specially when one want
to operate fast at the native resolution of the sensor while being able to
absorb event rates of the order of Giga Events per seconds. Deeper scales
however, are integrating information on larger times scales and do not need to
operate fast. This is very similar to biological brains where low level areas
need to process at the submillisecond the output of senses while higher areas
and prefrontal cortex operate on slower time scales. We can then consider
temporal integration steps above the order or 15-20 ms to be easily and
efficiently performed by conventional architectures. Specially that when one
integrates cues on larger time scales the amout of information tend to
stabilize and decrease, the last layers are naturally compressing information.

AI has drastically improved during the last decade with the advent of more
adapted hardware such as GPU with also the advent of increasingly more
optimized hardware. As we will show it is possible to take the best out of the
conventional narrow AI and the neuromorphic world. The solution lays in
rethinking temporal and spatial integration to be able to operate on GPU like
hardware while maintaining its concept of temporal and spatial scale
integration. We coined this Methodology Deep Temporal Nets to be more in
adequation with the current trend of deep AI technolgy to which it can add a
better use of time. DTNets rely on a new hierachical architecture relying on
the use of autoencoders to compress Time surfaces. We can then also add
pooling possiblities to full comply with existing deep structures.

The DTNets machine learning architecture is the same line of thought as deep
neural network. It is based on two main computational features that constitute
the complete architecture: the autoencoders for features extraction and the
classifiers to determine classes. DTNets take as input any type of
asynchronous signals based on existing techniques to sparsify signals such as
sigma-delta, level crossing or relative change. But is can also use
conventional video sequences not requiring to encode graylevel values in time
by generating ’fake’ events. DTNets are hierarchical multilayer structures
composed of two computational blocks. The first block extracts features using
autoencoders. Autoencoders allow to extract features and abstract them into
more and more complex representation. The second block is a classification
layer that takes as input the output from the last layer and feeds it into a
classifier.
A wide variety of classifiers can be used, from Multi Layer Perceptron
currenlty in use, SVM, to more advanced techniques such as LSTM. Currently the
retained classifier is MLP based and uses backpropagation to train. It
important to emphasize that both layers rely on backpropagation for feature
extraction and classification.

##  1 Time Surfaces

Figure 1: Definition of a time-surface from the spatio-temporal cloud of
events from [1]. A time-surface describes the recent time history of events in
the spatial neighborhood of an event. This figure shows how the time-surface
for an event happening at pixel $\mathbf{x_{0}}=[x_{0},y_{0}]^{T}$ at time
$t_{0}$ is computed. The event-driven time-based vision sensor (a) is filming
a scene and outputs events shown in (b) where ON events are represented on the
left hand picture and OFF events on the right hand one. For clarity, we
continue by only showing values associated to OFF events. When an OFF event
$ev_{i}=[\mathbf{x_{0}},t_{i},-1]$ arrives, we consider the times of most
recent OFF events in the spatial neighborhood (c) where brighter pixels
represent more recent events. Extracting a spatial receptive field allows to
build the event-context $\mathcal{T}_{i}(\mathbf{x},p)$ (d) associated with
that event. Exponential decay kernels are then applyed to the obtained values
(e) and their values at $t_{i}$ constitute the time-surface itself. (f) shows
these values as a surface. This representation will be used in the following
figures and the label of the axes will be removed for better clarity.

The process of building a time-surface from the output of an event-driven
time-based vision sensor as introduced in [1] is illustrated in Fig. 1 and
described hereafter.

Consider a stream of visual events (Fig. 1(b)) which can be mathematically
defined as

$\begin{array}[]{l
r}ev_{i}=[\mathbf{x_{i}},t_{i},p_{i}]^{T},&i\in\mathbb{N}\\\ \end{array}$ (1)

where $ev_{i}$ is the $i^{th}$ event and consists of a location
($\mathbf{x_{i}}=[x_{i},y_{i}]^{T}$), time ($t_{i}$) and polarity ($p_{i}$),
with $p_{i}\in\\{-1,1\\}$, where $-1$ and $1$ represent OFF and ON events
respectively. When an object (or the camera) moves, the pixels asynchronously
generate events which form a spatio-temporal point cloud representing the
object’s spatial distribution and dynamical behavior. Fig. 1(b) shows such
events generated by an object rotating in front of the sensor ((Fig. 1(a))
where ON and OFF events are represented respectively by white and black dots.

Because the structure of this point cloud contains information about the
object and its movement, we introduce the time-surface $\mathcal{S}_{i}$ of
the $i^{th}$ event $ev_{i}$ to keep track of the activity surrounding the
spatial location $\mathbf{x_{i}}$ just before time $t_{i}$. We can then define
$\mathcal{T}_{i}(\mathbf{u},p)$ a time-context around an incoming event
$ev_{i}$ as the array of most recent events times at $t_{i}$ for the pixels in
the $(2R+1)\times(2R+1)$ square neighborhood centered at
$\mathbf{x_{i}}=[x_{i},y_{i}]^{T}$ as:

$\mathcal{T}_{i}(\mathbf{u},p)=\max_{j\leq
i}\left\\{t_{j}\,|\,\mathbf{x_{j}}=(\mathbf{x_{i}}+\mathbf{u}),\,p_{j}=p\right\\},$
(2)

where $\mathbf{u}=[u_{x},u_{y}]^{T}$ is such that $u_{x}\in\\{-R,\dots,R\\},$
$u_{y}\in\\{-R,\dots,R\\}$ and $p\in\\{-1,1\\}$
$\mathcal{T}_{i}(\mathbf{x},p)$ is shown in Fig. 1(d) where intensity is
coding for time values: bright pixels show recent activity whereas dark pixels
received events further away in the past (only time values corresponding to
OFF events are represented in the figure for clarity).

Let $\mathcal{S}_{i}(\mathbf{u},p)$ be the time-surface around an event
$ev_{i}$ (shown in Fig. 1(e)), it is defined by applying an exponential decay
kernel with time constant $\tau$ on the values of
$\mathcal{T}_{i}(\mathbf{u},p)$.

$\mathcal{S}_{i}(\mathbf{u},p)=e^{-(t_{i}-\mathcal{T}_{i}(\mathbf{u},p))/\tau}.$
(3)

$\mathcal{S}_{i}$ provides a dynamic spatiotemporal context around an event,
the exponential decay expands the activity of passed events and provides
information about the history of the activity in the neighborhood. The
resulting surface $\mathcal{S}_{i}(\mathbf{u},p)$ is shown in Fig. 1(f) for
the OFF events represented all along Fig. 1. In the following sections
$\mathcal{S}_{i}(\mathbf{u},p)$ will be referred to directly as
$\mathcal{S}_{i}$ to simplify notations. In the figures it will be represented
as a surface showing the values of each of its element at their corresponding
spatial positions.

## Hierarchy of Temporal and Spatial Integration

###  1.1 Creating a Hierarchical Model

Figure 2: View of the proposed hierarchical model as introduced in [1]. From
left to right, a moving digit (a) is presented to the event based camera (b)
which produces ON and OFF events (c) which are fed into Layer $1$. The events
are convolved with exponential kernels (d) to build event contexts from
spatial receptive field of sidelength $(2R_{1}+1)$. These contexts are
clustered into $N_{1}$ features (e). When a feature is matched, it produces an
event (f). Events from the $N_{1}$ features are merged into the output of the
layer (g). Each layer $k$ (gray boxes) takes input from its previous layer and
feeds the next by reproducing steps (d)-(g). The output of Layer $k$ is
presented between Layer $k$ and $k+1$ ((g),(h),(i)). To compute event
contexts, each layer considers a receptive field of sidelength $(2R_{k}+1)$
around each pixel. The event contexts are compared to the different features
(represented as surfaces in the gray boxes) and the closest one is matched.
The images next to each features show the activation of their associated
features in each layer. These activations are merged to obtain the output of
the layer. The output (i) of the last layer is then fed to the classifier (j)
which will recognize the object.

Fig. 2 illustrates the hierarchical model we introduce in this paper. Steps
(a) to (g) sum up the process described in the previous sections. As shown in
Fig. 2, a moving digit (a) is presented to the event based camera (b) which
produces ON and OFF events (c). Time-surfaces are built by convolving them
with an exponential kernel of time constant $\tau_{1}$ (d) and considering
spatial receptive fields of sidelength $(2R_{1}+1)$. These time-surfaces are
then clustered into $N_{1}$ prototypes represented as surfaces (e) in the
Layer $1$ box. When a cluster center is matched, an event is produced,
resulting in the activations shown in (f). These events are merged to form the
output of Layer $1$ (g). One can see that each incoming event from the
observed pattern is associated with the most representative prototype surface.
The nature of the output of Layer $1$ is exactly the same as its input: Layer
$1$ outputs timed events. Once a prototype matches the temporal surface around
the incoming event it immediately emits an event. Thus, the same steps used in
Layer 1 (from (d) to (g)) can be applied in Layer $2$. However the emitted
event is now representing the temporal activity of a prototype surface, it
thus carries more meaning than the initial camera event. The prototype
surfaces of Layer $2$ represent the temporal signature of the activity of
complex features. Layer 2 uses different constants for space-time integration
of features ($R_{2}$, $N_{2}$ and $\tau_{2}$). The goal is to introduce
stability of the perceptual representation and sensitivity to the accumulation
of sensory evidence over time. This integration over longer and longer time
period will thus be able to accumulate evidence in favor of alternative
propositions in a recognition process. When alternatives with a barely
discernible difference in their sensory inputs are presented over an extended
period of time, longer time and spatial integration scales can accumulate the
small differences over time until it becomes eventually possible to
discriminate the alternatives through its ever growing output. This
accumulation dynamics is at the heart of the HOTS model, the difference
between time scales can be substantial and can start from 50ms for Layer $1$
to 250ms for Layer $2$ to finally reach 1.25 s for Layer $3$.
Layer $3$ receives input from Layer $2$, it is the last layer of the system
and it provides the highest level information integration, as shown in Fig.
2(i) time-surface prototypes are also larger both spatially and temporally.
The output of the temporal activity of Layer $3$ can finally be used for
object recognition by being fed to a classifier (shown in Fig. 2(j)).

As stated above, each layer is then defined by only a few parameters (we add
an index $l$ for the $l$th layer of the system):

  * •

$R_{l},$ which defines the size of the time-surface neighborhood

  * •

$\tau_{l},$ the time constant of the exponential kernel applied to events

  * •

$N_{l},$ the number of cluster centers (prototypes) learnt by the clustering
algorithm.

To increase the information extracted by each subsequent layer, we make these
parameters evolve between subsequent layer. For each layer, we define the
parameters $K_{R}$, $K_{\tau}$, $K_{N}$ so that:

$\displaystyle R_{l+1}$ $\displaystyle=$ $\displaystyle K_{R}\cdot R_{l}$ (4)
$\displaystyle\tau_{l+1}$ $\displaystyle=$ $\displaystyle
K_{\tau}\cdot\tau_{l}$ (5) $\displaystyle N_{l+1}$ $\displaystyle=$
$\displaystyle K_{N}\cdot N_{l}$ (6)

The obtained architecture consists in a Hierarchy Of Time-Surfaces (HOTS)
which is building and extracting a set of features (the prototypes from the
final layer) out of a stream of input events. The time-surface prototypes will
then be called time-surface features in the rest of the paper.

###  1.2 Replacing Clustering by Auto-Encoders, compressive coding and more

Clustering has been used to determine from a set of computed temporal contexts
which ones are the most used, namely representing dominant information of
scenes. Several other techniques could be used according to what type of
information one would want to extract. In [2] we used a sparse coding approch
such as introduced in [3]. One could consider many other strategies to extract
valuable time surfaces. Here we use autoencoders as compressive mean of
selecting time surfaces. This has the advantage of providing a compact
representation of features and more important a higher generalization of
features. An autoencoder always consists of two parts, the encoder and the
decoder (see Fig.3), which can be defined as transitions $\phi$,$\psi$ such
that:

$\phi:\mathcal{X}\rightarrow\mathcal{F}$
$\psi:\mathcal{F}\rightarrow\mathcal{X}$
$\phi,\psi=\arg\min_{\phi,\psi}\|\hat{\mathcal{S}}_{i}(\mathbf{u},p)-(\psi\circ\phi)\mathcal{S}_{i}(\mathbf{u},p)\|^{2}$

We consider a single hidden layer, the encoder stage of the autoencoder takes
the input volume patch $\mathcal{S}_{i}(\mathbf{u},p),$ and maps it to
$\mathbf{z}\in\mathbb{R}^{p}=\mathcal{F}$ :

$\mathbf{z}=\sigma(\mathbf{W\mathcal{S}_{i}(\mathbf{u},p)}+\mathbf{b})$

$\mathbf{z}$ is usually referred to as code, latent variables, or latent
representation. Here, $\sigma$ is an element-wise activation function such as
a sigmoid function or a rectified linear unit.

$\mathbf{W}$ is a weight matrix and $\mathbf{b}$ is a bias vector. After that,
the decoder stage of the autoencoder maps $\mathbf{z}$ to the reconstruction
$\hat{\mathcal{S}}_{i}(\mathbf{u},p)$ of the same shape as :

$\hat{\mathcal{S}}_{i}(\mathbf{u},p)=\sigma^{\prime}(\mathbf{W^{\prime}z}+\mathbf{b^{\prime}})$

where $\mathbf{\sigma^{\prime}},\mathbf{W^{\prime}},\text{ and
}\mathbf{b^{\prime}}$ for the decoder may differ in general from the
corresponding $\mathbf{\sigma},\mathbf{W},{\text{ and }}\mathbf{b}$ for the
encoder, depending on the design of the autoencoder.

The hidden layer $\mathbf{z}$ encode a compressed version of time surfaces and
can be used as an output feature to the next layer. Several possiblities can
considered. The content can be sent as it is to the next layer. One can also
consider as introduced in [2] to encode values of $\mathbf{z}$ as time delay,
as shown in Figure.3 by considering the new time as the time of arrival of the
time surface $t_{in}$ plus an additional delay with regards to the value of
value of $\mathbf{z}_{i}$ weighted by a scalar $\alpha$. One can also consider
more crude version by simply thresholding the values of value of
$\mathbf{z}_{i}$ to operate in the next layer on the meaninglfull components
of $\mathbf{z}$. if raw values are used one could use a pooling layer
similarly using a wide variety of pooling strategies similarly to what is
being done in conventional deep convolution networks

Figure 3: Schematics of the single layer autoencoder used the select features
from incoming time surfaces.

The general architecture of the system is then as follows: as show in Figure.

##  2 General Architecture

Figure 4: General architecture of DTNets.

A moving pattern observed by the event based camera (Figure.4(a)). For each
incoming event at a spatial location $\mathbf{p}$, we define a neighbourhood
as shown in yellow in Figure.4(c) for which we which we compute a time surface
$\mathcal{S}^{1}(\mathbf{u},p)$ centered on the event (the subscript 1 stands
for the layer 1). The time surface is fed into the autocoder that will
compress the time surface and reduce dimensionality (4(c)). Training is
performed layer per layer, therefore at an initial stage, the autoencoder is
trained for each incoming time surface, one the reconstruction erros is lower
that a certain percentage, the training is considered complete, and the output
of the encoder is sent to the rest of the processing chain. The autoencoder
hidder layer $\mathbf{z}^{1}$ contains the compressed information and its size
$N_{z}$ is smaller than the size of the input time surface $(2R_{1}+1)^{2}$.
The output if the encoder is stored in a feature space (or features volume) at
the location $\mathbf{p}$. The size of the volume is equal to the size of the
focal plane $N_{l}*N_{c}$ (minus $2*R_{1}$ because of borders, but for clarity
this precision is omitted in the figure), and depth $N_{z}$, the length of the
hidden layer $z^{1}$. We can then define a volume in the feature space of size
$R_{2}*R_{2}*N_{z}$ as show in purple in Figure.4(f). The content of the
volume can processed in several ways. Many options are possible, the content
of the volume can be kept as it is and pooled, or it can be encocded in the
time domain on a larger time scale $\tau_{2}$, then clustered. All these
strategies have both advantages and disadvantages that will be discussed
further. The output of these two stage shown in Figure.4(g-h) is then sent to
another time encoder that will compress the information. The content of the
hidder layer $\mathbf{z}^{2}$ can then be sent to another volume and so on.
Finally, the last layer sends its ouput to a classifier.

##  3 Experiments

we will first evluate a network’s configuration that uses the direct ouput of
the autoencoders wihtout adding any temporal concept. The network is simply
encoding time surfaces at the inpout layer and then similarly to what
conventional convolution networks operate. The idea of the network in to
create feature exploiting the combinatory of the input without adding more
temporal integration scales.The network is basically working of learning
sequences without adding more temporal information into the network.

The perfomances are evaluated on two learning tasks: 1) the recognition of
asynchronous handwritten characters from the N-MNIST database, used by the
neuromorphic community. 2) the recognition of the presence or absence of a car
in an N-CARS database. Once the auto-encoders are running with the best
reconstruction rate, we go out of the last layer to realize the final overall
volume of the example and then vectorize it to finally send it to a
classifier, here a Multi Layer Perceptron (MLP).

We can see in Figure.5 and Figure.6 that the overall recognition rate is 93.8%
based on N-MINIST data. The parameters used are:

  * •

Single layer of auto-encoder

  * •

A context radius $R_{1}$= 2

  * •

A time constant $\tau_{2}$ = 30ms

  * •

The number of neurons in the hidden layer of the auto-encoder N1 = 10 - The
number of neurons in the hidden layer of the MLP NMLP = 200

Figure 5: Confusion matrix on the NMNIST Database Figure 6: Recognition rates
on N-MNIST Figure 7: Recognition rates on N-CARS

The N-CARS database is more complex to classify, we are introducing
preliminary work, further tuning of the parameters will surely improve this
rate (see Figure.6). The parameters used are:

  * •

Single layer of auto-encoder

  * •

A context radius $R_{1}$ = 3

  * •

A time constant $\tau_{1}$1 = 40ms

  * •

The number of neurons in the hidden layer of the auto-encoder N1 = 20

  * •

The number of neurons in the hidden layer of MLP NMLP = 200 The use of several
layers greatly increases the dimensionality

##  4 Conclusions

These results although preliminary shows that the algorithm performs well
altough using only one layer and not integrating any additional temporal
scales that is the essence of these temporal networks. We will report more in
depth results and more advanced tuning strategies as soon as as our current
non optimized implementation is improved as it currently tends to be slow to
execute. Compared with other results of methods of classification of event
based machine learning these Deep Temporal Networks provide promissing results
and more imporant they can be easily implemented on GPU while maintaining an
event based input but also any type of frame based input too. The method
scales easily and can be added a variety of improvements to compress the data
from layer to layer. We beleive that type type of structure not only
integrates over space but can also integrzte easily over time while benefiting
form the conventional machinel learning advances.

## Aknowledgments

We are gratefull to the fondation Voir et Entendre and also Streetlab to have
funded this project. We also would like to thanks all member of the laboraotry
not mentioned in this paper that also contributed to develop temporal
networks.

## References

  * [1] X. Lagorce, G. Orchard, F. Galluppi, B. E. Shi, and R. B. Benosman, “Hots: A hierarchy of event-based time-surfaces for pattern recognition,” IEEE Transactions on Pattern Analysis Machine Intelligence, vol. 39, pp. 1346–1359, July 2017.
  * [2] G. Haessig and R. Benosman, “A sparse coding multi-scale precise-timing machine learning algorithm for neuromorphic event-based sensors,” CoRR, vol. abs/1804.09236, 2018.
  * [3] B. Olshausen and D. Field, “Emergence of simple-cell receptive field properties by learning a sparse code for natural images,” Nature, vol. 381, pp. 607–609, 1996.


# Unbalanced Job Approximation using Taylor Series Expansion and Review of
Performance Bounds

Alexander Thomasian
Thomasian & Associates
Pleasantville  NY 10570 USA
<EMAIL_ADDRESS>
###### Abstract

Unbalanced Job Approximation - UJA is a family of low-cost formulas to obtain
the throughput of Queueing Networks - QNs with fixed rate servers using Taylor
series expansion of job loadings with respect to the mean loading. UJA with
one term yields the same throughput as optimistic Balanced Job Bound - BJB,
which at some point exceeds the maximum asymptotic throughput. The accuracy of
the estimated throughput increases with more terms in the Taylor series. UJA
can be used in parametric studies by reducing the cost of solving large QNs by
aggregating stations into a single Flow-Equivalent Service Center - FESCs
defined by its throughput characteristic. While UJA has been extended to two
classes it may be applied to more classes by job class aggregation. BJB has
been extended to QNs with delay servers and multiple jobs classes by Eager and
Sevcik, throughput bounds by Eager and Sevcik, Kriz, Proportional Bound - PB
and Prop. Approximation Bound - PAM by Hsieh and Lam and Geometric Bound - GB
by Casale et al. are reviewed.

##  1 Introduction to Queueing Network Modeling

Queueing Networks - QNs are used to model the realistic situation where
customers (resp. jobs) require processing at more than one service station
(resp. device at computer system). QNs have been used in many application
domains such as manufacturing systems, computer communication networks, and
multiprogrammed computer systems.

While efficient solution methods have been developed for a subcategory of QN
models known as product-form, the solution cost increases with the number of
stations, number of job classes and the number of jobs.

QN models can be categorized to open, closed, and mixed. In open QNs there are
external job arrivals and completed jobs leave the system, while in closed QNs
completed jobs are immediately replaced with new jobs, as if there is an
infinite backlog of jobs. We consider product-form QNs, which lend themselves
to an exact solution and can be analyzed relatively efficiently Kleinrock 1975
[24].

The analysis of open product-form QNs with Poisson arrivals is trivial in that
each station can be analyzed separately Jackson 1957 [22]. Closed PF QNs
require the degree of concurrency and similarly to open QNs service demands or
loadings, which is the mean total time jobs spend being served at service
stations.

Unbalanced Job Approximation - UJA Thomasian and Nadji 1981 [38] is a family
of approximations based on Taylor series expansions of job loading with
respect to the average loading $X_{0}$. UJA starts with an upper bound on
throughput, which equals the optimistic throughput of Balanced Job Bounds -
BJBs Zahorjan et al. 1982 [43]. UJA obtains more accurate throughput estimates
by using a Taylor series expansion with respect to $X_{0}$. UJA was presented
with BJB at SIGMETRICS’81 Zahorjan et al. 1981 [42], but unlike BJB which was
publicized in Lazowska et al. 1984 [28], UJA is little known and hence this
publication.

UJA was applied in the analysis of the QN model of a computer system running
CAD applications originated by users of Time-Sharing Option - TSO on an IBM
mainframe with the MVS operating system. A performance study based on BEST/1
capacity planning tool [5] is reported by Silvester and Thomasian 1981 [36].
Given large number of disks measurement results were used to aggregate disks
with almost balanced utilizations, while treating disks with very low
utilizations as delay servers [28].

The paper is organized as follows. Product-form QNs are described in Section
2. The GF method for analyzing closed QNs is presented in Section 3.
Aggregation of balanced stations to a Flow Equivalent Service Center - FESC is
described in Section 3.1. This is followed UJA for a single and two job
classes in Section 4 and Section 5 based on Nadji and Thomasian 1981/1984 [38,
31] Conclusions are drawn in Section 6. Appendix I presents Performance Bound
Hierarchies - PBH by Eager and Sevcik 1983,86 [13, 14] Appendix II presents
extensions to BJB by Kriz 1984 [26]. Appendix III discusses asymptotic
expansions with multiple bottlenecks by George et al. 2012 [15]. Appendix IV
presents Geometric Bounds - GB by Casale et al. [CaSM08] and compares them
with Proportional Bounds - PB by Hsieh and Lam [19]. Appendix V present
Proportional Approximation Method - PAM for multichain QNs by Hsieh and Lam
19989 [21]. We have preserved the notation of the original papers in the
Appendices to make it easer for readers to refer to them.

##  2 Closed QNs with Product-Form Solutions

Product-form QNs allow stations with single and multiple servers with
exponential service times and FCFS scheduling Jackson 1957 [22]. Arrivals to
open QNs are Poisson and jobs are routed from station to station according to
probabilistic routing. The analysis of such QN models is inexpensive, since
each station can be analyzed separately as if it is subjected to Poisson
arrivals.

The arrival rates of jobs to the $N$ stations is given by:

$\underline{\gamma}=(\gamma_{1},\gamma_{2},\ldots,\gamma_{N})$

Job arrival rates to the $N$ stations, denoted by $\lambda_{n},1\leq n\leq N$,
are obtained by solving the set of linear equations:

$\displaystyle\underline{\lambda}=\underline{\gamma}+\underline{\lambda}P\mbox{
hence }\underline{\lambda}=[I-P]^{-1}\underline{\gamma}.$ (1)

Server utilization factor at the $n^{th}$ station is
$\rho_{n}=\lambda_{n}\bar{x}_{n}/m_{n}$, where $\bar{x}_{n}$ is the mean
service time and $m_{n}$ the number of servers. In the case $m_{n}=1$ and
exponential service times $\bar{x}_{n}=1/\mu_{n}$, the mean residence time at
the station is that of an M/M/1 queue.

$R_{n}=\bar{x}_{n}/(1-\rho_{n})=(\mu_{n}-\lambda_{n})^{-1}\mbox{ for
}\rho_{n}<1$

The formula for M/M/m queues is given in [24].

#### Central Server Model - CSM

CSM is a closed QN model of a multiprogrammed computer system proposed by
Buzen 1973 [4], which consists of a CPU and multiple disks. Jobs alternate
between CPU and disk processing until they are completed. An efficient
solution method known as the Convolution Algorithm - CA is presented in [4] to
analyze CSM and closed product-form QNs in general.

The CPU is designated as the central station ${\cal S}_{1}$ and the $N-1$
disks as peripheral stations ${\cal S}_{n},2\leq n\leq N$. Given the state
transition probabilities ${\cal S}_{i}\xrightarrow{p_{i,j}}{\cal S}_{j}$ the
following transitions are applicable to CSM.

$p_{1,n},\hskip 5.69054pt2\leq n\leq N,\hskip 8.53581ptp_{n,1}=1,\hskip
5.69054pt2\leq n\leq N$

The self-transition $p_{1,1}=1-\sum_{n=2}^{N}p_{1,n}$ implies the completion
of a job in a closed QN (or a job that leaves the system in an open QN). The
number of visits to the CPU is given by the geometric distribution [24].

$q_{k}=p_{1,1}(1-p_{1,1})^{k-1},k\geq 1\hskip
5.69054pt=\bar{k}=v_{1}=1/p_{1,1}.$

The relative number of visits to the stations is obtained by solving

$\underline{v}=\underline{v}{\bf P}$

It follows $v_{n}=p_{1,n}v_{1}=p_{1,n}/p(1,1),2\leq n\leq N$.

Given mean service time at ${\cal S}_{n}$ per visit is $\bar{x}_{n}$, the mean
loading per job is $X_{n}=v_{n}\bar{x}_{n},1\leq n\leq N$.

Inputs required for the analysis of closed QNs are the number of jobs or the
degree of concurrency, probabilistic routings among stations which lead to the
mean number of visits to the stations and mean service times per visit, which
yields loadings. Performance metrics of interest in modeling closed QNs are
system throughput, mean residence time, device utilizations and mean queue-
lengths.

Jackson theorem which only allowed service stations with exponential servers
with FCFS scheduling was extended by the BCMP theorem to four types of
stations Baskett et al. 1975 [3].

$\displaystyle F_{n}(k)=\begin{cases}X_{n}^{k}\hskip 28.45274pt\mbox{FCFS, PS,
LCFSPR}\\\ X_{n}^{k}/k\hskip 25.60747pt\mbox{Infinite-Server - IS}\\\
X_{n}^{k}/\prod_{j=1}^{k}a(j)\hskip 17.07164pt\mbox{Queue-
dependent}\end{cases}$ (2)

In a queue-dependent station the service rate varies with the number of jobs:

$F_{n}(k)=X_{n}^{k}/a(k),\mbox{ where }a(k)=\mu(k)/\mu(1)$

$a(k)$ is the ratio of the service rate with $k$ versus one job at the
station. Multiserver and IS or delay server queues are special cases with
$a(k)=k,1\leq k\leq m$ for $m$-server. For IS or delay servers $m=\infty$.
General service times are allowed for Processor Sharing - PS and IS service
stations.

Buzen 1973 [4] developed the Convolution - CA algorithm to analyze closed QNs
and applied it to CSM. We consider mostly fixed rate single server and delay
servers in our discussion.

The solution of open product-form QN models is trivial. Given the arrival rate
to each station its mean residence time can be obtained by solving the
corresponding station independently.

Consider the processing of $K$ jobs in a closed QN with a single jobs with $N$
stations. Job service demands or loadings are the products of mean service
time $\bar{x}_{n}$ and the mean number of visits ($v_{n}$) to station ${\cal
S}_{n}$, i.e., $X_{n}=v_{n}\bar{x_{n}}$.

The steady-state state probability in a closed QN with $K$ jobs with $N$
fixed-rate single server-stations with job distribution:

$\underline{k}=(k_{1},k_{2},\ldots,k_{N})\mbox{ with }\sum_{n=1}^{N}=K.$

$\displaystyle P[\underline{k}]=\frac{X_{1}^{k_{1}}X_{2}^{k_{2}}\ldots
X_{N}^{k_{N}}}{G(K)}\hskip 5.69054ptG(K)=\sum_{\underline{k}\in{\cal
K}}\prod_{n=1}^{N}X_{n}^{k_{n}}.$ (3)

The number of terms grows rapidly with the size of the network such that a
direct summation of multinomial expressions may be computationally
intractable. Multiple job classes or chains yield a more realistic
representations of QN models, where each class has its own routings
probabilities and service times per visit and hence different loadings at the
stations. Given $R$ job classes and $k_{r},1\leq r\leq R$ jobs in class $r$,
the number of states is.

$\prod_{r=1}^{R}\binom{K_{r}+N-1}{N-1}$

There are several efficient computational methods to analyze product-form
closed QNs. The Convolution Algorithm - CA was proposed in the context of CSM
with a single job class Buzen 1973 [4] and was extended as part of the BEST/1
capacity planning package by Buzen et al. [5].

IBM’s Reiser and Kobayashi 1975 [32] used generating functions, see e.g.,
[24], in extending BCA to multiple job classes. The application of GFs to
analyze closed QNs is presented in a tutorial manner in Williams and Bhandiwad
1976 [41], Thomasian and Nadji 1981 [38], and Trivedi 2002 [40],

##  3 Generating Functions for Analyzing Product-Form QNs with a Single Job
Class

The GF for $n^{th}$ QN station for $1\leq n\leq N$ is given by:

$\displaystyle x_{n}(t)=\sum_{k=0}^{\infty}F_{n}(k)t^{k},$ (4)

where $F_{n}(k)$ was given by Eq. (2).

For stations with fixed-rate single servers $F_{n}(k)=X_{n}^{k}$.

$\displaystyle x_{n}(t)$ $\displaystyle=\sum_{k\geq
0}(X_{n}t)^{k}=1+(X_{n}t)+(X_{n}t)^{2}+\ldots$ (5)
$\displaystyle=(1-X_{n}t)^{-1},|X_{n}t|<1.$

If the $n^{th}$ station is IS:

$\displaystyle x_{n}(t)=1+X_{n}t+\frac{(X_{n}t)^{2}}{2!}+\dots=e^{X_{n}t}.$
(6)

The GF for the QN is the product of the GFs for individual stations.

$\displaystyle g(t)=\prod_{n=1}^{N}x_{n}(t).$ (7)

which can be rewritten as:

$\displaystyle g(t)=1+G(1)t+G(2)t^{2}+\ldots+G(K)t^{K}+\ldots,$ (8)

In the case of single-server stations with $N=2$ stations as $K$ increases:

$\begin{cases}G(1)=X_{1}+X_{2},\\\ G(2)=X_{1}^{2}+X_{2}^{2}+X_{1}X_{2},\\\
G(3)=X_{1}^{3}+X_{1}^{2}X_{2}+X_{1}X_{2}^{2}+X_{2}^{3},\\\ \vdots\end{cases}$

Given

$g_{n}(t)=g_{n-1}(t)x_{n}(t),1\leq n\leq N$

and noting that $x_{n}(t)=(1-X_{n}t)^{-1}$ it follows

$g_{n}(t)=g_{n-1}(t)+X_{n}g_{n}(t),\hskip 14.22636pt1\leq k\leq K,1\leq n\leq
N$

Equating the coefficients of $t^{k}$ on both sides we have:

$\displaystyle\boxed{G_{n}(k)=G_{n-1}(k)+X_{n}G_{n}(k-1),}$ (9) $\displaystyle
1\leq k\leq K,1\leq n\leq N.$

The utilization of a single server at station $n$ is the sum of state
probabilities with at least one job at the station:

$\displaystyle\boxed{U_{n}(K)=X_{n}T(K),}$ (10)

The system throughput follows from Little’s formula [24]:

$\displaystyle\boxed{T(K)=\frac{U_{n}(k)}{X_{k}}=\frac{G(K-1)}{G(K)}.}$ (11)

If $X$s are expressed in seconds $T(K)$ with be jobs/second.

The mean queue length at the $n^{th}$ single-server station is:

$\displaystyle\boxed{Q_{n}(K)=\frac{1}{G(K)}\sum_{k=1}^{K}X_{n}^{k}G(K-k).}$
(12)

CA to compute the normalization constant $G(K)$ in a single class QN is:

Initializations:

$\begin{cases}G_{0}(k)=0\hskip 14.22636pt1\leq k\leq K.\\\ G_{n}(0)=1\hskip
14.22636pt0\leq n\leq N\end{cases}$

for $k=1$ to $K$ do
for $n=1$ to $N$ do
$G_{n}(k)=G_{n-1}(k)+X_{n}G_{n}(k-1)$
end
end

###  3.1 Aggregation of Balanced Stations

Given that a subset $M\leq N$ of the QN stations have balanced loads, i.e.,
$X_{n}=X,1\leq n\leq M$, it follows from Eq. (7) that:

$\displaystyle v(t)=[1+Xt+(Xt)^{2}+\dots]^{M}=(1-Xt)^{-M}$ (13)

Applying the Binomial Theorem we have:

$\displaystyle v(t)=\sum_{k=0}^{\infty}\binom{M+k-1}{k}(Xt)^{k}.$ (14)

The GF for the aggregate station can be written as:

$\displaystyle v(t)=\sum_{k=0}^{\infty}(MXt)^{k}/\prod_{j=1}^{k}a(j),\hskip
5.69054pta(j)=\frac{jM}{M+j-1}.$ (15)

The throughput of the using Eq. (11) is:

$\displaystyle T(k)=\frac{\binom{M+k-2}{k-1}X^{k-1)}}{\binom{M+k-1}{k}X^{k}}$
(16)

It follows that job throughout with $j$ jobs is:

$\displaystyle\boxed{T(k)=\frac{k}{M+k-1}\frac{1}{X}}$ (17)

As $k\rightarrow\infty$ the maximum throughput is $T_{max}=1/X$

We can use the symmetry due to balancedness to obtain the mean residence time
at single-server stations using the key equation in Mean Value Analysis - MVA
Reiser and Lavenberg 1980 [33] which is based on the Arrival Theorem Lavenberg
and Reiser 1980 [27].

In a QN with $k$ jobs a job arriving at a station sees a mean queue-length
with one job less. The $k-1$ jobs are equally distributed among $M$ stations,
hence:

$r_{n}(k)=\bar{x}_{k}\left(1+\frac{k-1}{M}\right)$

Multiplying both sides by $v_{k}$:

$\displaystyle\boxed{R_{n}(k)=X_{n}\left(1+\frac{k-1}{M}\right).}$ (18)

The mean residence time in the system is obtained by multiplying by $M$.

$R(k)=Mr(k)=(M+k-1)X,$

from which the throughput $T(k)=k/R(K)$ as given by Eq. (16) follows.

The aggregation of IS or delay station in a QN is particularly simple, so that
all such stations whose indices are in the set ${\cal I}$ can be replaced by a
single IS station with service demand

$\boxed{X_{I}=\sum_{i\in{\cal I}}X_{i}.}$

##  4 Aggregation of Unbalanced Stations in a Single Class

We use the throughput of an FESC with balanced stations as a starting point to
obtain the throughout of an FESC with unbalanced stations. Defining the mean
service demand over $M$ fixed-rate single-servers stations:

$X_{0}=\sum_{n=1}^{M}X_{n}/M.$

The deviation of service demands from the mean and its moments are as follows:

$\displaystyle e_{n}=\frac{X_{n}-X_{0}}{X_{0}},\hskip 5.69054pt1\leq n\leq
M\mbox{ and }E_{j}=\sum_{n=1}^{M}e_{n}^{j}.$ (19)

The GF of the $n^{th}$ station can be given as its Taylor series expansion as
follows:

$\displaystyle x_{n}(t)=x_{0}(t)+(X_{n}-X_{0})x_{0}^{(1)}(t)+$ (20)
$\displaystyle\frac{1}{2!}(X_{n}-X_{0})^{2}x_{0}^{(2)}(t)+\ldots,$

Making the substitution $(X_{n}-X_{0})=e_{n}X_{0}$ and noting that that the
$j^{th}$ derivative with respect to $t$:

$x_{0}^{(j)}(t)=j!t^{j}[x_{0}(t)]^{j+1}$

$\displaystyle x_{n}(t)=x_{0}(t)+(e_{n}X_{0}t)x_{0}^{(2)}(t)+$ (21)
$\displaystyle(e_{n}X_{0}t)^{2}x_{0}^{(3)}(t)+(e_{n}X_{0}t)^{3}x_{0}^{(4)}(t)+\ldots$

Substituting $x_{n}(t)$ into Eq. (13) and noting that $E_{1}=0$ we have

$\displaystyle
v(t)=x_{0}^{M}(t)+\frac{1}{2}(X_{0}t)^{M+2}E_{2}x_{0}^{M+2}(t)+$ (22)
$\displaystyle\frac{1}{3}(X_{0}t)^{M+3}E_{3}x_{0}^{M+3}(t)+\ldots$

In the above expression $x_{0}^{(}M)(t)$ represents the GF for $M$ balanced
stations. It follows from the definition of $E_{j}$’s that this series will
converge when the stations are not highly unbalanced and $e_{n}$s are small.

The coefficient of $t^{k}$ in Eq. (22) is given as:

$\displaystyle V(k)=V_{0}(k)+\frac{X_{0}^{2}}{2}V_{0}(k-2)E_{2}$ (23)
$\displaystyle+\frac{X_{0}^{3}}{3}V_{0}(k-3)E_{3}+\ldots$

$V_{0}(k=\binom{M+k-1}{k}X_{0}^{k}$ is the coefficient of $t^{k}$ in
$X_{0}^{M}(t)$.

We can rewrite Eq. (23) as

$\displaystyle V(k)=V_{0}(k)[1+R_{0}(k)],\mbox{ where }$ (24) $\displaystyle
R_{0}(k)\sum_{j=2}^{\infty}\frac{E_{j}}{j}\left[\prod_{i=0}^{j-1}\frac{k-i}{M+i}\right]$

The throughput of the subnetwork with $k$ jobs is given as

$\displaystyle
T(k)=\frac{V(k-1)}{V(k)}=T_{0}(K)\frac{1+R_{0}(k-1)}{1+R_{0}(k)}$ (25)

where $T_{0}(k)$ is the throughput of the $M$ balanced stations.

We can obtain a family of approximations for $T(k)$ according to the number of
terms in the summation for $R_{0}(k)$. Maintaining multiple terms in Eq. (23):

$\displaystyle V(k)=V_{2}(k)+V_{0}(k)R_{2}(k)$ (26)

$\displaystyle
V_{2}(k)=V_{0}(k)\left[1+\frac{1}{2}\frac{k(k-1)}{M(M-1)}E_{2}+\frac{1}{3}\frac{k(k-1)(k-2)}{M(M+1)(M+2)}E_{3}\right]$

$\displaystyle
R_{2}(k)=\sum_{j=4}^{\infty}\frac{E_{j}}{j}\prod_{i=0}^{j-1}\frac{k-i}{M+i}$

The expression for throughput with with one and two summation terms are:

$\displaystyle
T_{1}(k)=T_{0}(k)\left[\frac{1-\frac{k-1}{M+1}E_{2}}{1+\frac{k(k-1)}{2(M+1)}E_{2}}\right]$
(27)

$\displaystyle
T_{2}(k)=T_{0}(k)\left[1-\frac{\frac{k-1}{M(M+1)}\left(E_{2}+\frac{k-2}{M+2}E_{3}\right)}{1+\frac{k(k-1)}{M(M+1)}\left(\frac{1}{2}E_{2}+\frac{1}{3}\frac{k-2}{M+2}E_{3}\right)}\right]$
(28)

Given the coefficient of variation $c$, which is the standard deviation
divided by the mean and the coefficient of skewness, $\beta$ which is the
third central moment divided by $c^{3}$ we have $E_{2}=Mc^{2}$ and
$E_{3}=M\beta c^{3}$.
111https://grapherhelp.goldensoftware.com/WTOPICS/WKS_Skew.htm Eq. (28) and
Eq. (27) can be rewritten using this notation as:

$\displaystyle\boxed{T_{1}(k)=T_{0}(k)\left[1-\frac{\frac{k-1}{M+1}c^{2}}{1+\frac{k(k-1)}{2(M+1)}c^{2}}\right].}$
(29)

$\displaystyle\boxed{T_{2}(k)=T_{0}(k)\left[1-\frac{\frac{k-1}{M+1}c^{2}\left(1+\frac{k-2}{M+2}c\beta\right)}{1+\frac{k(k-1)}{M+1}c^{2}\left(\frac{1}{2}+\frac{k-2}{3(M+2)}c\beta\right)}\right]}$
(30)

In [39] we determine the relative errors: $T_{0}(k)/T_{2}(k)-1$ (resp.
$T_{1}(k)/T_{2}(k)-1$) for varying $c$ (resp. $c$ and $\beta$), with respect
to $T_{2}(k)$ since it is quite accurate.

$T_{0}(k)$ which is the throughput based on the average service demand $X_{0}$
yields the maximum throughput as shown in [43].

Given that there are $m_{n}$ servers at station $n$ with mean service demands
$X_{n}$ it follows from the fact the utilization factor should not exceed one:

$U_{n}(k)=\frac{T(k)X_{n}}{m_{n}}<1$

It is shown in [30] that the maximum throughput of a closed QN is given as the
minimum of service station throughputs, hence as $k$ is increased $T(k)$ is
bounded as:

$\displaystyle\boxed{T_{max}=\mbox{min}\left(\frac{m_{1}}{X_{1}},\frac{m_{2}}{X_{2}},\ldots,\frac{m_{N}}{X_{N}}\right)\geq
T(k),k\geq 1.}$ (31)

Since we are concerned with single-server queues $m_{n}=1,\forall{n}$ and
$T_{max}=1/D_{max}$.

In fact the throughput for lower values of $k$ the throughput is bounded by a
line that connects the origin and $T(1)=[\sum_{n=1}^{N}X_{n}]^{-1}$.
Throughput convexity is shown by Dowdy et al. 1984 [11] for QN’s with fixed
rate and delay servers.

Optimistic (resp. pessimistic) Balanced Job Bound - BJB analysis of QNs
utilize the average (resp. maximum) service demand yielding an upper (resp.
lower) bounds to throughput [43],[28] Since the optimistic BJB bound may
exceed $T_{max}$ it should only be used up to the point it intersects
$T_{max}$.

$\displaystyle\frac{K}{(K+N-1)D_{max}}\leq T(k)\leq$
$\displaystyle\mbox{min}\left(\frac{1}{D_{max}},\frac{K}{(K+N-1)D_{avg}}\right).$

It follows easily from Eq. (31) that when the sum of service demands is fixed
then the throughput is maximized when the service demands at the stations
equal the mean.

From a pragmatic viewpoint there is a drop from the maximum throughput as the
MultiProgramming Level - MPL is increased in virtual memory systems. This is
due to increased paging which may ensue in a severe performance degradation.
The frequency of paging is given by the lifetime curve, which is the time
between page faults and the number of page frames dedicated to a program (see
e.g., Figure 9.5 in [28]). As MPL is increased less page frames can be
allocated to each program A QN model is used to obtain the throughput
characteristic vs the MPL, which show a drop beyond a certain maximum
throughput. Increasing the MPL further may result in a severe degradation in
throughput referred to as thrashing Denning 1970 [9].

##  5 Aggregation of Balanced Stations with Two Classes

The aggregation of balanced stations for a single class is extended to two job
classes in the Appendix of [39] and Nadji and Thomasian 1984 in [31]. There
are $N$ stations and $K$ and $L$ jobs in two classes with service demands
given by $X_{n}=X,1\leq n\leq N$ and $Y_{n}=Y,1\leq n\leq N$. The equilibrium
state probability distribution is given as:

$\displaystyle P[{\bf
k},\underline{\ell}]=\frac{1}{G(K,L)}\prod_{n=1}^{N}F_{n}(k_{n},\ell_{n}).$
(32)

For single-server stations with FCFS, PS, and LCFSPR queueing disciplines: In
the first case both job classes should have the same exponentially distributed
service times.

$\displaystyle
F_{n}(k_{n},\ell_{n})=\frac{(k+\ell)!}{k_{n}!\ell!}X_{n}^{k_{n}}Y_{n}^{\ell_{n}}.$
(33)

The GF for each station is given by:

$\displaystyle h_{n}(t,u)$
$\displaystyle=\sum_{k=0}^{\infty}\sum_{\ell=0}^{\infty}\frac{(k+\ell)!}{k!\ell!}(X_{n}t)^{k}(Y_{n}u)^{\ell}$
(34) $\displaystyle=(1-X_{n}t+Y_{n}u)^{-1}\mbox{ provided }|X_{n}t+Y_{n}u|<1.$

In the case of IS stations

$F_{n}(k,\ell)=\frac{(X_{n}t)^{k}}{k!}\frac{(Y_{n}u)^{\ell}}{\ell!}$

$\displaystyle
h_{n}(t,u)=\sum_{k=0}^{\infty}\sum_{\ell=0}^{\infty}\frac{(X_{n}t)^{k}}{k!}\frac{(Y_{n}u)^{\ell}}{\ell!}=$
(35) $\displaystyle\sum_{m=0}^{\infty}(X_{n}t+Y_{n}u)^{m}=e^{X_{n}t+Y_{n}u}.$

The GF for the $M\leq N$ stations to be aggregated is

$\displaystyle v(t,u)=\prod_{i=1}^{M}h_{i}(t,u).$

When the QN is balanced $X_{n}=X,Y_{n}=Y,1\leq n\leq N$. The GF for the $M$
stations to be aggregated is:

$\displaystyle v(t,u)=[x(t,u)]^{M}=(1-Xt-Yu)^{-M}$ (36)
$\displaystyle\sum_{j=0}^{\infty}\binom{M+j-1}{j}(X_{0}t+Y_{0}u)^{j}.$

$\displaystyle
v(t,u)=\sum_{k=0}^{\infty}\sum_{\ell=0}^{\infty}V(k,\ell)(MXt)^{k}(MYu)^{\ell}.$
(37)

Equating the coefficients of $t^{k}u^{\ell}$ of $v(t,u)$

$\displaystyle V(k,\ell)=\binom{M+k+\ell-1}{k+\ell}\binom{k+\ell}{k}.$ (38)

The throughputs in the two classes are given as:

$\displaystyle
T^{(1)}(K,L)=\frac{V(K-1,L)}{V(K,L)}=\frac{K}{M+K+L-1}.\frac{1}{X_{0}}$ (39)
$\displaystyle
T^{(2)}(K,L)=\frac{V(K,L-1)}{V(K,L)}=\frac{L}{M+K+L-1}.\frac{1}{Y_{0}}$

If the QN is not balanced we expand the GF for ${\cal S}_{n}$ around the mean
loadings $X_{0}$ and $Y_{0}$:

$\displaystyle h_{n}(t,u)=\sum_{j=0}^{\infty}\frac{1}{j!}\hskip 128.0374pt$
(40) $\displaystyle\left[\left((X_{n}-X_{0})\frac{\partial}{\partial
X_{n}}+(Y_{n}-Y_{0})\frac{\partial}{\partial
Y_{n}}\right)^{j}h_{n}(t,u)\right]=$
$\displaystyle\sum_{j=0}^{\infty}\left[(X_{n}-X-0)t+(Y_{n}-Y_{0})u\right]^{j}h_{0}^{j+1}(t,u).$

$\displaystyle v(t,u)/h_{0}^{M}(t,u)=\hskip 128.0374pt$ (41) $\displaystyle
1+\frac{1}{2}h_{0}^{2}(t,u)\sum_{n=1}^{M}[(X_{n}-X_{0})t+(Y_{n}-Y_{0})u]^{2}+$
$\displaystyle\frac{1}{3}h_{0}^{3}(t,u)\sum_{n=1}^{M}[(X_{n}-X_{0})t+(Y_{n}-Y_{0})u]^{3}+\ldots$

Equating the coefficients of $t^{k}u^{\ell}$

$\displaystyle V(K,L)=V_{0}(K,L)+$ (42) $\displaystyle
X_{0}^{K}Y_{0}^{L}\sum_{i=2}^{\infty}\frac{1}{i}\sum_{j=0}^{i}\left[\binom{M+K+L-1}{K+L-1}\binom{K+L-i}{K-j}\binom{i}{j}E_{j,i-j}\right]$

$E_{i,j}=\sum_{n=1}^{M}\left(\frac{X_{n}-X_{0}}{X_{0}}\right)^{i}\left(\frac{Y_{n}-Y_{0}}{Y_{0}}\right)^{j}$

The first order approximation can be expressed as:

$\displaystyle V_{1}(K,L)=V_{0}(K,L)[1+\frac{K(K-1)}{2M(M+1)}E_{2,0}+$ (43)
$\displaystyle\frac{L(L-1)}{2M(M+1)}E_{0,2}+\frac{KL}{M(M+1)}E_{1,1}].$

Given that $c_{X}$, $c_{Y}$, and $c_{X,Y}$ are the coefficients of variation
and correlation for $X_{n}$ and $Y_{n}$ then:

$\displaystyle E_{2,0}=Mc^{2}_{X}\hskip 5.69054ptE_{0,2}=Mc^{2}_{Y}\hskip
5.69054ptE_{1,1}=Mc_{X,Y}c_{X}c_{Y}.$

$\displaystyle T_{1}^{(1)}(K,L)=T_{0}^{(1)}(K,L)$ (44)

$\displaystyle\frac{1-\frac{1}{2(M+1)}[K(K-1)(K-2)c^{2}_{X}+L(L-1)c^{2}_{Y}+2(K-1)Lc_{X}c_{Y}c_{X,Y}]}{1+\frac{1}{2(M+1)}[K(K-1)c^{2}_{X}+L(L-1)c^{2}_{Y}+2KLc_{X}c_{Y}c_{X,Y}]}$

### Accuracy of Unbalanced Job Approximation

A small experiment to assess the accuracy of UJA is given below. The service
demands at the $N=5$ stations are as follows:

$\underline{X}=\\{0.25,0.23,0.19.17,0.15\\},$

such that $X_{0}=0.2$. $c=0.179$, and $\beta=0.079$. Exact results were
obtained using CA.

$k$ | $T_{0}(k)$ | $T_{1}(k)$ | $T_{2}(k)$ | $T(k)$
---|---|---|---|---
1 | 1.0000 | 1.0000 | 1.0000 | 1.0000
2 | 1.6667 | 1.6578 | 1.6578 | 1.6578
3 | 2.1429 | 2,1204 | 2.1203 | 2.1203
4 | 2.5000 | 2.4612 | 2.4611 | 2.4610
5 | 2,7778 | 2,7215 | 2.7212 | 2.7209
6 | 3.0000 | 2.9259 | 2.9254 | 2.9245
Table 1: Comparison of approximate and exact throughputs for a single job
class.

A larger experiment with 100,000 random service demands in a QN with $N=18$
stations is reported in [31]. More than 99% of the results obtained by $T_{1}$
and $T_{2}$ throughputs are within 10% accuracy for the range of number of
jobs considered. $T_{0}$ is accurate within 15% in 97% of the examples when
average utilization was 40%. An experiment similar to Table 1 with two job
classes is reported in [31].

##  6 Conclusions

UJA - Unbalanced Job Approximation - UJA is a family of low-cost methods
starting with an upper bound for the throughput based on the mean service
demand $X_{0}$, but more accurate throughput estimates are obtained using
additional terms in a Taylor series expansion,

Given a large number of stations, instead of aggregating all stations at once,
aggregation can be applied to groups of single-server stations with close
utilizations.

Given that UJA is restricted to two classes and that it is difficult to extend
to more classes, more job classes can be dealt with by clustering multiple
chains into one, see e.g., deSouza et al. 1986 [10]. The error introduced by
applying such methods is discussed in Cheng and Muntz 1986 [8].

Performance bound hierarchies for single and multiple job classes are given by
Eager and Sevcik [13, 14]. Kerola 1896 [23] is a less influential work on this
topic.

## Appendix I: Performance Bound Hierarchies - PBH

Eager and Sevcik 1983 [13] extended BJB to Performance Bound Hierarchies - PBH
Consider a closed QN with $K$ stations, $N$ job and loadings $L_{k},1\leq
j\leq K$ According to the arrival theorem [27]

$\displaystyle R_{k}(N)=L_{k}[1+\bar{n}_{k}(N-1)]$ (45)

where $\bar{n}_{k}(N-1)$ is the mean number of jobs at station $k$. It follows
the mean residence time of jobs is $R(N)=\sum_{k=1}^{K}R_{k}(N).$ An
application of Little’s result yields.

$\displaystyle\bar{n}_{k}(N-1)=\frac{R_{k}(N_{1})}{Z+R(N-1)}(N-1)$ (46)

Substituting $n_{k}(N-1)$ in Eq. (45)

$\displaystyle
R_{k}^{(i)}(N)=L_{k}[1+\frac{R_{k}^{(i)}(N_{1})}{Z+R^{(i-1)}(N-1)}](N-1)$ (47)

The optimistic hierarchy starts with the initialization:

$\displaystyle R_{k\mbox{ opt}}^{(0)}=\frac{1}{K}\mbox{max}[NL_{B}-Z,1]$ (48)

where $b$ is the index of the bottleneck station. It follows:

$\displaystyle R_{k\mbox{ opt}}^{(0)}=\frac{1}{K}\mbox{max}[NL_{b}-Z,1]$ (49)

Using $a(N)=\mbox{max}[NL_{b}-Z,1]$

$\displaystyle
R_{\mbox{opt}}^{(1)}=1+\frac{1}{K}\left(\frac{a(N-1)}{Z+a(N-1)}\right)(N-1).$
(50)

When $Z=0$ and $S=\sum_{k=1}^{K}L^{2}_{k}$

$\displaystyle R^{(2)}_{\mbox{opt}}(N)=1+S(N-1).$ (51)

For the pessimistic hierarchy we have

$\displaystyle R^{(0)}_{k\mbox{pess}}=\begin{cases}N\mbox{ for }k=b\\\ 0\mbox{
for }k\neq b\end{cases}$ (52)

$\displaystyle R^{(0)}_{\mbox{pess}}(N)=N,$ (53)

For $k\neq b$ $R_{k}{\mbox{p}ess}^{(1)}=L_{a}$ and for $k=b$

$\displaystyle
R_{k}{\mbox{p}ess}^{(1)}(N)=1+L_{b}\left(\frac{N-1}{Z+N-1}\right)(N-1)$ (54)

When $Z=0$ the bound corresponds to BJB pessimistic bound

Level 2 pessimistic bound for $Z=0$ is

$\displaystyle
R^{(2)}_{pess}(N)=1+\left(\frac{L_{b}^{2}(N-2)}{1+L_{b}(N-2)}\right)(N-1)$
(55)

A measure of the error magnitude is

$\displaystyle\frac{R^{(i)}(N)-R^{(i)}(N)}{R^{(i)}(N)+R^{(i)}(N)}\times 100\%$
(56)

Eager and Sevcik 1986 [14] is an extension of the single class PBH bounding
algorithm to QNs with multiple classes of customers. The resulting algorithm
is applicable to all of the types of product-form QNs that are commonly used
in computer system and computer-communication network applications. Asymptotic
Expansions - AE algorithm proposed by McKenna and Mitra 1994 [29] is also
capable of providing reasonably tight bounds for a practically useful subset
of multiple class closed QNs. AE algorithm can only treat QNs in which each
class has a significant service demand at one or more IS servers. Both
algorithms offer a hierarchy of bounds with differing accuracy levels and
computational cost requirements.

## Appendix II: Kriz’s Extension to BJB

The discussion follows Kriz 1984 [26] and preserves its notation. The number
of jobs is $N$ Single-server or waiting stations are denoted by (${\cal W}$).
IS or delay stations are denoted ${\cal D}$ and the sum of their service
demands by $Z$. The set of all servers is denoted by ${\cal Q}$ with $M=|{\cal
Q}|$. The loading at server $m$ is $t_{m}=v_{m}s_{m}$, where $v_{m}$ is the
number of visits and $s_{m}$ the service time. The sum of all loading at the
${\cal W}$ servers is $T$ and the sum of all loadings is $r=Z+T$. The maximum
and average loading at ${\cal W}$ servers is $t^{\prime}$ and $t"$.

$L_{m}(N)$ is the queue-length, $R_{m}(N)$ the response time, and $X_{m}(N)$
the throughput at $m^{th}$ station. The system throughput is
$X(N)=X_{m}(N)/v_{m}$ and the mean residence time in the system is
$R(N)=N/X(N)$.

The throughput ABA and BJB can be expressed succinctly as:

$X(N)\leq\mbox{min}(N/r,1/t^{\prime})$

$\frac{N}{r+(N-1)t^{\prime}}\leq X(N)\leq\frac{N}{r+(N-1)t"}$

Given that that load at ${\cal W}$ stations is balanced with loadings given by
$t$,

$X(1)=1/r,\hskip 8.53581ptX(n)=n/[r+(n-1-ZX(n-1)t],n\geq 2$

With initialization $\underline{Y}_{0}(N)=0$ we have:

$\underline{Y}_{i}(N)=N[r+(N-1-ZY_{i-1}(N-1))t^{\prime}]^{-1},1\leq i\leq N$

Given $\bar{Y}_{0}(N)=\mbox{min}\\{N/r,1/t^{\prime}\\}$

$\bar{Y}_{i}(N)=\mbox{min}\left(N[r+(N-1-Z\bar{Y}_{i-1}(N-1)t"]^{-1},1/t^{\prime}\right)$

Then $\bar{Y}_{i-1}(N)\leq\bar{Y}_{i-1}(N)\leq
X(N)\leq\bar{Y}_{i}(N)\leq\bar{Y}_{i-1}(N)$

## Appendix III: Asymptotic Expansions with Multiple Bottlenecks

Considered in George et al. 2012 [15] is a closed QN with $N$ jobs and $M$
stations with ${\cal M}$ denoting the set of stations. ${\cal S}$ denotes
single-servers, ${\cal L}$ multiservers, and ${\cal I}$ IS stations. Given the
routing matrix $P=[p_{i,j}]$, the relative throughputs are:

$\pi=\sum_{j\in{\cal M}}\pi_{j}p_{j,i}.$

For any sequence $\\{f(n),n\geq 0\\}$ its z-transform in alternate form is:

$\tilde{f}(z)=\sum_{n=0}^{\infty}f(n)z^{-n}$

Given three theorems regarding the properties of $\tilde{f}(z)$ and an
assumption is followed by a theorem.

Assumption 1. We assume that the relative throughputs $\pi_{i},\in{\cal M}$
are chosen such that $\mbox{max}\\{\rho_{i}=1\\}$ stations are labeled by
their relative utilization, such that:

$1=\rho_{1}\geq\rho_{2}\geq\ldots\geq\rho_{M}$

The set of bottleneck stations is defined as ${\cal B}:=\\{i\in{\cal
M}:\rho_{i}=1\\}$. Clearly $|B|\geq 1$.

The main theorem that characterizes the asymptotic behavior of the
normalization constant $G(N)$ in exact order.

Theorem 4. For any $M$ station closed QN with $N$ jobs the normalization
constant $G(N)$ satisfies the exact asymptotics:

$\displaystyle G(N)\sim C_{B}N^{|B|-1}\mbox{ where }$ $\displaystyle
C_{B}=\frac{1}{(|B|-1)!}\prod_{i=1}^{|B|}\frac{\gamma_{i}^{s^{i}}}{s_{i}!}\prod_{j=|B|+1}^{M}\tilde{f}_{j}(1)$
(57)

As $N$ grows the rate at which $G(N)$ grows is on the order $N^{|B|-1}$. The
actual throughput $\Lambda_{i}(N)$ and utilization $U_{i}(m)$ satisfy the
following exact-order asymptotics:

$\Lambda_{i}(N)\sim\pi_{i}(1-\frac{1}{N})^{|B|-1},\forall{i}\in{\cal M}$

$U_{i}(N)\sim\rho_{i}(1-\frac{1}{N})^{|B|-1},\forall{i}\in{\cal M^{\prime}}$

where ${\cal M}$ is the set of all stations and ${\cal M^{\prime}}$ excludes
IS stations (${calI}$).

An upper bound to system throughput is given as:

$\Lambda(N)\leq\mbox{min}\left[\frac{N}{D+Z},s_{1}\mu_{1}\right]:=\Lambda^{ABA}(N)$

where $D=\sum_{i\in\cal M^{\prime}}\pi_{i}/\pi_{1}\mu_{i}$ and
$Z=\sum_{i\in{\cal I}}\pi_{i}/\pi_{1}\mu_{i}$.

As $N$ increases a saturation point is reached where the bottleneck resource
reaches 100% utilization.

Noting that $\Lambda_{i}(N)=\pi_{i}/\pi_{1}\lambda(N)$

$\displaystyle\mbox{min}\left[\frac{N}{D+Z},s_{1}\mu_{1}(1-\frac{1}{N})^{|B|-1}\right]$
(58)

Numerical examples are given why the AE bound is more accurate than others.

We note that in the case of a network with a single bottleneck, Eq. (58)
reduces to the previous ABA bound and the proposed approximation provides the
greatest improvements

### Appendix IV: Proportional and Geometric Bounds

The Geometric Bounds - GBs proposed by Casale et al. 2006/08 [6, 7] are fast
and accurate noniterative bounds on closed QN metrics. Compared to BJB GB
achieves higher accuracy at similar computational cost, limiting the worst-
case bounding error typically within 5-13% while for BJB the error is usually
in the range of 15-35%. 222Application of GB as an approximation to fork-join
processing in closed QNs is beyond the scope of this discussion.

We restate the steps for MVA [33] in a closed QN with $M$ queues with service
demands or loadings $L_{i},1\leq i\leq M$ with $L=\sum_{i=1}^{M}L_{i}$. The
delay at the delay server is $Z$. The network throughout with $N$ jobs is
$X(N)$. The utilization at the $i^{th}$ queue is $U_{i}(N)$ and the queue-
length $Q_{i}(N)$. It follows from Little’s result

$U_{i}(N)=L_{i}X(N),1\leq i\leq M$

According to the arrival theorem [27] the mean residence time at the $i^{th}$
station is:

$W_{i}(N)=L_{i}[1+Q_{i}(N-1)],1\leq i\leq M$

The residence time in the network is:

$R(N)=\sum_{i=1}^{M}W_{i}(N)$

The network throughput follows from Little’s result:

$X(N)=N/(Z+R(N)$

The mean number of jobs at the $i^{th}$ queue is:

$Q_{i}=X(N)W_{i}(N)=U_{i}(N)[1+Q_{i}(N-1)],1\leq i\leq M$ 333The paper uses
$R_{i}(N)$ which is undefined, it should be $W_{i}(N)$.

By recursively expanding the above equation:

$\displaystyle Q_{i}(N)=U_{i}(N)+U_{i}(N)U_{i}(N-1)+$ $\displaystyle
U_{i}(N)U_{i}(N-1)U_{i}(N-2)+$ $\displaystyle\prod_{i=0}^{N-1}U_{i}(N-n)$

The geometric summation yields:

$\displaystyle\sum_{i=1}^{N}y^{i}=\frac{y(1-y^{N})}{1-y}$ (59)

The throughput is bounded as follows where $X_{max}$ is the maximum throughput
obtained by ABA:

$\displaystyle N/(Z+LN)\leq X(N)\leq$ (60)
$\displaystyle\mbox{min}(N/(Z+L),X_{Max})\mbox{ ABA }$

The throughput for BJBs is:

$\displaystyle N/(Z+L+L_{max}(N-1-ZX^{-1})\leq X(N)\leq$ (61) $\displaystyle
N/(Z+L+L(N-1)ZX^{+})/M)\mbox{ BJB }$

The throughput for Proportional Bounds - PB Hsieh and Lam 1987 [19].

$\displaystyle
N/\left(Z+L+\frac{\sum_{i=1}^{M}L_{i}^{N}(N-1-ZX^{-})}{\sum_{j=1}^{M}L_{j}^{N-1}}\right)\leq
X(N)$ (62) $\displaystyle\leq
N/\left(Z+L+\frac{\sum_{i=1}^{N}L_{i}^{2}(N-1-ZX^{+})}{\sum_{j=1}^{M}L_{j}}\right)\mbox{
PB }$

BJB always offers greater accuracy than the ABA. The bounds hold true for any
$X^{+}$ and $X^{-}$ such that

$X(N-1)\leq X^{+}\leq X_{max},\hskip 14.22636pt0\leq X^{-}\leq X(N-1)$

The following formula is exact:

$\displaystyle X(N)=N/[Z+L+L_{max}(N-1)-ZX(N-1))-D(N)]$ (63)

where

$D(N)=\sum_{i=1}^{N}(L_{max}-L_{i})Q_{i}(N-1).$

Theorem 4. $X(n)$ for $Z>0$ is bounded by

$\displaystyle 2N/\left(b(N)+\sqrt{b^{2}(N)-4ZL_{max}(N-1)}\right)\leq
X(N)\leq$ (64) $\displaystyle 2N/\left(b(N)+\sqrt{b^{2}(N)-4ZL_{max}N}\right)$

$\displaystyle b(N)=Z+L+L_{max}(N-1)-\sum_{i:L_{i}<L_{max}}$ (65)
$\displaystyle(L_{max}-L_{i})Q_{i}(N-1).$

$Q_{i}$ is bounded by $Q_{i}^{-}$ and $Q_{i}^{+}$ given by Eq. (Appendix IV:
Proportional and Geometric Bounds and Eq. (Appendix IV: Proportional and
Geometric Bounds), respectively.

Theorem 1. The queue length of station is bounded from below by

$\displaystyle
Q_{i}^{-}(N)=\begin{cases}\frac{y_{i}(N)-{y_{i}(N)}^{N+1}}{1-y_{i}(N)}\mbox{
if }L_{i}<L_{max},\\\
\frac{1}{m_{max}}\left(N-ZX^{+}-\sum_{k:L_{k}<M_{max}}Q^{+}_{k}(N)\right)\\\
\mbox{ if }L_{i}=L_{max}\end{cases}$

For any $X^{+}$ such that $X(N)<X^{+}$ and where

$y_{i}(N)=L_{i}N/(Z+L+L_{max}N)$

is the ratio of the underlying sum in Eq. (59), $M_{max}$ is the number of
queues with service demand $L_{max}$, and $Q^{+}_{k}$ is the upper bound in
Theorem 2 for $L_{i}<L_{max}$

Theorem 2. The queue-length $Q_{i}(N)$ is bounded above as follows:

$\displaystyle
Q^{+}_{i}(N)=\begin{cases}\frac{Y_{i}(N)[1-{Y_{i}(N)}^{N}]}{1-Y_{i}(N)}\mbox{
if }L_{i}<L_{max}\\\
\frac{1}{M_{max}}\left(N-ZX^{-}-\sum_{k:L_{k}<L_{max}}Q^{+}_{k}(N)\right)\\\
\mbox{ if }L_{i}=L_{max}.\end{cases}$

For any $X^{+}$ and $X^{-}$ such that $X(N)\leq X^{+}\leq X_{max}$ and $0\leq
X^{-}\leq X(N)$ and where $Y_{i}(N)=L_{i}^{+}$ is the ratio $y$ of the
underlying geometric sum.

$\frac{N-1}{N}X(N)\leq X(N-1)\leq X(N)$

$X^{-}$ and $X^{+}$ are used to represent lower and upper bounds.

The it Geometric Square Bound - GSB is obtained by replacing queue=lengths by
$Q_{i}^{-}(N-1)$.

Two classes of performance bounds one for single chain and the other for
multi-chain QNs are proposed in [19]. Proportional Bounds - PBs assume that
mean queue-lengths are proportional to server loads. PBs are more accurate
than BJBs in that individual server loads are retained as parameters in bounds
formula. Several tables are given in [7] comparing the accuracy of various
methods. GSB is shown to be the most accurate even when one or two iterations
are allowed for BJB and PB.

## Appendix V: Proportional Approximation Method - PAM

Proportional Approximation Method - PAM is that latest bounding method for
closed multichain QNs by Hsieh and Lam 1988/89 [20, 21].

Approximate MVA algorithms for separable queueing networks are based upon an
iterative solution of a set of modified MVA formulas. Let $M$ denotes the
number of queues and $K$ the number of chains or job classes. Each iteration
has a computational cost of O($MK^{2}$) or less, many iterations are typically
needed to attain convergence.

Presented are faster approximate noniterative solution algorithms, which are
suitable for the analysis and design of communication networks requiring
thousands of closed chains.

Three PAM algorithms of increasing accuracy are presented. Two of them have
time and space requirements of O($MK$). The third algorithm has a time
requirement of O($MK^{2}$) and a space requirement of O($MK$).

Three PAMs developed in this study consider QNs with fixed-rate and delay
servers.

Let $\tau_{mk}$ denote the loading of chain $k$ at queue $m$.

Algorithm PAM_BASIC

Step 1. Calculate proportional approximations of mean queue lengths.

$\gamma_{mk}=\tau_{mk}/\sum_{i=1}^{M}\tau_{ik}$

$\mbox{for }m=1,\ldots,M,\mbox{ and }k=1,\ldots,K$

${q^{\prime}}_{mh}(\underline{N})=\gamma_{nh}N_{h}$

${q^{\prime}}_{mh}(\underline{N}-\underline{1}_{K})=\begin{cases}{q^{\prime}}_{mh}(\underline{N})\mbox{
if }h\neq k.\\\ {q^{\prime}}_{mh}(\underline{N})-\gamma_{mh}\mbox{ if }h=k.\\\
\end{cases}$

$\mbox{for }m=1,\ldots,M,\mbox{ and }k=1,\ldots,K$

Step 2: Calculate approximate mean delay of chain $k$ at server $m$ and
approximate throughput of chain $k$

$D_{mk}(\underline{N})=\begin{cases}\tau_{mk}[1+\sum_{k=1}^{K}{q^{\prime}}_{mh}(\underline{N})-\underline{1}_{k})]\mbox{
if fixed rate}\\\ \tau_{mk}\mbox{ if delay server}\end{cases}$

$T_{k}(\underline{N})=\frac{N_{k}}{\sum_{m=1}^{M}D_{mk}(\underline{N})}\mbox{
for }k=1,\ldots K$

The accuracy improvement of PAM_IMPROVED over PAM_BASIC is obtained by a
simple scaling operation to ensure that server utilizations do not exceed one.

Algorithm PAM_Improved calculates server utilizations and finds the the server
with maximum utilization visited by chain $k$.

Step 3.

$U_{m}(\underline{N})=\sum_{k=1}^{K}T_{k}(\underline{N}),\hskip
8.53581ptm=1,\ldots,M$

Step 4. Find the largest utilization $S_{k}$ visited by chain $k$.

$S_{k}=\mbox{max}\\{U_{m}(\underline{N})\\}\hskip 8.53581ptk=1,\ldots,K$

Step 5. Scale down the throughputs of individual chains in necessary

for $k=1,\ldots,K$

$\mbox{ If }S_{k}>1\mbox{ then
}T_{k}(\underline{N})=T_{k}(\underline{N})/S_{k},$

Step 6. Calculate total throughput and utilizations

$T(\underline{N})=\sum_{k=1}^{K}T_{k}(\underline{N})$

$U_{m}(\underline{N})=\sum_{k=1}^{K}\tau_{ml}T_{k}(\underline{N}),\hskip
8.53581ptm=1,\ldots,M$

Algorithm PAM_II

The accuracy of the PAM algorithm is improved by executing the last two steps
of the MVA recursion instead of just the last step using the proportional
approximation to get initial mean queue length estimates. The reader is
referred to [20, 21] Relative errors in chain throughputs for 500 networks
calculated by PAM_Improved and PAM_TWO are 2.3% vs 0.8% for average error and
40.3 vs 30.8% for the maximum.

The additional accuracy of PAM_TWO is obtained by executing the final two
steps of the MVA recursion instead of just the last step. The computational
time requirements are O(MK) for PAM_BASIC and PAM_IMPROVED, and O($M^{2}$) for
PAM_TWO. All three PAM algorithms have space requirements of O(MK).

### Acknowledgement

The discussion of UJA is based on joint publications with Dr Behzad (Brad)
Nadji at the EE-Systems Dept. at USC.

## References

  * [1] G. Balbo and G. Serazzi. Asymptotic analysis of multiclass closed queueing networks: Common bottleneck. Performance Evaluation 26, 1 (1996), 51-72.
  * [2] G. Balbo and G. Serazzi. Asymptotic analysis of multiclass closed queueing networks: Multiple bottlenecks. Performance Evaluation 30, 3 (1997), 115-152.
  * [3] F. Baskett, K. M. Chandy, R. R. Muntz, and F. G. Palacios. Open, closed, and mixed networks of queues with different classes of customers. J. ACM 22, 2 (1975), 248-260.
  * [4] J. P. Buzen. Computational algorithms for closed queueing networks with exponential servers. Commun. ACM 16, 9 (1973), 527-531.
  * [5] J. P. Buzen, R. P. Goldberg, A. M. Langer, E.S. Lentz, H. S. Schwenk, D. A. Sheetz and A. W. Shum. BEST/1 - Design of a tool for computer system capacity planning. In Proc. AFIPS National Computer Conf. 1978: 447-455
  * [6] G. Casale, R. R. Muntz, and G. Serazzi: A new class of non-iterative bounds for closed queueing networks. In Proc. 14th Int’l Symp. on Modeling, Analysis, and Simulation of Computer and Telecomm. Systems \- MASCOTS 2006: 69-76.
  * [7] G. Casale, R. R. Muntz, and G. Serazzi. Geometric bounds: A noniterative analysis technique for closed queueing networks. IEEE Trans. Computers 57, 6 (June 2008), 780-794.
  * [8] W. C. Cheng and R. R. Muntz. Bounding errors introduced by clustering of customers in closed product-form queueing networks. J. ACM 43, 4 (1996), 641-669.
  * [9] P. J. Denning, Virtual Memory. ACM Computing Surveys 2, 3 (1970), 153-189.
  * [10] E. de Souza e Silva, S. S. Lavenberg, and R. R. Muntz. A clustering approximation technique for queueing network models with a large number of chains. IEEE Trans. Computers 35, 5,(May 1986), 419-430.
  * [11] L. W. Dowdy, D. L. Eager, K. D. Gordon, and L. V. Saxton. Throughput concavity and response time convexity. Information Processing Letters 19, 4 (1984), 209-212.
  * [12] L. W. Dowdy, B. M. Carlson, A. T. Krantz, and S. K. Tripathi, Single-class bounds of multi-class queuing networks. J. ACM 39, 1 (1992), 188-213.
  * [13] D. L. Eager and K. C. Sevcik. Performance bound hierarchies for queueing networks. ACM Trans. Computer Systems 1, 2 (May 1983), 99-115.
  * [14] D. L. Eager and K. C. Sevcik. Bound hierarchies for multiple class queueing Networks. J. ACM 33, 1 (May 1986) 178-185.
  * [15] F. George, C. Xia, and M. Squillante. Exact-order asymptotic analysis for closed queueing networks. J. Applied Probability 49, 2 (2012), 503-520.
  * [16] B. M. Gijsen, R .D.van der Mei, P. Engelberts, J.L.van den Berga, and K.M. C.van Wingerden. Sojourn time approximations in queueing networks with feedback. Performance Evaluation 63, 8 (Aug. 2996), 743-758.
  * [17] W. J. Gordon and G. F. Newell. Closed queueing systems with exponential servers. Operations Research 15, 2 (April 1967), 254-265.
  * [18] P. G. Harrison et al. On the asymptotic behavior of closed multiclass queueing networks, Performance Evaluation 47, 2/3 (2002), 131-138.
  * [19] C.-T. Hsieh and S. S. Lam. Two classes of performance bounds for closed queueing networks. Performance Evaluation 7 (1987), 3-30.
  * [20] C.-T. Hsieh and S. S. Lam: PAM - A noniterative approximate solution method for closed multichain queueing networks. ACM SIGMETRICS Performance Evaluation Review 16, 1 (May 1988), 261–-269.
  * [21] C.-T. Hsieh and S. S. Lam. PAM - A noniterative approximate solution method for closed multi-chain queueing networks. Performance Evaluation 9, 2 (1989), 119-133.
  * [22] J. R. Jackson. Networks of waiting lines. Operations Research 5, 4 (1957), 518-521.
  * [23] T. Kerola. The composite bound method for computing throughput bounds in multiple class environments. Performance Evaluation 6, 1 (1986), 1-9.
  * [24] L. Kleinrock. Queuing Systems Vol. I: Theory. Wiley-Interscience, 1975.
  * [25] L. Kleinrock. Queuing Systems Vol. II: Computer Applications. Wiley-Interscience, 1976.
  * [26] J. Kriz. Throughput Bounds for Closed Queueing Networks Performance Evaluation 4(1): 1-10 (Feb. 1984),
  * [27] S. S. Lavenberg and M. Reiser. Stationary state probabilities of arrival instants for closed queueing networks with multiple customer types. J. Applied Probability 17 (1980), 1948-1961.
  * [28] E. D. Lazowska, J. Zahorjan, G. Scott Graham, and K. C. Sevcik: Quantitative System Performance: Computer System Analysis Using Queueing Network Models Prentice-Hall 1984.
  * [29] J. McKenna and D. Mitra: D Asymptotic expansions and integral representations of moments of queue lengths in closed Markovian networks. J. ACM 31, 2 (1984), 346-360.
  * [30] R. R. Muntz and J. Wong. Asymptotic Properties of Closed Queueing Network Models. In Proc. 8th Annual Princeton Conf. on Informations Sciences and Systems, 1974: 348-352.
  * [31] B. Nadji and A. Thomasian. Throughput estimation in queueing networks. Computer Performance 5, 4 (Dec. 1984), 197-206.
  * [32] M. Reiser and H. Kobayashi. Queuing networks with multiple closed chains: Theory and computational algorithms. IBM J. Research & Development 19, 3 (1975), 283-294.
  * [33] M. Reiser and S. S. Lavenberg: Mean-value analysis of closed multi-chain queuing networks. J. ACM 27, 2 (1980), 313-322.
  * [34] P. J. Schweitzer. Approximate analysis of multiclass closed networks of queues. In Proc. Int’l Conf. on Stochastic Control and Optimization, 1979.
  * [35] P. J. Schweitzer, G. Serazzi, and M. Broglia. A survey of bottleneck analysis in closed networks of queues. Performance/SIGMETRICS Tutorials, Lecture Notes in Computer Science - LNCS 729, Springer 1993: 491-508
  * [36] J. Silvester and A. Thomasian. Performance modeling of a large-scale multiprogrammed computer using BEST/1. In Proc. 3rd Int’l. Conf. on Computer Capacity Management - ICCCM, 1981, pp. 197-211.
  * [37] R. Suri, S. Sahu, and M. Vernon. TBD Approximate mean value analysis for closed queuing networks with multiple-server stations.. Proc. 2007 Industrial Engineering Research Conf.
  * [38] A. Thomasian and B. Nadji. Aggregation of stations in queueing network models of multiprogrammed computers. In Proc. ACM SIGMETRICS Conf. on Measurement and Modeling of Computers 1981, 86-96.
  * [39] A. Thomasian and B. Nadji. Algorithms for queueing network models of multiprogrammed computers using generating functions. Computer Performance 2(3): 100-123 (Sept. 1981), IPC Press, UK.
  * [40] K. S. Trivedi. Probability and Statistics with Reliability, Queueing, and Computer Science Applications, 2nd ed. Wiley-Interscience 2002.
  * [41] A. C. Williams and R. A. Bhandiwad. A generating function approach to queueing network analysis of multiprogrammed computers. Networks 6, 1 (1976), 1-22.
  * [42] J. Zahorjan, K. Sevcik, D. Eager, and B. Galler. Balanced job bound analysis of queueing net- works. Abstract in SIGMETRICS’81, 58.
  * [43] J. Zahorjan, K. Sevcik, D. Eager and B. Galler. Balanced job bound analysis of queueing networks. Commun. ACM, 25 2 (Feb. 1982), 134-141.