Modélisez un processus de naissance et de mort

Dans une file d’attente, les clients arrivent un par un et repartent un par un. Dans un atelier, les machines tombent en panne une à une et sont réparées une à une. Dans une population, les individus naissent et meurent un à la fois. Tous ces systèmes ont un point commun : l’état est un nombre, et il ne bouge que de plus un ou de moins un.

Cette structure très particulière rend les calculs remarquablement simples : le régime stationnaire s’obtient par une formule explicite, sans résoudre de système. Ce chapitre vous la donne, et elle servira de moteur à toute la partie 3 sur les files d’attente. Notre fil rouge : un petit atelier de trois machines et un seul réparateur.

Le fil rouge : trois machines, un réparateur

Un atelier de fabrication possède trois machines identiques. Chacune, quand elle fonctionne, tombe en panne au taux $\lambda = 0{,}1$ par heure (10 heures de fonctionnement en moyenne). Un seul technicien les répare, une à la fois, en une durée exponentielle de taux $\mu = 0{,}5$ par heure (2 heures en moyenne). Les durées sont indépendantes.

On note $X_t$ le nombre de machines en panne à l’instant $t$ (qu’elles soient en réparation ou en attente). Les états possibles sont 0, 1, 2 et 3. Que peut-il se passer à partir de l’état $n$ ?

• Une panne fait passer de $n$ à $n + 1$. Il y a $3 - n$ machines en marche, chacune avec une horloge exponentielle de taux $\lambda$ ; la première panne survient au taux $(3 - n)\lambda$ (les taux s’additionnent, chapitre 5).
• Une réparation fait passer de $n$ à $n - 1$, au taux $\mu$ dès qu’au moins une machine est en panne (le technicien ne répare qu’une machine à la fois).

Deux événements simultanés ont une probabilité négligeable : l’état ne saute jamais de deux unités d’un coup.

Définition et générateur

Un processus de naissance et de mort est une chaîne de Markov à temps continu à valeurs dans $\mathbb N$ (ou dans $\lbrace 0, \dots, K \rbrace$) dont les seules transitions sont $n \to n+1$ et $n \to n-1$. Sur un petit intervalle de durée $h$ :

$$p_{n,n+1}(h) = \lambda_n h + o(h), \qquad p_{n,n-1}(h) = \mu_n h + o(h), \qquad p_{nm}(h) = o(h) \text{ si } |n - m| \geq 2.$$

• $\lambda_n$ est le taux de naissance dans l’état $n$ (arrivée d’un client, panne d’une machine) ;
• $\mu_n$ est le taux de mort dans l’état $n$ (départ d’un client, fin d’une réparation), avec $\mu_0 = 0$ puisqu’on ne peut pas descendre sous zéro.

Attention

La mort fait passer de $n$ à $n - 1$ : c’est bien $p_{n,n-1}(h) = \mu_n h + o(h)$, avec l’indice de départ en premier. Écrire $p_{n-1,n}$ reviendrait à décrire une naissance.

Le graphe des taux est une simple ligne :

      λ0        λ1        λ2               λ(n-1)       λn
  0 ─────► 1 ─────► 2 ─────► 3  ···  n-1 ───────► n ───────► n+1  ···
    ◄─────   ◄─────   ◄─────             ◄───────   ◄───────
      μ1        μ2        μ3                μn         μ(n+1)

et le générateur est tridiagonal : sur la ligne $n$, on trouve $\mu_n$ juste avant la diagonale, $\lambda_n$ juste après, et $-(\lambda_n + \mu_n)$ sur la diagonale.

$$A = \begin{pmatrix} -\lambda_0 & \lambda_0 & 0 & 0 & \cdots \\ \mu_1 & -(\lambda_1 + \mu_1) & \lambda_1 & 0 & \cdots \\ 0 & \mu_2 & -(\lambda_2 + \mu_2) & \lambda_2 & \cdots \\ \vdots & & \ddots & \ddots & \ddots \end{pmatrix}$$

Pour l’atelier : $\lambda_0 = 0{,}3$, $\lambda_1 = 0{,}2$, $\lambda_2 = 0{,}1$, $\lambda_3 = 0$ (plus aucune machine ne peut tomber en panne) et $\mu_1 = \mu_2 = \mu_3 = 0{,}5$.

Le graphe de l'atelier est une ligne : on monte d'un cran à chaque panne, à un taux qui baisse avec le nombre de machines encore en marche, et on descend d'un cran à chaque réparation, toujours au taux 0,5.

Le régime transitoire : des équations en chaîne

L’équation de Kolmogorov $p'(t) = p(t)\,A$ prend ici une forme très lisible :

