# Optimisation sous contrainte de fiabilité d'une coque
imparfaite

Vincent Dubourg, Bruno Sudret, Jean-Marc Bourinet, Mikaël Cazuguel

► To cite this version:

Vincent Dubourg, Bruno Sudret, Jean-Marc Bourinet, Mikaël Cazuguel. Optimisation sous contrainte de fiabilité d'une coque imparfaite. 10e colloque national en calcul des structures, May 2011, Giens, France. hal-00592795

HAL Id: hal-00592795

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

Submitted on 3 May 2011

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
# Optimisation sous contrainte de fiabilité d'une structure en treillis

V. Dubourg<fnref n="1" />,<fnref n="2" />, B. Sudret<fnref n="1" />,<fnref n="2" />, J.-M. Bourinet<fnref n="1" />, M. Cazuguel<fnref n="3" />

**Résumé** — On s'intéresse dans ce papier au problème d'optimisation sous contrainte de fiabilité (RBDO, pour reliability-based design optimization) qui vise à minimiser le coût d'une structure en modifiant la valeur de ses paramètres de conception (e.g. dimensions caractéristiques) tout en garantissant un niveau de fiabilité. Celui-ci est évalué en prenant en compte les incertitudes inhérentes au processus de fabrication (e.g. tolérances sur les dimensions, variabilité des propriétés matériaux) ou à l'environnement dans lequel sera placé la structure (chargements appliqués). Dans l'optique de résoudre des problèmes industriels faisant intervenir des modèles numériques fins et coûteux, on propose ici une méthode de résolution qui s'appuie sur la construction d'un méta-modèle de la fonction de performance de la structure. On utilise pour ce faire le krigeage. On présente la construction adaptative de ce méta-modèle et son intégration dans l'algorithme d'optimisation global, dans lequel la fiabilité de la structure et sa sensibilité par rapport aux paramètres de conception sont évalués à chaque étape par la méthode de subset simulation. L'ensemble de la démarche est appliquée au dimensionnement optimal d'un treillis isostatique.

**Mots clés** — optimisation sous contrainte de fiabilité, krigeage, subset simulation, espace de fiabilité augmenté

## 1 Introduction

L'optimisation sous contrainte de fiabilité cherche à déterminer une conception à coût minimal tout en assurant une performance avec un haut degré de fiabilité. Précisément, on vise à minimiser une fonction objectif $c$ en agissant sur des variables de dimensionnement caractéristiques $\theta$ (e.g. des valeurs moyennes) sous des contraintes de fiabilité venant remplacer ou compléter les contraintes déterministes usuelles. Le dimensionnement optimal d'intérêt $\theta^*$ est donc défini par :

$$ \theta^* = \arg \min_{\theta \in \mathcal{D}_{\theta}} c(\theta) : \begin{cases} f_i(\theta) \le 0, & i = 1, \dots, n_c \\ \mathbb{P}(g_l(\mathbf{X}(\theta)) \le 0) \le P_{fl}^0, & l = 1, \dots, n_p \end{cases} \quad (1) $$

où :

- $\mathcal{D}_{\theta}$ est l'espace de conception admissible,

- $\{f_i, i = 1, \dots, n_c\}$ sont des contraintes déterministes ne mettant pas en jeu l'intégrité du système,

- $\{\mathbb{P}(g_l(\mathbf{X}(\theta)) \le 0) \le P_{fl}^0, l = 1, \dots, n_p\}$ sont les contraintes dites de fiabilité faisant intervenir le vecteur aléatoire $\mathbf{X}$ modélisant les incertitudes inhérentes au système à concevoir.

Dans ces dernières équations, $\{P_{fl}^0, l = 1, \dots, n_p\}$ sont les probabilités de défaillance maximales tolérées pour chaque mode de défaillance. La notation $\mathbf{X}(\theta)$ signifie que les paramètres de conception sur lesquels va porter l'optimisation interviennent dans la définition du vecteur aléatoire des paramètres d'entrée du problème : les $\theta_i$ seront en général les valeurs moyennes de variables aléatoires $X_i$, dont la variabilité (supposée non nulle) est liée à la qualité du process de fabrication. Le vecteur $\mathbf{X}$ regroupe par ailleurs les variables d'environnement (e.g. propriétés matériaux, chargements) dont l'incertitude va intervenir dans le calcul de fiabilité.

