<figure><img src="image_1.png" /></figure>

HAL
open science

# Discrete event simulation modelling of computer system for performance evaluation

Jacques Leroudier, Michel Parent

► To cite this version:

Jacques Leroudier, Michel Parent. Discrete event simulation modelling of computer system for performance evaluation. [Research Report] IRIA-RR-177, IRIA. 1976, pp.69. hal-04716348

HAL Id: hal-04716348

https://inria.hal.science/hal-04716348v1

Submitted on 1 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
IRIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B.P. 5 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

<figure><img src="image_3.png" /></figure>

ARCHIVES

DISCRETE EVENT
SIMULATION MODELLING
OF
COMPUTER SYSTEM
FOR
PERFORMANCE EVALUATION

Jacques LEROUDIER
Michel PARENT

Rapport de Recherche N° 177

Juin 1976
DISCRETE EVENT SIMULATION MODELLING

OF COMPUTER SYSTEM FOR PERFORMANCE EVALUATION

Jacques LERODIER - Michel PARENT

**ABSTRACT :**

Performance evaluation through computer simulation consists of three separate phases once one has decided what questions must be answered. First of all, a model of the system, that is a simplification of its structure and its operation in time, has to be described in terms of a computer program. The basic principles of this description will be presented in chapter II while the implementation languages will be reviewed in Chapter III. The second phase consists in collecting measurements during the execution of the program. This operation is very similar to the measurements of real systems and can be performed with a number of statistical techniques discussed in Chapter IV. Since "the purpose of a simulation experiment is to predict some aspect of reality" (Naylor), it is important to test this function to gain confidence in the results obtained. This problem and its various implications such as the modelling assumptions and random number generation are considered in the last chapter.

**RESUME :**

Une fois que les questions auxquelles l'on veut répondre ont été clairement déterminées, l'évaluation de performance par simulation se déroule en 3 phases. Tout d'abord un modèle du système, c'est-à-dire une simplification et de sa structure et de son fonctionnement dans le temps, doit être décrit en termes de programme d'ordinateur. Les principes de base d'une telle description sont présentés au chapitre II tandis que les langages de simulation le sont au chapitre III. La deuxième phase consiste à recueillir les mesures intéressantes durant l'exécution de la simulation. Cette opération, qui est très semblable à la collecte de données sur un système réel, nécessite un nombre de techniques statistiques qui sont exposées au chapitre IV. Puisque "le but d'une simulation est de prédire quelques aspects de la réalité" (Naylor), il est important de valider toute simulation pour avoir une certaine confiance en ses résultats. Ce problème et ses implications, en particulier sur les hypothèses de modélisation ou la génération de variables aléatoires, sont considérés au dernier chapitre.
DISCRETE EVENT SIMULATION MODELLING
---
OF COMPUTER SYSTEM FOR PERFORMANCE EVALUATION
---

Jacques LERODIER

Michel PARENT

IRIA - LABORIA

# I - INTRODUCTION

In the past few years, computer performance evaluation has led to the development and improvement of measurement and modelling techniques. Due to the growing complexity of computer systems it has become more and more difficult to understand system behavior and to predict performances. Therefore designers and many users are now starting to use modelling as a standard design and evaluation tool.

Computer performance modelling consists in the use of a formalism to express relations between parameters and variables of interest. In this formalism, the variables are often stochastic and it is usually possible to distinguish the input variables (environment) from the output variables (performances) or equivalently the independent and the dependent variables.

Let us consider for example a multiprogramming virtual memory computer system where we wish to study the relations between the multiprogramming level, the program characteristics, the CPU utilization and the various channel activities. One formalism is to express the relations between these variables by a number of equations obtained, for instance, by regression analysis from direct measurements on an existing system. This formalism does not consider the structure of the computer. If we take this structure into account and simplify it in the queueing system of fig. 1, we can develop a formalism in terms of probabilistic equations under certain assumptions. If these equations cannot be obtained or solved due to the complexity of the model, it remains to express the operation of the system (that is the changes
in time of the input or output variables needed to describe the system),
in the form of an algorithm which can be executed on a machine.

<figure><img src="image_4.png" /><figcaption>Figure 1 : Queueing model of a multiprogrammed virtual memory computer.</figcaption></figure>

This last approach is called computer simulation if the machine
used to execute the algorithm is a digital computer and we will look at it
in some detail in this paper. Performance evaluation through computer
simulation consists of three separate phases once one has decided what
questions must be answered. First of all, a model of the system, that
is a simplification of its structure and its operation in time, has to be des-
cribed in terms of a computer program. The basic principles of this descrip-
tion will be presented in chapter II while the implementation languages will
be reviewed in chapter III. The second phase consists in collecting measure-
ments during the execution of the program. This operation is very similar to
the measurements of real systems and can be performed with a number of statis-
tical techniques discussed in chapter IV. Since "the purpose of a simulation
experiment is to predict some aspect of reality" [Nay 67], it is important to
test this function to gain confidence in the results obtained. This problem
and its various implications such as the modelling assumptions and random
number generation are considered in the last chapter.
## II - DISCRETE EVENT SIMULATION

We will now present the basic principle of discrete event simulation and the various techniques needed for its implementation on a digital computer. The approaches presented here are used in simulation languages but have to be coded by the user if he wishes to use any general purpose language.

### 2.1 - Basic principle

The basic principle of simulation with respect to time is to follow all the changes in the system modelled as time goes by. Therefore the operation of any simulator of this type consists in changing the values of the variables which define the state of the system as the simulation time is advanced.

The state of the system modelled is sometimes defined by a (fixed or variable) number of entities and the values of their attributes, but although this terminology is well accepted in simulation, it appears to be too restrictive. We define the state of the system in a more general way, as the set of variables defined by the user to describe the system at any time t.

> In discrete event simulation, the state space of the system must be finite or denumerable and changes of state, or events, occur at discrete instant of time $t_i$, i ∈ N called occurrence times of the events.

It is not possible in discrete event simulation to consider in the system state, variables which can vary continuously as it can be done with continuous time simulation [For 68]. Hence the user will have to be very careful in defining the system state so that the basic principle of discrete event simulation can be applied. Furthermore, this is only a necessary condition : the user must also be able to define the evolution of his system in terms of the events defined on the system state changes. It is however possible to use variables which may or may not be continuous if they are not included in the state space and considered only at occurrence times of events.

Since the events occur in a discrete fashion, they can be sequenced in time according to their chronological order and hence simulated one after another on a sequential computer. The parallelism which appears at the functional level of the system under study does not exist when changes of state are considered, and is thus referred to as quasi-parallelism.
We shall assume here that all the state changes are well defined even when they are the result of "composite events", that is, when several logical actions take place at the same time. This is not always true as can be seen in the cases of "race condition" where two or more processes can modify simultaneously the same element of a system (semaphore, memory cell, flip-flop,...). To remove these ill defined state changes the true parallelism must be replaced by an explicit or implicit sequencing of the elementary operations. This is always the case in simulations performed on a sequential machine, so we shall not come back to this problem (for further details see [Par 74]).

In discrete event simulation, the occurrence of an event results in a state change. Two methods are used to advance the time from one occurrence to the next occurrence time: the *event driven method* and the *clock driven method*.

### 2.1.1 - Event driven simulation

The only instances of time considered during the simulation are occurrence times, and the time variable is incremented at the occurrence of each event (the time can be expressed as a real variable which may take any value). The time management is based on a data structure which can be thought of as an ordered linear list of events called the *event list* or *time table*: at the head of the list is the current event, i.e. the event being processed and the last element of the list is the event whose occurrence time is the most remote in the future. However, all the future events in the list are only *potential*: their occurrence can be cancelled or displaced in time by any previous event. The different possible implementations and management techniques of this data structure in event driven simulation are presented in section 2.2.

It should be noticed that if two events have the same occurrence time, the order in which they will be processed depends on the way the time table is updated and scanned. Therefore, in most simulation languages, there exist special statements such as BEFORE, AFTER or PRIOR in SIMULA-67, which force the desired sequencing.
### 2.1.2 - Clock driven simulation

Time is discretized with an appropriate unit of time, and the simulation is driven by a clock which is periodically incremented by this unit of time. At each incrementation step the potential events are examined and occurrence times of events are compared to the clock time in order to find whether an event occurs at this time or not [Pee 74]. Therefore, a list of events does not have to be maintained and the implementation is simpler than in the event driven simulation case. The main drawback of the clock driven technique is in the number of times the scanning of the potential events does not produce an event time. This is measured by the *event scanning miss ratio* which depends on the unit of time selected to increment the clock time [Bab 75].

## 2.2 - The structure of a synchronization kernel

As it was pointed out above, an event driven simulation requires the implementation of an ordered list of events. The data structure and the set of programs which are used to manipulate this list are called the *synchronization kernel* or the *time-flow mechanism* and is to be found in any event driven simulator or simulation language such as SIMULA-67, SIMSCRIPT, GPSS.

### 2.2.1 - Events and activities representation

An event has been defined in section 2.1 as a change of state of the system under study : during a simulation, the occurrence of an event causes the execution of a program which actually realizes this change of state. Each execution of such a program is refered to as an *activity* of the simulation and can be implemented in the two following ways :

- a - The synchronization kernel provides elementary functions which allow one to schedule events by statements such as : "EVENT X AT time t". When the event under consideration arrives at the head of the time table, the *subroutine* associated with the event is executed. This subroutine, called the
**event routine**, describes the activity i.e. the state change of the system. The event-routine may itself reschedule the event which called it or other events. Such synchronisation kernels are said to be **event oriented**: it is for instance the case of SIMSCRIPT.

b. Each **activity** of the simulation (which remains defined as the description of a state change) and is thus associated with an **event**, is seen at a logical level as the "activity" of a process. The model then consists of a set of processes which are simultaneously active. A process represents an entity which evolves in time, as opposed to other entities of the simulation which remain passive. The parallelism between the different **activities** implies that they are described by coroutines [Knu 1] instead of **subroutines**. Hence the processes are executed in quasi-parallelism. The synchronisation kernel provides statements such as: "ACTIVATE PROCESS X AT time t", and is said to be **process oriented**: it is the case for SIMULA.

The two approaches are indeed equivalent but differ by the facilities they provide to describe a given problem. The first one is well suited for problems set in terms of events as for instance the interruptions handling in a computer system. The second approach is better adapted to problems defined in terms of parallel activities and allows a formulation of the system close to its logical description. Moreover, the **process-oriented** approach allows a duality in the description of the simulated system by allowing the active entities to become passive entities and vice versa. For example, in a computer system, programs can be thought of either as passive entities which receive services from the active entities which are the different units of the computer system (C.P.U., channel), or as active entities which serve the units of the system. This duality is important since it allows one to reduce the cost of writing the simulator and to improve the performance of the simulation by choosing the alternative which minimizes the number of processes in the simulation. Consider for instance a car
assembly line where a single car passes through 50 stations : it will be
more efficient to consider the car as the active entity which serves the
stations than the dual approach.

2.2.2 - Performance of a synchronization kernel

The performance of one simulation run is defined as the ratio of
the simulated time to the execution time. This ratio is greatly influenced
by the structure of the synchronization kernel and in particular depends
on :

- the structure of the time-table,

- the updating routines of the time-table,

- the number of events in the time-table,

- the distribution of the intervals of time between events.

Since the last two parameters depend on the environment, a robust
rather than an optimal solution must be sought for. In this section, we
present the results reported by VAUCHER and DUVAL in [Vau 75] where the
efficiency of the synchronization kernel is measured by the execution time
of the function which schedules the current event. This frequently used function
- HOLD in SIMULA-67 - takes the head of the event list and scans it
in order to find its new position according to its new occurrence time.
Several data structures can be considered for the time-table : a linear list,
an indexed linear list, a binary tree, a balanced binary tree etc... With a
linear list, the average seek time is of the order of the number *N* of events in
the list, and it is of the order of $\log_2 N$ in the case of a balanced binary tree.
For large *N*, a tree structure is thus more efficient than a linear list.
Moreover, MYHRAUG [Myr 75] reports that no significant discrepancy is
observed between a balanced tree and an unbalanced one for as long as *N* does
not exceed 400.
When a linear list is used, it can be indexed in order to avoid scanning too many elements. The list is then divided in as many sublists as indices and the sublist to be scanned is determined by the index value.

The results cited in [Vau 75] are illustrated in figures 2, 3 and 4 where the execution time of the HOLD function in SIMULA-67 is plotted as a function of the number of events in the list for different distributions of interval of time between events.

<figure><img src="image_5.png" /><figcaption>Fig. 2 : Linear list</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4 : Binary tree</figcaption></figure>
The following time distributions have been considered :

1 - negative exponential

2 - uniform

3 - bimodal

4 - discrete (taking the values 0, 1 and 2 with equal probabilities).

It appears from these results that the best structure for the time table seems to be the *indexed linear list*. However this structure can lead to some difficulties because to be efficient, all the sublists should have in the average the same number of entries. The indexing step should hence depend on the time distribution between events. VAUCHER and DUVAL [Vau 75] propose a formula to compute the indexing step from the mean time between events and they mention that further investigations should be done in order to implement an adaptive algorithm.
# III - SIMULATION LANGUAGES

We have seen that simulation modelling consists in representing the operation of the system modelled in form of a computer program. This computer program can be obtained by using one of the following programming tools :

## 3.1 - General purpose languages

At the level closest to the simulation programs in binary form, one can use an assembly language. Although this seems to be a quite painful programming approach, it has been used to program large systems such as OSS [Sea 69] or special purpose simulators like the Waterloo network simulator [Fee 74] where space and efficiency are of prime importance.

It is also possible to program in any general purpose high level language. For instance FORTRAN was used to build OSSL, a programmable simulator for computer system [Dew 72] and numerous other large or small simulations [Reh 68, Noe 72] ; PL/1 was used in one case by IBM to simulate in detail OS under LASP [Wal 72] ; MAD was selected to study various scheduling algorithms in CTSS [Sch 65] ; etc... These languages require the programming of the same simulation functions (random variable generator, clock, quasi-parallelism, list handling,...) and such a programming effort should be kept for simulations where flexibility and closeness to the application are necessary to insure execution speed and the best use of memory space. Considerations of portability are also important and may impose the selection of a standard high level language if the simulator is to be run on different machines.