$$p_0'(t) = -\lambda_0\, p_0(t) + \mu_1\, p_1(t)$$ $$p_n'(t) = \lambda_{n-1}\, p_{n-1}(t) - (\lambda_n + \mu_n)\, p_n(t) + \mu_{n+1}\, p_{n+1}(t) \quad (n \geq 1)$$

La probabilité d’être en $n$ augmente grâce aux naissances depuis $n - 1$ et aux morts depuis $n + 1$, et diminue de tout ce qui quitte $n$. Pour un espace infini, c’est un système infini d’équations différentielles couplées ; sauf cas particuliers (comme le processus de Poisson du chapitre suivant), on ne sait pas le résoudre à la main. En pratique, on le résout numériquement sur un espace tronqué, ou on se contente du régime stationnaire, beaucoup plus simple.

Le système différentiel de l'atelier, résolu numériquement : partant de 0 panne, les probabilités se stabilisent en une dizaine d'heures sur les valeurs stationnaires (pointillés), que la suite du chapitre obtient sans aucun calcul différentiel.

La balance locale et la formule du régime stationnaire

Écrivons la balance globale (flux sortant = flux entrant) de l’état 0 : on ne sort de 0 que par une naissance, on n’y entre que par une mort depuis 1.

$$\lambda_0\, \Pi_0 = \mu_1\, \Pi_1.$$

Pour l’état 1, la balance globale s’écrit $(\lambda_1 + \mu_1)\,\Pi_1 = \lambda_0\,\Pi_0 + \mu_2\,\Pi_2$. En y retranchant l’égalité précédente, il reste $\lambda_1\,\Pi_1 = \mu_2\,\Pi_2$. Par récurrence, on obtient pour tout $n$ :

$$\lambda_n\, \Pi_n = \mu_{n+1}\, \Pi_{n+1}.$$

C’est l’équation de balance locale. Elle a une interprétation géométrique simple : coupez le graphe entre $n$ et $n + 1$. Le seul moyen de franchir la coupure vers la droite est une naissance depuis $n$, et le seul moyen de la franchir vers la gauche est une mort depuis $n + 1$. En régime stationnaire, il passe autant de flux dans chaque sens.

En déroulant la balance locale, $\Pi_n = \dfrac{\lambda_{n-1}}{\mu_n}\,\Pi_{n-1} = \cdots$, d’où la formule générale :

$$\Pi_n = \gamma_n\, \Pi_0 \quad \text{avec} \quad \gamma_n = \frac{\lambda_0\, \lambda_1 \cdots \lambda_{n-1}}{\mu_1\, \mu_2 \cdots \mu_n}, \qquad \Pi_0 = \left(1 + \sum_{n \geq 1} \gamma_n\right)^{-1}.$$

• $\gamma_n$ est le produit des taux de naissance de 0 à $n - 1$, divisé par le produit des taux de mort de 1 à $n$ ;
• $\Pi_0$ vient de la normalisation $\sum_n \Pi_n = 1$.

À retenir

Pour un processus de naissance et de mort, le régime stationnaire se lit directement : $\Pi_n = \gamma_n\,\Pi_0$ avec $\gamma_n = \dfrac{\lambda_0 \cdots \lambda_{n-1}}{\mu_1 \cdots \mu_n}$, et $\Pi_0 = (1 + \sum_{n \geq 1} \gamma_n)^{-1}$. Il existe si et seulement si la série $\sum \gamma_n$ converge (condition toujours vérifiée si l’espace est fini).

Piège

Les indices ne sont pas alignés : en haut, les naissances vont de $\lambda_0$ à $\lambda_{n-1}$ ; en bas, les morts vont de $\mu_1$ à $\mu_n$. Pour vous en souvenir, pensez au chemin de 0 à $n$ : on monte $n$ marches, en utilisant les naissances des états $0, \dots, n-1$ et les morts des états d’arrivée $1, \dots, n$.

Exemple

Pour l’atelier ($\lambda = 0{,}1$, $\mu = 0{,}5$) :

• $\gamma_1 = \dfrac{0{,}3}{0{,}5} = 0{,}6$ ;
• $\gamma_2 = \dfrac{0{,}3 \times 0{,}2}{0{,}5^2} = 0{,}24$ ;
• $\gamma_3 = \dfrac{0{,}3 \times 0{,}2 \times 0{,}1}{0{,}5^3} = 0{,}048$.

Donc $\Pi_0 = 1/(1 + 0{,}6 + 0{,}24 + 0{,}048) = 1/1{,}888 \approx 0{,}530$, puis

