IRIA
laboria

laboratoire de recherche
en informatique
et automatique

Quelques problèmes liés
à la simulation de modèles
de systèmes informatiques

M. Badel

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

Rapport de Recherche n°144

Novembre 1975

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
78150 - Le Chesnay
Tél. 954.90.20
# QUELQUES PROBLEMES LIÉS A LA SIMULATION DE MODELES DE SYSTEMES INFORMATIQUES *

M. Badel

## RÉSUMÉ :

Le travail présenté décrit deux applications de la simulation de modèles de systèmes informatiques. Une première partie évoque les problèmes généraux posés par la simulation de systèmes informatiques et introduit les outils statistiques utilisés dans la suite de l'étude.

La première application exposée qui est la validation réciproque d'un modèle mathématique et d'une simulation, donne un exemple de mise en œuvre des méthodes précédentes et montre comment utiliser la simulation pour obtenir une précision statistique donnée.

La seconde application montre un exemple d'étude par simulation d'un algorithme de régulation adaptive d'un système, à partir des mesures de performances prises en sortie. Des techniques statistiques inspirées des précédentes et permettant des estimations avec oubli du passé y sont exposées.

## ABSTRACT :

Two applications of computer system simulation are presented in this work. An introduction presents the general problems in computer system modeling by simulation and deals with the statistical tools used in the rest of the study.

The first application is the mutual validation of a mathematical model and a simulation model, and it gives an example of the use of the statistical methods presented before to obtain a given statistical precision.

The second application shows an example for the study of an adaptive control scheme for a computer system based on performance monitoring. Statistical methods close to the previous ones and allowing the use of smoothed estimates are presented.

<footnote>* Thèse soutenue le 28 mai 1975 à l'Université de PARIS VI devant un jury composé de :
MM. J. ARSAC (Président), J.C. SIMON, E. GELENBE, S. KRAKOWIAK, E. DIDAY, J. LEROUDIER,</footnote>
PRESENTATION

L'étude objet de la thèse décrit quelques problèmes liés à la
simulation de modèles de systèmes informatiques et développe deux appli-
cations. Une attention toute particulière est portée aux mesures et à
quelques unes de leurs utilisations possibles.

Une première partie générale de l'étude situe le problème de la simulation de systèmes informatiques et précise les buts recherchés. Le but retenu par nous est l'évaluation de performances, ce qui suppose la mise en oeuvre de techniques de mesures et de prises de statistiques, et pose le problème de la qualité des résultats obtenus. Les différents problèmes liés à la discussion de la qualité des résultats sont d'abord passés en revue à travers la bibliographie. Le problème de la détermination d'une précision statistique sur les résultats (construction pratique d'intervalle de confiance) est étudié en détail. Après quelques rappels statistiques, la technique retenue pour la suite de l'étude (groupement des mesures afin de diminuer l'effet des corrélations) est entièrement exposée.

La suite de l'étude donne deux exemples d'application de la
simulation de modèles de systèmes informatiques, et des techniques décrites
dans le chapitre I. Nous précisons d'abord les modèles mis en oeuvre
(extensions de modèles mathématiques à serveurs et à files d'attente) et
leurs techniques particulières de simulation et de prises de mesures. Le
second chapitre en donne ensuite une première application au contrôle
réciproque d'une résolution mathématique approchée et d'une simulation.
Après présentation du modèle objet de la comparaison, la méthode de
résolution théorique, empruntée à la bibliographie, (approximation du
comportement du modèle par un processus de diffusion) est rappelée. La
simulation est ensuite utilisée pour la comparaison des résultats : taux
d'utilisation des serveurs et probabilités d'état des files d'attente. La
conclusion tirée après expériences est que, dans les hypothèses du modèle,
la concordance des taux d'utilisation des serveurs obtenus théoriquement
et par simulation est excellente. Nous n'avons que l'approximation

<figure><img src="image_2.png" /></figure>
faite pour la résolution théorique est très bonne dans ce cas. A partir
de là, nous proposons pour l'étude des probabilités d'état des files d'at-
tente une méthode d'intégration des équations théoriques, relatives à un
processus continu, de façon à retourner au processus discret donné au départ.

Le troisième chapitre montre une utilisation possible des mesures
dans un système informatique pour le pilotage du système. Nous présentons
d'abord le modèle du système objet de l'étude, représentation d'un système
paginé à partage de temps. Dans ce modèle, le comportement des programmes
soumis au système a été représenté par des lois tirées de mesures antérieures
prises sur des programmes réels.

Le comportement du modèle est étudié systématiquement dans ses
différentes plages de fonctionnement, ce qui met en évidence le phénomène
d'écroulement bien connu dans les systèmes (à mémoire virtuelle) et intro-
duit la nécessité d'une régulation.

Nous exposons ensuite un algorithme de régulation directe du
degré de multiprogrammation à partir de courbes de performances évaluées
en temps réel. Cette régulation met en œuvre une technique d'estimation
avec oubli (fondée sur le lissage exponentiel) des différents points des
courbes de performance. Une détermination de précision statistique est
requise pour le filtrage des fluctuation aléatoires. La régulation comprend
en outre un algorithme de recherche dynamique du maximum des courbes.

La mise en œuvre d'une telle régulation demande la définition
d'un bon critère de mesure de performances. Après avoir passé en revue et
écarté divers critères classiques, nous en définissons un nouveau, chiffrant
directement le gain apporté par la multiprogrammation. Cette partie de
l'étude est conclue par la recherche par simulation du meilleur réglage de
paramètres à introduire pour la régulation, puis par la présentation des
améliorations de performances apportées.
TABLE DES MATIERES

Pages

INTRODUCTION 1

CHAPITRE I : LA SIMULATION DES MODELES DE SYSTEMES INFORMATIQUES

PROBLEME DES MESURES ET DE LEUR PRECISION 3

I - GENERALITES SUR LA SIMULATION DE SYSTEMES INFORMATIQUES 3

1 - Généralités sur la simulation. 3

2 - Langages de simulation utilisables. 4

3 - Nature des résultats de simulation. 4

II - DISCUSSION DE LA QUALITE DES MESURES DANS UN SIMULATEUR 5

1 - Vérification de la bonne représentativité d'un modèle. 5

2 - Détermination d'une précision statistique sur les résultats. 6

2.1 - Première construction d'un intervalle de confiance. 7

2.2 - Intervalles de confiance sur des moyennes. 8

2.3 - Techniques de détermination d'intervalle de confiance en simulation. 9

2.4 - Description de la mise en oeuvre de la méthode du blocage des mesures. 13

3 - Longueur d'une simulation pour obtenir une précision donnée. 16

III - CONCLUSION 17

CHAPITRE II : EXEMPLE D'UNE VALIDATION SIMULATEUR D'UN MODELE MATHE-

MATIQUE ET D'UNE SIMULATION 18

I - MODELISATION ADOPTÉE POUR LES SYSTEMES INFORMATIQUES 18
1 - Structure et classification des modèles. 18

2 - Description des techniques de simulation utilisées. 19

## II - LE MODELE ETUDIE ET SA RESOLUTION MATHEMATIQUE 20

1 - Description du modèle. 20

2 - Résolution mathématique du modèle. 20

2.1 - Equations du modèle. 22

2.2 - Résolution pratique du modèle. 25

## III - RESOLUTION DU MODELE PAR SIMULATION 26

1 - Problème des mesures. 26

1.1 - Détermination des grandeurs mesurées. 27

1.2 - Détermination d'une précision sur les teaux d'utilisation des serveurs. 27

1.3 - Mise en œuvre pratique du calcul de précision. 28

? - Résultats obtenus par simulation et comparaison avec la théorie. 29

2.1 - Validation réciproque des deux résolutions. 29

2.2 - Etude de la file d'attente de l'unité centrale. 31

## VI - CONCLUSION 38

# CHAPITRE III : UTILISATION D'UN MODELE DE SIMULATION POUR L'ETUDE DE LA REGULATION ADAPTATIVE D'UN SYSTEME INFORMATIQUE 44

## I - INTRODUCTION 44

## II - PRESENTATION DU MODELE DE SYSTEME ETUDIE 45

1 - Description du modèle. 45

1.1 - Mécanisme des éléments du système. 45
<table><thead><tr><th></th><th>Pages</th></tr></thead><tbody><tr><td>1.2 - Comportement des programmes à l'unité centrale.</td><td>47</td></tr><tr><td>1.3 - Hypothèses sur les lois du modèle.</td><td>48</td></tr><tr><td>2 - Mesures prises au cours de la simulation.</td><td>49</td></tr><tr><td>3 - Validation du simulateur.</td><td>50</td></tr><tr><td>4 - Comportement du modèle.</td><td>54</td></tr><tr><td>III - ETUDE D'UN PRINCIPE DE REGULATION DU SYSTEME</td><td>61</td></tr><tr><td>1 - Présentation du problème de la régulation.</td><td>61</td></tr><tr><td>2 - Description d'ensemble de la régulation adaptative proposée.</td><td>63</td></tr><tr><td>3 - Description de l'estimateur EST.</td><td>65</td></tr><tr><td>4 - Description de l'optimiseur OPT.</td><td>71</td></tr><tr><td>4.1 - Comparaison de deux points.</td><td>71</td></tr><tr><td>4.2 - Algorithme d'estimation du degré optimal de de multiprogrammation.</td><td>73</td></tr><tr><td>5 - Conclusion.</td><td>74</td></tr><tr><td>IV - MISE EN OEUVRE ET ETUDE DU COMPORTEMENT DE LA REGULATION</td><td>75</td></tr><tr><td>1 - Choix d'un critère d'optimisation.</td><td>75</td></tr><tr><td>2 - Etude de la régulation par simulation.</td><td>81</td></tr><tr><td>2.1 - Hypothèses de fonctionnement.</td><td>81</td></tr><tr><td>2.2 - Etude des temps d'adaptation. Réglage des paramètres.</td><td>82</td></tr><tr><td>2.3 - Performances obtenues grâce à la régulation.</td><td>88</td></tr><tr><td>V - CONCLUSION</td><td>89</td></tr><tr><td>CONCLUSION</td><td>92</td></tr><tr><td>BIBLIOGRAPHIE</td><td>95</td></tr><tr><td>ANNEXES</td><td>99</td></tr></tbody></table>
# INTRODUCTION

Le comportement d'un système informatique peut s'étudier par
deux approches :

- Mesures dans un système réel soumis à l'exploitation,

- Représentation du système par un modèle rendant bien compte de
la réalité.

La première approche peut se conduire à l'aide de sondes câblées
[ 28 ] (sans modification ni perturbation du système dans ce cas, mais
toutes les mesures ne sont pas possibles) ou programmées [ 34 ]. Dans ce
dernier cas, des modifications du système sont nécessaires et celui-ci
peut-être perturbé. Dans les deux cas, une bonne connaissance du système
est nécessaire. Les essais de nouveaux algorithmes demandent un travail de
programmation système important. Enfin l'étude d'un système inexistant à
l'état d'avant-projet est impossible par cette méthode.

La seconde approche demande la construction d'un modèle représen-
tant bien la réalité. Le comportement de tels modèles doit être comparé,
directement ou indirectement, à celui d'un système réel. Le modèle construit
peut être résolu par des méthodes analytiques: au prix d'hypothèses simpli-
ficatrices, le problème peut se traduire par des équations que l'on sait
résoudre. Un cas particulier important est celui des modèles probabilistes
à serveurs et à files d'attente.

Les principales techniques de résolution de ces modèles sont :

- Utilisation des chaînes de Markov. Jackson [ 32 ] a résolu
le problème général d'un système markovien constitué d'un graphe de serveurs
et de files d'attente.

- Approximation à l'aide de processus de diffusion [29,33] dont le
principe sera brièvement rappelé au chapitre II. Les méthodes de résolution
analytique imposent des hypothèses restrictives sur le fonctionnement du
modèle. Elles sont peu coûteuses pour le calcul des solutions et décrivent
les variations du phénomène étudié par des équations. Mais certaines études
(mécanismes compliqués, comportements transitoires par exemple) ne peuvent
être effectués par de telles méthodes faute de savoir résoudre les équations,
ou même de savoir construire le modèle...

Le modèle construit peut se résoudre également par simulation, seule méthode dont il sera question dans cette étude. La simulation consiste à reproduire le fonctionnement du système pour en étudier des comportements possibles dans le temps. Le but d'une simulation peut être :

- De tester la correction fonctionnelle d'un mécanisme [39]

- D'évaluer des performances, seul cas qui nous intéressera.

Nous décrions deux expériences d'évaluation de performances par simulation. Dans l'une d'elles, toute résolution théorique du modèle est impossible, à la fois en raison de la prise en compte de certains mécanismes (dépassements, boucle d'asservissement) et parce que les phénomènes étudiés sont variables dans le temps. Nous étudierons un tel cas au chapitre III avec l'étude de l'adaptation d'un système à la variation de sa charge dans le temps. Dans l'autre, la simulation est utilisée pour valider une résolution mathématique toutes les fois que celui-ci comporte des approximations. Le chapitre II donnera l'exemple d'une telle étude.

Contrairement aux résolutions mathématiques, une simulation pour évaluer des performances pose le problème de la qualité, notamment de la précision, des résultats. En effet, ceux-ci sont des statistiques dont il faut apprécier le caractère significatif. Ce sera notamment le cas dans les deux exemples étudiés, d'une part pour permettre la comparaison de deux méthodes, d'autre part pour prendre une décision dans une boule d'asservissement avec une confiance statistique déterminée.

Nous présentons successivement d'autres problèmes généraux de mesures dans les simulations dans un premier chapitre et nous y développons plus particulièrement le problème de la précision statistique. Dans un second chapitre, la simulation nous servira de comparaison avec une méthode de résolution mathématique approchée. Le troisième chapitre donne un exemple d'utilisation des mesures dans un système pour adapter celui-ci à des variations de sa charge. La simulation y aboutit au test d'un algorithme susceptible d'être employé dans un système réel.
# CHAPITRE I

## LA SIMULATION DES MODELES DE SYSTEMES INFORMATIQUES

### PROBLEMES DES MESURES ET DE LEUR PRECISION

#### I - GENERALITES SUR LA SIMULATION DES SYSTEMES INFORMATIQUES

##### 1 - Généralités sur la simulation

La simulation est la reproduction du fonctionnement du système physique étudié, existant ou non, représenté par un modèle, dans le temps. Le fonctionnement des diverses parties du système peut être déterministe ou aléatoire. Dans ce dernier cas, il est reproduit par des procédures de tirage au sort au moyen de nombres au hasard. On parle de simulation de Monte-Carlo.

L'évolution d'un système physique dans le temps peut s'effectuer [ 35 ] d'une manière continue ou discrète. Dans le premier cas, le fonctionnement se décrit par des systèmes d'équations différentielles ou aux dérivées partielles résolues par des méthodes d'intégration numérique. Le temps varie d'une façon continue ou plutôt, pour les besoins du calcul, par petits intervalles de temps élémentaires. Nous ne nous intéresserons pas à ce cas.

Les systèmes auxquels nous nous intéresserons sont tous à évolution discrète dans le temps. Entre deux instants T0 et T1, il existe une suite finie d'instants $t_1, t_2, ..., t_i, ..., t_n$ tels qu'entre deux instants $t_{i-1}$ et $t_i$, le système n'évolue pas. Les instants $t_i$ de changement d'état du système seront appelés _occurences d'évènements_ et les changements d'états seront les _évènements_. Les évènements arrivant d'une manière discrète, il est possible de les ordonnancer dans le temps suivant un échéancier et de les simuler séquentiellement l'un après l'autre. Les simulations de systèmes informatiques présentées ci-après sont toutes à évènements discrets.
## 2 - Langages de simulations utilisables

Tous les langages algorithmiques généraux sont utilisables pour la simulation du fait de la possibilité d'exécution séquentielle. Il existe cependant différents langages disposant d'un certain nombre des outils nécessaires [37]. Citons notamment DYNAMO de Forrester [21] adapté à la simulation de phénomènes continus, et les langages adaptés à la simulation par évènements discrets. Ceux-ci comprennent notamment GPSS [31] diffusé par IBM, SIMSCRIPT [17, 38] qui est une extension de Fortran et SIMULA 67, choisi par nous, qui est une extension d'ALGOL 60 comprenant en outre les notions de classes d'objets et de coroutines. Les objets sur lesquels portent les évènements sont des processus. Une description sommaire de SIMULA 67 est donnée en annexe I. ASPOL [36] est un langage proche de SIMULA 67 dont il a la notion de processus, les synchronisations s'y opérant par des sémaphores.

## 3 - Nature des résultats de simulation

Le but d'une simulation peut être, soit de vérifier la correction fonctionnelle du système étudié ou de l'un de ses sous-ensembles, soit d'évaluer des performances. Dans le premier cas, l'on se borne à reproduire fidèlement le fonctionnement étudié [39] (exemple simulation d'une machine non existante ou d'une logique câblée sur lesquels des tests doivent être effectués), ou à mettre en évidence des phénomènes gênants tels l'interblocage dans les systèmes.

Nous ne nous intéresserons qu'aux évaluations de performances. Diverses grandeurs sont prélevées périodiquement dans le système, soit sur évènement (entrée ou sortie de client, début et fin d'utilisation d'une ressource, en particulier d'une unité centrale ou périphérique), soit sur top d'horloge. Les grandeurs prélevées sont considérées comme des variables aléatoires sur lesquelles nous effectuerons des statistiques. Les résultats donnés le plus couramment sont : moyennes, écarts-types, probabilités simples ou conditionnelles, parfois histogrammes ou corrélations. Nous mesurons essentiellement les taux d'utilisation des différentes unités, les tailles des files d'attente, les temps de réponse des usagers dans un système conversationnel (un temps de réponse est défini comme la différence des heures de sortie et d'entrée de l'usager considéré dans le système).
Les grandeurs mesurées peuvent se répartir en deux catégories : Les unes, telles que les temps de réponse des usagers (ou les différentes mesures globales particulières à chacun) ne peuvent être échantillonnées que lors de certains évènements précis dans le système (ici une sortie d'usager). Les mesures peuvent être assez rares si l'évènement considéré est peu fréquent. De telles mesures serviront à comparer qualitativement deux systèmes entre eux, ou deux réglages différents de paramètres. D'autres grandeurs, telles que les taux d'utilisation des différentes unités, les probabilités de trouver le système dans un certain état, peuvent se calculer entre deux instants $t_1$ et $t_2$ quelconques. De telles mesures peuvent être prélevées aussi souvent que l'on veut, sur top d'horloge par exemple, et sont calculées généralement par sommation de différents intervalles de temps. Elle permettent elles aussi des comparaisons qualitatives comme les précédentes, mais elles sont également utilisables pour le pilotage du système lui-même comme nous le verrons au chapitre III.

## II - DISCUSSION DE LA QUALITE DES MESURES DANS UN SIMULATEUR

Les mesures de performances effectuées dans une simulation peuvent être peu satisfaisantes pour plusieurs raisons [ 35 ].

Tout d'abord, le modèle proposé ne correspond pas suffisamment à la réalité, soit à cause d'un mauvais réglage des paramètres, soit à la suite d'hypothèses simplificatrices trop peu réalistes. Ces deux inconvénients peuvent exister pour tout modèle, qu'il soit résolu par simulation ou mathématiquement.

D'autre part, les résultats de simulation peuvent n'être pas statistiquement significatifs et nous serons donc amenés à calculer une précision statistique.

### 1 - Vérification de la bonne représentativité d'un modèle

Cette vérification se fait dans des cas particuliers par comparaison avec une autre méthode d'évaluation : mesures dans un système réel ou résultats d'un autre modèle déjà contrôlé, par exemple un modèle mathématique déjà bien confronté à la réalité. Les modèles simulés dans cette étude seront tous comparés à des modèles mathématiques, soit pour donner une vérification supplémentaire de ceux-ci, soit pour contrôler la simulation elle-même.
La comparaison d'un modèle avec une autre méthode se limite le plus souvent à la comparaison visuelle de tableaux de chiffres ou de courbes dans le plus grand nombre de cas possibles. Une comparaison plus rigoureuse comprend deux opérations [ 35 ] [ 2, 3 ] :

- L'étalonnage, c'est-à-dire le réglage des paramètres du simulateur jusqu'à l'accord des résultats obtenus dans quelques cas particuliers.

- La validation, qui comprend aussi le test de la robustesse du modèle, c'est-à-dire la détermination du domaine dans lequel on peut faire confiance aux résultats obtenus.

Dans le cas de la comparaison avec un modèle mathématique, ces opérations se réduisent à assurer la coïncidence des paramètres respectifs et à comparer les résultats obtenus. Dans le cas de systèmes réels, des méthodes statistiques ont été définies [ 3 ] pour l'étalonnage et la validation de simulateurs alimentés par une trace de mesures. Les mesures extraites de la trace sont comparés à celles extraites du simulateur sous forme de régressions multilinéaires. Une première partie de la trace sert à régler les paramètres du modèle et une autre à vérifier que le modèle étalonné se comporte convenablement sur un ensemble de points différents.

Nous nous limiterons dans cette étude à la comparaison avec des modèles mathématiques.

## 2 - Détermination d'une précision statistique sur les résultats

Un modèle de simulation étant correct et bien conforme à la réalité, les résultats peuvent être mauvais s'ils sont peu significatifs statistiquement. Les résultats obtenus sont en effet des estimations statistiques des grandeurs exactes recherchées et ces estimations sont elles-même des variables aléatoires. Le problème de la précision statistique d'une estimation conduit à déterminer, pour la valeur théorique exacte recherchée, un intervalle de confiance, c'est-à-dire à affirmer avec une certaine probabilité maximum d'erreur (en général 5 %) que la grandeur recherchée se trouve dans un certain intervalle (voir annexe II).

Nous nous limiterons à des études de précision sur les espérances mathématiques et nous nous ramènerons toujours à ce cas.
## 2.1. - Première construction d'un intervalle de confiance

Soit un échantillon $u_1, \dots, u_n$ de n observations d'un phénomène aléatoire et X une estimation sans biais (voir annexe II) d'une grandeur théorique G à partir de l'échantillon. X est elle-même une variable aléatoire et $E(X) = G$ par définition. Soit $\sigma_X$ l'écart-type de l'estimation X. L'inégalité de Bienaymé-Tchebycheff donne un premier intervalle de confiance sur G. En effet :

$$ (I-1) \quad \text{Prob} \left[ |X - E(X)| < \delta \sigma_X \right] > 1 - \frac{1}{\delta^2} \quad \text{quel que soit } \delta > 0. $$

Comme $E(X) = G$ par hypothèse, nous aurons bien, avec une confiance d'au moins $1 - \frac{1}{\delta^2}$

$$ (I-2) \quad X - \delta \sigma_X < G < X + \delta \sigma_X $$

Pour une confiance habituelle de 95 %, nous aurons $1 - \frac{1}{\delta^2} = 0,95$ soit $\delta = \sqrt{20} \approx 4,47$ soit $\delta = 4,5$ plus simplement.

L'inégalité de Bienaymé-Tchebycheff ne fait aucune hypothèse sur la loi suivie par l'estimation X et il sera souvent possible de déterminer un intervalle plus serré.

Notamment si notre estimation X suit une loi gaussienne (cas de la moyenne d'un nombre suffisant d'observations indépendantesd'après le théorème Central-limite) les tables de la loi de Gauss donnent pour une confiance de 95 % :

$$ (I-3) \quad X - 1,96 \sigma_X < G < X + 1,96 \sigma_X $$

D'une façon générale, la détermination de la précision statistique d'une estimation se ramènera à la détermination d'une variance (donc un écart-type) sur cette estimation.

Un cas particulièrement intéressant est celui d'un estimateur convergent [25, 27]. Dans ce cas, nous pourrons faire tendre vers 0 l'amplitude de l'intervalle déterminé en prolongeant la simulation, autrement dit nous parviendrons à atteindre une précision fixée à l'avance. Ce ne sera pourtant pas toujours le cas comme nous le verrons dans le chapitre III.
## 2.2 - Intervalles de confiance sur des moyennes

Soit un échantillon $X_1, \dots, X_n$ de n observations d'un même phénomène. Les n observations sont supposées identiquement distribuées et mutuellement indépendantes, ou plus simplement et d'une façon suffisante, de même moyenne et variance $E(X)$ et $V(X)$ et non corrélées deux à deux.

L'estimateur le plus usuel de $E(X)$ est la moyenne arithmétique $\bar{X}$ de l'échantillon (estimateur de variance minimale)

$$ (I-4) \quad \bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_i $$

Cet estimateur est sans biais car :

$$ E(\bar{X}) = \frac{1}{n} \sum_{i=1}^{n} V(X_i) = E(X) $$

Les observations $X_i$ étant non corrélées nous avons (annexe 2)

$$ (I-5) \quad V(\bar{X}) = \frac{1}{n^2} \sum_{i=1}^{n} V(X_i) = \frac{1}{n} V(X) $$

d'où si $\bar{\sigma}_X = \sqrt{V(\bar{X})}$ et $\sigma_X = \sqrt{V(X)}$ sont les écarts-types sur $\bar{X}$ et sur les $X_i$ :

$$ (I-6) \quad \sigma_{\bar{X}} = \frac{\sigma_X}{\sqrt{n}} $$

$V(X)$ s'estime par l'estimateur sans biais classique [ 25 ] :

$$ (I-7) \quad V(X) = \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] $$

valable si les observations $X_i$ ne sont pas corrélées entre elles.

La moyenne d'un nombre suffisant d'observations indépendantes étant approximativement gaussienne (théorème central-limite), la formule (I-3) s'applique et devient

$$ (I-8) \quad \bar{X} - \frac{1,96 \sigma_X}{\sqrt{n}} < E(X) < \bar{X} + \frac{1,96 \sigma_X}{\sqrt{n}} $$
Si l'effectif n de l'échantillon devient suffisamment grand, il est théoriquement possible d'atteindre toute précision fixée à l'avance, en raison du $\sqrt{n}$ au dénominateur.

Si les observations de l'échantillon sont corrélées entre elles, les considérations ci-dessus cessent d'être valables. D'après une formule générale (voir annexe 2) :

$$ (I-9) \quad \sum_{i=1}^{n} V(\sum_{i=1}^{n} X_i) = \sum_{i=1}^{n} V(X_i) + \sum_{i=1}^{n} \sum_{j=1, i \neq j}^{n} \text{Cov}(X_i, X_j) $$

et la formule (5) se transforme en :

$$ (I-10) \quad V(\bar{X}) = \frac{1}{n} V(X) + \frac{1}{n^2} \sum_{\substack{i=1 \\ i \neq j}}^{n} \sum_{i=1}^{n} \text{Cov}(X_i, X_j) $$

Si les observations $X_i$ sont fortement corrélées entre elles, ce qui est le cas en simulation pour les mesures proches les unes des autres dans le temps, le dernier terme devient prépondérant et l'application bruitale de la formule (I-8) risque d'être très optimiste dans le cas où les corrélances sont positives. Nous allons passer en revue diverses techniques possibles pour tourner cette difficulté.

## 2.3 - Techniques de détermination d'intervalle de confiance en simulation.

Pour déterminer en simulation une précision sur une grandeur déterminée, nous ferons en général apparaître cette grandeur sous la forme d'une espérance mathématique. Nous estimerons la grandeur par une moyenne et nous appliquerons la formule (I-8) (ou nous nous y ramènerons). Le problème est donc d'obtenir pour la grandeur recherchée n mesures de même espérance, même variance et non corrélées deux à deux.

### a) Première technique : faire n simulations

Il est théoriquement possible de faire n passages différents du simulateur, avec les mêmes paramètres et des données différentes mais statistiquement semblables (générées avec des suites de nombres pseudo-aléatoires différentes dans le cas de simulateurs probabilistes). Chaque simulation donne une
mesure et une seule de la grandeur étudiée. Les n simulations donnent n mesures identiquement distribuées et mutuellement indépendantes. La formule (I-8) s'applique bien. Aucune hypothèse n'est à faire sur l'évolution de la grandeur étudiée dans le temps. Mais cette méthode trop coûteuse est irréaliste. Elle conduit en effet à un très grand nombre d'expériences pour évaluer une précision. Aucune détermination de celle-ci n'est d'autre part possible en temps réel.

Les techniques utilisées dans la pratique consistent à effectuer les n mesures $X_i$ nécessaires au cours d'une seule simulation plus longue. Mais des hypothèses sont nécessaires sur la grandeur étudiée, que nous représenterons par un processus stochastique $\chi = \{X(t)\}$.

Tout d'abord, les mesures $X_i$ doivent être toutes de même espérance mathématique et de même variance. D'où la nécessité d'hypothèses de stationnarité sur le processus $\chi$.

**Définitions** : [11,27,22] – $\chi$ est un processus strictement stationnaire si la distribution des variables aléatoires $X(t_1)$, $X(t_2)$, ..., $X(t_m)$ ne dépend pas d'une translation dans le temps, c'est-à-dire si leur distribution est identique à celle de $X(t_1+h)$, $X(t_2+h)$, ..., $X(t_m+h)$. La distribution de $X(t)$ ne dépend donc pas de t.

$\_\chi$ est un processus faiblement stationnaire si

1. L'espérance de $X(t)$ est constante dans le temps soit

$$ (I-11) \quad E[X(t)] = E(X) $$

2. L'autocovariance de $X(t)$ sur un intervalle $h$ c'est-à-dire $\text{Cov}[X(t), X(t+h)]$ ne dépend pas d'une translation dans le temps autrement dit :

(I-12) $\text{Cov}[X(t), X(t+h)] = \gamma(h)$ fonction de h seul
$\gamma(h)$ est appelée fonction d'autocovariance du processus.
En particulier $V[X(t)] = \gamma(0)$ est indépendante du temps soit

$$ (I-13) \qquad V[X(t)] = V(X) $$

La forte stationnarité d'un processus est difficile à tester à moins d'hypothèses théoriques. La stationnarité faible peut par contre se tester et permet d'affirmer la stabilité dans le temps des deux premiers moments.

Nous n'étudierons en simulation que des processus faiblement stationnaires dans le temps ou du moins pendant des intervalles de temps aussi longs. Dans le cas d'un processus non stationnaire, nous devrons supposer que les caractéristiques du processus varient dans le temps d'une façon lente, ce qui nous permettra de négliger ces variations. Nous verrons un exemple de cette technique au chapitre III. Il est impossible en particulier, sauf en effectuant n simulations, d'étudier quantitativement des phénomènes transitoires rapidement variables.

Nous ferons d'autre part l'hypothèse, bien vérifiée dans la pratique, que des mesures de plus en plus éloignées dans le temps sont de moins en moins corrélées. Autrement dit, si $\gamma(h)$ est la fonction d'auto-covariance définie par (I-12)

$$ (I-14) \qquad \gamma(h) \to 0 \text{ si } h \to \infty $$

Sous ces hypothèses, nous pouvons déterminer entièrement la moyenne d'une grandeur avec une précision déterminée, tout en n'effectuant qu'une seule simulation. Nous effectuerons n mesures $X_1, \dots, X_n$ au cours du temps et nous pourrons affirmer que $E(X_n) = E(X)$ et $V(X_n) = V(X)$ indépendants de la mesure. Nous aurons en outre :

$$ (I-15) \qquad \operatorname{Cov}(X_{i+h}, X_i) \to 0 \text{ si } h \to \infty $$
D'où les méthodes :

b) Négliger l'autocorrélation des mesures

Nous supposons les n mesures peu corrélées entre elles et
nous appliquons brutalement la formule (I-8) . Cette méthode
sommaire risque de donner des résultats optimistes, mais
peut fournir des ordres de grandeur.

c) Faire une hypothèse sur la décroissance de la fonction d'autocovariance $\gamma(h)$ [8, 10]

