Voir la source de Méthode FORM-SORM

La '''méthode FORM-SORM''' est une méthode utilisée en [[fiabilité]] pour déterminer la [[loi de fiabilité|probabilité de défaillance]] d'un composant.
* 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 du maximum de vraisemblance. La solution est obtenue par une procédure d'optimisation.

== Présentation du problème ==

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 [[statistique]]s ; 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 α, choisie en fonction des incertitudes et des conséquences d'une situation de défaillance :
: 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 : La notion de durée de vie implique celle de processus, c'est-à-dire de variable aléatoire indexée sur le temps, ce que la méthode FORM présentée ne prend pas en compte.
* 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 [[Grandeur physique|grandeurs chiffrées]] ''x<sub>i</sub>''. Ces valeurs sont des réalisations de [[variable aléatoire réelle|variables aléatoires]] X<sub>''i''</sub> : 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''(''x''<sub>1</sub>, ''x''<sub>2</sub>, …, ''x<sub>n</sub>'') pour exprimer la condition sur la durée de vie :
: F(''t'') ≤ α ⇔ ''g''(''x''<sub>1</sub>, ''x''<sub>2</sub>, …, ''x<sub>n</sub>'') ≥ 0.
L'espace à ''n'' dimensions (''x''<sub>1</sub>, ''x''<sub>2</sub>, …, ''x<sub>n</sub>'') est donc coupé en deux demi-espaces :
* le domaine de défaillance F ''(failure)'', vérifiant l'équation ''g''(''x''<sub>1</sub>, ''x''<sub>2</sub>, …, ''x<sub>n</sub>'') &lt; 0 ;
* le domaine de bon fonctionnement, vérifiant l'équation ''g''(''x''<sub>1</sub>, ''x''<sub>2</sub>, …, ''x<sub>n</sub>'') &gt; 0 ;
la frontière entre les deux est une [[hypersurface]] d'équation ''g''(''x''<sub>1</sub>, ''x''<sub>2</sub>, …, ''x<sub>n</sub>'') = 0 qui définit l'[[États limites ultime et en service|état limite]].

Chaque variable aléatoire X<sub>''i''</sub> suit une loi de probabilité dont la fonction de [[densité de probabilité]] est notée ƒ<sub>''i''</sub>. La densité de probabilité d'avoir un ''n''-uplet donné (''x''<sub>1</sub>, ''x''<sub>2</sub>, …, ''x<sub>n</sub>'') vaut donc
: ''y'' = ƒ<sub>1</sub>(''x''<sub>1</sub>)×ƒ<sub>2</sub>(''x''<sub>2</sub>)×…×ƒ<sub>''n''</sub>(''x<sub>n</sub>'')
La fonction ''y''(''x''<sub>1</sub>, ''x''<sub>2</sub>, …, ''x<sub>n</sub>'') forme une hypersurface de l'espace (''x''<sub>1</sub>, …, ''x<sub>n</sub>'', ƒ). Le maximum de ''y'', le « sommet » de l'hypersurface, est le point le plus probable.

Pour déterminer la probabilité de défaillance P<sub>f</sub> du système, il faut intégrer ''y'' sur la zone de défaillance F
: <math>\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</math>
La valeur de P<sub>f</sub> 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 P<sub>f</sub> consiste à utiliser une [[méthode de Monte-Carlo]] : on génère un grand nombre de valeurs aléatoires (''x''<sub>1</sub>, ''x''<sub>2</sub>, …, ''x<sub>n</sub>'') 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 (''x''<sub>1</sub>, ''x''<sub>2</sub>, …, ''x<sub>n</sub>'') 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 P<sub>f</sub>.

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 P<sub>f</sub>.

== Méthode générale ==

La méthode FORM-SORM consiste à :
# Remplacer les variables aléatoires par des [[Loi normale#Loi normale centrée réduite|lois normales centrées réduites]].
# Trouver le ''n''-uplet (''x''<sub>1</sub>, ''x''<sub>2</sub>, …, ''x<sub>n</sub>'') 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.
# Faire une approximation linéaire (FORM) ou quadratique (SORM) de ''g'' pour déterminer la probabilité de défaillance P<sub>f</sub>.

