Passez au temps continu

Un serveur de calcul ne tombe pas en panne « à la fin de l’heure » : il tombe en panne à 14 h 37 min 12 s, puis il est réparé à 16 h 05. Découper le temps en pas réguliers, comme dans la partie 1, est commode mais artificiel. Et si l’on veut modéliser un système où les événements surviennent à n’importe quel instant ?

On passe alors aux chaînes de Markov à temps continu. Presque tout ce que vous savez reste vrai, à une différence près : la matrice de transition $P$ est remplacée par une matrice de taux, le générateur. Nous suivrons un serveur qui tombe en panne et qu’on répare, et nous calculerons sa disponibilité.

La loi exponentielle, l’horloge sans mémoire

Commençons par la durée de fonctionnement du serveur avant la panne, notée $D$. En fiabilité, on la modélise très souvent par une loi exponentielle de paramètre $\lambda > 0$, notée $\mathcal{E}(\lambda)$ :

$$f(t) = \lambda\, e^{-\lambda t}, \qquad P(D > t) = e^{-\lambda t}, \qquad E[D] = \frac{1}{\lambda} \qquad (t \geq 0)$$

• $f$ est la densité de $D$ ;
• $P(D > t)$ est la probabilité que le serveur fonctionne encore à l’instant $t$ (on l’appelle la fonction de survie) ;
• $1/\lambda$ est la durée moyenne de fonctionnement, et $\lambda$ s’interprète comme un taux : un nombre moyen de pannes par unité de temps de fonctionnement.

Pour notre serveur, prenons $\lambda = 0{,}02$ panne par heure, soit une durée moyenne de fonctionnement de 50 heures. La probabilité qu’il tienne plus de 100 heures vaut $e^{-0{,}02 \times 100} = e^{-2} \approx 0{,}135$.

La survie du serveur décroît exponentiellement. À la durée moyenne de 50 h, il ne reste que 37 % de serveurs en marche : la majorité tombe en panne avant la moyenne. À 100 h, il en reste 13,5 %.

La propriété essentielle de cette loi est d’être sans mémoire :

$$P(D > t + u \mid D > u) = P(D > t)$$

Sachant que le serveur fonctionne depuis $u$ heures, la probabilité qu’il tienne encore $t$ heures est la même que s’il venait d’être mis en route. Après 80 heures sans panne, la probabilité de tenir 50 heures de plus vaut toujours $e^{-1} \approx 0{,}37$. L’exponentielle est la seule loi continue qui ait cette propriété, et c’est précisément ce qui rend le processus markovien : l’état présent (« le serveur marche ») suffit, il n’y a pas besoin de savoir depuis quand.

Sachant que le serveur a tenu 80 h, sa courbe de survie repart de 1 avec exactement la même forme : tenir 50 h de plus a la même probabilité, 0,37, que pour un serveur neuf.

Attention

L’absence de mémoire est une hypothèse forte. Un disque dur qui s’use a au contraire une mémoire : plus il est vieux, plus il risque de lâcher. Le modèle exponentiel convient aux pannes « accidentelles » (taux constant), pas à l’usure. Si l’usure compte, il faut enrichir l’état (par exemple ajouter un état « dégradé »), comme à l’exercice 2.

Deux propriétés de calcul serviront tout au long de la partie. Si $D_1 \sim \mathcal E(\lambda_1)$ et $D_2 \sim \mathcal E(\lambda_2)$ sont indépendantes, alors :

• le premier des deux événements survient après un temps $\min(D_1, D_2) \sim \mathcal E(\lambda_1 + \lambda_2)$ : les taux s’additionnent ;
• c’est le premier qui gagne la course avec probabilité $\lambda_1 / (\lambda_1 + \lambda_2)$.

Du temps discret au temps continu

Une chaîne de Markov à temps continu est un processus $(X_t)_{t \geq 0}$ à valeurs dans un espace discret $E$ qui vérifie la propriété de Markov à tous les instants : pour tous instants $t_0 < t_1 < \dots < t_n$,

