# Éléments finis stochastiques non intrusifs par chaos polynomiaux creux et adaptatifs

Géraud Blatman, Bruno Sudret

► To cite this version:

Géraud Blatman, Bruno Sudret. Éléments finis stochastiques non intrusifs par chaos polynomiaux creux et adaptatifs. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01413783

HAL Id: hal-01413783

https://hal.science/hal-01413783v1

Submitted on 11 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Éléments finis stochastiques non intrusifs par chaos
polynomiaux creux et adaptatifs

G. Blatman¹,², B. Sudret²,³

¹ EDF R&D, Département Matériaux et Mécanique des Composants
Site des Renardières, F-77818 Moret-sur-Loing

² LaMI, Institut Français de Mécanique Avancée et Université Blaise Pascal
Campus des Cézeaux, F-63175, Aubière Cedex France

³ Phimeca Engineering S.A.
Centre d'affaires du Zénith, 34 rue de Sarliève, 63800 Cournon

**Résumé** — On propose un algorithme permettant de construire une approximation par chaos polynomial (CP) creux de la réponse d'un modèle mécanique, dont les paramètres d'entrée sont aléatoires. On obtient au final une représentation par CP contenant peu de termes, pouvant être estimés par régression à partir d'un faible nombre d'évaluations du modèle, autrement dit à un coût de calcul réduit. La méthode est appliquée à l'analyse de fiabilité d'une ossature de bâtiment.

**Mots clés** — Analyse de fiabilité, chaos polynomial creux, procédure itérative.

## 1 Introduction