=== Hypothèses ===

On suppose que les variables aléatoires ont une [[loi normale|distribution normale]].

Si les variables sont normales d'espérance μ<sub>''i''</sub> et d'écart type σ<sub>''i''</sub>, on les transforme également en variables de loi centrées réduites
: U<sub>''i''</sub> = (X<sub>''i''</sub> - μ<sub>''i''</sub>)/σ<sub>''i''</sub>
: X<sub>''i''</sub> = σ<sub>''i''</sub>U<sub>''i''</sub> + μ<sub>''i''</sub>

Si les variables aléatoires X<sub>''i''</sub> ne sont pas normales, on les transforme pour avoir des variables U<sub>''i''</sub> suivant lois normales centrées réduites ; cette opération est la « transformation de Nataf et Rosenblatt » :
: U<sub>''i''</sub> = Φ<sup>-1</sup>(F<sub>''i''</sub>(X<sub>''i''</sub>))
et réciproquement
: X<sub>''i''</sub> = F<sup>-1</sup>(Φ<sub>''i''</sub>(U<sub>''i''</sub>))
où Φ (lettre grecque Phi) est la [[fonction de répartition]] de la loi normale centrée réduite, et F<sub>''i''</sub> la fonction de répartition de X<sub>''i''</sub>.

=== Méthode générale ===

Plus on est proche du sommet de l'hypersurface ''y'', plus la densité de probabilité est élevée. Les fonctions de densité ƒ<sub>''i''</sub> é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''(''x''<sub>1</sub>, …, ''x<sub>n</sub>'') = 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''*(''u''<sub>1</sub>, …, ''u<sub>n</sub>'') ≥ 0.
La distance à l'origine est donc
: <math>d^* = \sqrt{u_1^2 + u_2^2 + \cdots + u_n^2}</math>
et ainsi,
: la méthode consiste à minimiser ''d''*(''u''<sub>1</sub>, …, ''u<sub>n</sub>'') avec pour contrainte ''g''*(''u''<sub>1</sub>, …, ''u<sub>n</sub>'') = 0.
C'est une méthode d'[[optimisation quadratique]] : en effet, minimiser ''d''* revient à minimiser ''d''*<sup>2</sup> (puisque la fonction carré est strictement croissante sur [0 ; +∞[), on doit donc minimiser une fonction quadratique ∑''x<sub>i</sub>''<sup>2</sup> 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 ==
On fait l'hypothèse que la fonction ''g'' est [[application linéaire|linéaire]]. Si ce n'est pas le cas, la méthode revient finalement à en prendre le [[développement limité]] du premier ordre.

La probabilité de défaillance vaut alors
: P<sub>f</sub> = Φ(-β).

{{démonstration
 | contenu = Plaçons-nous à deux dimensions ; nous considérons donc un problème à deux variables aléatoires X<sub>1</sub> et X<sub>2</sub>, et les variables normales centrées réduites correspondantes sont U<sub>1</sub> et U<sub>2</sub>. Le graphe de l'état limite ''g''*(''u''<sub>1</sub>, ''u''<sub>2</sub>) = 0 est donc une droite.

Comme la densité de probabilité ƒ = ƒ<sub>1</sub>×ƒ<sub>2</sub> est à symétrie centrale de centre (0 ; 0), on peut « tourner la figure » ; supposons donc que cette droite ''g''*(''u''<sub>1</sub>, ''u''<sub>2</sub>) = 0 est verticale et située du côté des abscisses négatives, elle a pour équation
: ''x'' = -β.
L'intersection de la surface ƒ avec le plan vertical d'abscisse ''x'' est une courbe gaussienne
: ''z<sub>x</sub>''(''y'') = ƒ<sub>1</sub>(''x'')×ƒ<sub>2</sub>(''y'') = φ(''x'')×φ(''y'')
où φ est le densité de probabilité de la loi normale centrée réduite. Le volume dV situé sous la portion de surface de ƒ située entre les abscisses ''x'' et ''x'' + d''x'' vaut :
: <math>\mathrm{dV} = \int_{y = -\infty}^{+\infty} \varphi(x) \times \varphi(y) \times \mathrm{d}x \times \mathrm{d} y
 = \varphi(x) \mathrm{d}x \int_{y = -\infty}^{+\infty} \varphi(y) \mathrm{d} y = \varphi(x) \mathrm{d}x</math>
et donc, le volume -β sous la courbe ƒ situé entre les abscisses -∞ et -β vaut :
: <math>\mathrm{V} = \int_{x = -\infty}^{-\beta} \varphi(x) \mathrm{d}x = \Phi(-\beta)</math>
on a donc bien
: Pf = Φ(-β).
}}

== Méthode SORM ==

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 P<sub>f</sub>.

On définit les [[courbure]]s principales κ<sub>''i''</sub> (lettre grecque kappa) :
: <math>\kappa_i = \frac{\partial^2 g^*}{\partial u_i^2}</math>
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 :
: P<sub>f</sub> = S<sub>1</sub> + S<sub>2</sub> + S<sub>3</sub>
avec :
: <math>\mathrm{S}_1 = \Phi(-\beta) \prod_{j = 1}^{n - 1} \frac{1}{\sqrt{1 + \beta \kappa_j}}</math>
: <math>\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 )</math>
: <math>\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 ]</math>
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 S<sub>1</sub> est l'approximation asymptotique de P<sub>f</sub> 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 S<sub>2</sub> et S<sub>3</sub> sont des termes correctifs.

