Méthode FORM-SORM

La méthode FORM-SORM est une méthode utilisée en fiabilité pour déterminer la 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 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 α, 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) = α.

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. 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 x_i. Ces valeurs sont des réalisations de variables aléatoires X_i : 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₁, x₂, …, x_n) pour exprimer la condition sur la durée de vie :

F(t) ≤ α ⇔ g(x₁, x₂, …, x_n) ≥ 0.

L'espace à n dimensions (x₁, x₂, …, x_n) est donc coupé en deux demi-espaces :

le domaine de défaillance F (failure), vérifiant l'équation g(x₁, x₂, …, x_n) < 0 ;
le domaine de bon fonctionnement, vérifiant l'équation g(x₁, x₂, …, x_n) > 0 ;

la frontière entre les deux est une hypersurface d'équation g(x₁, x₂, …, x_n) = 0 qui définit l'état limite.

Chaque variable aléatoire X_i 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é (x₁, x₂, …, x_n) vaut donc

y = ƒ₁(x₁)×ƒ₂(x₂)×…×ƒ_n(x_n)

La fonction y(x₁, x₂, …, x_n) forme une hypersurface de l'espace (x₁, …, x_n, ƒ). Le maximum de y, le « sommet » de l'hypersurface, est le point le plus probable.

Pour déterminer la probabilité de défaillance P_f 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 P_f 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_f consiste à utiliser une méthode de Monte-Carlo : on génère un grand nombre de valeurs aléatoires (x₁, x₂, …, x_n) 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₁, x₂, …, x_n) 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_f.

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_f.

Méthode générale

La méthode FORM-SORM consiste à :

Remplacer les variables aléatoires par des lois normales centrées réduites.
Trouver le n-uplet (x₁, x₂, …, x_n) 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_f.

Hypothèses

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

U_i = (X_i - μ_i)/σ_i

X_i = σ_iU_i + μ_i

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

U_i = Φ^-1(F_i(X_i))

et réciproquement

X_i = F^-1(Φ_i(U_i))

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

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é ƒ_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(x₁, …, x_n) = 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₁, …, u_n) ≥ 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*(u₁, …, u_n) avec pour contrainte g*(u₁, …, u_n) = 0.

C'est une méthode d'optimisation quadratique : en effet, minimiser d* revient à minimiser d*² (puisque la fonction carré est strictement croissante sur [0 ; +∞[), on doit donc minimiser une fonction quadratique ∑x_i² 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 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_f = Φ(-β).

Démonstration

Plaçons-nous à deux dimensions ; nous considérons donc un problème à deux variables aléatoires X₁ et X₂, et les variables normales centrées réduites correspondantes sont U₁ et U₂. Le graphe de l'état limite g*(u₁, u₂) = 0 est donc une droite.

Comme la densité de probabilité ƒ = ƒ₁×ƒ₂ est à symétrie centrale de centre (0 ; 0), on peut « tourner la figure » ; supposons donc que cette droite g*(u₁, u₂) = 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_x(y) = ƒ₁(x)×ƒ₂(y) = φ(x)×φ(y)

où φ est la 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 + dx vaut :

\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

et donc, le volume -β sous la courbe ƒ situé entre les abscisses -∞ et -β vaut :

\mathrm {V} =\int _{x=-\infty }^{-\beta }\varphi (x)\mathrm {d} x=\Phi (-\beta )

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_f.

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 :

P_f = S₁ + S₂ + S₃

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 S₁ est l'approximation asymptotique de P_f 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₂ et S₃ sont des termes correctifs.

Application à un problème à deux variables normales

Présentation du problème

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.

Résolution avec LibreOffice Calc

Le logiciel LibreOffice Calc est un tableur libre et gratuit.

La fonction ALEA (ou RAND en anglais)^[2] 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 LOI.NORMALE.INVERSE (NORMINV)^[3]. On compare simplement deux valeurs avec l'opérateur >, ce qui donne un booléen ; Calc calcule la moyenne (commande MOYENNE, MEAN) 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^[4].

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

	A	B	C	D
1		R	S
2	mu	250	230
3	sigma	17	14
4
5	n°	Valeur R	Valeur S	S>R
6	1	`=LOI.NORMALE.INVERSE(ALEA() ; B$2 ; B$3)`	`=LOI.NORMALE.INVERSE(ALEA() ; C$2 ; C$3)`	`=C6>B6`
…
1006				`=MOYENNE(D6:D1005)`
1007
1008	epsilon	`=RACINE(LN(2/0,005)/(2*1001))`

Résolution avec R

Le logiciel 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 rnorm()).
On crée une matrice booléenne Mbool, contenant des « vrai » si la contrainte S est supérieure à la résistance R (opérateur >).
On calcule la moyenne de cette matrice, les « vrai » étant comptés comme des 1 et les « faux » comme des 0 (commande mean()).

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)

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 rand() ; l'option "normal" 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 >. La fonction sum() calcule la somme des éléments de la matrice en considérant les « vrai » comme des 1 et les « faux » comme des 0.

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

