Propagation des incertitudes portées par la loi de comportement

Les méthodes de propagation d’incertitudes portées par la loi de comportement peuvent se diviser en 2 catégories : les méthodes non intrusives et les méthodes intrusives. Pour les méthodes de type non intrusives, on s’intéresse à une grandeur G qui peut être le flux traversant une surface, l’énergie du système, les pertes ou la valeur du potentiel en certains points… Cette grandeur est une variable aléatoire dépendant du vecteur des inconnues Ω′( ) ξ de (1.59) : G G ( ) ξ = ( )) (Ω′ ξ (1.61) La caractérisation de G(ξ) est réalisée grâce à un nombre d’évaluations de G : G(ξ k ), k=1 : n. Il faut alors réaliser n calculs déterministes de Ω′( ) ξ en résolvant le système (1.59) correspondant au ξ = ξ k . Puis les évaluations de G(ξ k ) sont obtenues par (1.61). Pour les méthodes non intrusives, on ajoute une « couche » au modèle déterministe permettant de prendre en compte les aspects probabilistes. Concernant les méthodes de type intrusive, le système (1.59) se ramène à un nouveau système à résoudre pour obtenir une expression de Ω′( ) ξ . Dans la suite, on présentera les méthodes de résolution non intrusives et intrusives d’un problème stochastique.

 Méthodes non intrusives

Dans cette partie, on présente trois méthodes non intrusives pour caractériser la grandeur G(ξ). 

Méthode de Monte Carlo

Une manière simple d’introduire la méthode de Monte Carlo est de considérer le théorème central limite (TCL). Il est donné par: Théorème 2 : Soit (Yi , i ≥1) une suite de variables aléatoires réelles indépendantes de carrée intégrable de même loi d’espérance m et d’écart-type σ. Propagation des incertitudes portées par la loi de comportement 32 Alors la suite ( ) n n Y m σ − , avec 1 1 n n i i Y Y n = = ∑ , converge en loi vers une variable de loi gaussienne centrée réduite Z,. Ce qui signifie que : lim ( ( ) ) ( ) n n n Y m a Z a σ →∞ P P − < = < (1.62) On peut déduire à partir de ce théorème que : 2 2 1 lim ( ) lim ( ( ) ) ( ) 2 b x n n n n n a n Y b m Y a a Y m b a Z b e dx n n σ σ σ π − →+∞ →+∞ − < < − = < − < = < < = ∫ P P P Par conséquent, pour estimer la moyenne m d’une variable aléatoire Y, une possibilité est d’obtenir un nombre n suffisamment grand de réalisations yi , i=1 : n, de Y et de supposer que ( ) n n Y m σ − suit effectivement une loi normale centrée réduite. L’intervalle de confiance à 95% par exemple de m est le suivant: 1 1 1 1 [ 1.96 , 1.96 ] n n i i i i y y n n n n σ σ = = ∑ ∑ − ⋅ + ⋅ Dans le cas où on ne connaît pas l’écart type σ, celui-ci peut être estimé par la formulation suivante : 2 2 1 1 1 1 ( ) 1 n n n i k i k y y n n σ = = = − − ∑ ∑ Cette estimation vient du fait que 2 2 1 1 1 1 ( ) 1 n n n i j i j Y Y n n σ = = = − − ɶ ∑ ∑ est un estimateur sans biais de σ 2 et 2 σ n ɶ converge presque surement vers σ 2 . Ce qui signifie que : 2 2 2 2 ( ) et (lim ) 1 σ σ σ σ n n n = = = →∞ E ɶ ɶ P En remplaçant σ par σ n , on peut approcher un intervalle de confiance à 95% de m : 1 1 1 1 [ 1.96 , 1.96 ] n n n n i i i i y y n n n n σ σ = = ∑ ∑ − ⋅ + ⋅ Maintenant, on s’intéresse à la grandeur physique aléatoire d’intérêt G(ξ) considérée précédemment. La méthode de Monte-Carlo permet de déterminer les différents moments de G(ξ) tels que son espérance, sa variance, etc… Comme cela a été montré ci-dessus, il suffit de calculer un nombre n suffisamment grand de réalisations de G : ( ) k G G k = ξ avec k =1 : n. Pour la mise en œuvre de la méthode de Monte-Carlo, on génère d’abord n réalisations indépendantes ξ k de ξ, k =1 : n. Pour chaque réalisation ξ k l’équation (1.63) est résolue: 33 k k k Λ Ω = ( ) ⋅ ( ) ( ) ξ ξ β ξ ′ (1.63) La réalisation Gk est obtenue par : ( k G G k = ( )) Ω′ ξ (1.64) L’intervalle de confiance de 95% de E( ( )) G ξ peut être obtenue par : [ 1.96 , 1.96 ] G G G G n n n n σ σ − ⋅ + ⋅ avec: 1 1 n n k k G G n = = ∑ et 2 σ G de la forme suivante : 2 2 1 1 ( ) 1 n G k n k G G n σ = = − − ∑ Cette méthode est non intrusive parce qu’il suffit de réaliser n calculs déterministes correspondant aux n réalisations des données d’entrée ξ k , k=1 : n. Cette méthode est facile à mettre en œuvre et on peut maîtriser la convergence de l’estimation. Par contre, la vitesse de convergence est assez lente (proportionnelle à 1/ n ). En effet, dans certain cas, on doit réaliser un grand nombre de calcul ce qui peut s’avérer numériquement lourd. On peut trouver dans la littérature la méthode de l’hypercube latin [25], la méthode « importance sampling » [38] permettant de réduire le nombre de calculs. Les méthodes d’échantillonnage précédentes permettent d’estimer certaines grandeurs statistiques liées à G(ξ). Néanmoins, il peut être utile d’avoir une expression explicite de G(ξ) en fonction de ξ. Dans la littérature, différentes bases de dimension finie sont utilisées pour approcher G(ξ) (le chaos polynomial, les ondelettes, etc. [39]). Dans la suite, nous allons nous concentrer sur le chaos polynomial sachant qu’une extension à d’autres bases est possible. La grandeur G(ξ) est approchée sous la forme : 1 ( ) ( ) ( ) P P i i i G G α = ξ ξ ξ ≈ = Ψ ∑ (1.65) avec ( ) Ψi ξ le chaos polynomial et αi des coefficients à déterminer. Nous allons présenter des méthodes qui permettent d’identifier les coefficients αi du développement (1.65). 

Méthode de régression

Dans la méthode de régression, les coefficients αi sont déterminés de la façon suivante [4] : 1 1 2 1 2 1 2 ( , ,… ) ( , ,… ) arg Min( ( ( ) ( )) ) P P P P G G α α α α α α ∈ = −   R   E ξ ξ (1.66) 34 En réalité, on ne dispose que de n évaluations de G(ξ) pour estimer 2 ( ( ) ( )) P   G G−   E ξ ξ . On obtient alors : 1 1 2 1 2 ( , ,… ) 1 2 ( , ,… ) arg Min ( ( , ,… )) P P P P r α α α α α α α α α ∈ =   R   (1.67) avec : 2 1 2 1 1 ( , ,… ) (G( ) ( )) n P k k P k i i k i r α α α ω α = = = − Ψ ∑ ∑ ξ ξ (1.68) En cherchant le point stationnaire de 1 2 ( , ,… ) P r α α α de (1.68), on obtient le système matriciel suivant: S ⋅ = α b (1.69) où : 1 1 ( ) ( ) ( )G( ) et avec , 1: S n ij k i k j k k n i k i k k k i i i j P ω ω α = = = Ψ Ψ = Ψ = = ∑ b ∑ ξ ξ ξ ξ α (1.70) La résolution de (1.69) permet d’obtenir les coefficients αi . Dans cette méthode, on doit effectuer n calculs G(ξ k ), k =1 : n. Rappelons que pour chaque ξ k , le problème (1.63) est résolu et l’évaluation G(ξ k ) est obtenue à partir de Ω’(ξ k ). Ces n calculs peuvent être réalisés en parallèle. Le choix des points de réalisation ξ k ainsi que celui des poids ωk est un point délicat. En effet, un choix non pertinent peut conduire à un système (1.69) singulier ou parfois mal conditionné. Dans [26], les points d’évaluation ξ k peuvent se construire à partir des racines du polynôme monodimensionnel ψi(ξ) et les poids peuvent être donnés par 1/ k ω = n . On peut constater aussi que pour un système (1.69) non singulier il faut prendre n P ≥ . Dans le cas où la dimension d est très grande et on utilise tous les polynômes d’ordre inférieur ou égal à p, la valeur de P (1.15) peut être très élevée. La réalisation au moins de P calculs peut être numériquement très lourde. Dans [15], une méthode est proposée consistant à choisir d’une façon adaptative des polynômes utilisés dans le cas d’une grande dimension d. Ce choix a pour but d’exclure les polynômes dont la contribution est faible dans l’approximation (1.65).