Avant-propos - La fiabilité, c’est…

\[R(t)=Pr(T>t)\] (où T est la variable aléatoire qui associe le temps de bon fonctionnement.)

L’analyse de survie est particulièrement utile pour diagnostiquer les défaillances, imprévues ou prématurées. De plus cette analyse permet de déterminer les périodes de garantie afin d’éviter les coûts de réparations ou de remplacements excessifs.

Waloddi Weibull a inventé la distribution qui porte son nom en 1937 et a livré son article emblématique sur ce sujet il y a 70 ans en septembre 51 devant l’ASME (American Society of Mechanical Engineers). Il illustra son exposé avec sept exemples allant de la résistance de l’acier à la taille des hommes adultes dans les îles Britanniques… [Waloddi Weibull] “A Statistical Distribution Function of Wide Applicability”, ASME Joumal Of Applied Mechanics Paper, 1951.

Waloddi Weibull (1887-1979)

L’avantage de sa distribution est qu’elle s’applique à un large éventail de problèmes y compris ceux avec des taux de défaillances croissants, constants ou décroissants. Les anomalies sont mises en évidence par des graphiques qui aident les ingénieurs à diagnostiquer les causes premières des défaillances.

L’étude des défaillances des systèmes ou des durées de vie a pour objectif de répondre à des questions du type :

Bien que dans les années 50 sa méthode fut contestée (Weibull lui même était assez scéptique… ) aujourd’hui, l’analyse de Weibull est très utilisée pour l’ajustement et l’analyse des données de durée de vie. L’U.S. Air Force a reconnu le mérite de la méthode de Weibull et a financé ses recherches.

Le principal avantage de l’analyse de Weibull est sa capacité à fournir un modèle raisonnablement précis avec de faibles échantillons.

Weibulll a mis l’accent en particulier sur la polyvalence de sa distribution, en effet elle dépend de trois paramètres, \(\beta\), \(\eta\) et \(\gamma\)

\[R(t)=e^{-(\frac{t-\gamma}{\eta})^{\beta}}\]

Ainsi de nombreuses formes de distributions sont décrites par Weibull (exactement ou par approximation hormis le log-normal qui est aussi très utilisé en analyse de survie) comme la loi Normale, exponentielle, Rayleigh et même Poisson

La MTTF (Temps moyen de fonctionnement avant panne) est égale à \(\eta\) si \(\beta=1\)

Dans le cas général ils sont reliés par : \(MTTF=\eta\cdot\Gamma(1+\frac{1}{\beta})\)

La MTBF, le temps moyen entre pannes (Mean Time Between Failures), ne doit pas être confondu avec le MTTF (Mean time to failure), temps moyen avant défaillance. La MTBF est estimée en divisant le temps de fonctionnement total de toutes les unités par le nombre de pannes observées. Ce sont des paramètres différents, mais égaux lorsqu’il n’y a pas de suspensions (en français on parle de données censurées). Sinon, ils peuvent être très différents.

La MTBF est utilisée avec des systèmes réparables.