Une première approche de résolution du problème d'optimisation formulé à l'équation 1 consisterait à utiliser un algorithme d'optimisation sous contrainte dans laquelle on viendrait évaluer, à chaque itération, les contraintes probabilistes par des simulations de Monte-Carlo. Cette approche dite brute-force n'est pas applicable en pratique car les fonctions de performances $\{g_l, l = 1, \dots, n_p\}$ font intervenir des variables issues de simulations numériques souvent coûteuses (e.g. ici, la flèche maximale d'un treillis

<footnote n="1">¹ Clermont Université, IFMA, EA 3867, Laboratoire de Mécanique et Ingénieries, BP 10448, F-63000 Clermont-Ferrand</footnote>

<footnote n="2">² Phimeca Engineering, Centre d’Affaires du Zénith, 34 rue de Sarliève, F-63800 Cournon d’Auvergne</footnote>

<footnote n="3">³ DCNS, Ingénierie Navires Armés, rue Choiseul, F-56311 LORIENT Cedex</footnote>
évalué par un calcul aux éléments finis). C'est pourquoi des stratégies alternatives ont été proposées dans la littérature RBDO (Reliability-Based Design Optimization), voir [1] pour une revue des méthodes existantes.

Pour évaluer le plus efficacement possible le niveau de fiabilité à chaque itération de l'algorithme d'optimisation, la méthode FORM (First Order Reliability Method) a été utilisée [2]. Celle-ci n'est cependant pas applicable aux cas industriels visés qui peuvent présenter des états-limites fortement non-linéaires et/ou des configurations défaillantes les plus probables multiples.

Ainsi, en vue d'étendre le champ d'application de l'optimisation sous contrainte de fiabilité à des problèmes plus complexes, une nouvelle stratégie de résolution a été proposée. Elle est basée sur l'utilisation :

- de méta-modèles par krigeage [3] substitués aux fonctions de performance originales afin d'en réduire le coût d'évaluation de façon significative ;

- d'une technique d'estimation de probabilités faibles plus efficace que les simulations de Monte-Carlo connue sous le nom de subset simulation [4].

# 2 Krigeage et apprentissage actif pour la représentation des fonctions de performance

## 2.1 Principe du krigeage

Dans le domaine de l'expérimentation numérique, on appelle méta-modèle une fonction peu coûteuse à évaluer qui reproduit le comportement d'un modèle mathématique $\mathcal{M}: \mathcal{D}_x \subseteq \mathbb{R}^n \to \mathcal{D}_y \subseteq \mathbb{R}$ coûteux en temps de calcul, comme par exemple un modèle aux éléments finis. Pour construire le méta-modèle, on définit un plan d'expérience numérique, c'est-à-dire un ensemble de points $X = \{x_i \in \mathcal{D}_x, i = 1, \dots, m\}$ sur lequel on évalue le modèle, soit $Y = \{y_i = \mathcal{M}(x_i), i = 1, \dots, m\}$.

La méthode du krigeage [3] est une technique particulière de construction de méta-modèles qui s'appuie sur l'hypothèse heuristique suivante : la réponse $y = \mathcal{M}(x)$ est supposée être une réalisation particulière d'un processus gaussien $Y(x)$ défini sur $\mathcal{D}_x$ comme suit :

$$Y(x) = f(x)^{\top} \beta + Z(x) \quad (2)$$

Dans cette équation, on a noté :