Les développements par chaos polynomial (CP) permettent de représenter de manière explicite la réponse d'un système mécanique dont les paramètres d'entrée sont des variables aléatoires. Les coefficients du CP peuvent être calculés efficacement au moyen de méthodes non intrusives telles que la régression des moindres carrés [1]. Cependant, le nombre d'évaluations du modèle (e.g. modèle aux éléments finis) - autrement dit le coût de calcul - augmente avec la taille du CP, qui s'accroît elle-même considérablement avec le nombre de variables d'entrée lorsque l'on applique le schéma usuel de troncature (i.e. on retient tous les polynômes dont le degré n'excède pas un *p* donné). Pour résoudre ce problème, on utilise une stratégie de troncature basée sur le principe de hiérarchie des effets, selon lequel la plupart des modèles dépendent principalement des effets principaux et des interactions simples. Les développements par CP correspondants contiennent un faible nombre de termes par rapport à une représentation pleine. De plus, une version anisotrope est développée pour réduire encore davantage la taille du métamodèle. Elle repose sur une sélection des variables aléatoires d'entrée associées à des indices de sensibilité importants [2].

À partir de ces stratégies de troncature, on propose une procédure iterative visant à détecter progressivement un faible nombre de termes importants dans la décomposition par CP, ce qui conduit à une représentation creuse. Par ailleurs, le plan d'expériences est systématiquement en-richi de sorte à ce que le problème de régression soit bien posé [3]. On utilise dans ce but des techniques d'échantillonnage séquentielles, telles que les suites de nombres quasi-aléatoires [4].
## 2 Développement par chaos polynomial de la réponse d'un modèle

On considère un système mécanique décrit par un modèle numérique $\mathcal{M}$ qui peut être analytique ou plus généralement algorithmique (e.g. modèle aux éléments finis). On suppose que ce modèle dépend de $\mathcal{M}$ paramètres incertains que l'on représente par des variables aléatoires indépendantes $\{X_1, \dots, X_M\}^T$, regroupées dans un vecteur aléatoire $X$. On suppose connue la densité jointe $f_X(x)$ de $X$. La réponse du modèle notée $Y = \mathcal{M}(X)$ est alors elle aussi aléatoire. On considère dans ce travail une réponse scalaire $Y$ de variance finie. On montre que $Y$ peut s'exprimer comme suit [5] :

$$Y = \mathcal{M}(X) = \sum_{\alpha \in \mathbb{N}^M} a_{\alpha} \psi_{\alpha}(X) \quad (1)$$

Ce développement est appelé représentation par chaos polynomial (CP) de $Y$. Les $a_\alpha$ sont des coefficients déterministes inconnus et les $\psi_\alpha$ sont des polynômes multivariés qui sont orthonormaux par rapport à la densité $f_X(x)$, i.e. $\mathbb{E}[\psi_\alpha(X)\psi_\beta(X)] = 1$ si $\alpha = \beta$ et 0 sinon.

On considère un ensemble de réalisations de $X$ noté $X = \{x^{(1)}, \dots, x^{(N)}\}$ et appelé plan d'expériences (PE). On note $\mathcal{Y} = \{\mathcal{M}(x^{(1)}), \dots, \mathcal{M}(x^{(N)})\}$ l'ensemble des évaluations du modèle correspondantes. Les coefficients du CP peuvent être estimés par régression des moindres carrés [1]. Cette méthode requiert le choix a priori d'une troncature de CP, autrement dit d'un ensemble fini non vide $\mathcal{A} = \{\alpha_0, \dots, \alpha_{p-1}\} \subset \mathbb{N}^M$, dit ensemble de troncature, qui contient les multi-indices de la base polynomiale retenue $\{\psi_{\alpha_0}, \dots, \psi_{\alpha_{p-1}}\}$. L'approximation par CP associée est notée $Y_{\mathcal{A}} \equiv \mathcal{M}_{\mathcal{A}}(X)$. En pratique, on tronque souvent le développement par CP en tenant les polynômes de degré inférieur ou égal à $p$, ce qui conduit à l'ensemble de troncature $\mathcal{A}^{M,p} = \{\alpha \in \mathbb{N}^M : \|\alpha\|_1 \equiv \sum_{i=1}^M \alpha_i \le p\}$, de cardinalité $P = \binom{M+p}{p}$. Le nombre de termes du CP augmente donc fortement avec $p$ et $M$. Par conséquent, la taille minimale du PE (i.e. le coût de calcul) requise pour obtenir une solution précise du problème de régression s'accroît considérablement.

## 3 Stratégie de troncature du chaos polynomial basée sur le principe de hiérarchie des effets

### 3.1 Schéma isotrope

On propose une stratégie de troncature du CP basée sur le principe de hiérarchie des effets [6], selon lequel la plupart des modèles dépendent principalement des effets principaux et des interactions simples. On définit ainsi des ensembles de troncature à partir de normes $q$, $0 < q < 1$ :

$$\mathcal{A}_{q}^{M,p} = \left\{ \alpha \in \mathbb{N}^{M} : \|\alpha\|_{q} = \left( \sum_{i=1}^{M} \alpha_{i}^{q} \right)^{1/q} \le p \right\} \quad (2)$$

Les normes $q$ pénalisent d'autant plus les interactions d'ordre élevé que $q$ est faible. On note que le choix $q = 1$ correspond au schéma de troncature usuel présenté dans la section précédente.

La stratégie proposée conduit ainsi à un développement par CP qui contient un faible nombre de coefficients inconnus, qui peuvent donc être calculés à partir d'un nombre $N$ réduit d'évaluations du modèle. Il est possible de réduire davantage le coût de calcul en prenant en compte le fait que toutes les variables d'entrée du modèle n'ont pas le même impact sur la réponse, comme le montrent les analyses de sensibilité [2].
## 3.2 Schéma anisotrope

Étant donné un développement par CP, on définit tout d'abord les indices de sensibilité totaux de la réponse du modèle à chacune des variables d'entrée $X_i$ par :

$$S_i^T = \sum_{\alpha \in I_i} a_{\alpha}^2 \qquad (3)$$

où $I_i$ est l'ensemble de tous les multi-indices de $\mathbb{N}^M$ dont la i-ème composante est non nulle, c'est-à-dire $I_i = \{\alpha \in \mathbb{N}^M : \alpha_i > 0\}$. Les quantités $S_i^T$ représentent approximativement la part de la variance de la réponse qui est expliquée par la variance de la variable $X_i$ (voir [2] pour plus de détails).

On cherche ici à tronquer le CP tout en favorisant les variables d'entrée qui correspondent à des $S_i^T$ importants. On définit dans ce but un ensemble de poids $w = \{w_1, \dots, w_M\}$ par :

$$w_i = \frac{\sum_{j=1}^{M} S_j^T + \max_{1 \le k \le M} S_k^T - S_i^T}{\sum_{j=1}^{M} S_j^T}, \quad i = 1, \dots, M \qquad (4)$$

On considère alors les normes anisotropes et les ensembles de troncature suivants :

$$\|z\|_{q,w} = \left( \sum_{i=1}^{M} |w_i z_i|^q \right)^{1/q}, \quad \mathcal{A}_{q,w}^{M,p} = \{ \alpha \in \mathbb{N}^M : \|\alpha\|_{q,w} \le p \} \qquad (5)$$

Le schéma de troncature anisotrope proposé peut ainsi conduire à une réduction notable du nombre de termes du CP, et donc à une forte diminution du nombre d'évaluations du modèle.

# 4 Chaos polynomial creux adaptatif

À partir de la stratégie de troncature décrite dans la section précédente, on définit un algorithme permettant une construction pas-à-pas d'un CP creux, en détectant automatiquement les termes significatifs [3]. Une description synthétique de la procédure est présentée en Figure 1. On sélectionne un PE initial $\mathcal{X}$ basé sur des suites de nombres quasi-aléatoires [4], et l'on rassemble les évaluations du modèle correspondantes dans un vecteur $\mathcal{Y}$. Si l'on choisit la stratégie de troncature anisotrope, on initialise les poids comme suit : $w = \{1, \dots, 1\}$. On approxime tout d'abord la réponse du modèle par un CP de degré $p=0$, i.e. $\mathcal{A} = \{0\}$, et l'on calcule des estimateurs de précision $R^2$ et $Q^2$. $R^2$ est le coefficient de détermination couramment utilisé en régression, et $Q^2$ est un estimateur basé sur une technique de validation croisée qui est relativement robuste au phénomène de surapprentissage (voir [3] pour plus de détails). Ensuite, on enrichit la base du CP en ajoutant des termes candidats de la base, de norme $(q,w)$ croissante. On ne retient au final que les termes qui conduisent à une augmentation significative du $R^2$, i.e. plus qu'une valeur-seuil $\varepsilon_1$ (étape ascendante). Par ailleurs, on retire tous les termes dont l'élimination de la base courante ne diminue pas significativement le $R^2$, i.e. moins qu'une valeur-seuil $\varepsilon_2$ (étape descendante). Si l'on choisit la stratégie de troncature anisotrope, on calcule les indices de sensibilité (3) et on actualise les poids (4).

Lorsque l'on enrichit la base du CP (i.e. l'ensemble de troncature $\mathcal{A}$), le nombre $P$ de termes retenus peut s'approcher de la taille $N$ du PE, d'où un mauvais conditionnement du problème de régression à résoudre. Pour éviter ce problème, on ajoute des points au PE de sorte à satisfaire la règle heuristique $N \ge 2P$. Dans ce but, il est nécessaire d'adopter une stratégie séquentielle d'échantillonnage. Lorsque $N > 2P$, l'ensemble de troncature $\mathcal{A}$ est réinitialisé au singleton $\{0\}$ et l'on redemarre la procédure d'enrichissement de la base du CP. On stoppe l'algorithme dès qu'une précision-cible $Q_{cible}^2$ est atteinte.
<figure><img src="image_2.png" /><figcaption>Figure 1 – Procédure itérative pour construire un chaos polynomial creux</figcaption></figure>