Si la forme de la fonction $\gamma(h)$ est connue ou si une majoration
est possible, nous pouvons appliquer la formule (I-10) au lieu
de la formule (I-5). D'où une variance de la moyenne $\bar{X}$ des me-
sures. Avec l'hypothèse que $\bar{X}$ suit une loi gaussienne la for-
mule (3) devient ici :

$$
(I-16) \quad \bar{X} - 1,96 \sqrt{V(\bar{X})} < E(X) < \bar{X} + 1,96 \sqrt{V(\bar{X})}
$$

où d'après (I-10), $V(\bar{X})$ tendra vers 0 si le nombre n de mesures augmente.

Fishman [22, 23, 24] a relié l'étude de la fonction d'autocor-
rélation $\rho(h)$ définie par

$$
(I-17) \quad \rho(h) = \frac{\gamma(h)}{V(X)}
$$

aux techniques de l'analyse harmonique.

d) Se ramener aux cas des mesures non corrélées

L'échantillon total des n mesures est décompé en p blocs de
mesures consécutives. Les p blocs constitueront un nouvel
échantillon non autocorrélé. Cela est possible de deux façons :

- En séparant les blocs par des points de régénération [12,23]. Un état du système définit un point de régénération si l'évolution future du système ne dépend plus du passé. Par exemple
dans un système aux arrivées de clients poissonniennes, les instants où le système est vide sont des points de régénération. Les blocs forment en quelque sorte p sous-simulations indépendantes.

- En prenant les p blocs égaux et de taille K suffisamment grande pour que l'autocorrélation soit négligeable d'un bloc à l'autre [9]. C'est la méthode que nous avons retenu dans la pratique.

## 2.4 - Description de la mise en œuvre de la méthode du blocage des mesures

Soit à répartir les n mesures prises au cours de la simulation en p blocs égaux de taille K suffisamment grande pour que les résultats sur chaque bloc (nous prendrons des moyennes partielles comme précisé plus loin) soient non corrélés entre eux. Les opérations complètes à effectuer sont :

- Détermination d'une taille minimum de bloc soit $K_0$. Cette détermination demande un premier passage du simulateur.

- Nouveau passage du simulateur, où une taille de bloc $K \ge K_0$ est donnée en paramètre. La taille de bloc $K_0$ doit avoir été choisie pour être valable d'une simulation à l'autre, avec des paramètres différents.

Sur chacun des p blocs, nous effectuerons la moyenne des mesures sur le bloc, ou moyenne partielle, que nous considérons comme la mesure sur l'ensemble du bloc. Nous effectuerons sur les moyennes partielles un nouveau calcul de moyenne, variance et intervalle de confiance.

### 2.4.1 - Détermination d'une taille minimum $K_0$ de bloc

Nous effectuerons un premier passage du simulateur dont nous extrairons l'échantillon complet de toutes les mesures (par exemple sur une bande magnétique). Nous traitons ensuite de nouveau l'échantillon et nous estimons la fonction d'autocovariance de l'échantillon des n mesures. Par définition pour un échantillon discret :

$$ (I-18) \quad \gamma(h) = \text{Cov}(Y_i, X_{i+h}) $$
où $X_i$ est la ième mesure de la mesure. En raison de la stationnarité faible nous considérons encore cette nouvelle grandeur comme fonction de h seul.

Soit n le nombre total de mesures. Nous estimerons $\gamma(h)$ par l'estimateur :

$$ (I-19) \quad \hat{\gamma}(h) = \frac{1}{n-h} \sum_{i=1}^{n-h} X_{i,h} \cdot X_{i+h} - \frac{1}{(n-h)^2} \left( \sum_{i=1}^{n-h} X_{i,h} \cdot \sum_{i=1}^{n-h} X_{i+h} \right) \quad (\text{avec } h \le n) $$
estimateur classique d'une covariance. Dans (I-19) h doit être petit devant n pour que l'estimation soit valable.

Nous estimerons successivement ainsi $\gamma(0), \gamma(1), \dots, \gamma(i)$ et nous choisirons $K_0$ pour que $\gamma(K_0)$ soit négligeable à côté de la variance $\gamma(0)$.

La taille minimum de bloc étant déterminée, nous la majorerons éventuellement par une valeur K susceptible de donner des résultats satisfaisant dans des simulations effectuées avec d'autres paramètres.

### 2.4.2 - Mise en oeuvre du blocage proprement dit

L'intérêt de déterminer préalablement une taille de bloc convenable K, valable d'une simulation à l'autre, est d'utiliser cette valeur pour différents passages du simulateur. D'autre part, la méthode pourra dorénavant être séquentielle, d'où la possibilité d'effectuer les calculs en temps réel.

Soient donc p blocs successifs de K mesures consécutives, tels que deux mesures appartenant à un bloc différent soient non corrélées. Nous noterons $X_{i,j}$ la jième mesure du bloc i, n le nombre total de mesures, $\bar{X}$ la moyenne des n mesures, $\bar{X}_i$ la moyenne (dite partielle) des mesures sur le ième bloc soit

$$ (I-20) \quad \bar{X}_i = \frac{1}{K} \sum_{j=1}^{K} X_{ij} $$
Les $\bar{X}_i$ sont non corrélés entre eux par construction des blocs. Leur moyenne est :

$$ \frac{1}{p} \sum_{i=1}^{p} \frac{1}{K} \sum_{j=1}^{K} X_{ij} = \frac{1}{pK} \sum_{i=1}^{p} \sum_{j=1}^{K} X_{ij} = \frac{1}{pK} \sum_{i=1}^{p} \sum_{j=1}^{K} = \bar{X} $$

La moyenne des $\bar{X}_i$ est bien la moyenne des mesures initiales. Nous calculerons donc successivement chaque $\bar{X}_i$, puis leur moyenne $\bar{X}$ et leur variance $\nu$ d'où leur écart-type $\sigma$. Soit $E(X)$ l'espérance mathématique des $X_{ij}$ (ou aussi des $m_i$). Sur le nouvel échantillon constitué par les $m_i$, la formule (I-8) s'applique (toutes ses hypothèses sont vérifiées) d'où :

$$ (I-21) \quad \bar{X} - \frac{1.96\sigma}{\sqrt{p}} < E(X) < \bar{X} + \frac{1.96\sigma}{\sqrt{p}} $$

avec une confiance statistique de 95 %.

**Remarque :** Il est de peu d'importance pour le calcul de la précision obtenue de majorer la taille des blocs. La réduction résultante du nombre $p$ de blocs sera compensée par une diminution parallèle de $\sigma$. Montrons le à titre d'exemple si $K$ est changé en $2K$, l'effectif total n des mesures restant inchangé :

