# Surfaces de réponse physiques et polynomiales

Frederic Duprat, Franck Schoefs, Bruno Sudret

► To cite this version:

Frederic Duprat, Franck Schoefs, Bruno Sudret. Surfaces de réponse physiques et polynomiales. Julien Baroth, Franck Schoefs, Denys Breysse. Fiabilité des ouvrages Sûreté, variabilité, maintenance, sécurité, pp.165-193, 2012, 9782746231443. hal-01007518

HAL Id: hal-01007518

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

Submitted on 16 Jun 2014

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
# Surfaces de réponse physiques et polynomiales

Frédéric DUPRAT, Franck SCHOEFS et Bruno SUDRET

## 7.1. Introduction

L'analyse de fiabilité d'une structure s'appuie classiquement sur la donnée de modèles mécanique et probabiliste, ainsi qu'une fonction d'état limite. On définit tout d'abord un modèle mécanique décrivant le comportement de la structure. C'est au sens large une fonction de transfert mathématique $M$ permettant d'évaluer les effets de sollicitations en fonction de paramètres d'entrée, ou stimuli, décrivant la structure et son environnement, paramètres regroupés dans un vecteur $x$. On note ici $y = M(x)$ la réponse du modèle.

On définit ensuite un modèle probabiliste pour les paramètres d'entrée qui sont considérés comme mal connus ou incertains, soit à partir du traitement statistique d'échantillons de données lorsque ceux-ci sont disponibles, soit à partir de jugement d'expert et de bases de données [JCS 02]. Ce modèle probabiliste est caractérisé par la densité jointe des variables aléatoires d'entrée $X$, notée $f_X(x)$.