# 5 Application à l'analyse de fiabilité d'une ossature de bâtiment

## 5.1 Présentation du problème

On considère le modèle d'ossature de bâtiment illustré en Figure 2. La réponse d'intérêt est le déplacement horizontal $u$ au noeud situé en haut à droite de la structure.

<figure><img src="image_3.png" /><figcaption>Figure 2 – Problème de l'ossature de bâtiment</figcaption></figure>

<table><thead><tr><th>Element</th><th>Young's modulus</th><th>Moment of inertia</th><th>Cross-sect. area</th></tr></thead><tbody><tr><td>B<sub>1</sub></td><td>E<sub>4</sub></td><td>I<sub>10</sub></td><td>A<sub>18</sub></td></tr><tr><td>B<sub>2</sub></td><td>E<sub>4</sub></td><td>I<sub>11</sub></td><td>A<sub>19</sub></td></tr><tr><td>B<sub>3</sub></td><td>E<sub>4</sub></td><td>I<sub>12</sub></td><td>A<sub>20</sub></td></tr><tr><td>B<sub>4</sub></td><td>E<sub>4</sub></td><td>I<sub>13</sub></td><td>A<sub>21</sub></td></tr><tr><td>C<sub>1</sub></td><td>E<sub>5</sub></td><td>I<sub>6</sub></td><td>A<sub>14</sub></td></tr><tr><td>C<sub>2</sub></td><td>E<sub>5</sub></td><td>I<sub>7</sub></td><td>A<sub>15</sub></td></tr><tr><td>C<sub>3</sub></td><td>E<sub>5</sub></td><td>I<sub>8</sub></td><td>A<sub>16</sub></td></tr><tr><td>C<sub>4</sub></td><td>E<sub>5</sub></td><td>I<sub>9</sub></td><td>A<sub>17</sub></td></tr></tbody></table>