Dans le premier cas, le nombre de blocs est $p$, dans la seconde, il est $p/2$. Chaque bloc dans le second cas résulte de la conaténation de deux blocs obtenus dans le premier cas. Soit $\bar{X}_i'$ la moyenne partielle des mesures sur le ième bloc dans le second cas, les moyennes partielles dans le premier cas étant notées $\bar{X}_1, \bar{X}_2, ..., \bar{X}_p$. Soit $\nu$ la variance des moyennes partielles $\bar{X}_i$ dans le premier cas, $\nu'$ la variance des $\bar{X}_i$ dans le second cas, $\sigma$ et $\sigma'$ les écarts-types. Les mesures seront toujours indicées comme pour le premier cas. Nous avons :

$$ \bar{X}_i' = \frac{1}{2K} \left[ \sum_{j=1}^{K} X_{2i,j} + \sum_{j=1}^{K} X_{2i+1,j} \right] \quad \text{soit} $$

$$ \bar{X}_i' = \frac{1}{2} (\bar{X}_{2i} + \bar{X}_{2i+1}) $$
d'où en raison de l'absence de corrélation :

$$ \gamma' = V(\bar{X}_i') = \frac{1}{2} V(\bar{X}_i) = \frac{1}{2} \gamma $$

$$ d'où \sigma' = \frac{\sigma}{\sqrt{2}} . $$

La demi-amplitude de l'intervalle de confiance déterminé vaut $\frac{1,96\sigma}{\sqrt{p}}$ dans le premier cas et $\frac{1,96\sigma}{\frac{\sqrt{p}}{2}}$ (nous avons maintenant $\frac{N}{2}$ blocs) dans le deuxième cas, soit $\frac{1,96\sigma}{\sqrt{2}\frac{\sqrt{p}}{2}} = \frac{1,96\sigma}{\sqrt{p}}$ de même que dans le premier cas.

Cette remarque montre qu'une taille de bloc grande n'est pas un inconvénient. Cette possibilité de majoration permet :

- D'utiliser une même taille de bloc pour des simulations assez différentes.

- D'éviter dans certains cas la première étape (détermination de la taille de blocs nécessaire) dans les simulations longues.

Nous choisirons à priori pour K une très grande valeur, de façon à disposer encore d'un nombre de blocs compris entre 50 et 100, suffisant pour que l'estimation d'une variance entre blocs soit possible.

### 3 - Longueur d'une simulation pour obtenir une précision donnée

Parmi les conséquences des méthodes exposées plus haut, il est possible de donner un ordre de grandeur de la durée de fonctionnement du modèle à simuler ou de la dépense de temps machine, pour obtenir une précision donnée sur une moyenne.

Nous pouvons faire les hypothèses suivantes, vérifiées approximativement dans la pratique :