$$P(X_{t_n} = i_n \mid X_{t_{n-1}} = i_{n-1}, \dots, X_{t_0} = i_0) = P(X_{t_n} = i_n \mid X_{t_{n-1}} = i_{n-1}).$$

On la suppose homogène : $p_{ij}(t) = P(X_{s+t} = j \mid X_s = i)$ ne dépend que de la durée $t$, pas de $s$. On range ces probabilités dans la matrice $G(t) = (p_{ij}(t))$, qui joue le rôle de $P^n$. Tout ce que vous connaissez se transpose :

Temps discret	Temps continu
$p(n) = p(0)\,P^n$	$p(t) = p(0)\,G(t)$
$P^{n+m} = P^n\,P^m$	$G(t+s) = G(t)\,G(s)$ (Chapman-Kolmogorov)
$P^0 = I$	$G(0) = I$

Mais il n’y a plus de « pas élémentaire » : comment décrire la chaîne par un seul objet, comme $P$ en temps discret ? En regardant ce qui se passe sur un intervalle de temps infiniment court.

Le générateur infinitésimal

Sur un petit intervalle de durée $h$, un serveur en marche tombe en panne avec une probabilité à peu près proportionnelle à $h$ : $p_{01}(h) \approx \lambda h$. C’est cette proportionnalité qu’on formalise. Pour $i \neq j$ :

$$p_{ij}(h) = \lambda_{ij}\, h + o(h), \qquad p_{ii}(h) = 1 + \lambda_{ii}\, h + o(h)$$

• $\lambda_{ij} \geq 0$ est le taux de transition de $i$ vers $j$ ;
• $\lambda_{ii} = -\sum_{j \neq i} \lambda_{ij} \leq 0$, et $-\lambda_{ii}$ est le taux de sortie de l’état $i$ ;
• $o(h)$ désigne un terme négligeable devant $h$ quand $h \to 0$.

La matrice $A = (\lambda_{ij})$ est le générateur infinitésimal de la chaîne. C’est la dérivée de $G$ en 0 :

$$A = G'(0) = \lim_{h \to 0^+} \frac{G(h) - I}{h}.$$

Comme chaque ligne de $G(h)$ somme à 1 et chaque ligne de $I$ aussi, chaque ligne de $A$ somme à 0. Pour notre serveur, avec l’état 0 = marche et 1 = panne, un taux de panne $\lambda$ et un taux de réparation $\mu$ :

$$A = \begin{pmatrix} -\lambda & \lambda \\ \mu & -\mu \end{pmatrix}.$$

Le graphe de transition porte maintenant les taux sur les arcs, et n’a pas de boucles (rester sur place n’est pas une transition).

Mais que fait concrètement la chaîne ? Grâce aux propriétés de l’exponentielle, on peut la décrire comme une suite de séjours et de sauts : dans l’état $i$, elle reste un temps de loi $\mathcal E(-\lambda_{ii})$, puis elle saute vers $j \neq i$ avec probabilité $\lambda_{ij} / (-\lambda_{ii})$. C’est exactement une course d’horloges exponentielles, une par transition possible, dont la première qui sonne décide du saut.

Le graphe porte des taux, pas des probabilités. Une trajectoire simulée alterne de longs séjours en marche (50 h en moyenne) et de brèves pannes (2 h en moyenne), chaque durée étant exponentielle.

À retenir

Une chaîne à temps continu est décrite par son générateur $A$ : hors diagonale, les taux de transition $\lambda_{ij} \geq 0$ ; sur la diagonale, l’opposé du taux de sortie. Les lignes de $A$ somment à 0 (celles de $P$ sommaient à 1). On reste dans $i$ un temps exponentiel de paramètre $-\lambda_{ii}$.

Les équations de Kolmogorov et le régime transitoire

Dérivons la relation de Chapman-Kolmogorov $G(t + s) = G(t)\,G(s)$ par rapport à $s$, puis par rapport à $t$, et prenons la dérivée en 0. On obtient les équations de Kolmogorov :