(Image de Seongwoo Woo https://www.researchgate.net/publication/334205633)

Les données

Il faut réfléchir sur les données, car ce sont elles qui détermineront les paramètres. Comme toute analyse, la précision dépend de la quantité et de la qualité des données analysées. Dans une analyse de Weibull il y a deux catégories de données : les temps de défaillance (TTF) et les données censurées ou suspendues.

Les données sont des intervalles de l’unité considérée (celle caractérisant la durée de vie du produit).

L’omission de données censurées produit des résultats qui sous-estiment la fiabilité réelle de l’élément analysé.

Un exemple (tiré du [Abernethy])

Nous inspectons des rivets d’étanchéité d’un compresseur. Nous avons constaté que sur certains rivets la partie évasée du rivet était manquante.

Un test accéléré en laboratoire utilisant les anciens rivets a permis d’établir une méthode. Les résultats sont présentés dans le tableau ci-dessous. Les numéros de série des rivets 3, 6 et 7 sont considérés comme des défaillances non représentatives (“Suspensions”) car ils ont été produits par deux autres modes de défaillance et non pas le mode de défaillance spécifique étudié.

Numéro Failure time (mn) Remarque
1 90 F (failure)
2 96 F
3 100 S
4 30 F
5 49 F
6 45 S
7 10 S
8 82 F

En première analyse ces trois défaillances seront ignorées (elles sont censurées, “suspensions” à droite) mais seront correctement prises en compte par la suite. Il reste donc cinq points d’échec pour la défaillance.

Pour estimer les fréquences des défaillances l’approximation de [Benard] dite méthode des rangs médians (Benard’s Median Rank) est la plus populaire (Voir [Johnson 1964]).

\[ Benard's Median Rank = \frac{i-0.3}{N+0.4} \]

Nous obtenons :

i Failure time t Median Rank %
1 30 12.96
2 49 31.48
3 82 50
4 90 68.52
5 96 87.04

Que nous pourrions retrouver avec R par :

ppoints(1:5,a=0.3)
## [1] 0.1296296 0.3148148 0.5000000 0.6851852 0.8703704

Précédemment nous avions ignoré certaines données (celles “censurées”). Or ces données ont une influence même si elles ne sont pas tracées. L’ordre de classement et les rangs médians doivent être ajustés en tenant compte de ces données censurées.

Leonard Johnson a donné une méthode qui tient compte des données censurées, puis elle fut améliorée et simplifiée par Drew Auth (voir page 52).

Formule d’ajustement des rangs :

On note \(RA(i)\) : Rang Ajusté au rang \(i\)

\(RA(i)=\frac{s(t)RA(i-1)+N+1}{s(t)+1}\) (où \(s\) donne le nombre de “survivants” à \(t\))

i t S(t) RA Rang médian
1 10S 8 Suspended
2 30F 7 (7x0+8+1)/(7+1)=1.125 9.8%
3 45S 6 S
4 49F 5 2.438 25.5%
5 82F 4 3.750 41.1%
6 90F 3 5.063 56.7%
7 96F 2 6.375 72.3%
8 100S 1 S

Là encore nous retrouvons facilement ces résultats avec le package WeibullR

library(WeibullR)
failures<-c(30,49,82,90,96)
suspensions<-c(10,45,100)
(getPPP(failures,suspensions))
##   time        ppp adj_rank
## 1   30 0.09756403   1.1250
## 2   49 0.25330810   2.4375
## 3   82 0.41023658   3.7500
## 4   90 0.56732480   5.0625
## 5   96 0.72429972   6.3750

L’effet majeur des suspensions est d’augmenter \(\eta\) ; \(\beta\) sera peu affecté. Par conséquent, ignorer les suspensions rendra le modèle plus pessimiste.

Mise en pratique avec R

Reprenons l’exemple de nos rivets.

failures<-c(30, 49, 82, 90, 96)
res<-MRRw2p(failures,bounds = T,show = T)

Avec les suspensions :

failures<-c(30,49,82,90,96)
suspensions<-c(10,45,100)
MRRw2p(x = failures,s = suspensions,bounds = T,show = T)->res

Comparons les deux cas :

obj1<-wblr(x = failures)
obj1<-wblr.conf(wblr.fit(obj1,col="red",dist="weibull2p"),lwd=1)
obj2<-wblr(x = failures,s = suspensions)
obj2<-wblr.conf(wblr.fit(obj2,col="blue",dist="weibull2p"),lwd=1)
plot.wblr(list(obj1,obj2))

La droite bleue a un R² supérieur à la rouge et nous avons pris en compte les suspensions ce qui donne un \(\eta\) d’environ 95.

Interprétation de \(\beta\)

Le taux d’avarie moyen par unité de temps est bien :

\[\frac{R(t)-R(t+\Delta t)}{\Delta tR(t)}\]

Donc le taux instantané est (\(\Delta t \rightarrow 0\)) :

\[\lambda(t)=\frac{-R'(t)}{R(t)}\]

L’étude de \(\lambda\) nous donne la célèbre courbe en baignoire, qui traduit la relation entre \(\beta\) est le taux de défaillance instantané :

Public Domain : https://commons.wikimedia.org/w/index.php?curid=7458336


Exemple : Courbe en S

Les données d’un diagramme de Weibull peuvent apparaître courbées ou mal ajustées, comme l’illustre la figure ci-dessous. Cela peut-être dû à plusieurs modes de défaillance ou un changement soudain de la cause de défaillance ou bien alors que l’environnement a évolué…

tf<-c(3.46,3.73,4.05,4.81,5.89,10,20,36.87, 53.91,80)
obj<-wblr(x = tf)
obj<-wblr.conf(wblr.fit(obj,col="blue",dist="weibull2p"),lwd=1,col="red")
plot.wblr(obj,main="Forme en S")

Cette courbe en S peut résulter du fait que la moitié de la population échoue prématurément (défauts de jeunesse) tandis que le reste échouera en raison d’une forte usure.

Une autre raison est peut être due à l’origine des temps qui se trouve au mauvais endroit. Le temps peut ne pas commencer à zéro.

Par exemple, il peut être physiquement impossible que le mode de défaillance se produise au début de la vie. Une défaillance de roulement à billes due à l’écaillage ou à un déséquilibre ne peut se produire sans que la rotation du roulement n’entraîne suffisamment de dommages pour que le roulement ne soit plus en mesure de fonctionner.

Le temps commence lorsque des défaillances sont possibles.

Le modèle de Weibull à 3 paramètres, détermine le paramètre \(\gamma\) comme origine des temps.

tf<-c(3.46,3.73,4.05, 4.81,5.89,10,20,36.87, 53.91,80)
MRRw3p(x = tf,show = T)

##        Eta       Beta         t0       Rsqr 
## 11.2514136  0.4353637  3.4373115  0.9883045

Avec un \(\gamma=3.44\) nous constatons une nette amélioration et un coefficient de corrélation R²=0.988 nettement meilleur.

Le \(\beta=0.435\) indique des défauts de jeunesse

Une courbe concave (\(\gamma >0\)) est plus fréquente qu’une courbe convexe \(\gamma<0\). Un \(\gamma\) négatif est assez difficile à interpréter physiquement (vieillissement prématuré, mixage entre plusieurs modes de défaillance…)

Prévoir les défaillances - Analyse des risques

Une analyse des risques a été développée en 83 chez Pratt & Whitney Aircraft [Abernethy 1983]. Aujourd’hui nous l’appelons “Analyse prédictive”…

Une analyse des risques prédit le nombre attendu d’incidents pouvant survenir sur une période de temps spécifique.

Quelques exemples :

I] Etude adaptée et tirée du [Abernethy]

Soit un parc de 20 pompes, l’une a échoué à 2000 heures et une autre à 3000 heures. Il y a 5 suspensions à 1000 heures et 4000 heures puis 4 à 2000 et 3000 heures.

Générons le tracé de Weibull

library(WeibullR)
failures<-c(2000,3000)
suspensions<-c(rep(c(1000,4000),each=5),rep(c(2000,3000),each=4))
MRRw2p(x = failures,s = suspensions,show = T)->res

Nous obtenons : \(\beta= 2,829\) et \(\eta = 5047\) heures.

La valeur de F(t) peut être lue directement à partir la fonction de densité cumulée de Weibull ou calculée à partir de l’équation de Weibull :

\[R(t)=e^{-(\frac{t}{\eta})^{\beta}}\]

Par exemple, au temps 1000 et 2000 :

F<-function(t){
1-exp(-(t/5047)^2.829)
}
c(F(1000),F(2000))
## [1] 0.01020688 0.07030700

La dernière colonne du tableau ci-dessous nous donne la probabilité qu’une pompe en état de marche à l’instant \(t\) ait une défaillance entre les instants \(t\) et \(t+u\)

t t+u F(t) F(t+u) \[ \frac{F(t+u)-F(t)}{1-F(t)} \]
1000 2000 0.010 0.070 0.061
2000 3000 0.070 0.205 0.145
3000 4000 0.205 0.404 0.251
4000 5000 0.404 0.622 0.366

Entre 3000 et 4000 heures le taux d’avarie et de 25.1%, le taux moyen par heure sera de 0.0251%

II] Nous étudions les cassures de cage de roulement à billes.

Sur 6 roulements. On nous donne le nombre de cycles avant rupture :

\(4\cdot 10^{5}, 1.3\cdot 10^{5},9.8\cdot 10^{5},2.7\cdot 10^{5},6.6\cdot 10^{5},5.2\cdot 10^{5}\)

Lançons une analyse :

cycles<-c(4,1.3,9.8,2.7,6.6,5.2)
MRRw2p(x = cycles,show = T)->res

Les fabricants de roulements nomment \(L_{10}\) la durée de vie nominale qui correspond à un seuil de fiabilité de 90%, c’est-à-dire 90% atteignent \(t=L_{10}\)

Dans notre cas avec \(R(t)=e^{-(\frac{t}{5.729})^{1.494}}\) on trouve \(L_{10}=1.27\cdot10^{5}\) cycles




Références