On définit enfin une fonction d'état limite traduisant mathématiquement le critère de défaillance par rapport auquel la structure doit être justifiée. Cette fonction s'écrit sous la forme générale $g(M(X), X')$ et porte sur les effets des sollicitations auxquels on fixe des limites, regroupées dans un vecteur $X'$.
En notant $f_{X,X'}$ la densité jointe des vecteurs $X$ et $X'$, l'objectif de l'analyse est d'évaluer la probabilité de défaillance $P_f$ de la structure définie par :

$$P_f = \int_{\{x: g(M(x), x') \le 0\}} f_{X,X'}(x, x') dx \quad [7.1]$$

Il existe diverses méthodes pour résoudre le problème [DIT 96, LEM 05]. Parmi celles-ci, dans ce chapitre, seules la simulation de Monte-Carlo et la méthode FORM sont présentées et utilisées.

## 7.2. Historique de la méthode des surfaces de réponse

De nombreux domaines scientifiques et courants de pensée ont contribué à l'élaboration de ce qui est appelé méthodologie des surfaces de réponse (MSR). Les prémices de l'approche par des surfaces de réponse apparaissent quelques années avant les développements de Box et Wilson [BOX 50]. Des domaines scientifiques variés étaient concernés :

- la biologie animale et végétale avec la construction de courbes de croissance [WIS 39, WIN 32, REE 29] ;

- les sciences humaines et l'analyse de la réponse d'une population à des stimuli [GAD 33, BLI 35a et b], travaux basés sur ceux du psychiatre Fechner vers 1860 ;

- l'agronomie et l'étude de fertilisation des sols [MIT 30, CRO 35, STE 51].

Ces démarches reposaient principalement sur l'hypothèse de base (justifiée mathématiquement par le théorème de Weierstrass) que, sous certaines conditions de régularité, la réponse peut être approchée par une forme déduite des polynômes. C'est donc naturellement dans ce cadre que les chimistes Box et Wilson ont développé en 1951 le concept de surface de réponse en s'appuyant sur des techniques de régression analytique et la conduite d'expériences. En particulier le soin à apporter au choix des variables stimuli et à l'affectation de leur poids relatif fut déjà largement décrit. Ces modèles expérimentaux se sont enrichis par la définition des périodes d'observation [BOX 55] et par un calcul d'erreur [BOX 57]. Avec l'accroissement du nombre de modèles potentiels, il a fallu en parallèle proposer des critères de choix [BOX 59a] comme la minimisation de la variance généralisée des estimations des variables [BOX 59b]. La période 1950-1970 a été marquée par l'apparition de trois courants majeurs d'inspiration probabiliste :

- la recherche d'une représentation fonctionnelle optimale par approximation stochastique et ce en présence de points non représentatifs, avec une extension à des configurations multivariables [ROB 51, KIE 52] ;

- la comparaison des courbes de croissance en biométrie [RAO 58, ELS 62], où la fonction de réponse résulte de projections sur une famille de polynômes orthogonaux. Leurs coefficients sont ensuite utilisés pour des études de prévision ;
– la théorie du modèle optimal sous contraintes dans le cas de modèles linéaires [KIE 60]. L’optimisation des fonctions de réponse est reliée à la minimisation de la variance généralisée des paramètres.

Ce dernier courant a constitué un apport théorique bien structuré qui sert encore de référence. De nombreux travaux d'optimisation sont venus ensuite [ANS 63, POW 65, NEL 65]. La fin de cette période a vu l'émergence de modèles non linéaires. Les travaux les plus significatifs ont été certainement ceux menés sur les polynômes inverses [NEL 66] et l'estimation statistique des paramètres à partir de processus expérimentaux [ATK 68]. Cependant, c'est plus récemment, grâce à la puissance croissante des moyens de calculs, que l'utilisation des modèles non linéaires a connu un réel essor. Les critères de validation des modèles y apparaissent comme très spécifiques à chaque domaine d'activité. On peut signaler les articles de synthèse en biométrie [MEA 75], et dans le secteur du nucléaire [HEL 93].

## 7.3. Concept de surface de réponse

### 7.3.1. Définitions de base

L'appellation de surface de réponse dénote le souci de visualiser par une représentation géométrique ce qu'est la réponse d'un processus physique à des stimuli. La propriété étudiée, ou réponse Y, résulte alors de l'application de la fonction de transfert caractérisant le système aux paramètres d'entrée. Elle varie par conséquent avec les variations de ces paramètres d'entrée, appelés aussi stimuli. Ces derniers sont modélisés par des variables ou champs aléatoires, notés $X_i$, $i=1,..,n$, et sont caractérisés par un ensemble d'informations statistiques disponibles, notées $\theta_j$, $j=1,..,p$, (fonctions de densité de probabilité indépendantes ou corrélées, moments normalisés, etc.). Ce transfert peut se représenter tel que :

<figure><img src="image_2.png" /><figcaption>**Figure 7.1.** Réponse d'une fonction de transfert à un stimuli, modélisé par des variables ou champs aléatoires</figcaption></figure>

Une modélisation probabiliste par des champs aléatoires est parfois nécessaire pour représenter les variations spatio-temporelles des paramètres incertains d'entrée. Une modélisation par des variables aléatoires, plus simple, peut néanmoins suffire dans un certain nombre de problèmes. Par commodité, nous utilisons par la suite des variables aléatoires pour modéliser les paramètres d'entrée.
Généralement, on n'accède pas à cette fonction de transfert sous une forme explicite des variables de base. Nous en cherchons une approximation, appelée fonction de réponse *M*, le plus souvent dans une famille de fonctions usuelles linéaires ou non, caractérisées par des paramètres $\chi_k$, $k=1,...,l$, aléatoires ou déterministes. Ces paramètres sont déduits de l'ajustement de la réponse rendu possible par l'existence de données expérimentales. La représentation géométrique sous la forme d'une courbe, d'une surface, ou d'une hypersurface est appelée surface de réponse. L'introduction de notions géométriques comme les lignes de niveau relatives à cette surface peut alors se traduire en termes de frontière du domaine de sécurité. Pour construire une surface de réponse, il faut fournir :

- $X = \{X_1, ..., X_n\}$, ensemble représentatif et ordonné de variables aléatoires ;

- $\theta = \{\theta_1, ..., \theta_p\}$, ensemble d'informations statistiques sur *X* (fonctions de densité de probabilité indépendantes ou corrélées, moments normalisés, etc.) ;

- $M(X/\theta)$, approximation de la réponse *Y*, formulée explicitement en fonction de *X* sachant $\theta$, et obtenue par ajustement de l'ensemble des paramètres $\chi$ ;

- $|\cdot|$, une métrique dans l'espace des variables de base et de la réponse. La qualité de l'ajustement de l'approximation *M* à la réponse *Y* peut être ainsi mesurée.

La fonction de réponse peut donc s'écrire formellement :

<figure><img src="image_3.png" /><figcaption>**Figure 7.2.** *Ecriture formelle d'une fonction de réponse*</figcaption></figure>

### 7.3.2. Différentes formulations

Le choix d'un type de formulation de la fonction de réponse se fait à partir de critères spécifiques résultant de la méthodologie scientifique retenue pour étudier le phénomène. Les premières questions sont le niveau de complexité envisageable, la possibilité d'une approche expérimentale, et l'état actuel des connaissances servant de références. On peut noter deux orientations possibles mais la construction de surfaces de réponse relève de plus en plus d'une solution mixte :

- l'ajustement d'un modèle approché du transfert en utilisant des fonctions mathématiques usuelles (notamment polynomiales) et une base de données sélectionnées. [SCH 96] présente et compare les modèles usuels ;
– l’utilisation des lois physiques déterministes sous-jacentes dans lesquelles sont introduites des variables aléatoires pour rendre compte des variabilités intrinsèques (comme la hauteur ou la période d’une vague) ou des incertitudes sur les paramètres des modèles [LAB 96]. Deux questions dirigent alors la construction de la surface de réponse, la représentativité physique des modèles déterministes et la sélection des variables de base, c’est-à-dire celles qui gouvernent les variations des quantités étudiées.

Enfin, à la sélection de la formulation analytique du transfert s’ajoute la difficulté de caractériser statistiquement les variables de base. Dans le cas des surfaces de réponse décrivant des états limites, cette question conditionne la mesure de la fiabilité. Par exemple, dans le domaine des matériaux composites laminés unidirectionnels et chargés dans le plan, on considère souvent qu’un état limite peut se déduire de trois critères que sont : les contraintes, les efforts externes plans, et les dimensions géométriques du matériau. Des hypothèses différentes sur le nombre de variables et leur type de loi de distribution peuvent conduire à des variations de plus de 30 % sur les valeurs des facteurs de sécurité correspondants [NAK 95].

### 7.3.3. Critères de construction

On peut penser que les critères de construction des surfaces de réponse sont spécifiques à chaque domaine d’utilisation. Aussi, dans ce qui suit, les critères sont abordés par ordre d’importance pour les problèmes relatifs à la sécurité des ouvrages et structures. Les critères ainsi mis en évidence soulèvent des questions qui sont néanmoins généralisables. Par ailleurs, il est illusoire de considérer qu’une solution unique dans l’absolu puisse exister. Le choix final résulte d’une optimisation sous contraintes dont on propose ici de spécifier la nature.

On cherche notamment à profiter de la puissance croissante des moyens de calcul qui offre des facilités en permettant bien souvent d’affiner la représentation mathématique et de contrôler les incertitudes inhérentes à l’ajustement du modèle. Ceci ne doit cependant pas masquer les impératifs de représentativité physique des phénomènes. Les critères considérés comme majeurs dans notre approche sont :

- le sens physique de la représentation ;

- les effets du choix de la modélisation probabiliste ;

- la mesure de la qualité de l’ajustement ;

- la réduction du niveau de complexité et l’adaptation au calcul numérique.
### 7.3.3.1. Sens physique de la représentation

La compréhension du mécanisme sous-tendant un phénomène physique est fondamentale dans le choix de l'ensemble des variables d'entrée et de la fonction d'approximation. Ce critère peut révéler la nécessité de baser la formulation de la fonction de réponse sur des lois déterministes. La variabilité intrinsèque s'exprime alors par l'introduction de variables aléatoires et l'incertitude sur le modèle est prise en compte par l'intermédiaire de paramètres aléatoires. L'utilisation de lois déterministes et une sélection soignée des variables de base, pour peu qu'une formulation analytique soit accessible, se montrent plus représentatives que le recours à des modèles d'ajustement par des fonctions d'approximation usuelles.

### 7.3.3.2. Effets du choix de la modélisation probabiliste

La loi de probabilité de la réponse du système dépend de la caractérisation probabiliste des paramètres d'entrée (loi de probabilité, dispersion, dissymétrie, aplatissement, etc.). Dans un cas simple où la fonction de transfert M est linéaire, en fonction de variables aléatoires gaussiennes, la loi de la réponse est aussi gaussienne. Dans un cas général, les variables aléatoires sont d'une loi non gaussienne et la fonction de transfert d'un certain degré de non linéarité.

Afin de contrôler l'effet de ce degré de non linéarité de la fonction de transfert, divers travaux proposent d'approcher cette fonction par une forme polynomiale (en général linéaire, quadratique ou cubique). Une approximation linéaire est souvent insuffisante [BOU 95] ; les approximations cubiques et de degré cinq permettent parfois d'évaluer de manière satisfaisante les moments respectivement jusqu'à l'ordre 3 et jusqu'à l'ordre 4. Les approximations polynomiales pour le type de transfert étudié peuvent être de degré élevé et donc coûteuses en temps de calcul.

En fait, pour garantir le transfert des lois de probabilité, on va montrer qu'il est nécessaire de veiller au bon ajustement de la matrice jacobienne $[D(X)/D(Y)]$ [LAB 95, SCH 08]. On considère pour cela une surface de réponse à une variable de la forme $Y = M(X)$ avec M dérivable, bijective et monotone, et X variable de base. Connaissant la densité de probabilité $f_X$ de X, on peut calculer $f_Y$, densité de probabilité de Y. Soit G la fonction de répartition associée à la densité de probabilité g, on sait par définition :

$$G_Y(y) = P(Y < y) \text{ et } f_Y(y) = \frac{d}{dy} P(Y < y)$$
$$f_Y(y) = P[X < M^{-1}(y)] \text{ soit enfin,} \qquad f_Y(y) = f_X[M^{-1}(y)] \left| \frac{1}{M'}[M^{-1}(y)] \right|$$

où $|\cdot|$ désigne la valeur absolue.

Dans le cas plus général d'un problème multivariable pour lequel X et Y sont des vecteurs, les fonctions de densité de probabilité des variables aléatoires d'entrée et des variables aléatoires de sorties sont donc liées par la relation :

$$f_Y(y) = f_X(x) [M^{-1}(y)] \left| \frac{D(X)}{D(Y)} \right|$$

Ainsi, pour obtenir une bonne approximation de $P_Y$, il faut bien ajuster la fonction $M$, mais aussi rechercher une adéquation avec les dérivées partielles. Les formes polynomiales linéaires, quadratiques ou cubiques évoquées précédemment ont respectivement des fonctions dérivées constantes, linéaires et quadratiques. L'écart entre ces trois fonctions est en général important sur les bords du domaine d'étude et les perturbations qu'elles apportent sur les queues des fonctions de densité de probabilité peuvent donc être significatives et modifier les résultats du calcul fiabiliste. Ainsi, le choix d'un polynôme de degré peu élevé, pratique pour la conduite des calculs et apparemment adapté à une approximation de la réponse, peut mener à des densités de probabilité erronées en sortie. Un tel choix s'avère très sensible aux effets du choix de la modélisation probabiliste des paramètres d'entrée.

### 7.3.3.3. Mesure de la qualité de l'ajustement

Il s'agit ici de définir une métrique fournissant un bon indicateur de la qualité de l'ajustement effectué. Les métriques usuelles sont dites du second ordre car elles permettent d'obtenir les variables réellement prépondérantes sur la réponse en évaluant leur influence sur l'écart de la variance en sortie. Elles sont donc peu performantes en présence d'événements singuliers externes aux lois de probabilité d'entrée et ne sont qu'indicatives en tant que mesure de l'incertitude [BIE 83, IMA 87]. Pour pallier ces insuffisances, il existe des métriques basées sur les écarts inter-quantiles [KHA 89] ou sur la mesure de l'entropie du système [PAR 94].

Une autre démarche consiste à utiliser les métriques de régression. Notons $f$ la fonction de réponse dont il faut ajuster les paramètres et $\varepsilon$ l'erreur. Considérons par exemple un modèle de régression du type :

$$Y = M(X/\theta) + \varepsilon$$
Lorsque l'erreur $\varepsilon$ est supposée suivre une distribution normale de moyenne zéro et de matrice de covariance diagonale, alors les ajustements de $f$ aux moindres carrés et au maximum de vraisemblance sont identiques. La métrique $L_2$ (intégration du carré du résidu $u$) est alors la plus efficace. La métrique $L_1$ (intégration de la norme du résidu) est quant à elle mieux adaptée à une distribution exponentielle de l'erreur, alors nettement plus dispersée.

Nous avons conclu au paragraphe précédent à la nécessité d'ajuster la matrice jacobienne des dérivées partielles. Ainsi, il paraît particulièrement intéressant de choisir une métrique de la forme rencontrée en théorie variationnelle :

$$ \|u\| = \sqrt{\|u\|_{L_2}^2 + \sum_{i=1,...,n} \left\| \frac{\partial u}{\partial X_i} \right\|_{L_2}^2} $$

avec $\|\cdot\|_{L_2}$, norme au sens $L_2$, dans l'espace de Sobolev $H_1$ :

$$ H_1 = \{ u \mid \frac{\partial^\alpha u}{\partial X_i^\alpha} \in L^2 \text{ (carré intégrable) quels que soient } \alpha = 0, 1 \text{ et } i = 1...,n \} $$

L'idée sous-jacente de ce choix de métrique est donc de préférer le contrôle des queues de distribution par transferts successifs que de privilégier la partie centrale par l'intermédiaire des premiers moments. Toute construction de surface de réponse doit être accompagnée d'un choix de métrique qui conditionne le sens de l'approximation et permet d'expliquer certaines limites dans la représentation.

### 7.3.3.4. Réduction du niveau de complexité et adaptation au calcul numérique

Pour plus de précision, nous pouvons introduire des ordres d'approximation élevés dans les procédures de calcul. Compliquer les procédures doit cependant être justifié, ne serait-ce que par l'impératif d'ingénierie de ne pas augmenter inutilement les temps des calculs.

Cette question est illustrée dans [SCH 07] au travers de plusieurs études portant sur l'interaction houle-structure. Celles-ci étudient les effets de l'ordre du modèle cinématique de Stokes, de la prise en compte du terme d'inertie des efforts et du nombre de tronçons de discrétisation de l'élément-poutre.
## 7.4. Méthodes de calcul de fiabilité

### 7.4.1. Problème de fiabilité et simulation de Monte-Carlo

La fiabilité des structures vise à calculer la probabilité de défaillance d'un système mécanique dont les paramètres sont mal connus et modélisés dans un cadre probabiliste, relativement à un scénario de défaillance. Les méthodes de fiabilité permettent non seulement de calculer la probabilité de défaillance $P_f$ d'un composant ou d'un système constitué de plusieurs composants en interaction, mais également de déterminer la sensibilité de cette probabilité à chaque variable aléatoire intervenant dans le problème [DIT 96, LEM 05].

Un critère de défaillance peut être représenté par une fonction de performance, notée $g: \mathbb{R}^M \rightarrow \mathbb{R}$ définie conventionnellement dans l'espace des paramètres de sorte que $D_f = \{\mathbf{x}: g(\mathbf{x}) \le 0\}$ corresponde au domaine de défaillance et $D_s = \{\mathbf{x}: g(\mathbf{x}) > 0\}$ au domaine de sûreté. La frontière $\partial D = \{\mathbf{x}: g(\mathbf{x}) = 0\}$ est la surface d'état limite. La probabilité de défaillance est alors définie par :

$$P_f = \mathbb{P}(g(\mathbf{X}) \le 0) = \int_{D_f} f_X(\mathbf{x}) d\mathbf{x} \qquad [7.2]$$

où $f_X$ est la densité de probabilité jointe de $X$.

Le domaine d'intégration $D_f$ étant défini implicitement et la fonction de performance analytique, dans le cas où la fonction g est définie à partir d'une surface de réponse, le calcul direct de cette intégrale peut être effectué numériquement par simulation de Monte-Carlo : $N_{sim}$ réalisations du vecteur d'entrée $X$ sont simulées, et pour chaque tirage, la fonction g est évaluée. On estime $P_f$ par le rapport $N_f / N_{sim}$, où $N_f$ est le nombre de tirages ayant conduit à la défaillance (signe négatif de g). Cette méthode simple dans son principe est très coûteuse en pratique : pour estimer une probabilité de l'ordre de $10^{-k}$ à 5 % près, il faut de l'ordre de $N_{sim} \approx 4.10^{k+2}$ tirages. Les probabilités calculées en fiabilité étant de l'ordre de $10^{-2}$ à $10^{-6}$, la méthode n'est souvent pas applicable directement à l'échelle industrielle, en particulier si l'on ne dispose pas d'une surface de réponse pertinente. Pour remédier à ce problème, diverses autres méthodes ont été développées parmi lesquelles la méthode FORM (First Order Reliability Method), très utilisée dans l'industrie.
### 7.4.2. Méthode FORM

La méthode FORM permet d'approcher la probabilité de défaillance à partir d'une reformulation du problème dans l'espace gaussien centré réduit, espace dans lequel toutes les variables aléatoires $\xi$ sont de loi gaussienne de moyenne nulle et d'écart type unitaire. Pour cela, on utilise une transformation isoprobabiliste $T: X \to \xi(X)$. Si les variables d'entrées sont indépendantes, de fonction de répartition marginale $F_{X_i}(x_i)$, cette transformation s'écrit simplement $\xi_i = \Phi^{-1}(F_{X_i}(x_i))$, où $\Phi$ est la fonction de répartition gaussienne centrée réduite. Dans le cas général, on peut utiliser les transformations de Nataf ou Rosenblatt [LEM 05, chapitre 4]. L'équation [7.2] se réécrit :

$$P_f = \int_{\{\xi: G(\xi)=g(T^{-1}(\xi))\le 0\}} \varphi_M(\xi) d\xi_1 \dots d\xi_M \quad [7.3]$$

où $G(\xi) = g(T^{-1}(\xi))$ est la fonction de performance dans l'espace réduit et $\varphi_M$ est la densité de probabilité multinormale centrée réduite de dimension $n$, définie par
$$\varphi_M(\mathbf{x}) = (2\pi)^{-\frac{M}{2}} \exp\left[-\frac{1}{2}(x_1^2 + \dots + x_M^2)\right].$$

Cette densité de probabilité est maximale à l'origine et décroît exponentiellement avec $\|\xi^2\|$ – cette décroissance diminue fortement avec le nombre de variables. Les points contribuant le plus à l'intégrale [7.3] sont donc les points appartenant au domaine de défaillance $D_f$ qui sont les plus proches de l'origine. La seconde étape de la méthode FORM consiste en la détermination du point de conception $P^*$, point du domaine de défaillance $D_f$ le plus proche de l'origine. Ce point est solution du problème d'optimisation :

$$P^* = \operatorname{Arg min}_{\xi \in \mathbb{R}^M} \left\{ \frac{1}{2} \|\xi\|^2 / G(\xi) \le 0 \right\} \quad [7.4]$$

Un algorithme approprié d'optimisation sous contraintes permet de le déterminer. L'indice de fiabilité $\beta$ est alors défini comme la distance algébrique de l'origine de l'espace réduit à la surface d'état-limite $\partial D : \beta = \operatorname{sign}(G(0)) \|\xi^*\|$. Une fois $\xi^*$ déterminé, on approche la surface d'état-limite $\partial D$ par un hyperplan tangent au point $P^*$ et le domaine de défaillance $D_f$ par le demi-espace associé.
L'intégrale [7.3] sur ce demi-espace conduit à l'approximation FORM de la probabilité de défaillance : $P_f \approx P_{f,FORM} = \Phi(-\beta)$. L'équation de l'état limite linéarisé (hyperplan) peut se mettre sous la forme $\tilde{G}(\xi) = \beta - \alpha^\top \cdot \xi$. Dans cette expression, le vecteur unitaire $\alpha$ des cosinus directeurs, qui est perpendiculaire à l'hyperplan, permet de définir les facteurs d'importance $\alpha_i^2$ de chaque variable aléatoire $\xi_i$, et par extension, de chaque $X_i$.

La méthode FORM présente donc l'avantage d'estimer $P_f$ pour un coût de calcul relativement faible (quelques dizaines à quelques centaines d'évaluations de la fonction $g$). L'approximation est souvent satisfaisante, et d'autant meilleure que $\beta$ est grand, à condition d'avoir bien identifié un point de conception $P^*$, sans extrema secondaires. De plus, cette méthode fournit les facteurs d'importance et des mesures de sensibilité (variation de la probabilité de défaillance en fonction des paramètres d'entrée) qui sont des indicateurs intéressants pour le concepteur.

Le lecteur est renvoyé à la lecture de [LEM 05] concernant la méthode SORM (Second Order Reliability Method), basée sur une approximation quadratique de l'équation de l'état limite. La méthode est a priori plus coûteuse que FORM, mais permet une approximation plus précise.

## 7.5. Surfaces de réponse polynomiales

Parmi les surfaces de réponse reposant sur l'ajustement de fonctions mathématiques sur des bases de données numériques, on a eu majoritairement recours en mécanique aux surfaces de réponse polynomiales. On se borne dans cette section aux surfaces de réponse polynomiales dont l'expression est un polynôme simple par opposition à celles construites sur la base d'un chaos polynomial qui font l'objet du chapitre 8.

### 7.5.1. Formulation générale

Si $X=\{X_i, i=1...,M\}$ dénote le vecteur des variables aléatoires, et si l'approximation $\hat{g}(X)$ de la fonction d'état limite $g(X)$ est quadratique :

$$ \hat{g}(X) = a_0 + \sum_{i=1}^{M} a_i X_i + \sum_{i=1}^{M} a_{ii} X_i^2 + \sum_{i=1}^{M} \sum_{j=1 \neq i}^{M} a_{ij} X_i X_j \quad [7.5] $$
où $\mathbf{a} = \{a_0, a_i, a_{ii}, a_{ij}\}^T$ est le vecteur des coefficients de la surface de réponse. Ces derniers sont déterminés par la méthode des moindres carrés à partir du plan d'expériences numériques choisi $\{\mathbf{x}^{(k)}, k=1,...,N\}$, où le nombre de points $N$ doit être au moins égal à la dimension de $\mathbf{a}$, soit :

$$ \mathbf{a} = \arg \min_{\mathbf{a}} \sum_{k=1}^{N} (y_k - \hat{g}(\mathbf{x}^{(k)}))^2 \quad [7.6] $$

où $y_k = g(\mathbf{x}^{(k)})$ est la valeur exacte de fonction d'état limite au point $\mathbf{x}^{(k)}$. En écrivant $\hat{g}(\mathbf{X}) = \{1, X_i, X_i^2, X_iX_j, \dots, a_0, a_i, a_{ii}, a_{ij}\}^T = \mathbf{B}(\mathbf{X})^T \mathbf{a}$, la détermination de $\mathbf{a}$ devient [FAV 89] :

$$ \mathbf{a} = (\mathbf{C}^T \mathbf{C})^{-1} \mathbf{C}^T \mathbf{y} \quad [7.7] $$

où $\mathbf{C}$ est la matrice dont les lignes sont les vecteurs $\mathbf{B}(\mathbf{x}^{(k)})^T$ et $\mathbf{y}$ le vecteur de composantes $y_k$.

La formulation générale de la méthode est simple. Elle soulève cependant des difficultés et des controverses lorsque l'on considère les points suivants, dont certains sont communs à tout type de surface de réponse, et d'autres plus particuliers aux surfaces polynomiales liées à des problèmes mécaniques :

- l'espace de travail : il plus aisé de contrôler le plan d'expériences numériques dans l'espace physique pour assurer une réponse cohérente du modèle pour chaque point du plan d'expériences. Mais il est plus commode de le construire dans l'espace standardisé. En effet, cet espace est adimensionnel, permet de s'assurer d'une distance entre les points, indispensable au bon conditionnement de $\mathbf{C}$, et facilite ainsi la liaison avec la détermination de l'indice de fiabilité ;

- le nombre de points du plan d'expériences : avec l'augmentation du nombre de points, la convergence de la régression est améliorée (statistique-F et variance des coefficients), mais pas nécessairement la qualité de l'approximation (coefficient de corrélation ajusté $R^2$ ou de validation croisée $Q^2$), laquelle est toujours satisfaisante avec un nombre de points égal au nombre de coefficients à calculer (car dans ce cas une interpolation est réalisée) ; la validité de la surface de réponse en dehors du plan d'expérience est souvent meilleure si elle est obtenue avec un grand nombre de points ; si la méthode de surface de réponse vise à réduire le coût par rapport à l'utilisation du modèle physique exact, le nombre de points doit être minimisé ;

- la topologie du plan d'expériences : un plan d'expériences compact permet un ajustement de meilleure qualité, mais concentre le domaine de validité de l'approximation ; le choix des directions de l'espace prospectées dépend idéalement
du comportement de la surface de réponse, qui n'est pas connue, et du modèle exact, dont on cherche à limiter le nombre de réponses à calculer ;

– la validité de la surface de réponse : en fonction de l'objectif fiabiliste (indice de fiabilité et/ou probabilité de défaillance estimée par simulations de Monte-Carlo), il peut être nécessaire de valider localement ou globalement la surface de réponse ;

– l'adaptabilité de la surface de réponse : la construction de la surface de réponse est rarement unique, en particulier si l'on cherche à la rendre performante pour un sous-domaine inconnu au départ (autour du point de défaillance le plus probable par exemple), on doit donc envisager une construction séquentielle couplée à un algorithme probabiliste, avec éventuellement une différentiation par sous-domaine ;

– le degré du polynôme doit être inférieur ou égal au degré inconnu de la fonction d'état limite, afin d'obtenir un système linéaire plus facile à résoudre ;

– la forme du polynôme est relative à la présence ou non de termes mixtes, qui contribuent à la prise en compte des interactions entre variables.

### 7.5.2. Choix de l'espace de travail

Dans le contexte de la fiabilité des structures deux espaces de travail sont possibles : l'espace physique X (variables aléatoires dimensionnelles du problème physique, de lois quelconques éventuellement corrélées), et l'espace standardisé U (variables aléatoires adimensionnelles gaussiennes centrées réduites et indépendantes). Le passage d'un espace à l'autre, via par exemple les transformations de Nataf, n'est requis pour la construction de la surface de réponse que lorsque celle-ci est opérée dans l'espace standardisé, car le modèle mécanique ne répond que dans l'espace physique. Ces transformations sont non linéaires, sauf pour des variables gaussiennes, et modifient la topologie du plan d'expérience ainsi que celle de la surface d'état limite [DEV 97] comme illustré sur la figure 7.3 où le plan d'expériences est composite centré dans l'espace standardisé. De nombreux auteurs ont ainsi choisi d'opérer dans l'espace physique [BUC 90, RAJ 93, MUZ 93, KIM 97, KAY 04, GAV 08]. Dans le cas où la détermination de l'indice de fiabilité guide de façon séquentielle la construction de la surface de réponse, il est préférable d'opérer directement dans l'espace standardisé, où l'adimensionnalité peut faciliter les aspects numériques [DEV 97, ENE 94, GAY 03, GUP 04, DUP 06, NGU 09].

Travailler en aveugle dans l'espace standardisé soulève la question de savoir si le modèle mécanique est en mesure de fournir une réponse pour la totalité des points du plan d'expériences, et impose de contrôler la validité des points retenus [DEV 97]. Néanmoins, si le choix des lois de distribution des variables aléatoires du modèle mécanique tient compte simultanément de l'information statistique disponible et des restrictions imposées par le modèle, quant au domaine que peuvent couvrir les variables, le passage de l'espace standardisé à l'espace physique ne doit
pas conduire à l'impossibilité pour le modèle de fournir une réponse. Il n'y a donc en principe pas d'inconvénients à opérer dans l'un ou l'autre des espaces.

<figure><img src="image_4.png" /><figcaption>Figure 7.3. Incidence du changement d'espace</figcaption></figure>

### 7.5.3. Expression de la surface de réponse

Quel que soit l'espace de travail retenu, l'objectif poursuivi est de pouvoir obtenir une estimation satisfaisante de la fiabilité mécanique, si possible à un coût numérique moindre que celui requis par l'emploi du modèle mécanique exact.

Dans cette optique, il apparaît naturel, si l'on se contente des estimations courantes FORM/SORM (paragraphe 7.4.2), d'orienter son choix vers des surfaces polynomiales linéaires ou quadratiques, nécessairement numériquement moins gourmandes que les surfaces d'ordre plus élevé. On peut également conjecturer que le degré de non linéarité de la surface d'état limite au voisinage du point de défaillance le plus probable n'est pas supérieur à trois la plupart du temps, ce qui permet d'espérer une bonne représentativité de la surface de réponse dans cette région. L'expression de la surface de réponse prend alors l'une des deux formes suivantes :

$$ 
\begin{aligned}
\hat{g}_I(\mathbf{X}) &= a_0 + \sum_{i=1}^{M} a_i X_i \\
\hat{g}_{II}(\mathbf{X}) &= a_0 + \sum_{i=1}^{M} a_i X_i + \sum_{i=1}^{M} a_{ii} X_i^2 + \sum_{i=1}^{M} \sum_{j=1 \neq i}^{M} a_{ij} X_i X_j
\end{aligned}
\qquad [7.8]
$$

Néanmoins, lorsque l'on projette d'affiner l'estimation de la probabilité de défaillance par l'emploi de simulations de Monte-Carlo, avec un domaine de prospection plus ou moins large autour du point de défaillance le plus probable, ou bien lorsque l'hypothèse de faible non-linéarité de la surface d'état limite n'est plus
acceptable, il peut devenir opportun d'augmenter le degré du polynôme [GUP 04, GAV 08]. L'expression de la surface de réponse prend alors la forme suivante :

$$ \hat{g}_p(\mathbf{X}) = a_0 + \sum_{i=1}^{M} \sum_{j=1}^{p_i} a_{ij} X_i^j + \sum_{q=1}^{M_m} b_q \prod_{i=1}^{M} X_i^{p_{iq}} \quad [7.9] $$

où les coefficients $a_{ij}$ portent sur la variable $X_i$ élevée à la puissance $j$, et les coefficients $b_q$ portent sur les $M_m$ produit de deux variables ou plus, chacune pouvant être élevée à une puissance différente $p_{iq}$ dans ce produit.

Gupta et Manohar [GUP 04] ne considèrent qu'une seule valeur $p_i=p$ et utilisent une expression quadratique sans terme mixte ($M_m=0$). S'appuyant sur la technique développée par Bucher & Bourgund [BUC 90], ils localisent approximativement le point de défaillance le plus probable $U_1^*$ et définissent une région de prospection de rayon $\beta_2 = -\Phi^{-1}(10^{-4}\Phi(-\beta_1))$ où $\beta_1=\|U_1^*\|$. D'autres points situés au voisinage de la surface d'état limite sont ensuite localisés de la même façon, chacun dans un nouveau repère déduit du repère initial par translation le long d'un des axes sur une distance égale à $(2\beta_2/p)$. Une collection d'au moins $N=(pM+1)$ points situés approximativement sur la surface d'état limite doit ainsi être constituée, nécessitant $(4M+3)(kM+1)$ valeurs de la fonction d'état limite exacte. Les auteurs ne précisent pas comment l'ordre $p$ est choisi (deux exemples sont traités avec $p=3$ puis $p=5$).

Gavin et Yau [GAV 08] ont développé une approche en deux étapes. La première étape consiste à déterminer les exposants $p_i$ axe par axe, en ignorant les termes mixtes ($M_m=0$). Pour ce faire une surface de réponse instrumentale est définie pour chaque axe à partir de monômes de Tchebychev :

$$ \hat{g}_i(\mathbf{X}_i) = \sum_{j=0}^{p} d_j T_j(\lambda_i) \quad [7.10] $$

où $T_j(\lambda_i)$ est le monôme de Tchebychev d'ordre $j$ en $\lambda_i$, définie sur [-1,1] :

$$ T_j(\lambda_i) = \cos(j \text{Arc cos}(\lambda_i)) \quad [7.11] $$

Les monômes de Tchebychev forment une base orthogonale et varient sur [-1,1], ce qui permet d'une part de calculer de façon unique la variance $\sigma_d^2$ des coefficients $d_j$ (pour $j > 0$), et d'autre part d'appliquer une statistique-T sur la variable $t_j=d_j/\sigma_d$. Les coefficients $d_j$ sont déterminés par la méthode des moindres carrés. Le plan d'expériences est constitué de $p$ points dont toutes les coordonnées sont fixées aux valeurs moyennes $x_{i \neq i}^{(j)} = \mu_i$, sauf pour l'axe $i$ pour lequel $x_i^{(j)} = \mu_i + h_{ord} \lambda_i^{(j)} \sigma_i$ où
$\lambda^{(j)}$, est la $j$ème racine du monôme d'ordre $j$, $\mu_i$ et $\sigma_i$ sont la moyenne et l'écart type de $X_i$. La valeur de l'ordre $p_i$ à retenir dans l'expression [7.9] est telle que l'hypothèse de nullité de $t_{pi}$ n'est pas rejetée au risque 90 % ($t_{pi} \le 3,5$), ou bien que le coefficient $d_{pi}$ est inférieur au centième de la plus grande valeur $max(d_j)$. Si la non-linéarité de la fonction d'état limite exacte est modifiée au voisinage du point moyen, l'ordre $p_i$ retenu pour chaque axe peut dépendre du coefficient $h_{ord}$, en fonction des changements de courbure rencontrés au fur et à mesure de l'extension de la région prospectée (valeurs de $h_{ord}$ croissantes). Dans les exemples traités [GAV 08], il apparaît que les ordres $p_i$ changent lorsque $h_{ord} \ge 3$, pouvant passer pour certaines variables de 2 ou 3 à 4 ou 5. Le nombre minimum de points nécessaires à la détermination des ordres $p_i$ est $(Mp) \ge (M \max(p_i))$.

La seconde étape consiste à recourir à deux règles simples pour déterminer les exposants $p_{iq}$, soit $p_{iq} \le p_i$ et $\sum_{i=1}^{M} p_{iq} \le \max(p_i)$. Le nombre de termes mixtes devient rapidement important dès que $\max(p_i) > 2$. Pour l'exemple précédent à 7 variables aléatoires, le nombre total de coefficients est 288 pour $h_{ord}=5$, dont 274 sur les termes mixtes.

Dans un exemple de comportement dynamique de structure [GUP 04], cas où la fonction d'état limite est fortement non linéaire, un ordre $p=3$ est retenu sans qu'une comparaison soit établie avec un ordre plus faible permettant de confirmer la nécessité du recours à une surface de réponse cubique. Pour des exemples similaires au précédent, traités dans [GAV 08], peu de variables interviennent avec un ordre $p_i > 2$ dans l'expression de la surface de réponse. Cela confirme l'idée que les surfaces de réponses quadratiques sont le plus souvent satisfaisantes, moyennant quelques précautions, telles que précisées ci-après.

### 7.5.4. Construction du plan d'expériences

La construction du plan d'expériences est un point crucial dans la méthode des surfaces de réponse : elle conditionne la qualité de la surface de réponse obtenue et agrège l'essentiel des efforts à produire, tout au moins dans le cas de surfaces de réponse linéaires ou quadratiques.

#### 7.5.4.1. Nombre, distribution et espacement des points

Le nombre minimal de points (nombre de coefficients à déterminer) est $N_c=(M+1)$ dans le cas d'une surface linéaire, $N_c=(2M+1)$ dans le cas d'une surface quadratique sans terme mixte et $N_c=(M+1)(M+2)/2$ dans le cas d'une surface quadratique avec termes mixtes. Avec ce nombre minimal, seule une interpolation peut être obtenue et la signification statistique des coefficients est alors nulle.
D'autre part, l'utilisation de la surface de réponse en dehors du plan d'expériences utilisé pour la construire est délicat. Même si l'économie numérique fait partie des objectifs, on peut rarement se contenter du nombre minimal de points, car la qualité de la surface de réponse s'en ressent.

Dans la phase initiale de construction de la surface de réponse, si plusieurs étapes sont envisagées, on opte généralement pour une distribution uniforme des points autour d'un point central (plan dit centré), comprenant au moins un point par axe, et possiblement des points en dehors des axes (plan dit composite centré). Le plan d'expériences de la figure 7.3 (dans l'espace standardisé) comprend $(2M+1+2^M)$ points ou $(3^M)$ points selon que l'on dispose un ou plusieurs points par quadrant. Les distributions non uniformes de points relèvent de procédures adaptatives où l'on cherche à obtenir de l'information dans la région de l'espace où devrait se situer la surface d'état limite (voir paragraphe 7.4.2).

La distance entre points du plan d'expérience est exprimée le long d'un axe i donné la plupart du temps en nombre d'écarts types :

$$x_i^{(k)} = X_{C,i} \pm h_i \sigma_i \quad [7.12]$$

où $x^{(k)}$ est le kème point courant du plan d'expériences et $X_C$ en est le point central. Similairement, si l'on opère dans l'espace standardisé :

$$u_i^{(k)} = U_{C,i} \pm h_i \quad [7.13]$$

Hormis dans le cas de procédures adaptatives, la valeur de $h_i$ est considérée constante ($h_i=h$). Des valeurs trop faibles ou trop grandes de $h$ compromettent la qualité de la surface de réponse quadratique, et conduisent à des résultats erronés. Le choix d'une valeur optimale de $h$ n'est pas immédiat [GUA 01]. Pour un exemple à 21 variables aléatoires [GUA 01] (état limite mécanique sur un portique multiple), la probabilité de défaillance passe ainsi de $10^{-4}$ à 0,94 lorsque $h$ varie de 0,1 à 5. Sur ce même exemple, la probabilité de défaillance et l'indice de fiabilité varient respectivement de $9,3 \cdot 10^{-5}$ à 0,37, et de 0,06 à 4,47 lorsque $h$ varie de 1,8 à 2,2.

La valeur de $h$ doit également être suffisamment large pour éviter un mauvais conditionnement du système $C$ [7.7]. Pour ce faire, on peut proposer des valeurs minimales pour $h$ variant de 0,2 à 0,5 [ENE 94, DUP 06], ou bien contrôler la proximité des points en calculant leur distance relative [DEV 97]. Divers auteurs fixent une valeur de $h$ dans la plage 1 à 3 [BUC 90, MUZ 93, ENE 94, DEV 97, KIM 97, KAY 04, WON 05, DUP 06, NGU 09].
### 7.5.4.2. Procédures adaptatives

Les procédures adaptatives de construction de la surface de réponse ont pour objet de combiner une approximation satisfaisante de la fonction d'état limite exacte au voisinage du point de défaillance le plus probable, avec une dépense numérique aussi limitée que possible. L'adaptation comporte plusieurs étapes portant sur les éléments qui suivent :

– la forme de la surface de réponse : il peut être intéressant de modifier une première surface de réponse linéaire ($\hat{g}_I(X)$) ou quadratique sans terme mixte ($\hat{g}_{II}(X)$ avec $a_{ij}=0$), permettant une estimation peu coûteuse du point de défaillance le plus probable, par une seconde surface, quadratique avec termes mixtes [GAY 03, NGU 09] ;

– le plan d'expériences initial : il est généralement centré. Si le rôle des variables aléatoires dans la fonction d'état limite est connu a priori, un point de départ $X_{C1}$ « meilleur candidat » au point de défaillance le plus probable, peut être localisé plus proche de la région de défaillance, en choisissant $X_{C1,i} = \mu_i - h_i\sigma_i$ si $X_i$ est une variable favorable à la performance, $X_{C1,i} = \mu_i + h_i\sigma_i$ si $X_i$ est une variable favorable à la défaillance. Si le rôle des variables est inconnu, on peut opter pour une phase de pré-prospection nécessitant quelques appels à la fonction d'état limite exacte [GAY 03] ;

– le second point central : il peut être fixé connaissant le point de défaillance le plus probable obtenu avec la première surface de réponse, à partir d'une interpolation [BUC 90], ou bien en utilisant directement ce dernier, moyennant son appartenance au premier plan d'expériences [ENE 94, DEV 97] ;

– le maillage du plan d'expériences : sa dimension est souvent réduite au fur et à mesure de la progression vers le point de défaillance le plus probable [MUZ 93, ENE 94, DEV 97, GAY 03, KAY 04]. Il peut également s'avérer intéressant de tenir compte du nombre de variables aléatoires dans cette dimension [KIM 97], ainsi que de la sensibilité de la surface de réponse aux différentes variables aléatoires [DUP 06, NGU 09]. Enfin, la dimension du maillage peut être guidée par des considérations statistiques sur les coordonnées du point de défaillance le plus probable [GAY 03] ;

– la pondération des points du plan d'expériences : les coefficients de la surface de réponse sont calculés en introduisant dans [7.7] une matrice de pondération favorisant les points proches de la surface d'état limite [KAY 04], et simultanément ceux proches du dernier point de défaillance le plus probable estimé [NGU 09] ;

– l'enrichissement du plan d'expériences : le plus grand nombre possible de points où la fonction d'état limite a été calculée aux précédentes itérations est conservé, cela permet de tirer parti du coût numérique consenti. Mais quelques précautions sont à prendre pour éviter de conserver les points risquant de nuire à la qualité de l'approximation [ENE 94, DEV 97, GAY 03, DUP 06, NGU 09].
### 7.5.4.3. Qualité de l'approximation obtenue

L'indicateur le plus commode pour juger de la qualité de l'approximation reste le coefficient de corrélation ajusté qui peut s'écrire, dans le cas d'un point scalaire $x^{(k)}$ :

$$R^2 = 1 - \frac{(N-1) \sum_{k=1}^{N} (\hat{g}(x^{(k)}) - y_k)^2}{(N - N_c - 1) \sum_{k=1}^{N} (y_{\text{mean}} - y_k)^2} \quad [7.14]$$

où $y_{\text{mean}}$ est la valeur moyenne de la fonction d'état limite sur l'ensemble du plan d'expériences. Si $R^2$ est supérieur à 0,9, on peut souvent considérer que l'approximation est satisfaisante au sein du plan d'expériences. Il est possible également d'utiliser un coefficient de corrélation croisée (voir le chapitre 8).

## 7.5.5. Exemple d'une méthode adaptative

### 7.5.5.1. Description d'ensemble de la méthode [NGU 09]

#### 7.5.5.1.1. Première itération

La première surface de réponse $\hat{g}_I(U)$ est linéaire dans le repère standardisé. Le plan d'expériences est centré sur $U_C$ et comprend $U_C$ et un point par axe défini par :

$$u_i^{(k)} = U_{C,i} + h_i \quad \text{avec} \quad h_i = - \frac{h_0}{\| \nabla \hat{g}(U_C) \|} \frac{\partial \hat{g}(U_C)}{\partial U_i} \quad [7.15]$$

où $u_i^{(k)}$ est le $k$-ème point courant du plan d'expériences. Pour le premier plan d'expériences, $U_{C,i}=U_0$ étant l'origine du repère, un gradient fictif de $\hat{g}_I(U_{C,i})$ est proposé en fonction de la connaissance *a priori* du rôle joué par les variables, de telle sorte d'avoir $h_i=\pm 1$. Pour favoriser les points situés plus près de la surface d'état limite exacte, une matrice de pondération diagonale $W$ est introduite dans l'équation [7.7], selon la proposition de Kaymaz & McMahon [KAY 04], qui devient :

$$a = (C^T W C)^{-1} C^T W z \quad [7.16]$$

où les coefficients de $W$ sont calculés par :

$$W_k = \exp \left( - \frac{g(u^{(k)}) - y_{\min}}{g(U_0)} \right) \quad [7.17]$$
où $z_{min}$ est la valeur minimale de la fonction d'état limite pour tous les points du plan d'expériences. Les coefficients de $\hat{g}_I(U)$ calculés par [1.16], le premier point de défaillance le plus probable $U^{*(1)}$ est déterminé selon une approche FORM.

### 7.5.5.1.2. Deuxième itération

A partir de cette itération, la forme de la surface de réponse est la forme quadratique $\hat{g}_{II}(U)$ à $(M+1)(M+2)/2$ coefficients inconnus. Le point central du second plan d'expériences est le point $U_{C2}$ défini, selon la proposition de Bucher & Bourgund [BUC 90], par :

$$U_{C2} = U_0 + (U^{*(1)} - U_0) \frac{g(U_0)}{g(U_0) - g(U^{*(1)})} \quad [7.18]$$

situé plus près de la surface d'état limite exacte que le point $U^{*(1)}$, pour lequel en général $g(U^{*(1)}) \neq 0$.

Tous les points du premier plan d'expériences sont conservés dans le second, auquel on adjoint $((M+1)(M+2)/2 - (M+1))$ points complémentaires disposés en demi-étoile autour de $U_{C2}$. Parmi les points complémentaires, dont $U_{C2}$ fait partie, M points sont placés sur les axes selon [7.15]. Ce choix assure que les points de calage sont situés vers la région de défaillance (par rapport au point central) et à une distance proportionnelle à la sensibilité locale du modèle par rapport à la variable.

Les points complémentaires restant sont engendrés à partir des points situés sur les axes, comme les M premiers l'ont été à partir de $U_{C2}$. Les axes sont parcourus dans l'ordre décroissant des valeurs absolues des composantes du gradient de $\hat{g}_I(U)$ calculées en $U_{C2}$, dans le but d'améliorer l'approximation selon les axes de plus forte sensibilité. Le nombre maximal de points ainsi engendrés vaut $\sum_{i=1}^{M-1} (M-i)$, ce qui est toujours supérieur au nombre minimal requis. Le point $U^{*(1)}$ est conservé si la condition suivante est respectée :

$$\|U_{C2} - U^{*(1)}\| \le h_0 \quad [7.19]$$

Les coefficients de la matrice de pondération W sont à présent modifiés pour tenir compte également de la distance entre les points du plan d'expériences et le dernier point de défaillance le plus probable estimé, soit :
$$W_k = \exp\left(-\frac{g(\mathbf{u}^{(k)}) - y_{\min}}{g(\mathbf{U}_0)}\right) \exp\left(-\frac{\|\mathbf{U}^{*(\text{iter}-1)} - \mathbf{u}^{(k)}\|^2}{2}\right) \quad [7.20]$$

Dans le cas de la deuxième itération il s'agit du point $\mathbf{U}^{*(1)}$ si la condition [7.19] est observée, sinon du point $\mathbf{U}_{C2}$. Une fois calculés les coefficients de $\hat{g}_{II}(\mathbf{U}, \mathbf{d})$ par [7.20], le second point de défaillance le plus probable $\mathbf{U}^{*(2)}$ est déterminé.

### 7.5.5.1.3. Itérations suivantes (iter>2)

Le plan d'expériences est enrichi du point $\mathbf{U}^{*(\text{iter}-1)}$. Les facteurs de pondérations sont mis à jour selon [7.20] puis les coefficients de $\hat{g}_{II}(\mathbf{U})$ et le nouveau point de défaillance le plus probable sont déterminés. Les itérations sont stoppées lorsque les critères de convergence sont satisfaits :

$$\left\{ \begin{array}{l} \left| \beta^{(\text{iter})} - \beta^{(\text{iter}-1)} \right| \le \varepsilon_{\beta} \\ \left\| \mathbf{U}^{*(\text{iter})} - \mathbf{U}^{*(\text{iter}-1)} \right\| \le \varepsilon_{P} \end{array} \right. \qquad [7.21]$$

### 7.5.5.2. Exemples d'application

Trois exemples sont présentés permettant de montrer l'intérêt mais aussi la limite des méthodes de surfaces polynomiales. La méthode décrite ci-avant est dénommée SRDP (surface de réponse à double pondération). Le premier exemple traite d'une fonction d'état limite simple non linéaire (exemple 1) :

$$g(\mathbf{U}) = \exp(0,4(\mathbf{U}_2 + 2) + 6,2) - \exp(0,3\mathbf{U}_1 + 5) - 200 \quad [7.22]$$

Toutes les variables aléatoires ont une distribution normale standard. Les résultats sont donnés dans le tableau 7.1. On constate que les valeurs de l'indice de fiabilité sont très proches les unes des autres pour toutes les méthodes. De même, les rapports $|g(\mathbf{U}^*) / g(\mathbf{U}_0)|$ sont très voisins de zéro et traduisent la proximité du dernier point de défaillance le plus probable estimé avec la surface d'état limite. La formation progressive de la surface de réponse est illustrée sur la figure 7.4. Enfin, les tableaux 7.2 et 7.3 traduisent l'influence de la valeur retenue pour $h_0$ sur la qualité et l'efficacité de la surface de réponse. On voit que l'utilisation d'un système de pondération des points du plan d'expériences permet de réduire cette influence et d'obtenir une meilleure proximité avec la surface d'état limite exacte.
Puis étudions une fonction d'état limite non linéaire plus complexe (exemple 2) :

$$g(\mathbf{X}) = X_1 - 10^4 \left[ \frac{X_2 (X_4 X_5)^{1,71}}{X_3} + \frac{(1-X_2) (X_4 X_5)^{1,188}}{X_6} \right] \quad [7.23]$$

Les lois de probabilité des variables aléatoires sont portées dans le tableau 7.4. L'indice de fiabilité obtenu est $\beta = 2,387$ nécessitant 6 itérations et 32 calculs de la fonction d'état limite exacte. La proximité à la surface d'état limite du dernier point de défaillance le plus probable est donnée par $|g(U^*) / g(U_0)| = 1,09 \times 10^{-4}$.

La qualité de la surface de réponse est satisfaisante ($R^2 = 0,999$). On constate encore sur cet exemple que l'influence de la valeur retenue pour $h_0$ est réduite lorsqu'un système de pondération est utilisé (tableaux 7.5 et 7.6).

Le troisième exemple (exemple 3) traite d'une fonction d'état limite non explicite où intervient le résultat d'un calcul aux éléments finis (ossature de bâtiment). Il s'agit du même exemple que celui traité dans le chapitre 8 (paragraphe 8.4.2), où l'on peut trouver le détail des lois de probabilité des 21 variables aléatoires corrélées. Les résultats sont présentés dans le tableau 7.7, où plusieurs valeurs du seuil en déplacement $u_{max}$ sont considérées.

Les résultats en termes de convergence et de proximité avec la surface d'état limite exacte sont satisfaisants. La comparaison avec le calcul de l'indice d'Hasofer-Lind obtenu par couplage direct selon une méthode FORM, montre d'une part que la SRDP produit une valeur correcte de l'indice de fiabilité, mais d'autre part que la dépense numérique consentie est beaucoup plus importante. La méthode FORM est donc ici plus économique.

<table><caption>**Tableau 7.1.** Comparaison de plusieurs méthodes de SR – Exemple 1</caption><thead><tr><th>Méthode</th><th>β</th><th>U*<sub>1</sub></th><th>U*<sub>2</sub></th><th>|g(U*)/g(U<sub>0</sub>)|</th><th>N<sub>c</sub></th></tr></thead><tbody><tr><td>MC adaptatif [KAY 04]</td><td>2,710</td><td>0,969</td><td>-2,531</td><td>3,36×10<sup>-5</sup></td><td>-</td></tr><tr><td>SR [KIM 97]</td><td>2,691</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SR [KAY 04]</td><td>2,686</td><td>0,820</td><td>-2,558</td><td>5,84×10<sup>-3</sup></td><td>8</td></tr><tr><td>SR [DUP 06]</td><td>2,710</td><td>0,951</td><td>-2,538</td><td>9,10×10<sup>-4</sup></td><td>21</td></tr><tr><td>SRDP (R<sup>2</sup> = 0,997)</td><td>2,707</td><td>0,860</td><td>-2,567</td><td>8,48×10<sup>-4</sup></td><td>12</td></tr></tbody></table>
<table><caption>**Tableau 7.2.** Influence de $h_0$ avec système de pondération ($W_k \neq 1$) – Exemple 1</caption><tr><td>h<sub>0</sub></td><td>1</td><td>2</td><td>3</td><td>4</td></tr><tr><td>β</td><td>2,707</td><td>2,724</td><td>2,713</td><td>2,714</td></tr><tr><td>|g(U*)/g(U<sub>0</sub>)|</td><td>8,48×10<sup>-4</sup></td><td>4,06×10<sup>-4</sup></td><td>8,02×10<sup>-4</sup></td><td>6,88×10<sup>-4</sup></td></tr><tr><td>R<sup>2</sup></td><td>0,997</td><td>0,999</td><td>0,999</td><td>0,930</td></tr><tr><td>N<sub>c</sub></td><td>12</td><td>10</td><td>11</td><td>13</td></tr></table>

<table><caption>**Tableau 7.3.** Influence de $h_0$ sans système de pondération ($W_k=1$) – Exemple 1</caption><tr><td>h<sub>0</sub></td><td>1</td><td>2</td><td>3</td><td>4</td></tr><tr><td>β</td><td>2,707</td><td>2,715</td><td>2,715</td><td>2,726</td></tr><tr><td>|g(U*)/g(U<sub>0</sub>)|</td><td>1,15×10<sup>-3</sup></td><td>1,45×10<sup>-3</sup></td><td>1,12×10<sup>-3</sup></td><td>3,16×10<sup>-3</sup></td></tr><tr><td>R<sup>2</sup></td><td>0,997</td><td>0,999</td><td>0,999</td><td>0,885</td></tr><tr><td>N<sub>c</sub></td><td>14</td><td>11</td><td>11</td><td>15</td></tr></table>

<table><caption>**Tableau 7.4.** Distribution des variables aléatoires – Exemple 2</caption><thead><tr><td>Variable</td><td>Distribution</td><td>Moyenne</td><td>Ecart type</td></tr></thead><tbody><tr><td>X<sub>1</sub></td><td>Lognormale</td><td>1,044</td><td>0,313</td></tr><tr><td>X<sub>2</sub></td><td>Normale</td><td>0,7</td><td>0,07</td></tr><tr><td>X<sub>3</sub></td><td>Lognormale</td><td>0,239</td><td>0,0956</td></tr><tr><td>X<sub>4</sub></td><td>Lognormale</td><td>1,011</td><td>0,1516</td></tr><tr><td>X<sub>5</sub></td><td>Gumbel Max</td><td>0,0005</td><td>0,00008</td></tr><tr><td>X<sub>6</sub></td><td>Lognormale</td><td>1,802</td><td>0,7208</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Figure 7.4. Formation cumulative de la surface de réponse – Exemple 1</figcaption></figure>

<table><caption>Tableau 7.5. Influence de h₀ avec système de pondération (W<sub>k</sub>≠1) – Exemple 2</caption><tr><td>h<sub>0</sub></td><td>1</td><td>2</td><td>3</td><td>4</td></tr><tr><td>β</td><td>2,387</td><td>2,386</td><td>2,392</td><td>2,389</td></tr><tr><td>|g(U*)/g(U<sub>0</sub>)|</td><td>1,09x10<sup>-4</sup></td><td>9,32x10<sup>-4</sup></td><td>4,05x10<sup>-4</sup></td><td>8,59x10<sup>-4</sup></td></tr><tr><td>R<sup>2</sup></td><td>0,999</td><td>0,999</td><td>0,968</td><td>0,922</td></tr><tr><td>N<sub>c</sub></td><td>32</td><td>31</td><td>35</td><td>15</td></tr></table>

<table><caption>Tableau 7.6. Influence de h₀ sans système de pondération (W<sub>k</sub>=1) – Exemple 2</caption><tr><td>h<sub>0</sub></td><td>1</td><td>2</td><td>3</td><td>4</td></tr><tr><td>β</td><td>2,389</td><td>2,387</td><td>2,387</td><td>2,433</td></tr><tr><td>|g(U*)/g(U<sub>0</sub>)|</td><td>4,88x10<sup>-4</sup></td><td>4,68x10<sup>-4</sup></td><td>1,15x10<sup>-3</sup></td><td>9,10x10<sup>-3</sup></td></tr><tr><td>R<sup>2</sup></td><td>0,997</td><td>0,999</td><td>0,893</td><td>0,762</td></tr><tr><td>N<sub>c</sub></td><td>32</td><td>32</td><td>40</td><td>43</td></tr></table>
<table><caption>Tableau 7.7. Résultats obtenus pour l'exemple 3</caption><thead><tr><th colspan="2">u<sub>maxn</sub> (cm)</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td rowspan="2">β</td><td>SRDP</td><td>2,19</td><td>2,86</td><td>3,28</td></tr><tr><td>FORM</td><td>2,20</td><td>2,86</td><td>3,40</td></tr><tr><td rowspan="2">|g(U*)/g(U<sub>0</sub>)|</td><td>SRDP</td><td>5,39x10<sup>-3</sup></td><td>2,80x10<sup>-3</sup></td><td>1,06x10<sup>-3</sup></td></tr><tr><td>FORM</td><td>6,03x10<sup>-6</sup></td><td>6,25x10<sup>-4</sup></td><td>1,72x10<sup>-5</sup></td></tr><tr><td rowspan="2">N<sub>c</sub></td><td>SRDP</td><td>256</td><td>258</td><td>255</td></tr><tr><td>FORM</td><td>80</td><td>75</td><td>104</td></tr><tr><td colspan="2">R<sup>2</sup></td><td>0,993</td><td>0,994</td><td>0,993</td></tr></tbody></table>

## 7.6. Conclusion

La méthode des surfaces de réponse a été utilisée en premier lieu dans des domaines où la réponse d'un système à des stimuli ne trouvait pas de modèles explicatifs satisfaisants. Elle a été portée vers la fiabilité mécanique avec l'idée de réduire les coûts numériques liés aux modèles de calcul mécaniques, comme les modèles aux éléments finis. Les exemples présentés dans ce chapitre montrent que cet objectif n'est pas toujours atteint, même dans l'emploi d'une méthode de surface de réponse polynomiale quadratique adaptative, dès lors que le nombre de variables aléatoires dépasse la vingtaine. La représentativité de ce type de surface de réponse est cependant satisfaisante : les résultats de l'analyse de fiabilité restent très proches de ceux obtenus par l'emploi de la fonction d'état limite exacte, lorsque la topologie du domaine de défaillance est favorable. Par ailleurs, disposer d'une fonction d'état limite explicite définie au voisinage du point de défaillance le plus probable facilite tous les post-traitements au calcul de l'indice de fiabilité : estimation de la probabilité de défaillance, calcul des sensibilités, etc. L'utilisation d'une surface de réponse quadratique paraît donc appropriée pour des problèmes à faible nombre de variables aléatoires ($M<20$).

## 7.7. Bibliographie

[ATK 68] ATKINSON A.C., HUNTER W.G., « The design of experiments for parameter estimation », Technometrics, vol. 10, p. 271-289, 1968.

[ANS 63] ANSCOMBE F.J., TUKEY F.W., « The examination and analysis of residuals », Technometrics, vol. 5, p. 141-160, 1963.
[BIE 83] BIER V.M., « A measure of uncertainty importance for components in fault trees », *ANS Trans.*, n° 45(1), p. 384-385, 1983.

[BLI 35a] BLISS C.I., « The calculation of the dosage-mortality curve », *Ann. Appl. Biol.*, vol. 22, p. 134-167, 1935.

[BLI 35b] BLISS C.I., « The comparison of the dosage-mortality data », *Ann. Appl. Biol.*, vol. 22, p. 307-333, 1935.

[BOU 95] BOUYSSY V., RACKWITZ R., « Approximation of non-normal responses for drag dominated offshore structures », *Proc. 6th IFIP WG 7.5, Reliability and Optimization of Structural Systems*, n° 7, Chapman & Hall, p. 161-168, 1995.

[BOX 50] BOX G.E.P., WILSON K.B., « On the experimental attainment of optimum conditions », *J. R. Statist. Soc.*, B13, p. 1-45, 1950.

[BOX 55] BOX G.E.P., YOLE P.V., « The exploration and exploitation of response surfaces : an example of the link between the fitted surface and the basic mechanism of the system », *Biometrics*, n° 11, p. 287-322, 1955.

[BOX 57] BOX G.E.P., HUNTER J.S., « Multifactor experimental designs for exploring response surfaces », *Ann. Math. Statist.*, n° 28, p. 195-241, 1957.

[BOX 59a] BOX G.E.P., DRAPER N.R., « A basis for the selection of a response surface design », *J. Amer. Statist. Ass.*, n° 54, p. 622-654, 1959.

[BOX 59b] BOX G.E.P., LUCAS H.L., « Design of experiments in nonlinear situations », *Biometrika*, n° 46, p. 77-90, 1959.

[BUC 90] BUCHER C.G., BOURGUND U., « A fast and efficient response surface approach for structural reliability problems », *Structural Safety*, vol. 7 (1), p. 57-66, 1990.

[CRO 35] CROWTHER E.M., YATES F., « Fertiliser policy in war-time », *Empire J. Exp. Agric.*, n° 9, p. 77-97, 1941.

[DEV 97] DEVICTOR N., MARQUES M., LEMAIRE M., « Utilisation des surfaces de réponse dans le calcul de la fiabilité des composants mécaniques », *Revue Française de Mécanique*, n° 1997-1, p. 43-51, 1997.

[DIT 96] DITLEVSEN O., MADSEN H., *Structural reliability methods*, John Wiley & Sons, Chichester, 1996.

[DUP 06] DUPRAT F., SELLIER A., « Probabilistic approach to corrosion risk due to carbonation *via* an adaptive response surface method », *Probabilistic Engineering Mechanics*, n° 21(4), p. 207-216, 2006.

[ELS 62] ELSTON R.C., GRIZZLE J.E., « Estimation of time response curves and their confidence bands », *Biometrics*, vol. 18, p. 148-159, 1962.
[ENE 94] ENEVOLDSEN I., FABER M.H., SORENSEN J.D., «Adaptive response surface technique in reliability estimation », dans G.I. Schuëller, M. Shinozuka et J.T.P. Yao (dir.), *Structural Safety and Reliability, Proceedings of the ICOSSAR'93*, Balkema Brookfield, Rotterdam, p. 1257-1264, 1994.

[FAV 89] FAVARELLI L., « Response surface approach for reliability analysis », *Journal of Engineering Mechanics*, vol. 115 (12), p. 2763-2781, 1989.

[GAD 33] GADDUM J.H., « Reports on biological Standards III : methods on biological assay depending on a quantal response », *Special Report n° 183*, Medical Research Council, H. M.S.O, Londres, 1933.

[GAV 08] GAVIN H.P., YAU S.C., « High-order limit state functions in the response surface method for structural reliability analysis », *Structural Safety*, vol. 30(2), p. 162-179, 2008.

[GAY 03] GAYTON N., BOURINET J.M., LEMAIRE M., « CQ2RS : a new statistical approach to the response surface method for reliability analysis », *Struct. Safety*, vol. 25(1), p. 99-121, 2003.

[GUA 01] GUAN X.L., MELCHERS R.E., « Effect of response surface parameter variation on structural reliability estimates », *Structural Safety*, vol. 23 (4), p. 429-444, 2001.

[GUP 04] GUPTA S., MANOHAR C.S., « An improved response surface method for the determination of failure probability and importance measures », *Structural Safety*, vol. 26 (2), p. 123-139, 2001.

[HEL 93] HELTON J.C., « Uncertainty and sensitivity analysis techniques for use in performance assessment for radioactive waste disposal », *Rel. Engrg. and Syst. Saf.*, vol. 42, p. 327-367, 1993.

[IMA 87] IMAN R.L., « A matrix-based approach to uncertainty and sensitivity analysis for fault trees », *Risk analysis n° 7* (1), p.21-33, 1987.

[JCS 02] JOINT COMMITTEE ON STRUCTURAL SAFETY, JCSS probabilistic model code, 2002, http://www.jcss.ethz.ch.

[KAY 04] KAYMAZ I., MCMAHON C., « A response surface method based on weighted regression for structural reliability analysis », *Probabilistic Engineering Mechanics*, vol. 20(1), p.1-7, 2004.

[KHA 89] KHATIB-RAHBAR M., « A probabilistic approach to quantifying uncertainties in the progression of severe accidents », *Nucl. Sci. Engrg.*, n° 10 (2), p. 219-259, 1989.

[KIE 52] KIEFER J., WOLFOWITZ J., « Stochastic estimation of the maximum of a regression function », *Ann. Math. Statist.*, n° 23, p. 462-466, 1952.

[KIE 60] KIEFER J., WOLFOWITZ J., « The equivalence of two extremum problems », *Canadian J. Math.*, n° 12, p. 363-366, 1960.

[KIM 97] KIM S.H., NA S.W., « Response surface method using vector projeted sampling points », *Structural Safety*, n° 19(1), p. 3-19, 1997.
[LAB 95] LABEYRIE J., SCHOEFS F., « A discussion on response surface approximations for use in structural reliability », 6th IFIP WG7.5, « Reliability and Optimisation of Structural Systems », n° 15, p. 161-168, Chapman and Hall, Londres, 1995.

[LAB 96] LABEYRIE J., SCHOEFS F., « Matrix response surfaces for describing environmental loads », vol. II, Safety and Reliability, *Proc. of 15th int. conf. on Offshore Mechanics and Arctic Engineering*, (O.M.A.E'96), p. 119-126, ASME, Florence, 1996.

[LEM 05] LEMAIRE M., *Fiabilité des structures – Couplage mécano-fiabiliste statique*, Hermès, Paris, 2005.

[MEA 75] MEAD R., PIKE D.J., « A review of response surface methodology from a biometric viewpoint », *Biometrics*, n° 31, p. 803-851, 1975.

[MIT 30] MITSCHERLICH E.A., *Die Bestimmung des Dungerbedurfnisses des Bodens*, Paul Paray, Berlin, 1930.

[MUZ 93] MUZEAU J.P., LEMAIRE M., EL-TAWIL K., « Méthodes fiabilistes des surfaces de réponse quadratiques et évaluation des règlements », *Construction Métallique*, vol. 3, p. 41-52, 1993.

[NAK 95] NAKAYASU H., « Relation between parameters sensitivities and dimensional invariance on Stochastic materials design of fibrous composite laminates », 6th IFIP WG 7.5, « Reliability and Optimisation of structural system », n° 21, p. 209-216, 1995.

[NEL 65] NELDER J.A., MEAD R., « A simplex method for function minimization », *Comp. J.*, n° 7, p. 308-313, 1965.

[NEL 66] NELDER J.A., « Inverse polynomials, a useful group of multi-factor response functions », *Biometrics* 22, p. 128-141, 1966.

[NGU 09] NGUYEN X.S., SELLIER A., DUPRAT F., PONS G., « Adaptive response surface method based on a double weighted regression technique », *Probabilistic Engineering Mechanics*, vol. 24 (2), p. 135-143, 2009.

[PAR 94] PARK C.K., ANH K.I., « A new approach for measuring uncertainty importance and distributional sensitivity in probabilistic safety assessment », *Rel. Engrg. and Syst. Saf.*, n° 46, p. 253-261, 1994.

[POW 65] POWEL M.J.D., « A method for minimizing a sum of squares of nonlinear functions without calculating derivatives », *Comp. J.*, n° 8, p. 303-307, 1965.

[RAJ 93] RAJASHEKHAR M.R., ELLINGWOOD B.R., « A new look at the response surface approach for reliability analysis », *Structural Safety*, vol. 12, p. 205-220, 1993.

[RAO 58] RAO C.R., « Some statistical methods for the comparison of growth curves », *Biometrics*, n° 14, p. 1-17, 1958.

[REE 29] REED L.J., BERKSON J., « The application of the logistic function for experimental data », *J. Phys. Chem.*, n° 33, p. 760-799, 1929.
[ROB 51] ROBINS H., MONRO S., « A stochastic approximation method », *Ann. Math. Statist.*, n° 22, p. 400-407, 1951.

[SCH 08] SCHOEFS F., « Sensitivity approach for modelling the environmental loading of marine structures through a matrix response surface », *Reliability Engineering and System Safety*, vol. 93, n° 7, juillet 2008, p. 1004-1017, 2008, doi: dx.doi.org/10.1016/j.ress.2007.05.006.

[SCH 07] SCHOEFS F., Méthodologie des surfaces de réponse physiques pour une analyse de risque intégrée des structures offshore fissurées, Mémoire d'habilitation à diriger les recherches, ED MTGC, Université de Nantes, 3 décembre 2007.

[SCH 96] SCHOEFS F., « Surface de réponse des efforts de houle dans le calcul de fiabilité des Ouvrages », Thèse de doctorat, ED 82-208, Université de Nantes, 21 novembre 1996.

[STE 51] STEVENS W.L., « Asymptotic regression », *Biometrics*, n° 7, p. 247-267, 1951.

[WIN 32] WINSOR C.P., « The Gompertz curve as a growth curve », *Proc. Natl. Acad. Sci.*, n° 18, p. 1-8, 1932.

[WIS 39] WISHART J., « Statistical treatment of animal experiment », *J. R. Statist. Soc.*, B6, p. 1-22, 1939.

[WON 05] WONG S.M., HOBBS R.E., ONOF C., « An adaptive response surface method for reliability analysis of structures with multiple loading sequences », *Structural Safety*, vol. 217 (3), p. 287-308, 2005.