Simulez un système à événements discrets

Les formules du chapitre 4 supposent des services exponentiels. Or le module de paiement de la billetterie a des temps de service presque constants : 20 ms, à peu de chose près, pour chaque commande. La formule $R = S/(1 - U)$ s’applique-t-elle encore ? Et si l’on ajoutait une règle métier compliquée, comme « abandon du panier après 3 secondes d’attente » ? Aucune formule simple ne répond à cela.

Quand les formules ne suffisent plus et que le système n’existe pas encore (ou qu’on ne peut pas le tester en vrai), on le simule. Dans ce chapitre, vous allez écrire un simulateur complet de file d’attente en une trentaine de lignes, puis apprendre à en tirer des résultats dignes de confiance.

Le principe : sauter d’événement en événement

Dans une file d’attente, l’état du système (le nombre de clients présents) ne change qu’à des instants précis : quand un client arrive et quand un client part. Entre deux de ces instants, il ne se passe rien d’intéressant. Une simulation à événements discrets (discrete event simulation) exploite cette propriété : au lieu d’avancer le temps milliseconde par milliseconde, elle saute directement au prochain événement.

L'état ne change qu'aux arrivées et aux départs (triangles). Entre deux, il est constant : le simulateur n'a aucune raison d'y passer du temps et saute directement à l'événement suivant.

Elle repose sur quatre ingrédients :

• une horloge de simulation, qui indique la date courante (rien à voir avec le temps réel d’exécution) ;
• un état : ici, la liste des clients présents et leurs dates d’arrivée ;
• un échéancier (event list) : la liste des événements futurs, triée par date ;
• des routines d’événement : pour chaque type d’événement, ce qu’il faut faire à l’état, et quels nouveaux événements programmer.

La boucle principale est toujours la même : retirer de l’échéancier l’événement le plus proche, avancer l’horloge à sa date, exécuter sa routine, recommencer.

La boucle d'un simulateur à événements discrets : l'échéancier fournit toujours l'événement le plus proche, et chaque routine y réinsère les événements qu'elle déclenche.

Pour une file à un serveur, les routines sont simples :

• arrivée : ajouter le client à la file ; si le serveur était libre, programmer sa fin de service ; programmer l’arrivée suivante ;
• départ : retirer le client servi et noter son temps de réponse ; si d’autres attendent, programmer la fin de service du suivant.

Le hasard, sous contrôle

Une simulation a besoin de tirer des temps entre arrivées et des temps de service aléatoires. On utilise un générateur pseudo-aléatoire : une suite de nombres qui ont l’air aléatoires, mais qui est entièrement déterminée par une valeur de départ, la graine (seed). C’est une bonne nouvelle : avec la même graine, la simulation redonne exactement les mêmes résultats, ce qui permet de déboguer et de reproduire une expérience.

