Méthode FORM-SORM

Un article de Wikipédia, l'encyclopédie libre.
Aller à : navigation, rechercher

La méthode FORM-SORM est une méthode utilisée en fiabilité pour déterminer la probabilité de défaillance d'un système.

  • FORM est l'acronyme de First Order Reliability Method, méthode de fiabilité du premier ordre ;
  • SORM est l'acronyme de Second Order Reliability Method, méthode de fiabilité du second ordre.

Il s'agit de trouver le point de conception, c'est-à-dire le cas de défaillance le plus probable, par des méthodes d'optimisation.

Présentation du problème[modifier | modifier le code]

Une conception « au plus juste » — compromis entre le coût de fabrication, le coût de service (maintenance, garantie) et la satisfaction de l'utilisateur (sécurité, disponibilité du système) — nécessite d'avoir recours aux statistiques ; c'est le domaine de la fiabilité. On veut connaître la probabilité que la durée de vie T d'un système soit supérieure à une valeur t, R(t) = P(T ≥ t), ce qui revient à connaître la probabilité de défaillance avant t,

F(t) = P(T ≤ t) ; F = 1 - R.

En particulier, on veut garantir que cette probabilité de défaillance est inférieure à une limite α, appelée « risque alpha », souvent prise à 5 ou 10 % :

F(t) ≤ α, c'est-à-dire R(t) ≥ 1 - α.

On se place en général à la limite

F(t) = α.

Cette durée de vie dépend de plusieurs facteurs :

  • des facteurs de charge, ou contraintes S (stress) : c'est ce que subit la machine, sollicitations mécaniques, humidité, température, intensité du courant électrique…
  • des facteurs de fabrication, ou résistance R : dimensions des composants, propriétés de la matière, …

Tous ces facteurs, contraintes ou résistances, peuvent s'exprimer par des grandeurs chiffrées xi. Ces valeurs sont des réalisations de variables aléatoires Xi : tous les systèmes, même produits à la chaîne, sont différents les uns des autres, et ils ne subissent pas non plus les mêmes sollicitations.

Dans un certain nombre de cas, on peut définir une fonction de défaillance g(x1, x2, …, xn) pour exprimer la condition sur la durée de vie :

F(t) ≤ α ⇔ g(x1, x2, …, xn) ≥ 0.

L'espace à n dimensions (x1, x2, …, xn) est donc coupé en deux demi-espaces :

  • le domaine de défaillance F (failure), vérifiant l'équation g(x1, x2, …, xn) < 0 ;
  • le domaine de bon fonctionnement, vérifiant l'équation g(x1, x2, …, xn) > 0 ;

la frontière entre les deux est une hypersurface d'équation g(x1, x2, …, xn) = 0 qui définit l'état limite.

Chaque variable aléatoire Xi suit une loi de probabilité dont la fonction de densité de probabilité est notée ƒi. La densité de probabilité d'avoir un n-uplet donné (x1, x2, …, xn) vaut donc

y = ƒ1(x1)׃2(x2)ׅ׃n(xn)

La fonction y(x1, x2, …, xn) forme une hypersurface de l'espace (x1, …, xn, ƒ). Le maximum de y, le « sommet » de l'hypersurface, est le point le plus probable.

Pour déterminer la probabilité de défaillance Pf du système, il faut intégrer y sur la zone de défaillance F

\mathrm{P_f} = \int_\mathrm{F} f_1(x_1) f_2(x_2) \cdots f_n(x_n) \mathrm{d}x_1 \mathrm{d}x_1 \cdots \mathrm{d}x_n

La valeur de Pf est le « volume » — l'hypervolume — délimité par les hypersurfaces y et g = 0.

La solution analytique de cette intégrale est en général complexe, voire impossible. Une manière d'estimer Pf consiste à utiliser une méthode de Monte-Carlo : on génère un grand nombre de valeurs aléatoires (x1, x2, …, xn) suivant les lois statistiques connues, et l'on compte le nombre de cas pour lesquels g est négatif.

Toutefois, la majeure partie de l'information est contenue dans une petite zone de l'espace, dite « zone critique », qui est l'ensemble des n-uplet (x1, x2, …, xn) pour lesquels la densité de probabilité de défaillance est importante ; le point où la densité de probabilité de défaillance est maximale est appelée « point de conception » et noté x* (c'est un n-uplet). Or, si la conception est robuste, cette zone critique est loin du sommet de l'hypersurface y, puisque l'on veut que la plupart des situations soient des situations de bon fonctionnement. Donc, la méthode de Monte-Carlo génère beaucoup de n-uplet proche du maximum de y, c'est-à-dire loin de la zone d'intérêt.