$$G'(t) = A\, G(t) \quad \text{(équation arrière)}, \qquad G'(t) = G(t)\, A \quad \text{(équation avant)}.$$

En multipliant l’équation avant à gauche par $p(0)$, on obtient l’équation d’évolution des probabilités d’état :

$$p'(t) = p(t)\, A.$$

Pour un espace d’états fini, ce système différentiel linéaire a pour solution $G(t) = e^{tA}$, l’exponentielle de matrice. C’est l’analogue continu de $P^n$.

Exemple

Pour le serveur, $p'(t) = p(t)\,A$ donne $p_0'(t) = -\lambda\,p_0(t) + \mu\,p_1(t)$. Comme $p_1 = 1 - p_0$, cela devient

$$p_0'(t) = -(\lambda + \mu)\,p_0(t) + \mu,$$

une équation linéaire du premier ordre. Si le serveur marche à $t = 0$ :

$$p_0(t) = \frac{\mu}{\lambda + \mu} + \frac{\lambda}{\lambda + \mu}\, e^{-(\lambda + \mu)t}.$$

Avec $\lambda = 0{,}02$ et $\mu = 0{,}5$ par heure (réparation en 2 heures en moyenne), la probabilité qu’il marche 2 heures plus tard vaut $0{,}9615 + 0{,}0385\,e^{-1{,}04} \approx 0{,}975$. Le terme exponentiel s’éteint à la vitesse $\lambda + \mu = 0{,}52$ par heure : au bout d’une dizaine d’heures, on a oublié l’état initial.

Quel que soit l'état de départ, p₀(t) rejoint la disponibilité 0,9615 en une dizaine d'heures, à la vitesse e^(−0,52 t). Le zoom montre la valeur transitoire 0,975 à t = 2 h, encore marquée par le départ en marche.

Pour la classification des états, rien ne change par rapport à la partie 1 : accessibilité, classes, irréductibilité et états absorbants ($\lambda_{ij} = 0$ pour tout $j \neq i$) se lisent sur le graphe des taux. Seule différence, la notion de période disparaît : en temps continu, on peut revenir dans un état après n’importe quelle durée.

Le régime stationnaire : ΠA = 0

À long terme, $p(t)$ se stabilise et sa dérivée s’annule. Comme $p'(t) = p(t)\,A$, la limite $\Pi$ doit vérifier $\Pi A = 0$.

À retenir

Si la chaîne à temps continu est irréductible et récurrente positive (automatique si elle est finie), alors $\Pi_j = \lim_{t \to \infty} p_{ij}(t) = \lim_{t \to \infty} p_j(t)$ existe, ne dépend pas de l’état initial, et $\Pi$ est l’unique solution de

$$\Pi\, A = 0, \qquad \sum_j \Pi_j = 1.$$

Pas besoin d’hypothèse d’apériodicité : elle n’a pas de sens en temps continu.

Chaque équation de $\Pi A = 0$ a une lecture physique, appelée équation de balance globale : en régime stationnaire, le flux de probabilité qui sort d’un état égale le flux qui y entre.

$$\Pi_j \sum_{k \neq j} \lambda_{jk} = \sum_{i \neq j} \Pi_i\, \lambda_{ij}$$

À gauche, la probabilité d’être en $j$ multipliée par le taux de sortie ; à droite, la somme des flux venant des autres états.

Exemple

Pour le serveur, la balance de l’état « marche » s’écrit $\lambda\,\Pi_0 = \mu\,\Pi_1$. Avec $\Pi_0 + \Pi_1 = 1$ :

$$\Pi_0 = \frac{\mu}{\lambda + \mu} = \frac{0{,}5}{0{,}52} \approx 0{,}9615, \qquad \Pi_1 = \frac{\lambda}{\lambda + \mu} \approx 0{,}0385.$$