## 3.2 - General purpose simulation languages

It was quickly realized that the techniques used in discrete event simulation were quite standard and hence numerous simulation languages appeared as either extensions of general purpose languages or as new languages. These languages can be used efficiently for the simulation of computer systems [Tei 66, Kiw 67, McD 70] although some are more appropriate to solve particular types of problem.
We shall simply mention some of the most commonly used
simulation languages together with their characteristics :

* SIMSCRIPT [Dim 64, Mar 64, Mar 66] is an event oriented extension of FORTRAN (it used to be a preprocessor) and the most widely available simulation language. The language provides event managing functions, and the synchronization is achieved by the event themselves. The occurrence of an event causes the call of an event routine which is a subprogram describing the changes of state caused by the event.

* GPSS [Her 64, IBM 67] is a language based on block diagrams. It was specifically developed for simulation. Hence it is not as flexible as extensions of general purpose high level languages but its ease of use by people having little programming skill makes it one of the best known simulation packages. The description of a system is done through the use of functional blocks to build a diagram. The entities, called *transactions*, are created by the GENERATE blocks and flow through the diagram according to the conditions and delays specified by the blocks.

* SIMULA 67 [Dah 66, Dah 68, Bir 69] is a true extension of ALGOL 60 (i.e. it contains ALGOL 60 as a subset, except for the I/O statements which differ in SIMULA 67). It is a powerful language allowing complex dynamic data structures through the use of *classes* and *ALGOL 60 blocks*. The language is process oriented, and the processes are programmed independently and can share data. They are activated or suspended by themselves or by other processes.

* Other simulation languages: many other general purpose simulation languages have been developed each one trying to offer new features either to facilitate programming or to allow more detailed and more complex simulations to be performed. Most of the time, a user will not decide to use a particular language on its characteristics but on its availability and on the knowledge he has
of it, on its documentation and on the amount of simulation programming
he has to perform. The following are a few additional simulation languages
which are described in the literature :

- ASPOL [McD 73] although designed for simulating computer systems
is a general process oriented simulation language,

- CSL [IBM 66] was first designed as source language to a FORTRAN
translator (like SIMSCRIPT) but is now compiled directly into
machine language and is available on several machines,

- DESPL/1 [Red 73] and SIMPL/1 [Ret 74] are preprocessors producing
PL/1 code and are available on any machine with a PL/1 compiler
(they are written in PL/1),

- GASP [Kiv 64, Prk 69] is an extension of FORTRAN implemented by a
set of FORTRAN written subroutines which makes it easily available
on almost any machine. GASP has the same programming philosophy as
GPSS through the use of GASP oriented flow-charts,

- SOL [Knu 64] is an ALGOL like language which has many of the
ideas implemented later in SIMULA 67.

3.3 - *Specialized languages for the simulation of computer systems*

Although the simulation of computer systems becomes much simpler
through the use of general purpose simulation languages, the similarity of
such systems made it attractive to develop highly flexible simulation models
or specialized simulation languages. The programming investment to develop
such systems was found by several manufacturers and research centers to be
justified when the number or the complexity of simulations become large
enough. The two possible approaches are to extend an existing simulation
language by adding several functions to represent the standard parts of a
computer system or to build a simulation model general enough to model a
wide range of architectures and environments. The following is a list of a
few such systems using one or the other of the two approaches :
- CSS (Computer System Simulator) [Sea 69, IBM 70] developed by IBM uses many of the techniques of GPSS but is a completely new simulation language designed to model computer systems. However it has been implemented with the idea of simulating specifically IBM 360/370 hardware and software [Ber 74] and hence its generality to model any system is not established. Furthermore this is an IBM proprietary program.

- ECSS (Extendable Computer System Simulator) [Nie 70] written at the RAND Corporation takes a different standpoint. It is an extension of SIMSCRIPT II (the prototype version is a translator) and hence has all the flexibility of a general purpose simulation language while providing facilities to describe the system configuration and the job load. Some tools are also provided to describe the software but not to the same extent as for the hardware.

- OSSL (Operating System Simulation Language) [Dew 71, Dew 72, Wya 72] developed in universities for teaching and research on operating systems is a language in which a system can be described by using the same terminology as system designers. The system description in OSSL consists of the hardware characteristics, the operating system and the environment and is used as an input to a simulator. Both the simulator and the language interface have been coded in FORTRAN IV and are hence easily available.

- QAL (Queueing Analysis Language) [Fos 74] is a "very high level language" to describe networks of queues in order to simulate them. Although it is not specifically designed to simulate computer systems, its convenience of use and the number of computer systems which can be modelled at a global level by a network of queues when parallelism and congestion are taken into account make it very attractive. Another interesting feature is the possibility to use a QAL subprogram in any FORTRAN program.
- Other packages such as SCERT [SCE 67] or S3 [S3 67] can simulate computer systems as long as they fit into a quite restricted framework and hence lack a desired flexibility to simulate new architectures.

## 3.4 - Graphical descriptions

It has always been convenient to represent simple or complex systems by the way of a graph. This has been done in a more or less informal way until several approaches have been formalised to provide a more powerful tool. One such approach is in the logical representation of computer systems through the use of Petri nets [Pet 62], Parallel Program Schemata [Kar 69], Flow Graph Schemata [Slu 68] which have seen an extension and an implementation as a graphical tool on a computer system with the LOGOS project [Ros 72]. Such systems can be further extended [Par 74] to provide integrated tools for the design of future system where both the logical consistency and performances can be evaluated throughout the design stage [Bae 73]. Petri nets have already been extended [Noe 73] to provide for better representation and evaluation [Nut 73, Nut 72] of computer system and a simulator using as input not a program but a graph (Macro Evaluation Net) has been developed at the University of Colorado [Nut 74]. Such approaches can be compared to the graphical representation of systems in GPSS or the queueing networks described with QAL but the programming burden is greatly facilitated by graphical input on display units. The possibility of having an integrated design tool may be also an enormous advantage for computer designers.

## 3.5 - Performance of simulation runs

The usual criticism made of simulation language and particularly of specialized packages is the poor performance of simulation runs when they are compared to the run times of the same models coded in assembly or high level languages. This criticism is not always justified [Kos 73] if we consider the cost in development time (ease of modelling, coding, debugging modifications...) and even in the worst case of languages producing interpretable code, the ratio rarely exceeds twenty if both simulators describe basically the same events and collect the same information. The run
time ratio of the actual computer time spend to the simulated interval of time does not depend so much on the language used but on the level of detail of the simulation, that is the average time between occurrences of events and the complexity of these events. The table presented by Teorey [Tec 73] can be cited as reference :

<table><caption>Table 1: Run time ratios in computer simulation</caption><thead><tr><th>Type</th><th>Average Occurrence Between Events</th><th>Run Time Ratios</th><th>Application</th></tr></thead><tbody><tr><td>1. Gross functional description</td><td>50-500 milliseconds</td><td>10-50 times faster than real time</td><td>Reconfiguration</td></tr><tr><td>2. Detailed functional description</td><td>1-10 milliseconds</td><td>20 times slower to 10 times faster than real time</td><td>System design, Computer Selection</td></tr><tr><td>3. Machine instruction execution</td><td>1-10 microseconds</td><td>More than 20 times slower than real time</td><td>Subsystem Analysis</td></tr></tbody></table>

This shows clearly that it is inconceivable to simulate with the same model events that occurs at the micro-second level with events that occur at the second or minute level due to run time constraints. For example it is unfeasible using classical digital simulation to model a complete memory hierarchy going from the cache level which operates at the nano-second down to the tapes or library level which run in a range $10^{-9}$ times superior.
# IV - MEASUREMENTS AND STATISTICAL ANALYSIS

## 4.1 - Measurement problems

A simulation allows the generation of a series of events which conform to the operation -on a more or less idealized level- of a real system. Sometimes the interest of the simulation lies only in the production of this series of events and no quantitative processing is needed : this is the case when one wants to study the possibility of occurrence of a special event (deadlock, fatal transmission error,...) or when one wants to prove the functional correctness of a system (such as the hardware simulators [Mer 73, Hay 69]).

In this section we shall be concerned with simulations performed to obtain measurements and to evaluate the performance of the system under study on a quantitative basis. In order to do this, we must compute performance criterion estimates (utilization ratios, page fault rate,...) from measurements obtained during the simulation runs.

The problems brought by the collecting and the analysis of these measurements are often similar in real system performance evaluation although the technology or the method used to obtain the samples are different [Bor 71, Hol 71, Ier 73]. A major difference exists in the planning of different experiments which can be much easier in simulation.

In simulation or in real system measurements, data are always collected through one of the following ways :

- change of state of the system under study,

- human intervention (seldom in simulation),

- clock signal.

This collecting can cause several technical problems :

- the variable to be measured is not always accessible or is not updated [Boi 71],

- the reading of the variable may affect its value (during the reading period or after) [Hol 71],

- the time scale between sampling points is too large (fast variations may not be observed) [Gau 73].

and statistical problems.
In order to analyze these measurements, the series of events is considered as a stochastic process which we want to study. Any performance is hence taken as moments (or functions of moments) of associated processes [Bar 71, Rod 72]. Therefore, the collecting of measurements is nothing but the sampling of a process [Boi 73, Knu 73, Ler 75].

From this viewpoint, the following problems associated with stochastic process statistics arise :

- stationarity (or non-stationarity) test for a sampled process [Cox 68, Cox 72, Pry 69], or evaluation of a transient phase,

- influence of the process auto-correlation [Cra 74, Fis 65, Fis 67],

- search for special events which may simplify the analysis of the sample -the knowledge of renewal points* for instance may be extremely valuable-,

- in the case where such event do not exist : what is the influence of the initial state on the time spent by the process in a "transient state" before reaching a "stable state" [Kle 75],

- simulation time in order to obtain a given confidence interval [Fer 75, Fis 67]...

From the reading of many papers on simulation or on real measurements, it appears that many of these problems are often not considered. However the rigor of a sound statistical approach is not all unfounded. The lack of homogeneity of different samples of a same variable that is often observed can lead to erroneous conclusions unless great care is taken.

Data collection and analysis operated in less than satisfactory conditions can often be traced back to a lack of cooperation between :

- operating system specialists who implement software or hardware probes in order to collect real system measurements (for calibration and for trace extraction),

<footnote>*For a stochastic process X(t), $t_0$ is a renewal point if the knowledge of X(t) for $t \le t_0$ does not bring more information about the behaviour of X(t) for $t > t_0$ than the knowledge of $X(t_0)$.</footnote>
- mathematicians and statisticians who develop statistical methods to analyze data and to test hypothesis such as stationarity,

- and simulation language specialists who are too much concerned with modelling and language problems to think about the quality and the validity of the results obtained.

It seems highly desirable to work on these problems in a global way through a better communication between these specialists. In particular it is necessary to implement sophisticated collection and analysis tools in simulators and in real systems to process data. Such tools are specially needed in feed-back systems where it is imperative to detect a significant change in the operating environment [Wul 69, Bad 74], which implies the use of efficient statistical methods.

In the next sections, we shall formalize the notion of estimate and we shall give the formulas to obtain usual estimations and their confidence intervals and the hypotheses to do so.

## 4.2 - Characterization of the underlying stochastic processes

We consider a simulation model from which a series of measurements of one or several variables are extracted. The variations in time of any of the measured variables can be mathematically represented by a stochastic process defined as follows :

$$ \begin{array}{ccc} T \times \Omega & \longrightarrow & R \\ (t, \omega) & \longrightarrow & X(t, \omega) \end{array} $$

where $t$ belongs to a time interval $T$ and $\omega$ is one outcome on the probability space $\Omega$.

Hence for a given $\omega$, $X(t, \omega)$, function of $t$, is one realization of the stochastic process and is called a **trajectory** or **sample path**. This function may be defined for any $t$ (continuous parameter process) or only at specific instants such as throughput which is defined only at service ends. We should also mention point processes where only the occurrence time of events are of interest.

Let us suppose for the moment that $X(t, \omega)$ is defined for any $t$ and any $\omega$. The process can be sampled basically in two different ways :
1. a sampling can be performed on the probability space $\Omega$, that is for a given $t$, we collect a sample of the random variable $X(t)$. This method can only be used through the replication of several runs with different seeds for the random number generator (cf : § 5.2). The sample obtained will be uncorrelated if the random number generator is "good" (see 5.2) and statistical analysis for independent samples can be applied for any $t$.

2. or the sampling can be performed on the time axis for a given trajectory ($\omega$ fixed). This method is the only one possible with real systems and is the most commonly used in simulation where the replication of runs is too cumbersome and also expensive in computer time. In this case, we observe at discrete sampling times $t_i$ ($i = 1, 2, ..., N$) with $t_i > t_{i-1}$ one trajectory $X(t, \omega)$ for a given $\omega$ and since the sample may be autocorrelated and non stationary, this leads to several technical difficulties.

<figure><img src="image_7.png" /><figcaption>Fig. 5: Stochastic process sample paths</figcaption></figure>
We shall now concentrate on these difficulties and the methods used to obtain estimates and confidence intervals in the case where we observed one single trajectory of a continuous parameter random process X(t) (of a real or a simulated system).

The stochastic process would be perfectly known if, for any set $t_1, t_2, \dots, t_m$, and any finite m, the joint distribution could be obtained :

$$F(t_1, t_2, \dots, t_m; x_1, x_2, \dots, x_m) = \text{Prob}[X(t_1) \le x_1, \dots; X(t_m) \le x_m]$$

It is obviously impracticable to estimate this function for any m. Most often only the marginal distribution function

$$F(t; x) = \text{Prob}[X(t) \le x]$$

is sought for and indeed in many cases only estimates of the mean $E[X(t)] = \mu(t)$ and the variance $E[(X(t) - \mu(t))^2] = \sigma^2(t)$ are computed.