Il faut donc un nombre considérable de simulations pour avoir une bonne estimation de Pf.

La méthode FORM-SORM est une méthode d'approximation consistant à trouver le point de conception, et de s'en servir pour déterminer la probabilité de défaillance Pf.

Méthode générale[modifier | modifier le code]

La méthode FORM-SORM consiste à :

  1. Remplacer les variables aléatoires par des lois normales centrées réduites.
  2. Trouver le n-uplet (x1, x2, …, xn) de l'hypersurface d'état limite le plus proche de l'origine ; c'est alors le point de conception x*, le n-uplet de la zone de défaillance ayant la densité de probabilité maximale.
  3. Faire une approximation linéaire (FORM) ou quadratique (SORM) de g pour déterminer la probabilité de défaillance Pf.

Hypothèses[modifier | modifier le code]

On suppose que les variables aléatoires ont une distribution normale.

Si les variables sont normales d'espérance μi et d'écart type σi, on les transforme également en variables de loi centrées réduites

Ui = (Xi - μi)/σi
Xi = σiUi + μi

Si les variables aléatoires Xi ne sont pas normales, on les transforme pour avoir des variables Ui suivant lois normales centrées réduites ; cette opération est la « transformation de Nattaf et Rosenblatt » :

Ui = Φ-1(Fi(Xi))

et réciproquement

Xi = F-1i(Ui))

où Φ (lettre grecque Phi) est la fonction de répartition de la loi normale centrée réduite, et Fi la fonction de répartition de Xi.

Méthode générale[modifier | modifier le code]

Plus on est proche du sommet de l'hypersurface y, plus la densité de probabilité est élevée. Les fonctions de densité ƒi étant des courbes en cloche symétriques, y forme une « colline » : plus on s'éloigne du sommet, plus la densité de probabilité y diminue. Donc, le point le plus critique, le point de conception x*, est le point de la zone de défaillance F le plus proche du sommet.

La méthode FORM-SORM consiste à trouver le point de la zone frontière g(x1, …, xn) = 0 le plus proche du sommet. Dans l'espace centré réduit, le sommet de y se trouve à l'origine, le point de coordonnées (0, 0, …, 0). La fonction de défaillance devient alors g* :

F(t) = α ⇔ g*(u1, …, un) ≥ 0.

La distance à l'origine est donc

d^* = \sqrt{u_1^2 + u_2^2 + \cdots + u_n^2}

et ainsi,

la méthode consiste à minimiser d*(u1, …, un) avec pour contrainte g*(u1, …, un) = 0.

