FR 2933792 A1 20100115 FR 0803874 A 20080708 PROCEDE DE CONSTRUCTION D'UN METAMODELE DE SIMULATION DE DONNEES TECHNIQUES. L'invention proposée concerne l'aide à l'analyse de données techniques obtenues par un procédé de simulation technique ou scientifique ou d'ingénierie. Il peut donc concerner une très large gamme d'applications dans de nombreux secteurs techniques et scientifiques, parmi lesquels on peut citer à titre d'exemple la microélectronique et les micro- et nano-technologies (avec par exemple les logiciels de "TCAD" (Technology Computer Aided Design), les logiciels de "SPICE" (Simulation Program with Integrated Circuit Emphasis)), la construction automobile, aéronautique ou spatiale, le bâtiment et les travaux publics. la chimie, l'extraction minière, pétrolière et gazière, l'énergie, l'agro-alimentaire, le transport, la finance, le traitement numérique de données audio et/ou vidéo. Dans tous ces domaines, l'homme de l'art désire souvent comprendre ou visualiser un phénomène, ou optimiser une pièce, un dispositif ou 15 procédé, ou faire des études de sensibilité, de variabilité ou de la modélisation inverse, ou faire des études statistiques ou de rendement, ou créer des programmes qui simulent certains résultats dans une gamme fixée de paramètres. Pour parvenir à son but, l'homme de l'art a souvent à sa disposition des logiciels techniques ou scientifiques ou d'ingénierie qui permettent de simuler le 20 dispositif, phénomène ou concept étudié, et plus particulièrement des procédés qui se caractérisent par les faits suivants : l'homme de l'art décrit le dispositif ou phénomène ou concept étudié dans un langage ou avec des jeux de données et avec une méthodologie propre. Pour cela, il fait intervenir un certain nombre de "paramètres", c'est-à-dire de 25 grandeurs numériques lui permettant de décrire son dispositif (par exemple la géométrie du dispositif, les propriétés physiques des matériaux ou des caractéristiques des parties ou des sous-ensembles du dispositif) ou phénomène ou concept étudié, ou encore qui conditionnent la mise en oeuvre du procédé (par exemple le maillage, des paramètres ou choix algorithmiques ou numériques). Ces 30 "paramètres" sont choisis par l'homme de l'art dans une certaine gamme possible. • les résultats sont caractérisés par une ou (plus généralement) plusieurs "réponses" numériques qui renseignent l'homme de l'art sur le dispositif ou pièce ou procédé ou phénomène ou concept étudié. • le procédé est déterministe, c'est-à-dire que si l'homme de l'art 35 utilise plusieurs fois le même procédé avec les mêmes paramètres sur le même ordinateur, il aura à chaque fois les mêmes "réponses" numériques. très proche, c'est-à-dire que le métamodèle soit un "interpolateur" ou un "quasi-interpolateur". Un interpolateur est un métamodèle qui passe exactement par tous les points d'apprentissage, alors qu'un quasi-interpolateur est un métamodèle qui passe très près de tous les points d'apprentissage. Par "très près de tous les points de la base d'apprentissage", il faut comprendre que la variance résiduelle doit être petite par rapport à la variance totale, par exemple inférieure à 0,01 fois la variance totale ou, mieux, par exemple inférieure à 0,001 fois la variance totale ou, encore mieux, par exemple . nférieure à 0,0001 fois la variance totale. La variance résiduelle est ici définie (très classiquement) par : n Var résid = W • . ( R 2 simulation, j û R méta mod èle ( Xi %% j=l et la variance totale est ici définie (classiquement) par : n Var _ totale = 1 col . (Rsimulation, j û ) j=1 Les o sont des poids tous positifs ou nuls (par défaut et si il n'y a 15 pas de poids prédéfinis, tous les poids co sont pris égaux à 1). n est le nombre de simulations dans la base d'apprentissage. Xi est le vecteur représentant l'ensemble des paramètres pour la jème simulation, Rsi,,,,,iatio,,, i est la valeur de la réponse pour la jème simulation, R111èta,,,odèie(Xi) est le métamodèle associé à cette réponse. est la moyenne pondérée par les poids coi de toutes les simulations de la base 20 d'apprentissage. Plus le paramètre de lissage est grand, plus la variance résiduelle est petite. Le paramètre de lissage p est considéré comme "grand" si le métamodèle ainsi obtenu est un quasi interpolateur. p peut être par exemple choisi entre 103 et 1012 et plus particulièrement entre 104 et 1010. 25 Une première manière connue de construire ce métamodèle et de le chercher sous forme d'un polynôme dépendant des divers paramètres. Dans ce cas, le métamodèle est souvent nommé "surface de réponse", quoique cette appellation regroupe parfois, chez certains auteurs, un ensemble plus vaste de métamodèles. Par exemple, la réponse R représentée par un polynôme de degré total q dépendant de d 30 variables x1 , x2, ... , Xk, ... ,Xd s'écrit: 10 2 2 Les procédés connus de simulation technique ou scientifique ou d'ingénierie sont souvent assez lourds, c'est-à-dire qu'ils demandent beaucoup de temps de processus (CPU) et/ou de ressources mémoire et/ou de place d'archivage. Ceci implique qu'il est souvent difficile ou impossible de le mettre en oeuvre sur un grand nombre de cas, par exemple quelques milliers, ou dizaines de milliers, ou centaines de milliers, ou millions de cas, chaque cas étant défini par la valeur des différents paramètres. Or il est souvent nécessaire à l'homme de l'art de pouvoir estimer ce que seraient les réponses pour un nombre élevé de cas. Il lui faut donc un moyen d'avoir ces estimations. Pour arriver à son objectif (qui est souvent de comprendre ou visualiser un phénomène, ou d'optimiser un dispositif ou procédé, ou de faire des études de sensibilité, de variabilité ou de la modélisation inverse, ou de faire des études statistiques ou de rendement, ou de créer des programmes qui simulent certains résultats dans une gamme fixée de paramètres), il est connu de construire des "métamodèles". Un "métamodèle" est ici défini comme une fonce ion analytique, que l'on souhaite rapide à calculer par ordinateur, qui représente une réponse et qui dépend de plusieurs paramètres numériques. Ce "métamodèle", parce qu'il est très rapide à calculer, remplace le logiciel complet pour de nombreuses applications, par exemple la visualisation des phénomènes, l'optimisation, les études de sensibilité, de variabilité ou de modélisation inverse, les études statistiques ou de rendement, etc ... Un tel métamodèle est construit à partir d'un certain nombre de simulations et pour lesquelles les paramètres varient dans un certain domaine. Il est conseillé de choisir ces simulations suivant des techniques de plans de simulations numériques, par exemple décrites par Thomas J. Santner, Brian J. Williams, William 1. Notz, "Design and Analysis of Computer Experiments", Springer Series in Statistics, 2003 On appelle "base d'apprentissage" un tel ensemble de simulations faites avec le procédé et qui servent à bâtir le métamodèle. Typiquement, cette base d'apprentissage contient de quelques dizaines à plusieurs centaines ou. plusieurs milliers de simulations représentés par des vecteurs d'apprentissage, mais il n'y a aucun inconvénient à ce qu'elle en comporte moins ou plus. Vu que le procédé que nous considérons est déterministe, il y a une parfaite reproductibilité des réponses, et il est très fortement préférable que le métamodèle, appliqué à l'un des points de la base d'apprentissage, redonne soit exactement la valeur de réponse obtenue par le logiciel en ce point, soit une valeur 4 R = ail i2 ... ik ... id ' xll • x2 • ... • xk • .... xd avec, pour tous les ik, 0 ik S q et le degré total est q , c'est-à-dire que les 5 coefficients a;1 '2 ...'k ...'d sont nuls si la somme il + i:2 + ... + ik + ... + id est supérieure à q. Les coefficients al1 12 ... ik ... définissent le polynôme. Pour les déterminer, la méthode la plus classique est de faire une régression aux moindres carrés. Cela revient à minimiser la somme S des 15 carrés (éventuellement pondérés) des écarts entre réponse réellement simulée (issue du logiciel scientifique ou technique), c'est-à-dire la variance résiduelle : n S = j ' (Rsimulation, j ù Rméta mod èle (Xi ) )2 j=1 Les poids sont tous positifs ou nuls. XI est le vecteur 20 représentant l'ensemble des paramètres pour la jeme simulation, Rs mutation, j est la valeur de la réponse pour la je' simulation, Rmetamodel, est le métamodèle associé à cette réponse. Les détails sur les régressions aux moindres carrés et sur les méthodes utilisées pour déterminer les valeurs numériques des coefficients 25 ai, ;2 'k id sont connus de l'homme de l'art et sont par exemple décrits dans G.E.P. Box & N.R. Draper, "Empirical Model-Building and Response Surfaces", Wiley series in Probability and Mathematical Statistics, 1987 ou bien encore W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, "Numerical Recipes, The Art of Scientific Computing, third edition", Cambridge University Press, 2007. 30 Le principal inconvénient de la construction d'un métamodèle par des polynômes est le dilemme suivant : 10 • soit on impose que le métamodèle soit un interpolateur et alors on est obligé d'avoir un polynôme de degré total relativement élevé avec souvent un nombre total de coefficients a;, '2 ... ;k ... ;d non nuls égal (ou à peine inférieur) au nombre de simulations dans la base d'apprentissage. Dans ce cas il est connu de 5 l'homme de l'art que le comportement du polynôme est extrêmement chaotique en interpolation et encore plus en extrapolation, rendant ces interpolations et extrapolations peu fiables. • soit on impose un nombre plus limité de coefficients non nuls dans le polynôme pour s'assurer d'un comportement peu chaotique en interpolation et extrapolation. L'homme de l'art connaît plusieurs manière de choisir ces coefficients et de déterminer leurs valeurs: ces manières sont souvent basées sur l'analyse de la variance et sont décrites par exemple dans l'ouvrage précité de G.E.P BOX et N. DRAPER. Mais cette manière de faire présente plusieurs gros inconvénients : les fondements théoriques statistiques sont basés sur un modèle où la réponse obtenue peut être considérée comme le résultat d'un tirage d'une fonction aléatoire, ce qui n'est pas du tout le cas ici. De plus le métamodèle ainsi obtenu n'est pas un interpolateur ni un quasi-interpolateur, alors que nous avons vu qu'il était nettement préférable qu'il le soit. Pour pallier ces inconvénients, l'homme de l'art utilise de plus en plus fréquemment les méthodes diverses de krigeage ("kriging" en anglais) pour construire son métamodèle. Une bonne introduction à ces méthodes de krigeage (issues de la géostatistique) peut être trouvée par exemple dans l'ouvrage précité de Thomas J. SANTNER et collaborateurs. Le krigeage est une méthode d'in:erpolation spatiale, parfois considérée comme la plus juste d'un point de vue statistique, qui permet une estimation linéaire basée sur l'espérance mathématique et aussi sur la variance de la donnée spatialisée. Il fait intervenir la fonction de corrélation entre la valeur de la fonction aux points X; et X; et Xj représentant chacun un jeu de valeurs numériques des d paramètres, donc par exemple chacun une simulation par le logiciel. Mathématiquement, ce sont des éléments de Vd, c'est-à-dire des vecteurs constitués de d termes réels, et l'hypothèse de base est que cette fonction de corrélation ne dépend que d'une "distance" au sens mathématique du terme entre les points X; et Xi , notée X; - Xj Il et qui peut être définie par exemple 'far : d x(k) _x(k) '2 J k=1 Ok Xi X~ où x;(k) est la kème composante de Xi ( le kème paramètre). Les grandeurs Ok sont souvent nommées les "portées" du paramètre k. Elles doivent être non nulles et sont 5 usuellement choisies positives. Un exemple classique (mais non unique, il en existe bien d'autres) de fonction de corrélation r(X Xi) entre X et Xi est: r(X, Xi ) = exp X-X~ = exp i 2 Bk i Io Le métamodèle donné par le krigeage est de la forme : p n .Î(X)=Qi gi(X) +12, r(X,X1) i=l 15 où : • X représente un point courant quelconque, c'est-à-dire un jeu de valeurs numériques des d paramètres. • Les X;, i variant de 1 à n, sont des vecteurs représentant les n points (ou simulations, effectuées avec le logiciel) de la base d'apprentissage, chacun 20 de ces points étant un jeu de d valeurs numériques. • Les p fonctions gi(X) sont choisies a priori et sont nommées "termes de tendances" • r(X , X;) est la fonction de corrélation entre X et X; • Les p coefficients (3i et les n coefficients 2.; sont calculés 25 d'après les réponses et paramètres de la base d'apprentissage. Dans une variante du krigeage qui est le krigeage ordinaire, on a p=1 (1 seul terme de tendance) et gi(X) = 1 (terme constant et unitaire). 7 Dans le krigeage universel qui est une autre variante du krigeage, il y a plusieurs termes de tendance. Un choix classique de krigeage universel est de prendre p = d + 1 et de choisir les d effets linéaires : gi(X) = x0) et un terme constant gd+1 (x) = 1 Un autre choix classique de krigeage universel est de choisir les p termes les plus influents tels que déterminés par une analyse de la variance faite sur la base d'apprentissage avec comme fonctions candidates les différents monômes issu d'un polynôme de degré total 2 ou 3. La manière d'obtenir les p coefficients Pi et les n coefficients Xi est connue de l'homme de l'art et est décrite dans de nombreux ouvrages, par exemple dans l'ouvrage de Thomas J. SANTNER et Collaborateurs. Elle nécessite la résolution de systèmes linéaires de tailles n*n avec plusieurs seconds membres. Ce document décrit aussi des méthodes statistiques pour déterminer les portées Ok optimales, mais ce sont souvent des méthodes itératives qui nécessitent des quantités de calcul assez importantes. Le krigeage est une méthode classique de construire un interpolateur. On peut également utiliser les calculs intermédiaires du krigeage pour estimer (au prix de calculs un peu lourds) en chaque point un écart-type autour du méta-modèle et à partir des intervalles de confiances. Ces écart-types sont nuls pour les points d'apprentissage (c'est un interpolateur) et augmentent au fur et à mesure qu'on s'éloigne de ces points d'apprentissage. La figure 1 schématise ce comportement dans un cas monodimensionnel (d=1): La manière d'obtenir ces écart-types est connue de l'homme de l'art. 25 Elle est décrite dans de nombreux ouvrages, par exemple dans l'ouvrage précité de Thomas J. SANTNER et collaborateurs. Le krigeage a cependant un certain nombre d'inconvénients: • Le système linéaire n*n à résoudre est parfois peu stable, rendant sa résolution impossible ou très imprécise, en particulier lorsque n est 30 relativement élevé (par exemple de l'ordre du millier ou plus) et que les points de la base d'apprentissage ne sont pas régulièrement répartis. Pour pallier cet inconvénient, on peut utiliser ce que l'homme de l'art appelle des "pépites" (technique décrite par exemple dans le document précité de Thomas J. SANTNER et collaborateurs) mais de ce fait le krigeage perd sa qualité première d'interpolateur. 35 • L'expérience prouve que l'aptitude du krigeage à bien reproduire le comportement du logiciel sur des points tests choisis dans le même domaine que les points d'apprentissages mais distincts de ceux-ci est parfois décevante, c'est-à-dire que l'erreur quadratique moyenne entre valeurs de la réponse prédite par le métamodèle et valeur de la réponse obtenue par simulation avec le logiciel est dans de tels cas anormalement grande. L'invention concerne ainsi un procédé de construction d'un métamodèle pour interpoler des données à partir de simulations constituant une base d'apprentissage comportant des vecteurs X;, i variant de 1 à n, présentant chacun des composantes x;(k), k variant de 1 à d, caractérisé en ce que le métamodèle est une spline d'interpolation ou de quasi-interpolation multidimensionnelle pseudo-cubique de type plaque mince ayant une fonction analytique f(X) f(X) = n 2i X - Xi l 3 + d akx(k) + ao i=1 k=1 d ( (k) (k) 2 2 _ dilk -xi k-1 \ ak / V - X; p étant une constante par exemple supérieure ou égale à 103 ak désignant l'écart type des ketl1e composantes dilk étant d'autant plus grand que les effets non linéaires associés à ce paramètre sont grands et dilk une dilatation d'échelle éventuellement associée à ladite keme composante (dilk est d'autant plus grand que les effets non linéaires associés à ce paramètre sont grands) 2, et ah. étant les solutions d'un système linéaire symétrique de dimension (n + d + Xi -Xi • x(.k) = 0 J 0 n Â'jd 3 / ( k ) + -+ lak xi +ao =Yi pour i E {1, 2, 3, Poi k=1 pour kE {1, 2, 3, ... , d} 10 15 9 y; désignant la valeur de la réponse pour la ieme simulation pour la valeur X; des paramètres, w étant un poids associé à la ie"e simulation. Si les simulations ne sont pas pondérées, c'est-à-dire si elles ont toute la même importance, alors w = 1 quelque soit i. p est un paramètre de lissage, qui est infini pour une spline d'interpolation (vraie interpolation) Pour déterminer les dilatations d'échelle dilk, le procédé peut être caractérisé en ce qu'on effectue les opérations suivantes : a) partir d'une estimation initiale des d dilatations d'échel les dilk b) partitionner de manière aléatoire ou pseudo-aléatoire les n résultats de simulation Xi en np parties disjointes (np 1), c) itérer sur P parties avec P Plus particulièrement, on peut avoir p = 1, np = n et nq = 1 quelque soit q. Selon une variante préférée, ce procédé est caractérisé en ce que, t o pour calculer les dilatations d'échelle dilk, il comporte : a') Partir d'une estimation initiale de toutes les dilatations a d'échelles, qui respecte la condition fl dilk = 1, notamment dilk = 1 quelque soit k. k=1 b') Construire la fonction spline multidimensionnelle pseudocubique type plaque mince f(X) avec ces dilatations d'échelles. 15 c') Calculer les d dérivées secondes de la fonction spline multidimensionnelle pseudo-cubique type plaque mince par rapport au d paramètres, a2f(X; ) ax(k)2 pour k=l à d, et ce pour chacun des n points de la base 20 d'apprentissage (i variant de 1 à n). d') Calculer, pour chacun des d paramètres, la moyenne quadratique, sur l'ensemble des points de la base d'apprentissage, des dérivées secondes calculées à l'étape c, cette moyenne quadratique étant définie par : 25 SQD2F(k) n a2f(Xi ) 2 ax(k)2 i pour k variant de 1 à d n e') Calculer, pour ce jeu de dilatations d'échelles, une variance de validation croisée ou une variance de validation croisée généralisée. f) A partir de la 2eme itération du processus itératif et si cet écart quadratique de validation croisée est plus grand que celui obtenu à l'itération précédente, sortir du processus itératif avec les valeurs des dilatations d'échelles de l'itération précédente et conserver la fonction spline multidimensionnelle pseudo-5 cubique type plaque mince pour ces dilatations d'échelles de l'itération précédente. g') Calculer la moyenne géométrique C des quantités SQD2F(k) définies plus haut et pondérées par le carré de ak : (d C = fl 6k' SQD2F(k) k=1 h') Calculer de nouvelles dilatations d'échelles notamment par la newdilk = 6 _ k 10 formule : SQD2F(k) C i') Eventuellement, tester la convergence du processus itératif. Ceci peut être fait par exemple en testant le maximum des valeurs absolues des variations 15 relatives entre dilk et new dilk, c'est-à-dire en comparant la quantité maxk=làn dilk -new_dilk dil k ( avec par exemple, un seuil de convergence prédéfini. Si cette quantité est inférieure au seuil de convergence, la convergence est réputée atteinte et on sort du processus 20 itératif On conserve les nouvelles dilatations d'échelles obtenue en h'. il') (Optionnellement) Vérifier le bon déroulement du processus itératif, fait par exemple en vérifiant au moins une condition, à savoir si le nombre d'itérations est inférieur à un nombre maximum d'itération autorisé, et/ou si le rapport 25 max(new dilk , k=1 à d) / min(new dilk , k=1 à d) est inférieur à un autre seuil prédéfini, et si une dite condition n'est pas vérifiée, sortir du processus itératif, max et min désignant le maximum et le minimum de la valeur des d paramètres entre parenthèses. j') Actualiser les dilatations d'échelles, c'est-à-dire remplacer les 30 dilk , k=1 à d , par les new dilk , k=1 à d, et retourner à l'étape b' pour la prochaine itération, à partir de ces dilatations d'échelles actualisées. 12 Pour calculer les dilatations d'échelle, il est avantageux dans le cas où le métamodèle est un quasi-interpolateur d'utiliser une variance de validation croisé généralisée GCV. A cet effet le procédé est caractérisé en ce que ladite variance de validation croisée est une variance généralisée GCV, avec : n 1 co (f(xi ) - Yi )2 GCV=i=1 (n ù trace(Ap ))2 = n , et trace(Ap) est la trace de la matrice de lissage Ap qui fait passer du vecteur Y des n réponses y; au vecteur F des n valeurs quasi-interpolées f(X;), : soit F= Ap •Y avec F = (f(X1), f(X2), f(X), ... , f(Xä))` et Y = (y>, y2, ... , y;, yn)` l'indice t désignant la matrice transposée, la spline étant calculée avec un paramètre de lissage p non infini, mais tel que le métamodèle soit considéré comme un quasi-interpolateur. 15 Le paramètre trace (Ap) peut être calculé notamment par une méthode de Monte-Carlo. L'invention sera mieux comprise à la lecture de la description qui va suivre, la figure 1 illustrant les intervalles de confiance d'une méthode de krigeage selon l'art antérieur, dans le cas monodimensionnel, c'est-à-dire qu'il y a 20 seulement d = 1 paramètre, avec en abscisse ce paramètre explicatif et en ordonnée la valeur du métamodèle, et les figures 2 et 3 étant des organigrammes correspondant à deux procédés de calcul des dilatations d'échelle. L'idée de base de l'invention est donc d'utiliser comme métamodèles des splines d'interpolation ou de quasi-interpolation 25 multidimensionnelles pseudo-cubiques type plaque mince. Ces splines, connues mais relativement rarement utilisées dans la littérature, font partie de la famille des splines multidimensionnelles invariantes par rotation et utilisent des fonctions radiales de base et sont connues de l'homme de l'art par le document suivant : - J. Duchon, "Splines minimizing rotation-invariant semai-norms in 30 Sobolev space", Lecture Notes in Mathematics, vol 571., pp85-100, 1977, also in the 10 et 10, proceedings of a conference "Constructive Theory of Functions of Several Variables", held at Oberwolfach, april 25- May 1, 1976 and edited by W. Schemp and K. Zellers (Eds.), Berlin: Springer-Verlag pp.85-100. La manière de résoudre un système linéaire symétrique est classique, très bien connue de l'homme de l'art, décrite dans de nombreux ouvrages, programmée dans de nombreux logiciels ou bibliothèques mathématiques. Les familles de splines multidimensionnelles utilisées dans le cadre de la présente invention sont des fonctions de Vd dans V, où d est la dimension de l'espace, c'est-à-dire le nombre de paramètres. De manière générale, une spline de lissage est une fonction f qui minimise une énergie totale, somme d'une énergie de courbure et d'une variance résiduelle multipliée par le paramètre de lissage p : n Et =Ec(f)+P.1coi(f(Xi)-yi)2 i=1 où n est le nombre de points, y; la valeur de la réponse pour la ième simulation et X; le 15 vecteur des variables explicatives, c'est-à-dire des valeurs des d paramètres pour la i'n'e simulation, w; le poids affecté au ième point ou simulation. Chaque famille de spline se caractérise par l'expression de son énergie de courbure. L'énergie de courbure pour les splines multidimensionnelles pseudo cubiques type "plaque mince" est , d'après le document précité de J. 20 DUCHON: d d E~(f)= fll Rd k=1 P=1 82f atk atP , du, du, .... dud Four 2 d-1 u (u) d est la dimension de l'espace des variables explicatives et Four la transformée de Fourier dans Vd , le carré de la transformée de Fourier étant le carré de la norme. 25 La régression aux moindres carrés correspond au cas limite de p La spline d'interpolation correspond au cas limite de p ù> : la spline doit passer par tous les points, c'est-à-dire que la variance résiduelle est nulle, et on en minimise l'énergie de courbure. Une spline d'interpolation est un cas limite 30 d'une spline de lissage. La spline de quasi-interpolation correspond au cas où le paramètre de lissage p est grand, c'est-à-dire que la variance résiduelle est petite par rapport à la 5 1015 variance totale des y;, c'est-à-dire que la variance résiduelle est inférieure à par exemple 0,01 fois la variance totale, ou, mieux, que la variance résiduelle est inféreire à par exemple 0,001 fois la variance totale, ou mieux encore, que la variance résiduelle est inférieure à par exemple 0,0001 fois la variance totale. La variance résiduelle est ici définie (très classiquement) par : n Var résid = (y j ùf (X j ) 2 j=1 et la variance totale est ici définie (classiquement) par : n Var totale= .(yjù )2 j =1 ou est la valeur moyenne pondérée par les w; des n valeurs yi. L'utilisation de splines de lissage multidimensionnelles pseudocubiques type plaque mince pour analyser des résultats expérimentaux est déjà connue de l'homme de l'art, par exemple de la publication de F. DE CRECY "Pseudo-cubic thin-plate-spline method for analysing experimental data", paru dans Nuclear Engineering and Design, vol. 149, pp. 459-465 (1994). La forme analytique d'une spline multidimensionnelle pseudocubique type plaque mince, qu'elle soit d'interpolation ou de lissage, est : x-x n f(X)=1.1 i=1d l3 + ak x") + ao k=1 où x(k) est la kè°1e composante (ou paramètre) de X , Xi est le vecteur des d valeurs 20 des paramètres pour la ième simulation, et où la norme est définie par : 2 d (x(k) ùx (k) 2 Il X - X; _ dilk k=1 ak / 6k étant une grandeur caractéristique positive associé aux kemes composantes des X; , par exemple l'écart-type des valeurs du kèmes composantes des X1 (c'est-à-dire des 25 valeurs du keme paramètre) et ou dilk est la dilatation d'échelle associée au keme paramètre. Ces dilatations d'échelle traduisent, le cas échéant, l'importance relative des effets non linéaires des divers paramètres. Elles sont généralement choisies positives. Plus elles sont grandes, plus les effets non linéaires associés à ce paramètre sont importants. Il est fréquent et conseillé, mais non obligatoire, d'imposer une condition supplémentaire sur ces dilatations d'échelles, par exemple: d Ud11 Une telle spline est définie sur tout V7d. Quand on s'éloigne de la 5 zone où il y a des points X1 , elle tend localement vers un hyperplan. Pour chaque valeur du paramètre de lissage p et des dilatations d'échelle dilk, les coefficients 2 et îk sont solutions d'un système linéaire symétrique de dimension (n+d+l) x (n+d+l). Ce système est : 3 z. d + 1 + l ak • xi k) + ao =Yi J ~1 k=1 Xj-Xi 10 pour iE {1,2,3,...,n} n l 'l1 • x(k) = 0 1=1 pour kE {1, 2, 3, ... , d} n 12i = 0 i=1 15 L'homme de l'art sait résoudre (avec un ordinateur) ce système linéaire. De nombreuses méthodes existent, dont les plus classiques sont décrites dans le document précité de W.H. PRESS et collaborateurs. Les notations sont les mêmes que plus haut au sein de ce chapitre: en particulier, la norme fait intervenir les dilatations d'échelle. Les y; sont les 20 grandeurs à lisser ou à interpoler (résultats des mesures ou réponse de simulations: y; est la valeur de la réponse pour la ieme simulation) pour la valeur X1 des paramètres. Une spline de lissage de quasi-interpolation est obtenue lorsque le paramètre de lissage p est suffisamment grand pour que le métamodèle soit une quasi-interpolation. 25 Dans le cas d'une spline d'interpolation, on a p = oD et ce système se simplifie directement en : Xj ù 16 3 d k + ak • XI + ao =y; pour i e {1, 2, 3, ... , n} j=1 k=l n 12 •x~k) =0 j=1 pour kE 11, 2, 3, ... , d} n 12j=0 j=1 La diagonale principale de la matrice associée à ce système linéaire symétrique ne comporte donc que des zéros dans le cas de l'interpolation. Il est conseillé de choisir le cas échéant des dilatations d'échelle dilk associées aux divers paramètres, si on veut tenir compte des non-linéarités. Une première variante est d'utiliser une méthode de type "Validation croisée"simple VVC. Les techniques de validation croisée sont classiques et sont décrites notamment dans l'article de Ron KOHAVI "A Study of Cross-validation and Bootstrap for Accuracy Estimation and Model Sélection" International Jont Conference on Artificial Intelligence ISCA, 1995. Elle comporte schématiquement dans la procédure suivante (voir figure 2) : 1. Partir d'une estimation initiale des d dilatations d'échelles, dilk ,kE {1, 2, 3, ...,d}. 2. Partitionner de manière aléatoire ou pseudo-aléatoire l'ensemble des n résultats de simulation de la base d'apprentissage en np parties disjointes. Un choix classique est de prendre np entre 5 et 50. Soit nq le nombre de résultats dans la genre partie. Un choix classique est d'avoir le même ordre de grandeur de nq quelque soit la partie. Un choix conseillé est d'avoir 1 nq égal à n ou proche de n. Si on est très large en temps de calcul, un choix ccîriseillé, sûr et robuste (mais long) est de prendre np=n et nq=1, quelque soit q. 3. Choisir P parties sur lesquelles on va itérer. Il est préférable d'avoir P = np , mais cela entraîne souvent des calculs très longs. Aussi est-il parfois acceptable de prendre P Une difficulté peut parfois survenir, à savoir que la matrice associée au système linéaire (n+d+1)*(n+d+l) à résoudre est mal conditionnée, même avec des paramètres centrés réduits. Cela peut arriver par exemple quand la répartition des simulations dans Od est très irrégulière, avec quelques simulations beaucoup plus proches les unes des autres que la moyenne. La man:ère de résoudre t 0 cette difficulté est d'approximer la spline d'interpolation par une spline de lissage avec un paramètre de lissage p très grand, par exemple dans la gamme [ 103 à 1010] si tous les paramètres et réponses ont été centrés et réduits. Avec des grandes valeurs du paramètre de lissage p , la spline de lissage est très proche de la spline d'interpolation: le métamodèle obtenu est un "quasi-interpolateur". 15 Une deuxième variante pour choisir les dilatations d'échelle (uniquement dans le cas du quasi-interpolateur, pas dans le cas du vrai interpolateur) est la validation croisée généralisée GCV. Cette autre manière est basée sur le fait qu'une bonne approximation de la variance de validation croisée simple VVC, définie plus haut, est GCV, nommée variance de validation croisée généralisée. 20 définie par: n 1 coi (f la sommation se faisant sur les n points de la base d'apprentissage exprimés en U. On peut calculer une quantité similaire à partir de la spline f(X) exprimée en X : ( 2 n a2f(Xi) ax(k)2 n La relation liant u(k) et x(k) implique que: a'h(U.) 82f(X;) *' 6k au(k)2 ax(k)2 dllk et que donc : SQD2U(k) = SQD2F(k)*(6k/dilk)2 La minimisation de VVC ou de GCV se fait sur le chemin qui vise à avoir le même SQD2U(k) pour tous les paramètres k, tout en respectant la contrainte que le produit de toutes les dilatations d'échelles est égal à 1. Pour cela, une méthode proposée est le processus itératif suivant : SQD2U(k) = SQD2F(k) _ 1. On part d'une estimation initiale de toutes les dilatations d'échelles, qui respecte la condition : d n dilk =1 k=1 Si on n'a pas de connaissance préalable de ces dilatations d'échelles initiales, les choisir toutes égales à 1 est généralement un bon choix. 2. Construire la fonction spline multidimensionnelle pseudo- cubique type plaque mince avec ces dilatations d'échelles, comme expliqué plus haut. 3. Calculer les d dérivées secondes de la fonction spline multidimensionnelle pseudo-cubique type plaque mince par rapport au d paramètres, a2f(X ax(k)2 pour k=1 à d, et ce pour chacun des n points de la base d'apprentissage (i = 1 à n 4. Calculer, pour chacun des d paramètres, la moyenne quadratique (sur l'ensemble des points de la base d'apprentissage) des dérivées secondes calculées à l'étape précédente, cette moyenne quadratique étant définie par : SQD2F(k) = pour k = 1 à d 5. Calculer, pour ce jeu de dilatations d'échelles, soit la variance de validation croisée VVC, soit la variance de validation croisée généralisée GCV, ces deux grandeurs ayant été définies plus haut et la manière de les calculer ayant été indiquée. 6. Si on n'est pas à la lere itération du processus itératif et si cette 25 quantité VVC ou GCV est plus grande que celle obtenue à l'itération précédente, sortir du processus itératif avec les valeurs des dilatations d'échelles de l'itération précédente et en recalculant ou en restaurant la fonction spline multidimensionnelle pseudo-cubique type plaque mince pour ces dilatations d'échelles de l'itération précédente. 20 7. Calculer la moyenne géométrique C des quantités SQD2F(k) définies plus haut et pondérées par le carré de 6k : (d \id C = flo SQD2F(k) \ k=1 8. Conceptuellement, cette quantité représente la valeur moyenne 5 géométrique des SQD2U(k) . 9. Calculer des nouvelles dilatations d'échelles par une formule new_dilk = 6k telle que: SQD2F(k) C 10. Tester la convergence du processus itératif. Ceci peut être fait 10 par exemple en testant le maximum des valeurs absolues des variations relatives entre dilk et new_dilk, c'est-à-dire en comparant la quantité maxk=1àn t dilk -new_dilk à, par exemple, un seuil de convergence prédéfini. Si cette quantité est inférieure au seuil de convergence, la convergence est réputée atteinte et on doit sortir du 15 processus itératif. 11. Vérifier le bon déroulement du processus itératif. Ceci peut être fait par exemple en vérifiant si le nombre d'itérations est inférieur au nombre maximum d'itération autorisé, et par exemple aussi en vérifiant si le ratio max(new_dilk , k=1 à d) / min(new dilk , k=1 à d) est inférieur à un 20 autre seuil prédéfini. Si par exemple une de ces deux conditions n'est pas vérifiée, sortir du processus itératif 12. Actualiser les dilatations d'échelles, c'est-à-dire remplacer les dilk , k=1 à d , par les newdilk , k=1 à d, et retourner à l'étape 2 pour la prochaine itération. 25 L'utilisation, à l'étape 5 du processus ci-dessus, de la variance de validation croisée VVC ou de la variance de validation croisée généralisée GCV correspondent à deux variantes possibles de ce mode de réalisation de l'invention. L'utilisation de la variance de validation croisée généralisée GCV est généralement plus rapide en temps de calcul, niais est parfois moins robuste et ne s'applique pas 30 aux vrais interpolateurs mais aux quasi-interpolateurs. L'utilisation de la variance de 23 validation croisée VVC est plus robuste et est conseillée chaque fois qu'on peut se permettre le surcoût en temps de calcul tout en conservant un nombre de parties dans la partition suffisamment élevé, par exemple supérieur à 10 ou 15, et en utilisant toutes ces parties dans le calcul de VVC, c'est à dire que P=np dans l'étape 5 de la description du mode de calcul de la méthode "validation croisée". L'utilisation de VVC et non pas de GCV est aussi fortement conseillée quand on utilise un vrai interpolateur et non pas un quasi-interpolateur. Lorsque les simulations sont pondérées, c'est-à-dire affectées de poids w; non tous égaux entre eux, une variante conseillée du mode de réalisation de l'invention est de prendre SQD2F(k), aux étapes 4 et suivantes du processus itératif ci-dessus, comme moyenne quadratique pondérée (sur l'ensemble des points de la base d'apprentissage) des dérivées secondes calculées à l'étape précédente, cette moyenne quadratique pondérée étant définie par : SQD2F(k) = , n coi(a2f(XI) \ z =1 n pour k = 1 à d ax,k)2 1 ;=1 Les principaux avantages apportés par cette invention sont les suivants: • Le métamodèle construit avec la spline d'interpolation multidimensionnelle pseudo-cubique type plaque mince est bien un interpolateur ou 20 un quasi-interpolateur, ce qui est un avantage considérable par rapport aux polynômes. • L'expérience montre que le système linéaire à résoudre pour la spline d'interpolation multidimensionnelle pseudo-cubique type plaque mince est plus stable numériquement que le système linéaire à résoudre pour le krigeage. • II n'est pas nécessaire de faire le choix, toujours un peu arbitraire, des termes de tendance du krigeage universel. • La spline d'interpolation multidimensionnelle pseudo-cubique type plaque mince est le seul interpolateur qui minimise l'énergie de courbure susmentionnée. • La spline d'interpolation multidimensionnelle pseudo-cubique type plaque mince est une fonction qui devient rapidement très lisse quand on s'éloigne des points expérimentaux. Au fur et à mesure qu'on s'éloigne des zones où il y a les points expérimentaux, la spline tend localement vers un hyperplan. 25 30 24 • L'expérience montre que la spline d'interpolation multidimensionnelle pseudo-cubique type plaque mince arrive mieux à retrouver (par rapport au krigeage) des points tests qui n'ont pas servi à bâtir la spline et qui sont dans la même gamme de paramètres. • La spline d'interpolation multidimensionnelle pseudo-cubique type plaque mince est un peu plus rapide à utiliser que le krigeage. L'ordinateur a à calculer des racines carrées au lieu d'exponentielles. • Il n'y a aucune contrainte sur l'emplacement des simulations dans l'espace des paramètres \7d. Si plusieurs simulations sont au même endroit dans l'espace des paramètres Vd , il suffit de n'en retenir qu'une seule de ces simulations (qui donnent le même résultat, puisque le logiciel est supposé déterministe). • Si le logiciel n'est pas parfaitement déterministe et qu'il fournit des réponses avec une certaine variabilité, on peut passer continûment de la spline d'interpolation à une spline de quasi-interpolation en jouant sur le paramètre de lissage p . • Il est possible de dériver indéfiniment la spline d'interpolation multidimensionnelle pseudo-cubique type plaque mince presque partout, sauf sur les points expérimentaux où elle est de classe C2, c'est-à-dire qu'on peut la dériver 2 fois sur les points expérimentaux et que toutes les dérivées partielles d'ordre total inférieur ou égal à 2 sont définies et continues. • La spline d'interpolation multidimensionnelle pseudo-cubique type plaque mince permet d'utiliser un nombre d de paramètres assez élevé. L'expérience montre qu'il n'y a pas d'inconvénient à ce que ce nombre d de paramètres soit de l'ordre de 30 ou même supérieur. Il faut par contre qu'il reste nettement inférieur au nombre n de simulations faites. L'invention s'applique à des procédés de simulations déterministes dans tous les secteurs techniques et industriels. Elle peut donc concerner une très large gamme d'applications, parmi lesquels on peut citer à titre d'exemple : la microélectronique et les micro- et nano-technologies (avec par exemple les logiciels de "TCAD" (Technology Computer Aided Design), les logiciels de "SPICE" (Simulation Program with Integrated Circuit Emphasis)), la construction mécanique et automobile, les télécommunications, les techniques aéronautique et/ou spatiale, le bâtiment et les travaux publics, la chimie, l'extraction minière, pétrolière et gazière, l'énergie, l'agro-alimentaire, le transport, etc...35 REVENDICATIONS1) Procédé d'interpolation de données techniques en construisant un métamodèle à partir de simulations constituant une base d'apprentissage comportant des vecteurs de données techniques X;, i variant de 1 à n, présentant chacun des composantes x;(k), k variant de 1 à d, caractérisé en ce que le métamodèle est une spline d'interpolation ou de quasi-interpolation multidimensionnelle pseudo-cubique de type plaque mince ayant une fonction analytique f(X) n f(X)=12i i=1 d 3 + Eakx(k) +ao k=1 d (x 2 (k) ù. _ dilk k=1 6k crk désignant l'écart type des keme composantes X;(k) des vecteurs X, de données techniques et dilk la dilatation d'échelle associée à ladite keme composante 2, et ak étant les solutions d'un système linéaire symétrique de dimension (n + d + 1)2 : X ù xi d 3 /li ( k ) + + E ak ' xi + ao =Y i Pua k=1 j=1 X j ùXi pour i E {1, 2,3,...,n} n pour k E {1, 2, 3, ... , d} 26 y; désignant la valeur de la réponse pour la ieme simulation pour la valeur X; des paramètres de données techniques, p étant un paramètre de lissage dont la valeur est infinie pour une vraie interpolation, 6); étant un poids associé à la ième simultion. 2) Procédé selon la revendication 1, caractérisé en ce que pour déterminer des dilatations d'échelles dilk, on effectue les opérations suivantes : a) partir d'une estimation initiale des d dilatations d'échelles dilk b) partitionner de manière aléatoire ou pseudo-aléatoire les n 1 o résultats de simulation X; en np parties disjointes (p > 1), c) itérer sur P parties avec P np d) initialiser à zéro une variance de validation croisée e) calculer une variance de validation croisée à partir des valeurs de dilatations d'échelle dilk 15 f) choisir de nouvelles valeurs de dilatations d'échelle minimisant ladite variance de validation croisée obtenue lors de e, et g) itérer à partir de e avec lesdites nouvelles valeurs de dilatation d'échelle g') éventuellement entre e et f ou entre f et g tester la convergence 20 selon un critère de convergence, et si ce critère est satisfait, sortir du processus itératif avec les valeurs de dilatations d'échelle en cours. 3) Procédé selon la revendication 2, caractérisé en ce que P 4) Procédé selon une des revendications 2 ou 3, caractérisé en ce que 25 e comporte, pour p variant entre 1 et P el) supprimer provisoirement tous les résultats de simulation de cette partie e2) déterminer la fonction f(X) sans ces points. e3) utiliser la fonction f(X) uniquement pour recalculer les valeurs 30 prédites sur les nq points provisoirement supprimés. e4) incrémenter la variance de validation croisée par la somme des carrés des écarts entre la valeur prédite sans les points provisoirement supprimés et la valeur de la simulation sur les points provisoirement supprimées, cette somme se faisant sur tous les points provisoirement supprimés. 35 e5) réintégrer les points provisoirement supprimés. 5) Procédé selon une des revendications 2 à 4, caractérisé en ce que la valeur de np est comprise entre 5 et 50. 6) Procédé selon une des revendications 2 à 5, caractérisé en ce que nq étant le nombre de résultats dans la qme partie, nq est sensiblement le même pour chaque partie. 7) Procédé selon la revendication 6, caractérisé en ce que np 1nq = n q=1 8) Procédé selon la revendication 7, caractérisé en ce que p = 1, et donc np = n, et en ce que nq = 1 quelque soit q. 9) Procédé selon la revendication 1., caractérisé en ce que pour calculer les dilatations d'échelle dilk, il comporte : a') Partir d'une estimation initiale de toutes les dilatations d d'échelles, qui respecte la condition 11 dilk = 1, notamment dilk = 1 quelque soit k. k=I b') Construire la fonction spline multidimensionnelle pseudocubique type plaque mince f(X) avec ces dilatations d'échelles. c') Calculer les d dérivées secondes de la fonction spline 15 multidimensionnelle pseudo-cubique type plaque mince par rapport aux d a2f(X;) paramètres, pour k=1 à d, et ce pour chacun des n points de la base ax(k)- d'apprentissage (i variant de 1 à n). d') Calculer, pour chacun des d paramètres, la moyenne quadratique sur l'ensemble des points de la base d'apprentissage, des dérivées secondes calculées 20 à l'étape c', cette moyenne quadratique étant définie par 2 1 a2f(Xt) i=1, ax(k)2 , n e') Calculer, pour ce jeu de dilatations d'échelles, une variance de validation croisée ou une variance de validation croisée généralisée. 25 f) A partir de la 2e1e itération du processus itératif et si ladite variance est plus grande que celle obtenue à l'itération précédente, sortir du processus itératif avec les valeurs des dilatations d'échelles de l'itération précédente et conserver la fonction spline multidimensionnelle pseudo-cubique type plaque mince pour ces dilatations d'échelles de l'itération précédente. SQD2F(k) = pour k compris entre 1 à dg') Calculer la moyenne géométrique C des quantités SQD2F(k) définies plus haut et pondérées par le carré de 6k : C = flakz SQD2F(k) k=l cette quantité représentant la valeur moyenne géométrique des SQD2U(k) . h') Calculer des nouvelles dilatations d'échelles notamment par la formule suivante: new_dilk = 6k SQD2F(k) C i') éventuellement tester la convergence du processus intéractif et sortir du processus itératif avec lesdites nouvelles dilatations d'échelles obtenues en 10 h'. j') Actualiser les dilatations d'échelles, c'est-à-dire remplacer les dilk , k=l à d , par les new dilk , k=1 à d, et retourner à l'étape b' pour la prochaine itération. 10) Procédé selon la revendication 9, caractérisé en ce que le test 15 de convergence à l'étape i') prévoit de tester le maximum des valeurs absolues des variations relatives entre dilk et new dilk, c'est-à-dire en comparant la quantité maxk=làn dilk -new dilk dilk / avec un seuil de convergence prédéfini, et si cette quantité est inférieure au seuil de 20 convergence, sortir du processus itératif. 11) Procédé selon une des revendications 9 ou 10, caractérisé en ce qu'il comporte entre i' et j', une étape i'1) pour vérifier au moins une condition, à savoir que le nombre d'itérations est inférieur au nombre maximum d'itération autorisé, et/ou que le rapport 25 max(newdilk , k=l à d) / min(new_dilk , k=1 à d) est inférieur à une valeur seuil prédéfinie, et si une dite condition n'est pas vérifiée, sortir du processus itératif. 12) Procédé selon une des revendications 2 à 11, caractérisé en ce que le métamodèle est un métamodèle de quasi-interpolation et en ce que ladite 30 variance de validation croisée est une variance généralisé GCV, avec :n oi(.f(xi-Yi)2 GCV = i =1 (n ù trace(Ap ))2 il et =n, et trace(Ap) est la trace de la matrice de lissage Ap qui fait passer du vecteur Y des n réponses y, au vecteur F des n valeurs quasi-interpolées f(X,), : soit F= Ap •Y avec F = (f(X 1), f(X2), f(X,), ... , f(Xä ))t et Y=(Y1,Y2, ... ,y;, Yä) r l'indice t désignant la matrice transposée, la spline étant calculée avec un paramètre de lissage p compris entre 104 et 101°. 13) Procédé selon la revendication 12, caractérisé en ce que trace (Ap) est déterminée par une méthode de Monte-Carlo.