Le serveur est disponible 96,15 % du temps. On retrouve la formule classique des ingénieurs fiabilistes : disponibilité = MTTF / (MTTF + MTTR), c’est-à-dire durée moyenne de fonctionnement sur durée moyenne d’un cycle, ici $50 / (50 + 2)$.

Piège

Ne confondez pas la limite $\lim_{t \to \infty} p(t)$, qui décrit le régime stationnaire (ou permanent), avec le régime transitoire, qui décrit $p(t)$ à un instant $t$ fini. La disponibilité 0,9615 est une valeur stationnaire ; la probabilité 0,975 calculée à $t = 2$ h est une valeur transitoire, qui dépend de l’état initial.

Le code suivant calcule $G(t) = e^{tA}$ avec SciPy, résout $\Pi A = 0$, et simule la chaîne par sa description « séjour exponentiel puis saut ».

import numpy as np
from scipy.linalg import expm

lam, mu = 0.02, 0.5                          # taux de panne et de réparation (par heure)
A = np.array([[-lam, lam],
              [mu,  -mu]])
assert np.allclose(A.sum(axis=1), 0)          # les lignes d'un générateur somment à 0

G2 = expm(2 * A)                              # G(2) = exp(2A)
print(G2[0, 0])                               # P(marche à t=2 | marche à t=0) ≈ 0.975

# Régime stationnaire : Pi A = 0 et somme = 1
M = A.T.copy()
M[-1, :] = 1
pi = np.linalg.solve(M, [0, 1])
print(pi)                                     # ≈ [0.9615 0.0385]

# Simulation : séjour exponentiel dans l'état courant, puis saut
rng = np.random.default_rng(0)
t, x, temps_marche, horizon = 0.0, 0, 0.0, 1_000_000
while t < horizon:
    taux_sortie = -A[x, x]
    sejour = rng.exponential(1 / taux_sortie)  # NumPy attend la moyenne 1/taux
    if x == 0:
        temps_marche += min(sejour, horizon - t)
    t += sejour
    x = 1 - x                                  # deux états : on saute vers l'autre
print(temps_marche / horizon)                  # proche de 0.9615

Piège

rng.exponential prend en argument la moyenne $1/\lambda$, pas le taux $\lambda$. Écrire rng.exponential(lam) simule des pannes 2 500 fois trop rapprochées. C’est l’erreur la plus fréquente en simulation de files d’attente.

En résumé

• La loi exponentielle $\mathcal E(\lambda)$ a pour moyenne $1/\lambda$ et vérifie $P(D > t) = e^{-\lambda t}$ ; elle est sans mémoire, ce qui rend le processus markovien. Le minimum d’exponentielles indépendantes est exponentiel, de taux la somme des taux.
• Une chaîne à temps continu homogène vérifie $p(t) = p(0)\,G(t)$ et $G(t+s) = G(t)\,G(s)$.
• Le générateur $A = G'(0)$ contient les taux $\lambda_{ij}$ ; ses lignes somment à 0. On séjourne en $i$ un temps $\mathcal E(-\lambda_{ii})$, puis on saute vers $j$ avec probabilité $\lambda_{ij}/(-\lambda_{ii})$.
• Équations de Kolmogorov : $G'(t) = A\,G(t) = G(t)\,A$, d’où $p'(t) = p(t)\,A$ et $G(t) = e^{tA}$.
• Pas de période en temps continu. Si la chaîne est irréductible (et récurrente positive), $\Pi$ est l’unique solution de $\Pi A = 0$, $\sum \Pi_j = 1$ : flux sortant = flux entrant.
• Disponibilité d’un équipement : $\mu / (\lambda + \mu)$ = MTTF / (MTTF + MTTR).

Exercices

Exercice · Une ampoule sans mémoire

facile

La durée de vie d’une LED suit une loi exponentielle de moyenne 1 000 heures.

1. Quelle est la probabilité qu’elle dure plus de 1 000 heures ?
2. Elle fonctionne depuis 800 heures. Quelle est la probabilité qu’elle tienne encore 500 heures ?

Voir la solution

Le taux vaut $\lambda = 1/1000$ par heure.