- $f(x)^{\top} \beta$ la moyenne du processus, qui est représentée sur une base de fonctions choisies $\{f_i, i = 1, \dots, p\} \in L_2(\mathcal{D}_x, \mathbb{R})$ (en pratique, un polynôme complet d'ordre donné) ;

- $Z(x)$ un processus gaussien stationnaire de moyenne nulle, de variance $\sigma_Y^2$, dont on définit la fonction d'auto-corrélation $(\ell)$ rassemble les paramètres de la fonction) :

$$C_{YY}(x, x') = \sigma_Y^2 R(|x - x'|, \ell), \quad (x, x') \in \mathcal{D}_x \times \mathcal{D}_x \quad (3)$$

En pratique, on utilise le plus souvent un modèle d'auto-corrélation de type "exponentiel carré", qui s'écrit sous la forme :

$$R(|x - x'|, \ell) = \exp\left(\sum_{k=1}^{n} - \left|\frac{x_k - x'_k}{\ell_k}\right|^2\right) \quad (4)$$

L'estimation par krigeage suppose qu'au point courant $x$ le méta-modèle est la réalisation d'une variable gaussienne notée $\hat{Y}(x) \sim \mathcal{N}(\mu_{\hat{Y}}(x), \sigma_{\hat{Y}}(x))$. En cherchant le meilleur estimateur linéaire sans biais (BLUE en anglais, pour best linear unbiased estimator [3]), on arrive tous calculs faits à l'expression suivante pour $\mu_{\hat{Y}}(x)$ :

$$\mu_{\hat{Y}}(x) = f(x)^{\top} \hat{\beta} + r(x)^{\top} \mathbf{R}^{-1} (\mathbf{Y} - \mathbf{F}\hat{\beta}) \quad (5)$$

Dans cette équation on a introduit les notations classiques $r$, $\mathbf{R}$ et $\mathbf{F}$ suivantes :

$$r_i(x) = R(|x - x_i|, \ell), \quad i = 1, \dots, m \quad (6)$$

$$R_{ij} = R(|x_i - x_j|, \ell), \quad i = 1, \dots, m, \quad j = 1, \dots, m \quad (7)$$

$$F_{ij} = f_i(x_j), \quad i = 1, \dots, p, \quad j = 1, \dots, m \quad (8)$$
En plus de la “prédiction moyenne” donnée à l’éq. (5), le krigeage permet d’estimer, à travers la variance de krigeage $\sigma_{\hat{Y}}(x)$, une incertitude épistémique de prédiction, liée au nombre fini de points utilisés dans le plan d’expériences. Celle-ci vaut :

$$ \sigma_{\hat{Y}}^2(x) = \sigma_Y^2 \left( 1 - \langle \mathbf{f}(x)^\top \mathbf{r}(x)^\top \rangle \begin{bmatrix} \mathbf{0} & \mathbf{F}^\top \\ \mathbf{F} & \mathbf{R} \end{bmatrix}^{-1} \begin{bmatrix} \mathbf{f}(x) \\ \mathbf{r}(x) \end{bmatrix} \right) \quad (9) $$

On peut montrer que le méta-modèle donné par l’éq. (5) est interpolant, c’est-à-dire $\mu_{\hat{Y}}(x_i) = \mathcal{M}(x_i)$ d’une part, et $\sigma_{\hat{Y}}(x_i) = 0$ d’autre part.

Jusqu’ici on n’a pas précisé comment calculer à partir du plan d’expériences les différents paramètres du modèle de krigeage, à savoir les coefficients $\beta$ de la partie “régression” qui décrivent la tendance moyenne, la variance du processus $\sigma_Y^2$ ainsi que les longueurs de corrélation de la fonction d’auto-corrélation regroupés dans $\ell = \{\ell_k, k = 1, \dots, n\}$.

Cette opération est réalisée en appliquant le principe du maximum de vraisemblance bien connu en inférence statistique : les évaluations du modèle disponibles, regroupées dans $\mathbf{Y}$, correspondent à une unique réalisation d’un vecteur gaussien dont les propriétés dépendent du jeu de paramètres $\{\beta, \sigma_Y^2, \ell\}$ à déterminer. La maximisation de la fonction de vraisemblance s’effectue en deux temps : pour $\ell$ fixé (ce qui définit complètement la matrice $\mathbf{R}$), on a une expression analytique des coefficients de la régression :

$$ \beta = (\mathbf{F}^\top \mathbf{R}^{-1} \mathbf{F})^{-1} \mathbf{F}^\top \mathbf{R}^{-1} \mathbf{Y} \quad (10) $$

ainsi que de la variance du processus $Z(x)$ :

$$ \sigma_Y^2 = \frac{1}{m} (\mathbf{Y} - \mathbf{F}^\top \beta)^\top \mathbf{R}^{-1} (\mathbf{Y} - \mathbf{F}^\top \beta) \quad (11) $$

Il reste ensuite à optimiser une fonction qui ne dépend plus que des paramètres de la fonction d’auto-corrélation, c’est-à-dire $\ell$, que l’on reporte ensuite dans les équations (10) et (11).

## 2.2 Apprentissage actif

Dans la section précédente, on a montré comment la méthode du krigeage permet de construire un méta-modèle $\hat{Y}(x)$ à partir d’un plan d’expériences $\mathcal{X}$ sur lequel est évalué le modèle mathématique $\mathcal{M}$ d’intérêt. La variance de krigeage donnée par l’éq. (9) permet de quantifier en chaque point la “précision” du modèle ainsi construit. Elle peut être utilisée dans un processus adaptatif dans lequel on cherche à raffiner itérativement le modèle de krigeage en enrichissant de façon intelligente le plan d’expériences $\mathcal{X}$.

Pour l’application qui nous intéresse ici, à savoir représenter par krigeage chaque surface d’état-limite $S_l = \{g_l = 0, l = 1, \dots, n_p\}$ dans le problème RBDO, on souhaite avoir une grande précision dans la prédiction pour les points de l’espace $\mathcal{D}_x$ qui conduisent à des valeurs nulles des fonctions $g_l$. Dans un processus d’apprentissage adaptatif, il faut donc être capable de déterminer quelles zones de l’espace $\mathcal{D}_x$ correspondent aux valeurs proches des zéros de $g$, et d’enrichir le plan d’expériences dans ces zones.

On définit pour ce faire la marge de prédiction suivante (on se limite à une seule fonction $g$ dans la suite de l’exposé par souci de clarté) :

$$ \mathcal{M} = \{ x : -k \sigma_{\mu_{\hat{G}}} (x) \le \mu_{\hat{G}} (x) \le +k \sigma_{\mu_{\hat{G}}} (x) \} \quad (12) $$

où $k$ est en pratique pris égal à $k = \Phi^{-1}(97.5\%) = 1.96$. Ce choix signifie que du point de vue de l’incertitude épistémique de prédiction du krigeage, on retient un intervalle de confiance à 95% autour de la moyenne. La probabilité qu’un point $x$ (par rapport à l’incertitude épistémique liée au krigeage) appartienne à cette marge de prédiction s’écrit de façon analytique, puisque par construction la variable $\hat{G}(x)$ est gaussienne :

$$ \mathbb{P}(x \in \mathcal{M}) = \Phi \left( \frac{k \sigma_{\mu_{\hat{G}}}(x) - \mu_{\hat{G}}(x)}{\sigma_{\mu_{\hat{G}}}(x)} \right) - \Phi \left( \frac{-k \sigma_{\mu_{\hat{G}}}(x) - \mu_{\hat{G}}(x)}{\sigma_{\mu_{\hat{G}}}(x)} \right) \quad (13) $$
où $\Phi(\bullet)$ désigne la fonction de répartition de la loi normale centrée réduite. Pour raffiner le plan d'expérience dans la région proche de la surface d'état-limite réelle définie par $S = \{x : g(x) = 0\}$, on doit chercher des points maximisant la probabilité définie à l'éq. (13). Ce concept est largement répandu dans la littérature sur les plans d'expérience adaptatifs appliqués aux prédictions par krigeage – voir e.g. [5, 6]. Néanmoins, la plupart des stratégies proposées cherchent le point $\boldsymbol{x}$ maximisant un critère d'enrichissement tel que celui défini à l'éq. (13). Or, ces critères sont généralement hautement multi-modaux ce qui signifie que les candidats à l'enrichissement ne sont pas uniques et il arrive même qu'il en existe une infinité le long du contour identifié.

C'est pourquoi on préfère appliquer ici la méthode d'enrichissement par échantillonnage présentée dans [7]. On utilise le critère défini précédemment comme une densité de probabilité à une constante de normalisation près que l'on vient pondérer par une densité de probabilité $w$ :

$$C(\boldsymbol{x}) \propto \mathbb{P}(\boldsymbol{x} \in \mathbb{M}) w(\boldsymbol{x}) \quad (14)$$

Un bon choix pour cette densité de pondération $w$ est la densité de probabilité $f_X$ du vecteur aléatoire $X$, ce qui permet de concentrer l'enrichissement dans les régions les plus probables pour les paramètres incertains. Si celle-ci n'est pas disponible, on peut éventuellement définir une densité uniforme sur une région bornée qui recouvre une région de confiance suffisamment large (cette idée est développée dans la Section 3.2). On utilise ensuite un algorithme de simulation par chaînes de Markov [8] pour échantilloner une population suivant la densité $C$ concentrée au voisinage de la surface d'état-limite approximée. Plus précisément, on utilise ici la technique baptisée *slice sampling* [9]. Plutôt que de tirer seulement le nombre $N_{\text{enr}}$ de points destinés à enrichir le plan d'expériences à l'itération courante, on simule par ce procédé une population plus grande (typiquement $10^3 - 10^4$ points), que l'on réduit ensuite à $N_{\text{enr}}$ points représentatifs de la population initiale par une technique de *K-means clustering* [10]. La fonction de performance est évaluée en ces $N_{\text{enr}}$ nouveaux points pour mettre à jour le méta-modèle. Les détails de la mise en œuvre de cette technique sont plus largement explicités dans [7].

L'opération est répétée de manière itérative jusqu'à ce que l'état-limite approximé soit suffisamment précis. La façon dont ce degré de précision est mesuré dans le cadre d'une application fiabiliste est détaillé ci-après.

## 2.3 Critère de précision sur la probabilité de défaillance

Dans le cadre d'une analyse de fiabilité, on vient évaluer la probabilité de défaillance sur la prédiction moyenne du krigeage par simulations – c'est ce que l'on appelle *la substitution*. Dans le but de réduire encore le temps de calcul total, on propose d'utiliser une technique de réduction de variance robuste connue sous le nom de *subset simulation* [4] en lieu et place des simulations de Monte-Carlo classiques qui s'avèrent trop coûteuses lorsque la probabilité recherchée est faible.

Cependant, la probabilité de défaillance évaluée sur la prédiction moyenne du krigeage peut être erronée par rapport à la probabilité de défaillance réelle que l'on aurait pu estimer sur le modèle réel. Et ce *biais* s'avère difficile à estimer en pratique. C'est pourquoi on propose de calculer les probabilités de défaillance sur les surfaces d'état-limite "limites" associées aux bornes de la marge de prédiction définie $\mathbb{M}$ à l'éq. (13), en plus de celle évaluée sur la prédiction moyenne, en vue de propager l'incertitude épistémique de méta-modélisation jusqu'à la grandeur d'intérêt. Les trois probabilités sont ainsi définies :

$$p_f^i = \mathbb{P}(\mathbf{X} \in \hat{\mathbf{F}}^i), \quad i = -1, 0, +1 \quad (15)$$

pour les 3 régions de défaillances approchées définies par :

$$\hat{\mathbf{F}}^i = \{ \boldsymbol{x} : \mu_{\hat{\mathbf{G}}} (\boldsymbol{x}) + i k \sigma_{\hat{\mathbf{G}}} (\boldsymbol{x}) \leq 0 \}, \quad i = -1, 0, +1. \quad (16)$$

Etant donné que l'écart-type de krigeage est toujours positif, la proposition suivante est valide :

$$\hat{\mathbf{F}}^{+1} \subseteq \hat{\mathbf{F}}^{0} \subseteq \hat{\mathbf{F}}^{-1} \Rightarrow p_f^{+1} \le p_f^{0} \le p_f^{-1}. \quad (17)$$

Il est néanmoins important de signaler ici qu'il n'existe aucune preuve que la probabilité de défaillance satisfasse cet encadrement de par les hypothèses heuristiques faites jusqu'à présent. Pourtant,
l'étendue de l'intervalle $[p_f^+, p_f^-]$ s'avère être une bonne mesure pratique de la précision de l'estimation de la probabilité de défaillance sur la surface krigée. Il est donc proposé d'enrichir le métamodèle tant que l'étendue logarithmique de l'intervalle $[p_f^+, p_f^-]$ est plus grande qu'un seuil de tolérance fixé. La condition d'arrêt pour la procédure d'enrichissement itérative détaillée ci-avant s'écrit donc :

$$ \varepsilon_{p_f} = \log \left( \frac{p_f^-}{p_f^+} \right) \leq \varepsilon_{p_f 0}. \qquad (18) $$

En pratique, on visera des étendues logarithmiques $\varepsilon_{p_f 0}$ inférieures à 1, ce qui signifie que l'estimation de la probabilité de défaillance est au moins dans l'ordre de grandeur de la probabilité réelle.

# 3 Algorithme RBDO

## 3.1 Résolution du problème d'optimisation sous contrainte fiabiliste

Le problème d'optimisation sous contrainte est résolu de manière itérative au moyen de l'algorithme de Polak-He [11]. Cet algorithme procède itérativement à partir d'un dimensionnement initial $\theta^{(0)}$. La direction d'optimisation est déterminée à chaque itération par résolution d'un sous-problème quadratique approximé. La matrice hessienne intervenant dans le développement limité du Lagrangien est approximée à partir des gradients des fonctions objectifs et contraintes. Une fois la direction déterminée, le pas d'optimisation est obtenu par une recherche linéique approximée suivant la règle de Goldstein-Armijo. Le lecteur interessé pourra se référer à [11] pour plus de détails.

Le méta-modèle de chaque fonction de performance est construit et enrichi dans un espace de fiabilité dit “augmenté” à chaque itération de la boucle d'optimisation. Cet espace défini ci-après englobe l'ensemble des choix de conception admissibles $\theta \in \mathcal{D}_{\theta}$ et le support des lois de probabilités conditionnelles $f_X(x | \theta)$ des variables d'environnement. Le méta-modèle ainsi construit est ensuite supposé être suffisamment précis pour l'estimation des contraintes de fiabilité et de ses gradients pendant toute l'itération d'optimisation (calcul de la direction et recherche linéique). Les gradients des probabilités de défaillances sont estimées par un simple post-traitement de l'analyse par subset simulation tel que présenté dans [12].

## 3.2 Espace de fiabilité augmenté

Selon [13], l'espace de fiabilité augmenté correspond à l'espace support décrit par les variables aléatoires d'environnement $X$ dont l'incertitude est augmentée artificiellement en supposant le vecteur des paramètres de conception $\theta$ uniforme ($\pi(\theta) = \frac{1}{|\mathcal{D}_{\theta}|} \mathbf{1}_{\mathcal{D}_{\theta}}(\theta)$) sur le domaine de conception admissible $\mathcal{D}_{\theta}$. On note le vecteur décrivant cet espace augmenté $\mathcal{D}_V$ par $v$, et sa densité de probabilité est définie comme suit :

$$ h(v) = \int_{\mathcal{D}_{\theta}} f_X(x|\theta) \pi(\theta) d\theta = \frac{1}{\mathcal{V}_{\mathcal{D}_{\theta}}} \int_{\mathcal{D}_{\theta}} f_X(x|\theta) d\theta \propto \int_{\mathcal{D}_{\theta}} f_X(x|\theta) d\theta \qquad (19) $$

en notant $\mathcal{V}_{\mathcal{D}_{\theta}}$ le volume de l'espace de conception admissible. On pourrait ensuite utiliser cette densité de probabilité comme densité de pondération pour définir le critère d'enrichissement $C$ présenté à la Section 2.3. Cependant, son calcul explicite, même à la constante de normalisation $\mathcal{V}_{\mathcal{D}_{\theta}}$ près, pour une loi jointe $f_X$ quelconque n'est pas trivial.

L'objectif est de construire les méta-modèles des fonctions de performance sur un domaine $\mathcal{D}_V$ suffisamment large permettant une estimation précise de la fiabilité pour n'importe quel choix de conception dans le domaine de conception admissible. Ainsi, une approximation hyperrectangulaire d'une région de confiance suffisamment large de la densité augmentée $h$ est suffisante. Cette approximation est construite comme le produit tensoriel des intervalles de confiance sur ses marginales. Les bornes inférieures et supérieures de ces intervalles de confiance sont définis comme suit :

$$ q_{V_i}^{-} = \min_{\theta \in \mathcal{D}_{\theta}} F_{X_i}^{-1}(\Phi(-\beta_0) | \theta), \quad i = 1, \dots, n \qquad (20) $$

$$ q_{V_i}^{+} = \max_{\theta \in \mathcal{D}_{\theta}} F_{X_i}^{-1}(\Phi(+\beta_0) | \theta), \quad i = 1, \dots, n \qquad (21) $$
où $\beta_0$ correspond à l'indice de fiabilité maximal que l'on souhaite pouvoir estimer au long de l'optimisation. En pratique, on fixe $\beta_0 = 8$ ce qui correspond à une région de confiance relativement large par rapport aux indices de fiabilité visés dans les applications (entre 3 et 5). $q_{\bar{V}_i}$ et $q_{V_i}^+$ peuvent être obtenus par optimisation numérique.

Finalement, on utilise la densité uniforme sur l'approximation rectangulaire de la région de confiance comme densité de pondération $w$ dans le critère $C$ (éq.(14). La densité uniforme, à sa constante de normalisation près, est simplement définie comme la fonction indicatrice suivante :

$$ \mathbb{1}_{\mathcal{D}_V, \beta_0}(v) = \begin{cases} 1 & \text{si } q_{\bar{V}_i} \le v_i \le q_{V_i}^+, \quad i = 1, \dots, n \\ 0 & \text{sinon} \end{cases} \qquad (22) $$

# 4 Optimisation sous contrainte de fiabilité d'un pont en treillis

Les modèles mécanique et probabiliste de cet exemple sont inspirés de [14].

## 4.1 Modèle mécanique

Il s'agit d'étudier le dimensionnement de la structure en treillis illustrée sur la Figure 1. Elle est composée de 11 barres horizontales et 12 barres obliques. La partie horizontale supérieure est soumise à un chargement vertical. Le déplacement vertical maximum de la structure (repéré $V_1$ sur la Figure 1) est obtenu par résolution d'un modèle aux éléments finis de barres plan et linéaire. Le modèle est ainsi constitué de 23 éléments et 13 noeuds comportant chacun 2 degrés de libertés avant application des conditions aux limites.

<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Pont en treillis – Modèle mécanique</figcaption></figure>

## 4.2 Modèle probabiliste

Le modèle probabiliste initial [14] est spécifié dans le Tableau 1. Il comprend 10 variables aléatoires indépendantes modélisant les modules d'Young et les sections des deux groupes de barres (horizontales et obliques), ainsi que les chargements aux 6 noeuds de la partie supérieure. Les moyennes initiales des sections sont respectivement fixées à $\mu_{A_1} = 2 \times 10^{-3} \text{ m}^2$ et $\mu_{A_2} = 1 \times 10^{-3} \text{ m}^2$.

<table><caption>TABLE 1 – Pont en treillis – Modèle probabiliste</caption><thead><tr><th>Variable</th><th>Distribution</th><th>Moyenne</th><th>Écart Type</th></tr></thead><tbody><tr><td>$E_1, E_2$ (Pa)</td><td>Lognormale</td><td>2, 10 &times; 10<sup>11</sup></td><td>2, 10 &times; 10<sup>10</sup></td></tr><tr><td>$A_1$ (m<sup>2</sup>)</td><td>Lognormale</td><td>$\mu_{A_1}$</td><td>2 &times; 10<sup>-4</sup></td></tr><tr><td>$A_2$ (m<sup>2</sup>)</td><td>Lognormale</td><td>$\mu_{A_2}$</td><td>1 &times; 10<sup>-4</sup></td></tr><tr><td>$P_1, \dots, P_6$ (N)</td><td>Gumbel</td><td>5 &times; 10<sup>4</sup></td><td>7, 5 &times; 10<sup>3</sup></td></tr></tbody></table>

Le déplacement à la moyenne du modèle probabiliste initial détaillé ci-dessus vaut $v_1(\mu) = -7,78 \text{ cm}$. La règle de dimensionnement veut que le déplacement maximal ne dépasse pas $v_{\text{max}} = 10 \text{ cm}$, de sorte que la fonction de performance s'écrit :

$$ g(\boldsymbol{x}) = v_{\text{max}} - |v_1(\boldsymbol{x})| \qquad (23) $$
## 4.3 Formulation du problème d'optimisation

Le problème d'optimisation consiste à minimiser le volume total de la structure en agissant sur la moyenne des sections des deux groupes de barres $\mu_{A_1}$ et $\mu_{A_2}$ (exprimées en mètres et regroupées dans le vecteur $\theta$) tout en cherchant à garantir un indice de fiabilité minimal $\beta_0$ égal à 3. Le problème d'optimisation ainsi formulé s'écrit :

$$ \theta^* = \arg \min_{\theta \in \mathcal{D}_{\theta}} L_1 \mu_{A_1} + L_2 \mu_{A_2} \quad \text{telque} \quad \beta_0 \le \beta(\theta) \quad (24) $$

où $L_1$ et $L_2$ sont les sommes des longueurs des barres par groupe (barres horizontales et obliques respectivement). Le domaine de conception admissible $\mathcal{D}_{\theta}$ est défini comme suit :

$$ \mathcal{D}_{\theta} = \left\{ \theta \in \mathbb{R}^2 : \begin{array}{l} 6 \times 10^{-4} \le \mu_{A_1} \le 6 \times 10^{-3} \\ 3 \times 10^{-4} \le \mu_{A_2} \le 3 \times 10^{-3} \end{array} \right\} \quad (25) $$

L'indice de fiabilité correspond au choix intial de conception ($\mu_{A_1} = 2 \times 10^{-3} \text{ m}^2$ et $\mu_{A_2} = 1 \times 10^{-3} \text{ m}^2$) vaut $\beta_{ini} = 1,68$, pour un coût associé de 0,1219 $\text{m}^3$ de matière.

## 4.4 Résultats

La Figure 2 illustre la bonne convergence de la stratégie proposée. Les traits pointillés qui encadrent la courbe des indices de fiabilité par itération correspondent aux indices de fiabilités évalués sur les bornes de l'intervalle à 95% sur la prédiction de l'états-limite.

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Pont en treillis – Convergence de la stratégie proposée</figcaption></figure>

Les résultats obtenus par la stratégie RBDO présentée ci-avant sont confrontés aux résultats obtenus par deux autres stratégies double-boucle basées sur FORM et sur subset simulation – cette dernière analyse constituant une référence. Les résultats sont données dans le Tableau 2. Tous les résultats sont présentés en termes d'indices de fiabilité généralisés, et on donne également le coefficient de variation (CoV) lorsque la probabilité associée est estimée par simulation. Les dimensionnements optimaux sont exprimés en mètres et le coût en mètres cubes. La stratégie proposée fournit des résultats équivalents à ceux obtenus par subset simulation sur le modèle réel – i.e. sans méta-modèle, tout en permettant de limiter considérablement le coût de calcul exprimé dans le Tableau 2 en nombre d'appels au modèle aux éléments finis.

FORM semble fournir ici des résultats légèrement non conservatifs. Une analyse SORM Point Fitting [15] au dimensionnement optimal donné par FORM confirme l'unicité du point de conception, mais l'état-limite présente un léger repli sur le domaine de sûreté qui diminue l'indice de fiabilité : $\beta_{SORM-PF} =$
<table><caption>TABLE 2 – Pont en treillis – Résultats de l'optimisation sous contrainte de fiabilité</caption><thead><tr><th rowspan="2">Méthode</th><th colspan="2">Dimensionnement</th><th rowspan="2">Coût</th><th rowspan="2">Indice de fiabilité (CoV)</th><th rowspan="2">Nb d'appels à g</th></tr><tr><th>μ<sup>*</sup><sub>A1</sub></th><th>μ<sup>*</sup><sub>A2</sub></th></tr></thead><tbody><tr><td>Krigeage & sub. sim.</td><td>2,53 × 10<sup>-3</sup></td><td>8,13 × 10<sup>-4</sup></td><td>0,1388</td><td>3,05 (2%)</td><td>350</td></tr><tr><td>FORM</td><td>2,43 × 10<sup>-3</sup></td><td>7,99 × 10<sup>-4</sup></td><td>0,1341</td><td>2,99 (N/A)</td><td>1621</td></tr><tr><td>Sub. sim.</td><td>2,53 × 10<sup>-3</sup></td><td>7,98 × 10<sup>-4</sup></td><td>0,1383</td><td>3,05 (8%)</td><td>317060</td></tr></tbody></table>

$2,74 < β_{HL} = 2,99$. On notera néanmoins, que l'approche proposée reste compétitive par rapport à une stratégie double-boucle basée sur FORM grâce, notamment, à l'espace global dans lequel le méta-modèle est construit.

## 5 Conclusions

On propose dans cette communication un algorithme d'optimisation sous contrainte de fiabilité basé sur l'utilisation de méta-modèles par krigeage pour approximer les fonctions de performance définissant les contraintes fiabilistes. Le calcul de la probabilité de défaillance courante à chaque itération de l'optimisation peut alors se faire par une méthode de simulation avancée (subset simulation) ce qui permet de s'affranchir des approximations comme FORM. On montre sur l'exemple d'application (optimisation des sections des barres d'un treillis) que le coût de calcul, en termes de nombre d'appels au code éléments finis, peut être divisé par 1000 du fait de l'utilisation du méta-modèle, sans dégradation de la précision.

## Références

[1] Y. Tsompanakis, N.D. Lagaros, and M. Papadrakis, editors. *Structural design optimization considering uncertainties*. Taylor & Francis, 2008.

[2] I. Enevoldsen and J. D. Sorensen. Reliability-based optimization in structural engineering. *Structural Safety*, 15(3):169–196, 1994.

[3] T.J. Santner, B.J. Williams, and W. Notz. *The design and analysis of computer experiments*. Springer series in Statistics. Springer, 2003.

[4] S.K. Au and J.L. Beck. Estimation of small failure probabilities in high dimensions by subset simulation. *Prob. Eng. Mech.*, 16(4):263–277, 2001.

[5] V. Picheny, D. Ginsbourger, O. Roustant, and R. Haftka. Adaptive designs of experiments for accurate approximation of a target region. *J. Mech. Des.*, 132(7), 2010.

[6] B. Echard, N. Gayton, and M. Lemaire. AK-MCS: an Active learning reliability method combining Kriging and Monte Carlo Simulation. *Structural Safety*, In press, 2011.

[7] V. Dubourg, B. Sudret, and J.-M. Bourinet. Reliability-based design optimization using kriging and subset simulation. Submitted to Struct. Multidisc. Optim., 2011.

[8] C.P. Robert and G. Casella. *Monte Carlo statistical methods (2nd Ed.)*. Springer Series in Statistics. Springer Verlag, 2004.

[9] R.M. Neal. Slice sampling. *Annals Stat.*, 31:705–767, 2003.

[10] C. Elkan. Using the triangle inequality to accelerate K-means. In *Proceedings of the Twentieth International Conference on Machine Learning (ICML-2003)*, Washington D.C., 2003.

[11] E. Polak. *Optimization algorithms and consistent approximations*. Springer, 1997.

[12] S. Song, Z. Lu, and H. Qiao. Subset simulation for structural reliability sensitivity analysis. *Reliab. Eng. Sys. Safety*, 94(2):658–665, 2009.

[13] A.A. Taflanidis and J.L. Beck. Stochastic subset optimization for reliability optimization and sensitivity analysis in system design. *Computers & Structures*, 87(5-6):318–331, 2009.

[14] G. Blatman and B. Sudret. Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach. *Comptes Rendus Mécanique*, 336(6):518–523, 2008.

[15] A. Der Kiureghian, H.-Z. Lin, and S.-J. Hwang. Second order reliability approximations. *Journal of Engineering Mechanics*, 113(8):1208–1225, 1987.