== Application à un problème à deux variables normales ==

=== Présentation du problème ===

[[Fichier:Contrainte resistance 3d proche.svg|vignette|300px|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 &gt; S).]]

[[Fichier:Contrainte resistance niveaux proche.svg|vignette|300px|Même situation vue de dessus.]]

[[Fichier:Contrainte resistance 2d proche.svg|vignette|300px|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é de probabilité|densités de probabilité]]) ƒ<sub>R</sub> et ƒ<sub>S</sub>, 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]) = ƒ<sub>S</sub>×ƒ<sub>R</sub>×dS×dR.
On peut tracer la surface ''y'' = ƒ<sub>S</sub>×ƒ<sub>R</sub>, 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 μ<sub>S</sub> = {{unité|230|MPa}} et un écart type σ<sub>S</sub> = {{unité|14|Mpa}} ;
* la résistance R suit une loi normale d'espérance μ<sub>R</sub> = {{unité|250|MPa}} et un écart type σ<sub>R</sub> = {{unité|17|Mpa}}.
Graphiquement, on constate que le maximum de l'état limite se trouve à {{unité|238|MPa}}. Ce point est en général proche de l'intersection des courbes ƒ<sub>S</sub> et ƒ<sub>R</sub>.

Une méthode de Monte-Carlo avec un million de tirages et un intervalle de confiance de 5 %<ref>l'erreur pour un risque α avec N tirages vaut<br /> <math>\varepsilon = \sqrt{\frac{1}{2(\mathrm{N} + 1)} \cdot \ln \left ( \frac{2}{\alpha}\right )}</math> <br /> soit <br /> <math>\varepsilon = \sqrt{\frac{1}{2(10^6 + 1)} \cdot \ln \left ( \frac{2}{0,005}\right )} \simeq 1,7 \cdot 10^{-3} = 0,17\ \%</math></ref> donne 18,1 % ± 0,2 % de défaillance.

{{boîte déroulante/début | titre = Résolution avec LibreOffice Calc}}
Le logiciel [[LibreOffice Calc]] est un tableur libre et gratuit.