1. $P(D > 1000) = e^{-1} \approx 0{,}368$. Remarquez que ce n’est pas 1/2 : la loi exponentielle n’est pas symétrique autour de sa moyenne, et la majorité des LED lâchent avant la durée moyenne.
2. Par l’absence de mémoire, $P(D > 1300 \mid D > 800) = P(D > 500) = e^{-0{,}5} \approx 0{,}607$. Les 800 heures déjà écoulées ne changent rien.

En code

Python

import numpy as np
from scipy.stats import expon

duree = expon(scale=1000)                # SciPy attend la moyenne 1/lambda
print(duree.sf(1000))                    # 0.3679 = P(D > 1000)
print(duree.sf(1300) / duree.sf(800))    # 0.6065 = P(D > 1300 | D > 800)
print(duree.sf(500))                     # 0.6065 = P(D > 500) : sans mémoire

# Vérification par simulation
rng = np.random.default_rng(0)
d = rng.exponential(1000, 1_000_000)
print((d > 1000).mean(), (d[d > 800] > 1300).mean())   # 0.3677 0.6065

R

taux <- 1 / 1000                        # R attend le taux lambda
pexp(1000, taux, lower.tail = FALSE)    # 0.3679 = P(D > 1000)
pexp(1300, taux, lower.tail = FALSE) /
  pexp(800, taux, lower.tail = FALSE)   # 0.6065 = P(D > 1300 | D > 800)
pexp(500, taux, lower.tail = FALSE)     # 0.6065 = P(D > 500) : sans mémoire

# Vérification par simulation
set.seed(0)
d <- rexp(1e6, taux)
c(mean(d > 1000), mean(d[d > 800] > 1300))   # 0.3683 0.6071

Exercice · Un serveur qui se dégrade

moyen

On raffine le modèle du serveur avec trois états : marche (M), dégradé (D), panne (P). Un serveur en marche se dégrade au taux 0,1 par heure. Un serveur dégradé tombe en panne au taux 0,2 par heure, ou est remis en état par une maintenance préventive au taux 0,3 par heure. Un serveur en panne est réparé (et revient en marche) au taux 0,5 par heure.

1. Écrivez le générateur $A$ et vérifiez ses lignes.
2. Combien de temps un serveur reste-t-il dégradé en moyenne ? Avec quelle probabilité quitte-t-il cet état par une panne ?
3. Calculez la loi stationnaire. Quelle proportion du temps le serveur est-il en panne ?

Voir la solution

1. Dans l’ordre M, D, P :

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

Chaque ligne somme bien à 0.

2. Le taux de sortie de D vaut $0{,}3 + 0{,}2 = 0{,}5$, donc le séjour moyen est de $1/0{,}5 = 2$ heures. La sortie se fait par une panne avec probabilité $0{,}2/0{,}5 = 0{,}4$.

3. Balances globales : pour D, $0{,}5\,\Pi_D = 0{,}1\,\Pi_M$, donc $\Pi_D = 0{,}2\,\Pi_M$ ; pour P, $0{,}5\,\Pi_P = 0{,}2\,\Pi_D$, donc $\Pi_P = 0{,}08\,\Pi_M$. On vérifie la balance de M : $0{,}1\,\Pi_M = 0{,}3 \times 0{,}2\,\Pi_M + 0{,}5 \times 0{,}08\,\Pi_M = 0{,}1\,\Pi_M$. La normalisation donne $\Pi_M\,(1 + 0{,}2 + 0{,}08) = 1$, soit

$$\Pi = (0{,}78125 \,;\, 0{,}15625 \,;\, 0{,}0625).$$

Le serveur est en panne 6,25 % du temps.

En code

Python

import numpy as np

# États M (marche), D (dégradé), P (panne)
A = np.array([[-0.1, 0.1, 0.0],
              [0.3, -0.5, 0.2],
              [0.5, 0.0, -0.5]])
print(A.sum(axis=1))                     # [0. 0. 0.]