$$\Pi \approx (0{,}530 \,;\, 0{,}318 \,;\, 0{,}127 \,;\, 0{,}025).$$

On en déduit les grandeurs utiles au chef d’atelier :

• nombre moyen de machines en panne : $\sum_n n\,\Pi_n \approx 0{,}318 + 2 \times 0{,}127 + 3 \times 0{,}025 \approx 0{,}648$ ;
• nombre moyen de machines en marche : $3 - 0{,}648 \approx 2{,}35$ ;
• taux d’occupation du technicien : $1 - \Pi_0 \approx 47\,\%$ ;
• probabilité que tout l’atelier soit à l’arrêt : $\Pi_3 \approx 2{,}5\,\%$.

La balance locale sur les chiffres de l'atelier : à chaque coupure, le flux qui monte (taux de panne × probabilité de l'état de gauche) égale le flux qui descend (0,5 × probabilité de l'état de droite).

Et quand l’espace d’états est infini ? Tout dépend de la convergence de $\sum \gamma_n$. Deux cas fondamentaux, que vous retrouverez dans la partie 3 :

• Taux constants $\lambda_n = \lambda$, $\mu_n = \mu$ : alors $\gamma_n = \rho^n$ avec $\rho = \lambda/\mu$. La série géométrique converge si et seulement si $\rho < 1$, et alors $\Pi_n = (1 - \rho)\,\rho^n$. C’est la file M/M/1. Si $\rho \geq 1$, il n’y a pas de régime stationnaire : la file grossit indéfiniment (ou erre sans jamais se stabiliser).
• Morts proportionnelles $\lambda_n = \lambda$, $\mu_n = n\,\mu$ : alors $\gamma_n = u^n / n!$ avec $u = \lambda/\mu$, la série converge toujours (vers $e^u$), et $\Pi_n = e^{-u}\,u^n/n!$ est une loi de Poisson. C’est la file M/M/∞.

Même rapport λ/μ = 0,8, deux comportements : avec un seul serveur (taux de mort constant), la loi géométrique a une longue queue et une moyenne de 4 ; quand chaque client a son serveur (morts proportionnelles à n), la loi de Poisson reste concentrée près de 0.

Voici une fonction générique, qui accepte les taux sous forme de fonctions de $n$ et tronque l’espace si nécessaire :

import numpy as np

def stationnaire_nm(taux_naissance, taux_mort, n_max):
    """Régime stationnaire d'un processus de naissance et de mort sur {0, ..., n_max}.
    taux_naissance(n) = lambda_n, taux_mort(n) = mu_n (pour n >= 1)."""
    gamma = np.ones(n_max + 1)
    for n in range(1, n_max + 1):
        gamma[n] = gamma[n - 1] * taux_naissance(n - 1) / taux_mort(n)   # balance locale
    return gamma / gamma.sum()

lam, mu = 0.1, 0.5
pi = stationnaire_nm(lambda n: (3 - n) * lam, lambda n: mu, n_max=3)
print(np.round(pi, 3))                    # [0.53  0.318 0.127 0.025]
print((np.arange(4) * pi).sum())          # nombre moyen de machines en panne ≈ 0.648

# M/M/1 tronquée : si rho < 1, la troncature a un effet négligeable
pi_mm1 = stationnaire_nm(lambda n: 8.0, lambda n: 10.0, n_max=200)
print(pi_mm1[:3])                         # ≈ [0.2  0.16  0.128] = (1 - rho) rho^n

En entretien

« Dans un processus de naissance et de mort, on ne passe que de $n$ à $n \pm 1$. En coupant le graphe entre $n$ et $n+1$, le flux montant égale le flux descendant : $\lambda_n \Pi_n = \mu_{n+1}\Pi_{n+1}$. Cela donne le régime stationnaire par un simple produit de rapports de taux, à condition que la série converge. »

En résumé

• Un processus de naissance et de mort est une chaîne à temps continu sur les entiers qui ne saute que de $+1$ (taux $\lambda_n$) ou de $-1$ (taux $\mu_n$, avec $\mu_0 = 0$).
• Son générateur est tridiagonal ; les équations du régime transitoire couplent chaque $p_n$ à ses deux voisins et sont en général difficiles à résoudre.
• En régime stationnaire, la balance locale $\lambda_n \Pi_n = \mu_{n+1}\Pi_{n+1}$ exprime l’égalité des flux à travers la coupure entre $n$ et $n+1$.
• Formule générale : $\Pi_n = \gamma_n \Pi_0$, $\gamma_n = \dfrac{\lambda_0\cdots\lambda_{n-1}}{\mu_1\cdots\mu_n}$, $\Pi_0 = (1 + \sum_{n\geq1}\gamma_n)^{-1}$, si la série converge.
• Taux constants : loi géométrique $\Pi_n = (1-\rho)\rho^n$ si $\rho = \lambda/\mu < 1$. Morts proportionnelles à $n$ : loi de Poisson de paramètre $\lambda/\mu$.

Exercices

Exercice · Un calcul direct

facile

Un processus de naissance et de mort sur $\lbrace 0, 1, 2 \rbrace$ a pour taux $\lambda_0 = 2$, $\lambda_1 = 1$, $\mu_1 = 1$, $\mu_2 = 2$. Écrivez son générateur et calculez sa loi stationnaire.

Voir la solution

$$A = \begin{pmatrix} -2 & 2 & 0 \\ 1 & -2 & 1 \\ 0 & 2 & -2 \end{pmatrix}$$

(chaque ligne somme à 0). On a $\gamma_1 = \lambda_0/\mu_1 = 2$ et $\gamma_2 = \lambda_0\lambda_1/(\mu_1\mu_2) = 2/2 = 1$. Donc $\Pi_0 = 1/(1 + 2 + 1) = 1/4$ et $\Pi = (\tfrac14, \tfrac12, \tfrac14)$. On vérifie la balance locale : $\lambda_0\Pi_0 = 2 \times \tfrac14 = \tfrac12 = \mu_1\Pi_1$ et $\lambda_1\Pi_1 = \tfrac12 = \mu_2\Pi_2$.

En code

Python

import numpy as np

A = np.array([[-2, 2, 0],
              [1, -2, 1],
              [0, 2, -2]], dtype=float)
print(A.sum(axis=1))                     # [0. 0. 0.]

# Formule produit : gamma_n = (lambda_0 ... lambda_{n-1}) / (mu_1 ... mu_n)
lam = [2, 1]
mu = [1, 2]
gamma = np.cumprod([1] + [l / m for l, m in zip(lam, mu)])
pi = gamma / gamma.sum()
print(pi)                                # [0.25 0.5  0.25]
print(pi @ A)                            # [0. 0. 0.] : Pi A = 0

R

A <- matrix(c(-2,  2,  0,
               1, -2,  1,
               0,  2, -2), nrow = 3, byrow = TRUE)
rowSums(A)                              # 0 0 0

# Formule produit : gamma_n = (lambda_0 ... lambda_{n-1}) / (mu_1 ... mu_n)
lam <- c(2, 1); mu <- c(1, 2)
gamma <- cumprod(c(1, lam / mu))
pi_ <- gamma / sum(gamma)
pi_                                     # 0.25 0.50 0.25
pi_ %*% A                               # 0 0 0 : Pi A = 0

Exercice · Embaucher un second technicien ?

moyen

Dans l’atelier du chapitre, on envisage d’embaucher un second technicien : deux machines peuvent alors être réparées en même temps.

1. Quels taux changent ? Donnez les nouveaux $\mu_n$.
2. Calculez la nouvelle loi stationnaire et le nombre moyen de machines en panne.
3. Le gain vous paraît-il justifier l’embauche ?

Voir la solution

1. Les taux de panne ne changent pas. Les réparations se font au taux $\mu_n = \min(n, 2)\,\mu$ : $\mu_1 = 0{,}5$, $\mu_2 = \mu_3 = 1$.
2. $\gamma_1 = 0{,}3/0{,}5 = 0{,}6$ ; $\gamma_2 = 0{,}3 \times 0{,}2/(0{,}5 \times 1) = 0{,}12$ ; $\gamma_3 = 0{,}3 \times 0{,}2 \times 0{,}1/(0{,}5 \times 1 \times 1) = 0{,}012$. Donc $\Pi_0 = 1/1{,}732 \approx 0{,}577$ et

$$\Pi \approx (0{,}577 \,;\, 0{,}346 \,;\, 0{,}069 \,;\, 0{,}007).$$

Le nombre moyen de machines en panne vaut $0{,}346 + 2 \times 0{,}069 + 3 \times 0{,}007 \approx 0{,}506$, contre 0,648 avec un seul technicien. 3. On gagne en moyenne $0{,}648 - 0{,}506 \approx 0{,}14$ machine en marche, soit environ 6 % de production en plus ($2{,}49$ machines en marche au lieu de $2{,}35$). La probabilité d’arrêt total passe de 2,5 % à 0,7 %. Le premier technicien n’est occupé que 47 % du temps : l’embauche ne se justifie que si une machine arrêtée coûte très cher.

En code

Python

import numpy as np

def stationnaire_nm(taux_naissance, taux_mort, n_max):
    gamma = np.ones(n_max + 1)
    for n in range(1, n_max + 1):
        gamma[n] = gamma[n - 1] * taux_naissance(n - 1) / taux_mort(n)
    return gamma / gamma.sum()

lam, mu = 0.1, 0.5
panne = lambda n: (3 - n) * lam
for techniciens in (1, 2):
    pi = stationnaire_nm(panne, lambda n: min(n, techniciens) * mu, 3)
    en_panne = (np.arange(4) * pi).sum()
    print(techniciens, np.round(pi, 3), round(en_panne, 3), round(3 - en_panne, 2))
# 1 [0.53  0.318 0.127 0.025] 0.648 2.35
# 2 [0.577 0.346 0.069 0.007] 0.506 2.49

R

stationnaire_nm <- function(taux_naissance, taux_mort, n_max) {
  gamma <- cumprod(c(1, sapply(1:n_max, function(n)
    taux_naissance(n - 1) / taux_mort(n))))
  gamma / sum(gamma)
}
lam <- 0.1; mu <- 0.5
panne <- function(n) (3 - n) * lam
for (techniciens in 1:2) {
  pi_ <- stationnaire_nm(panne, function(n) min(n, techniciens) * mu, 3)
  en_panne <- sum(0:3 * pi_)
  cat(techniciens, round(pi_, 3), round(en_panne, 3), round(3 - en_panne, 2), "\n")
}
## 1 0.53 0.318 0.127 0.025 0.648 2.35
## 2 0.577 0.346 0.069 0.007 0.506 2.49

Exercice · Des clients découragés

difficile

Devant un guichet unique, les clients arrivent au taux $\lambda$, mais se découragent quand il y a du monde : lorsqu’il y a déjà $n$ personnes, un client qui se présente n’entre qu’avec probabilité $1/(n+1)$. Le service est exponentiel de taux $\mu$.

1. Justifiez que le nombre de personnes présentes est un processus de naissance et de mort, et donnez ses taux.
2. Montrez que la loi stationnaire existe toujours et identifiez-la.
3. Application : $\lambda = 3$ clients par heure, $\mu = 2$ par heure. Quelle est la probabilité que le guichet soit vide ? Combien de personnes sont présentes en moyenne ?

Voir la solution

1. Les arrivées effectives se font au taux $\lambda_n = \lambda/(n+1)$ (un processus d’arrivées au taux $\lambda$ dont on ne garde chaque arrivée qu’avec probabilité $1/(n+1)$), et les départs au taux $\mu_n = \mu$ pour $n \geq 1$. L’état ne varie que de $\pm 1$ : c’est un processus de naissance et de mort.
2. $\gamma_n = \dfrac{\lambda \cdot \tfrac{\lambda}{2} \cdots \tfrac{\lambda}{n}}{\mu^n} = \dfrac{u^n}{n!}$ avec $u = \lambda/\mu$. La série $\sum u^n/n!$ converge toujours (vers $e^u$), donc le régime stationnaire existe quels que soient $\lambda$ et $\mu$, et $\Pi_n = e^{-u}\,u^n/n!$ : c’est une loi de Poisson de paramètre $u$, la même que pour la file M/M/∞.
3. $u = 1{,}5$ : $\Pi_0 = e^{-1{,}5} \approx 0{,}223$, et le nombre moyen de personnes présentes est l’espérance de la loi de Poisson, $u = 1{,}5$.

En code

Python

import numpy as np
from scipy.stats import poisson

lam, mu, n_max = 3.0, 2.0, 60            # espace tronqué à 60 personnes
gamma = np.ones(n_max + 1)
for n in range(1, n_max + 1):
    gamma[n] = gamma[n - 1] * (lam / n) / mu   # lambda_{n-1} = lam / n, mu_n = mu
pi = gamma / gamma.sum()

print(pi[0], (np.arange(n_max + 1) * pi).sum())   # 0.2231 1.5
print(np.abs(pi - poisson.pmf(np.arange(n_max + 1), lam / mu)).max())  # 5.6e-17 : Poisson(1.5)

R

lam <- 3; mu <- 2; n_max <- 60          # espace tronqué à 60 personnes
gamma <- cumprod(c(1, (lam / (1:n_max)) / mu))   # lambda_{n-1} = lam / n, mu_n = mu
pi_ <- gamma / sum(gamma)

c(vide = pi_[1], moyenne = sum(0:n_max * pi_))   # 0.2231 1.5
max(abs(pi_ - dpois(0:n_max, lam / mu)))         # 5.6e-17 : c'est une Poisson(1.5)