Les 3 charges appliquées, les 2 modules d'Young, les 8 moments d'inertie et les 8 sections droites des poutres constituant la structure sont supposées aléatoires. Ces variables sont réunies dans le vecteur aléatoire $X = (P_1, P_2, P_3, ..., I_6, ..., I_{13}, A_{14}, ..., A_{21})$ de taille $M = 21$. Les propriétés des variables aléatoires sont indiquées dans le Tableau 5.1.

De plus les variables aléatoires sont corrélées. Soit $F_i(x_i)$ la fonction de répartition marginale de la variable $X_i$. On pose $\hat{\xi} = \{\hat{\xi}_i = \Phi^{-1}(F_i(X_i)), i = 1, ..., M\}$. Les $\hat{\xi}_i$ sont des variables aléatoires gaussiennes standard, corrélées comme suit :
<table><thead><tr><th>Variable</th><th>Distribution</th><th>Moyenne</th><th>Écart-type</th></tr></thead><tbody><tr><td>P<sub>1</sub> (kN)</td><td>Lognormale</td><td>133,454</td><td>40,04</td></tr><tr><td>P<sub>2</sub> (kN)</td><td>"</td><td>88,97</td><td>35,59</td></tr><tr><td>P<sub>3</sub> (kN)</td><td>"</td><td>71,175</td><td>28,47</td></tr><tr><td>E<sub>4</sub> (kN/m<sup>2</sup>)</td><td>Gaussienne tronquée sur [0, +∞[</td><td>2,1738×10<sup>-7</sup></td><td>1,9152×10<sup>-6</sup></td></tr><tr><td>E<sub>5</sub> (kN/m<sup>2</sup>)</td><td>"</td><td>2,3796×10<sup>-7</sup></td><td>1,9152×10<sup>-6</sup></td></tr><tr><td>I<sub>6</sub> (m<sup>4</sup>)</td><td></td><td>8,1344×10<sup>-3</sup></td><td>1,0834×10<sup>-3</sup></td></tr><tr><td>I<sub>7</sub> (m<sup>4</sup>)</td><td></td><td>1,1509×10<sup>-2</sup></td><td>1,2980×10<sup>-3</sup></td></tr><tr><td>I<sub>8</sub> (m<sup>4</sup>)</td><td></td><td>2,1375×10<sup>-2</sup></td><td>2,5961×10<sup>-3</sup></td></tr><tr><td>I<sub>9</sub> (m<sup>4</sup>)</td><td></td><td>2,5961×10<sup>-2</sup></td><td>3,0288×10<sup>-3</sup></td></tr><tr><td>I<sub>10</sub>(m<sup>4</sup>)</td><td>"</td><td>1,0812×10<sup>-2</sup></td><td>2,5961×10<sup>-3</sup></td></tr><tr><td>I<sub>11</sub>(m<sup>4</sup>)</td><td></td><td>1,4105×10<sup>-2</sup></td><td>3,4615×10<sup>-3</sup></td></tr><tr><td>I<sub>12</sub>(m<sup>4</sup>)</td><td></td><td>2,3279×10<sup>-2</sup></td><td>5,6249×10<sup>-3</sup></td></tr><tr><td>I<sub>13</sub>(m<sup>4</sup>)</td><td></td><td>2,5961×10<sup>-2</sup></td><td>6,4902×10<sup>-3</sup></td></tr><tr><td>A<sub>14</sub>(m<sup>2</sup>)</td><td></td><td>3,1256×10<sup>-1</sup></td><td>5,5815×10<sup>-2</sup></td></tr><tr><td>A<sub>15</sub>(m<sup>2</sup>)</td><td></td><td>3,7210×10<sup>-1</sup></td><td>7,4420×10<sup>-2</sup></td></tr><tr><td>A<sub>16</sub>(m<sup>2</sup>)</td><td></td><td>5,0606×10<sup>-1</sup></td><td>9,3025×10<sup>-2</sup></td></tr><tr><td>A<sub>17</sub>(m<sup>2</sup>)</td><td></td><td>5,5815×10<sup>-1</sup></td><td>1,1163×10<sup>-1</sup></td></tr><tr><td>A<sub>18</sub>(m<sup>2</sup>)</td><td>"</td><td>2,5302×10<sup>-1</sup></td><td>9,3025×10<sup>-2</sup></td></tr><tr><td>A<sub>19</sub>(m<sup>2</sup>)</td><td></td><td>2,9117×10<sup>-1</sup></td><td>1,0232×10<sup>-1</sup></td></tr><tr><td>A<sub>20</sub>(m<sup>2</sup>)</td><td></td><td>3,7303×10<sup>-1</sup></td><td>1,2093×10<sup>-1</sup></td></tr><tr><td>A<sub>21</sub>(m<sup>2</sup>)</td><td></td><td>4,1860×10<sup>-1</sup></td><td>1,9537×10<sup>-1</sup></td></tr></tbody></table>

Les moyennes et les écart-types des sections droites, moments d'inertie et modules d'Young sont ceux des distributions gaussiennes originelles avant la troncature

<caption>Tableau 1 – Propriétés des variables aléatoires</caption>

- les $\tilde{\xi}_i$ associées à la section droite et au moment d'inertie d'un même élément sont corrélées par $\rho = 0,95$;

- les autres $\tilde{\xi}_i$ correspondant aux propriétés géométriques sont corrélées par $\rho_{A_i,A_j} = \rho_{I_i,I_j} = \rho_{A_i,A_j} = 0,13$;

- les $\tilde{\xi}_i$ reliées aux 2 modules d'Young sont corrélées par $\rho = 0,9$;

La réponse du modèle peut finalement s'exprimer comme une fonction de variables gaussiennes standard indépendantes $\xi_i$ au moyen d'une transformation de Nataf [7].

## 5.2 Analyse de fiabilité

On étudie la fiabilité de l'ossature vis-à-vis de la fonction d'état-limite suivante :

$$g(X) = u_{max} - \mathcal{M}(X) \quad , \quad u_{max} = 6,1 \text{ cm} \qquad (6)$$

La probabilité de défaillance associée s'écrit $P_f = \mathbb{P}[g(X) < 0]$. On estime cette quantité en remplaçant la réponse par des approximations par différents types de CP sur des bases de polynômes d'Hermite normalisés, d'où une fonction d'état-limite *analytique*. Cette dernière est post-traitée en utilisant la méthode FORM suivie d'un tirage d'importance (50 000 tirages sont effectués). On utilise des PE basés sur des suites quasi-aléatoires de Sobol' [4] pour estimer les coefficients du CP par régression. Pour une juste comparaison, on utilise des PE de taille $N = 2\text{card}(\mathcal{A})$ dans le cas du CP plein. Les valeurs-seuils $\varepsilon_1$ et $\varepsilon_2$ sont fixées respectivement à 0,005(1 - $Q_{cible}^2$) et 0,001(1 - $Q_{cible}^2$). En effet, ces valeurs donnent de bons résultats d'après l'expérience des auteurs. À titre de validation, on estime $P_f$ par simulation de Monte Carlo directe du problème (54 438 évaluations du modèle sont effectuées de sorte à avoir un coefficient de variation inférieur à 1%). Les résultats en termes d'indices de fiabilité généralisés $\beta = \Phi^{-1}(P_f)$ sont donnés dans le Tableau 5.2.

Il apparaît que toutes les approximations par CP donnent des estimateurs de $\beta$ avec des erreurs relatives inférieures à 2% par rapport au résultat de référence. Le CP creux isotrope associé à la
<table><thead><tr><th rowspan="2"></th><th rowspan="2">Référence</th><th rowspan="2">CP plein<br/>p = 3</th><th colspan="3">CP creux - Q<sup>2</sup><sub>cible</sub> = 0,999</th></tr><tr><th>ISO q = 1</th><th>ISO q = 0,75</th><th>ANIS q = 0,75</th></tr></thead><tbody><tr><td>Indice de fiabilité</td><td>3,56</td><td>3,66</td><td>3,66</td><td>3,64</td><td>3,66</td></tr><tr><td># évaluations du modèle</td><td>54 438</td><td>4 048</td><td>1 500</td><td>900</td><td>500</td></tr><tr><td># termes du CP</td><td>-</td><td>2 024</td><td>213</td><td>149</td><td>144</td></tr></tbody></table>

ISO et ANIS signifient respectivement *isotrope* and *anisotrope*

<caption>Tableau 2 – Estimations de l'indice de fiabilité généralisé par les différentes approximations par chaos polynomial</caption>

norme $q = 1$ permet un gain en coût de calcul d'un facteur 3 par rapport au CP plein. Le coût est encore réduit en utilisant une norme $q = 0,75$, avec un gain de facteur 4,5. Le CP anisotrope utilisé avec la norme $q = 0,75$ se révèle le plus efficace, avec un gain de facteur 8. On note que la procédure itérative proposée fournit des CP très creux, avec un nombre relatif de termes du CP de 10% (resp. 7%) par rapport au CP plein quand on utilise la norme $q = 1$ (resp. $q = 0,75$).

## 6 Conclusion

On propose une méthode pour construire de manière itérative un développement creux par CP de la réponse d'un modèle mécanique dont les paramètres sont aléatoires. Elle repose sur un schéma original de troncature basé sur l'utilisation de normes $q$ ($0 < q < 1$) isotropes ou anisotropes. Un algorithme est présenté pour réduire le nombre de coefficients inconnus du CP, ce qui conduit à une représentation creuse. Par ailleurs, le PE est systématiquement enrichi de sorte à ce que les problèmes de régression soient bien posés. Les coefficients retenus du CP peuvent finalement être efficacement évalués par régression au moyen d'un faible nombre d'évaluations du modèle, comparé à celui qui serait requis pour calculer un développement plein par CP. L'exemple de l'ossature montre que la stratégie proposée peut être appliquée à des analyses de fiabilité, conduisant à un gain d'un facteur 8 en coût de calcul par rapport aux CP usuels. La sensibilité de la procédure à la valeur $q$ de la norme sera étudiée dans de prochains travaux.

## Références

[1] M. Berveiller, B. Sudret, and M. Lemaire. Stochastic finite elements : a non intrusive approach by regression. *Eur. J. Comput. Mech.*, 15(1-3):81–92, 2006.

[2] B. Sudret. Global sensitivity analysis using polynomial chaos expansions. *Reliab. Eng. Sys. Safety*, 93:964–979, 2008.

[3] G. Blatman and B. Sudret. Use of sparse polynomial chaos expansions in adaptive stochastic finite element analysis. *Prob. Eng. Mech.*, 2008. accepted for publication.

[4] H. Niederreiter. *Random number generation and quasi-Monte Carlo methods*. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 1992.

[5] C. Soize and R. Ghanem. Physical systems with random uncertainties : chaos representations with arbitrary probability measure. *SIAM J. Sci. Comput.*, 26(2):395–410, 2004.

[6] D.C. Montgomery. *Design and analysis of experiments*. 2004.

[7] A. Nataf. Détermination des distributions dont les marges sont données. *C. R. Acad. Sci. Paris*, 225:42–43, 1962.