sortie_D = -A[1, 1]
print(1 / sortie_D, A[1, 2] / sortie_D)  # 2.0 heures en D, sortie par panne : 0.4

# Pi A = 0 et somme = 1
M = A.T.copy()
M[-1, :] = 1
pi = np.linalg.solve(M, [0, 0, 1])
print(pi)                                # [0.78125 0.15625 0.0625 ]

R

# États M (marche), D (dégradé), P (panne)
A <- matrix(c(-0.1,  0.1,  0,
               0.3, -0.5,  0.2,
               0.5,  0,   -0.5), nrow = 3, byrow = TRUE)
rowSums(A)                              # 0 0 0

sortie_D <- -A[2, 2]
c(1 / sortie_D, A[2, 3] / sortie_D)     # 2 heures en D, sortie par panne : 0.4

# Pi A = 0 et somme = 1
M <- t(A)
M[3, ] <- 1
solve(M, c(0, 0, 1))                    # 0.78125 0.15625 0.06250

Exercice · Combien de temps pour oublier une panne ?

difficile

On reprend le serveur à deux états ($\lambda = 0{,}02$, $\mu = 0{,}5$ par heure). À l’instant 0, il vient de tomber en panne.

1. Résolvez $p_0'(t) = -(\lambda + \mu)\,p_0(t) + \mu$ avec $p_0(0) = 0$.
2. Quelle est la probabilité qu’il fonctionne 2 heures plus tard ?
3. Au bout de combien de temps $p_0(t)$ est-elle à moins de 1 % de sa limite (en écart relatif) ?

Voir la solution

1. La solution générale est $p_0(t) = \dfrac{\mu}{\lambda + \mu} + C\,e^{-(\lambda + \mu)t}$. La condition $p_0(0) = 0$ impose $C = -\dfrac{\mu}{\lambda + \mu}$, donc

$$p_0(t) = \frac{\mu}{\lambda + \mu}\left(1 - e^{-(\lambda + \mu)t}\right).$$

2. $p_0(2) = 0{,}9615 \times (1 - e^{-1{,}04}) \approx 0{,}9615 \times 0{,}6465 \approx 0{,}622$. C’est bien moins que la disponibilité stationnaire : 2 heures après une panne, la réparation (de durée moyenne 2 heures) a une chance sur trois environ de ne pas être terminée.

3. L’écart relatif à la limite vaut $e^{-(\lambda + \mu)t}$. Il passe sous 1 % quand $0{,}52\,t > \ln 100 \approx 4{,}605$, soit $t > 8{,}86$ heures environ. Au-delà, l’état initial n’a pratiquement plus d’influence.

En code

Python

import numpy as np
from scipy.linalg import expm

lam, mu = 0.02, 0.5
A = np.array([[-lam, lam],
              [mu, -mu]])

p0 = np.array([0.0, 1.0])                # le serveur vient de tomber en panne
print((p0 @ expm(2 * A))[0])             # 0.6217 = P(marche à t = 2 h)
print(mu / (lam + mu) * (1 - np.exp(-(lam + mu) * 2)))   # 0.6217, formule fermée

t_oubli = np.log(100) / (lam + mu)       # écart relatif e^{-(lam+mu)t} < 1 %
print(t_oubli)                           # 8.86 heures

R

lam <- 0.02; mu <- 0.5
A <- matrix(c(-lam, lam,
              mu, -mu), nrow = 2, byrow = TRUE)

# exp(tA) par diagonalisation (A est diagonalisable)
expm_diag <- function(M) {
  e <- eigen(M)
  Re(e$vectors %*% diag(exp(e$values)) %*% solve(e$vectors))
}
p0 <- c(0, 1)                           # le serveur vient de tomber en panne
(p0 %*% expm_diag(2 * A))[1]            # 0.6217 = P(marche à t = 2 h)
mu / (lam + mu) * (1 - exp(-(lam + mu) * 2))   # 0.6217, formule fermée

log(100) / (lam + mu)                   # 8.86 heures pour un écart relatif < 1 %