C'est une méthode d'optimisation quadratique : en effet, minimiser d* revient à minimiser d*2 (puisque la fonction carré est strictement croissante sur [0 ; +∞[), on doit donc minimiser une fonction quadratique ∑xi2 avec une contrainte linéaire.

La distance minimale trouvée est appelée indice de fiabilité et est noté β (lettre grecque bêta) :

β = d*(x*).

Méthode FORM[modifier | modifier le code]

On fait l'hypothèse que la fonction g est linéaire. Si ce n'est pas le cas, la méthode revient au final à en prendre le développement limité du premier ordre.

La probabilité de défaillance vaut alors

Pf = Φ(-β).

Méthode SORM[modifier | modifier le code]

La méthode SORM est identique à la méthode FORM, mais si g* n'est pas linéaire, alors on utilise une approximation du second ordre pour déterminer Pf.

On définit les courbures principales κi (lettre grecque kappa) :

\kappa_i = \frac{\partial^2 g^*}{\partial u_i^2}

Elles sont au nombre de n - 1, puisque l'équation g* = 0 définit une hypersurface de dimension n - 1.

La probabilité de défaillance s'écrit alors :

Pf = S1 + S2 + S3

avec :

\mathrm{S}_1 = \Phi(-\beta) \prod_{j = 1}^{n - 1} \frac{1}{\sqrt{1 + \beta \kappa_j}}
\mathrm{S}_2 = (\beta\Phi(-\beta) - \varphi(\beta)) \left (
 \prod_{j = 1}^{n - 1} \frac{1}{\sqrt{1 + \beta \kappa_j}}
 - \prod_{j = 1}^{n - 1} \frac{1}{\sqrt{1 - (\beta + 1)\kappa_j}} \right )
\mathrm{S}_3 = (\beta + 1)(\beta\Phi(-\beta) - \varphi(\beta)) \left [
 \prod_{j = 1}^{n - 1} \frac{1}{\sqrt{1 + \beta \kappa_j}}
 - \mathrm{Re} \left (\prod_{j = 1}^{n - 1} \frac{1}{\sqrt{1 - (\beta + i)\kappa_j}} \right ) \right ]

où Re est la partie réelle, i est l'imaginaire, et φ est la densité de probabilité de la loi normale centrée réduite.

La fonction S1 est l'approximation asymptotique de Pf lorsque β tend vers +∞ — c'est donc le terme prépondérant si la conception est robuste (le point de conception est loin de l'origine) —, et S2 et S3 sont des termes correctifs.

Application à un problème à deux variables normales[modifier | modifier le code]

Présentation du problème[modifier | modifier le code]

Distribution de la contrainte S et de la résistance R : toutes les situations (R, S) ont une densité de probabilité (surface en niveaux de gris). La zone où la marge m = R - S est positive est la zone où le système survit (R > S).
Même situation vue de dessus.
Point de conception : maximum de la courbe m = 0.

Dans la méthode contrainte-résistance, la condition sur la durée de vie

F(t) ≤ α

s'exprime par une résistance R supérieure à une contrainte S (stress) :

R ≥ S

ou plutôt en définissant la marge m = R - S :

m ≥ 0.

Cette marge m est donc une fonction de défaillance :

g(R, S) = m.

Si l'on dispose des distributions (densités de probabilité) ƒR et ƒS, on peut alors déterminer cette probabilité :

F(t) = P(R ≤ S) = P(m ≤ 0)

La probabilité d'avoir une couple de valeurs dans un intervalle ([S, S + dS], [R, R+dR]) est donné par

P([S, S + dS]∩[R, R+dR]) = P([S, S + dS])×P([R, R+dR]) = ƒS׃R×dS×dR.

On peut tracer la surface y = ƒS׃R, ainsi que la ligne de séparation m = 0 : cette ligne marque l'état limite entre la situation de survie et la situation de défaillance. Une méthode simple d'estimer la probabilité P(m ≤ 0) consiste à utiliser une méthode de Monte-Carlo : on effectue N tirages aléatoires de couples (R, S), et l'on dénombre les n situations pour lesquelles m ≤ 0 :

P(m ≤ 0) ≃ n/N.

Si l'on veut un risque α faible, cela signifie que le sommet de la surface doit être loin de l'état limite m = 0.

Le point que l'on veut le mieux connaître est le maximum de la courbe m = 0 ; c'est notre « point de conception ». Si le risque α est faible, alors le point de conception est loin du sommet de la surface. Cela signifie que lors du tirage aléatoire, peu de points seront proches du point de conception ; il faudra donc de très nombreux tirages — N devra être très grand — si l'on veut caractériser ce point.

Exemple

Nous avons représenté ci-contre une situation pour laquelle :

  • la contrainte S suit une loi normale d'espérance μS = 230 MPa et un écart type σS = 14 Mpa ;
  • la résistance R suit une loi normale d'espérance μR = 250 MPa et un écart type σR = 17 Mpa.

Graphiquement, on constate que le maximum de l'état limite se trouve à 238 MPa. Ce point est en général proche de l'intersection des courbes ƒS et ƒR.

Une méthode de Monte-Carlo avec un million de tirages et un intervalle de confiance de 5 %[1] donne 18,1 % ± 0,2 % de défaillance.

Mise en œuvre de la méthode FORM[modifier | modifier le code]

Détermination du point de conception par la méthode FORM-SORM.

La méthode FORM est une méthode permettant d'obtenir une valeur approchée de la probabilité de défaillance à partir du point de conception. Elle suppose que les distributions de S et de R sont normales.

La première étape consiste à passer dans l'espace des lois centrées réduites : on transforme les coordonnées pour avoir des probabilités de densité d'espérance nulle et de variance 1. Cette étape est appelée « transformation de Nattaf et Rosenblatt ». Les nouvelles coordonnées sont appelées S* et R* :

\mathrm{R}^* = \frac{\mathrm{R} - \mu_{\mathrm{R}}}{\sigma_{\mathrm{R}}} ;
\mathrm{S}^* = \frac{\mathrm{S} - \mu_{\mathrm{S}}}{\sigma_{\mathrm{S}}}.

Le sommet de la surface ƒS*׃R* se trouve donc aux coordonnées (0 ; 0). La droite d'état limite m = 0 (R - S = 0) devient donc :

\mathrm{R} - \mathrm{S} = 0 \Longleftrightarrow
\sigma_{\mathrm{R}} \times \mathrm{R}^* + \mu_{\mathrm{R}} - (\sigma_{\mathrm{S}} \times \mathrm{S}^* + \mu_{\mathrm{S}}) = 0

ce que l'on peut noter

m* = 0

avec

m*(R*, S*) = σR×R* + μR - (σS×S* + μS).

Le point de conception est le point de densité de probabilité maximale sur la droite m* = 0. Du fait des propriétés de symétrie de la surface ƒS*׃R*, c'est aussi le plus proche de l'origine (0 ; 0) — sommet de la surface —, c'est-à-dire celui pour lequel la quantité

d*2 = R*2 + S*2

est minimale.

Pour trouver ce point, on part d'un couple (R*, S*) arbitraire, et l'on applique un algorithme d'optimisation (itératif) pour avoir le minimum de d* avec comme contrainte m* = 0.

La distance d* atteinte est alors l'indice de fiabilité β, et la probabilité de défaillance Pf s'obtient avec la fonction de répartition Φ de la loi normale centrée réduite :

Pf = Φ(-β) ;

c'est la méthode FORM.

Notons que la surface obtenue avec la transformation de Nataff et Rosenblatt est toujours la même ; c'est la droite m* = 0 qui change d'un problème à l'autre.

Exemple

La transformation de Nattaf-Rosenblatt est :

\mathrm{R}^* = \frac{\mathrm{R} - 250}{17} ;
\mathrm{S}^* = \frac{\mathrm{S} - 230}{14}.

On a

m* = 17×R* + 250 - (14×S* + 230) = 17×R* - 14×S* + 20.

L'optimisation nous donne le point de conception suivant :

R* = -0,701 ; S* = 0,577 ; R = S = 238 MPa
β = 0,908 ;
Pf = 0,182 = 18,2 %.

Comme dans notre exemple les lois sont effectivement normales et g* est effectivement linéaire, nous avons un résultat « exact » aux arrondis près : l'approximation provient uniquement du fait que l'on a une résolution numérique.

Les logiciels ont en général deux types de solveurs. Certains utilisent un solveur non-linéaire « généraliste » ; il faut donc leur définir la fonction g, ainsi que la contrainte

h = 0, avec h = σR×R* - σS×S* + (μR - μS).

Le solveur cherche à minimiser une fonction qui est la distance à laquelle on ajoute une fonction permettant de prendre en compte la contrainte (par exemple une pénalité qui augmente vers +∞ lorsque l'on s'éloigne des conditions de contrainte). C'est le cas des solveurs des tableurs, et du solveur constrOptim() du langage R.

Mais nous sommes ici en présence d'un problème d'optimisation quadratique, ce qui permet d'utiliser un solveur spécifique (Goldfarb et Idnani) ; c'est le cas de la fonction solve.QP du langage R, et de la fonction qpsolve de Scilab. Le problème s'écrit alors de manière matricielle

g(\mathbf{x}) = \frac{1}{2}d^{*\,2} = \frac{1}{2} {}^\mathrm{t}\mathbf{x}
\times \mathrm{Q} \times \mathbf{x}
+ {}^\mathrm{t} p \times \mathbf{x}

avec

\mathbf{x} = \begin{pmatrix} \mathrm{R}^* \\ \mathrm{S}^* \end{pmatrix}\text{, }
\mathrm{Q} = \begin{pmatrix} 1 & 0 \\ 0 & 1\end{pmatrix}\text{, }
p = \begin{pmatrix} 0 \\ 0 \end{pmatrix}

avec une contrainte d'égalité unique :

x = b

avec

C = (σR ; -σS) ; b = (μS - μR).

Notes et références[modifier | modifier le code]

  1. l'erreur pour un risque α avec N tirages vaut
    \varepsilon = \sqrt{\frac{1}{2(\mathrm{N} + 1)} \cdot \ln \left ( \frac{2}{\alpha}\right )}
    soit
    \varepsilon = \sqrt{\frac{1}{2(10^6 + 1)} \cdot \ln \left ( \frac{2}{0,005}\right )} \simeq 1,7 \cdot 10^{-3} = 0,17\ \%
  2. Fonctions mathématiques (aide lLibreOffice)
  3. Fonctions statistiques - Quatrième partie (aide LibreOffice 3.3)
  4. Quel est le nombre maximum de cellules dans une feuille de calcul LibreOffice  ?, FAQ LibreOffice Calc

Voir aussi[modifier | modifier le code]

Bibliographie[modifier | modifier le code]

  • (en) Shankar Sankararaman et Kai Goebel, « Uncertainty Quantification in Remaining Useful Life of Aerospace Components using State Space Models and Inverse FORM », IEEE Aerospace Conference, IEEE,‎ mai 2013 (lire en ligne)
  • (en) L. Cizelj, B. Mavko et H. Riesch-Oppermann, « Application of first and second order reliability methods in the safety assessment of cracked steam generator tubing », Nuclear Engineering and Design, Elsevier, no 147,‎ 1994 (lire en ligne)