Pour tirer selon une loi exponentielle de taux $\lambda$, on peut utiliser la méthode de l’inversion : si $u$ est uniforme sur $]0, 1[$, alors $-\ln(u)/\lambda$ suit une loi exponentielle de moyenne $1/\lambda$. En Python, random.expovariate(lam) le fait pour vous.

La transformation −ln(u)/λ envoie les petits u vers de longues durées et les u proches de 1 vers des durées courtes ; 5 000 tirages reproduisent bien la densité exponentielle de moyenne 1/λ = 25 ms.

Voici le simulateur complet d’une file M/M/1 FCFS, avec le paiement de la billetterie ($\lambda = 40$ commandes/s, $S = 20$ ms en moyenne) :

import heapq
import random
import statistics

def simuler_file(lam, s_moy, n_clients, graine):
    """Simule une file M/M/1 FCFS ; renvoie les temps de réponse dans l'ordre des départs."""
    rng = random.Random(graine)                       # générateur à graine : reproductible
    echeancier = [(rng.expovariate(lam), "arrivee")]  # tas de (date, type d'événement)
    file = []                                         # dates d'arrivée des clients présents
    reponses = []
    while len(reponses) < n_clients:
        horloge, evenement = heapq.heappop(echeancier)  # l'horloge saute à l'événement suivant
        if evenement == "arrivee":
            file.append(horloge)
            if len(file) == 1:                        # serveur libre : service immédiat
                heapq.heappush(echeancier, (horloge + rng.expovariate(1 / s_moy), "depart"))
            heapq.heappush(echeancier, (horloge + rng.expovariate(lam), "arrivee"))
        else:                                         # fin de service
            reponses.append(horloge - file.pop(0))
            if file:                                  # le client suivant entre en service
                heapq.heappush(echeancier, (horloge + rng.expovariate(1 / s_moy), "depart"))
    return reponses

r = simuler_file(lam=40, s_moy=0.020, n_clients=20_000, graine=1)
print(f"R moyen : {statistics.mean(r) * 1000:.1f} ms")
print(f"p95     : {statistics.quantiles(r, n=100)[94] * 1000:.0f} ms")

L’échéancier est un tas (module heapq), qui donne toujours l’événement de plus petite date en temps logarithmique. Cette exécution affiche un temps de réponse moyen de 106,8 ms et un p95 de 307 ms. La théorie du chapitre 4 prévoit 100 ms et 300 ms. Proche, mais pas identique : pourquoi ? Et d’ailleurs, 106,8 ms, c’est « égal » à 100 ms ou non ? Il nous manque un outil pour conclure.

Vérifier, puis valider

Avant d’interpréter le moindre chiffre, deux questions distinctes se posent.

• La vérification : le programme fait-il ce que je crois ? C’est une question de débogage. On la traite en testant des cas dont on connaît la réponse : une file M/M/1, justement, dont on connaît $R$, $Q$ et la loi du temps de réponse. On peut aussi faire tourner le simulateur pas à pas avec des durées fixées à la main, ou vérifier la loi de Little sur les sorties.
• La validation : le modèle ressemble-t-il au système réel ? C’est une question de modélisation. Les arrivées sont-elles vraiment poissoniennes ? Les services, exponentiels ? On compare les sorties du simulateur à des mesures du système, quand elles existent, ou à l’avis d’experts.

À retenir

Vérifier, c’est s’assurer que le programme implémente bien le modèle. Valider, c’est s’assurer que le modèle représente bien le système. Un simulateur parfaitement vérifié peut être invalide (mauvaises hypothèses), et on vérifie toujours avant de valider.

Le régime transitoire

Notre simulation démarre avec une file vide, ce qui n’est pas représentatif d’un système en régime de croisière. Les premiers clients trouvent un serveur libre et passent vite. Pour le voir, relançons une petite simulation de 200 clients (graine 3) : les 20 premiers ont un temps de réponse moyen d’environ 26 ms, les 180 suivants d’environ 109 ms.

Cette phase de démarrage s’appelle le régime transitoire (warm-up). Si on la garde, elle biaise la moyenne vers le bas. On la supprime en jetant les premières observations. Combien ? Il n’y a pas de règle universelle : on trace la moyenne glissante des temps de réponse et on coupe là où elle se stabilise, en gardant une marge.

La file démarre vide : les 20 premiers clients passent presque sans attendre (26 ms en moyenne, contre 109 ms ensuite). Les garder tirerait la moyenne vers le bas.

Attention

Le même problème existe pour les mesures réelles : caches froids, compilation à la volée, connexions à ouvrir. On y reviendra aux chapitres 7 et 10 : on ne mesure jamais sans phase de chauffe.

Un intervalle de confiance pour une simulation

Une exécution de simulation donne une valeur, qui dépend de la graine. Pour estimer la vraie moyenne et sa précision, il faut plusieurs valeurs indépendantes, puis un intervalle de confiance (sa formule est détaillée au chapitre 7). Or les temps de réponse successifs d’une même exécution ne sont pas indépendants : si un client a attendu longtemps, le suivant aussi, très probablement. On ne peut donc pas calculer un intervalle de confiance directement sur les 20 000 temps de réponse. Deux méthodes classiques contournent le problème.

Les réplications indépendantes

On lance $n$ exécutions avec des graines différentes, on supprime le transitoire de chacune, et on calcule une moyenne par exécution. Ces $n$ moyennes sont indépendantes, et on peut leur appliquer la formule de l’intervalle de confiance :

moyennes = []
for graine in range(10):
    r = simuler_file(40, 0.020, 21_000, graine)
    moyennes.append(statistics.mean(r[1_000:]))       # on jette les 1 000 premiers clients

m = statistics.mean(moyennes)
s = statistics.stdev(moyennes)
t = 2.262                                             # quantile t(0,975 ; 9 degrés de liberté)
demi = t * s / len(moyennes) ** 0.5
print(f"R = {m * 1000:.1f} ms ± {demi * 1000:.1f} ms (IC à 95 %, 10 réplications)")

On obtient $R = 102{,}9 \pm 5{,}0$ ms, soit l’intervalle $[97{,}9 \, ; \, 107{,}9]$ ms. Il contient la valeur théorique de 100 ms : le simulateur est vérifié sur ce cas. Remarquez l’écart entre les réplications : les moyennes vont de 94,6 à 119,6 ms. Une seule exécution, même de 20 000 clients, pouvait donc tomber à 20 % de la vérité.

Les moyennes par lots

Si le transitoire est long, le refaire à chaque réplication coûte cher. On lance alors une seule longue exécution, on supprime le transitoire une fois, et on découpe la suite en lots (batches) consécutifs assez longs pour que leurs moyennes soient à peu près indépendantes. Avec une exécution de 200 000 clients utiles découpée en 20 lots de 10 000, on obtient $99{,}1 \pm 3{,}8$ ms.

Chaque point est une moyenne presque indépendante des autres ; leur dispersion donne l'intervalle de confiance (bande verte). Les deux méthodes contiennent la valeur théorique de 100 ms.

Quand s’arrêter ?

Le critère d’arrêt d’une simulation n’est pas « quand j’ai assez attendu », mais « quand l’intervalle de confiance est assez étroit ». On fixe une précision cible (par exemple ±2 % de la moyenne) et on ajoute des réplications ou des lots jusqu’à l’atteindre.

Piège

Présenter le résultat d’une seule exécution de simulation comme « le » résultat. Une simulation est une expérience aléatoire : sans réplications et sans intervalle de confiance, on ne sait pas si deux configurations diffèrent vraiment ou si l’on observe seulement l’effet de la graine.

En entretien

« Une simulation à événements discrets fait sauter l’horloge d’un événement au suivant grâce à un échéancier trié. Pour en tirer un résultat fiable, je la vérifie sur un cas connu comme M/M/1, je supprime le régime transitoire et je donne un intervalle de confiance calculé sur des réplications indépendantes. »

En résumé

• Une simulation à événements discrets gère une horloge, un état, un échéancier d’événements triés par date et des routines par type d’événement ; l’horloge saute d’un événement au suivant.
• Le hasard vient d’un générateur pseudo-aléatoire à graine : même graine, mêmes résultats. Une loi exponentielle se tire par inversion, $-\ln(u)/\lambda$.
• Vérifier (le programme implémente le modèle) avant de valider (le modèle représente le système).
• Le régime transitoire biaise les résultats : on supprime les premières observations.
• Les observations d’une même exécution sont corrélées : on construit l’intervalle de confiance sur des réplications indépendantes ou des moyennes par lots.
• On arrête la simulation quand l’intervalle de confiance atteint la précision voulue.

Exercices

Exercice · Dérouler un échéancier

facile

Une file FCFS à un serveur est vide à $t = 0$. Les clients arrivent aux dates 1, 2 et 6 (en secondes), et leurs services durent respectivement 3, 2 et 1 s. Déroulez la simulation événement par événement, et donnez le temps de réponse de chaque client.

Voir la solution

Horloge	Événement	Action	Échéancier après
1	arrivée du client 1	serveur libre : service jusqu’à 1 + 3 = 4	arrivée 2 (t = 2), départ 1 (t = 4)
2	arrivée du client 2	serveur occupé : il attend	départ 1 (t = 4), arrivée 3 (t = 6)
4	départ du client 1	$R_1 = 4 - 1 = 3$ s ; le client 2 entre en service jusqu’à 6	départ 2 (t = 6), arrivée 3 (t = 6)
6	départ du client 2	$R_2 = 6 - 2 = 4$ s ; file vide	arrivée 3 (t = 6)
6	arrivée du client 3	serveur libre : service jusqu’à 7	départ 3 (t = 7)
7	départ du client 3	$R_3 = 7 - 6 = 1$ s	vide

Les temps de réponse sont 3 s, 4 s (dont 2 s d’attente) et 1 s. Remarquez que deux événements ont lieu à $t = 6$ : l’ordre de traitement des événements simultanés est une convention à fixer (ici, le départ d’abord), et il peut changer le résultat dans d’autres cas.

En code

Python

import heapq

arrivees, services = [1, 2, 6], [3, 2, 1]
# (date, priorité, type, client) : à date égale, le départ (priorité 0) passe avant l'arrivée
echeancier = [(a, 1, "arrivée", i) for i, a in enumerate(arrivees)]
heapq.heapify(echeancier)
file = []                                        # clients présents, dans l'ordre FCFS
while echeancier:
    horloge, _, evenement, client = heapq.heappop(echeancier)
    if evenement == "arrivée":
        file.append(client)
        if len(file) == 1:                       # serveur libre : service immédiat
            heapq.heappush(echeancier, (horloge + services[client], 0, "départ", client))
    else:
        file.pop(0)
        print(f"t = {horloge} : départ du client {client + 1}, R = {horloge - arrivees[client]} s")
        if file:                                 # le suivant entre en service
            suivant = file[0]
            heapq.heappush(echeancier, (horloge + services[suivant], 0, "départ", suivant))
# t = 4 : départ du client 1, R = 3 s
# t = 6 : départ du client 2, R = 4 s
# t = 7 : départ du client 3, R = 1 s

R

arrivees <- c(1, 2, 6); services <- c(3, 2, 1)
debut <- depart <- numeric(3)
for (i in 1:3) {
  # FCFS à un serveur : le service commence à l'arrivée, ou au départ du précédent
  libre <- if (i == 1) 0 else depart[i - 1]
  debut[i] <- max(arrivees[i], libre)
  depart[i] <- debut[i] + services[i]
}
data.frame(client = 1:3, debut, depart,
           attente = debut - arrivees, R = depart - arrivees)
#   client debut depart attente R
# 1      1     1      4       0 3
# 2      2     4      6       2 4
# 3      3     6      7       0 1

Exercice · Réplications et précision

moyen

Cinq réplications d’une simulation donnent les temps de réponse moyens suivants (en ms) : 52, 48, 55, 50, 45.

1. Donnez l’intervalle de confiance à 95 % de la moyenne (le quantile de Student à 4 degrés de liberté vaut 2,776).
2. On voulait une précision de ±2 ms. Environ combien de réplications faudrait-il ?

Voir la solution

1. Moyenne : $\bar{x} = 250/5 = 50$ ms. Écarts à la moyenne : 2, −2, 5, 0, −5 ; somme des carrés : 4 + 4 + 25 + 0 + 25 = 58 ; écart type $s = \sqrt{58/4} \approx 3{,}81$ ms. Demi-largeur : $2{,}776 \times 3{,}81/\sqrt{5} \approx 4{,}73$ ms. L’intervalle vaut environ $[45{,}3 \, ; \, 54{,}7]$ ms.
2. La demi-largeur décroît comme $1/\sqrt{n}$. En approchant le quantile par 2 (on s’attend à un $n$ assez grand), il faut $2 \times 3{,}81/\sqrt{n} \leq 2$, soit $\sqrt{n} \geq 3{,}81$ et $n \geq 14{,}5$. Il faut donc environ 15 réplications ; avec le quantile de Student à 14 degrés de liberté (2,145), la demi-largeur vaudrait $2{,}145 \times 3{,}81/\sqrt{15} \approx 2{,}1$ ms, et on ajouterait une ou deux réplications pour passer sous 2 ms.

En code

Python

import numpy as np
from scipy import stats

x = np.array([52, 48, 55, 50, 45])                 # moyennes des réplications (ms)
n, m, s = len(x), x.mean(), x.std(ddof=1)
demi = stats.t.ppf(0.975, n - 1) * s / np.sqrt(n)
print(m, round(s, 2), round(demi, 2), round(m - demi, 1), round(m + demi, 1))
# 50.0 3.81 4.73 45.3 54.7

# Plus petit nombre de réplications pour une demi-largeur <= 2 ms (même s)
for n in range(5, 30):
    if stats.t.ppf(0.975, n - 1) * s / np.sqrt(n) <= 2:
        break
print(n, round(1.96 * s / np.sqrt(15), 2), round(stats.t.ppf(0.975, 14) * s / np.sqrt(15), 2))
# 17 1.93 2.11 : ~15 avec l'approximation normale, 17 avec Student

R

x <- c(52, 48, 55, 50, 45)                         # moyennes des réplications (ms)
n <- length(x); m <- mean(x); s <- sd(x)
demi <- qt(0.975, n - 1) * s / sqrt(n)
c(m, s, demi, m - demi, m + demi)                  # 50  3.81  4.73  45.3  54.7
t.test(x)$conf.int                                 # même intervalle : 45.27 54.73

# Plus petit nombre de réplications pour une demi-largeur <= 2 ms (même s)
n_ok <- 5:30
n_ok <- n_ok[qt(0.975, n_ok - 1) * s / sqrt(n_ok) <= 2][1]
c(n_ok, qt(0.975, 14) * s / sqrt(15))              # 17  2.11

Exercice · Quand les formules ne suffisent plus

difficile

Le paiement a en réalité un temps de service constant de 20 ms : c’est une file M/D/1. Modifiez le simulateur pour qu’il accepte une fonction de tirage du temps de service, relancez 10 réplications à $\lambda = 40$ et comparez avec la formule de Pollaczek-Khinchine pour une file M/G/1 :

$$R = S + \frac{U \, S \, (1 + C_s^2)}{2(1 - U)}$$

où $C_s$ est le coefficient de variation du temps de service (écart type divisé par la moyenne).

Voir la solution

On remplace s_moy par une fonction tirer_service(rng), appelée aux deux endroits où l’on programme un départ : lambda rng: rng.expovariate(50) pour M/M/1, et lambda rng: 0.020 pour M/D/1.

Avec 10 réplications (graines 0 à 9, 1 000 clients de chauffe jetés sur 21 000), la simulation donne environ $60{,}7 \pm 1{,}8$ ms pour M/D/1, contre $102{,}9 \pm 5{,}0$ ms pour M/M/1.

La formule : pour un service constant, $C_s = 0$ et $R = 20 + \frac{0{,}8 \times 20 \times 1}{2 \times 0{,}2} = 20 + 40 = 60$ ms. Pour un service exponentiel, $C_s = 1$ et on retrouve $R = 20 + \frac{0{,}8 \times 20 \times 2}{0{,}4} = 100$ ms, soit $S/(1 - U)$.

Les deux intervalles contiennent la valeur théorique. La leçon : à utilisation égale, rendre le service régulier divise ici le temps de réponse par près de deux. C’est encore la variabilité qui crée l’attente.

En code

Python

import heapq, random, statistics

def simuler_file(lam, tirer_service, n_clients, graine):
    """File FCFS à un serveur, arrivées de Poisson, service tiré par tirer_service(rng)."""
    rng = random.Random(graine)
    echeancier, file, reponses = [(rng.expovariate(lam), "arrivee")], [], []
    while len(reponses) < n_clients:
        horloge, evenement = heapq.heappop(echeancier)
        if evenement == "arrivee":
            file.append(horloge)
            if len(file) == 1:
                heapq.heappush(echeancier, (horloge + tirer_service(rng), "depart"))
            heapq.heappush(echeancier, (horloge + rng.expovariate(lam), "arrivee"))
        else:
            reponses.append(horloge - file.pop(0))
            if file:
                heapq.heappush(echeancier, (horloge + tirer_service(rng), "depart"))
    return reponses

def intervalle(tirer_service):
    moy = [statistics.mean(simuler_file(40, tirer_service, 21_000, g)[1_000:]) for g in range(10)]
    return statistics.mean(moy) * 1000, 2.262 * statistics.stdev(moy) / 10 ** 0.5 * 1000

print("M/D/1 : %.1f ± %.1f ms" % intervalle(lambda rng: 0.020))              # 60.7 ± 1.8 ms
print("M/M/1 : %.1f ± %.1f ms" % intervalle(lambda rng: rng.expovariate(50)))  # 102.9 ± 5.0 ms
U, S = 0.8, 20                                                   # Pollaczek-Khinchine (ms)
print([round(S + U * S * (1 + cs2) / (2 * (1 - U)), 1) for cs2 in (0, 1)])  # [60.0, 100.0]

R

# File FCFS à un serveur par la récurrence de Lindley :
# départ_i = max(arrivée_i, départ_{i-1}) + service_i
simuler_file <- function(lam, tirer_service, n_clients) {
  arrivees <- cumsum(rexp(n_clients, lam))
  services <- tirer_service(n_clients)
  depart <- numeric(n_clients); libre <- 0
  for (i in 1:n_clients) {
    depart[i] <- max(arrivees[i], libre) + services[i]
    libre <- depart[i]
  }
  depart - arrivees                               # temps de réponse
}

intervalle <- function(tirer_service) {
  moy <- replicate(10, mean(simuler_file(40, tirer_service, 21000)[-(1:1000)]))
  1000 * c(moyenne = mean(moy), demi = qt(0.975, 9) * sd(moy) / sqrt(10))
}
set.seed(0)
intervalle(function(n) rep(0.020, n))             # M/D/1 : 60.6 ± 1.1 ms
intervalle(function(n) rexp(n, 50))               # M/M/1 : 99.0 ± 4.6 ms
U <- 0.8; S <- 20                                 # Pollaczek-Khinchine (ms)
S + U * S * (1 + c(0, 1)) / (2 * (1 - U))         # 60  100