- le nombre de mesures effectuées dans une simulation est à peu près proportionnel au temps simulé.
- Le temps machine dépensé est à peu près proportionnel lui aussi au temps simulé (il dépend en fait du nombre d'évènements générés, aux opérations initiales près).

Nous avons vu (formules (I-8), (I-21) que la précision obtenue est inversement proportionnelle, soit à la racine carrée du nombre de mesures, soit du nombre de blocs indépendants, donc dans ce cas encore à la racine carrée du nombre total de mesures. Donc :

La précision obtenue sur une moyenne est à peu près inversement proportionnelle à la racine carrée du temps simulé, ou à la racine carrée du temps machine dépensé pour l'exécution du programme de simulation.

Divers calculs de précision (ainsi que la comparaison visuelle avec des résultats obtenus par des méthodes mathématiques) pour différents temps simulés nous ont montré que cette loi était à peu près vérifiée.

Nous en déduisons le coût important d'une simulation destinée à obtenir une grande précision. On se contente souvent d'une précision relative de 10 à 15 % pour un simulateur. Certaines applications (voir Chapitre II) imposent cependant des précisions inférieures à 1 % (sur des grandeurs peu dispersées cependant). Une bonne précision courante est de l'ordre de quelques pour cent.

### III - CONCLUSION

Les techniques de calcul de précision décrites ci-dessus pour des mesures dans un phénomène stochastique autocorrélé sont valables pour tout échantillon de mesures d'un phénomène étudié au cours du temps. Elles peuvent être mises en oeuvre hors de tout contexte de simulation. Nous allons maintenant mettre en oeuvre ces techniques dans deux applications : la comparaison de deux modèles et le pilotage d'un système simulé en fonction des mesures de ces performances.
# CHAPITRE II

## EXEMPLE D'UNE VALIDATION MUTUELLE D'UN MODELE MATHEMATIQUE ET D'UNE SIMULATION

### I - MODELISATION ADOPTEE POUR LES SYSTEMES INFORMATIQUES

Dans la suite de cette étude, les systèmes informatiques seront représentés par des modèles à serveurs et à files d'attente. Une telle modélisation va nous permettre à chaque fois (au moins dans des cas particuliers) d'entreprendre une résolution à la fois par des méthodes mathématiques et par simulation, ce qui permettra de contrôler l'une par l'autre chacune de ces deux méthodes.

#### 1 - Structure et classification des modèles

Dans les modèles étudiés, les usagers du système sont représentés par des entités appelées clients parcourant un graphe de serveurs et de files d'attente. A chaque serveur est associé sa file et réciproquement. Chaque file d'attente est gérée suivant le principe : premier arrivé-premier servi (l'extension au cas de plusieurs files d'attentes avec des priorités différentes ne présente pas de difficulté supplémentaire). Les serveurs représentent les divers éléments du système tels qu'unité centrale, mémoire secondaire (tambour de pagination) pour un système paginé, périphérique d'entrées-sorties. Le fonctionnement d'un servaur est le suivant : au démarrage du système ou lors du départ d'un client, le serveur prend en compte le premier client contenu dans sa file d'attente et le sort pendant un temps de service constant ou aléatoire suivant une distribution donnée. Si la file est vide, le serveur se met en attente de l'arrivée d'un nouveau client. A la fin du service, le client est introduit dans la file d'attente d'un autre serveur choisi suivant une procédure que nous appellerons aiguillage. Un aiguillage
peut être suivant le cas, soit probabiliste, soit dépendant d'un mécanisme propre au serveur.

On distingue les modèles ouverts où les clients entrent et sortent du système et les modèles fermés contenant un nombre de clients invariable sans arrivée ni départ.

Lorsque toutes les distributions des temps de service sont exponentielles et que tous les aiguillages en fin de service sont probabilistes avec des probabilités données, et si de plus les arrivées éventuelles de nouveaux clients sont poissonniennes, nous retrouvons les modèles résolus mathématiquement par Jackson [ 32 ]. Nous en verrons un exemple au chapitre III.

## 2 - Description des techniques de simulation utilisées

Les simulateurs utilisés dans nos études sont écrits en langage SIMULA 67. Les serveurs sont décrits comme des processus. Les files d'attente sont représentées par des listes. Les clients sont des objets passifs circulant sur le réseau de serveurs et de files d'attente et porteurs de paramètres propres. Ils sont décrits en SIMULA 67 comme des liens. Dans les réseaux ouverts, les clients sont générés par une station particulière du réseau nommée naissance ; la sortie du système, ou mort, est considérée elle aussi comme une station particulière.

La mise à jour des mesures se fait lors des changements d'état du système, notamment arrivée d'un nouveau client, début ou fin de service à un serveur. Le prélèvement des grandeurs resurées se fait, soit globallement à la fin du temps simulé, soit périodiquement sur horloge. (Cas des calculs de précision comme nous le reverrons plus loin).

Nous allons maintenant étudier un premier exemple d'un modèle fermé à deux serveurs récluté à la fois mathématiquement et par simulation. Dans un premier temps, nous indiquerons une résolution mathématique approchée
du modèle. Dans un second temps, une partie des résultats théoriques servira
à contrôler la simulation. Enfin, l'ensemble des résultats de simulation
guidera dans l'interprétation des formules théoriques de façon à obtenir par
celles-ci tous les résultats.

II - LE MODELE ETUDIE ET SA RESOLUTION MATHEMATIQUE

1 - Description du modèle

Le modèle étudié (figure 1) est un modèle fermé à deux serveurs
que nous nommerons par la suite unité centrale et tambour (ou mémoire secon-
daire). Toutes les fonctions de lors de service sont admissibles à chacun des
deux serveurs, mais nous nous limiterons au cas particulier : temps de service
au tambour constant et pris égal à l'unité de temps ; temps de service à
l'unité centrale distribué exponentiellement avec la moyenne $1/μ$. Le modèle
sera étudié en fonction de deux paramètres : le nombre M de clients dans
le système et l'inverse $μ$ du temps de service moyen à l'unité centrale.
Les résultats à obtenir sont dans un premier temps les taux d'utilisation
t à l'unité centrale, et dans un second temps les probabilités $P(n)$ de n
clients dans la file d'attente de l'unité centrale ($0 \le n < M$).

2 - Résolution mathématique du modèle

Nous allons brièvement décrire la résolution approchée du modèle
par la méthode des équations de diffusion [11], [29], [33]. L'idée de
la méthode est l'approximation du processus stochastique discret qu'est le
nombre de clients dans la file d'attente de l'unité centrale par un processus
continu $\{X(t), t \ge 0\}$ avec $0 \le X \le M$, fonction du temps. Les équations
différentielles et aux différences décrivant le modèle dans le cas discret
deviennent dans le cas continu des équations aux dérivées partielles du
type équations aux diffusions de la chaleur [30] résolubles analytiquement
dans des cas plus généraux. La résolution indiquée ici est le rappel et la
reproduction partielle des résultats de [30].
<figure><img src="image_3.png" /><figcaption>Figure 1</figcaption></figure>
2.1 – Équations du modèle

Considérons le processus stochastique $\{X(t), t \ge 0\}$ sur l'intervalle fermé $[0, M]$. Sur l'intervalle ouvert $[0, M]$, le processus étudié est un processus de diffusion [6]. Lorsque l'un des points extrêmes 0 ou M est atteint (file de l'unité centrale ou du tambour vide), le processus y reste pendant un temps distribué exponentiellement jusqu'à une arrivée ou départ d'un client. A ce moment là, il saute instantanément à 1 ou M-1 suivant le cas. Soient $m_1(t)$ et $m_2(t)$ les masses de probabilités concentrées entre 0 et M, $1/\lambda_1$ et $1/\lambda_2$ les temps moyens avant saut à chacun de ces points. Avant de poser les équations, définissons :

$$ (II-1) \quad b(x, t) = \lim_{\Delta t \to 0} \frac{E\{X(t+\Delta t) - X(t) | X(t) = x\}}{\Delta t} $$

$$ (II-2) \alpha(x, t) = \lim_{\Delta t \to 0} \left[ \frac{E\{[X(t+\Delta t) - X(t)]^2 | X(t) = x\}}{\Delta t} - \frac{E\{X(t+\Delta t) - X(t) | X(t) = x\}}{\Delta t} \right]^2 $$

et les opérateurs :

$$ (II-3) \quad A_{x,t} f = - \frac{\partial}{\partial t} f(x, t) - \frac{\partial}{\partial x} [b(x, t) f(x, t)] + \frac{1}{2} \frac{\partial^2}{\partial x^2} [\alpha(x, t) f(x, t)] $$

$$ (II-4) \quad B_{x,t} f = - b(x, t) f(x, t) + \frac{1}{2} \frac{\partial}{\partial x} [\alpha(x, t) f(x, t)] $$

où $f = f(x, t)$ est la fonction de distribution au temps $t$. Les équations du processus avec retour instantané sont $[6]$, $[19]$, $[20]$ :

$$ (II-5) \quad A_{x,t} f' + \lambda_1 m_1(t) f_1(x) + \lambda_2 m_2(t) f_2(x) = 0 $$

$$ (II-6) \quad \frac{d}{dt} m_1(t) = B_{0,t} f - \lambda_1 m_1(t) $$

$$ (II-7) \quad \frac{d}{dt} m_2(t) = - B_{M,t} f - \lambda_2 m_2(t) $$
où $f_1(x)$ et $f_2(x)$ sont les densités de probabilité immédiatement après un saut depuis O ou M respectivement.

La conservation des masses s'écrit :

$$ (II-8) \quad \frac{\partial}{\partial t} \left[ \int_{0}^{M} f(x, t) dx + m_{1}(t) + m_{2}(t) \right] = 0 $$

$$ \text{Comme } \int_{0}^{M} f_{1}(x)dx = 1 \text{, } \int_{0}^{M} f_{2}(x)dx = 1 $$

une conséquence de la formule (II-8) et des formules précédentes est :

$$ (II-9) \quad \int_{0}^{M} f(x, t) dx + m_{1}(t) + m_{2}(t) = 1 $$

Nous avons de plus les conditions aux limites

$$ (II-10) \quad \lim_{x \to 0} f(x, t) = \lim_{x \to M} f(x, t) = 0 $$

Nous allons maintenant appliquer ces équations générales à notre modèle particulier et en régime permanent (au bout d'un temps infini) c'est-à-dire que nous annulerons toute dérivée partielle par rapport au temps.

Nous considérons ici le cas plus général où les distributions des temps de service de l'unité centrale et du tambour sont arbitraires et caractérisées par leurs moyennes respectives $\frac{1}{\mu}$ et $\frac{1}{\lambda}$ et leurs variances $V_s$ et $V_a$.

Nous supposons aussi que ces mêmes distributions deviennent exponentielles, tout en conservant les mêmes moyennes, lorsque le nombre de clients $x$ dans la file de l'unité centrale devient 0 ou M.

Dans les équations précédentes, nous sommes conduits à faire $\lambda_1 = \lambda$ , $\lambda_2 = \mu$ et $[40]$, $[26]$ :

$$ (II-11) \quad b = \lambda - \mu \text{ et } \alpha = \lambda^3 V_a + \mu^3 V_s $$
Soient $m_1$ et $m_2$ les masses de probabilité en 0 et $\mu$ respectivement en régime stationnaire. Les équations (II-5), (II-6), (II-7) deviennent (compte tenu du régime permanent)

$$ (II-12) \quad \frac{1}{2} \alpha f'' - b f' = - \lambda m_1 \delta(x-1) - \mu m_2 \delta(x-M+1) $$

$$ (II-13) \quad \lim_{x \to 0} \left[ \frac{1}{2} d f' - b f \right] = \lambda m_1 $$

$$ (II-14) \quad \lim_{x \to M} \left[ \frac{1}{2} \alpha f' - b f \right] = \mu m_2 $$

$$ \text{avec } \lim_{x \to 0} f' = \lim_{x \to M} f = 0. $$

L'équation (II-12) doit être intégrée séparément dans les trois régions $0 \le x \le 1$, $1 \le x \le M-1$ et $M-1 \le x \le M$ (en assurant la continuité en 1 et $M-1$). La résolution donne :

$$ (II-15) \qquad f = \begin{cases} \frac{-\lambda m_1}{b} [1 - e^{\gamma x}], & 0 \le x \le 1 \\ \frac{-\lambda m_1}{b} [e^{-\gamma} - 1] e^{\gamma x}, & 1 \le x \le M-1 \\ \frac{-\mu m_2}{b} [e^{\gamma(x-M)} - 1], & M-1 \le x \le M \end{cases} $$
$$ \text{où (II-16)} \qquad \gamma = \frac{2b}{\alpha} $$

La continuité en $M-1$ permet d'obtenir $m_2$. Il vient

$$ (II-17) \qquad m_2 = \frac{\lambda m_1}{\mu} e^{\gamma(M-1)} $$

devient :

$$ (II-18) \quad \int_0^M f \, dx + m_1 + m_2 = 1 $$
d'où l'on tire $m_1$ grâce à (II-15) et (II-17). Posons :

$$ (II-19) \qquad \rho = \frac{\lambda}{\mu} \quad \text{d'où après calculs :} $$

$$ (II-20) \qquad m_1 = \frac{1 - \rho}{1 - \rho^2 e^{\gamma(M-1)}} $$

(II-20) et (II-17) permettent d'obtenir entre $m_1$ et $m_2$ la relation :

$$ (II-21) \qquad (1 - m_2) \lambda = (1 - m_1) \mu $$

## 2.2 - Résolution pratique du modèle

Les données de départ du problème sont $\lambda = 1$, $1/\mu$ donné en para-
mètre ainsi que le nombre M de clients. Comme la loi du tambour a été prise
constante et celle de l'unité centrale exponentielle, $V_s = \frac{1}{2\mu}$ et $V_a = 0$.
Nous calculerons b et a par les formules (II-11), puis $\gamma$ et $\rho$ par (II-16)
et (II-19). D'où $m_1$ d'après (II-20) et $m_2$ d'après (II-21). Le premier résul-
tat obtenu est le taux d'utilisation $\tau$ de l'unité centrale. Ce taux est
défini comme la probabilité que la file de l'unité centrale ne soit pas vide.
Or $m_1$ et $m_2$ sont par définition les probabilités de 0 ou M clients respec-
tivement dans la file de l'unité centrale. D'où :

$$ (II-22) \qquad \tau = 1 - m_1 $$

D'autres résultats à obtenir sont les probabilités $\rho(n)$ de n
clients dans la file d'attente de l'unité centrale. D'après la définition
de $m_2$, nous avons $\rho(M) = m_2$. Pour $1 \le n \le M-1$, les probabilités $\rho(n)$
doivent se déduire de la fonction de distribution f définie par (II-15).
Les probabilités $\rho(n)$ sont définies physiquement pour n entier alors que
la fonction f définit la distribution d'une variable continue en raison de
l'approximation de n par un processus de diffusion. Les résultats obtenus
par simulation guideront sur la façon de retourner au cas discret.
Le tableau de la figure 2 donne les taux d'utilisation de l'unité centrale calculés pour $\frac{1}{\mu} = 0.25, 0.50, 0.75, 0.80, 0.90$ et $M = 1, 3, 5, 7, 9$ à côté des mêmes résultats obtenus par simulation.

### III - RESOLUTION DU MODELE PAR SIMULATION

Le simulateur utilisé pour l'étude du système précédent écrit en SIMULA 67, obéit aux principes généraux exposés au paragraphe I-2. Nous allons examiner successivement le problème des mesures, puis la comparaison des résultats avec le modèle théorique.

#### 1 - Problème des mesures

Les résultats édités sont :

- Les taux d'utilisation des ressources : unité centrale et tambour,

- Les probabilités expérimentales (histogramme) des nombres de clients dans les files d'attente.

Le simulateur devant contrôler et être contrôlé par la résolution théorique exposée précédemment, la comparaison portera sur les taux d'utilisation de l'unité centrale obtenus pour les différentes valeurs des paramètres du système. Cette comparaison pour être valable demandera le calcul d'une précision statistique (cf Chapitre I, §II-2.2) sur ces taux d'utilisation. Nous serons en outre conduits à allonger les temps de simulation pour que cette précision soit au moins de l'ordre de $1 \text{ } \%$ (l'expérience nous a d'ailleurs montré, comme nous le reverrons au chapitre III, que les taux d'utilisation des serveurs sont des variables aléatoires peu dispersées. Une telle précision peut donc être atteinte à un prix raisonnable dans la pratique). La détermination de la précision sur les taux d'utilisation nous conduira à transposer à une telle variable la méthode de blocage exposée au chapitre I, §II-2.4.
## 1.1 - Détermination des grandeurs mesurées

Le taux d'utilisation d'un serveur entre deux instants $t_1$ et $t_2$ est estimé comme le rapport du temps total d'activité du serveur au temps total de mesure soit $t_2 - t_1$. La probabilité de n clients dans une file d'attente est estimée de la même façon comme le rapport de la somme des intervalles de temps où n clients se trouvaient dans la file d'attente au temps total, La simulation est effectuée pendant un temps total TMAX. Les probabilités de n clients dans les files d'attentes sont mesurées sur la durée totale de la simulation. Les taux d'utilisation des serveurs sont prélèvés périodiquement entre deux tops d'horloge pour les raisons étudiées ci-après.

## 1.2 - Détermination de la précision sur les taux d'utilisation des serveurs

Un taux d'utilisation est une variable que l'on sait mesurer entre deux instants $t_1$ et $t_2$ (tops d'horloge par exemple). Transposons au cas de telles variables mesurées comme rapports de temps les méthodes exposées au chapitre I § II-2.2.4 pour les séries de mesures corrélées. La table K de bloc sera remplacée par la notion de période T de prise de mesures qu'il est naturel de prendre comme période d'horloge entre deux prélèvements de mesures. De la même façon qu'au chapitre I, le temps T sera pris suffisamment grand pour que les mesures prises sur deux intervalles de temps différents soient très peu corrélées entre elles. La moyenne partielle $\bar{a}_i$ définie par (I-20) sera remplacée par la mesure $\tau_i$ faite sur le ième intervalle de temps. Soit p le nombre total d'intervalle de temps, pris tous égaux entre eux, TMAX le temps total de simulation, $a_i$ le temps total d'activité du serveur étudié sur le ième intervalle de temps, A son temps total d'activité, $\tau$ le taux d'utilisation mesuré sur le temps total. Nous avons :

$$ (II-22) \qquad TMAX = pT $$

$$ (II-23) \qquad A = \sum_{i=1}^{p} a_i $$

$$ (II-24) \qquad \tau = \frac{A}{TMAX} \quad \text{et} \quad \tau_i = \frac{a_i}{T} $$
par définition d'un taux d'utilisation.

Soit $\bar{\tau}$ la moyenne des $\tau_i$. soit :

$$ (II-25) \qquad \bar{\tau} = \frac{1}{p} \sum_{i=1}^{p} \tau_i $$

Nous avons :

$$ \bar{\tau} = \frac{1}{p} \sum_{i=1}^{n} \frac{a_i}{T} = \frac{1}{pT} \sum_{i=1}^{n} a_i = \frac{A}{TMAX} = \tau $$

Les $\tau_i$ seront appelés taux d'utilisation partiels sur le ième intervalle. $\tau$ sera appelé taux d'utilisation global.

D'après la définition de T (analogue à celle de la taille de bloc K du chapitre I § 2.2.4), les $\tau_i$ sont très peu corrélés entre eux. Nous calculerons leur moyenne $\bar{\tau}$, leur variance $V(\bar{\tau})$ et nous obtiendrons un intervalle de confiance sur $E(\bar{\tau}) = E(\bar{\tau}_0)$ d'une façon analogue au cas de la formule (I-20). Soit $\sigma_{\bar{\tau}} = \sqrt{V_{\bar{\tau}}}$. Nous aurons avec une confiance statistique de 95 % :

$$ (II-26) \qquad \tau - \frac{1,96\sigma_{\bar{\tau}}}{\sqrt{p}} < E(\tau) < \tau + \frac{1,96\sigma_{\bar{\tau}}}{\sqrt{p}} $$

De la même façon que précédemment, la précision obtenue sera inversement proportionnelle à la racine carrée du temps simulé.

## 1.3 - Mise en oeuvre pratique du calcul de précision

Les taux d'utilisation partiels sont relevés sur tops d'horloge loge de période T à déterminer. Il est possible de déterminer une période minimum $T_0$ par des méthodes similaires à celles exposées précédemment. Dans la pratique, il suffit de prendre $T \ge T_0$ minimum pour les mêmes raisons qu'au chapitre 1. Cette remarque permet de déterminer empiriquement une valeur convenable de T par quelques simulations préalables. Pour T croissant et
pour une même durée de simulation, les amplitudes IC(T) des intervalles
de confiance calculés croissent (à mesure que les taux partiels successifs
deviennent de moins en moins corrélés) jusqu'à se stabiliser à une valeur
IC(∞) qui est celle à adopter. Nous pouvons choisir la valeur minimum T₀ de
T de sorte que IC(T₀) ≈ IC(∞) et prendre alors T ≥ T₀.

Dans la pratique, nous avons choisi une période T grande et
égale à 600 fois la somme des temps de service moyens du tambour et de l'unité
centrale. Le temps total simulé TMAX a été choisi égal à 60 000 unités de temps
pour obtenir uniformément la précision voulue (inférieure à 1.0/°) pour toutes
valeur des paramètres. Nous obtenons ainsi de 50 à 100 mesures de taux partiels.

## 2 - Résultats obtenus par simulation et comparaison avec la théorie

Les résultats sont toujours obtenus avec les hypothèses : tambour
à temps de service constant égal à 1 et unité centrale à temps de service
exponentiel de moyenne $1/μ$. Les résultats sont donnés pour $1/μ = 0.25, 0.50,$
$0.75, 0.80, 0.90$ et un nombre M de clients dans le système prnent les valeurs :
1, 3, 5, 7, 9.

### 2.1 - Validation réciproque des deux résolutions

La validation a été effectuée par comparaison des taux d'utilisation de l'unité centrale obtenus mathématiquement et par simulation. Cette comparaison est résumée par le tableau de la figure 2 où les résultats de simulation sont donnés avec indication de la précision. L'examen de ce tableau montre que :

- Les résultats des deux méthodes de résolution différent au plus de quelques unités sur la troisième décimale.

- Les précisions calculées en simulation sont de l'ordre de grandeur des différences avec les résultats mathématiques. En particulier pour M ≥ 3 et $1/μ ≥ 0.50$, la plupart des résultats mathé-
Figure 2

Comparaison des taux d'utilisations de l'unité centrale obtenus par la méthode des équation de diffusion et par simulation.

D : diffusions

S : simulation (sous la forme : résultat ± précision).

<table><thead><tr><th>1/μ</th><th>0.25</th><th>0.50</th><th>0.75</th><th>0.80</th><th>0.90</th></tr></thead><tbody><tr><th>M</th><td></td><td></td><td></td><td></td><td></td></tr><tr><th>1</th><td>D 0.200<br>S 0.198<sup>+</sup>2x10<sup>-3</sup></td><td>D 0.333<br>S 0.330<sup>+</sup>2x10<sup>-3</sup></td><td>D 0.429<br>S 0.426<sup>+</sup>3x10<sup>-3</sup></td><td>D 0.444<br>S 0.442<sup>+</sup>3x10<sup>-3</sup></td><td>D 0.474<br>S 0.472<sup>+</sup>3x10<sup>-3</sup></td></tr><tr><th>3</th><td>D 0.248<br>S 0.247<sup>+</sup>2x10<sup>-3</sup></td><td>D 0.482<br>S 0.483<sup>+</sup>4x10<sup>-3</sup></td><td>D 0.685<br>S 0.684<sup>+</sup>5x10<sup>-3</sup></td><td>D 0.719<br>S 0.715<sup>+</sup>5x10<sup>-3</sup></td><td>D 0.781<br>S 0.775<sup>+</sup>5x10<sup>-3</sup></td></tr><tr><th>5</th><td>D 0.250<br>S 0.247<sup>+</sup>2x10<sup>-3</sup></td><td>D 0.498<br>S 0.495<sup>+</sup>4x10<sup>-3</sup></td><td>D 0.729<br>S 0.730<sup>+</sup>6x10<sup>-3</sup></td><td>D 0.770<br>S 0.770<sup>+</sup>6x10<sup>-3</sup></td><td>D 0.843<br>S 0.839<sup>+</sup>6x10<sup>-3</sup></td></tr><tr><th>7</th><td>D 0.250<br>S 0.247<sup>+</sup>2x10<sup>-3</sup></td><td>D 0.500<br>S 0.495<sup>+</sup>4x10<sup>-3</sup></td><td>D 0.743<br>S 0.740<sup>+</sup>6x10<sup>-3</sup></td><td>D 0.788<br>S 0.786<sup>+</sup>6x10<sup>-3</sup></td><td>D 0.868<br>S 0.864<sup>+</sup>6x10<sup>-3</sup></td></tr><tr><th>9</th><td>D 0.250<br>S 0.247<sup>+</sup>2x10<sup>-3</sup></td><td>D 0.500<br>S 0.495<sup>+</sup>5x10<sup>-3</sup></td><td>D 0.747<br>S 0.743<sup>+</sup>7x10<sup>-3</sup></td><td>D 0.795<br>S 0.791<sup>+</sup>7x10<sup>-3</sup></td><td>D 0.885<br>S 0.877<sup>+</sup>7x10<sup>-3</sup></td></tr></tbody></table>
matiques sont exactement dans l'intervalle déterminé par la simulation.

Nous pouvons conclure que :

- Dans le cas des hypothèses faites sur le modèle, l'approximation du phénomène par un processus de diffusion est très bonne (résultats approchés à mieux que $1 \, \%$).

- Les méthodes utilisées dans la simulation sont correctes.

Les calculs de précision donnent bien l'ordre de grandeur de celle-ci.

## 2.2 - Etude de la file d'attente de l'unité centrale

Nous allons étudier maintenant les probabilités $P(n)$ de n clients dans la file d'attente de l'unité centrale. Ces probabilités ont été étudiées en simulation comme dit plus haut. La méthode théorique a permis d'oolenir jusqu'à présent :

$$P(0) = m_1 \text{ et } P(M) = m_2$$

Ces résultats s'accordent convenablement avec la simulation. Comme $P(0) = 1 - \tau$, $\tau$ étant le taux d'utilisation de l'unité centrale, la comparaison effectuée figure 2 s'applique. Pour $P(M)$, les ordres de grandeur des résultats correspondent (du fait que $P(M)$ est toujours faible pour $M$ supérieur à 3, les déterminations expérimentales deviennent imprécises). L'accord est cependant visible sur les queues de courbes des figures 4 à 8.

Pour $1 \le n \le M-1$, la théorie (équations (II-15) ne fournit de résultats que pour n continu, alors que la simulation a fourni ces mêmes résultats pour n discret. Nous sommes conduits à la recherche d'un algorithme permettant de retourner du cas continu au cas discret.
<figure><img src="image_4.png" /></figure>
(II-15) se transforme alors en :

$$ (II-33) \quad f = \begin{cases} C_1 [1 - e^{\gamma^x}] , & 0 \le x \le 1 \\ C_2 e^{\gamma^x} , & 1 \le x \le M-1 \\ C_3 [e^{\gamma(x-M)} - 1], & M-1 \le x \le M \end{cases} $$

dont l'intégration sur x fournit :

$$ (II-34) \quad F = \begin{cases} m_1 + C_1 \left[x + \frac{1}{\gamma} (1 - e^{\gamma x})\right] , & 0 \le x \le 1 \\ m_1 + C_1 \left[1 + \frac{1}{\gamma} (1 - e^{\gamma})\right] + \frac{C_2}{\gamma} [e^{\gamma x} - e^{\gamma}] , & 1 \le x \le M-1 \\ m_1 + C_1 + \frac{C_2}{\gamma} e^{\gamma(M-1)} + C_3 \left[M - x - 1 + \frac{1}{\gamma} [e^{\gamma(x-M)} - e^{-\gamma}]\right] , & M-1 \le x < M \end{cases} $$

(II-27) s'écrira donc :

$$ (II-35) \quad \rho(n) = F(\beta_n) - F(\beta_{n-1}) $$

Le problème se ramène à la détermination des bornes $\beta_n$ de façon
à retrouver convenablement les résultats de simulation.

### 1er algorithme : Découpage de l'intervalle $]0, M[$ en intervalles constants :

Nous ferons quel que soit $n$ :

$$ \beta_n - \beta_{n-1} = \frac{M}{M-1} $$

ce qui avec (II-28) $\beta_0 = 0$, (II-35) et (II-34) permet de calculer $\rho(n)$. Les résultats obtenus sont médiocres (sauf en 0 et M, cas traités à part comme vu plus haut). Les probabilités pour les petites valeurs de n sont très sous-estimées au profit des autres qui sont exagérées. La distribution obtenue tend à être uniformisée. Deux courbes comparatives types sont représentées
<caption>Figure 3</caption>

Exemple d'application de l'algorithme des intervalles égaux pour le calcul de la file d'attente de l'unité centrale

$$1/\mu = .5$$

M = 5

Calcule

Simulation

<figure><img src="image_5.png" /></figure>
schématiquement sur la figure 3.

## 2ème algorithme :

P(1) étant sous-estimé, recherchons un algorithme qui favorise le premier intervalle. Nous remarquons que la concentration en 1 de la densité de probabilité comprise entre 0 et 1 (plus un petit intervalle au delà de 1) a pour effet d'augmenter l'espérance mathématique puisque :

Pour $0 \le x \le 1$, on a $x \le f \le f$ d'où :

$$ \int_{0}^{X} x f(x) \, dx < \int_{0}^{X} f(x) \, dx $$

pour $0 \le X \le 1$

Nous chercherons à remplacer notre distribution continue par une distribution discrète de même espérance mathématique. Nous tenterons de respecter l'espérance mathématique terme à terme, le n ième terme dans le cas continu correspondant à l'intégrale de $x f(x)$ sur le n ième intervalle. Nous poserons donc :

$$ (II-36) \qquad n \mathcal{P}(n) = \int_{\beta_{n-1}}^{\beta_n} x f(x) \, dx $$

Définissons $F_e(x) = \int_0^x u f(u) \, du$ (l'espérance mathématique du nombre $x$ de clients dans la file d'attente de l'untié centrale est donnée par $f_e(M) + m_2 M$ en raison de la masse ponctuelle de probabilité $m_2$ au point $M)$.

La fonction $F_e$ est donnée par (d'après (II-33)):
$$
(II-37) \quad F_e = \begin{cases}
C_1 \frac{x^2}{2} - \frac{C_1}{\gamma} \left[e^{\gamma x} \left(x - \frac{1}{\gamma}\right) + \frac{1}{\gamma}\right] , & 0 \le x \le 1 \\
\frac{C_1}{2} - \frac{C_1}{\gamma} \left[e^{\gamma} \left(1 - \frac{1}{\gamma}\right) + \frac{1}{\gamma}\right] + \frac{C_2}{\gamma} \left[e^{\gamma x} \left(x - \frac{1}{\gamma}\right) - e^{\gamma} \left(1 - \frac{1}{\gamma}\right)\right] , & 1 \le x \le M-1 \\
\frac{C_1}{2} - \frac{C_1}{\gamma} \left[e^{\gamma} \left(1 - \frac{1}{\gamma}\right) + \frac{1}{\gamma}\right] + \frac{C_2}{\gamma} \left[e^{\gamma(M-1)} \left(M-1-\frac{1}{\gamma}\right) - e^{\gamma} \left(1-\frac{1}{\gamma}\right)\right] + \\
\qquad + C_3 \left[ \frac{(M-1)^2 - x^2}{2} + \frac{1}{\gamma} \left[ \left(x - \frac{1}{\gamma}\right) e^{\gamma(x-1)} - (M-1-\frac{1}{\gamma}) e^{-\gamma} \right] \right] , & M-1 \le x \le M
\end{cases}
$$

L'équation (II-36) devient :

$$
(II-38) \quad n \mathcal{F}(n) = F_e(\beta_n) = F_e(\beta_{n-1})
$$

d'où d'après (II-35) :

$$
(II-39) \quad n [F(\beta_n) - F(\beta_{n-1})] = F_e(\beta_n) - F_e(\beta_{n-1})
$$

équation fournissant $\beta_n$ à partir de $\beta_{n-1}$, donc par récurrence, $\beta_n$ quel que
soit n puisque $\beta_0 = 0$. $\beta_{n-1}$ étant connu, l'équation (II-39) en $\beta_n$ est résolue
numériquement par approximations successives. Une solution parasite de (II-39)
est $\beta_n = \beta_{n-1}$. Cette solution devant être écarté, l'algorithme ne sera pra-
tiquable que s'il existe une autre solution $\beta_n$ supérieure à $\beta_{n-1}$. Cela est
bien le cas pour n=1 (puisque $\beta_0 = 0$). Malheureusement, la résolution numé-
rique pour certaines données a fait apparaître des valeurs de $\beta_1$ de l'ordre
de 2.5. Or si $\beta_{n-1} > n$, nous avons pour $x \ge \beta_{n-1}$ :

$$
n f(x) \leq x f(x)
$$

qui donne par intégration :

$$
(II-40) \quad n [F(x) - F(\beta_{n-1})] < F_e(x) - F_e(\beta_{n-1})
$$
ce qui montre que l'équation (II-39) dans ce cas n'a pas de solution supérieure à $\beta_{n-1}$. Dans le cas précité, on ne peut déterminer $\beta_2$. L'application de l'algorithme est donc impossible dans ce cas : la seule possibilité de remplacement de la distribution continue par une distribution discrète, tout en respectant l'espérance mathématique terme à terme, conduirait à donner la valeur 0 à certaines probabilités $\mathcal{P}(n)$.

**Algorithme retenu :**

Nous avons finalement retenu un algorithme intermédiaire entre les deux précédents. Un calcul à la main des bornes $\beta_n$ successives à partir des résultats de simulation nous a en effet montré que $\beta_1 - \beta_0 = \beta_1$ devait être pris très supérieur à l'amplitude $\beta_n - \beta_{n-1}$ (pour $n \ge 2$) des autres intervalles, ceux-ci étant à peu près constants. Cette constatation nous suggère donc l'algorithme suivant :

- Appliquer le second algorithme pour déterminer $\beta_1$, c'est-à-dire résoudre l'équation (tirée de (II-39) en faisant $\beta_0 = 0$,

$$F(0) = m_1 \text{ et } F_e(0) = 0$$

$$(II-41) \qquad F(\beta_1) - m_1 = F_e(\beta_1)$$

- Déterminer les bornes $\beta_n$ pour $n \ge 2$ en découpant l'intervalle $[\beta_1, M]$ en $M-2$ intervalles égaux soit

$$(II-42) \qquad \beta_n - \beta_{n-1} = \frac{M - \beta_1}{M-2} \quad \text{pour } n \ge 2.$$

Tous les $\beta_n$ étant déterminés, la formule (II-35) donne les valeurs successive de $\mathcal{P}(n)$.

Les valeurs calculées des probabilités $\mathcal{P}(n)$ ont été représentées sous forme de courbes pour les différentes valeurs des paramètres à côté des
mêmes courbes déterminées à partir des simulations. Les figures 4 à 8 donnent quelques exemples de cette série de courbes. Les résultats de simulation sont retrouvés par le calcul avec une bonne précision.

IV - CONCLUSION

L'étude précédente a donné un exemple de validation d'une approximation mathématique par une simulation. une première partie de l'étude a permis la vérification réciproque des deux méthodes de résolution.

Les techniques de comparaison exposées dans cette première partie sont applicables à d'autres types de modèles. Dans une seconde partie, la facilité d'obtenir par simulation des résultats de nature très diverse a permis l'utilisation d'équations mathématiques pour l'obtention par le calcul de résultats plus généraux.

L'étude du modèle entre les points sur lesquels la comparaison a pu porter sera désormais effectuée entièrement à partir des équations. Celles-ci, une fois vérifiées et interprétées à l'aide de l'expérience, permettent en effet des résolutions numériques moins coûteuses ainsi que l'étude mathématique des variations du comportement du modèle, lorsque les paramètres de celui-ci varient.
<figure><img src="image_6.png" /><figcaption>Figure 4</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6</figcaption></figure>

Nombre de clients 3

Temps de service moyen UC .80
<figure><img src="image_9.png" /><figcaption>Figure 7</figcaption></figure>

Nombre de clients 3

Temps de service moyen UC .80
<figure><img src="image_10.png" /><figcaption>Figure 8</figcaption></figure>
# CHAPITRE III

## UTILISATION D'UN MODELE DE SIMULATION POUR L'ETUDE DE LA REGULATION ADAPTATIVE D'UN SYSTEME INFORMATIQUE

### I - INTRODUCTION

La partie précédente de cette étude a traité des techniques de mesures et des méthodes de calculs de précision statistique utilisées par nous en simulation. Nous en avons montré ensuite une application à la validation réciproque de deux méthodes de résolution d'un modèle. Nous allons maintenant donner un exemple d'application de la simulation à l'étude d'une régulation adaptative d'un système informatique lorsque les caractéristiques de la charge du système varient dans le temps. Cette étude montrera l'utilisation possible des mesures de performances à l'interaction sur le système lui-même de façon à constituer une boucle d'asservissement. Nous aurons également un exemple de mise en oeuvre de techniques statistiques adaptées à la mesure d'un système non stationnaire. D'autre part, le comportement de la régulation dans le temps fournira un exemple d'étude par simulation d'un phénomène trés itoire.

Nous présèntrons d'abord le modèle de système étudié et la validation du simulateur. Nous étudierons ensuite le comportement du modèle suivant les conditions de fonctionnement et nous en déduirons l'intérêt d'une régulation adaptative. Nous proposerons un algorithme applicable à une telle régulation à partir de mesures de performances sur le système. Nous définirons alors une bonne mesure physique des performances et nous étudierons le système régulé du double point de vue de la bonne adaptation de la régulation suivant les algorithmes définis et de l'amélioration des performances du système.
# II - PRESENTATION DU MODELE DE SYSTEME ETUDIE

## 1 - Description du modèle

Le modèle de simulation étudié [ 1 ] est représenté sur la figure 1 et modélise un système de temps partagé à mémoire virtuelle. Il comprend le système proprement dit et un ensemble de N consoles à partir desquelles sont générés des programmes. Le système comprend trois serveurs : unité centrale UC, mémoire secondaire (ou tambour de pagination) MB et disque de fichiers DF modélisant les entrées-sorties.

### 1.1 - Mécanisme des éléments du système

a) **Consoles** : Les consoles modélisent une suite d'interactions entre un utilisateur main et le système. Une console émet une interaction (ou un programme) au bout d'un temps de réflexion aléatoire de loi de distribution exponentielle négative donnée par sa moyenne $\tau$.

La console attend la sortie du programme qu'elle a émis et émet un nouveau programme au bout d'un nouveau temps de réflexion.

Le temps moyen de réflexion $\tau$ est un paramètre du modèle. Il est possible en particulier de faire $\tau = 0$ pour maintenir dans le système un nombre de programmes constant égal à N.

b) **Unité centrale** : Les programmes entrant dans le système sont mis en file d'attente derrière l'unité centrale. L'exécution d'un programme sur l'UC se poursuit jusqu'à l'un des trois évènements suivants :

- **Fin d'exécution** : Le temps total de calcul d'un programme est généré par la console et décrémenté à chaque service à l'unité centrale. A l'épuisement de ce temps, le programme quitte le système et retourne à la console qui l'avait généré.

- **Défaut de page, auquel cas le programme est mis en file d'attente derrière la mémoire secondaire**.

- **Entrée-sortie** : dans ce cas, le programme est mis en file d'attente derrière le disque de fichiers.
<figure><img src="image_11.png" /><figcaption>Figure 1</figcaption></figure>
Dans chacun de ces trois cas, nous tiendrons compte d'un temps d'activité monitée (é constant et fixé pour chaque cas) afin de tenir compte du temps d'unité centrale dépensé par le superviseur.

c) mémoire secondaire et disque de fichiers : Les demandes de services à la mémoire secondaire et au disque de fichiers sont traitées dans l'ordre "Premier-arrivé-premier-servi". Le fonctionnement des unités est représenté par le temps total demandé pour un service. Ce temps est égalatoire et sa distribution (donnée par son type et ses paramètres) pour chaque unité est un par mètre du simulateur.

## 1.2 - Comportement des programmes à l'unité centrale

Le comportement d'un programme, fixé par un ensemble de paramètres générés lors du départ de sa console, commande le déroulement d'un service à l'unité centrale. Dans notre modélisation, tous les comportements possibles peuvent être caractérisés par les distributions de trois variables :

- Temps total de calcul du programme,

- Temps ininterrompu à l'unité centrale entre deux demandes d'entrées-sorties,

- Temps ininterrompu à l'unité centrale entre deux défauts de page.

Les deux premières distributions ne dépendent que du programme considéré. La troisième dépend également de l'état du système. Cette dernière loi supposée exponentielle a sa moyenne q donnée par la formule de Belady [ 5 ]

$$ (III - 1) \qquad q = a \ s^k $$

où s est la taille de la mémoire allouée au programme k un coefficient dépendant du programme et que nous appelerons sa localité et a un coefficient dépendant de la machine et du système. Soit M la mémoire totale réservée à l'exécution des programmes et n le degré de multiprogrammation (nombre de programmes dans le système par définition), à un instant donné, Nous supposerons la mémoire également répartie entre tous les programmes, c'est-à-dire
que nous ferons $s = \frac{M}{n}$ dans la formule (III-1). D'où :

$$ (III-2) \qquad q = a \left(\frac{M}{n}\right)^k $$

Nous avons pris $a = 0.01$. Dans la pratique, $1.5 \le k \le 2.5$ [ 4 ].

Les lois du temps total de calcul et du temps entre entrées-sorties ont été prises exponentielles de moyennes c et r. Les paramètres d'un programme sont donc :

- Son temps total de calcul, tiré au hasard lors de sa création

- Son temps moyen r entre entrées-sorties et sa localité k.

Nous supposerons (au moins pendant des intervalles de temps longs) les programmes générés toujours statistiquement semblables : même localité k, même temps moyen r entre entrées-sorties, et temps total de calcul généré toujours avec la même loi et la même moyenne c. Dans une première partie de l'étude, ces caractéristiques resteront constantes tout au long de la simulation. Ces hypothèses, avec celle du partage équitable de la mémoire, signifient que l'on s'intéresse à un instant donné qu'aux caractéristiques moyennes de l'ensemble des programmes dans le système.

## 1.3 - Hypothèses sur les lois du modèle

Nous avons déjà mentionné que les lois liées au comportement des programmes étaient supposées exponentielles. Les lois de service au disque de fichiers et à la mémoire secondaire sont prises elles aussi exponentielles de moyennes respectives $\frac{1}{\mu_2}$ et $\frac{1}{\mu_1}$. Ces hypothèses sont pour but de faciliter la validation du modèle par des méthodes de résolution mathématique (voir ci-après). Des essais ont cependant été effectués en supposant constants les services au disque de fichiers et à la mémoire secondaire. Les figures 3a et 3b montrent les différences de résultats obtenus entre ce cas et le cas des lois exponentielles. Les différences obtenues sont relativement peu importantes. Ce fait est dû à ce que les résultats dépendent
essentiellement des moyennes des lois du modèle, beaucoup plus que de la
forme de ces lois. D'où la légitimité du choix de lois exponentielles,
choix destiné à faciliter des calculs théoriques servant de base de
comparaison.

## 2 - Mesures prises au cours de la simulation

Diverses grandeurs sont prélevées à divers instants sur le
système au cours de la simulation afin d'élaborer des statistiques pour la
mesure des performances. Ce sont :

- Les temps de réponse des programmes, c'est-à-dire la différence entre les heures de sortie et d'entrée dans le système.

- Le degré n de multiprogrammation, c'est-à-dire le nombre de programmes dans le système.

- Les taux d'utilisations des unités : unité centrale, mémoire secondaire et disque de fichiers, prélevés sur l'ensemble de la simulation ou entre deux instants $t_1$ et $t_2$.

- La dilatation et l'efficacité des clients que nous définirons à propos de la régulation.

Les statistiques effectuées sur ces grandeurs ont un double but,
suivant leur nature. Les unes valuent globalement les performances de
l'ensemble du système et sont prises sur l'ensemble de la simulation. Elles
comprènnent les moyennes (et parfois les variances) de toutes les grandeurs
décrites précédemment. C'est ainsi qu'ils ont été construits les tableaux des
figures 17 et 18. D'autres statistiques évaluent les variations
de performances obtenues suivant l'état de certaines variables. Nous en
avons fait l'application à l'étude des variations des taux d'utilisation (et
de la dilatation qui leur est liée, voir ci-après) en fonction du degré n de
multiprogrammation (courbes 3a à 4 pour les taux d'utilisation de
l'unité centrale, courbes 11 et 12 pour la dilatation). Deux techniques
sont utilisables pour obtenir ces résultats (dans l'exemple de notre appli-
cation, mais l'extension à d'autres cas est facile) :
- Effectuer plusieurs simulations à n constant pour les différentes valeurs de n (par exemple en fixant le nombre de consoles égal à n et en annulant les temps de réflexion). Nous sommes ramenés à l'étude par simulation d'un système fermé (cf. chapitre II § I-1).

- Obtenir tous les résultats en une seule simulation (à n n figé) par prélèvement conditionnel suivant la valeur de n. Pour une valeur de n donnée, les statistiques sont mises à jour toutes les fois que la valeur de n donnée se réalise dans le système

Sans être entièrement équivalentes du point de vue théorique, ces deux méthodes conduisent en fait à des résultats comparables dans le cas de notre application.

L'intérêt principal de la seconde technique est la possibilité d'évaluer en temps réel (avec des techniques d'oubli du passé appropriées) les courbes de variation des performances d'où son utilisation possible à une interaction sur le système (voir plus loin). Nous l'avons appliquée aussi à la validation du simulateur par comparaison des chiffres obtenus avec des résultats mathématiques, ce qui valide à la fois le simulateur et la technique de mesures. Nous avons cependant utilisé la première technique (simulation d'un réseau fermé) pour certaines études systématiques des variations (courbes 5 et 6).

## 3 - Validation du simulateur

La validation du modèle de simulation exposé ici a été entièrement décrite dans [ 1 ] que nous reproduisons partiellement pour l'exposé de la méthode et les calculs.

La validation s'opère par comparaison avec des résultats obtenus mathématiquement. La comparaison a porté sur les taux d'utilisation de l'unité centrale en fonction du degré de multiprogrammation. Le modèle de la figure 1 est remplacé pour cela par le modèle fermé de la figure 2, où n est le nombre de clients dans le système. Outre les hypothèses déjà faites
<figure><img src="image_12.png" /><figcaption>Figure 2</figcaption></figure>
sur l'exponentialité des lois, nous supposerons nuls les temps d'activité moniteur lorsqu'un programme quitte l'unité centrale. Nous allons montrer qu'avec ces hypothèses, le modèle est un réseau de Jackson [ 32 ] à trois serveurs et nous en déduirons les formules de résolution.

Montrons d'abord que l'unité centrale est équivalente à un serveur exponentiel : son temps de service $t$ est le minimum de trois temps exponentiels $t_1, t_2$ et $t_3$ de fonctions de répartition respectives :

$$ (III-3) \quad F_1(t_1) = 1 - e^{-\frac{t_1}{c}} \quad \text{pour le temps } t_1 \text{ avant fin de calcul} $$

$$ (III-4) \quad F_2(t_2) = 1 - e^{-\frac{t_2}{q}} \quad \text{pour le temps } t_2 \text{ avant défaut de page} $$

$$ (III-5) \quad F_3(t_3) = 1 - e^{-\frac{t_3}{r}} \quad \text{pour le temps } t_3 \text{ avant demande d'entrée-sortie} $$

Soit $F(t)$ la fonction de répartition du temps de service ininterrompu $t$. On a :

$$ F(t) = \text{Prob} (\tau \le t) \text{ avec } \tau = \min (t_1, t_2, t_3) \text{ d'où} $$

$$ F(t) = 1 - \text{Prob} (t_1 > t) \cdot \text{Prob} (t_2 > t) \cdot \text{Prob} (t_3 > t) \text{ d'où} $$

$$ F(t) = 1 - [1 - F_1(t)] [1 - F_2(t)] [1 - F_3(t)] \text{ soit} $$

$$ (III-6) \quad F(t) = 1 - e^{-t(\frac{1}{c} + \frac{1}{q} + \frac{1}{r})} = 1 - e^{-\frac{U}{o}} t $$

en posant :

$$ (III-7) \quad \frac{U}{W_o} = \frac{1}{c} + \frac{1}{q} + \frac{1}{r} $$

Soient $p_0, p_1, p_2$ les probabilités respectives de sortie du système, de défaut de page et de demande d'entrée sortie du programme à la fin du service à l'unité centrale. Nous avons :
$$p_1 = \int_0^\infty \text{Prob} \left[ (t - dt \le t_1 \le t) \text{ et } (t_2 > t) \text{ et } (t_3 > t) \right] dt$$

d'où

$$p_1 = \int_0^\infty dt \frac{1}{q} e^{-\frac{t}{c}} e^{-\frac{t}{q}} e^{-\frac{t}{r}} = \frac{1}{q} \int_0^\infty e^{-\frac{t}{U_0}} e^{-\frac{t}{r}} dt \text{ soit}$$

$$(III-8) \qquad p_1 = \frac{1}{U_0 q}$$

De la même façon

$$(III-9) \qquad p_2 = \frac{1}{U_0 r} \quad \text{d'où}$$

$p_0 = 1 - p_1 - p_2$ soit après calculs :

$$(III-10) \qquad p_0 = \frac{1}{U_0 c}$$

Le modèle vérifie bien les hypothèses de Jackson [32] et se réscut en régime stationnaire (au bout d'un temps infini). Soit $p(n_0, n_1, n_2)$ la probabilité de $n_0$ clients dans la file d'attente de l'unité centrale, $n_1$ dans celle de la mémoire secondaire, $n_2$ dans celle du disque de fichiers et $A_0, A_1, A_2$ les taux d'utilisation de ces trois unités. Porons :

$$(III-11) \qquad G(N) = \sum_{\substack{n_1, n_2 \\ n_1 + n_2 \le N}} (p_1 \frac{U_0}{U_1})^{n_1} (p_2 \frac{U_0}{U_2})^{n_2}$$

Les formules de Jackson [32] donnent :

$$(III-12) \qquad p(n_0, n_1, n_2) = \frac{1}{G(n)} (p_1 \frac{U_0}{U_1})^{n_1} (p_2 \frac{U_0}{U_2})^{n_2}$$

Comme $A_0(n) = \sum_{i \le n_0 \le n} p(n_0, n_1, n_2)$, nous avons après sommation :
$$ (III-13) \qquad A_0(n) = \frac{G(n-1)}{G(n)} $$

Nous avons ainsi :

$$ (III-14) \qquad A_1(n) = p_1 \frac{U_0}{U_1} \quad \text{ou} \quad A_1(n) = \frac{1}{U_1 q} A_0(n) \quad \text{et} $$

$$ (III-15) \qquad A_2(n) = p_2 \frac{U_0}{U_2} \quad \text{ou} \quad A_2(n) = \frac{1}{U_2 r} A_0(n) $$

La comparaison des résultats théoriques et de simulation est illustrée par les courbes de la figure 3a avec les paramètres :

- Taille mémoire M de 128 pages

- Localité k des programmes prise à 1.5

- Temps moyen r entre entrées-sorties de 10 ms et 30 ms

- Temps de service moyens $\frac{1}{J_1}$ et $\frac{1}{U_2}$ à la mémoire secondaire pris à 5 ms et 30 ms respectivement

- Moyenne du temps total de calcul à l'unité centrale pris à 800 ms.

Les courbes donnent la comparaison des taux d'utilisation de l'unité centrale obtenus par les deux méthodes lorsque le degré n de multi-programmation varie.

## 4 - Comportement du modèle

Divers essais préalables ont été effectués en faisant varier le nombre de consoles et leur rapidité. Lorsque le système est très chargé, on constate qu'il s'écroule : le temps de réponse des programmes croît très rapidement avec la charge du système ; les taux d'utilisation de l'unité centrale et du disque de fichiers décroissent jusqu'à des valeurs faibles ; le taux d'utilisation de la mémoire secondaire devient très proche de 1 (l'accès à celle-ci est donc saturé). Les courbes des figures 3a à 6
<caption>Figure 3 a</caption>

Utilisation
unite centrale

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

Depuis de
multi-ouïranisation
<figure><img src="image_14.png" /><figcaption>Figure 3 b</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 3. c</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 4</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 5</figcaption></figure>
Figure 6

Taux d'utilisation effective UC (sans les temps moniteur) en fonction du degré de multiprogrammation pour différentes valeurs du temps de service moyen $1/U_1$ de la mémoire secondaire

$k = 2$

$r = 30$ ms

<figure><img src="image_18.png" /></figure>
montrent que ce phénomène se produit lorsque le degré n de multiprogramma-
tion est trop élevé : des usagers de plus en plus nombreux dans le système
s' voient allouer de moins en moins de mémoire et font de plus en plus de
léfauts de page. L'apparition du phénomène d'écroulement lorsque n augmente
est liée essentiellement à la localité moyenne k des programmes et au temps
de service $1/\eta_1$ de la mémoire secondaire (constante de la machine). Une
localité k élevée est favorable à un bon fonctionnement du système
(courbes 3c, 4, 5, 6). Les courbes 5 et 6 montrent l'intérêt d'une
mémoire secondaire la plus rapide possible. Avec une mémoire secondaire trop
lente (temps de service moyen 100 ms sur la figure 5 ), aucune multiprogram-
mation avec pagination n'est plus possible.

La bibliographie sur le phénomène d’écroulement est abondante
[15], [16], [45]. Il apparaît que le degré de multiprogrammation est un
paramètre important dont dépendent les performances du système. Ceci pose
le problème de la recherche d’un degré optimal de multiprogrammation :
combien de programmes peut-on admettre à se partager les ressources ? Ce
problème est l’objet de la régulation que nous allons maintenant exposer.

### III - ÉTUDE D'UN PRINCIPLE DE REGULATION DU SYSTEME

#### 1 - Présentation du problème de la régulation

Pour empêcher l'écroulement, les programmes demandant à entrer
lors de leur émission par les consoles sont mis en file d'attente devant
l'entrée du système jusqu'à leur admission par le régulateur de travaux.

Nous laisserons délibérément de côté les méthodes de régulation
fondées sur la notion d'ensemble de travail [14]. Nous nous attach erons à
montrer la possibilité d'une régulation directe du degré de multiprogramma-
tion par contrôle de l'entrée des programmes (nous ne prenons pas en compte
leur évacuation sur mémoire secondaire). D'où le nouveau schéma de système
représenté par la figure 7 où la commande d'admission des programmes (dans
leur ordre d'arrivée dans nos essais) est symbolisée par l'interrupteur
d'entrée K. Divers algorithmes de gestion de l'interrupteur K sont possibles :
<figure><img src="image_19.png" /><figcaption>Figure 7</figcaption></figure>
- Entrée des programmes dans le système jusqu'à concurrence d'un degré de multiprogrammation maximum $\bar{n}_0$ choisi une fois pour toutes. Or le meilleur $\bar{n}_0$ possible (par exemple la valeur de n donnant le meilleur taux d'utilisation de l'unité centrale, ou le maximum sur tout autre critère, voir plus loin) dépend essentiellement des caractéristiques moyennes des clients dans le système (localité et temps moyen entre entrées-sorties dans notre modèle) (courbes (3c, 4, 5 et 6). Ces caractéristiques moyennes pouvant varier dans le temps, nous serons conduits à choisir prudemment un $\bar{n}_0$ faible correspondant à des caractéristiques pessimistes. Le fonctionnement du système ne pourra être que sous-optimal.

- Estimation en temps réel d'un degré optimal de multiprogrammation. Cette méthode permet de tenir compte des variations dans le temps des caractéristiques moyennes de la charge du système. Elle est pratiquable si ces variations sont lentes par rapport aux temps d'estimation nécessaires. Une telle régulation adaptative fera l'objet de la suite de cette étude.

## 2 – Description d'ensemble de la régulation adaptative proposée

La régulation que nous allons décrire entretient à tout instant une estimation $\bar{n}_0$ du degré optimal de multiprogrammation et admet des programmes dans le système jusqu'à concurrence de $\bar{n}_0$. La régulation s'effectue uniquement par l'admission des usagers à l'exclusion de toute évacuation sur mémoire secondaire.

La définition et l'estimation d'un optimum sur le degré de multiprogrammation demandent le choix préalable d'un critère bien représentatif des performances du système et que nous chercherons à maximiser.

Laissant de côté pour l'instant un tel choix, nous allons développer d'abord les techniques de mesures et d'adaptation utilisées. Nous ferons les hypothèses suivantes sur le critère utilisé :
<figure><img src="image_20.png" /><figcaption>Figure 8</figcaption></figure>
- Le critère est défini par l'espérance mathématique d'une grandeur aléatoire X mesurable périodiquement dans le système.

- Le critère est une fonction du degré de multiprogrammation $n$. Cette fonction présente un maximum et un seul et dépend en outre des caractéristiques moyennes (localité et temps moyen entre entrées-sorties) des programmes dans le système.

Le schéma fonctionnel du système régulé est représenté figure 8.

La régulation proposée comprend trois modules :

- Un estimateur EST appelé à chaque nouvelle mesure tenant à jour les estimations du critère pour chaque valeur $n$ du degré de multiprogrammation, tout en tenant compte des variations dans le temps consécutives aux variations de caractéristiques des programmes.

- Un optimiseur OPT qui estime en permanence la valeur optimale du degré de multiprogrammation.

- Une commande K d'admission des usagers dans le système. K n'est autre que l'interrupteur d'admission des programmes dans le système représenté figure 7. K est activé soit lors d'une arrivée de programme (pour le faire éventuellement entrer) soit lors d'un départ du système, soit lors d'un réajustement de $\bar{n}_o$.

Nous allons maintenant décrire plus particulièrement EST et OPT et développer les techniques de mesures et d'estimations.

## 3 - Description de l'estimateur EST

L'estimateur EST est activé toutes les fois qu'une mesure est prise dans le système et met à jour les estimations de performances. Au bout d'un nombre suffisant de mesures, EST appelle l'optimiseur OPT pour un nouvel ajustement de l'estimation $\bar{n}_o$. Nous noterons $X_n^i$ la ième mesure effectuée
pour le degré de multiprogrammation n. Par définition, la valeur du critère
au point n est l'espérance des prestures faites pour la valeur n. Cette
espérance varie dans le temps avec les caractéristiques de la charge du
système, mais son estimation sera possible si ces variations sont lentes.
Nous rechercherons dans ce but un estimateur avec subli du passé lointain,
d'où l'utilisation du lissage exponentiel. Nous calculerons, par définition :

$$
(III-16) \qquad Y_n^i = \alpha X_n^i + (1-\alpha) Y_{n-1}^{i-1}
$$

où $\alpha$ est un coefficient constant compris entre 0 et 1 appelé coefficient
de lissage.

En posant par définition $Y_n^C = 0$, quel que soit n, la formule
(III-16) permet le calcul itératif de $Y_n^i$. Lions $E(Y_n^i)$ à $E(X_n^i)$ lorsque
$E(X_n^i)$ est stationnaire (ou peut être considérée comme tel) c'est-à-dire :

$$
(III-17) \quad E(X_n^1) = E(X_n^2) = \dots = E(X_n^i) = E(X_n)
$$

La formule (III-16) donne :

$$
E(Y_n^1) = \alpha E(X_n^1) + \alpha(1-\alpha) E(X_n^{1-1}) + \dots + \alpha(1-\alpha)^{1-1} E(X_n^1) \text{ soit}
$$

$$
E(Y_n^i) = \alpha E(X_n) [1 + (1-\alpha) + \dots + (1-\alpha)^{i-1}] \text{ soit}
$$

$$
E(Y_n^i) = \frac{\alpha E(X_n) [1 - (1-\alpha)^i]}{1 - (1-\alpha)} \text{ soit}
$$

$$
(III-18) \quad E(Y_n^i) = E(X_n) [1 - (1-\alpha)^i]
$$

posons

$$
(III-19) \quad \lambda_i = \frac{1}{1 - (1-\alpha)^i}
$$

Pour faciliter le calcul, $\frac{1}{\mu}$ s'obtient par récurrence pour :

$$
(\text{III-19'}) \quad \frac{1}{\lambda_i} = (1-\alpha) \cdot \frac{1}{\lambda_{i-1}} + \text{d avec } \frac{1}{\lambda_0} = 0 \quad (\lambda_0 = \infty)
$$
d'où :

$$ (III-20) \quad E(Y_n^i) = \frac{1}{\lambda_i} E(X_n) $$

Si i est infini, $Y_n^i$ est un estimateur sans biais (cf annexe 2) de $E(X_n)$, estimateur avec oubli du passé (d'autant plus rapide que $\alpha$ est grand) en raison du facteur $(1-\alpha)$ dans la formule (III-16). Pour tenir compte du cas où i est faible (notamment au démarrage du système), nous remplacerons $Y_n^i$ par $Z_n^i$ défini par :

$$ (III-21) \quad Z_n^i = \lambda_i Y_n^i $$

D'après (III-20) et comme $E(Z_n^i) = \lambda_i E(Y_n^i)$ nous aurons $E(Z_n^i) = E(X_n)$. L'estimateur $Z_n^i$ estime donc sans biais notre critère quel que soit i.

La technique de recherche du degré optimal de multiprogrammation mettant en jeu des comparaisons de valeurs moyennes du critère pour différentes valeurs de $n$, il est utile de faire ces comparaisons avec une confiance statistique déterminée. On limite de cette façon les effets du bruit (variations aléatoires dans le temps). Nous sommes amenés ainsi à estimer une variance pour notre estimateur $Z_n^i$, cette nouvelle estimation oubliant elle-même le passé l'oinctain. En négligeant la corrélation entre $X_n^i$ et $Y_n^{i-1}$, l'équation (III-16) donne :

$$ (III-22) \quad V(Y_n^i) = \alpha^2 V(X_n^i) + (1-\alpha)^2 V(Y_n^{i-1}) $$

qui fournit $V(Y_n^i)$ par récurrence si $V(X_n^i)$ est connu. Nous sommes ramenés à l'étude de $V(X_n^i)$. Or :

$$ (III-23) \quad V(X_n^i) = E\left[\left(X_n^i\right)^2\right] - \left[E\left(X_n^i\right)\right]^2 $$

Il est naturel de chercher une estimation de $E\left[(X_n^i)^2\right]$ par un nouveau lissage exponentiel sur $(X_n^i)^2$. Nous définirons de même que précé-
donnent :

$$ (III-24) \qquad u_n^i = \beta (x_n^i)^2 + (1-\beta) u_n^{i-1} \qquad 0 < \beta < 1 $$

où $\beta$ est un nouveau coefficient de lissage

$$ (III-25) \qquad \mu_i = \frac{1}{1 - (1-\beta)^i} \quad \text{et} $$

$$ (III-26) \qquad v_n^i = \mu_i u_n^i $$

Nous avons bien (de la même façon que pour $z_n^i$ et $x_n^i$)

$$ (III-27) \qquad E(v_n^i) = E \left[ (x_n^i)^2 \right] $$

et $v_n^i$ estime sans biais $E [(x_n^i)^2]$.

Etudions maintenant le terme $[E(x_n^i)]^2$ dans l'équation (III-25). Comme $z_n^i$ estime sans biais $E(x_n^i)$, nous avons

$$ E(Z_n^i) = E(X_n^i) \quad \text{où} \quad [E(Z_n^i)]^2 = [E(X_n^i)]^2 $$

et compte tenu de (III-27), l'équation (III-23) devient :

$$ (III-28) \qquad V(X_n^i) = E(Y_n^i) - [E(Z_n^i)]^2 = E(Y_n^i) - \lambda_i^2 [E(Y_n^i)]^2 $$

Or $V(Y_n^i) = E [(Y_n^i)^2] - [E(Y_n^i)]^2$, formule qui peut s'écrire :

$$ [E(Y_n^i)]^2 = E[(Y_n^i)^2] - V(Y_n^i) $$

et la formule (III-28) devient :

$$ (III-29) \qquad V(X_n^i) = E(Y_n^i) - \lambda_i^2 E[(Y_n^i)^2] + \lambda_i^2 V(Y_n^i) $$
qui avec l'équation (III-22) forme un système de deux équations à deux inconnues $f(Y_n^i)$ et $V(X_n^i)$. Résolvant en $V(Y_n^i)$ il vient (en remplaçant $Y_n^i$ par $\mu_i u_n^i$ d'après (III-26)) :

$$ (\text{III-30}) \quad V(Y_n^i) = \frac{1}{1-\alpha^2 \lambda_i^2} \left\{ \alpha^2 \left[ \mu_i E(u_n^i) - \lambda_i^2 E\left[(Y_n^i)^2\right] \right] + (1-\alpha)^2 V(Y_n^{i-1}) \right\} $$

De la formule (III-30) et compte tenu que $V(Z_n^i) = \lambda_i^2 V(Y_n^i)$, nous déduisons deux estimations sans biais $\hat{V}(Y_n^i)$ et $\hat{V}(Z_n^i)$ pour $V(Y_n^i)$ et $V(Z_n^i)$ soit :

$$ (\text{III-31}) \quad (Y_n^i) = \frac{1}{1-\alpha^2 \lambda_i^2} \left\{ \alpha^2 \left[ \mu_i u_n^i - \lambda_i^2 (Y_n^i)^2 \right] + (1-\alpha)^2 \hat{V}(Y_n^{i-1}) \right\} \text{ et} $$

$$ (\text{III-32}) \quad \hat{V}(Z_n^i) = \lambda_i^2 \hat{V}(Y_n^i) $$

La variance de notre estimateur $Z_n^i$ ayant été calculée, nous effectuerons tous les calculs de précision nécessaires en faisant l'hypothèse que $u_n^i$ suit approximativement une loi gaussienne (cette affirmation résulte d'ailleurs du théorème central limite si toutes les mesures sont considérées comme indépendantes).

**Calculs pratiques** : A chaque nouvelle mesure $X_n^i$ obtenue pour le degré n de multiprogrammation, on entretient le nombre i de mesures effectuées pour le degré n, $Y_i^i$ d'après (III-16), $u_n^i$ d'après (III-24) et $\hat{V}(Y_n^i)$ d'après (III-31). Chaque fois que nous aurons à comparer les résultats obtenus pour deux valeurs diérentes de n, nous calculerons $Z_n^i$ d'après (III-21) et $\hat{V}(Z_n^i)$ d'après (III-32). Dans la pratique, dans les formules (III-16) et (III-24) nous faisons $\beta = \alpha$, c'est-à-dire que nous prenons deux coefficients de lissage identiques.

Remarque sur la variance $V(Y_n^i)$ :

La formule (III-22) permet d'obtenir une relation entre $V(Y_n^i)$
et $V(X_n^i)$ si $V(X_n^i)$ est stationnaire, c'est-à-dire si :

$$ (III-33) \quad V(X_n^i) = V(X_{n-1}^i) = \dots = V(X_1^i) = V(X_n). $$

En effet :

$$ V(Y_n^i) = \alpha^2 V(X_n^i) + \alpha^2 (1-\alpha) V(X_{n-1}^i) + \dots + \alpha^2 (1-\alpha)^{2(i-1)} V(X_1^i) $$

soit :

$$ V(Y_n^i) = \alpha^2 V(X_n) \left[ 1 + (1-\alpha)^2 + \dots + (1-\alpha)^{2(i-1)} \right] \text{d'où} $$

$$ V(Y_n^i) = \alpha^2 V(X_n) \frac{1 - (1-\alpha)^{2i}}{1 - (1-\alpha)^2} \quad \text{d'où :} $$

$$ (III-34) \quad V(Y_n^i) = \alpha V(X_n) \frac{1 - (1-\alpha)^{2i}}{2-\alpha} $$

Si $i$ est très grand, nous avons :

$$ (III-35) \quad V(Y_n^\infty) = \frac{\alpha V(X_n)}{2-\alpha} $$

La formule (III-35) montre qu'un $V(Y_n^i)$ faible demande un $\alpha$ faible, d'où un oubli faible ; autrement dit, les estimations sont d'autant moins dispersées que les mesures prises en compte (ici le passé) sont plus nombreuses, ce qui s'explique.

### Mesure de l'oubli :

L'oubli se mesure par ce que nous appellerons le seuil $\alpha$ : importance relative du passé dans la moyenne exponentielle après f newelles mesures. Par définition, le nombre $f$ sera nommé fenêtre. En supposant que pour les $f$ dernières mesures :
$$ (III-36) \quad E(X_n^i) = E(X_n^{i-1}) = \dots = E(X_n^{i+f+1}) = E(X_n) $$

la somme des termes relatifs aux f dernières mesures dans $Y_n^2$ (ou $Z_n^1$ si grand) vaut (après sommation) $E(X_n)$ $[1 - (1-\alpha)^f]$. L'importance relative du passé ancien vaut donc $(1-\alpha)^f$ et nous aurons par définition du seuil s :

$$ (III-37) \quad s = (1-\alpha)^{\frac{f}{n}} \quad \text{ou} \quad f = \frac{\text{Log } s}{\text{Log} (1-\alpha)} $$

Pour un seuil s fixé, la fenêtre f croît donc si $\alpha$ décroît. Le seuil s est un paramètre du modèle. Il doit être choisi de façon que l'importance relative du présent et du passé récent soit statistiquement significative pour prendre une décision.

## 4 - Description de l'optimiseur OPT

L'optimiseur OPT est appelé périodiquement par l'estimateur EST toutes les fois que f nouvelles mesures (f fenêtre d'oubli définie ci-dessus) ont été prises dans le système. OPT recherche le maximum du critère en fonction du degré de multiprogrammation n par comparaison des différentes estimations pour les diférentes valeurs de n. Définissons d'abord la technique de comparaison de deux points, avec une confiance statistique déterminée. Nous donnerons ensuite un algorithme de mise à jour de l'estimation $\bar{n}_0$.

### 4.1 - Comparaison de deux points

Soit à comparer les estimations $Z_{n1}$ et $Z_{n2}$ du critère pour les degrés de multiprogrammation $n_1$ et $n_2$, et $\hat{V}(Z_{n1})$ et $\hat{V}(Z_{n2})$. Les variances estimées de ces estimations (figure 9).
<figure><img src="image_21.png" /><figcaption>Figure 9</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 1</figcaption></figure>
En supposant toujours que $Z_{n_1}$ et $Z_{n_2}$ suivent approximativement des lois gaussiennes, la différence $Z_{n_1} - Z_{n_2}$ suit elle même une loi gaussienne dont la variance est estimée par $\hat{V}(Z_{n_1}) + \hat{V}(Z_{n_2})$ et l'écart-type $\sigma$ par $\hat{V}(Z_{n_1}) + \hat{V}(Z_{n_2})$. Nous devrons situer le zéro avec une confiance statistique déterminée par rapport à l'espérance de cette loi. Après les tables de la loi de Gauss, le signe $Z_{n_1} - Z_{n_2}$ sera significatif avec une confiance de 84 % si $|Z_{n_1} - Z_{n_2}| \ge \sigma$. Nous noterons $Z_{n_1} < Z_{n_2}$ si $Z_{n_1} - Z_{n_2} \ge \sigma$, c'est-à-dire si $Z_{n_1}$ est significativement inférieur à $Z_{n_2}$. Nous noterons aussi $Z_{n_1} > Z_{n_2}$ dans le cas où $Z_{n_1}$ est significativement supérieur à $Z_{n_2}$ et $Z_{n_1} \approx Z_{n_2}$ si la comparaison n'est pas significative.

Le but d'une telle technique de comparaison est le filtrage des fluctuations statistiques de l'agon à limiter les risques de mauvaises approximations.

## 4.2 - Algorithme d'estimation du degré optimal de multiPROGRAMMATION

L'algorithme envisagé doit à la fois effectuer une première estimation de l'optimum $n_0$ de n lors du démarrage du système et suivre les variations de cet optimum dans le temps. Dans le but de détecter un déplacement éventuel vers les n faibles de l'optimum $n_0$, et de ne faire inter-venir que des valeurs de n pour lesquelles des mesures ont été prises, l'ajustement de $\bar{n}_0$ se fait grâce aux comparaisons des performances aux points $\bar{n}_0$, $\bar{n}_0-1$ et $\bar{n}_0-2$, et la mise à jour consiste à faire +1 ou -1 sur l'ancienne estimation $\bar{n}_j$ (au démarrage, une première valeur (faible) de $\bar{n}_j$ est feurnie au système). Notons A, B, C les performances en $\bar{n}_0-2$, $\bar{n}_0-1$, $\bar{n}_0$ où $\bar{n}_0$ est pris à son ancienne valeur (figure 10). L'algorithme suivant a donné de bons résultats en simulation :

$$1 \rightarrow \text{Si } A \le C \text{ et si } B \le C \text{ alors } \bar{n}_0 := \bar{n}_0 + 1$$

$$\text{sinon } B > C \text{ alors } \bar{n}_0 := \bar{n}_0 - 1$$

$$\text{sinon } B \approx C \text{ alors } \bar{n}_0 := \bar{n}_0 + 1$$

$$2 \rightarrow \text{Si } A \le C \text{ alors } \bar{n}_0 := \bar{n}_0 - 1$$

$$3 \rightarrow \text{Si } A \approx C \text{ alors les mêmes actions que pour le cas 1.}$$
Cet algorithme a pour but :

- D'accroître l'estimation $\bar{n}_0$ lorsque la localité moyenne des programmes augmente (1'augmentation de l'optimum $n_0$ dans le temps s'accompagne de l'aplatissement des courbes).

- De diminuer $\bar{n}_0$ dès que les caractéristiques des programmes deviennent noirs favorables (l'estimation des performances en $\bar{n}_0$ décroît alors que les estimations en $\bar{n}_0-2$ et $\bar{n}_0-1$ sont inchangées : la comparaison donnera alors $A > C$ ou $B > C$).

Diverses précautions doivent être prises pour éviter les situations d'un manque trusque de programmes dans le système, s'ét-à-dire, un trou de charge. En effet, si le nombre de clients venus des consoliers a décroît brusquement, la commande K ne pourra plus maintenir dans le système le degré de multiprogrammation $\bar{n}_0$. Dans ce cas, nous pourrions ne plus disposer de nouvelles mesures aux points $\bar{n}_0, \bar{n}_0-1$ et $\bar{n}_0-2$. Les estimations du critère en ces trois points ne seront plus mises à jour et l'algorithme restera bloqué. Pour éviter cet inconvénient, on calcule la moyenne (entière par excès) p du degré de multiprogrammation entre deux appels de OFT et on remplace $\bar{n}_0$ par $\min(\bar{n}_0, p)$ avant d'appliquer l'algorithme. Le la sorte, l'estimation $\bar{n}_j$ sera au plus égale au nombre de clients demandant à entrer dans le système, ce qui n'est nullement un frein.

## 5 - Conclusion

Nous avons décrit une méthode de régulation possible d'un système informatique lorsque les caractéristiques des programmes varient dans le temps. Il nous reste pour en faire l'application à choisir un bon critère d'optimisation, puis à étudier par simulation le comportement du système régulé. Nous étudierons d'abord le temps nécessaires aux estimations de l'optimum, puis nous c mparerons les performances avec et sans régulation.
Si le système contient plusieurs programmes, la somme de tous les services peut devenir supérieure au temps total en raison du parallélisme dans le système, c'est-à-dire :

$$ (III-36) \quad \rho = \frac{\sum_{j=1}^{n} \sum_{k=1}^{K} t_{s_j}^{k}}{t_2 - t_1} \ge 1 $$

Dans le cas d'un recouvrement parfait de l'utilisation des ressources :

$$ \rho = \frac{\sum_{j=1}^{n} \sum_{k=1}^{K} t_{s_j}^{k}}{t_2 - t_1} = K $$

Le rapport $\rho$ défini par (III-36) sera appelé par définition dilatation globale. Sa valeur maximum théorique est de K, nombre de ressources. La dilatation globale est une mesure du gain apporté par la multiPROGRAMMATION, en raison du parallélisme dans le système.

Soit $\tau_k$ le taux d'utilisation de la kième ressource entre les instants $t_1$ et $t_2$. Par définition :

$$ (III-37) \quad \tau_k = \frac{\sum_{j=1}^{n} t_{s_j}^{k}}{t_2 - t_1} \quad \text{d'où :} $$

$$ (III-38) \quad \rho = \frac{\sum_{j=1}^{n} \sum_{k=1}^{K} t_{s_j}^{k}}{t_2 - t_1} = \sum_{k=1}^{K} \tau_k $$

La formule (III-38) montre que la dilatation globale est la somme des taux d'utilisation des ressources du système (unité centrale et disque de fichiers dans notre modèle).

La dilatation globale est liée physiquement aux temps de réponse des programmes. Nous définirons pour le programme j l'efficacité partielle entre $t_1$ et $t_2$ par
$$ (III-39) \quad \text{ep}_j = \frac{\sum_{k=1}^{K} t_{s_j}^k}{t_2 - t_1} $$

ep_j est par définition le rapport de la somme des services donnés au jième programme entre t_1 et t_2 au temps total. En particulier, si t_1 et t_2 sont les instants d'entrée et de sortie du programme du système, le somme $\sum_{k=1}^{K} t_{s_j}^k$ sera appelée par définition temps de réponse théorique TRH_j du programme, ep_j deviendra e_j : efficacité globale du programme par définition (efficacité partielle durant tout son séjour dans le système) et TR_j le temps de réponse physique du programme (III-39) devient

$$ (III-40) \quad e_j = \frac{\text{TRH}_j}{\text{TR}_j} $$

e_j est le rapport du temps de réponse théorique (durée totale des services demandés) au temps de réponse physique.

La dilatation globale entre t_1 et t_2 est aussi la somme des efficacités partielles des programmes entre t_1 et t_2 car

$$ (III-41) \quad \rho = \frac{\sum_{j=1}^{n} \sum_{k=1}^{K} t_{s_j}^k}{t_2 - t_1} = \sum_{j=1}^{n} \text{ep}_j $$

ce qui montre le lien entre dilatation globale et efficacités, d'où physiquement les temps de réponse.

La dilatation globale a été étudiée par simulation en fonction du degré n de multiprogrammation et différentes caractéristiques de programmes (figures 1 et 12). C'est un critère peu dispersé statistiquement (de même que les taux d'utilisation dont elle dérive) et à un seul maximum en n. La mesure fréquente est possible sur horloge. Ce critère tient compte de l'ensemble du système. Ces qualités font retenir son choix pour la régulation adaptative envisagée.
<figure><img src="image_23.png" /><figcaption>Figure 11</figcaption></figure>

Légende de multiprogramme.
<figure><img src="image_24.png" /><figcaption>Figure 12</figcaption></figure>

Moyenne

écart-type

autour de la moyenne

H = 512 pages

k = 2

r = 2 ms

t = 5 ms

t1F = 5 ms
# 2 - Etude de la régulation par simulation

## 2.1 - Hypothèses de fonctionnement

Pour étudier le comportement de la régulation dans le temps, nous étudierons la façon dont l'estimation $\bar{n}_o$ de l'optimum $n_o$ de n varie et atteint la bonne valeur, d'une part au démarrage du système, d'autre part lorsque les caractéristiques moyennes des programmes varient brusquement. Nous ferons l'hypothèse que pendant une période de temps déterminée, les consoles génèrent des programmes tous semblables statistiquement (même localité et même taux d'entrées-sorties) et nous ferons périodiquement varier brusquement ces caractéristiques. Nous avons défini à cet effet deux charges types que nous ferons alterner dans le temps :

- Une charge I avec une localité de 1.5 (pessimiste) et un temps moyen entre entrées-sorties de 20 ms.

- Une charge II avec une localité de 2 (optimiste) et un temps moyen entre entrées-sorties de 30 ms.

Nous ferons sur le système les hypothèses complémentaires suivantes :

- La taille sera prise à 512 pages dans ces expériences.

- Le temps superviseur nécessaire au fonctionnement de la régulation est compté à raison d'un temps d'activité moniteur fixe de 0.5 ms lors de chaque activation d'un module quelconque (EST, OPT, K) de la régulation.

- Les mesures de la dilatation globale (donc les activations de EST) se font, soit lors des variations de n (entrée ou sortie de programmes), soit sur horloge de garde au bout d'un temps maximum. La période P de l'horloge est un paramètre du modèle.

- Les consoles, au nombre de 20, ont été prises très rapides (temps moyen de réflexion de 4 secondes) afin de maintenir une charge très forte.
- Les temps de service moyens à la mémoire secondaire et au disque de fichiers ont été pris respectivement à 5 ms (tambour rapide) et 30 ms.

Pour l'étude du comportement de la régulation dans le temps, nous appliquerons la charge suivante :

- La charge II (favorable) pendant 200 secondes.

- La charge I (défavorable) pendant 300 secondes.

- De nouveau la charge II pendant 200 secondes pour un total de 700 secondes simulées.

La première estimation $\bar{n}_0$ de l'optimum $n_0$ a d'autre part été prise égale à 3 lors du démarrage du système.

Les outils nécessaires à l'étude de l'adaptation (évolution de l'estimation $\bar{n}_0$) sont d'une part la courbe de $\bar{n}_0$ en fonction du temps (courbes des figures 13 à 16 établies pour diverses valeurs des paramètres), d'autre part, dans un but d'explication des phénomènes observés ainsi que pour la mise au point, la trace, à chaque appel de OPT, des estimations tenues à jour par EST.

Nous allons d'abord chercher à optimiser les temps d'adaptation (temps nécessaires pour que l'estimation $\bar{n}_0$ atteigne la bonne valeur après chaque variation de charge) ; des temps courts sont en effet nécessaires à la mise en oeuvre de la régulation dans un système réel. Dans un second temps, nous montrerons l'intérêt de la commande par l'amélioration des performances moyennes obtenues.

## 2.2 - Etude des temps d'adaptation. Réglage des paramètres

Les paramètres de la régulation, sur lesquels on peut agir sont la période P de l'horloge de garde, le coefficient $\alpha$ de lissage exponentiel, et le seuil d'oubli s, défini par (III-37), qui commande la périodicité des appels à OPT. Notons que dans la pratique, presque toutes les mesures sur le système (donc les activations de EST) se font sur l'horloge, donc avec une période P.
Les courbes des figures 13 à 16 montrent le processus d'adaptation (variations de $\bar{n}_0$ dans le temps) pour différents réglages des paramètres. Chaque courbe comporte un premier front de montée, correspondant à l'estimation du $n_0$ optimal pour la charge II au démarrage du système, puis un front de descente lors de l'application de la charge I, enfin un nouveau front de montée lors du retour à la charge II.

Il apparaît que la cause principale de mauvaise montée est la dispersion statistique des estimations. Ce phénomène se traduit par de nombreuses oscillations de la courbe de $\bar{n}_0$ en fonction du temps, elles-mêmes dues à de mauvaises actions de OPT. Ce phénomène gênant est illustré par la courbe 13 obtenue pour un mauvais réglage des paramètres, et se traduit dans ce cas extrême par le non fonctionnement de l'adaptation. Ce phénomène se produit surtout lors du premier front de montée, au démarrage du système, ce qu'explique la formule (III-34) qui montre que la variance des estimations décroît jusqu'à une limite finie lorsque le nombre de mesures augmente. Une autre cause de montée lente est l'appel trop rare de l'optimiseur OPT (comparaison des courbes successives 14, 15 et 16).

Nous devrons donc à la fois chercher à réduire les dispersions statistiques et à appeler plus souvent l'optimiseur.

La seule façon raisonnable de diminuer la dispersion consiste à diminuer le coefficient $\alpha$ de lissage exponentiel (d'après les formules (III-34) et (III-35)). Malheureusement, d'après (III-37), la diminution de $\alpha$ s'accompagne de l'augmentation de la fenêtre $f$, nombre de mesures nécessaire avant l'appel de OPT. Nous accroîtons donc le nombre des mesures en diminuant la période $P$ de l'horloge de garde (avec une augmentation de la dispersion puisque nous diminuerons la période de temps sur laquelle sont prises les mesures). Pratiquement, nous diminuerons simultanément $\alpha$ et $P$ (courbes (13), (14) et (15) montrant des résultats toujours meilleurs). Enfin (d'après (III-37)), nous diminuerons $f$, donc nous appellerons plus souvent l'optimiseur, en augmentant le seuil $s$ tant que le passé récent reste significatif statistiquement (courbe 16) (pratiquement, on augmente s de façon à conserver un bon front de descente lors du
<figure><img src="image_25.png" /><figcaption>Figure 13</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure 14</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure 15</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Figure 16</figcaption></figure>
passage de la charge II à la charge I).

Les meilleurs temps d'adaptation (courbe 16) ont été obtenus pour $\alpha = 0.05$, $P = 100$ ms et $s = 0.25$. De la sorte, 20 à 30 secondes sont déjà suffisantes pour atteindre une bonne plage de fonctionnement lors du démarrage du système (d'après la courbe 12 représentant la dilatation en fonction de n pour la charge II). Le front de descente n'est plus limité que par la sortie des programmes et le front de montée suivant est très rapide.

## 2.3 - Performances obtenues grâce à la régulation

Les résultats nous permettant de juger à posteriori l'intérêt de la régulation seront les taux d'utilisation des serveurs, le temps de réponse moyen et l'efficacité moyenne des programmes, et le degré moyen de multiprogrammation.

Des essais préalables avec et sans régulation ont été conduits avec la charge I, puis la charge II seules, afin d'examiner l'intérêt de la régulation. Les résultats sont résumés dans le tableau de la figure 17. Il apparaît qu'avec la charge II, très favorable à la multiprogrammation, la régulation est inutile. Néanmoins, la baisse de performances due à la régulation, traduite par une légère diminution du taux d'unité centrale sans temps moniteur, du taux d'utilisation disque, de l'efficacité moyenne des programmes, et une légère augmentation du temps de réponse moyen, peut être considérée comme négligeable. Cette baisse n'affecte en rien l'utilité de la régulation comme sécurité contre l'écroulement lorsque les caractéristiques des programmes sont moins optimistes. Avec la charge I, beaucoup moins favorable à la multiprogrammation, la régulation double pratiquement les taux d'utilisation disque et unité centrale sans temps moniteur (l'utilisation totale unité centrale e. l'absence de régulation met d'ailleurs en évidence le temps exagéré utilisé par le superviseur lors de l'écroulement). L'utilisation du tambour est très diminuée. Le temps de réponse moyen est pratiquement divisé par 2 et l'efficacité moyenne des programmes est presque doublée. Dans le cas de l'alternance des charges I et II
telle que décrite ci-dessus, la régulation apporte globalement des améliorations résumées par le tableau de la figure 18. La régulation augmente l'utilisation de l'unité centrale, en particulier en mode problèmes (sans les temps moniteurs) et du disque de fichiers. L'utilisation moyenne du tambour diminue. Le temps de réponse moyen et le degré de multiprogrammation sont diminués, l'efficacité moyenne des clients augmentée.

V - CONCLUSION

Nous avons présenté une méthode de régulation possible du degré de multiprogrammation dans un système à mémoire virtuelle. Cette méthode est utilisable en l'absence de contrôle d'allocation mémoire. La régulation est surtout une sécurité contre l'écroulement dû à une pagination excessive lorsque le degré de multiprogrammation est trop important. Le contrôle du nombre de programmes qui se partagent les ressources est effectué ici par la mesure de l'activité générale du système. Nous ne nous sommes intéressés qu'à l'influence des caractéristiques moyennes de l'ensemble des programmes. Des schémas plus complexes de régulation peuvent être réalisés en prenant en compte plusieurs classes de programmes. Cette extension fera l'objet de développements futurs. Les fonctions de la régulation décrite ici ont été simulées telles qu'elles seraient mises en oeuvre dans la réalité. Elles peuvent donc être mises en place dans un système réel.
<table><caption>Figure 17</caption><thead><tr><th></th><th>Paris régulation</th><th>Avec régulation</th></tr></thead><tbody><tr><td>Utilisation</td><td>.59</td><td>.72</td></tr><tr><td>Utilisation sans les temps majeurs</td><td>.35</td><td></td></tr><tr><td>Utilisation</td><td>.97</td><td>.99</td></tr><tr><td>Util. disque de finition</td><td>.53</td><td>.9</td></tr><tr><td>Temps de répense moyen</td><td>33.12</td><td>18.79</td></tr><tr><td>Temps moyen de multi-grammation</td><td>10.9</td><td>4.98</td></tr><tr><td>Efficacité m.y.nr.</td><td></td><td>.11</td></tr><tr><td>Utilisation</td><td>.78</td><td>.73</td></tr><tr><td>Utilisation sans les temps majeurs</td><td>.44</td><td>.55</td></tr><tr><td>Utilisation</td><td>.24</td><td>.11</td></tr><tr><td>Util. disque de finition</td><td>.99</td><td>.97</td></tr><tr><td>Temps de répense moyen</td><td>15.79</td><td>10.55</td></tr><tr><td>Temps moyen de multi-grammation</td><td>13.55</td><td>8.98</td></tr><tr><td>Efficacité m.y.nr.</td><td>.13</td><td>.12</td></tr></tbody></table>
<table><caption>Figure 18</caption><thead><tr><th></th><th>Sans<br>régulation</th><th>Avec<br>régulation</th></tr></thead><tbody><tr><td>Utilisation UC</td><td>0,65</td><td>0,89</td></tr><tr><td>Utilisation UC sans<br>les temps moniteur</td><td>0,49</td><td>0,68</td></tr><tr><td>Utilisation tambour de<br>pagination</td><td>0,54</td><td>0,45</td></tr><tr><td>Utilisation disque de<br>fichiers</td><td>0,67</td><td>0,88</td></tr><tr><td>Temps de réponse moyen<br>(secondes)</td><td>25,35</td><td>18,67</td></tr><tr><td>Degré moyen de multipro<br>grammation</td><td>14,17</td><td>11,08</td></tr><tr><td>Efficacité moyenne des<br>clients</td><td>0,09</td><td>0,19</td></tr></tbody></table>
Nous avons décrit ici deux exemples d'études de modèles de systèmes informatiques menées par simulation. Les deux exemples décrits montrent bien deux catégories importantes d'application de la simulation de modèles.

Le premier type d'application porte sur des modèles résolubles mathématiquement de manière exacte ou approchée. Dans le dernier cas, la simulation permet de contrôler la validité des approximations faites (cf. Chapitre II). La simulation s'affranchissant des hypothèses simplificatives nécessaires à l'écriture et à la résolution des équations, il est également possible de faire varier ces hypothèses le façon à déterminer dans quelles limites la résolution mathématique donne encore des résultats acceptables ; ce contrôle s'appelle test de robustesse de la résolution mathématique. Certaines comparaisons évoquées au chapitre III (paragraphe II-1-3) en donnent un exemple (cas de la généralisation des lois du modèle).

Un autre type d'application correspond aux cas où l'on ne sait pas résoudre, ou même écrire, les équations du modèle, en raison de la trop grande complexité du problème. Nous trouvons dans ce cas les études de phénomènes transitoires, et l'extension des mécanismes des modèles mathématiques (après lesquelles tout calcul même approché est impossible). Un exemple de ce dernier point est donné au chapitre III par l'introduction des temps moniteurs. Un autre exemple est fourni par tous les essais d'algorithmes de gestion dans les systèmes, particulièrement lorsqu'on étudie le comportement de l'algorithme (exemple de la régulation au chapitre III).

Dans un dernier cas, on ne sait même pas formaliser correctement le problème par des objets mathématiques et l'on ne sait que reproduire un fonctionnement physique possible. Ce sera le cas en informatique pour la simulation d'un système particulier ou d'un avant-projet de système
en cours de conception ou de développement. Dans ce cas, la simulation permet d'expérimenter sans avoir à construire le système en vraie grandeur. Ce cas, déjà signalé, sort du cadre de cette étude.

Dans les deux cas où résolution mathématique et simulation sont pratiquables simultanément (au moins avec des hypothèses particulières) on peut opposer entre elles ces deux méthodes. La résolution mathématique (éventuellement approchée) nous semble recommandable dans tout domaine d'application où sa validité a été vérifiée d'une façon acceptable. Outre la possibilité de décrire les phénomènes par des équations, elle offre l'avantage d'un coût de calcul numérique et de programmation bien moindre. Son domaine d'application s'étend à tous les cas où les résultats, après expérience, se sont révélés peu sensibles aux hypothèses faites (par exemple la forme des lois de probabilités dans de nombreux calculs). A côté de ses avantages de grande souplesse et de généralité, la simulation présente en effet l'inconvénient d'être une méthode coûteuse. Elle est coûteuse d'abord en temps machine dépensé, surtout si une certaine précision des résultats est nécessaire (cf chapitre I et II). Coûteuse aussi en temps de conception, programmation et mise au point. A côté des difficultés habituelles d'écriture d'un programme assez souple, afin d'obtenir la généralité souhaitée, la simulation par évèvements par exemple impose une autre conception de la programmation qui l'apparente à la programmation système, particulièrement pour les prises de mesures. Dans les cas extrêmes (simulation d'un système complet), le simulateur peut devenir aussi complexe que le système lui-même, avec malgré tout l'avantage de pouvoir être écrit en langage évolué.

La mise au point des programmes simulant des phénomènes aléatoires est délicate, si l'on ne peut se ramener à des cas déterministes ; les méthodes habituelles de mise au point ne pouvant être utilisées dans ce cas, c'est habituellement la validation du simulateur par d'autres résultats qui permet d'affirmer que le programme est au point.
Devant les difficultés d'écriture des simulateurs, le problème s'est posé de disposer de systèmes de description de modèles de systèmes informatiques, capables de couvrir une large gamme d'applications. Nous pouvons citer quelques projets de description de modèles par un graphe [41], [42], [43] où les objets du graphe peuvent permettre des contrôles de validité [47] ou des prises automatiques de mesures. A notre connaissance, ce problème n'a pas reçu à ce jour de solution correcte dans toute la généralité souhaitée.

Dans les applications pratiques, les méthodes mathématiques et la simulation seront des méthodes complémentaires.
BIBLIOGRAPHIE

[ 1 ] BADEL M., GELENBE E., LEROUDIER J., POTIER D. - "Adaptive optimization of the performance of a virtual memory computer".
IRIA International Workshop proceedings published by North-Holland Company - August 12-14, 1974.

[ 2 ] BEILNER H. - "Problems in calibrating and validating simulation models 2nd Seminar on Experimental simulation and solving of probability problems".
Liblice (CSSR) 28-5 au 1-6 1973.

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

[ 4 ] BELADY L.A. - "A study of replacement algorithms for a virtual storage computer".
IMB Systems Journal, Vol 5 n°2, 1966.

[ 5 ] BELADY L., KUEHNER C.J. - "Dynamic Space Sharing in Computer Systems".
Comm. ACM, Vol 12, n° 5, May 1969, 282-288.

[ 6 ] BHARUCHA-REID A.T. - "Elements of the theory of Markov Processes and their Applications".
Mc Graw-Hill, New York, 1960.

[ 7 ] BIRTWISTLE G. - "Notes on the SIMULA system classes".
Norwegian Computing Center Oslo 3 Norway. April 1969.

[ 8 ] BROWN R.G. - "Smoothing, forecasting and prediction of discrete time series".
Prentice Hall Inc Englewood Cliffs N.J. 1963.

[ 9 ] CONWAY R.W. - "Some tactical problems in digital simulation".
Management Science Vol 10, n°1, Oct 1963.

[10] COX D.R., LEWIS P.A.W. - "The statistical analysis of series of events".
Methuen Co Ltd London 1968.

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

[12] CRANE M.A., IGLEHART D.L. - "Simulating stable stochastic system I : general multiserver queues".
JACM Vol 21, n°1, Jan. 1974.
[13] DAHL O.J., MYHRHANG B., NYGAARD K. - "SIMULA 67 Common base language". Norwegian Computing Center Oslo 3 Norway, April 1969.

[14] DENNING P.J. - "The working-set model for program behaviour". CACM, May 1968.

[15] DENNING P.J. - "Virtual memory". Computing Surveys n°3, September 1970.

[16] DENNING P.J. - "Thrashing: Its causes and preventions". AFIPS Conference Proceedings, FJCC 33 (1968).

[17] DIMSDALE B., MARKOWITZ H. - "A description of the STMSCRIPT language". IBM Systems Journal 3-1, 1964.

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

[19] FELLER W. - "The parabolic differential Equations and the associated semi-groups of transformations". Ann. Math. 55, 468-519 (1952).

[20] FELLER W. - "Diffusion processes in one dimension". Trans. Som. Math. Soc. 77, 1-31 (1954).

[21] FORRESTER J.W. - "Industrial dynamics". M.I.T. Press Student's Edition 1968.

[22] 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, Feb. 1967.

[23] FISHMAN G.S. - "Statistical analysis for queuing simulations". Management Science, Vol 20, n°3, Nov. 1973.

[24] FISHMAN G.S., KIVIAT P.J. - "Spectral analysis of time series generated by simulation models". Rand Corp. Memorandum RM 4393-PR. Feb. 1965.

[25] FOURGEAUD C., FUCHS A. - "Statistiques". Dunod 1967.

[26] GAVER D.P., SHEDLER G.S. - "Processor utilization in multiprogramming systems via diffusion approximations". Opns. Res. 21, 569-576 (1973).

[27] LE GALL P. - "Convergences des simulations et applications aux réseaux téléphoniques". AFIRO Monographies de recherche opérationnelle n°7, 1968.
[28] GAUDEUL B. - "Evaluation des performances des systèmes par méthode "Hardware". Compte rendu des Journées d'Etudes sur les recherches en mesures, simulation, modélisation des systèmes informatiques. Toulouse, Juin 1973, IRIA.

[29] GELENBE E. - "Modèles de systèmes informatiques". Thèse de Doctorat d'Etat ès Sciences Mathématiques, Université Paris VI (1973).

[30] GELENBE E. - "On approximate computer system models". IRIA International Workshop proceedings published by North-Holland Company. August 12-14, 1974.

[31] IBM Corp. - "General Purpose Simulation System/360". User's Manual H2C D326-0 New York 1967,

[32] JACKSON J.R. - "Jobshop-like Queueing Systems". Management Science, 10, 1, October 1963, 131-142.

[33] KOBAYASHI H. - "Application of the diffusion approximation to queueing networks" : Part I - Equilibrium queue distributions. ACM-SIGME First Symposium on Measurement and Evaluation, 54-62, Palo-Alto, Calif. (1973).

[34] LEROUDDIER J. - "Une analyse de système". Thèse de Docteur de 3ème cycle. Grenoble, avril 1973.

[35] LEROUDDIER J., PARENT M. - "Quelques aspects de la modélisation par simulation à évènements discrets". IRIA-LABORIA. Rapport de recherche n°98. Février 1975.

[36] Mac DOUGALL M.H., Mac ALPINE J.S. - "Computer system simulation with ASPOL". Symposium on the simulation of computer systems. June 1973.

[37] Mac DOUGALL M.H. - "Computer system simulation : an introduction". Computing Surveys, Vol 2, n°3, Sept. 1970.

[38] MARKOWITZ H., HANONER E., KARR H. - "SIMSCRIPT". Prentice Hall Englewood Cliffs, N.J. 1964.

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

[40] NEWELL G.F. - "Applications of Queueing Theory". Chapman and Hall. Ltd. London 1971.
[41] NCE J.D., NUTT G.J. - "Macro E.nets for representation of parallel systems". IEE Transactions of computers, Vol 022, n°8, August 1973.

[42] NUTT G.J., - "The computer system representation problem". Symposium on the simulation of the computer systems. June 1973.

[43] PARENT M.R., - "Presentation of the control graph model". Colloques IRIA, Aspects théoriques et pratiques des systèmes d'exploitations, Avril 1974.

[44] PUGH A. - "DYNAMO User's Manual". MIT Press. Cambridge. Mass. 1963.

[45] RODRIGUEZ-ROSSEL J., DUPUY J.F. - "The evaluation of a time-sharing page demand system". AFIPS Spring Joint Computer Conference, 1972, 759-765.

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

[47] WATTERS L.J., VASIL'K M.V., - "A stochastic approach to test and checkout". Rand corporation, Report P 4486, Oct. 1970.
# ANNEXE 1

## DESCRIPTION SOMMAIRE DU LANGAGE SIMULA 67 [7], [13].

Le langage SIMULA 67 est une extension de l'ALGOL 60. De même que l'ALGOL 60, il possède la notion de structure de blocs, et les instructions ALGOL (affectations, IF...THEN...ELSE..., boucles FOR et WHILE, déclarations de types de variables, de procédures). Le langage SIMULA est donc un langage algorithmique universel. Les principales extensions par rapport à l'algol comprennent les notions d'objet, de classe d'objet et de coroutine qui font de SIMULA 67 un langage général puisant.

Un objet est un ensemble de données, éventuellement de procédures et de code exécutable, dont le format est décrit par une déclaration de classe. Les déclarations de classes peuvent se faire hiérarchiquement (toujours en arbre) suivant le schéma classe et sous-classes. Il est possible ainsi de décrire le format d'un type d'objet général, puis d'objets de plus en plus particuliers.

Un objet SIMULA 67 est pointé par un nouveau type de variables appelées références, déclarées par une directive REF (classe) ou REF(classe) ARRAY. Un objet formant lui-même un bloc (au sens ALGOL), on accèdera à ses données ou procédures par notation pointée (éventuellement en cascade si les variables d'un objet comportent elles-mêmes des références). à partir d'une référence à cet objet.

### Exemple :

```pascal
CLASS POINT (X, Y)
REAL X; REAL Y;
CLASS SEGMENT (P1, P2)
REF (POINT) P1; REF (POINT) P2;
BEGIN REAL XMLIEU, YMILIEU;
XMILIEU := (P1.X + P2.X)/2.
YMILIEU := (P2.Y + P2.Y)/2.
END;
```
```pascal
REF(SEGMENT) S; REF(POINT) A, B;

A: - NEW POINT (1., 2.);
B: - NEW POINT (3., 4.);
S: - NEW SEGMENT (A, B);
```

Ce petit programme décrit d'abord le format d'un objet "point" par ses deux coordonnées, puis celui d'un objet "segment". Les variables de ce dernier sont les références aux deux extrémités et les coordonnées du milieu du segment. L'effet des ordres NEW est le suivant :

Pour les deux premiers (création de A et B), un nouveau point est créé à chaque fois avec les coordonnées indiquées. Ces nouveaux points créés sont pointés respectivement par les références A et B. L'ordre de création du segment S, dont la classe de description contient du code exécutable, sera le suivant : un nouvel objet segment (pointé par S) sera créé à partir des références aux points A et B. Le code exécutable décrit par la classe s'exécutera alors pour l'objet nouvellement créé, ce qui ici calculera les coordonnées du milieu du segment (le code indiqué ici donne un exemple de notation pointée).

Nous accèderons ensuite aux coordonnées du milieu du segment S par S.XMLMILIEU et S.YMILIEU.

Un ordre élaboré (INSPECT) permet d'éviter les notations pointées lorsqu'on travaille continuellement sur le même objet. Le code exécutable d'un objet est considéré comme un sous programme de type coroutine. Son exécution, commencée lors de la création de l'objet, peut être interrompue, puis reprise grâce aux ordres DETACH et RESUME. Ces derniers ordres sont utilisés implicitement dans les opérations sur les objets processus.

Le langage SIMULA 67 est présenté ainsi comme un langage tout à fait général. Divers éléments (sous-programmes et classes) prédéclarés rendent cependant faciles les applications à la simulation. Ce sont notamment les manipulations de listes (classes d'objets lien et tête de liste prédéclarées avec leurs procédures de mise en oeuvre) et processus. Un objet
processus n'est utilisable qu'avec la mise en oeuvre de la gestion du temps simulé (le programme est alors préfixé par la directive SIMULATION) et donne son plein sens à la notion de coroutine exposée précédemment. Les ordres de manipulation des processus agissent sur un échéancier, ce qui fournit la possibilité de simulation par évènements. Un processus peut être actif, suspendu, passif ou terminé. Dans les deux premiers cas, il possède dans l'échéancier son heure de mise en activité à laquelle son exécution sera relancée en séquence. Les opérations sur les processus se font par les ordres ACTIVATE ou REACTIVATE pour la programmation ou la reprogrammation dans l'échéancier du réveil d'un processus, HOLD pour la suspension du processus pendant un temps simulé donné, PASSIVATE et CANCEL pour rendre passif le processus actif ou un processus suspendu respectivement, TERMINATE pour terminer un processus.

Lorsqu'un processus se suspend ou devient passif, le contrôle est donné en séquence au prochain processus en tête d'échéancier et ainsi de suite.

Pour faciliter l'écriture des programmes de simulation, le langage SIMULA 67 possède en bibliothèque, outre des sous-programmes scientifiques ou utilitaires, des procédures de tirage aléatoire pour les lois de probabilité les plus usuelles.

Les ordinateurs sur lesquels le langage SIMULA 67 est disponible sont à notre connaisance les IBM série 360, UNIVAC 1107/1108, CDC 6600 et CII 10070 et Iris 80.
# ANNEXE II

## RAPPELS ET COMPLEMENTS DE STATISTIQUE [25], [27]

### I - DEFINITIONS RELATIVES AUX ESTIMATIONS

Nous poserons par hypothèse l'existence d'un échantillon $\{X_1, ..., X_n\}$ de n observations d'une variable aléatoire X. Nous rechercherons une valeur approchée d'une grandeur théorique G liée à la loi de probabilité de X par une fonction $\varphi$ de l'échantillon d'observation. Par définition, la fonction $\varphi$ sera appelée estimateur de la grandeur théorique G à partir de l'échantillon et le résultat :

$Y = \varphi\{X_1, ..., X_n\}$ sera appelé estimation de la grandeur G. Dans le cas général, l'estimation Y sera elle-même une variable aléatoire, puisqu'elle dérive l'observations aléatoires. Une telle estimation n'est utilisable que si elle possède certaines des propriétés que nous allons définir :

**Définition 1** : Une estimation Y de la grandeur théorique G est dite sans biais si $E(Y) = G$, et-st-à-dire si l'estimation a pour espérance mathématique la grandeur recherchée.

Nous ne nous intéressons dans cette étude qu'à des estimations sans biais.

**Exemple** : La moyenne

(1) $\bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_i$ de l'échantillon considéré précédemment est une estimation sans biais de l'espérance $E(X)$ car

$$ (2) \quad E(\bar{X}) = \frac{1}{n} \sum_{i=1}^{n} E(X_i) = \frac{1}{n} \sum_{i=1}^{n} E(X) = E(X) $$
Un cas intéressant dans les estimations est l'utilisation d'estimateurs convergents que nous allons définir. Dans de tels cas, il est possible d'atteindre par accroissement de l'échantillon toute précision (cf. chapitre I) fixée à l'avance comme nous allons le voir. Nous allons définir deux types intéressants de convergence aléatoire:

**Définition 2** : Y tend vers $E(Y)$ en moyenne quadratique, ou, ce qui revient au même, $[Y - E(Y)]$ tend vers 0 si la variance $\sigma^2$ tend vers 0, lorsque la taille de l'échantillon, ou la longueur de la simulation, croît indéfiniment.

**Définition 3** : Y tend vers $E(Y)$ presque sûrement quand l'échantillon d'observation croît indéfiniment s'il y a une probabilité de 1 pour que $[Y - E(Y)]$ tende vers 0 au sens de l'analyse. L'un de ces modes de convergence n'entraîne pas obligatoirement l'autre. Le cas le plus facilement vérifié dans les modèles de simulation (cf chapitre I) est la convergence en moyenne quadratique. Dans ce cas, en supposant par exemple $E(Y) > 0$, l'inégalité de Bienaymé-Tchebycheff donne :

$$ (3) \quad \text{Prob} \left\{ |Y - E(Y)| \ge \delta \sigma \right\} \le \frac{1}{\delta^2} $$

$\frac{1}{\delta^2}$ est appelé seuil de confiance et l'intervalle $[Y - \delta \sigma, Y + \delta \sigma]$ est un intervalle de confiance pour $E(Y)$ (cf. chapitre I). La convergence en moyenne quadratique, en assurant $\sigma \to 0$, permet d'affirmer que l'amplitude de l'intervalle de confiance donné tend vers 0 lorsque le nombre d'observations (la longueur de la simulation) croît indéfiniment.

## II - RAPPELS ET COMPLÉMENTS DE FORMULES USUELLES

### Propriétés de l'espérance mathématique

L'espérance mathématique est linéaire, c'est-à-dire

$$ (4) \quad E(X + Y) = E(X) + E(Y) \text{ et } E(\lambda X) = \lambda E(X) $$
Pour deux variables X et Y indépendantes

$$ (5) \qquad E(X.Y) = E(X) \cdot E(Y) $$

### Variances et covariances :

Par définition, (6) $V(X) = E[(X - E(X))^2]$. Dans la pratique, on applique :

$$ (7) \qquad V(X) = E(X^2) - [E(X)]^2 $$

La covariance de deux variables X et Y se définit par

$$ (8) \qquad \mathrm{Cov}(X, Y) = E[(X - \mathrm{E}(X)).(Y - \mathrm{E}(Y))] $$

On a évidemment (d'après la définition) $\mathrm{Cov}(X, Y) = \mathrm{Cov}(Y, X)$. Dans la pratique on utilise :

$$ (9) \qquad \mathrm{Cov}(X, Y) = E(X.Y) - E(X) \cdot E(Y) $$

La variance est bilinéaire (linéaire en X et Y à la fois). Par exemple pour X :

$$ (10) \qquad \mathrm{Cov}(X + Y, Z) = \mathrm{Cov}(X, Z) + \mathrm{Cov}(Y, Z) $$

$$ \text{et (11) } \mathrm{Cov}(\lambda X, Y) = \lambda \mathrm{Cov}(X, Y) $$

formules qui résultent de (9).

Les covariances interviennent dans les calculs des variances de sommes de variables aléatoires.
Pour deux variables, on a (d'après (7)) :

$$V(X + Y) = E[(X + Y)^2] - [E(X + Y)]^2 \text{ soit}$$

$$V(X + Y) = E(X^2) + E(Y^2) + 2E(X.Y) - [E(X)]^2 - 2E(X).E(Y) - [E(Y)]^2$$

d'où en regroupant les termes et en appliquant (7) et (9)

$$(12) \qquad V(X + Y) = V(X) + V(Y) + 2 \operatorname{Cov}(X, Y)$$

Dans le cas particulier de deux variables X et Y indépendantes, d'après (9) et (5), on a $\operatorname{Cov}(X, Y) = 0$ et la formule (12) se ramène à :

$$(13) \qquad V(X + Y) = V(X) + V(Y)$$

Dans le cas d'une somme de n variables, la formule (12) se généralise par récurrence sur le nombre des variables (et compte tenu de la formule (10)).

On a :

$$(14) \qquad V\left(\sum_{i=1}^{n} X_i\right) = \sum_{i=1}^{n} V(X_i) + \sum_{\substack{i=1 \\ i \neq j}}^{n} \sum_{j=1}^{n} \operatorname{Cov}(X_i, X_j)$$

$$\text{ou bien compte tenu de } \operatorname{Cov}(X_i, X_j) = \operatorname{Cov}(X_j, X_i)$$

$$(14 \text{ bis}) \qquad V\left(\sum_{i=1}^{n} X_i\right) = \sum_{i=1}^{n} V(X_i) + 2 \sum_{i=1}^{n-1} \sum_{j=i+1}^{n} \operatorname{Cov}(X_i, Y_j)$$

Supposons en effet la formule vraie pour $n=1$ variables (pour deux variables, elle se réduit à la formule (12)). On a (d'après (12)) :

$$V\left(\sum_{i=1}^{n} X_i\right) = V\left(\sum_{i=1}^{n-1} X_i\right) + V(X_n) + 2\operatorname{Cov}(X_n, \sum_{i=1}^{n-1} X_i) \text{ soit}$$
(d'après le formule (14) à l'ordre $n-1$ et (10)) :

$$V\left(\sum_{i=1}^{n-1} V(X_i) + V(X_n) + \sum_{i=1}^{n-1} \sum_{j=1, j \neq i}^{n-1} \operatorname{Cov}(X_i, Y_j) + 2 \sum_{i=1}^{n-1} \operatorname{Cov}(X_i, X_i)\right)$$

et en réécrivant :

$$2 \sum_{i=1}^{n-1} \operatorname{Cov}(X_n, X_i) \text{ sous la forme } \sum_{i=1}^{n-1} [\operatorname{Cov}(X_n, X_i) + \operatorname{Cov}(X_i, X_n)]$$

on a aboutit bien à la formule (14) à l'ordre n.