If more information is desired about the behavior of the stochastic process the correlation function between the two random variables X(t) and X(t+τ) :

$$\rho(t, \tau) = \frac{E[(X(t) - \mu(t))(X(t+\tau) - \mu(t+\tau))]}{\sigma(t) \sigma(t+\tau)} \quad \text{where } -\infty < \tau < +\infty$$

can be estimated. However, due to practical and phenomenological considerations, this function is estimated when the process is supposed to be weakly stationary, that is, when we have for any t :

$$\left. \begin{array}{l} \mu(t) = \mu \\ \text{Cov}(X(t), X(t+\tau)) = \gamma(\tau) \end{array} \right\} \quad \forall t$$

which implies

$$\rho(t, \tau) = \frac{\gamma(\tau)}{\gamma(0)} = \rho(\tau) \quad \forall t$$
Here $\rho(t)$ is called the _autocorrelation function_ which can now be
estimated through a single simulation run if the process is ergodic (4.2.1)

Under the same conditions, it is known that the spectral density :

$$
f(\lambda) = \frac{2}{\pi} \int_{0}^{\infty} \rho(\tau) \cos \lambda\tau \, d\tau \qquad 0 \le \lambda \le \pi
$$

which is the Fourier cosine transform of the autocorrelation function
[Jen 61, Par 61, B1o 76] can also describe the stochastic process and give
useful information in the study of the variations of X(t) around its mean.
It is also useful in the detection of cyclic trends [Fis 65, B1o 76].

We will see that the problems are quite different in complexity
according to the nature of the stochastic process X(t).

4.2.1. - X is a strictly stationary process

This means that for any $t_1, t_2, \dots, t_m$, the joint distribution function of the random variables

$$
X(t_1), X(t_2), \dots, X(t_m),
$$

is the same as the joint distribution function of

$X(t_1+h), X(t_2+h), \dots, X(t_m+h)$, for any $h$ that is

$$
\begin{array}{ll}
F_{X(t_1), \dots, X(t_m)} = F_{X(t_1+h), \dots, X(t_m+h)} & \forall h, t_1, \dots, t_m \\
& \forall m = 1, 2, \dots
\end{array}
$$

This implies in particular (by setting m=1) that the distribution
function F(t,u) of X(t) does not depend on t. Hence this distribution function
can be estimated in one simulation run, for example by drawing the histogram
or better, the empirical cumulative distribution function [Cox 68]
of the measured series X₁, X₂, ..., Xₙ where Xᵢ = X(tᵢ,ω) is the measured
value at time tᵢ of one realization (ω) of the stochastic process.
From the definition of a stationary process, we have :

$$ (1) \quad E[X(t)] = \mu \quad \forall t $$

$$ (2) \quad \rho(t, \tau) = \frac{\gamma(\tau)}{\gamma(0)} = \rho(\tau) \quad \forall t, \forall \tau $$

$$ (3) \quad \mathrm{Var}[X(t)] = \sigma^2 = \gamma(0) \quad \forall t $$

The estimators of the moments will be different if we can observe the complete trajectory $X(t, \omega)$ (for instance if $X(t)$ changes in a discrete fashion at times $t_i$), or if we sample it at arbitrary times $t_i$.

Since the process is stationary, we have the following relations* :

$$ (4) \quad E\left[\frac{1}{T} \int_{0}^{T} X^n(t) dt\right] = E[X^n] \quad \forall T $$

$$ (5) \quad E\left[\frac{1}{T} \int_{0}^{T} X(t) \cdot X(t+\tau) dt\right] = E[X(t) \cdot X(t+\tau)] = \gamma(\tau) + \mu^2 \quad \forall T $$

Hence, if it is possible to record a complete trajectory of the process $X(t)$, convergent estimates of any moment or of the autocorrelation function can be obtained by taking the time averages which intuitively collect the "maximum amount of information" from the run.

For example, if the value of $X(t)$ changes in a discrete fashion at times $t_1, t_2, \dots, t_N$, then a time estimate of an n-order moment can be obtained by :

$$ (6) \quad \hat{E}(X^n) = \frac{1}{t_N - t_0} \sum_{i=1}^{N} (t_i - t_{i-1}) X_i^n $$

if $t_N - t_0$ is a fixed time interval and not a random variable.

When an outcome of $X(t)$ cannot be known for all $t$ ($t_0 < t < t_N$), either because $X(t)$ is not discrete or because the time $t_i$ of variation of the outcome are not available, it is necessary to sample the outcome at arbitrary instants $t_i$. The method generally used consists in taking regular time intervals :

$$ \Delta t = t_i - t_{i-1} = \text{cst} \quad \forall i $$

<footnote>*More generally if the expressions themselves (and not their expected value) converge when $T \to \infty$, the process is said to be ergodic with respect to the corresponding moment or auto-covariance function.</footnote>
The following estimates for the mean and the variance are based on the standard *unbiased estimators* for independent samples (an estimator $\tilde{Y}$ of parameter Y is unbiased if $E(\tilde{Y}) = Y$):

$$ (8) \quad \hat{E}(X(t)) = \hat{\mu} = \frac{1}{N} \sum_{i=1}^{N} X_i, $$

$$ (9) \quad \hat{\mathrm{Var}}(X(t)) = \hat{\sigma}^2 = \frac{1}{N-1} \left[ \sum_{i=1}^{N} X_i^2 - \frac{1}{N} \left( \sum_{i=1}^{N} X_i \right)^2 \right]. $$

If $\Delta t$ is taken as the unit of time, the autocorrelation of lag k (which is time independent because of the stationarity) can be estimated by:

$$ (10) \quad \hat{\rho}(k) = \frac{\hat{\gamma}(k)}{\hat{\gamma}(0)} $$

where

$$ (11) \quad \hat{\gamma}(k) = \frac{1}{N-k} \sum_{i=1}^{N-k} X_i X_{i+k} - \frac{1}{(N-k)^2} \left( \sum_{i=1}^{N-k} X_i \sum_{i=1}^{N-k} X_{i+k} \right) $$

for $k = 0, 1, \dots, K$ are the estimates of $\gamma(k)$. K must be chosen so that $K < N$, because the sample size to estimate $\gamma(k)$ decreases as K increases. Experiments reported in [Box 70] have shown that satisfactory results were obtained up to $K = N/4$ for $N = 40$.

We can notice that $\gamma(0) = \mathrm{Var}(X(t))$.

A "smoothed estimate" of the spectral density can also be obtained with the following formula:

$$ (12) \quad \hat{f}(\lambda) = \frac{1}{\Pi} \left[ \hat{\gamma}(0) + 2 \sum_{k=1}^{N-1} W_k \hat{\gamma}(k) \cos \lambda k \right] $$

where the $W_k$ are suitably chosen weights (for example $W_k = 1 - \frac{k}{N}$ as reported in [Box 70]).

<footnote>* We have adopted the following notation throughout this chapter: if $Y$ is a parameter to be estimated, $\tilde{Y}$ is an *estimator* of this quantity (i.e. a random variable) and $\hat{Y}$ is an *estimate* of the same quantity (i.e. an outcome of the estimator $\tilde{Y}$).</footnote>
If the time intervals between the sampling points are not equal and depend on the stochastic process sampled, the estimates can be arbitrarily biased (consider the case where a queue length would be sampled each time the server becomes idle!). However, even when the intervals are equal, the estimates of the autocorrelation function or the spectral density can be misleading because of aliasing [B10 76].

### 4.2.2 - X is a weakly stationary process

Unfortunately without knowing theoretical properties of the model, it is not possible to know if a stochastic process is strictly stationary or not. However one can test [Prs 69] for the following two conditions :

- the expected value of $X(t)$ is constant :

$$ (13) \qquad E[X(t)] = \mu \quad \forall t $$

- the auto-covariance of $X(t)$ does not depend on $t$ but only on the considered interval :

$$ (14) \qquad \mathrm{Cov}[X(t), X(t+\tau)] = \gamma(\tau) \quad \forall t $$

A process satisfying these two conditions is said to be weakly stationary. Notice that only the two first moments of $X(t)$ have to be constant in time, since from (14), we have :

$$ (15) \qquad \mathrm{Var}(X(t)) = \gamma(0) = \sigma^2 $$

It is clear (since it is characterized by its two first moments) that a Gaussian process is _strictly stationary_, if and only if it is weakly stationary, that is in this case, the two definitions are equivalent [Cox 72].

Using the methods previously described in 4.2.1., it is possible to estimate only the first two moments, the autocorrelation function and the density spectrum of a _weakly stationary_ process.
### 4.2.3 - X is not a stationary process

In this case all the information we may want to gather on the stochastic process, in particular all its moments, will be time dependent. Hence they can be estimated for each time $t_1, t_2, ..., t_i, ...$ by performing several simulations (or measurement runs). The use of different seeds for the random number generators will insure the independence of the outcomes. Hence if we have N outcomes $X_1(t), ..., X_N(t)$ of $X(t)$ at time t, it will be possible to estimate all the moments of $X(t)$ at this time :

$$ (16) \qquad \hat{E}[X(t)^n] = \frac{1}{N} \sum_{i=1}^{N} (x_i(t))^n $$