Mise en œuvre de la méthode FORM

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*² = R*² + S*²

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_f s'obtient avec la fonction de répartition Φ de la loi normale centrée réduite :

P_f = Φ(-β) ;

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 ;

P_f = 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 :

C×x = b

avec

C = (σ_R ; -σ_S) ; b = (μ_S - μ_R).

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

expression littérale		application numérique
$\sigma _{\mathrm {R} }\times \mathrm {R} ^{}-\sigma _{\mathrm {S} }\times \mathrm {S} ^{}+(\mu _{\mathrm {R} }-\mu _{\mathrm {S} })$		17×R* - 14×S* + 20 = 0
$\Leftrightarrow \mathrm {S} ^{}={\frac {\sigma _{\mathrm {R} }}{\sigma _{\mathrm {S} }}}\mathrm {R} ^{}+{\frac {\mu _{\mathrm {R} }-\mu _{\mathrm {S} }}{\sigma _{\mathrm {S} }}}$		⇔ S* = (17/14)×R* + 20/14

On cherche à minimiser d*² = R*² + S*², qui s'écrit

expression littérale		application numérique
$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}$		d² = R² + ((17/14)×R* + 20/14)²

soit

expression littérale		application numérique
$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}}}$		d² ≃ 2,474×R² + 3,469×R* + 2,041

Le minimum est là où la dérivée s'annule : (d*²)' = 0 ⇔

expression littérale		application numérique
$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$		4,948×R* + 3,469 = 0
$\Leftrightarrow \mathrm {R} ^{*}=-{\frac {\sigma _{\mathrm {R} }(\mu _{\mathrm {R} }-\mu _{\mathrm {S} })}{\sigma _{\mathrm {S} }^{2}+\sigma _{\mathrm {R} }^{2}}}$		⇔ R* = -0,7010

et donc

expression littérale		application numérique
$d^{*}={\frac {\mu _{\mathrm {R} }-\mu _{\mathrm {S} }}{\sqrt {\sigma _{\mathrm {S} }^{2}+\sigma _{\mathrm {R} }^{2}}}}$		d* = 0,908

La 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.

Résolution avec LibreOffice Calc

Le tableur libre LibreOffice Calc dispose d'un solveur.

On construit le tableau suivant :

	A	B	C	D
1	muR	250	sigmaR	17
2	muS	230	sigmaS	14
3
4	R	250
5	S	230
6	R*	`=(B4 - B1)/D1`
7	S*	`=(B5 - B2)/D2`
8	d*²	`= B6^2 + B7^2`
9	m	`= B4 - B5`

Puis on va dans le menu Outil > Solveur, et dans la boîte de dialogue Solveur qui apparaît, on met :

Cellule cible : B8
Optimiser le résultat à : minimum
Par modification de cellules : B4:B5
Conditions de limitation :
- Référence de cellule : B9
- Opérateur : =
- Valeur : 0

Résolution avec R

Le solveur d'optimisation utilisé est la commande constrOptim. 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*² = R*² + S*²

La contrainte linéaire est sous la forme

u₁×x - c₁ ≥ 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* > R*. Il faut donc prendre pour inégalité

- σ_R×R* + σ_S×S* - (μ_R - μ_S) ≥ 0

donc u₁ est le vecteur ligne (-σ_R ; σ_S), x est le vecteur colonne (R* ; S*) et c₁ vaut μ_R - μ_S.

On lui passe en paramètres :

le point de départ theta, qui doit être dans la « zone de faisabilité » u₁×x - c₁ ≥ 0 ; on choisit arbitrairement le point (0 ; 2) ;
la fonction à optimiser ;
NULL, 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 a est sous la forme d'une liste de vecteurs :

a$par est le vecteur (R*, S*) trouvé, c'est-à-dire le point de conception ;
a$value est la valeur de la fonction en ce point, β².

# 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))

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

#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))

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_e = 1) ; il n'y a pas de limite inférieure ni supérieure. Le script utilisé est donc :

// 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));

Notes et références

↑ 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\ \%$
↑ Fonctions mathématiques (aide lLibreOffice)
↑ Fonctions statistiques - Quatrième partie (aide LibreOffice 3.3)
↑ Quel est le nombre maximum de cellules dans une feuille de calcul LibreOffice ?, FAQ LibreOffice Calc

Voir aussi

Bibliographie

(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, n^o 147,‎ 1994 (lire en ligne)