La fonction <code>ALEA</code> (ou <code>RAND</code> en anglais)<ref>[https://help.libreoffice.org/Calc/Mathematical_Functions/fr#ALEA Fonctions mathématiques] (aide lLibreOffice)</ref> donne un nombre aléatoire compris entre 0 et 1, mais avec une loi uniforme. Pour avoir une loi normale, il faut appliquer la fonction <code>LOI.NORMALE.INVERSE</code> (<code>NORMINV</code>)<ref>[https://help.libreoffice.org/3.3/Calc/Statistical_Functions_Part_Four/fr#LOI.NORMALE.INVERSE Fonctions statistiques - Quatrième partie] (aide LibreOffice 3.3)</ref>. On compare simplement deux valeurs avec l'opérateur <code>&gt;</code>, ce qui donne un booléen ; Calc calcule la moyenne (commande <code>MOYENNE</code>, <code>MEAN</code>) de booléens en prenant 1 pour « vrai » et 0 pour « faux ». Il faut ensuite étirer une boîte de… un million de lignes, ce qui n'est pas raisonnable, et touche les limites du logiciel<ref>[https://wiki.documentfoundation.org/FR/FAQ/Calc/022 Quel est le nombre maximum de cellules dans une feuille de calcul LibreOffice ?], FAQ LibreOffice Calc</ref>.

On se limite ici à mille lignes, ce qui donne un résultat à ± 5 %.

{| border = "1"
|-
!
! scope = "col" | A
! scope = "col" | B
! scope = "col" | C
! scope = "col" | D
|-
! scope = "row" | 1
| || R || S || 
|-
! scope = "row" | 2
| mu || 250 || 230 || 
|-
! scope = "row" | 3	
| sigma || 17 || 14 || 
|-
! scope = "row" | 4	
| ||  ||  || 
|-
! scope = "row" | 5
| n° || Valeur R || Valeur S || S>R
|-
! scope = "row" | 6
| 1
| <code>=LOI.NORMALE.INVERSE(ALEA() ; B$2 ; B$3)</code>
| <code>=LOI.NORMALE.INVERSE(ALEA() ; C$2 ; C$3)</code>
| <code>=C6>B6</code>
|-
| colspan="5" | …
|-
! scope = "row" | 1006
| ||  ||  || <code>=MOYENNE(D6:D1005)</code>
|-
! scope = "row" | 1007
| ||  ||  || 
|-
! scope = "row" | 1008
| epsilon || <code>=RACINE(LN(2/0,005)/(2*1001))</code> || || 
|}

{{boîte déroulante/fin}}

{{boîte déroulante/début | titre = Résolution avec R}}
Le logiciel [[R (langage de programmation et environnement statistique)|R]] est un logiciel libre et gratuit, sous la forme d'un langage de programmation interprété, qui permet d'effectuer des calculs statistiques. La méthode de Monte-Carlo se fait de la manière suivante :
# On crée deux matrices R et S de un million de nombres tirés aléatoirement avec une loi normale (commandes <code>rnorm()</code>).
# On crée une matrice booléenne Mbool, contenant des « vrai » si la contrainte S est supérieure à la résistance R (opérateur <code>&gt;</code>).
# On calcule la moyenne de cette matrice, les « vrai » étant comptés comme des 1 et les « faux » comme des 0 (commande <code>mean()</code>).

<syntaxhighlight lang="rsplus">
N = 1e6 # nombre d'événements (tirs)
alpha = 0.005 # risque

# paramètres des lois
muR = 250
sigmaR = 17
muS = 230
sigmaS = 14

# calculs
R = rnorm(N, muR, sigmaR)
S = rnorm(N, muS, sigmaS) 
Mbool = S>R

moyenne = mean(Mbool)

epsilon = sqrt(log(2/alpha)/(2*(N + 1)))

# affichage des résultats
print(moyenne)
print(epsilon)
</syntaxhighlight>
{{boîte déroulante/fin}}

{{boîte déroulante/début | titre = Résolution avec Scilab}}
[[Scilab]] est un logiciel libre et gratuit de calcul scientifique. Il se présente comme un langage interprété.

Le tirage aléatoire de R et de S se fait avec la fonction <code>rand()</code> ; l'option <code>"normal"</code> indique qu'il utilise une loi normale centrée réduite. On détermine une matrice booléenne Mbool contenant des « vrai » si S est plus grand que R, avec l'opérateur <code>&gt;</code>. La fonction <code>sum()</code> calcule la somme des éléments de la matrice en considérant les « vrai » comme des 1 et les « faux » comme des 0.

<syntaxhighlight lang="scilab">
N = 1e6; // nombre d'événements (tirs)
alpha = 0.005; // risque

// Paramètres des lois
muR = 250;
sigmaR = 17;
muS = 230;
sigmaS = 14;

// Calculs
R = muR + sigmaR*rand(N, 1, "normal"); // résistance
S = muS + sigmaS*rand(N, 1, "normal"); // contrainte
Mbool = S>R; // comparaison des valeurs (matrice de booléens)
moyenne = sum(Mbool)/N; // les "vrai" comptent pour 1 et les "faux" pour 0
epsilon = sqrt(log(2/alpha)/(2*(N + 1)));

// Affichage des résultats
disp(string(100*moyenne)+" % +/- "+...
    +string(100*epsilon)+" %") // affichage du résultat
</syntaxhighlight>
{{boîte déroulante/fin}}

{{clr}}

=== Mise en œuvre de la méthode FORM ===

[[Fichier:Form-sorm contours noir.svg|vignette|300px|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 [[Loi normale#Loi normale centrée réduite|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* :
: <math>\mathrm{R}^* = \frac{\mathrm{R} - \mu_{\mathrm{R}}}{\sigma_{\mathrm{R}}}</math> ;
: <math>\mathrm{S}^* = \frac{\mathrm{S} - \mu_{\mathrm{S}}}{\sigma_{\mathrm{S}}}</math>.
Le sommet de la surface ƒ<sub>S*</sub>×ƒ<sub>R*</sub> se trouve donc aux coordonnées (0 ; 0).
La droite d'état limite ''m'' = 0 (R - S = 0) devient donc :
: <math>\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</math>
ce que l'on peut noter
: ''m''* = 0
avec
: ''m''*(R*, S*) = σ<sub>R</sub>×R* + μ<sub>R</sub> - (σ<sub>S</sub>×S* + μ<sub>S</sub>).

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 ƒ<sub>S*</sub>×ƒ<sub>R*</sub>, c'est aussi le plus proche de l'origine (0 ; 0) — sommet de la surface —, c'est-à-dire celui pour lequel la quantité
: ''d''*<sup>2</sup> = R*<sup>2</sup> + S*<sup>2</sup>
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 P<sub>f</sub> s'obtient avec la fonction de répartition Φ de la loi normale centrée réduite :
: P<sub>f</sub> = Φ(-β) ;
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 :
: <math>\mathrm{R}^* = \frac{\mathrm{R} - 250}{17}</math> ;
: <math>\mathrm{S}^* = \frac{\mathrm{S} - 230}{14}</math>.
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 = {{unité|238|MPa}}
: β = 0,908 ;
: P<sub>f</sub> = 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 [[optimisation non linéaire|solveur non-linéaire]] « généraliste » ; il faut donc leur définir la fonction ''g'', ainsi que la contrainte
: ''h'' = 0, avec ''h'' = σ<sub>R</sub>×R* - σ<sub>S</sub>×S* + (μ<sub>R</sub> - μ<sub>S</sub>).
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énalisation (optimisation)|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 <code>constrOptim()</code> 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 <code>solve.QP</code> du langage R, et de la fonction <code>qpsolve</code> de Scilab. Le problème s'écrit alors de manière matricielle
: <math>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}</math>
avec
: <math>\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}</math>
avec une contrainte d'égalité unique :
: C×'''x''' = ''b''
avec
: C = (σ<sub>R</sub> ; -σ<sub>S</sub>) ; ''b'' = (μ<sub>S</sub> - μ<sub>R</sub>).

{{boîte déroulante/début| titre = Résolution directe}}
Le problème est ici relativement simple, nous pouvons donc le résoudre « à la main ». Nous évoluons sur une droite ''m''* = 0, soit d'équation

{| border="0"
|-
! scope = "col" | expression littérale
| &nbsp;
! scope = "col" | application numérique
|-
| <math>\sigma_\mathrm{R} \times \mathrm{R}^* - \sigma_\mathrm{S} \times \mathrm{S}^*
+ (\mu_\mathrm{R} - \mu_\mathrm{S})</math>
|
| 17×R* - 14×S* + 20 = 0
|-
| <math>\Leftrightarrow
\mathrm{S}^* = \frac{\sigma_\mathrm{R}}{\sigma_\mathrm{S}} \mathrm{R}^*
+ \frac{\mu_\mathrm{R} - \mu_\mathrm{S}}{\sigma_\mathrm{S}}</math>
|
| ⇔ S* = (17/14)×R* + 20/14
|}

On cherche à minimiser ''d''*<sup>2</sup> = R*<sup>2</sup> + S*<sup>2</sup>, qui s'écrit

{| border="0"
|-
! scope = "col" | expression littérale
| &nbsp;
! scope = "col" | application numérique
|-
| <math>d^{*\,2} = \mathrm{R}^{*\, 2}
+ \left ( \frac{\sigma_\mathrm{R}}{\sigma_\mathrm{S}} \mathrm{R}^*
+ \frac{\mu_\mathrm{R} - \mu_\mathrm{S}}{\sigma_\mathrm{S}} \right )^2</math>
|
| ''d''*<sup>2</sup> = R*<sup>2</sup> + ((17/14)×R* + 20/14)<sup>2</sup>
|}

soit

{| border="0"
|-
! scope = "col" | expression littérale
| &nbsp;
! scope = "col" | application numérique
|-
| <math>d^{*\,2} = \left ( 1 + \frac{\sigma_\mathrm{R}^2}{\sigma_\mathrm{S}^2} \right )
\mathrm{R}^{*\,2}
+ \frac{2 \sigma_\mathrm{R}(\mu_\mathrm{R} - \mu_\mathrm{S})}{\sigma_\mathrm{S}^2}\mathrm{R}^*
+ \frac{(\mu_\mathrm{R} - \mu_\mathrm{S})^2}{\sigma_\mathrm{S}^2}</math>
|
| ''d''*<sup>2</sup> ≃ 2,474×R*<sup>2</sup> + 3,469×R* + 2,041
|}

Le minimum est là où la dérivée s'annule : (''d''*<sup>2</sup>){{'}} = 0 ⇔
{| border="0"
|-
! scope = "col" | expression littérale
| &nbsp;
! scope = "col" | application numérique
|-
| <math>2 \left ( 1 + \frac{\sigma_\mathrm{R}^2}{\sigma_\mathrm{S}^2} \right )\mathrm{R}^*
+ \frac{2 \sigma_\mathrm{R}(\mu_\mathrm{R} - \mu_\mathrm{S})}{\sigma_\mathrm{S}^2}
= 0</math>
|
| 4,948×R* + 3,469 = 0 
|-
| <math>\Leftrightarrow \mathrm{R}^*
= -\frac{\sigma_\mathrm{R}(\mu_\mathrm{R} - \mu_\mathrm{S})}{\sigma_\mathrm{S}^2 + \sigma_\mathrm{R}^2}</math>
|
| ⇔ R* = -0,7010
|}

et donc
{| border="0"
|-
! scope = "col" | expression littérale
| &nbsp;
! scope = "col" | application numérique
|-
| <math>d^* = \frac{\mu_\mathrm{R} - \mu_\mathrm{S}}{\sqrt{\sigma_\mathrm{S}^2 + \sigma_\mathrm{R}^2}}</math>
|
| ''d''* = 0,908
|}

La [[Loi normale#Tables numériques et calculs|table de valeurs de la fonction de répartition de la loi normale]] donne
: Φ(-0,908) = 1 - Φ(0,908) ≃ 1 - Φ(0,91) = 1 - 0,81859 ≃ 0,181.
{{boîte déroulante/fin}}

{{boîte déroulante/début| titre = Résolution avec LibreOffice Calc}}
[[Fichier:Solveur libreoffice calc 3 fr.png|vignette|Interface du solveur de LibreOffice Calc 3.]]

Le tableur libre [[LibreOffice Calc]] dispose d'un solveur.

On construit le tableau suivant :

{| border="1px"
|-
|
! scope = "col" | A
! scope = "col" | B
! scope = "col" | C
! scope = "col" | D
|-
! scope = "row" | 1
| muR || 250 || sigmaR || 17
|-
! scope = "row" | 2
| muS || 230 || sigmaS || 14
|-
! scope = "row" | 3
| || || ||
|-
! scope = "row" | 4
| R || 250 || ||
|-
! scope = "row" | 5
| S || 230 || ||
|-
! scope = "row" | 6
| R* || <code>=(B4 - B1)/D1</code> || ||
|-
! scope = "row" | 7
| S* || <code>=(B5 - B2)/D2</code> || ||
|-
! scope = "row" | 8
| d*² || <code>= B6^2 + B7^2</code> || ||
|-
! scope = "row" | 9
| m || <code>= B4 - B5</code> || ||
|}
Puis on va dans le menu '''Outil''' &gt; '''Solveur''', et dans la boîte de dialogue '''Solveur''' qui apparaît, on met :
* '''Cellule cible''' : <code>B8</code>
* '''Optimiser le résultat à''' : <code>minimum</code>
* '''Par modification de cellules''' : <code>B4:B5</code>
* '''Conditions de limitation''' :
** '''Référence de cellule''' : <code>B9</code>
** '''Opérateur''' : <code>=</code>
** '''Valeur''' : <code>0</code>
{{boîte déroulante/fin}}

{{boîte déroulante/début| titre = Résolution avec R}}
Le solveur d'optimisation utilisé est la commande <code>constrOptim</code>. C'est un solveur non-linéaire qui applique les contraintes par la méthode de la « barrière logarithmique ». La fonction optimisée est
: ''d''*<sup>2</sup> = R*<sup>2</sup> + S*<sup>2</sup>
La contrainte linéaire est sous la forme
: ''u''<sub>1</sub>×''x'' - ''c''<sub>1</sub> ≥ 0
Bien que la méthode requiert une égalité dans la contrainte, l'inégalité est ici suffisante : le point le plus proche de l'origine est nécessairement sur la frontière ''g''* = 0, il n'est pas nécessaire de contraindre la solution à rester sur cette frontière, elle s'y retrouvera nécessairement. Il faut par contre choisir le bon sens pour l'inégalité, pour être du côté des distances « élevées », donc la valeur calculée doit être positive pour S* &gt; R*. Il faut donc prendre pour inégalité
: - σ<sub>R</sub>×R* + σ<sub>S</sub>×S* - (μ<sub>R</sub> - μ<sub>S</sub>) ≥ 0
donc ''u''<sub>1</sub> est le vecteur ligne (-σ<sub>R</sub> ; σ<sub>S</sub>), ''x'' est le vecteur colonne (R* ; S*) et ''c''<sub>1</sub> vaut μ<sub>R</sub> - μ<sub>S</sub>.

On lui passe en paramètres :
* le point de départ theta, qui doit être dans la « zone de faisabilité » ''u''<sub>1</sub>×''x'' - ''c''<sub>1</sub> ≥ 0 ; on choisit arbitrairement le point (0 ; 2) ;
* la fonction à optimiser ;
* <code>NULL</code>, indiquant que l'on ne donne pas l'expression du gradient ;
* les paramètres de la contrainte linéaire, ui et ci ;
* un paramètre de précision, mu.
Le résultat <code>a</code> est sous la forme d'une liste de vecteurs :
* <code>a$par</code> est le vecteur (R*, S*) trouvé, c'est-à-dire le point de conception ;
* <code>a$value</code> est la valeur de la fonction en ce point, β<sup>2</sup>.

<syntaxhighlight lang="rsplus">
# paramètres de la loi
muR <- 250
sigmaR <- 17
muS <- 230
sigmaS <- 14

# point de départ (R, S) de l'optimisation
theta <- rbind(0, 2)

# fonction à optimiser
detoile2 <- function(x) { # carré de la distance à l'origine, d*²
  Retoile <- x[1]
  Setoile <- x[2]
  Retoile^2 + Setoile^2
  # aperm(x) %*% x
}

# contrainte linéaire
u1 <- c(-sigmaR, sigmaS) 
c1 <- muR - muS

# calcul
a <- constrOptim(theta, detoile2, NULL, u1, c1, mu=1e-3) # optimisation

x <- a$par # (R*, S*)
R <- x[1]*sigmaR + muR
#S <- x[2]*sigmaS + muS
detoile <- sqrt(a$value)
Pf <- pnorm(-detoile)

# affichage des résultats
print(paste("R* = ", x[1], " ; S* = ", x[2]))
print(paste("R = S = ", R))
print(paste("beta = ", detoile))
print(paste("Pf = ", Pf))
</syntaxhighlight>

S'agissant d'un problème de programmation quadratique, on peut également utiliser la fonction <code>solve.QP</code> de l'extension <code>quadprog</code> (qui est normalement installé par défaut).

<syntaxhighlight lang="rsplus">
#install.packages("quadprog")
library("quadprog")

# paramètres de la loi
muR <- 250
sigmaR <- 17
muS <- 230
sigmaS <- 14

# Optimisation
Dmat <- rbind(c(1, 0), c(0, 1))
dvec <- c(0, 0)
Cmat <- cbind(c(-sigmaR, sigmaS))
bvec <- c(muR - muS)
a <- solve.QP(Dmat, dvec, Cmat, bvec)

x <- a$solution # (R*, S*)
R <- x[1]*sigmaR + muR
#S <- x[2]*sigmaS + muS
detoile <- sqrt(2*a$value)
Pf <- pnorm(-detoile)

# affichage des résultats
print(paste("R* = ", x[1], " ; S* = ", x[2]))
print(paste("R = S = ", R))
print(paste("beta = ", detoile))
print(paste("Pf = ", Pf))
</syntaxhighlight>
{{boîte déroulante/fin}}

{{boîte déroulante/début| titre = Résolution avec Scilab}}
Le logiciel libre [[Scilab]] dispose de fonctions d'optimisation.

C'est un cas de [[programmation quadratique]] avec une contrainte d'égalité unique (''m''<sub>e</sub> = 1) ; il n'y a pas de limite inférieure ni supérieure. Le script utilisé est donc :

<syntaxhighlight lang="scilab">
// Données

muR = 250;
sigmaR = 17;
muS = 230;
sigmaS = 14;

// Optimisation

Q = [1, 0 ; 0, 1];
p = [0 ; 0];
C = [sigmaR, -sigmaS];
b = muS - muR;
ci = [];
cs = [];
me = 1;

[x] = qpsolve(Q, p, C, b, ci, cs, me); // Solveur

// Exploitation des données

R = x(1)*sigmaR + muR;
d = norm(x);
Pf = cdfnor("PQ", -d, 0, 1); // fonction Phi

// Affichage

disp("R* = "+string(x(1))+" ; S* = "+string(x(2)));
disp("R = S = "+string(R));
disp("beta = "+string(d));
disp("Pf = "+string(Pf));
</syntaxhighlight>
{{boîte déroulante/fin}}

== Notes et références ==

{{références}}

== Voir aussi ==

=== Bibliographie ===

* {{article
 | prénom1 = Shankar | nom1 = Sankararaman
 | prénom2 = Kai | nom2 = Goebel
 | titre = Uncertainty Quantification in Remaining Useful Life of Aerospace Components using State Space Models and Inverse FORM
 | langue = en
 | mois = mai | année = 2013
 | périodique = IEEE Aerospace Conference
 | éditeur = IEEE | lien éditeur = Institute of Electrical and Electronics Engineers
 | lire en ligne = http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20130013373_2013013151.pdf
}}
* {{article
 | prénom1 = L. | nom1 = Cizelj
 | prénom2 = B. | nom2 = Mavko
 | prénom3 = H. | nom3 = Riesch-Oppermann
 | titre = Application of first and second order reliability methods in the safety assessment of cracked steam generator tubing
 | langue = en
 | périodique = Nuclear Engineering and Design
 | éditeur = Elsevier | lien éditeur = Elsevier (éditeur)
 | année = 1994
 | numéro = 147
 | lire en ligne = http://www.ijs.si/~cizelj/papers/formsorm.pdf
}}

{{Portail|probabilités et statistiques}}

{{DEFAULTSORT:Methode form sorm}}
[[Catégorie:Fiabilité]]
[[Catégorie:Optimisation quadratique]]