If we also have N outcomes $X_1(t'), X_2(t'), ..., X_N(t')$ at time t', it is possible to estimate the covariance between the two random variables $X(t)$ and $X(t')$ by :

$$ (17) \qquad \hat{\text{Cov}}(X(t), X(t')) = \frac{1}{N} \left[ \sum_{i=1}^{N} X_i(t) X_i(t') - \frac{1}{N} \sum_{i=1}^{N} X_i(t) \sum_{i=1}^{N} X_i(t') \right] $$

It is always very costly to study a non-stationary process through independent runs and, in some applications such as real time evaluation, the method cannot be used.

There are two main approaches to the reduction of non-stationary process to a stationary one : one tries either to decompose the process into a deterministic trend plus a stationary stochastic process [Cox 66] or to divide the time axis into slices during which the process may be considered as stationary. The second method is the most commonly used in computer simulation modelling and is very fruitful. It means that we may distinguish two time scales. For instance, in a model of computer system, we may have :

- a microscopic scale at the level of events (CPU time service between the page faults or the I/O)

- a macroscopic scale at the level of the time variations of probability distributions (changes in the user or program behaviors)

Without any theoretical assumption, it is very difficult to find the appropriate time division. One way might be to consider an arbitrary time slicing and use a non parametric recursive test to find such intervals where the process has a given probability (confidence level) to be stationary, but
such a test has not yet been developed.

In some cases the appropriate time slicing can be justified on theoretical grounds when it is known that the characteristics of the stochastic process change at particular instant (for example, a number of customers or programs in the system). If these instants are distant enough, it is possible to estimate the process which is stationary during each time slice. Systems having such a property are said to be quasi-decomposable and have been extensively studied in the case of Markov chains (variable aggregation method). Applications of this theory may be found in economics [Sim 61, And 63] and in computer systems [Cou 71, Avi 73].

### 4.3 - Confidence level of estimates

Once we have obtained an estimate of a performance variable (usually a mean value), it is important to have an indication of its "quality". This is usually done by measuring the variability of the estimator for example by comparing the estimate with the highest and the lowest values of the sample or by computing its variance. This last criterion is used in particular to compare two estimators $\tilde{X}$ and $\tilde{Y}$ of a same quantity (e.g. a time average and an arithmetic average) by defining the relative efficiency :

$$Q = \frac{\text{Var}(\tilde{X})}{\text{Var}(\tilde{Y})}$$

The variance of the estimator may also be used to compute a confidence interval around the estimate but one has to remember that this is simply one way to express the variability of the measure.

Let us consider in detail the problem of finding a confidence interval around the arithmetic average of n outcomes $X_1, \dots, X_n$ of a stochastic process $X(t)$ and the variance of this estimator of the mean value of $X(t)$ :
$$ (18) \quad \hat{E}(X(t)) = \hat{\mu} = \frac{1}{N} \sum_{i=1}^{N} X_i $$

which is an outcome of the random variable (the estimator) :

$$ (19) \quad \tilde{\mu} = \frac{1}{N} \sum_{i=1}^{N} X(t_i) $$

(recall that $X_i$ is one outcome of the random variable $X(t_i)$).

In this case, $\hat{E}(X(t)) = \hat{\mu}$ is an *unbiased estimate* of the expected value $E(X(t)) = \mu$ of $X(t)$, because the expected value of the estimator $\tilde{\mu}$ is equal to the parameter estimated :

$$ (20) \quad E(\tilde{\mu}) = \frac{1}{N} \sum_{i=1}^{N} E(X(t_i)) = \frac{1}{N} \sum_{i=1}^{N} \mu = \mu $$

The problem is to find a *confidence interval* [$\mu - \epsilon, \mu + \epsilon$] that is, an interval such that the estimator $\tilde{\mu}$ of $\mu$ has a given probability $\delta$ (or *confidence level*) to lie within :

$$ (21) \quad \text{Prob}[|\tilde{\mu} - \mu| < \epsilon] \ge \delta $$

The first method to compute a confidence interval is to use Tchebycheff inequality which requires only the weak stationarity of $\chi^2$ :

$$ (22) \quad \text{Prob}[|\tilde{\mu} - \mu| < \epsilon] \ge 1 - \frac{\sigma^2}{\epsilon^2} $$

where $\sigma^2 = \text{Var}[\tilde{\mu}] = E[(\tilde{\mu} - \mu)^2]$ if $\tilde{\mu}$ is an unbiased estimator ($E(\tilde{\mu}) = \mu$).

Since the inequality can be rewritten :

$$ (23) \quad \text{Prob}[|\tilde{\mu} - \mu| < \alpha\sigma] \ge 1 - \frac{1}{\alpha^2} $$

to obtain a 90 % confidence interval, this implies that

$$ \alpha = \sqrt{10} \approx 3.16 $$

and hence gives a confidence interval of 3.16 standard deviation of $\tilde{\mu}$ each direction around the mean.
Better confidence intervals can be obtained, if $\tilde{\mu}$ can be supposed
Normally distributed. This is indeed the case if $\tilde{\mu}$ is obtained through the
sum of independent random variables identically distributed (central limit
theorem). Then the formula giving the confidence interval is :

$$
(24) \quad \text{Prob}[|\tilde{\mu} - \mu| < \alpha\sigma] = \frac{1}{\sqrt{2\pi}} \int_{-\alpha}^{+\alpha} e^{-\frac{x^2}{2}} dx = 2\Phi(\alpha) - 1
$$

where $\Phi(\alpha) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\alpha} e^{-\frac{x^2}{2}} dx$ and is given by Normal distribution tables.

Typical values for $\alpha$ are 1.65 and 1.96 which give respectively
90 % and 95 % confidence levels.

It is easily seen that the lenght of the confidence interval is
thus considerably reduced compared to the previous method since for a confidence
level of 90 % it falls from 6.3 σ down to 3.3 σ. The underlying assumption
to have this reduction is the strict stationarity of X and the independence
of the sample used to estimate E(X(t)) = μ.

In all cases, to obtain a confidence interval, we are led to compute
the variance Var($\tilde{\mu}$) of the estimator $\tilde{\mu}$ (this variance should not be
confused with the variance of the stochastic process X(t)).

Anyway, this variance by itself is a good indication of the
variability of the estimator and it may not be necessary to compute a
confidence interval specially if the estimator is non Normally distributed.
The variance depends very much on the estimator selected but it will be
interesting to try to have a convergent estimator :

$$
(25) \quad \lim_{N \to \infty} \mathrm{Var}(\tilde{\mu}) = 0
$$

to have an arbitrarily small variance (and confidence interval).. In this case, we
can reach a given precision by increasing the simulation lenght. The
estimation of Var($\tilde{\mu}$) is critical because it is sensitive to the sample auto-
covariance :

$$
(26) \quad \mathrm{Var}(\tilde{\mu}) = \mathrm{Var}\left(\frac{1}{N} \sum_{i=1}^{N} X(t_i)\right) = \frac{1}{N^2} \left[ \sum_{i=1}^{N} \mathrm{Var}(X(t_i)) + \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \mathrm{Cov}_X(X(t_i), X(t_j)) \right]
$$
If X is weakly stationary, we have :

$$
(27) \quad \operatorname{Var}(\bar{x}) = \frac{\operatorname{Var}(X(t))}{N} + \frac{1}{N^2} \sum_{\substack{i,j=1 \\ i \neq j}}^{N} \operatorname{Cov}(X(t_i), X(t_j))
$$

If the variables X(t₁),...,X(tN) are independent we have ∀tᵢ, tⱼ and
tᵢ ≠ tⱼ Cov(X(tᵢ), X(tⱼ)) = 0, thus

$$
(28) \quad \mathrm{Var}(\bar{x}) = \frac{\mathrm{Var}(X(t))}{N}
$$

where Var(X(t)) is estimated by (9). We may notice that (28) is convergent.

If the random variables X(tᵢ) are not independent, from the assumption
of weak stationarity, we have :

$$
\operatorname{Cov}(X(t_i), X(t_j)) = \gamma(|t_i - t_j|).
$$

Since the estimation of the last terms in (27) is a difficult and
time consuming task, the auto-covariance γ can be estimated from the sample
by using (11) and fitted by an arbitrary function [Cox 68, Bro 63] to be
used in other simulation runs of the same model. The function chosen is
often linear or exponential.

The method we suggested takes into account the possible autocorrelation
of the sample, now we are going to study methods which attempt to eliminate
it to estimate Var(μ).

4.3.1 - Using blocked measurements [Con 63]

The sample X₁, X₂, ..., XN is divided into p blocks of equal length l :

<figure><img src="image_8.png" /></figure>

and a partial average is computed on each block j :
$$ (29) \qquad \hat{E}(X(t))^{(j)} = \hat{\mu}^{(j)} = \frac{1}{\ell} \sum_{i=(j-1)\ell}^{j\ell} X_i $$

The idea is to have a block partition such that $(\hat{\mu}^{(1)}, \hat{\mu}^{(2)}, \dots, \hat{\mu}^{(p)})$ is an uncorrelated sample of the stochastic variable $\tilde{X}$. Then the variance of $\tilde{\mu}$ is estimated by :

$$ (30) \qquad \mathrm{Var}(\tilde{\mu}) = \frac{1}{p-1} \left[ \sum_{i=1}^{p} (\hat{\mu}^{(i)}) - \frac{1}{p} \left( \sum_{i=1}^{p} (\hat{\mu}^{(i)})^2 \right) \right] $$

We may notice that after having determined a maximum block size $\ell_0$, this method can be *iterative*, that is, used to estimate in real time $\mathrm{Var}(\tilde{\mu})$ in order to stop a simulation run when reaching a given precision.

The approach used in a simulation will consist of two steps :

1. determining a minimum block size $\ell_0$ from a first trajectory of $X(t)$.

2. during all the other simulation runs the block size $\ell > \ell_0$ is given as a parameter.

The minimum block size $\ell_0$ must be chosen to be consistent from a simulation run to another one, it will be determined from the autocorrelation function in such a way that :

$$ (31) \qquad |\rho(\ell_0)| = \left| \frac{\gamma(\ell_0)}{\gamma(0)} \right| \ll 1 $$

Theoretically, we should choose : $\ell > \ell_0$ and consider only the blocks separated by a distance of $\ell_0$. In practice, it is sufficient to choose $\ell_0 < \ell \le 2\ell_0$ and consider all the blocks, but to be very careful, $\ell_0$ must be re-estimated periodically when different experiments are performed with the simulation model.

### 4.3.2 - Using spectral analysis [Fis 65, Fis 67]

In fact, this method differs from the previous one (cf. : 4.3.1) by a different approach used to estimate $\ell_0$. The spectral density function
(cf. : 4.2) is estimated by (12) and is used to define the correlation time of the stochastic process [Fis 67] *:

$$ (32) \qquad \tau^* = \frac{\Pi}{2} f(0) $$

This correlation time corresponds to the ideal case of a stochastic process which is entirely autocorrelated on the distance $\tau^*$ and not correlated beyond. Fishman defines an equivalent number of independent measurements which is :

$$ (33) \qquad N^* = \frac{N}{2\tau^*} $$

This means that two measurements which are separated by a distance of $2\tau^*$ can be considered as independent. Fishman considers only such measurements in his statistical analysis which leads to a great loss of information. We suggest to use the previous approach (4.3.1.) with :

$$ (34) \qquad \ell_0 = 2\tau^* $$

### 4.3.3 - Using regeneration points [Cra 74, Cra 75, Lav 75]

We now present a method which solves both the problem of autocorrelation and the problem of the initial state (cf. : 4.5), it is the technique of _regenerative simulation_ using _regeneration points_ [Cra 74, Cra 75, Lav 75]. Informally the process $X(t)$ is said to be regenerative, if when the system starts from a given state (called regeneration state) at time $t = 0$, there exists with probability one a time $t \neq 0$ at which the system will be found in exactly the same state. The times at which the system is found in the regeneration state are called regeneration points. Hence we can decompose a simulation run into a number of independent and statistically identical subruns, each subrun starting from a regeneration point. The _regeneration state_ is chosen as the _initial state_ of the simulation run.

<footnote>*However, the estimate of $f(0)$ obtained this way is very sensitive and it is better to estimate it directly from the data ($f(0) \approx [\text{Var}(X_1 + \dots + X_n)/n]$)</footnote>

<footnote>**Including past history</footnote>
If the time intervals between regeneration points have a finite expectation, it is possible to obtain any given number of subruns (also called *cycles* or *tours*) if we allow the simulation to run for a sufficiently long time.

This method can be used each time the process is known to have a regenerative behavior and when the mean time between points is finite.

These conditions are satisfied, for example, in the case of ergodic Markov chains (any state can be taken indifferently as a regeneration state) or for certain open queueing systems (e.g. : an M/G/1 queue) if the traffic intensity is less than one and if the service times have finite fourth moments [Cra 75a]. However, computer systems are often modelled as a complex network of queues.

In an open system with a Poisson arrival process, each instant where a customer arrives when the system is empty is a regeneration point. But it is in general impossible to prove that the time between such points has a finite expectation. So, in order to use the method, this has to be taken as conjecture [Lav 75a]. Approximate regeneration states can also be used in order to guarantee the finiteness of the cycle lenght [Cra 75b].

Since all cycles are independent and stochastically identical, it is possible to use standard statistics to compute confidence intervals on the estimates obtained from the sample.

However one must be careful to take into account the duration of each cycle since it is a random variable correlated with the sample. For this, the following theoretical results have to be used (for a detailed description see [Cra 75a]):

Let $t_0, t_1, \dots, t_k, \dots$ be the regeneration points of a stochastic process $X(t)$ with general state space. Let $\varphi[X(t)]$ a measurable function on the state space such that $E\{|\varphi(X)|\} < \infty$. The function $\varphi$ can be selected to
produce the moments, density or autocorrelation function etc...

Let $Y_k$ be the time average of $\varphi(X)$ on the k-th cycle (which has a random length $T_k = t_k - t_{k-1}$), that is:

$$ (34) \quad Y_k = \int_{t_{k-1}}^{t_k} \varphi(X(t)) \, dt \quad . \text{Then we have :} $$

$$ (35) \quad E[\varphi(X)] = \frac{E[Y_k]}{E[T_k]} \quad \forall k = 1, 2, \dots $$

Hence if $\hat{Y}_k$ and $\hat{T}_k$ are estimates of $E[Y_k]$ and $E[T_k]$ obtained from the sample on the k-th cycle and if :

$$ (36) \quad \hat{Y}(n) = \sum_{k=1}^{n} \hat{Y}_k $$

$$ (37) \quad \hat{T}(n) = \sum_{k=1}^{n} \hat{T}_k $$

we can apply the central limit theorem where the variance is replaced by an estimate [Chu 68] :

$$ (38) \quad \lim_{n \to \infty} \text{Prob} \left\{ \left| E[\varphi(X)] - \frac{\tilde{Y}(n)}{\tilde{T}(n)} \right| < \frac{\sqrt{n}}{\hat{T}(n)} \right\} = 2 \Phi(t) - 1 $$

where $\Phi(t) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{t} e^{-\frac{x^2}{2}} dx$ and $\tilde{Y}(n), \tilde{T}(n)$ are the estimators whose outcomes are $\hat{Y}(n)$ and $\hat{T}(n)$ and $\hat{\sigma}^2(n)$ is an estimate of the variance of $Y_k - E[\varphi(X)]$. $T_k$ which is obtained by :

$$ (39) \quad \hat{\sigma}^2(n) = v_1(n) - 2 \frac{\hat{Y}(n)}{\hat{T}(n)} v_2(n) + \left[ \frac{\hat{Y}(n)}{\hat{T}(n)} \right]^2 v_2(n) $$

with

$$ (40) \quad v_1(n) = \sum_{k=1}^{n} \frac{(\hat{Y}_k - \hat{Y}(n))^2}{n-1} $$
$$ (41) \qquad V_2(n) = \sum_{k=1}^{n} \frac{(\hat{T}_k - \hat{T}(n))^2}{n-1} $$

$$ (42) \qquad V_{12}(n) = \sum_{k=1}^{n} \frac{(\hat{Y}_k - \hat{Y}(n))(\hat{T}_k - \hat{T}(n))}{n-1} $$

An example where the method is used to estimate average waiting time in a queueing system is presented in [Lav 75a-b].

## 4.4 - Variance reduction techniques

We have seen in the previous paragraph how to compute a confidence interval on an estimate obtained from a simulation run. Such confidence intervals are often in the range of 10 % or more of the estimate for a confidence level of 95 % even with long simulation runs (500 or 1000 cycles) [Lav 75b, Pot 75]. Clearly the width of the confidence interval depends on the stochastic process sampled, on the size of the sample but also on the estimator selected. In [Lav 75c], the computation of confidence intervals for different work rate estimators in closed queueing network shows large differences.

Once the stochastic process and the estimator are given, the only way to reduce the confidence interval is through the sample length [Fis 68].

For this two methods may be used : one can execute several simulation run or increase the length of a single run.

The first method is slightly more difficult to implement (two or more samples are collected), and sometimes even impossible to consider such as in real time measurements and evaluation (for instance for feedback systems [Bad 74]). However, it usually gives better results [Gaf 66] than the second one if a long transient period (see 4.5) does not have to be eliminated before measurement are collected. This comes from the fact that since the two runs are independent (if different seeds are taken for the random number generators, see 5.2), there is no correlation between the two samples and for a same total simulation time, the first method gives smaller confidence intervals. It is sometimes possible to introduce a
negative correlation between several runs reducing further the confidence interval with antithetic or stratified sampling techniques [Kle 75].

With the antithetic technique, each time a random variable X is generated from a pseudo-random number r(0 < r < 1) through the inverse of the distribution function of X (see 5.2.2) :

$$X_i = F^{-1}(r_i) \quad (\text{at the i-th drawing})$$

then in a companion run (which actually runs in parallel with the first simulation), an antithetic variable X' is generated from the pseudo-random number 1-r

$$X'_i = F^{-1}(1-r_i)$$

Hence a large value of X will lead to a small one for X' (with respect to the mean value) and vice-versa, introducing therefore a negative correlation between the input variables. If several random variables with different distributions have to be generated during the simulation, the method can be applied to only one, to several or to all of the random variables. However the method can work only if the antithetic random variables can be produced at the same time on the two runs, otherwise the lack of synchronisation will eventually destroy the negative correlation which may be observed between the two samples [Fer 75].

### 4.5 - Initial state of a simulation run

It is often said that a simulation run should start from a state close to the "steady state of the system". Another method often used is to let the simulation run during a certain time before taking measurements to eliminate the "transient phase". These notions exist only in a probabilistic sense and even in this sense, they sometimes cannot be defined.
Let us try to clarify these intuitive notions in the light of the
formalism developed before.

In what follows, we shall consider that the state of the simulation
is defined by X(t,ω) although in general this vector should be extended to
consider all the variables needed to run the simulation and the measurements.
We shall also assume, to clarify the notations, that the vector X(t,ω) is
reduced to a single component and we shall hence consider the initial state
as the value of X(t₀).

The initial value X(t₀) depends on the way the simulation has been
coded: it may be fixed (X(t₀,ω) = X₀, ∀ω) or randomly drawn from an arbitrary
distribution (using pseudo-random numbers). In any case the random variable
X(t), t > t₀ usually depends (if the process is autocorrelated) on some
<u>known</u> information introduced in the simulator for the initial condition.

If we want to study the transient behavior (defined by the joint
distribution) of X(t) (preferably through several simulations runs) it is
important to remember that the process is conditioned on the initial distribu-
tion of X(t₀).

If on the other hand, we want to study the stationary behavior
of X(t) (assuming that it exists), the methodology is quite different when
we use a single simulation run or several. If we want to have at a specific
instant t a representative sample of the limiting random variable $\lim_{t \to \infty} X(t)$
through a number of simulation runs, the initial value
of each run must be drawn from the stationary distribution. Otherwise
the joint distributions of the random variable X(t) will be different from
the stationary joint distributions and the estimates obtained from the
sample collected at time t will be biased. Of course, the influence of the
initial distribution of X(t₀) on the joint distributions of X(t) will
diminish as $t - t_0 \to \infty$ but there is no strict rule to find when this
influence can be neglected. It depends in general on the initial distribution
and on the process studied. For example, in Palm's problem, the expected
number of busy lines $\xi(t)$ in a telephone exchange with Poisson arrivals
($\frac{1}{\lambda}$) and exponential service time ($\frac{1}{\mu}$) is given by :

$$E[\xi(t)/\xi(0) = k] = \frac{\lambda}{\mu} + e^{-\mu t} \left(k - \frac{\lambda}{\mu}\right)$$

In this case, the influence of the initial state decreases exponentially with the time and is proportional to $|E[\xi(0)] - \frac{\lambda}{\mu}|$.

The problem is different if the stationary behavior has to be estimated in a single simulation run. In this case and if the method of the regeneration points is used, there is no problem with the initial state (which is a regeneration point) : theoretical considerations (4.3.4) allow us to have proper estimates for the stationary stochastic process. Otherwise if standard statistic analysis is applied on a single sample path obtained from one simulation run, the initial value $X(t_0)$ could be arbitrarily chosen as any value which could be drawn from the limiting distribution.

However, consider as an example the number of customers $\xi(t)$ in a simple queueing system (e.g. $M/M/1$) : if we let the simulation start with a huge number of customers (e.g. $\xi(0) = 1000$ a very unlikely although possible value if $\frac{\rho}{1-\rho} \ll 1000$) then the confidence interval on the mean number of customers in the system will not include the limiting mean number at least during a long simulation period. This simply means that the sample path observed is a very unlikely one : we are in one of the cases which occurs with probability $\delta$ for a $1-\delta$ confidence interval if we assume that the initial state has been drawn from the steady state distributions (even if this is not done in the simulation run).

Hence, starting the simulation in state close to the most probable state is justified but this comes from the fact that we perform a single run which should be representative of all the possible runs. The idea of separating the simulation run into a transient phase and a steady state phase is not justified on theoretical grounds when a single run is performed unless the initial state belongs to a transient class the system will certainly leave and never come back to it.
V - ROBUSTNESS AND VALIDATION OF A SIMULATION MODEL

As it was said earlier, a discrete event simulation (used for performance evaluation) is a mechanism which should act, with respect to the measurements obtained and under the same conditions, as would the physical system which is modelled. Therefore it is of prime importance if we want to obtain a meaningful evaluation of the physical system -whether real or hypothetical- to make sure that this function of the simulator is reliable. We will see at the end of this chapter how to define and measure this but we will first discuss some modelling techniques which play an important role in the quality of the simulation.

There are several ways to look at performance predictions. One can observe a real system and by performing numerous measurements on the variables of interest, try to obtain a mathematical relation between all these variables through statistical analysis (factorial analysis, regressions,...). Although this approach which has been used mostly in computing centers can provide meaningful results, it is entirely different from the simulation approach. With simulation, we try to obtain performance prediction not only through the analysis of raw data, but through a knowledge and a modelling of the operation of the system itself. Using a mathematical model of the system (for instance with queueing theory), it may be possible to obtain exact formulas expressing the relations between several variables. With simulation, the relations are obtained experimentally through the execution of runs with different data.

Although the results desired are similar with those obtained from pure statistical analysis -we always get statistical relation between variables of interest- the approach by simulation can lead to more accurate modelling because information on the structure of the system is taken into account. However, for this approach to be fruitful, one depends on two conditions :
- a correct modelling of the system at the structural and probabilistic level

- and the availability of appropriate data.

While the first point is purely a question of modelling methodology (although it can be tested by statistical analysis, see 5.3), the second point depends on the way the data are collected, processed and given to the simulator. This data - which may be used to describe the environment (e.g. interarrival times of programs) or the operation of the system itself (e.g. service time at the disk) - can be measured directly on a real system and used in the simulator or extrapolated from known data if the measurements cannot be performed.

Although in some cases complete traces of data are used to drive the simulation runs, the available data are often reduced to a few parameters by using a model to generate artificial samples. The real and the artificial samples should look alike on a statistical basis and therefore there is no real difference between a trace driven and a model driven simulator except when one wants to control the experiments by changing the inputs. In particular, one should not consider a trace driven simulation as a "deterministic" simulation : the sample has been generated outside the model but is still of a probabilistic nature in the cases which we are considering.

Now we will see how to consider and model the environment of the system we wish to evaluate, how random variables are generated and the influence of mathematical and structural assumptions (robustness).

## 5.1 - Input variables of a simulator

Each time a physical system is modelled it is important, although not always simple, to distinguish the system itself from its environment. In a computer system for example, the environment can consist in the number of users, the characteristics of the program load, the response times of the operator, the computer room temperature... In view of this, a simulator can be considered as a black box transforming a number of input variables - the environment - into a number of output variables - the measurements of interest. Obviously, all the variables which should be considered to describe the environment,

<footnote>* Here, the term "variable" means either a stochastic variable, the characteristics of a distribution (mean, variance,...) or some fixed parameter.</footnote>
are the ones which can have an influence on the output variables of the simulator. This in turn, depends on the way the system has been modelled and more precisely, on the level of abstraction of the model.

Once the variables describing the environment have been selected (we suppose that at least one variable describing the environment has an influence on the measurements although in some cases this may not be the case), two approaches can be taken to execute a simulation run under specific conditions, that is, under a given environment :

- a trace (or time recording) of all the input variables can be obtained directly by measuring a real environment (if this is possible) and used as the input to the simulator,

- or a model of the input can be built and included in the simulator to produce a statistically identical input.

### 5.1.1 - Real input (trace)

In this case, the real environment of the system under study is available and the input variables of the simulator can be measured directly. The environment may be an existing one (for example programs taken in the job stream of a computer), or created artificially (for example synthetic jobs). In any case, measures are taken on this environment (for example the addresses referenced at each instruction or the times between two I/O) and input directly (the measurement procedure is included in the simulator), or indirectly (via some physical medium like a tape) to the simulator. One should notice that the input variables are not necessarily time dependent : one or several variables may take a single fixed value during one simulation run (for example a number of users in the system during the simulated time period).

Typical examples of computer systems using real input are the so called "trace driven simulators" [She 73, Bei 72, Ber 74, Noe 72, She 72,...], where the environment consists usually in the job load. The input variables are monitored on an existing system processing the loads of interest and are recorded on a tape which is read during the execution of the simulation. Trace data
usually consists in supervisor calls and sometimes in references to virtual pages. If the paging activity is modelled at a very fine level, the trace should contain the page number of all the memory references at each instruction.

### 5.1.2 - Modelled input

The previous approach can only be used when the real environment of the system under study is available or when it can be created. It also implies that the measurement of the needed input variables can be performed and this during a period equal in length to the one we wish to simulate. Often these conditions are not satisfied but even when they are, the trace approach does not allow to have a good understanding of the influence of the input on the measurements. Unless the trace data consists in a few fixed parameters, it does not provide a good characterization of the environment at a global level and we are left with the comparison of results obtained with a few "known" environments.

If we want to use another approach in the case where the input variables are stochastic variables varying in time, we must build a model of the environment. This model - included in the simulator - will furnish the required input variables as the simulated time advances. Although they can be issued from a deterministic procedure, the variables are usually derived from probabilistic models through the use of pseudo-random number generators. The main assumption usually made with these models is to consider that all the input variables are generated from stationary distribution functions (although not necessarily independent nor non-autocorrelated).

In computer system simulation, when the model is based on a network of queues and servers, one of the most commonly used input model is the life time function. This model expresses the virtual time spent by a program on the CPU between two events of the same type, as random variables whose distribution are known. For example in [Bad 74], the CPU service time is taken as the minimum of three independent random variables which express the virtual times until :

- the end of the execution of the program,

- the next I/O,

- the next page fault.
These three random variables are described by their distribution functions which may depend on several parameters (program parameters, number of pages allocated,...). These distributions can be estimated by analyzing real program trace data [Bur 76] with statistical tools [Ler 75].

If the simulation model is described at a deeper level of detail and needs program references to virtual pages, a page reference model has to be used if one does not want to use a real program trace. Among the models available we can mention the independent reference model, the LRU stack model, the working set model [Len 75] which have been validated on program trace data.

If we want the model to produce input variables as would do the real environment, it must be derived from a statistical analysis of a sufficiently large sample obtained by measuring the environment to be modelled. This approach can prove if the assumptions are valid (dependence, stationarity,...), but it also provides the model with the distribution functions which best fit the sample [Bur 76, Len 75, Ler 75].

However this statistical analysis still cannot be performed if the environment whether real or artificially created is not available or cannot be measured. Nevertheless, in this case mathematical modelling is the only possibility and the user will have to depend upon his intuition or upon extrapolation of known data and assumptions to build a suitable model. Anyway, we will see in 5.3 that if the model is robust, the output variables depend very little on the assumptions made in order to have a simple mathematical model [BuP 76]. In this case, the user is left only with the validation of a few critical parameters which he can anyway experiment on a wide range.

## 5.2 - Random variable generation

We have seen previously (5.1) that random variables may be used to model the operation of a system or its environment. We will now present in detail how to generate these variables according to a given probability distribution. One of the following approaches has to be used depending on the knowledge one has of the probability distribution :
1. if it is known through a real sample, this sample may be used directly in the simulation to give the random variables one by one. This is the trace driven approach and we already have seen its main drawbacks (5.1.1).

2. the distribution function may be given by its histogram (or a table of its distribution function). This histogram may have been obtained from a trace if we suppose the distribution function to be stationary. We will see (5.2.2) that we need random numbers uniformly distributed in [0,1] to use this method.

3. the distribution function (or the density), is given by a formula. In this case it may be possible to obtain the random variable through the composition of available random variables (e.g. uniformly distributed) or through the inverse $F^{-1}(U)$ of the distribution function :

$$F(x) = \text{Prob}\{X < x\} = \int_{\infty}^{x} f(t) dt$$

In this last case, a random variable uniformly distributed in [0,1] is also needed.

From this brief introduction, it can be noticed that uniformly distributed random variables play an important role and we will now focus on their generation.

### 5.2.1 - Pseudo random number generation [Knu 2, Mir 72]

The first random number generators were built on physical phenomena whose stochastic nature had been demonstrated (white noise, particle counting,...). Therefore, to be used by a computer program, special hardware had to be connected to the machine. Besides this difficulty, the lack of reproducibility made it quite difficult for programmers to test and correct their programs.
For these two main reasons, algorithms were soon developed to
produce in a deterministic way, sequences of numbers which appear to be
random. These numbers are called *pseudo-random numbers* and are considered
to simulate a sequence of independent uniformly distributed random
variables if their sequence can pass a number of statistical and structu-
ral tests (in practice, perhaps a dozen). This concept of pseudo-
randomness is still ill-defined but it is at the moment a subject of
active research.

The principle of these algorithms is to define a machine state
by one or several numbers or bit configurations and to change this state
in a deterministic way by some algorithm each time a new number is required.
A pseudo-random number is deduced from the machine state by some simple
rule. Since the number of states is necessarily finite, a major property
of pseudo-random number generators is that they produce periodic sequences :
as soon as the state is identical to some previous state, the random number
sequence will reproduce itself. Therefore, it is important to try to obtain
the longest possible period. This depends on the definition of the state
(the maximal period is equal to the number of possible states) but also on
the characteristics of the algorithm (not all algorithms give the maximal
period).

Since the first generation of computers, many algorithms have
been proposed to generate sequences of uniformly distributed random numbers.
The first methods such as the mid-square technique proposed by Von Neumann
(1948) or the mid-product technique (Von Neumann and Metropolis) did not
produce very good results with respect to the statistical tests. Furthermore,
it was not possible to study the characteristics of the sequences in a
theoretical way (period, autocorrelation function,...). In 1951, Lehmer
proposed the congruential technique which is still in use. However, it is
now challenged by a new method based on feedback shift registers and
by methods combining two random number generators (shuffling) which
can lead to very large periods [Lea 73].

We will now present these two techniques with their properties
and possible implementations.

### 5.2.1.1 - Congruential technique

Several variations on this technique exist, the most widely
used being the multiplicative congruential technique. With this method
the k-th pseudo-random number is given by the following recurrence
formula :

$$U_k = aU_{k-1} \pmod m \quad k = 1, 2, 3, \dots$$

$U_0$, the first value given to the generator is called the *seed*. Investigations
have been conducted to find the best values for a and m. These should insure
that the period is maximum but also they must lead to good statistical
properties of the series which is uniformly distributed between 1 and m-1
(a division by m gives numbers uniformly distributed in $]0, 1[$).

The computer word size usually gives the value for m if the
modulo operation is performed by truncation (i.e. $m=2^{32}$ for a 32 bit
machine). Actually if we want the maximum period, m must admit primitive
roots and a should be selected among them [Sto 73]. Therefore m should
be of the following form [Kur 76] :

$$m = \begin{cases} p^r & r = 1, 2, 3, \dots \\ 2p^r & p \text{ any prime number } \neq 2. \end{cases}$$

However this implies that m cannot be of the form $2^\alpha$. In order
to have a fast algorithm to perform the modulo operation with no division,
an algorithm due to Lehmer [Pay 69] requires that m is of the form :

$$m = 2^{q} \pm 1 \pmod{p}$$
The previous constraint (maximum period) leads to the following values for q :

$$q = 2, 3, 5, 7, 13, 17, 19, 31, 61, 89, \dots$$

Therefore for a 32 bit machine, m can be selected as :

$$m = 2^{31} - 1$$

An extension of this method can be found in [Kur 76] where m is of the form :

$$m = 2^{\alpha} \pm 2^{\beta} \pm 1$$

These two methods only insure that the period is maximum (and equal to m-1). The value of a must be selected among the primitive roots of m to give good statistical and structural properties to the sequence. Two severe tests are usually performed to check the quality of the generator [Knu 2, Kur 76]. These are the serial correlation test (first order autocorrelation coefficient) and the spectral test. The first one can be performed on a theoretical basis if the period is maximum [Knu 2]. However other statistical tests have to be performed before any generator is accepted.

We give below a list of values for a and m which lead to satisfactory results for the statistical tests. Two cases where m=32 do not give the maximum period but are simpler to implement.

<table><caption>Table 2 : Some selected parameters for the congruential method</caption><thead><tr><th>Multiplier a</th><th>Modulo m</th><th>Period T</th><th>Source</th></tr></thead><tbody><tr><td>2100005341</td><td>2<sup>31</sup>-1</td><td>2<sup>31</sup>-2</td><td>[Kur 76]</td></tr><tr><td>7<sup>5</sup>=16807</td><td>2<sup>31</sup>-1</td><td>2<sup>31</sup>-2</td><td>[Lea 73]<br>Assembler 360</td></tr><tr><td>2<sup>22</sup>-2<sup>14</sup>+4</td><td>2<sup>31</sup>-1</td><td>2<sup>31</sup>-1</td><td>[Kur 76]</td></tr><tr><td>69069</td><td>2<sup>32</sup></td><td>2<sup>29</sup></td><td>[Lea 73]<br>Assembler 360</td></tr><tr><td>162571342215<sub>a</sub></td><td>2<sup>32</sup></td><td>?</td><td>[Cov 67]</td></tr><tr><td>2<sup>27</sup>-2<sup>18</sup>-1</td><td>2<sup>35</sup>-2<sup>5</sup>+1</td><td>2<sup>35</sup>-2<sup>5</sup></td><td>[Kur 76]</td></tr><tr><td>2<sup>44</sup>-2<sup>28</sup>+1</td><td>2<sup>47</sup>-2<sup>7</sup>+1</td><td>2<sup>47</sup>-2<sup>7</sup></td><td>[Kur 76]</td></tr></tbody></table>
There exist several variations on this technique, two of the best known being [Knu 2, Mir 72] :

- the additive congruential method where

$$U_k = (U_{k-1} + U_{k-2}) \pmod{m}$$

- and the mixed congruential method :

$$U_k = (aU_{k-1} + c) \pmod{m}$$

However the congruential techniques give often relatively mediocre results to the spectral test [Mar 68]. Another method to generate pseudo-random numbers is now being used more and more frequently. We will present now this technique (called the feed-back shift register method) which gives much better results to the statistical test. It also has the advantage of giving very long periods and being very easy to implement either in software or in hardware.

### 5.2.1.2 - Feed-back shift register (FSR) method [Tau 65, Lew 73, Lea 73, Sto 73]

This method is based on feed-back shift registers. These registers give series of bits from which numbers are extracted. The operation of the register consists in performing a circular right shift and to change the leftmost bit as an XOR combination of selected bits of the register. The bit sequence is obtained by the successive values of one of the bits.

<figure><img src="image_9.png" /><figcaption>Fig. 6 : Feed back shift register (FSR)</figcaption></figure>
The maximal period of the bit sequence is $2^n - 1$ for an n bit register and is obtained by carefully selecting the bits to perform the feed-back. A theory based on Gallois Fields allow us to find the right combinations [Sto 73]. The content of the register is considered as the coefficients of a n-1 degree polynomial (mod. 2) and each step is equivalent to a division by an n degree polynomial (the content of the register represents the polynomial rest). The non-zero coefficients of this dividing polynomial indicate the bits to be selected. It is shown that the period is maximum (that is all the polynomial will be formed by successive divisions) if and only if the dividing polynomial is irreducible. Figure shows an example of maximum length cycle on a four bit register obtained with $X^4 + X + 1$ (i.e. bits 1 and 0)

<figure><img src="image_10.png" /><figcaption>Fig. 7 : 4 bit FSR based on $X^4+X+1$</figcaption></figure>
In this example, the bit sequence B consists in the repetition of any row (they are all the same if we consider that they are circular) :

$$ B = \{ \dots 0 \underbrace{1 \ 1 \ 1 \ 1 \ 0 \ 0 \ 0 \ 1 \ 0 \ 0 \ 1 \ 1 \ 0 \ 1 \ 0 \ 1 \ 1 \ 1 \ 1 \ 0 \dots} \}_{\text{1 cycle}} $$

The problem is now to extract combinations of bits from the sequence to interpret them as numbers. The first method consisted in reading the content of the register after each shift. This is equivalent to considering n consecutive bits in the sequence and moving one bit to the left along the sequence each time a new number is required. A better scheme proposed by Kendall [Ken 38] is to select consecutive but disjoint n-tuples along the sequence (this can be obtained by operating the register n steps between two readings). A variant of this method has been proposed by Tausworthe [Tau 65]. He considers L consecutive bits in the sequence ($L \le n$) and successive reading are spaced q digits apart (q must be prime with respect to $2^n-1$ to keep the maximum period property). However, a new method called GFSR (Generalized Feedback Shift Register) has been proposed recently [Lew 73] and gives good results to statistical tests. In this method, L bits are selected in the sequence ($L \le n$) in a non consecutive fashion. The L bits are picked out in the sequence regularly every p bits. Each time a new number is needed, the next L bits are picked out by moving all the selection one bit up the sequence. Very good statistical characteristics were obtained for $p \ge 100 n$ [Lew 73].

<figure><img src="image_11.png" /><figcaption>Fig. 8 : Bit selection in a FSR sequence</figcaption></figure>
Although this method seems to be complex to realize, its implementation is quite simple in software (or in firmware) and furthermore is machine independent. The number of bits L selected to form a number will certainly depend on the machine word size (e.g. L = 32 on a 32 bit machine) but the number n of bits in the register may be arbitrary. Hence very large periods may be obtained. The idea behind the implementation is to consider L registers which generate identical but delayed sequences, the delay between two consecutive register being equal to p. Therefore it is not necessary to memorize the sequence : the simultaneous operation of the registers generates the required L bits at each iteration.

At the implementation level, the i-th bits ($1 \le i \le n$) of the L registers are stored in one machine word. Hence n words are needed to simulate the L registers (each register consists of n bits).

<figure><img src="image_12.png" /><figcaption>Fig. 9 : Implementation of the GFSR method</figcaption></figure>
For simplicity we are looking for polynomial of the form $x^n + x^m + 1$ to perform the division. Therefore the simultaneous operation of the L registers consists in performing an XOR on words 0 and m to obtain word n-1 and to obtain all the other words (0, 1, ..., n-2) by a general shift up one word. In the implementation, the shift is replaced by the operation of a pointer to word 0 which runs in a circular fashion around the list.

To initialize the table one may start with any combination of bits (but all zeros) in register 1. The registers are then run p iterations and the initial configuration of register 1 is copied in register 2. This operation is repeated until the L registers are filled. After the initialization, a pseudo random number is obtained at each iteration by reading word 0. Its value is equidistributed between 1 and $2^n - 1$.

The following table gives possible values for n and m (m should not be too small nor close to n/2 [Too 71]):

<table><thead><tr><th>n</th><th>m</th></tr></thead><tbody><tr><td>20</td><td>5</td></tr><tr><td>30</td><td>9</td></tr><tr><td>39</td><td>14</td></tr><tr><td>49</td><td>9, 12, 15</td></tr><tr><td>60</td><td>11, 15, 17, 23</td></tr><tr><td>71</td><td>18, 20</td></tr><tr><td>81</td><td>16</td></tr><tr><td>98</td><td>11, 27</td></tr><tr><td>110</td><td>33</td></tr><tr><td>150</td><td>53</td></tr><tr><td>201</td><td>59</td></tr><tr><td>300</td><td>73, 75, 91</td></tr><tr><td>401</td><td>152</td></tr><tr><td>500</td><td>75, 95, 125</td></tr></tbody></table>

Table 3 : Selected primitive polynomials $x^n + x^m + 1$ [Zib 68]

The major problem with feedback shift register techniques is that they have not been extensively tested. Two other drawbacks are their requirement for larger memory space and the rather cumbersome and long initialization period.
### 5.2.1.3 - Shuffling random generators [McL 65]

The last method we present now is based on the operation of two generators. Its main advantages are to produce good quality pseudo random numbers even if the two generators are mediocre, and to give very large periods (the product of the periods of the two generators if the two periods are respectively prime). On the other hand, this method requires usually a larger computer time for the generation of each number.

The principle is to fill a table with entries generated by the best generator and to use the other generator to access this table. Each time a number is selected in the table, it is replaced by a new one issued from the first generator. Good results have been obtained with 128 entries in the table [McL 65]. In this case, the table is initialized with 128 numbers obtained from a generator $U(u_1, u_2, ..., u_{128})$. After the initialization, at the k-th drawing ($k = 1, 2, ...$), the second generator $V$ gives an index $I$ ($0 \le I \le 127$) for example by taking the seven high order bits of the pseudo random number $V_k$. This index is used to access a number $X_k$ in the table. $X_k$ is then replaced by a new pseudo random number $U_{128+k}$ obtained from $U$. The result of the k-th drawing is $X_k$.

In [Lea 73], very good results have been obtained with

- $U$ = FSR method where the dividing polynomial is:

$$X^{17} + X^{15} + 1 \text{ and the period } T_U = 2^{17} - 1$$

- $V$ = congruential method with :

$$a = 69069 \quad m = 2^{32} \quad \text{Period: } T_V = 2^{29}$$

$$\text{The resulting period is } T = T_U \cdot T_V = 2^{46} - 2^{29}$$

### 5.2.2 - Random variable generation

We have seen previously how to generate uniformly distributed random variables in $[1, m]$ or in $]0, 1[$ if we place a decimal point in front of
the binary number generated. We will now see how to generate random variables with arbitrary distribution functions.

In a few cases good approximations can be obtained directly from several uniformly distributed numbers. For example a Normal variable can be obtained by summing several uniformly distributed variable (Central Limit Theorem). Similarly Gamma and Beta variables can be obtained through the ratio of uniformly distributed variables [Kin 76]. We will not insist on these specific cases. A collection of such recipies can be found in [Knu 2] and [Mih 72].

The next method can be used whenever we can have an explicit formula for the inverse $F^{-1}(U)$ of the distribution function $F(x)$. At the k-th drawing, a pseudo random number $U_k$ uniformly distributed in $[0, 1]$ is generated. The result of the drawing is the variable $X_k$ obtained by :

$$X_k = F^{-1}(U_k)$$

Let us show that this variable has the desired distribution $F(x)$ :

$$
\begin{align*}
\text{Prob}\{X_k \le x\} &= \text{Prob}\{F^{-1}(U_k) \le x\} \\
&= \text{Prob}\{F[F^{-1}(U_k)] \le F(x)\} \quad \text{since } F(x) \text{ is strictly increasing} \\
&= \text{Prob}\{U_k \le F(x)\} \\
&= F(x) \quad \text{since } U_k \text{ is uniformly distributed in } [0, 1].
\end{align*}
$$

We give below a few examples of classical distribution laws with their inverse functions :

$$\text{exponential law} \quad : F^{-1}(U) = \left(-\frac{1}{\lambda}\right) \log U$$

$$\text{uniform law} \quad : F^{-1}(U) = a + U(b-a)$$

$$\text{Pareto's law} \quad : F^{-1}(U) = (U)^{-1/\lambda}$$

$$\text{Rayleigh's law} \quad : F^{-1}(U) = \sqrt{-\frac{1}{\lambda} \log U}$$

$$\text{Weibull's law} \quad : F^{-1}(U) = \left(-\frac{1}{\lambda} \log U\right)^{1/b}$$
The last method supposes that we have a tabulation of the distribution function (or a cumulative histogram) F(X).

<figure><img src="image_13.png" /><figcaption>Fig. 10 : Tabulated distribution function method</figcaption></figure>

Hence we have in a table, couples [xᵢ, F(xᵢ)] xᵢ < xᵢ₊₁
i = 1, 2, ..., n.

The algorithm consists in generating at the k-th drawing a
pseudo random number U_k uniformly distributed in [0, 1] and to find the
couple [x_j, F(x_j)] such that

$$
F(x_j) < U_k \le F(x_{j+1})
$$

The result of the drawing is

$$
X_k = x_j
$$

This method is particularly will suited for discrete random variables
especially if the number of values of x is limited. For continuous random
variables, its precision will depend on the number of entries in the table.
## 5.3 - Calibration and validation

However carefully a simulation model has been designed and programmed, one may never be completely sure of the accuracy of its predictions. To gain some confidence in its results, it must be tested against some aspects of the reality it is supposed to represent. These comparisons may be performed with respect to measurements done in real life situations or with respect to a dependable (validated) model.

The first step is the calibration and consists in the "turning" of the simulator. This step is executed each time some parameters or even the model itself have to be adjusted. This is done by comparing the output of the simulator and real life (or another model) for a number of different operating conditions. When all the comparisons are favorable, the simulator is said to be calibrated. The second step, called validation, consists in checking that the comparisons are still valid for different operating conditions not included in the calibration. This whole procedure (calibration + validation) is often referred to as the validation of the model but the two steps have to be distinguished each time the simulator needs some adjustments [Bei 73, Ern 71].

We will now present in detail how the comparisons can be performed. In the first case, the simulator is validated against a mathematical model. In the second case, the output of the simulator is compared to another sample obtained either from another simulator or from real life.

### 5.3.1 - Validation with a mathematical model

The solutions of mathematical models (e.g. queueing networks) are obtained by mathematical formulas which usually give mean values, sometimes variances and seldom complete distribution functions of the stochastic variables used for performance evaluation (queue lengths, waiting times, resource utilizations,...). Therefore, the validation has to be performed with respect to the numerical values obtained from the model.
The procedure goes as follows for a given environment :

1 - the simulation model is run under the same hypotheses as the
mathematical model each time this is possible (Poisson arrivals,
exponential service times,...),

2 - the performance variables which can be computed from the mathematical model are estimated with the samples obtained from the simulation runs,

3 - these estimates are compared to the values predicted by the mathematical model under the same environment conditions. For this, the confidence intervals have to be used (cf. 4.3) : the comparison is good if the predicted value falls within the confidence interval at least n°/o of the times (for an n°/o confidence level).

Notice that those three steps are only conceptually different but
in practice, they can be merged together in a simulator for efficiency
purposes.

5.3.2 - Validation with another simulator or with real life

In these cases we are comparing several samples of the performance
variables. For each variable we must have two or more samples collected
from the simulator to be validated and from another simulator or from
real life but in the same environment conditions.

Several methods can be used to compare the samples [Kle 73]. We
present now two of these methods based on variance analysis and examples
where they have be used successfully.

The first method consists in representing, in the form of a
multilinear regression, the differences measured between two samples of a
same output variable (in general a performance variable because the model
is supposed to predict such variables) [Lin 63, Bei 72]. The difference is
expressed with respect to time :

$$
\epsilon(t_i) = Y(t_i) - Y'(t_i) \quad (i = 1, 2, \dots, N)
$$
where $Y(t)$ is obtained from the simulator being validated and $Y'(t)$ is the reference sample (from a simulation or from real life).

This difference $\epsilon$ is the variable to be explained by the explicative variables $X_j(t)$ ($j = 1,...,p$). These p variables are the ones on which it is possible to have an influence in the simulation model during the calibration stage.

In the multilinear regression model, $\epsilon$ is of the form :

$$ \epsilon = \sum_{j=1}^{p} a_j X_j $$

where the $a_j$ are computed to minimize a square sum norm, that is, such that :

$$ \sum_{i=1}^{N} \left[ \epsilon(t_i) - \sum_{j=1}^{p} a_j X_j(t_i) \right]^2 \text{ is minimal.} $$

The values can be easily obtained as the solution of a linear system whose matrix is the covariance matrix of the variables $X_j$. A detailed description of the theory can be found in [Lin 63].

Although the variables $X_j$ are rarely completely independent (this can be tested by computing the correlation coefficients $\rho(X_i, X_j)$), $a_j$ represents the influence of the variable $X_j$ on the error $\epsilon$. More complex models for example with a quadratic expression for $\epsilon$ can take into account the correlation between $X_i$ and $X_j$ but do not seem to lead to better results.

If $a_j$ is close to zero, it means that the values of $X_j$ are correct and do not participate to the error on the performance variable $Y$. Otherwise, one must try to modify the values of $X_j$ in order to decrease the value of $|a_j|$. This must be done by modifying the simulator (either some parameter values or its structure) in order to obtain the desired effect. When all the $a_j$ are sufficiently small, the simulator is considered to be calibrated and it can be validated with the same technique.
This method has been used in [Bei 72] to validate a simulator model of OS/MVT-LASP with respect to real life measurements. The calibration was performed on the parameters of the model operation (error on response time is expressed as a function of several parameters such as tape set up time) and on the load parameters (errors on waiting or execution times are expressed as functions of the number of steps per job, memory space required, number of tape or disk mounts per job stream).

The principle of the second approach is quite different. One tries to explain the performance variable Y not as a function of explicative variables which can be adjusted but as a function of the environment and of the operating conditions represented by variables $X_1, ..., X_p$

$$Y = f(X_1, ..., X_p) + e$$

where e is an error term expressed by a random variable. The variables $X_i$ need not be continuous variables such as in the previous method. Therefore it is possible to represent the influence of a particular operating condition or different types of environment. For example a variable $X_i$ may take the value 1 if a control scheme is in action and 0 if it is not.

The models are computed by measuring Y for different values of the variables $X_1, ..., X_p$ in the simulator and in the reference system (simulator or real life). If both models are identical, the simulator is calibrated (or validated).

The number of experiments to perform in order to compute the model depends on P and on the number of values each variable may take. If all the possible combinations are examined, this may lead to an impractical number of runs. This is why in practice the design of the experiment may not include all possibilities (factorial design) and the method is used with variables which may take only two values. In this case, if p=3 and $X_i = 1$ or 2 ($i = 1, 2, 3$), the model can be chosen as a linear expression of the form :
$$Y_{ijk} = \mu + \alpha_i + \beta_j + \gamma_k + (\alpha\beta)_{ij} + (\beta\gamma)_{jk} + (\alpha\gamma)_{ik} + e_{ijk} \quad i, j, k = 1, 2$$

where $\alpha_i$ is the influence of $X_1$ at level i, $\beta_j$, the influence of $X_2$ at level j, ..., $(\alpha\beta)_{ij}$ the interaction of $X_1$ and $X_2$ when $X_1$ is at level i and $X_2$ at level j and $e_{ijk}$ is the error term.

The main difference with the previous method is that instead of testing the simulator for one specific environment, the simulator is tested for all the conditions specified in the experiment design. A simulator of CP-CMS has been validated with this method [Sch 75]. Five two level variables have been used to model different operating conditions (operation of the multiprogramming level control, single or double level CPU queue,...). The performance variables considered were percentages of time spent by the CPU in problem state and supervisor state and the number of page faults per second.

The model is calibrated if all the significant terms of the model are the same for the simulator and for the reference system. If this is not the case, the variable or variables involved in the discordance can indicate which part of the simulator has to be modified to obtain better results.

## 5.4 - Robustness

The results of the second part of the validation process (time validation by opposition to the calibration step) are positive when the model is robust. This is a property of the simulator which indicates that the performance results are not critically dependent on the modelling hypotheses. This means for example that the results sought for are not greatly influenced by small changes in the level of detail selected or by the form of the distributions used to describe the environment or the operation of the system.
A robust and validated model will allow us to have a good confidence in the results predicted, even if the operating conditions are distant from the ones tested during validation. Therefore, although this is not simple, it may be interesting to have some indications on the robustness of a simulator.

At the moment, there is no general method to check the robustness. A few experimental results exist on the mathematical assumptions in queueing networks [Buz 74, BuP 76]. We know for example that in closed queueing network, server utilization ratios and service rates are not very sensitive to service time distributions except for the mean value (on the other hand, the waiting times are very dependent on the variation coefficients of these distributions). We also know that it is possible to replace a finite source of customers by an infinite source with a good accuracy.

Hence the robustness must be tested for each individual simulation. The methods presented before in 5.3 can be used to test some of the hypotheses one at a time and to check their influence on the performance variables.

Such an approach has been taken in [Bei 72] to test the influence of the system load parameters on the performance variables.
## VI - CONCLUSION

In this paper we have presented some of the problems which may be encountered in developing discrete event simulation models. We also have mentioned or described in some details the techniques currently available to solve these problems.

Although the discrete event simulation technique is described and even formalized in some details, this paper is addressed to users already familiar with simulation and who are facing a particular problem. More specifically, we have investigated the problem of writing a synchronization kernel, selecting a simulation language, and analysing the results. This last problem is one of the most important and often one of the most neglected topics, therefore, we have placed more emphasis on it.

Many areas are still open for further research in simulation. In particular, the techniques for data analysis and random number generation need to be further investigated. New tools should also be developed to simplify the writing and the testing of the complex simulators which are required to design and evaluate large scale systems in computer and other fields.

**Acknowledgements** : We would like to mention that this survey has been greatly influenced by many discussions and editorial comments. Among others, we want to thank specially many friends working on related subjects and in particular Y. Kurita, J. Lenfant, D. Potier and A. Schroeder. We are also grateful to R. Muntz and P. Lewis for the valuable discussions we have had with them and for all the informations they gave us to improve this paper.
[And. 65] ANDO A., FISHER F.M. - "Near decomposability, partition and aggregation, and the relevance of stability discussions" - International economic review - Vol. 4, n° 1, January 1963.

[Avi. 73] AVI-ITZEAK B., HEYMAN D.P. - "Approximate queueing models for multiprogramming computer systems" - Operations Research. - Vol. 21, n° 6, Nov.-December 1973.

[Bab. 75] BABICH A.F., GRASON J., PARNAS D.L. - "Significant event simulation" - ACM - Vol. 18, n° 6 - June 1975.

[Bad. 75] BADEL M., GELENBE E., LEROUDIER J., POTIER D. - "Adaptive optimization of a time-sharing system's performance" - IEEE Proceedings on Interactive Computer Systems - June 1975.

[Bad. 75] BADEL M., SHUM A.V.Y. - "Accuracy of an approximate computer system model" - Rapport IRIA-LABORIA n° 130 - Juillet 1975.

[Bae. 74] BAER J.L. - "Models for the Design, Simulation, and Performance of Distributed Function Architecture" - Computer - March 1974.

[Bar. 71] BARD Y. - "CP-67 Measurement and Analysis - Overhead and Throughput" - ACM-SIGOPS - Workshop on System Performance Evaluation - April 1971 - Harvard University.

[Bei. 73] BEILNER H. - "Problems in calibrating and validating simulation models" - 2nd Seminar on Experimental simulation and solving of Probability problems - Liblice (CSSR) - June 1973.

[Bei. 72] BEILNER H., WALDBAUM G. - "Statistical methodology for calibrating a trace driven simulator of a batch computer system" - IBM Research - RC-3855 - May 1972.

[Ber. 74] BERETEVAS T. - "A simulation model representing the OS/VS2 release 2 control program" - Proceedings of an International Symposium on Operating Systems - Rocquencourt - Springer-Verlag - April 1974.

[Bir. 69] BIRTWISTLE G. - "Notes on the SIMULA language" - Publication n° S.8 - Norwegian Computing Center - OSLO 3 NORWAY - April 1969.

[Blo. 76] BLOOMFIELD P. - "Fourier analysis of time series : an introduction" - John Wiley and Sons - 1976.

[Boi. 73] BOI L., DRUCBERT J.P. - "OSSYOCSCOPE : System on auxiliary processors for measuring operating systems" - International computing symposium 1973 - Davos - Switzerland - September 1973 -

[Bor. 71] BORDSEN D.T. - "Univac 1108 instrumentation system" - SIGOPS-ACM - Workshop on system performance evaluation - April 1971 - Harvard University.

[Box. 70] BOX G.E., JENKINS G.M. - "Time Series Analysis Forecasting and Control" - Holden-Day - 1970.
[Bro. 63] BROWN R.G. - "Smoothing, forecasting and prediction of discrete time series" - Prentice-Hall Inc. - Englewood Cliffs - N.J. 1963.

[Bur. 76] BURGEVIN P., LERODIER J. - "Characteristics and models of program behavior" - National Conference ACM 76 - Houston, Texas - Oct. 1976.

[Buz. 74] BUZEN J.P., GOLBERG P.S. - "Guidelines for the use of infinite source queuing models in the analysis of computer system performance" - National Conference 1974.

[BuP. 76] BUZEN J.P., POTIER D. - "Accuracy of the exponential approximation" - Rapport IRIA-LABORIA (to appear).

[Buz. 76] BUZEN J.P. - "Fundamental laws of computer system performance" - Proceedings of International Symposium on Computer Performance Modeling Measurement and Evaluation - Harvard University - March 1976.

[Cha. 74] CHANDY K.M., HOWARD J., KELLER T.W., TOWSLEY D.J. - "Local balance, Robustness, Poisson departure and the product form in queueing networks" - Internal Report of the University of Texas at Austin - 1974.

[Chu. 68] CHUNG K.L. - "A course in probability theory" - Harcourt, Brace and World Inc. - New York - 1968.

[Con. 63] CONWAY R.W. - "Some tactical problems in digital simulation" - Management Sciences - Vol 10, n° 1 - October 1963.

[Cou. 71] COURTOIS P. - "On the near-complete decomposability of networks of queues and of stochastic models of multiprogramming computing systems" - NBLE Report - Belgium - November 1971.

[Cov. 67] COVEYOU R.R., MACPHERSON R.D. - "Fourier analysis of uniform random number generator" - JACM, Vol. 14, n° 1, January 1967.

[Cox. 68] COX D.R., LEWIS P.A. - "The Statistical Analysis of Series of Events" - M.S. Bartlett's - Methuen monographs - London 1968.

[Cox. 72] COX D.R., MILLER H.D. - "The theory of stochastic processes" - Chapman and Hall - London 1972.

[Cra. 74] CRANE M.A., IGLEHART D.L. - "Simulating stable stochastic systems I : General multiserver queues - II : Markov chains" - JACM - Vol. 21, n° 1 - January 1974.

[Cra. 75a] CRANE M.A., IGLEHART D.L. - "Simulating Stable Stochastic Systems : III Regenerative Process and Discrete Event Simulations" - Operations Research - Vol. 23, N° 1 - February 1975.

[Cra. 75b] CRANE M.A., IGLEHART D.L. - "Simulating stable stochastic systems IV : approximation techniques" - Management Science - Vol. 21, n° 11 - July 1975.
[Dah. 66] DAHL O.J., NYGAARD K. - "Simula - an Algol based simulation language" - CACM - Vol. 9, n° 9 - September 1965.

[Dah. 68] DAHL O.J., MYHRHAUG B., NYGAARD K. - "SIMULA 67 Common Base Language" - Publication n° S-2 - Norwegian Computing Center - OSLO 3 - NORWAY - May 1968.

[Dew. 71] DEWAN P.B., WYATT J.B. - "OSSL : Operating System Simulation Language, a user's guide" - Report RS I-71 - November 1971 - University of Houston (Texas).

[Dew. 72] DEWAN P.B., DONAGHEY C.E., WYATT J.B. - "OSSL - A specialized language for simulating computer systems" - Spring Joint Computer Conference 1972.

[Dim. 64] DIMSDALE B., MARKOWITZ H.M. - "A description of the SIMSCRIPT language" - IBM Syst. Journal 3, 1 - 1964.

[Ern. 71] ERNST H.A. - "Problems in computing systems performance modelling : characterization, calibration and validation" - IBM Research - RC 3319 - April 1971.

[Fer. 75] FERRARI D., LINDSAY D.S., PEAK P.M. - "Variance Reduction Techniques in Computer System Performance Measurements" - International Computing Symposium 1975 - Antibes - France - North-Holland Publishing Company - June 1975.

[Fis. 65] FISHMAN G.S., KIVIAT P.J. - "Spectral analysis of time series generated by simulation models" - Rand Corporation - Memorandum - RM - 4395 - PR - February 1965.

[Fis. 67] FISHMAN G.S. - "Problems in the statistical analysis of simulation experiments : the comparison of means and the length of sample records" - CACM - Vol. 10, n° 2 - February 1967.

[Fis. 68] FISHMAN G.S. - "Digital computer simulation : the allocation of computer time in comparing simulation experiments" - Operation Research - Vol. 16, n° 2 - erratum in Vol. 16, n° 5 - March 1968.

[For. 68] FORRESTER J.W. - "Industrial Dynamics" - M.I.T. Press. Student's edition - 1968.

[Fos. 74] FOSTER D.U., Mac GEHEARTY, SAUER C.H., WAGGONER C.N. - "A language for analysis of queuing models" - Report TR33 - University of Texas at Austin - April 1974.

[Fou. 67] FOURGEAUD D.C., FUCHS A. - "Statistique" - Dunod 1967.

[Gaf. 66] GAFARIAN A.V., ANKER C.J. - "Mean value estimation from digital computer simulation" - Operation Research - Vol. 14, n° 25 - 1966.

[Gau. 73] GAUDEUL B. - "Evaluation des performances des systèmes par méthode Hardware" - Compte rendu des Journées d'Etudes sur les recherches en mesure, simulation, modélisation des systèmes informatiques - Toulouse - Juin 1973.
[Knu. 2] KNUTH D.E. - "The Art of computer programming - Seminumerical algorithms - Random numbers" - Vol. 2, Chapter 3 - Addison-Wesley Publishing Company.

[Kur. 74] KURITA Y., MORIMURA M. - "Generation and testing of random numbers" - Proceedings of Symposium on Pseudo-Random Number Generators - Computing Center of Tokyo University - Tokyo - February 1974.

[Kur. 76] KURITA Y. - "Consideration of pseudo-random number generators : choice of parameters of multiplicative congruential methods and its statistical tests" - Journal of Information Processing Society of Japan - 1976.

[Kos. 73] KOSY D.W. - "An Interim Empirical Evaluation of ECSS for Computer System Simulation Development" - Symp. on the Simulation of Comp. System - 1973.

[Lav.75a] LAVENBERG S.S., SLUTZ D.R. - "Introduction to regenerative simulation" - IBM Journal on Res. and Develop. - September 1975.

[Lav.75b] LAVENBERG S.S., SLUTZ D.R. - "Regenerative simulation of a queuing model of an automated tape library" - IBM Journal on Res. and Develop. - September 1975.

[Lav.75c] LAVENBERG S.S., SHEDLER G.S. - "Derivation of confidence intervals for work rate estimators in a closed queuing network" - SIAM Journal of comp. - Vol. 4, n° 2 - June 1975.

[Lea. 73] LEARMONTH G.P., LEWIS P.A.W. - "Statistical tests of some widely used and recently proposed uniform random number generators" - Report MP555LW73III1 - Naval Postgraduate School Monterey - California - Nov. 1973.

[LeG. 68] LE GALL P. - "Convergence des simulations et application aux réseaux téléphoniques" - AFIRO monographies de recherche opérationnelle - n° 7 - 1968.

[Leh. 59] LEHMANN E.L. - "Testing statistical hypotheses" - John Wiley and Sons, Inc. - New York - 1959.

[Len. 75] LENFANT J., BURGEVIN P. - "Empirical data on program behavior" - International Computing Symposium 1975 - Antibes - France - North-Holland Publishing Company - 1975.

[Ler. 73] LEROUDIER J. - "Analyse d'un système à partage de ressources" - RAIRO - October 1973.

[Ler. 75] LEROUDIER J., SCHROEDER A. - "A statistical approach to the estimation of service time distributions for operating system modelling" - International Computing Symposium 1975 - Antibes - North-Holland Publishing Company - France - June 1975.

[Lew. 69] LEWIS P.A.W., GOODMAN A.S., MILLER J.M. - "A pseudo-random number generator for the system/360" - IBM System Journal n° 2 - 1969.

[Lew. 73] LEWIS T.G., PAYNE W.H. - "Generalized feedback shift register pseudo-random number algorithm" - Journal of ACM, Vol. 20, n° 3 - July 1973.
[Lin. 63] LINNIK Y.V. - "Méthode des moindres carrés" - Dunod - Paris - 1963.

[Mar. 64] MARKOWITZ H.M., HANSNER B., KARR H. - "Simscript" - Prentice-Hall - Englewood Cliffs - N.J. 1964.

[Mar. 66] MARKOWITZ H.M. - "Simulating with SIMSCRIPT" - Management Science - Vol. 12, n° 10 - June 1966.

[Mar. 68] MARSAGLIA G. - "Random number fall mainly on the planes" - Proc. Nat. Acad. Sci. 61, 1 - Sept. 1968.

[Mar. 72] MARSAGLIA G. - "The structure of linear congruential sequences" - in S.K. Zaremba, Ed., Application of Number Theory to Numerical Analysis, New-York : Academic Press - 1972.

[McD. 73] Mac DOUGALL M.H., Mac ALPINE J.S. - "Computer System Simulation with ASPOL" - Symposium on the Simulation of Computer System - June 1973.

[McD. 70] Mac DOUGALL M.H. - "Computer System Simulation : an Introduction" - Comp. Surveys - Vol. 2, n° 3 - September 1970.

[McL. 65] Mac LAREN M.D., MARSAGLIA G. - "Uniform random number generators" - Journal of ACM, Vol. 12, n° 1 - January 1965.

[Mer. 73] MERMET J. - "Etude méthodologique de la conception assistée par ordinateur des systèmes logiques" - Thèse d'état - Grenoble - Avril 1973.

[Mih. 72] MIHRAM G.A. - "SIMULATION : statistical foundations and methodology" - Academic Press. - 1972.

[Nay. 67] NAYLOR T.H., FINGER J.M. - "Verification of computer simulation models" - Management Science - Vol. 14, n° 2 - Oct. 1967.

[Nee. 66] NEELY P.M. - "Comparison of several algorithms for computation of means, standard deviations and correlations coefficients" - CACM - Vol. 9, n° 7 - July 1966.

[Nie. 70] NIELSEN N.R. - "ECSS : An extendable computer system simulator" - Rand Corporation - Memorandum - RM-6132-NASA - February 1970.

[Noe. 72] NOE J.D., NUTT G.J. - "Validation of a trace driven CDC 6400 simulation" - SJCC 72.

[Noe. 73] NOE J.D., NUTT G.J. - "Macro E-Nets for Representation of Parallel Systems" - IEEE Transactions on computer 73 - Vol. C-22, n° 8 - August 1973.

[Nut. 72] NUTT G.J. - "Evaluation nets for computer system performance analysis" - Fall Joint Computer Conference - 1972.

[Nut. 73] NUTT G.J. - "The computer System Representation Problem" - Symp. on the Simulation of Comp. System - June 1973.
[Nut. 74] NUTT G.J. - "Evaluation net simulation system reference manual" - University of Colorado - Dept. of Comp. Sci. - Tech. Report n° CU-CS-042-74 - April 1974.

[Pag. 64] PAGE E.S. - "On Monte-Carlo methods in congestion problems : II - Simulation of queueing systems" - Operation Research 13 - n° 2 - Avril 1964.

[Par. 74] PARENT M.R. - "The control graph models : a unified approach to performance evaluation" - Ph.D. Thesis - Case Western Reserve - University - Cleveland - Ohio - June 1974.

[Par. 61] PARZEN E. - "Mathematical considerations in the estimation of Spectra" - Technometrics - Vol. 3, n° 2 - May 1961.

[Pee. 74] PEEBLES R.W., BUCK D.L., MANNING E. - "Simulation studies of a homogeneous computer network for data sharing" - CCNG Report E-19 - University of Waterloo - Ontario - Canada - March 1974.

[Pet. 62] PETRI C.A. - "Kommunikation mit automaten" - Thèse de doctorat - Université de Bonn - 1962.

[Pot. 76] POTIER D., LEROUDIER J., BADEL M. - "Un modèle d'analyse des performances d'ordinateurs multiprogrammés à mémoire virtuelle" - Rapport Laboria n° 152 - IRIA - France - Jan. 1976.

[PrS. 69] PRIESTLEY M.B., SUBBARAO T. - "A test for non-stationarity of time series" - J. of the Royal Statistical Society - Vol. B31, n° 1 - 1969.

[PrK. 69] PRITSKER A.A., KIVIAT P.J. - "Simulation with GASP II" - Prentice Hall - Englewood Cliffs - N.J. - 1969.

[Red. 73] REDDY Y.V., BRYAN R.H. - "DESPL/1 : A discrete simulation language based on PL/1" - SCSC - Summer Computer Simulation Conference - Montréal - Canada - July 1973.

[Reh. 68] REHMANN S.L., GANGWERE S.G. - "A simulation study of resource management in a time-sharing system" - FJCC 68.

[Ret. 74] RETTENMAYER J.W. - "SIMPL/1 : a simulation programming language" - Winter Simulation Conference - Washington D.C. - January 1974.

[Ric. 73] RICHARDS F.M., WEATHERBEE J.E., WILLIAMS T. - "A validation procedure for discrete digital simulation models" - 1973 SCSC - Summer Computer Simulation Conference - Montréal - Canada - July 1973.

[Rod. 72] RODRIGUEZ-ROSELL J., DUPUY J.P. - "The evaluation of a time sharing page demand system" - Spring Joint Computer conference - 1972.

[Ros. 72] ROSE C.W. - "LOGOS and the software engineer" - Proceedings of the FJCC - Vol. 41 - 1972.

[Sce. 67] "SCERT - System and Computers Evaluation Review and Technique" - Compress Inc. - Washington DC - 1967.

[Sch. 65] SCHERR A.L. - "An analysis of time-share computer system" - Thesis Project MAC - M.I.T (1965).
[Sch. 75] SCHATZOFF M., TILLMAN C.C. - "Design of experiments in simulator validation" - IBM Journal of Research and Development - May 1975.

[Sea. 69] SEAMAN P.H., SOUCY R.C. - "Simulating operating systems" - IBM System Journal - n° 4 - 1969.

[She. 72] SHERMAN S.W., BASKETT III F., BROWNE J.C. - "Trace driven modeling and analysis of CPU scheduling in a multiprogramming system" - Communications of ACM - Vol. 15 - n° 12 - Dec. 1972.

[She. 73] SHERMAN S.W., BROWNE J.C. - "Trace driven modeling : review and overview" - Symposium on the Simulation of Computer Systems - Gaithersburg, Md - June 1973.

[Sim. 61] SIMON H., ANDO A. - "Aggregation of variables in dynamic systems" - Econometrica - Vol 20, n° 2 - April 1961.

[Slu. 68] SLUTZ D.R. - "The flow graph schemata model of parallel computation" - Ph. D. Thesis - Massachusetts Institute of Technology - Sept. 1968.

[S3. 67] S3 : System and software simulator - Defense Document Clearinghouse - AD 679269 - Sept. 1967.

[Sto. 73] STONE H.S. - "Discrete mathematical structures and their applications" - Science Research Associates - Palo Alto - 1973.

[Tau. 65] TAUSWORTHE R.C. - "Random numbers generated by linear recurrence modulo two" - Math. Comp. 19 - 1965.

[Tei. 66] TEICHROEN D., LUBIN J.F. - "Computer Simulation - Discussion of the Technique and Comparison of Languages" - Communications of ACM - Vol. 9, n° 10 - October 1966.

[Teo. 73] TEOREY T.J., MERTEN A.G. - "Considerations on the Level of Detail in Simulation" - Symp. on the Simulation of Comp. System 1973.

[Too. 71] TOOTILL J.P.R., ROBINSON W.D., ADAMS A.G. - "The runs up and down performance of Tausworthe pseudo-random number generators" - Journal of ACM - Vol. 18, n°3 - July 1971.

[Van. 71] VAN HORN R.L. - "Validation of simulation results" - Management Science - Vol 17, n° 5 - January 71.

[Vau. 75] VAUCHER J.G., DUVAL P. - "A comparison of simulation event list algorithms" - CACM - Vol. 18, n° 4 - April 1975.

[Wal. 72] WALBAUM G., BEILNER H. - "SOUL : a simulation of OS under LASP" - IBM Research RC 3810 - March 30 - 1972.

[Wat. 70] WATTERS L.J., VASILIK M.V. - "A Stochastic Network Approach to Test and Chackout" - Fourth Conference on Applications of Simulation - New York - Dec. 1970.

[Wul. 69] WULF W.A. - "Performance Monitors for Multi-programming System" - Second Symposium on Operating System Principles - Princeton Univ. - Oct. 1969.

[Wya. 72] WYATT J.B., DEWAN P.B. - "A teaching Tool for Interactive Computer System Design" - ONLINE 72 - Conference Proceedings.
[Zib. 68] ZIELER N., BRILHART J. - "On primitive trinomials" - Information and Control - n° 13 - 1968.