Durant Ma Thèse, En Particulier

Amine Ammar

Elías Cueto

Chady Ghnatios

Xavier Kestelyn

Fodil Meraghni

J'aimerais tout d'abord remercier Paco Chinesta pour toutes les belles opportunités qu'il m'a apportées, pour la conance qu'il m'a accordée, pour la quantité incroyable de choses que j'ai apprises à son contact et pour tous les bons moments que nous avons passés ensemble. J'espère que nous pourrons continuer à collaborer encore longtemps et accomplir de belles choses. Je remercie également Charbel Farhat, bien que je n'aie pas pu aller en Californie, j'ai beaucoup apprécié et appris de nos échanges.

Merci à ESI Group de m'avoir fait conance pour travailler au sein de la chaire et collaborer avec les équipes. Je remercie Anoop, Morgan, Satya, Krishna, Bhushan et tous les autres avec qui j'ai travaillé pendant ces trois années. Mention spéciale à Fatima et Mustapha, pour leur gentillesse innie et tous les agréables moments et superbes projets que nous avons partagés.

Je remercie aussi l'équipe de CNRS@CREATE et du projet DesCartes de m'avoir accueilli à Singapour, que ce soit pour l'accueil chaleureux, l'aide administrative ou pour les travaux que nous avons

Réduction de modèle avancée et approches basées sur les données pour la construction de jumeaux numériques augmentés par la physique Résumé Au 20e siècle, l'ingénierie a fait des progrès remarquables dans divers domaines, tandis que d'autres disciplines se sont tournées vers les données à des ns de diagnostic et de pronostic. Reconnaissant le potentiel des données et de l'IA, les sciences de l'ingénieur ont adopté ces technologies pour faire de meilleures prédictions, améliorer les performances et mieux comprendre les systèmes complexes. C'est ainsi qu'est né le paradigme du jumeau numérique.

Le dé consiste à développer des modèles précis capables de prédire des résultats en fonction de données d'entrée. Dans le domaine de l'ingénierie, il peut être dicile de choisir entre les approches basées sur la physique et les approches basées sur les données. Les approches basées sur la physique orent des avantages, mais elles sont Réduction de modèle avancée et approches basées sur les données pour la construction de jumeaux numériques augmentés par la physique ix Résumé étendu Contexte et motivations Au 20e siècle, l'ingénierie s'est fortement appuyée sur des modèles pour concevoir des composants et des systèmes de manière ecace. Ces modèles utilisaient des équations algébriques ou des équations aux dérivées partielles pour décrire les relations entrée/sortie dans l'ingénierie. Alors que les modèles algébriques étaient faciles à manipuler et à résoudre, les modèles basés sur des équations étaient de plus grands dés, nécessitant des techniques de discrétisation telles que celle des éléments nis pour donner au problème une forme résoluble.

La discrétisation consiste à résoudre le problème en un nombre limité de points, aussi appelés n÷uds, et d'instants, au lieu d'essayer de le résoudre pour un nombre inni de points et d'instants. Le nombre d'inconnues du problème augmente avec la dimension de l'espace, ce qui entraîne une croissance exponentielle de la complexité. L'utilisation des ordinateurs et calculateurs a rendu la discrétisation possible, mais les problèmes à grande échelle exigent toujours un temps de calcul important, nécessitant des jours, des semaines, voire des mois pour être résolus.

Ces méthodes ont permis l'avènement du paradigme de l'ingénierie basée sur la simulation [18,19]. Ces techniques [START_REF] Peiró | Finite dierence, nite element and nite volume methods for partial dierential equations[END_REF][START_REF] Mattiussi | An analysis of nite volume, nite element, and nite dierence methods using some concepts from algebraic topology[END_REF] regroupent notamment la méthode des éléments nis [START_REF] Zienkiewicz | The nite element method in engineering science[END_REF][START_REF] Bathe | Finite element method[END_REF], la méthode des diérences nies [START_REF] Jordán | Calculus of nite dierences[END_REF][START_REF] Eh | The calculus of nite dierences[END_REF] et la méthode des volumesnis [START_REF] Versteeg | An introduction to computational uid dynamics: the nite volume method[END_REF][START_REF] Eymard | Finite volume methods, Handbook of numerical analysis[END_REF].

L'ingénierie basée sur la simulation requiert de multiples calculs pour explorer l'espace paramétrique et identier les paramètres de conception optimaux. Toutefois, le coût des simulations et la dimension élevée de l'espace paramétrique restent des dés permanents. En outre, la qualité des modèles paramétriques en ingénierie est aectée par des facteurs tels que la taille du système analysé, le grand nombre et la variabilité des paramètres impliqués et l'incertitude inhérente à l'écart entre le modèle et la réalité.

Pour relever ces dés, les ingénieurs utilisent diverses techniques et stratégies pour réduire les sources d'incertitude et garantir des prédictions précises pour des conceptions ables.

De l'ingénierie aux jumeaux numériques

Le 20e siècle a été le témoin de réalisations remarquables [START_REF] Constable | A century of innovation: Twenty engineering achievements that transformed our lives[END_REF] dans le domaine de l'ingénierie, notamment de progrès dans l'aéronautique, l'exploration spatiale, les transports, l'énergie et les infrastructures civiles. Pendant ce temps, des domaines tels que le marketing, l'économie, les sciences sociales et la médecine se sont orientés vers l'utilisation des données à des ns de diagnostic et de pronostic, motivant de fortes améliorations technologiques liées à la science des données. C'est ainsi que l'ingénierie a connu trois évolutions majeures : L'extension de l'ingénierie à l'utilisation : Avec l'avènement de technologies telles que l'Internet des objets (Internet of Things -IoT) et l'abondance des données, les fastidieuses méthodes traditionnelles d'ingénierie en amont ont été en partie remplacées par des opérations en temps réel. L'intérêt est passé d'une ingénierie basée sur le produit à une ingénierie basée sur le service, accentuant les prévisions en continu et le maintien de la performance tout au long de la durée de vie d'un produit.

x Réduction de modèle avancée et approches basées sur les données pour la construction de jumeaux numériques augmentés par la physique L'ingénierie à l'ère des systèmes connectés : L'ingénierie a élargi son champ d'action des composants individuels aux systèmes interconnectés, ce qui a conduit à des concepts tels que les industries, les villes et les nations intelligentes. Cette évolution a marqué la quatrième révolution industrielle, dans laquelle les données se sont intégrées aux révolutions précédentes, qui étaient alimentées par la vapeur, l'électricité, l'électronique et l'automatisation.

La prochaine révolution de l'ingénierie immersive centrée sur l'homme : La cinquième révolution, qui est déjà en marche, place l'homme au c÷ur de l'ingénierie, envisageant un métavers sensible à la physique. Les progrès réalisés dans le domaine de la vision par ordinateur, des dispositifs multimodaux et du traitement du langage naturel pourraient permettre de dépasser le test de Turing. Cependant, l'intégration de l'homme dans le système pose des dés (notamment d'interfaces homme-machine et de calcul temps réel) aux approches de modélisation traditionnelles.

En réponse à ces développements, les sciences de l'ingénieur se sont intéressées aux possibilités oertes par les données et l'IA, dans le but d'obtenir de meilleures prédictions, d'améliorer les performances et de mieux comprendre les systèmes complexes [START_REF] Madni | Leveraging digital twin technology in model-based systems engineering[END_REF].

L'apprentissage automatique, en particulier, ore un processus d'apprentissage général dans lequel la relation entre les entrées et les sorties peut être établie à partir des données collectées, ce qui permet l'application en temps réel des modèles appris. Ces nouvelles pratiques se placent dans le paradigme du jumeau numérique.

Le principal dé consiste à obtenir des modèles de régression précis qui prédisent la sortie en fonction de l'entrée. Il existe diérentes techniques, qui dépendent de facteurs tels que la quantité de données, le temps d'entraînement et la qualité des données.

Le choix entre ces techniques permet de répondre à des problématiques telles que l'obtention de données [4,5], le coût des données entraînant un besoin de travailler avec peu de données [6,7], l'interprétabilité du modèle [10,11] et l'extrapolation du modèle [8,9]. Bien que la modélisation basée sur les données présente des avantages, de nombreux dés persistent, notamment la diculté et le coût de la collecte de données, l'interprétabilité des modèles et l'extrapolation en dehors des conditions observées.

Pour relever ces dés, il faut examiner attentivement les stratégies de collecte de données, l'interprétabilité et l'utilisation appropriée des modèles dans leurs domaines de validité.

Entre la physique et les données : les modèles hybrides

Dans le domaine de l'ingénierie, le choix entre les approches fondées sur la physique et les approches fondées sur les données est dicile. Les approches fondées sur la physique orent des avantages, mais elles sont gourmandes en calculs et ont du mal à gérer les systèmes à grande échelle et l'incertitude. Les approches fondées sur les données sont prometteuses en l'absence de modèles précis, mais elles se heurtent à des dicultés telles que le coût élevé des données, les risques d'extrapolation et le manque d'explications et de certications.

La combinaison des deux approches [START_REF] Karniadakis | Physics-informed machine learning[END_REF][START_REF] Liu | Physics-augmented learning: A new paradigm beyond physics-informed learning[END_REF] semble donc être le choix optimal, puisqu'il ore un bon compromis entre les avantages et inconvénients de chacune. En combinant les approches fondées sur la physique et sur les données, nous pouvons exploiter les points forts de chaque approche et obtenir de meilleurs résultats en matière d'ingénierie.

Réduction de modèle avancée et approches basées sur les données pour la construction de jumeaux numériques augmentés par la physique xi Cette alliance présente l'avantage de réduire considérablement les besoins en données pour la construction de modèles. Cette réduction est obtenue en exploitant les lois de la physique connues, ou en se concentrant sur l'écart entre le modèle et la réalité, qui est censé être plus simple à approximer que le comportement réel. En outre, l'utilisation de modèles basés sur la physique permet d'expliquer les aspects fondamentaux du modèle et les prédictions qui en résultent, ce qui permet de certier les résultats, et parfois même d'atténuer les problèmes d'extrapolation.

La réduction des besoins en données n'est pas seulement permise par l'utilisation des modèles physiques, mais aussi au fait que la physique peut informer le choix des emplacements et instants optimaux de collecte des données. Cela est particulièrement évident dans le cadre de l'apprentissage actif, que l'incorporation des connaissances existantes fondées sur la physique peut rendre plus ecace.

Pour illustrer cela, prenons l'exemple de la mesure de la température dans une ville.

Au lieu de placer de nombreux thermomètres dans chaque rue et d'enregistrer la température en permanence, la connaissance de la physique et le bon sens nous indiquent que quelques thermomètres placés stratégiquement, par exemple un dans chaque quartier, et des mesures à quelques moments précis comme le matin, le midi, l'après-midi et la nuit susent généralement à décrire la distribution spatiale et l'évolution de la température.

En conclusion, il est bénéque de tirer le meilleur parti des connaissances existantes en intégrant la compréhension de la physique et les approches basées sur les données dans les projets d'ingénierie.

Réduction de modèle : physique en temps réel

Les techniques visant à améliorer la résolution des modèles basés sur la physique an d'obtenir des solutions en temps réel appartiennent au domaine de la réduction de modèle (Model Order Reduction -MOR) [13].

Il en existe deux déclinaisons :

La première consiste à construire un modèle d'ordre réduit, dit intrusif, en appliquant des techniques de réduction de la dimensionnalité à la solution du modèle physique. Ce mécanisme est intégré au processus de résolution, et permet au solveur de fonctionner rapidement.

La seconde solution consiste à construire un modèle de substitution, qui établit une relation directe entre la solution du modèle physique et les paramètres du modèle. Une fois le modèle de substitution construit, il permet de déduire la solution quasiment en temps réel lorsque les paramètres du modèle sont fournis.

Les modèles d'ordre réduit et les solutions paramétriques jouent un rôle crucial dans la réalisation de diverses tâches, notamment la simulation, l'optimisation, l'analyse inverse, le contrôle basé sur la simulation et la propagation d'incertitudes, tout en respectant des contraintes strictes en matière de temps réel. Ces solutions paramétriques orent un moyen exceptionnellement ecace et sans précédent pour eectuer de l'ingénierie basée sur des modèles en temps réel.

xii Réduction de modèle avancée et approches basées sur les données pour la construction de jumeaux numériques augmentés par la physique

Données et apprentissage

Les techniques de régression, qui permettent d'extraire des modèles des données, peuvent être appliquées à diérents types de données. Parmi les types de données les plus courants, on peut citer :

Listes (tableaux) : données organisées en tableaux, constitués de valeurs numériques continues ou discrètes, et parfois des caractéristiques catégorielles ;

Images : images réelles en 2D ou 3D ou résultats de simulations numériques ;

Formats tensoriels : images compressées ou données de champs spatiaux décomposées en produits tensoriels ;

Graphes : caractéristiques de n÷uds et d'arêtes représentées par une structure de graphe ;

Courbes : courbes géométriques ou fonctions à valeur scalaires ;

Séries temporelles : données ordonnées chronologiquement impliquant une causalité ;

Nuages de points : densités ou distributions de probabilité.

Plusieurs techniques d'apprentissage automatique sont couramment utilisées pour la régression :

Régressions polynomiales régularisées [START_REF] Wan | The eect of regularization coecient on polynomial regression[END_REF] : Elles sont utiles lorsque les données sont limitées et que la solution peut être exprimée à l'aide de bases polynomiales. Elles peuvent intégrer des connaissances existantes par la construction de variables adéquates.

Réseaux neuronaux articiels (NN) : Ils donnent de bons résultats lorsque les données sont abondantes et que les hyperparamètres sont bien réglés. Ils peuvent approximer n'importe quelle fonction [START_REF] Gorban | The general approximation theorem[END_REF], ce qui les rend extrêmement polyvalents.

Auto-encodeurs (AE) [START_REF] Tschannen | Recent advances in autoencoderbased representation learning[END_REF] : Ils encodent les données dans un espace latent, en approximant la dimensionnalité intrinsèque et en supprimant les corrélations.

D'autres techniques de réduction de la dimensionnalité telles que l'ACP [START_REF] Ringnér | What is principal component analysis?[END_REF], la kPCA [START_REF] Schölkopf | Kernel principal component analysis[END_REF], le LLE [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] et le tSNE [START_REF] Hinton | Stochastic neighbor embedding[END_REF] peuvent également être utilisées.

Réseaux antagonistes génératifs (GAN) [START_REF] Creswell | Generative adversarial networks: An overview[END_REF] : Ils génèrent des échantillons de données synthétiques permettant d'augmenter l'ensemble de données, ou de générer des échantillons réalistes.

Réseaux neuronaux convolutifs (CNN) [START_REF] Aloysius | A review on deep convolutional neural networks[END_REF] : Ils sont bien adaptés à l'analyse d'images et aux tâches de reconnaissance et extraction de motifs. Les réseaux neuronaux en graphes (GNN) [START_REF] Zhou | Graph neural networks: A review of methods and applications[END_REF] étendent les CNN pour qu'ils fonctionnent sur des maillages ou des graphes non structurés.

Réseaux neuronaux récurrents (RNN) [START_REF] Medsker | Recurrent neural networks[END_REF] et long short-term memory (LSTM) [START_REF] Hochreiter | Long short-term memory[END_REF] : Puissants pour traiter les données de séries temporelles. Ils peuvent capturer les dépendances temporelles grâce à des capacités de mémoire.

Réduction de modèle avancée et approches basées sur les données pour la construction de jumeaux numériques augmentés par la physique xiii Systèmes dynamiques : Diverses techniques sont utilisées pour modéliser les systèmes dynamiques, notamment les RNN, les LSTM, ResNet [START_REF] Elhassouny | Trends in deep convolutional neural networks architectures: A review[END_REF], NeuralODE [START_REF] Chen | Neural ordinary dierential equations[END_REF], DMD [START_REF] Tu | Dynamic mode decomposition: Theory and applications[END_REF], l'opérateur de Koopman [START_REF] Mauroy | Koopman operator in systems and control[END_REF], le calcul par réservoir [START_REF] Nakajima | Reservoir computing[END_REF] et les DeepONets [START_REF] Lu | Learning nonlinear operators via deeponet based on the universal approximation theorem of operators[END_REF].

Réseaux neuronaux informés par la physique (PINN) [13], réseaux neuronaux informés par la thermodynamique (TINN) [START_REF] Hernandez | Deep learning of thermodynamics-aware reduced-order models from data[END_REF] et modèles hybrides : Ils intègrent des connaissances fondées sur la physique dans le processus d'apprentissage, garantissant ainsi la satisfaction des équations fondamentales et combinent des approches fondées sur la physique et des approches fondées sur les données.

Ces techniques orent diverses options pour les tâches de régression, en fonction de la nature des données et des exigences de modélisation souhaitées.

Modélisation hybride

Dans le contexte de la modélisation hybride [START_REF] Glassey | Hybrid modeling in process industries[END_REF], il existe deux approches : l'enrichissement de la solution et l'enrichissement du modèle.

Dans le premier cas, la solution du modèle basée sur la physique, obtenue par des techniques spéciques, est combinée avec une correction basée sur les données. Le modèle est calibré en ligne en identiant les paramètres opérationnels qui minimisent la diérence entre les prédictions du modèle et les mesures disponibles. La prédiction hybride est la combinaison de la solution calibrée basée sur la physique et de la correction basée sur les données. Cette procédure s'applique également aux problèmes transitoires.

Dans le second cas, l'objectif est d'améliorer la capacité prédictive d'un modèle nominal connu. Les mesures sont comparées aux prédictions du modèle nominal et un écart notable est observé. Pour y remédier, un terme d'enrichissement du modèle est introduit an de mieux représenter les données collectées. Le modèle discret enrichi, ainsi que la prédiction associée, doivent satisfaire à la condition d'équilibre tout en représentant précisément les données. Une paramétrisation de la correction du modèle est choisie, ainsi qu'une technique de régularisation appropriée, pour calculer l'enrichissement du modèle et compléter les données de manière ecace.

Contenu du mémoire

Ce mémoire est divisée en trois chapitres : tout d'abord, quelques outils et méthodologies qui sont des éléments clés dans la construction de jumeaux hybrides sont décrits, puis les contributions personnelles dans le domaine de la réduction non intrusive de l'ordre des modèles sont présentées et enn, les contributions personnelles à la modélisation hybride sont décrites.

Chapitre 1 : Revisite de méthodes permettant la construction de jumeaux hybrides Section 1 : Méthodes de réduction de modèle intrusives et nonintrusives Les problèmes d'ingénierie nécessitent souvent des simulations rapides et précises utilisant des modèles complexes et de grandes quantités de données. Bien que les dévelop-xiv Réduction de modèle avancée et approches basées sur les données pour la construction de jumeaux numériques augmentés par la physique pements du hardware aient permis d'améliorer la qualité des simulations, la résolution de problèmes en grande dimension reste dicile en raison de la complexité algorithmique et des contraintes de traitement en temps réel. Les techniques de réduction de modèle (MOR), qui peuvent être intrusives ou non-intrusives, ont attiré l'attention pour résoudre ces problèmes.

La décomposition orthogonale aux valeurs propres (POD) [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF][START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent ows[END_REF] extrait les caractéristiques essentielles du système et crée une base d'approximation réduite, ce qui diminue considérablement le temps de calcul. La Proper Generalized Decomposition (PGD) [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF][START_REF] Ladevèze | The latin multiscale computational method and the proper generalized decomposition[END_REF] est une méthode intrusive qui utilise la séparation de variables pour réduire la complexité des algorithmes de résolution. La décomposition orthogonale aux valeurs propres avec interpolation (PODI) [START_REF] Tezzele | Shape optimization through proper orthogonal decomposition with interpolation and dynamic mode decomposition enhanced by active subspaces[END_REF][START_REF] Rama | Proper orthogonal decomposition with interpolation-based real-time modelling of the heart[END_REF] est une méthode non intrusive qui construit des solutions paramètriques à partir de la base réduite sans nécessiter de modier le solveur. La sparse PGD (sPGD) [START_REF] Ibáñez | A multidimensional datadriven sparse identication technique: the sparse proper generalized decomposition[END_REF] construit une solution paramétrique à l'aide d'un ensemble peu dense de solutions calculées par un solveur d'équations aux dérivées partielles en utilisant la séparation de variables pour traiter les problèmes avec beaucoup de paramètres.

Dans l'ensemble, les techniques de réduction de modèle fournissent des solutions ecaces pour des problèmes d'ingénierie complexes nécessitant une application en temps réel.

Section 2 : Réseaux neuronaux informés par la physique Les réseaux neuronaux informés par la physique (PINN) [13] sont apparus comme un cadre prometteur qui intègre dans les réseaux neuronaux les principes de la physique pour résoudre des problèmes scientiques et d'ingénierie complexes. Contrairement à l'apprentissage automatique traditionnel, les PINN intègrent des connaissances préalables sur les lois physiques an de permettre des prédictions précises et de saisir la physique sous-jacente, même avec un nombre limité de données annotées.

Par ailleurs, ces architectures sont un cadre propice aux stratégies hybrides. En eet, il est naturel d'ajouter dans la fonction de perte du réseau neuronal toutes les contraintes que doit respecter la solution. Ces contraintes peuvent provenir de lois physiques, souvent sous la forme d'équations aux dérivées partielles, mais aussi de données mesurées ou encore de mécanismes de régularisation. Pendant l'entraînement, la solution est guidée par le gradient de la fonction de perte vers un compromis entre toutes les contraintes.

Il est important de remarquer que combiner un modèle physique totalement déterminé avec des données expérimentales ne donne pas réellement un modèle hybride mais une solution hybride, qui est très utile dans de nombreux cas d'utilisation mais n'apporte pas de possibilité d'explicabilité ou de généralisation à des problèmes similaires.

Chapitre 2 : Progrès dans la réduction de modèle non intrusive Section 1 : Regularisation de représentations séparées Le dé de la régression avec des données limitées consiste à traiter les non-linéarités dans des contextes multiparamétriques, ce qui entraîne la malédiction de la dimensionnalité [START_REF] Köppen | The curse of dimensionality[END_REF][START_REF] Verleysen | The curse of dimensionality in data mining and time series prediction[END_REF]. Le principe du rasoir d'Occam [START_REF] Udrescu | Ai feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity[END_REF][START_REF] Brunton | Discovering governing equations from data by sparse identication of nonlinear dynamical systems[END_REF] suggère de privilégier les modèles les plus simples aux modèles les plus complexes, ce qui est possible grâce à la Engineering in the 20th century primarily relied on the utilization of models to eectively design components and systems. Generally, these models comprised mathematical operators that facilitated the transmission of input to its corresponding output.

To provide a clearer explanation, let's consider a deformable solid that undergoes shape modication due to the application of a force on its surface. This deformation results in an internal mechanical condition within the solid, leading to reversible or irreversible changes in its geometry. Such changes are often observed in manufacturing processes like stamping or forging, among others. Hence, the mechanical model mentioned earlier can be perceived as the operator that deduces the nal mechanical state in the structural system based on a given load.

There are two main typologies of models, referred later as Type I and Type II:

1. The rst one (Type I) expresses the input/output relationship in engineering in an algebraic form, whether linear or nonlinear. For instance, consider the example of a rod experiencing tension. In this case, the applied stress on the rod is directly proportional to the resulting deformation, which is dened as the relative change in its length. This relationship remains linear for small deformations but can transition to a nonlinear relationship for larger deformations.

Furthermore, for even larger deformations, it may become irreversible, leading to inelastic behavior.

These algebraic models are easily manipulated and oer signicant advantages in larger multi-component systems. Their algebraic nature greatly facilitates ecient manipulation and problem-solving.

2. The second type (Type II) involves an input/output relationship that gives rise to a mathematical problem involving derivatives in space and/or time. This situation arises when the solution at a particular point and time depends not only on the present conditions but also on the past history and the solution in all other points within the relevant domain. In such cases, models are expressed using more complex mathematical entities known as partial dierential equations.

These equations capture the interdependence of variables and their derivatives, enabling a comprehensive representation of the system dynamics.

Models of Type II generally pose greater challenges when it comes to nding ecient solutions. To illustrate this, let's consider a one-dimensional model, such as a rod, where we aim to determine the internal state at all points within it. Initially, we 2 0.1 Context and Motivations encounter a diculty: there is an innite number of unknowns, corresponding to the numerous points along the rod.

However, if the problem can be solved analytically, this challenge can be overcome. In such cases, the solution takes the form of a closed mathematical expression, such as the displacement u(x) in the example mentioned. This displacement at each point x in the rod can then be determined based on the forces f (x ′ ) applied at each corresponding point x ′ in the rod, considering the given boundary conditions.

It is important to note that most models cannot be solved analytically due to various reasons. They are often dened in complex geometries, involve intricate boundary conditions, and may exhibit strong nonlinearity and coupling eects. As a result, discretization techniques have emerged as prominent tools in Simulation-Based Engineering (SBE) [18,19]. Examples of such techniques [START_REF] Peiró | Finite dierence, nite element and nite volume methods for partial dierential equations[END_REF][START_REF] Mattiussi | An analysis of nite volume, nite element, and nite dierence methods using some concepts from algebraic topology[END_REF] include nite elements [START_REF] Zienkiewicz | The nite element method in engineering science[END_REF][START_REF] Bathe | Finite element method[END_REF],

nite dierences [START_REF] Jordán | Calculus of nite dierences[END_REF][START_REF] Eh | The calculus of nite dierences[END_REF], and nite volumes [START_REF] Versteeg | An introduction to computational uid dynamics: the nite volume method[END_REF][START_REF] Eymard | Finite volume methods, Handbook of numerical analysis[END_REF].

These discretization techniques aim to nd solutions at a limited number of points, referred to as nodes, and specic time instances, rather than attempting to solve the problem for an innite number of points and times. They eectively transform the continuous model into a discrete representation. Among these techniques, the nite element method has gained widespread popularity and is widely adopted in industrial practices today. It has proven to be a valuable approach for solving complex engineering problems.

However, the journey from the conceptual formalization of these discretization techniques to their widespread adoption in the industry has been quite extensive. Taking the previously mentioned one-dimensional problem as an example, the chosen discretization technique transforms the initial problem into a task of calculating the solution at only N points (referred to as nodes) distributed along the rod, forming a grid or mesh.

To illustrate, let us consider the case where N = 10. Once the solution is obtained at these 10 points, it can be interpolated to estimate the solution at any other location within the domain. However, in two dimensions (such as a square domain), the grid would contain 10x10 points that need to be calculated, and in three dimensions (e.g., a cube), it would involve 10x10x10 points. As a result, the algebraic problem associated with one-dimensional discretization has 10 unknowns, while two-dimensional and three-dimensional discretization problems have 100 and 1000 unknowns, respectively.

The number of unknowns scales with the number of nodes in the grid or mesh, which increases exponentially with the dimension of the space where the problem is dened.

Even if an algebraic problem of size 10 can be solved relatively quickly by hand, the one involving the solution of a problem with 1000 unknowns could require a signicant amount of time, potentially exceeding the lifespan of a human being.

Discretization, which is necessary for solving problems involving Type II models, became feasible thanks to the advent and contribution of computers. While computers can solve algebraic systems of equations quickly, problems with larger dimensions pose signicant challenges. For example, solving an algebraic problem with 1000 unknowns may only take fractions of a second on a standard laptop. However, engineering problems are increasingly growing in size and complexity. They often involve hundreds of millions of unknowns and require many iterations due to nonlinearity or transient behavior (evolving solutions over time).

Even with the most powerful computational platforms available today, simulations of such large-scale problems still demand substantial time. They may take days, weeks, or even months to complete. This extended computational time is necessary to accurately capture the intricate details of the system and obtain reliable results.

In summary, while most models can be solved with the aid of powerful computers, they often require a signicant amount of time to perform the necessary computations.

Engineering practices in the 20th century primarily focused on component design, and industries such as energy, automotive, and aerospace employed programs long enough to ensure the dialog between designers and simulation tools.

It is important to emphasize that the design process often necessitates multiple simulations to explore the parametric space and identify a quasi-optimal design. However, the cost associated with these simulations and the dimensionality of the parametric space, which can consist of hundreds of parameters, present ongoing challenges in achieving optimal certied designs even in modern times.

Another challenging aspect of Simulation-Based Engineering (SBE) relates to the quality of the models themselves. While models have evolved over centuries of scientic advancements, there are several factors that make it dicult to fully grasp and express the true behavior of a system:

The size of the analyzed system in terms of space and time poses the rst challenge. A model may be considered generally accurate, but it can exhibit noticeable local deviations. Similarly, a model that accurately represents the system at present may develop biases over time.

Another diculty lies in the large number of parameters involved in some models, making it challenging to calibrate and identify the appropriate values for these parameters.

Additionally, even when using the correct model, the parameters themselves exhibit variability. For example, specimens made of the same material from the same supplier, produced by the same machine on the same day, and processed under the same nominal conditions can still exhibit slight dierences in their mechanical properties. This leads to statistical distributions of model parameters, introducing uncertainty that propagates throughout the system.

Lastly, there is always an inherent epistemic uncertainty or ignorance between the model and reality. Dierences can exist at various scales in space and time, and these dierences can potentially grow as time progresses.

In practice, engineers strive to mitigate the impact of these challenges to ensure accurate predictions and maintain the quality of their designs. Various techniques and strategies are employed to address these issues and minimize the uncertainties introduced by the models used in SBE.

From virtual to digital twins

The 20th century engineering was prolic, with many impressive successes [START_REF] Constable | A century of innovation: Twenty engineering achievements that transformed our lives[END_REF]: aeronautics, space conquest, transport, energy, civil infrastructures, ...

In other domains such as marketing, economy, social sciences, and medicine, where traditional models were less established and had limited prognostic capabilities, researchers began shifting their focus towards the use of data for diagnosis and prognosis. This approach quickly led to signicant achievements. However, in the eld of engineering, three major new developments disrupted the status quo: 0.1 Context and Motivations 3. The upcoming human-centric immersive engineering revolution: The imminent fth revolution will place humans at the core of engineering endeavors. It envisions a physics-aware metaverse, where advancements in computer vision, multimodal devices, and natural language processing (NLP) have the potential to defy the Turing test. However, integrating humans into the system poses challenges for traditional modeling approaches, as humans may escape or resist conventional modeling frameworks.

In response to these developments, engineering science turned its attention to the promising possibilities oered by data and Articial Intelligence (AI). The expectations for leveraging these technologies in engineering are high, with the aim of achieving better predictions, improved performance, and enhanced understanding of complex systems [START_REF] Madni | Leveraging digital twin technology in model-based systems engineering[END_REF].

The learning process in machine learning seems general and agnostic to the specic physics of the problem. Once input and output data (represented as (p i , u i ), i = 1, . . . , n s , n s being the number of available data) are collected, one can attempt to establish the relationship between the inputs (p) and outputs (u), denoted as u = F(p).

By learning and establishing the functional relation F from the available data, the model can be applied to new input data to obtain the corresponding output. This allows for real-time application of the learned model to dierent scenarios.

The main challenge in this process is to obtain the regression model u = F(p) that accurately predicts the output u given the input p. There are various techniques available for constructing such regression models, and the choice depends on factors such as the quantity of available data (n s ), the time available for training (online versus oine), and the quality of the data.

In summary, if the quantity of interest (u) and the variables aecting it (p) are well identied, and if there is a sucient quantity (n s ) and quality of data, there are nu- merous machine learning techniques capable of extracting the functional relationship (regression) u = F(p).

Imagine for a while the set of data: (u 1 = 2, p 1 = 1), (u 2 = 4, p 2 = 2), ... (u n = 2n, p n = n), n = 3, . . . , n s . In that case the regression becomes linear and almost obvious: u = 2p, i.e. F(p) = 2p.

While data-driven modeling has its advantages, there are still challenges and diculties that persist in the process:

Data collection: Determining which data to collect, at what scale, and when and where to collect it is crucial. The collected data should be relevant and contribute to the desired output. It is important not to discard any features that may be useful. Techniques exist to identify and remove irrelevant features and discover missing features, but collecting data can be costly and constrained by various factors such as available resources, technological limitations, regulations, and environmental considerations [4,5].

Cost of data: Data collection and processing can be expensive in engineering.

Unlike in some domains where big data is readily available, engineering often faces limitations in data collection. The cost includes expenses related to sensing devices, instrumentation, data communication, data mining, processing, and addressing environmental concerns. The availability and cost of data can impact the feasibility and practicality of certain approaches [6,7].

Interpretability: While machine learning models can extract regressions and make predictions, explaining the learned model can be challenging. Understanding the underlying factors and relationships that the model has learned from the data may not be straightforward. It is easier to explain well-established scientic concepts like gravity, but explaining data-driven models may require additional eort and analysis [10,11].

Extrapolation and generalization: Models learned from data can be applied within the domain of validity of the collected data. When applying the learned model to situations or regions where data was not collected, extrapolation occurs. Extrapolation carries inherent risks as the model may not accurately predict outcomes in unobserved or extreme conditions. Care must be taken to understand the limitations and potential risks associated with extrapolating the model beyond its observed data range [8,9].

Addressing these challenges requires careful consideration of data collection strategies, cost-eectiveness, interpretability of models, and the appropriate use of models within their valid domains.

In between models and data

As discussed earlier, the option that relies primarily on physics-based approaches has its advantages, but it also has drawbacks such as longer computing time and the need for signicant computational resources. Additionally, it becomes challenging to handle large-scale systems with notable variability and uncertainty. On the contrary, the almost fully data-driven option appears promising, especially when accurate models are unavailable or lacking. However, this approach also faces diculties when applied in engineering practices, including high data costs, risks associated with extrapolation, and challenges in providing explanations and certications.

Perhaps the optimal choice lies in combining both approaches instead of choosing one over the other. Within this context, we will explore two possibilities:

1. Physics-informed learning [START_REF] Karniadakis | Physics-informed machine learning[END_REF],

0.1 Context and Motivations 2. Physics-augmented learning [START_REF] Liu | Physics-augmented learning: A new paradigm beyond physics-informed learning[END_REF].

In both cases, we leverage the advantage that learning specic parts can be more costeective than learning the entire system. In the realm of physics, it is widely accepted that conservation laws are universal, and therefore, there is no need to relearn them from data. However, there are other phenomenological relationships that can benet from a data-driven perspective, oering an opportunity for improvement.

Physics-informed learning

Let us consider a scenario where we are seeking to determine the output u at each point x of a domain Ω, i.e. u(x), x ∈ Ω. The regression model used to express u(x)

is expected to be complex and exhibit strong nonlinearity. Instead of employing the nite element method discussed earlier, we opt to approximate u using a nonlinear machine-learning-based regression, such as a neural network (NN).

However, according to established physics principles that we assume to be applicable here, u(x) adheres to a state-of-the-art model dened by a partial dierential equation in the form: L(u(x)) = f (x). In this equation, L(•) represents a linear or nonlinear dierential operator.

To address this, we turn to a physics-informed neural network (PINN), which seeks to nd a neural network representation of u from x that approximates u(x) while ensuring adherence to its governing equation. This is achieved by minimizing the residual R = ∥L(u(x))-f (x)∥, which enforces the verication of the partial dierential equation during the construction of the regression for u(x).

PINN can be seen as a type of collocation method, where the unknown eld, u(x), is approximated by a neural network, and the construction of the regression is based on minimizing the residual. This approach allows for easy assimilation of data in a straightforward and transparent manner.

Physics-augmented learning

The main concept here revolves around the assumption that the reference solution u(x), obtained from a given loading term f (x), can be decomposed into two contributions. The rst contribution represents our existing knowledge, which originates from the physics-based model and is denoted as u P (x).

Subsequently, we calculate the discrepancy between the actual behavior and the prediction of the physics-based model, u(x) -u P (x). This discrepancy is commonly referred to as "ignorance."

To address this gap, machine learning techniques are employed to model this discrepancy, resulting in the data-driven model denoted as u D (x). Consequently, the fundamental approximation is expressed as: u(x) = u P (x) + u D (x).

This relationship forms the basis of the hybrid modeling paradigm or the hybrid twin approach when applied to specic systems or assets.

Discussion

Both approaches, the physics-based model and the data-driven model, oer the advantage of signicantly reducing the amount of data required for constructing the models. This reduction is achieved either by leveraging the provided information or by focusing solely on modeling the discrepancy, which is assumed to be simpler to approximate than the actual behavior itself.

Furthermore, the use of the physics-based model enables us to explain the foundational aspects of the model and its resulting predictions, thereby alleviating extrapolation issues to some extent.

The reduction in data requirements is not solely due to the utilization of physics but also because physics informs us about the optimal locations and timing for data collection. This is particularly evident when considering active learning, which becomes highly eective when incorporating existing physics-based knowledge.

To illustrate this point, let us consider the example of temperature measurement in a city. To understand the temperature in a city, it is unnecessary to place thousands of thermometers on every street and record the temperature every millisecond. With knowledge of the physics that shapes our experience and common sense, we know that placing a few thermometers (e.g., one in each district) and measuring the temperature in the morning, at noon, in the afternoon, and at midnight would generally suce.

In general, we should make the most of the existing knowledge available to us, incorporating both physics-based understanding and data-driven approaches.

Physics in real-time: Model Order Reduction

In this section we focus on the physics-based model, and tackle the methodologies to enhance its solution procedure to obtain real-time performance. Such procedures are known as Model Order Reduction (MOR) [13].

For this purpose, we consider two alternatives:

The construction of an intrusive Reduced Order Model by applying dimensionality reduction techniques to the solution of the physics-based model, which subsequently enables the solver to operate very quickly.

The construction of a surrogate, which refers to a direct relation between the solution of the physics-based model and the parameters involved in the model.

Once this surrogate is built, it enables us to infer the solution almost in real-time by simply providing the model parameters.

Reduced order models and parametric solutions play a crucial role in enabling various tasks, such as simulation, optimization, inverse analysis, simulation-based control, and uncertainty propagation, all while operating under strict real-time constraints.

These parametric solutions facilitate real-time engineering based on models in an exceptionally ecient and unprecedented manner.

To illustrate the concept, let's consider a simple parametric solution denoted as u(x; p 1 , p 2 ), where the unknown u depends on two parameters, p 1 and p 2 , at each position x. It's important to note that while we are using two parameters for simplicity, the approach can be extended to include more parameters as needed.

The construction of the surrogate follows these steps:

1. A Design of Experiments (DoE) is dened, employing a sampling strategy. Several possibilities exist, such as Latin Hypercube [START_REF] Loh | On latin hypercube sampling[END_REF], quadrature methods (e.g.,

Context and Motivations

Gauss-Lobatto-Chebyshev [START_REF] Herman | Direct optimization using collocation based on high-order gauss-lobatto quadrature rules[END_REF] or Smolyak [START_REF] Smolyak | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF], or approaches associated with active learning [12] (such as Gaussian processes [START_REF] Williams | Gaussian processes for regression[END_REF] or Fisher information matrix [START_REF] Fisher | The statistical utilization of multiple measurements[END_REF]). The DoE provides a sampling of the parametric space: (p 1 1 , p 2 1 ), . . . , (p 1 D , p 2 D ).

2. Using the physics-based model and suitable software capable of solving it, the solution associated with each point in the DoE is calculated: ResNet [START_REF] Elhassouny | Trends in deep convolutional neural networks architectures: A review[END_REF] or NeuralODE [START_REF] Chen | Neural ordinary dierential equations[END_REF] can learn the forcing terms of dynamical systems.

u 1 (x) = u(x; p 1 1 , p 2
Techniques like Dynamic Mode Decomposition (DMD) [START_REF] Tu | Dynamic mode decomposition: Theory and applications[END_REF] and the Koopman operator [START_REF] Mauroy | Koopman operator in systems and control[END_REF] are ecient for addressing nonlinear behaviors. Reservoir computing [START_REF] Nakajima | Reservoir computing[END_REF] and DeepONets [START_REF] Lu | Learning nonlinear operators via deeponet based on the universal approximation theorem of operators[END_REF] are emerging techniques in this eld.

Physics-Informed Neural Networks (PINN) [13], Thermodynamics-Informed Neural Networks (TINN) [START_REF] Hernandez | Deep learning of thermodynamics-aware reduced-order models from data[END_REF] and Hybrid Models: These techniques incorporate physics-based knowledge into the learning process. PINNs, TINNs, and SPNNs [START_REF] Hernández | Structure-preserving neural networks[END_REF] ensure the satisfaction of governing equations. Hybrid models combine physics-based and data-driven approaches, often used in transfer learning scenarios.

Concurrently, learning modalities play a crucial role in machine learning. Here are some of the commonly used modalities:

Supervised Learning [START_REF] Nasteski | An overview of the supervised machine learning methods[END_REF]: In supervised learning, the model learns from labeled data, where the input samples are associated with corresponding target labels.

The model aims to learn the mapping between inputs and outputs based on the provided labeled examples.

Unsupervised Learning [START_REF] Celebi | Unsupervised learning algorithms[END_REF]: Unsupervised learning involves learning from unlabeled data. The model explores the inherent structure or patterns in the data without any specic target labels. Clustering, dimensionality reduction, and generative models are examples of unsupervised learning techniques.

Self-supervised Learning [START_REF] Liu | Self-supervised learning: Generative or contrastive[END_REF]: Self-supervised learning is a variation of unsupervised learning where the model learns to predict certain aspects of the data from the data itself. It formulates learning tasks using the available unlabeled data, creating proxy supervisory signals.

Semi-supervised Learning [START_REF] Reddy | Semi-supervised learning: A brief review[END_REF]: In semi-supervised learning, the model learns from a combination of labeled and unlabeled data. It leverages the limited labeled data and the additional unlabeled data to improve the learning process.

This approach is benecial when labeled data is expensive or time-consuming to obtain.

Context and Motivations

Transfer Learning [START_REF] Zhuang | A comprehensive survey on transfer learning[END_REF]: Transfer learning involves leveraging knowledge or models learned from one domain and applying it to another related domain. The pre-trained models on large-scale datasets can be ne-tuned or used as a feature extractor for a target task with limited data.

Reinforcement Learning [START_REF] Li | Deep reinforcement learning: An overview[END_REF]: Reinforcement learning is used when an agent interacts with an environment and learns to take actions that maximize a reward signal. The agent explores the environment and learns through trial and error, receiving feedback in the form of rewards or penalties.

The choice of learning modality depends on factors such as the availability of labeled data, the complexity of the data structure, and the need for online or real-time learning. Each modality has its strengths and limitations, and the selection depends on the specic problem and context.

Hybridization

In the context of hybridization [START_REF] Glassey | Hybrid modeling in process industries[END_REF], there are two main routes that can be followed:

enriching the solution or enriching the model and thus obtaining the enhanced solution.

Solution enrichment

The enriched or corrected solution is obtained by combining the physics-based model solution, u P (x, p) (computed by using the techniques discussed in Section 0.1.3), with the data-driven correction, u D (x, p) (learned by using one of the technologies introduced in Section 0.1.4). Here, p represents the vector of parameters associated with the model.

To calibrate the model online, as soon as data is collected, the operational parameters p * are identied by minimizing the dierence between the model predictions and the available measurements. This calibration process ensures that the model is aligned with the observed data and can accurately capture the behavior of the system.

Finally, and then the hybrid prediction yields the corrected or enriched solution from

u P (x, p * ) + u D (x, p * ).
The same procedure applies to transient problems, where the model parameters may vary over time. The online calibration is performed at each time step, and the datadriven model is particularized accordingly to provide an updated hybrid prediction.

Model enrichment

In the model enrichment route, the objective is to improve the predictive capability of a nominal model that is assumed to be known. The discrete form of the nominal model, for simplicity, can be represented as MU = F, where M is the matrix representing the system properties, U and F are the nodal vectors of displacements and forces, respectively.

In practice, measurements of displacements, denoted as Ũexp , are collected at certain locations. When comparing these measured displacements with the predictions of the nominal model in the same locations, denoted as Ũ, a noticeable discrepancy, or gap, is observed. The norm of the dierence between the measured and predicted displacements ∥ Ũexp -Ũ∥ remains larger than the acceptable error threshold ϵ, i.e. ∥ Ũexp -Ũ∥ > ϵ.

To address this discrepancy, an enrichment matrix M * is introduced. The goal is to nd an enriched discrete model, represented as M + M * , that can better represent the collected data. The corrected model, along with the associated displacement prediction, denoted as U+U * , should satisfy the equilibrium condition (M+M * )(U+ U * ) = F while accurately representing the collected data, i.e. ∥ Ũexp -( Ũ + Ũ * )∥ < ϵ.

To compute the model enrichment and perform data completion U+U * , a parametrization of the model correction M * is chosen, along with an appropriate regularization technique. This parametrization allows for the representation of the model correction in a computationally ecient manner.

Objectives of the thesis

The objective of the work performed in this thesis is to apply, develop, and investigate novel tools that enable the use of virtual, digital, and hybrid twins in order to tackle present-day industrial challenges. Specically, the focus has been placed on exploring and expanding the Model Order Reduction (MOR) and Machine Learning (ML) frameworks and the combination of both. The contributions tackle three major diculties in the context of Hybrid Twins:

Performing model reduction and regression at the low-data limit.

Employing adequate mathematical representations to empower model reduction.

Collecting relevant measurements for data-driven and hybrid modelling.

Using data to enrich physics-based models.

In particular, the dierent topics addressed concern:

Contributions to the development of machine learning and non-intrusive reduced modeling techniques [9395].

Development of general MOR methods which address the localization of features [START_REF] Champaney | Parametric curves metamodelling based on data clustering, data alignment, pod-based modes extraction and pgd-based nonlinear regressions[END_REF][START_REF] Torregrosa | Surrogate parametric metamodel based on optimal transport[END_REF].

Application of MOR and ML methods to industrial cases [98100].

Adaption of MOR techniques to solve dicult problems [START_REF] Vermiglio | Parametric electromagnetic analysis of radar-based advanced driver assistant systems[END_REF][START_REF] Jacot | Parametric damage mechanics empowering structural health monitoring of 3d woven composites[END_REF].

Development of tools to reduce component models in engineering applications [START_REF] Loreau | Learning the parametric transfer function of unitary operations for real-time evaluation of manufacturing processes involving operations sequencing[END_REF][START_REF] Reille | Learning data-driven reduced elastic and inelastic models of spot-welded patches[END_REF].

Employment of adequate representations and metrics to apply machine learning to materials science [START_REF] Loreau | Parametric analysis and machine learningbased parametric modeling of wire laser metal deposition induced porosity[END_REF][START_REF] Runacher | Describing and modeling rough composites surfaces by using topological data analysis and fractional brownian motion[END_REF].

Theoretical study on the feasibility of modeling from partial observations [START_REF] Champaney | Modeling systems from partial observations[END_REF].

Development of hybrid modeling techniques [START_REF] Kestelyn | Towards a hybrid twin for infrastructure asset management: Investigation on power transformer asset maintenance management[END_REF][START_REF] Torregrosa | Hybrid twins based on optimal transport[END_REF].

Application of hybrid modeling to model correction and damage identication [START_REF] Lorenzo | Data completion, model correction and enrichment based on sparse identication and data assimilation[END_REF][START_REF] Lorenzo | Physics informed and data-based augmented learning in structural health diagnosis[END_REF].

Construction of Digital and Hybrid Twin methodologies [START_REF] Champaney | Hybrid twins-a highway towards a performance-based engineering. part i: Advanced model order reduction enabling real-time physics[END_REF][START_REF] Champaney | Engineering empowered by physics-based and data-driven hybrid models: A methodological overview[END_REF]. Physics-Informed Neural Networks (PINN) [13] have emerged as a promising framework that integrates neural networks with physics principles to solve complex scientic and engineering problems. Unlike traditional machine learning, PINNs incorporate prior knowledge of physical laws to enable accurate predictions and capture underlying physics, even with limited labeled data. This section introduces PINN, explaining how it combines governing equations and constraints with a neural network architecture.

Chapter 2: Advances in non-intrusive Model Order Reduction Section 1: Regularization of separated representations

The challenge in regression with limited data is handling nonlinearities in multiparametric settings, leading to the curse of dimensionality [START_REF] Köppen | The curse of dimensionality[END_REF][START_REF] Verleysen | The curse of dimensionality in data mining and time series prediction[END_REF]. Occam's razor [START_REF] Udrescu | Ai feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity[END_REF][START_REF] Brunton | Discovering governing equations from data by sparse identication of nonlinear dynamical systems[END_REF] principle suggests favoring simpler models over complex ones, achieved through sparsity promotion in regression [START_REF] Brunton | Discovering governing equations from data by sparse identication of nonlinear dynamical systems[END_REF][START_REF] Ibañez | Some applications of compressed sensing in computational mechanics: model order reduction, manifold learning, data-driven applications and nonlinear dimensionality reduction[END_REF]. This section proposes robust and efcient regression methodologies for separated representation settings. Two techniques are presented: rs-PGD, combining L2 and L1 norm regularization, and s 2 -PGD, a doubly sparse regression method utilizing Lasso regularization. Predicting whole curves from given input features in parametric surrogates presents challenges. Non-alignement in curves, which can result from dierent sampling or ending points, and shifted critical points in curves with common patterns (e.g., elasticplastic transition in mechanics), make classical interpolation methods fail to provide physics-consistent results, requiring appropriate pre-processing. When bifurcations occur in the parametric space, coupling with clustering [START_REF] Nielsen | Hierarchical clustering, Introduction to HPC with MPI for Data Science[END_REF][START_REF] Murtagh | Algorithms for hierarchical clustering: an overview[END_REF] and classication [START_REF] Kotsiantis | Supervised machine learning: A review of classication techniques[END_REF] algorithms enhances the surrogate's accuracy. This section introduces methodologies to address these issues and create surrogates used for uncertainty quantication [START_REF] Lee | A comparative study of uncertainty propagation methods for black-box-type problems[END_REF],

providing statistical bounds for predicted curves. Section 3: Hybrid modeling by identifying properties

In the pursuit of rening physics-based models through hybrid modeling, an alternative approach arises: correcting uncertain physical parameters associated with the system. This correction enables the development of a model that accurately captures the behavior of a damaged system and identies the nature and location of the damage. This section presents a novel methodology focused on parameter correction, rening the model's predictive capabilities and providing insights into the damage.

Sparse regularization techniques are key to this approach, allowing for the selection of the most probable damage scenario and pinpointing its location and severity within the system. 0. Learning data-driven stable corrections of dynamical systems: Application to the simulation of the top-oil temperature evolution of a power transformer. Energies. iments. The full order solutions are obtained using any commercial software, and the PGD algorithm combined with one-dimensional approximation bases is then applied to build a parametric solution in which the separation of variables is applied to the parametric space.

Proper Orthogonal Decomposition

Proper Orthogonal Decomposition (POD) [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF][START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent ows[END_REF] is a method which consists in reducing the complexity of a problem by using a low-dimensional representation of its state instead of the initial high-dimensional representation. It is the application of Principal Component Analysis (PCA) [START_REF] Ringnér | What is principal component analysis?[END_REF] to Model Order Reduction.

Let us consider a general system of N ordinary dierential equations (ODE) which can in particular originate from the spatial discretization of a partial dierential equation (PDE):

E(p) du dt + A(p)u = f (p)
where u : (t; p) → u(t; p) ∈ R N is the state of the system being studied,

p ∈ D ⊂ R d is a vector of d parameters, E(p) ∈ R N ×N , A(p) ∈ R N ×N and f (p) ∈ R N .
Using any method of integration, we can compute a set of n s "snapshots" {u 1 , u 2 , ..., u n s } where u i = u(t i ; p i ) is the solution of the problem at t = t i and for the parameters p = p i ∈ D. The snapshots can also come from prior knowledge of the system or measurements.We can construct a matrix containing all the snapshots, the "snapshot matrix" U ∈ R N ×ns which has u i as its i-th column.

The Singular Value Decomposition (SVD) of U can be written U = LΣR T , L ∈ R N ×N and R ∈ R ns×ns being respectively the left and right singular vectors and Σ ∈ R N ×ns the matrix containing the singular values on its diagonal. The POD basis Φ = (ϕ 1 , ..., ϕ M ) is dened as the M ∈ N left singular vectors of U which correspond to the M largest singular values. This is an orthogonal basis which gives a good approximation of the snapshots provided M is chosen large enough. In practice, M is chosen by considering the singular values since they give a good idea of how much "energy" or "information" each basis function represents.

The POD reduced model is then written:

E M (p) du M dt + A M (p)u M = f M (p) where E M (p) = Φ T E(p)Φ ∈ R M ×M , A M (p) = Φ T A(p)Φ ∈ R M ×M and f M = Φ T f ∈ R M .
Computing the reduced solution u M is much faster than computing u because M << N . The approximation of the solution is then given by u ≈ Φu M .

One of the drawbacks encountered in POD is that it requires solving a large eigenvalue problem to obtain the reduced basis since the snapshot matrix has the size of the number of nodes or elements in the considered mesh, which can easily reach millions in industrial applications. This cost can be alleviated with the so-called snapshot POD [START_REF] Hilberg | The application of classical pod and snapshot pod in a turbulent shear layer with periodic structures[END_REF] by considering the fact that the number of snapshots is much smaller than the size of the discretization in space, hence solving for the eigenvalues of U T U is much faster than solving for those of U U T .

In a non-linear setting, this approach is not as ecient as in linear cases, because the computation of non-linear terms involves recalculating the full solution at each time step, increment or iteration of the non-linear solver, which almost completely negates the positive impact of POD on the computational complexity. Dierent methods have been proposed in the literature to circumvent this issue, such as the Empirical Interpolation Method (EIM) [START_REF] Chaturantabut | Nonlinear model reduction via discrete empirical interpolation[END_REF][START_REF] Peherstorfer | Localized discrete empirical interpolation method[END_REF] or hyper reduction [START_REF] Ryckelynck | Hyper-reduction of mechanical models involving internal variables[END_REF][START_REF] Farhat | Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear nite element dynamic models[END_REF].

Proper Generalized Decomposition

Proper Generalized Decomposition (PGD) [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF][START_REF] Ladevèze | The latin multiscale computational method and the proper generalized decomposition[END_REF] is a method which consists in computing a solution to a PDE in a reduced basis compared to an ordinary nite element basis. However, unlike POD, the reduced basis is not computed beforehand but at the same time as the solution. The numerical solution to a PDE using a nite element scheme can be written:

u(x, t) ≈ N m=1 α m (t)w m (x)
where the w i are the basis functions, each related to one node of the mesh and the α i are the unknown coecients. Independently of the problem at hand, this approximation involves a separation of variables between time and space.

Inspired by this idea, we can consider a more general space-time separation:

u(x, t) ≈ M m=1 X m (x)T m (t)
which is similar to the previous notation but this time both the X i and T i are unknown and will be computed by the solver. The number of terms in the sum M is also chosen by the solver based on a convergence criterion.

The solution is computed by iteratively solving a succession of one-dimensional problems. Before the M -th enrichment step, we have:

u M -1 (x, t) = M -1 m=1 X m (x)T m (t) The M -th term X M (x)T M (t) is computed such that u M (x, t) = u M -1 (x, t) + X M (x)T M (t)
X M (x) and T M (t) are calculated at the same time: starting from an arbitrary initial guess T 0 M (t), the alternating algorithm computes X j M (x) from T j-1 M (t), and then T j M (t) from X j M (x) by injecting everything into the PDE and solving the 1D problem obtained, until the condition

∥X j M (x)T j M (t) -X j-1 M (x)T j-1 M (t)∥ < ϵ is reached.
The same principle can be translated to separated representations in higher dimension, for instance separation of the dierent dimensions of the physical space, the phase 1.1 Intrusive and non-intrusive Model Order Reduction space or the parametric space:

u(x, y, z, t) ≈ M m=1 X m (x)Y m (y)Z m (z)T m (t) u(r, θ, t) ≈ M m=1 R m (r)Θ m (θ)T m (t) u(x, ẋ) ≈ M m=1 X m (x) Xm ( ẋ) u(x, t, p 1 , ..p d ) ≈ M m=1 X m (x)T m (t) d k=1 f k m (p k ) 1.1.

Proper Orthogonal Decomposition with interpolation

Proper Orthogonal Decomposition with Interpolation (PODI) [START_REF] Tezzele | Shape optimization through proper orthogonal decomposition with interpolation and dynamic mode decomposition enhanced by active subspaces[END_REF][START_REF] Rama | Proper orthogonal decomposition with interpolation-based real-time modelling of the heart[END_REF] is the nonintrusive counterpart of POD. The reduced basis is extracted through the same procedure as in intrusive POD, however the reduced coecients for a new set of parameters are obtained by using an interpolation or regression technique instead of projecting the PDE and associated solution in the lower-dimensional space.

From the matrix of snapshots U , the reduced basis Φ = (ϕ 1 , ..., ϕ r ) is calculated using SVD or the so-called snapshot-POD as described in 1.1.1. The reduced coordinates are similarly collected by extracting the rst M rows of matrix ΣR, where the mth column of this extracted matrix contains the reduced coordinates α i related to snapshot u i which allow reconstructing the full solution in an approximated manner:

u i ≈ Φα i .
Then, an interpolation or regression method is used to construct the function α(p) such that for all i, α i ≈ α(p i ). This can be performed with interpolation algorithms such as polynomial approximation, piecewise-polynomial interpolation [START_REF] Ly | Modeling and control of physical processes using proper orthogonal decomposition[END_REF], kriging [START_REF] Hamdaoui | Pod surrogates for real-time multi-parametric sheet metal forming problems[END_REF], Radial Basis Function (RBF) interpolation [START_REF] Nguyen | An ecient podi method for real-time simulation of indenter contact problems using rbf interpolation and contact domain decomposition[END_REF] ... or regression algorithms such as Moving Least Squares (MLS) [START_REF] Rama | Proper orthogonal decomposition with interpolation-based real-time modelling of the heart[END_REF], Dynamic Mode Decomposition (DMD) [START_REF] Tezzele | Shape optimization through proper orthogonal decomposition with interpolation and dynamic mode decomposition enhanced by active subspaces[END_REF],

Neural Networks (NN) [START_REF] Kovárnová | Shifted proper orthogonal decomposition and articial neural networks for time-continuous reduced order models of transport-dominated systems[END_REF], ...

Sparse Proper Generalized Decomposition

The sparse Proper Generalized Decomposition (sPGD) [START_REF] Ibáñez | A multidimensional datadriven sparse identication technique: the sparse proper generalized decomposition[END_REF] is a non-intrusive reduction method making use of the separation of variables to construct the solution of a parametric problem, which can be eciently obtained in real time once the training phase is complete. The problem consists in constructing the parametric solution

u(p 1 , . . . , p d ) : Ω ⊂ R d → R N , (1.1.1)
that depends on d parameters p k , k = 1, . . . , d, belonging to the parametric space Ω, in which a sparse sample of n s points and their corresponding solutions have been collected from high-delity simulations.

The sPGD expresses the function u from a low-rank separated representation

u(p 1 , . . . , p d ) ≈ ũM (p 1 , . . . , p d ) = M m=1 u m d k=1 ψ k m (p k ), (1.1.2)
constructed from greedy rank-one updates. In the previous expression ũM refers to the approximation, M the number of employed modes, u m ∈ R N are vectors which can be read as spatial modes and ψ k m are the one-dimensional functions concerning mode m and dimension k.

ψ k m (p k ) = D j=1 N k j,m (p k )a k j,m = (N k m ) T a k m , (1.1.3)
where D represents the number of degrees of freedom (nodes) of the chosen approximation and N k m is the vector collecting the shape functions.

In the context of usual regression the approximation ũM results from ũM = arg min

u * ∥u -u * ∥ 2 2 = arg min u * ns i=1 ∥u(p i ) -u * (p i )∥ 2 , (1.1.4)
where ũM takes the separated form of Eq.(1.1.2), n s is the number of sampling points to train the model and p i the vectors that contain the input data points of the training set. Notice that, to avoid overtting, the number of basis functions D must be D < n s .

The approximation coecients of each one-dimensional function are computed by employing a greedy algorithm, such that, once the approximation up to order M -1 is known, the M -th order term reads

ũM = M -1 m=1 u m d k=1 ψ k m (p k ) + u M d k=1 ψ k M (p k ).
(1.1.5)

The computed function is expected to approximate u not only in the training set but in any point p ∈ Ω. The main issue is how to ally rich approximations and scarce available data, while avoiding overtting. For that purpose a modal adaptivity strategy (MAS) was associated to the sPGD, however, it has been observed that the desired accuracy is not achieved before reaching overtting or the algorithm stops too early when using MAS in some cases. This last issue implies a parametric solution composed of low order approximation functions, thus not getting an as rich as desired function. Some papers describing the just referred techniques are [START_REF] Ibáñez | A multidimensional datadriven sparse identication technique: the sparse proper generalized decomposition[END_REF][START_REF] Borzacchiello | Non-intrusive sparse subspace learning for parametrized problems[END_REF].

In addition, in problems where just a few terms of the interpolation basis are present (that is, there are just some sparse non-zero elements in the interpolation basis to be determined), the strategy fails in recognizing the true model and therefore lacks accuracy.

This methodology is presented as a Model Order Reduction technique, but it can also be considered as a regression technique, and it will be used as such in the following chapters.

Physics-Informed Neural Networks

Physics-Informed Neural Networks (PINN) [13] have emerged as a promising framework that combines the power of neural networks with the fundamental principles of physics to address complex scientic and engineering problems. Traditional machine

       L(u)(x, t) = f (p(x, t)) ∀(x, t) ∈ Ω × [0, T ] u(x, t) = f d (x, t) ∀(x, t) ∈ Γ d × [0, T ] ∂ n u(x, t) = f n (x, t) ∀(x, t) ∈ Γ n × [0, T ] u(x, 0) = u 0 (x) ∀x ∈ Ω (1.2.1)
where L is a dierential operator, f is the loading or source term, p(x, t)R d is the vector which contains the properties of the loading or source at point x and time t,

Ω is the physical domain, f d (resp. f n ) is the Dirichlet (resp Neumann) boundary condition taking place on Γ d (resp. Γ n ) and u 0 is the initial condition.
When no measurement data is available and all the requirements coming from the eld of PDE analysis [START_REF] Evans | Partial dierential equations[END_REF] are met, the problem is usually well-posed in the sense that it has a unique solution that satises completely all the constraints.

In the PINN methodology, the solution is approximated by a neural network û(x, t) which is trained by minimizing the loss function L Φ :

L Φ = λ Ω L Ω + λ d L d + λ n L n + λ 0 L 0 (1.2.2)
where

L Ω = 1 N Ω N Ω i=1 L(û)(x i Ω , t i Ω ) -f (p(x i Ω , t i Ω )) 2 (1.2.3) L d = 1 N d N d i=1 û(x i d , t i d ) -f d (x i d , t i d ) 2 (1.2.4) L n = 1 N n Nn i=1 ∂ n û(x i n , t i n ) -f d (x i n , t i n ) 2 (1.2.5) L 0 = 1 N 0 N 0 i=1 û(x i 0 , 0) -u 0 (x i 0 ) 2 (1.2.6)
and λ Ω , λ d , λ n , λ 0 ∈ R * + are the weights associated to each term. Finding appropriate values for these weights is far from an easy task, but several works have addressed this issue, including [125127].

Neural Networks for hybrid solutions

PINNs are a very appealing tool when trying to combine physics and data. Indeed, a neural network can be built by adding one by one to its loss function all the dierent constraints it needs to fullll, whether they come from physical knowledge in the form of equations and boundary conditions, or from measurement data.

In section 1.2.1, only the physics-based model is considered, and the solution obtained with this methodology is the physics-based solution. Now, consider that measurements have been performed on the real system and data has been collected in various locations and at dierent times: u m (x i m , t i m ), i = 1, . . . , N m . To obtain the hybrid solution taking into account this additional information, the loss function becomes:

L H = λ Ω L Ω + λ d L d + λ n L n + λ 0 L 0 + λ m L m (1.2.7)
where

L m = 1 N m Nm i=1 û(x i m , t i m ) -u m (x i m , t i m ) 2 (1.2.8)
and λ m is its associated weight.

Minimizing this loss yields a solution which is a compromise between the physics-based solution and the measurements. However, while this solution may be very useful in many cases, there is no possibility to learn from this what was missing in the modeling which could explain the discrepancy between the physics-based model and reality.

2

Advances in non-intrusive Model Order Reduction

Model Order Reduction (MOR) techniques involve expressing the solution of a problem, typically a partial dierential equation (PDE), in a reduced basis with a strong physical or mathematical foundation. This basis is often obtained oine from the results of a high-delity solver, such as through the proper orthogonal decomposition (POD) or reduced basis method (RB) [2]. By using a reduced basis, the complexity of the solution scales with the size of this basis, which is typically much smaller than the multi-purpose approximation basis associated with the nite element method (FEM).

While a reduced basis may result in some loss of generality, it can greatly reduce computing time while maintaining accuracy as long as the problem solution falls within the space spanned by the reduced basis. However, accuracy may suer if the desired solution cannot be accurately approximated by the reduced basis. To improve generality and accuracy, the Proper Generalized Decomposition (PGD) [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial dierential equations encountered in kinetic theory modeling of complex uids[END_REF][START_REF] Chinesta | The proper generalized decomposition for advanced numerical simulations: a primer[END_REF] can construct the reduced basis and solve the problem simultaneously, although this approach can be more intrusive.

To address intrusiveness, non-intrusive procedures have been developed. These procedures construct the parametric solution of a parametric problem from a number of high-delity solutions performed oine. These solutions can be computationally expensive, as dierent choices of the model parameters are used to constitute the design of experiments (DoE).

There are various techniques to approximate solutions to parametric problems, including standard polynomial approximations on sparsely sampled domains. However, caution must be taken when using these methods. Orthogonal polynomial bases with Gauss-Lobatto points as the Design of Experiments (DoE) can produce very accurate approximations, but the number of samples required increases exponentially with the number of dimensions and polynomial degree. Using randomly sampled DoE or an overly complex approximation can result in overtting. To avoid this, one option is to use a basis that avoids these spurious oscillations, such as kriging approximations [START_REF] Papritz | Spatial prediction by linear kriging[END_REF],

or to restrict the polynomial approximation to a low degree.

Another approach is to use Proper Orthogonal Decomposition with Interpolation 28 2.1 Regularization of separated representations (PODI) [START_REF] Tezzele | Shape optimization through proper orthogonal decomposition with interpolation and dynamic mode decomposition enhanced by active subspaces[END_REF][START_REF] Rama | Proper orthogonal decomposition with interpolation-based real-time modelling of the heart[END_REF], where usual regressions are used to express the dependence of the modal coecients on the parameters. Sparse Subspace Learning (SSL) [START_REF] Borzacchiello | Non-intrusive sparse subspace learning for parametrized problems[END_REF] interpolates pre-computed solutions related to the DoE over the entire parametric space using a hierarchical approximation basis, which ensures the separated representation of the interpolated parametric solution. However, the volume of data required for SSL increases exponentially with the number of parameters involved, and the use of higher degree approximations with very little data increases the risk of overtting.

To address this issue, a sparse PGD (sPGD) [START_REF] Ibáñez | A multidimensional datadriven sparse identication technique: the sparse proper generalized decomposition[END_REF] was proposed, which uses a sparse sampling and adaptive approximation bases to limit over-oscillating behaviors. The sPGD is a nonlinear regression that uses the separation of variables, making it useful for multi-parametric settings. The Modal Adaptive Strategy (MAS) is used by the authors in [START_REF] Ibáñez | A multidimensional datadriven sparse identication technique: the sparse proper generalized decomposition[END_REF] to keep the degree of the approximation basis to a minimum for the rst PGD modes and increase it progressively for higher level modes. Other choices of the approximation bases, such as kriging, can also be used to limit overtting. This method is described in 1.1.4.

Besides, Model Order Reduction often focuses on the prediction of temporal, spatial or spatio-temporal responses, which could be a lot more complex than scalar quantities of interest and which are often unsuited to "point by point" interpolation or interpolation in a Euclidean space. A classical example of such an issue in the context of MOR is the Kolmogorov barrier [131133], a phenomenon related to the slow decay of the Kolmogorov n-width [START_REF] Kolmogoroff | Uber die beste annaherung von funktionen einer gegebenen funktionenklasse[END_REF][START_REF] Pinkus | N-widths in Approximation Theory[END_REF] in advection-dominated problems, which greatly reduces the eciency of classical MOR techniques.

Valuable strategies to circumvent these diculties often rely on curve alignment operations [START_REF] Wang | Alignment of curves by dynamic time warping[END_REF][START_REF] Bork | Chromatographic peak alignment using derivative dynamic time warping[END_REF] or registration methods [START_REF] Taddei | A registration method for model order reduction: data compression and geometry reduction[END_REF][START_REF] Blickhan | A registration method for reduced basis problems using linear optimal transport[END_REF] which act as a sort of pre-processing step on the data to map the solution to a reference domain where the computations are performed.

This chapter is divided into four sections which explore and attempt to give solutions to the aformentioned challenges. First, the use of separated representations in MOR is combined with regularization methodologies [START_REF] Bickel | Regularization in statistics[END_REF] to enhance the power of expression of the technique while reducing the potential overtting. Subsequently, the ANOVA-PGD is presented, a novel strategy which makes use of the so-called ANOVA decomposition [15] to enable the representation of complex parametric behaviors in MOR. Afterwards, the prediction of curves and temporal responses is tackled through the use of feature alignment and an adapted representation of the data. Lastly, challenges related to advection and localization of features in two or three dimensions are dealt with through the use of Optimal Transport [16].

Regularization of separated representations

Regressions are widely utilized in the eld of articial intelligence in general, and specically in supervised scientic machine learning [START_REF] Udrescu | Ai feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity[END_REF][START_REF] Brunton | Discovering governing equations from data by sparse identication of nonlinear dynamical systems[END_REF][START_REF] Hernández | Structure-preserving neural networks[END_REF] and the development of cognitive or hybrid digital twins [141143], as well as in neuroscience. Regression serves as a crucial component in automatically constructing models that represent the physical reality around us, which is essential for enabling articial intelligence to operate in the physical world [START_REF] Moya | Physically sound, self-learning digital twins for sloshing uids[END_REF][START_REF] Moya | Learning slosh dynamics by means of data[END_REF].

When implementing regression with limited data, the main challenge that arises is dealing with nonlinearities in multi-parametric settings. This leads to the curse of Advances in non-intrusive Model Order Reduction 29 dimensionality, where the number of degrees of freedom or sampling points required in the parametric space increases exponentially in order to achieve accurate results [START_REF] Laughlin | The theory of everything[END_REF].

When constructing models, it is always preferable to keep them as simple as possible.

In other words, it is better to have parsimonious models rather than complex ones.

This principle, known as Occam's razor [START_REF] Udrescu | Ai feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity[END_REF][START_REF] Brunton | Discovering governing equations from data by sparse identication of nonlinear dynamical systems[END_REF], suggests that simpler explanations should be favored when explaining any physical phenomenon among the available options. In the literature, this is achieved by promoting sparsity in regression [START_REF] Brunton | Discovering governing equations from data by sparse identication of nonlinear dynamical systems[END_REF][START_REF] Ibañez | Some applications of compressed sensing in computational mechanics: model order reduction, manifold learning, data-driven applications and nonlinear dimensionality reduction[END_REF]. To obtain parsimonious models that address sparsity, it is benecial to perform regression by combining L2 and L1 norm regularization. Numerical experiments are then provided to show the performance of these two methods compared to the regular sPGD algorithm.

rs-PGD

For the sake of simplicity of representation but without loss of generality, let us continue by assuming that the unknown sought function u is scalar-valued,

u(p 1 , . . . , p d ) : D ⊂ R d → R,
and that it is to be recovered from scarce data.

The goal is therefore to nd a function ũM which minimizes the distance to the function to be recovered:

ũM = arg min u * ns i=1 |u(p i ) -u * (p i )| 2 2 ,
and that takes the separated form

ũM (p 1 , . . . , p d ) = M m=1 d k=1 ψ k m (p k ) = M m=1 d k=1 (N k m ) T a k m , (2.1.1)
where n s is the number of sampling points employed to train the model (training set).

Here, the superscript M is employed to highlight the rank of the sought function. How to determine the precise value of M will be detailed hereafter.

In the PGD framework, an greedy enrichment algorithm combined with iterative scheme based on an alternating direction strategy is usually used to solve the resulting non-linear problem. At enrichment step m, this strategy initializes randomly all the a k m , and then proceeds to update them iteratively. At iteration l of the alternating direction strategy, a l m is updated by xing all the a k m for k ̸ = l and minimizing 2.1 Regularization of separated representations the residual. The iterations proceed until reaching a xed point according to a userspecied tolerance.

The following matrix equations dene the systems needed to be solved to perform these updates:

M l a l m = r, (2.1.2)
where:

r =    u(p 1 ) -ũm-1 (p 1 ) . . . u(p ns ) -ũm-1 (p ns )   , M l =         k̸ =l N k m (p k 1 ) T a k m N l m (p l 1 ) T . . . k̸ =l N k m (p k ns ) T a k m N l m (p l ns ) T        
In the context of sPGD, Eq.(2.1.2) is solved in the Ordinary Least Squares (OLS) setting:

a l m = (M T l M l ) -1 M l r, (2.1.3)
The rs-PGD is based on putting a penalty term when solving (2.1.2) in order to reduce overtting.

It is important to recognize that in the sPGD context, the problem of overtting can easily arise when employing high-order approximations, which are exacerbated by separated representations. This is primarily due to the typical usage of unstructured low data during the model training process. Overtting signicantly hampers the model's performance when faced with new and unseen datasets. Hence, the concept of utilizing a penalty term aims at enhancing the model's capability to perform well on new samples, even if it results in an increase in bias or error within the training set for a specic set of basis functions.

Dierent regularizations can be chosen depending on the properties of the problem such as the Tikhonov regularization [START_REF] Golub | Tikhonov regularization and total least squares[END_REF] or the Elastic Net regularization [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF][START_REF] De Mol | Elastic-net regularization in learning theory[END_REF].

In this work, the Elastic Net is employed, including in one of its extreme forms, the Ridge regularization (which is also a special case of Tikhonov regularization) and which is here presented in the rst place.

For this purpose, Eq.(2.1.3) is rewritten:

a l m = (M T l M l -λI) -1 M l r, (2.1.4)
where λ is the penalty factor and I is the identity matrix. In this case, all the dimensions are equally penalized but dierent penalty factors could be considered depending on the considered dimension.

The regularized problems associated to Eq.(2.1.4) is:

a l m = arg min a * ∥r -M l a * ∥ 2 2 + λ∥a * ∥ 2 2 , (2.1.5)
where the problem is divided in solving a ridge regression problem for each dimension when computing a l m during the alternating direction xed point strategy.

While a null intercept term is assumed in the deduction of equations (2.1.4) and (2.1.5) it could easily be includedand treated as in standard ridge procedures when solving the corresponding linear regularized regression problem for each dimension during the alternating direction strategy.

Since one is generally looking for the mode with best predictive abilities in each enrichment, the proposed criterion to choose λ is to perform a k-fold cross-validation [START_REF] Browne | Cross-validation methods[END_REF][START_REF] Fushiki | Estimation of prediction error by using k-fold cross-validation[END_REF] and select the value of λ that minimizes the cross-validated sum of squared residuals (or some other measure). It is also possible to use the one-standard error rule (heuristic) with cross-validation, in which the chosen model is the most penalized model which has an error that is no more than one standard error above the error of the best model. Such a rule acknowledges the fact that the tradeo curve is estimated with error, and hence takes a conservative approach [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF].

As the terminology used in this section shows, a regularization problem is formulated at each enrichment step. Thus, the best penalty factor is obtained at each updating stage, adapting the regularization whenever the approach is enriched. Other options have been investigated but numerical experiments suggest that this option is the most promising one.

If enough data is available, the split of the training set in two subgroups is equally a reasonable option to select λ and in addition, computationally less demanding. In this case, one subgroup is employed to construct the model and the other one to evaluate the predictive ability and then to select λ accordingly.

The Elastic Net regularization is obtained by adding an additional L1-norm regularization, from which Eq.2.1.5) becomes:

a l m = arg min a * ∥r -M l a * ∥ 2 2 + αλ 2 ∥a * ∥ 2 2 + (1 -α)λ 1 ∥a * ∥ 1 , (2.1.6)
where α ∈ [0, 1] and λ 2 and λ 1 are the penalty coecients that aect the L2 and L1 norm penalization terms respectively. Once again, these coecients could also be dierent for the dierent dimensions. The limit cases α = 0 and α = 1 result in the Ridge and Lasso regressions respectively.

s 2 -PGD

For the sake of simplicity of representation but without loss of generality, let us continue by assuming the same scalar unknown function discussed in section 2.1.1.

In this case, the sought solution admits a sparse solution for a certain basis using the PGD separated form. The goal is therefore to identify the correct non-zero coecients at each enrichment step in order to guide the approach to the correct separated representation.

Without a roadmap to select these nonzero coecients, the traditional sPGD fails to capture the true relationship between the model's features as well as its nal response.

Furthermore, if high-order terms appear in the searched function, this issues become even worse leading to serious overtting issues.

Let us consider the theory discussed in the previous section but now considering the L1 regularization with the aim of promoting sparsity in the coecients of the solution.

The L1 regularization is convenient because the nonlinear problem can be solved using the PGD constructor, with an alternating direction xed point strategy, where just a LASSO regression problem is considered in each dimension. Therefore, the regression problems for the iterative scheme will be:

a l m = arg min a * ∥r -M l a * ∥ 2 2 + λ∥a * ∥ 2 1 , (2.1.7)
that consists in solving a LASSO regression problem for each dimension when computing a l m within the alternate direction xed point strategy. Moreover, as previously discussed, in the present case again, all the dimensions are equally penalized but dierent penalty factors could be used.

As the LASSO problem are iteratively solved in each dimension, each one-dimensional function has a sparse representation with respect to its base choosing the right penalty factor. Again, a null intercept term is assumed in the formulation, but the penalty term could be non zero in a more general setting.

In the case of looking for sparsity only in some of the dimensions, Eq.(2.1.7) only applies in the computation of the considered dimensions, whereas the other dimensions are treated by invoking the standard sPGD or its ridge regularization counterpart, rs-PGD, addressed in the previous section.

To determine λ, the discussion of the previous section still applies. However, the following considerations can also be applied in the case of the s 2 -PGD:

Before selecting a model according to the predictive criterion, a lter is considered taking only the models with a minimum sparsity criterion ∥⃗ a l m ∥ 0 ≤ χ lim for the dimensions in which the sparsity is desired. Note: ∥ • ∥ 0 is dened by ∥x∥ 0 = #{k : x k ̸ = 0} and is refered to as "L0 norm", although it is actually not a norm.

Once model selection is performed, the OLS methodology is employed with the detected non-zero elements to obtain the correct update. The reason of this step is that LASSO regression terms are in general not accurate, and so it may be necessary to de-bias the obtained values. This is because the LASSO shrinkage causes the estimates of the non-zero coecients to be biased towards zero and in general they are not consistent [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF][START_REF] Brunton | Data-driven science and engineering: Machine learning, dynamical systems, and control[END_REF].

Numerical experiments

This section presents the results obtained from employing the aforementioned techniques in various scenarios. Initially, the reduction in error is demonstrated by comparing the rs-PGD approach with the classical procedure (sPGD). Then, the identication of sparsity and error reduction is showcased by contrasting the s 2 -PGD method with the standard sparse procedure (sPGD).

Results for the rs-PGD approach

The following examples consider the utilization of Elastic Net Regularization. To accomplish this, the parameter α that cotronls the relative importance between Ridge and Lasso regression is employed. α is determined by executing the algorithm multiple times with various α values, and subsequently selecting the value that yields superior predictive performance.

Advances in non-intrusive Model Order Reduction 33 In the rst example, the method is applied to try to approximate the ve-dimensional function

f (x 1 , x 2 , x 3 , x 4 , x 5 ) = (8x 3 1 -6x 1 -0.5x 2 ) 2 +(4x 3 3 -3x 3 -0.25x 4 ) 2 +0.1(2x 2 5 -1). (2.1.8)
The objective is to reconstruct the aforementioned function within the domain D = [-0.51, 0.51] 5 . The training set consists of 160 points, and thus, only these points are utilized to construct the model using either the sPGD or the rs-PGD methodology.

The Latin hypercube sampling (LHS) technique is employed to generate this set of data.

In contrast, a testing set comprising 54000 untrained points is considered to compare the results between the techniques when predicting unseen scenarios. This second set enables the evaluation of the predictive capability of both models once they are constructed.

Concerning the sPGD, a standard MAS is employed, reaching 4th degree polynomials.

To measure the error of both methodologies in the testing set, the following error criterion is used:

err pgd = ∥f -f pgd ∥ 2 ∥f ∥ 2 ; err rpgd = ∥f -f rpgd ∥ 2 ∥f ∥ 2
where f is the vector containing the values of f (x 1 , x 2 , x 3 , x 4 , x 5 ) in the testing set, f pgd and f rpgd are the vectors containing the prediction in the testing set of both methodologies (sPGD and rs-PGD, respectively).

After employing the discussed techniques in the above conditions, in this example the error is reduced by 52.38 % using the rs-PGD with α = 0.1.

To perceive the improvements and the overtting reduction, Concerning the sPGD, again a standard MAS is used up-to fourth-degree polynomials.

An error reduction of about 47% is accomplished with α = 0.5.

Results for the s 2 -PGD approach In the rst numerical example for s 2 -PGD, the considered function is: On the other hand, a testing set consisting of 27000 untrained points is available to compare the results between techniques when predicting unseen scenarios. This second set enables the evaluation of the predictive ability of both models once they are constructed.

f (x 1 , x 2 , x 3 ) = (sin(2x 1 ) -3.14)T 5 (x 2 ) + exp(x 3 ) cosh(x 1 ), ( 2
In this example, the conditions for employing the s 2 -PGD are as follows: a basis that encompasses eighth-degree polynomials is selected for the sparse dimensions.

Advances in non-intrusive Model Order Reduction 35 Additionally, a standard MAS-based s-PGD is employed, utilizing fourth-degree polynomials along the non-sparse dimensions. To detect sparsity, three simulations of the s 2 -PGD are conducted, penalizing a dierent dimension at each iteration. Consequently, the model with the best predictive ability (outside the training set) is selected.

As anticipated, the chosen model corresponds to the one obtained when penalizing the x 2 dimension.

Figure 2.4 presents the outcomes obtained from the standard sPGD methodology.

In this case, it is evident that the predictions are poor, as this methodology fails to identify the desired sparse solutions. This is one of the challenges faced by the sPGD approach, which s 2 -PGD aim to address.

Moreover, upon observing the sPGD solution, it is apparent that all possible elements are non-zero, indicating a failure in identifying sparsity. As a second test case, the approximation of the following function is tackled:

f (x 1 , x 2 , x 3 , x 4 , x 5 ) = T 5 (x 1 ) + 2T 1 (x 1 ) T 2 (x 2 ) + 2T 4 (x 2 )
sin(2x 3 ) -3.14 log(3x 4 + 1.5) cos(x 5 )

+ exp(x 4 ) cosh(x 3 ) sinh(x 5 ) Regarding the s 2 -PGD methodology, a basis consisting of sixth-degree polynomials is chosen for the sparse dimensions. Moreover, a standard MAS (Model-Adapted Sparse) approach is utilized, employing up to fourth-degree polynomials in the nonsparse dimensions.

Figure 2.6 illustrates the outcomes obtained from the standard sPGD methodology.

In this case, the predictions are unsatisfactory due to the incorrect identication of non-zero elements in the separated representation, leading to overtting issues.

To identify sparsity, ve separate simulations of the s 2 -PGD are conducted, penalizing one dimension at a time. Consequently, the model with the best predictive ability (outside the training set) is selected. As expected, the chosen model corresponds to the one obtained when penalizing the x 1 dimension. This choice is made because it correctly identies the non-zero terms for x 1 and x 2 when penalizing x 1 alone. Finally, the errors for the sPGD and s 2 -PGD solutions are reported as err pgd = 46.39% and err s2pgd = 2.4% respectively. This section is divided in three parts. To begin with, the ANOVA decomposition

[6466] and its low-cost approximation, the anchored ANOVA decomposition [6769], are introduced. Then, two versions of the ANOVA-PGD strategy are presented, the rst one for the prediction of scalar qunatity of interests and the second one for vector and full-eld outputs in the context of non-intrusive MOR. Finally, numerical experiments are performed to examine the performance of this method compared to other Machine Learning and non-intrusive MOR methods.

Advances in non-intrusive Model Order Reduction 39

ANOVA decomposition and anchored ANOVA

The ANOVA decomposition of a function u(p 1 , . . . , p d ) :

D = D 1 × • • • × D d ⊂ R d → R,
where the D k are intervals of R, is an orthogonal decomposition based on the analysis of variance [14,15], a statistical model designed for data analysis. Thus, the function u(p) can be written as a sum of orthogonal functions:

u(p) = u 0 + d k=1 u k (p k ) + d k 1 =1 d k 2 =k 1 +1 u k 1 ,k 2 (p k 1 , p k 2 ) + • • • + u 1,2,...d (p 1 , p 2 , . . . , p d ), (2.2.1)
verifying

E k (u k 1 ,...,k b (p k 1 , . . . , p k b )) = 0, (2.2.2)
where E k refers to the expectation with respect to any coordinate k in the set (k 1 , . . . , k b ), 1 ≤ b ≤ d. This property results in the orthogonality of functions involved in the previous decomposition.

To prove it, consider for example a simple 2D case with, p = (p 1 , p 2 ), u(p) = u(p 1 , p 2 ).

Thus, with E 1 (u 1 (p 1 )) = 0, E 1 (u 1,2 (p 1 , p 2 )) = 0 and E 2 (u 1,2 (p 1 , p 2 )) = 0, it results E 1,2 (u 1,2 (p 1 , p 2 )u 1 (p 1 )) = E 1 [E 2 (u 1,2 (p 1 , p 2 )) u 1 (p 1 )] = 0.
The number of functions involved in this decomposition (without considering the constant term) is 2 d -1, and they can be parametrized by the integer n, n = 1, . . . , 2 d -1. The dierent functions involved in the ANOVA decomposition are linked to the conditional expectations according to:

         E(u) = u 0 E(u|p k ) = u k (p k ) + u 0 E(u|p k , p l ) = u k,l (p k , p l ) + u k (p k ) + u l (p l ) + u 0 . . . (2.2.3)
hence the functions involved in the ANOVA decomposition can be expressed from the expectations in the following manner:

         u 0 = E(u) u k (p k ) = E(u|p k ) -u 0 u k,l (p k , p l ) = E(u|p k , p l ) -u k (p k ) -u l (p l ) -u 0 . . . (2.2.4)
where E(u|p k , p l ) refers to the integration with respect to all the variables except p k and p l .

Variance-based sensitivity analysis

The variance of u, Var(u), taking into account the orthogonality of the functions involved in the ANOVA decomposition, reads

Var(u) = 2 d -1 n=1 E u 2 n = 2 d -1 n=1
Var(u n ), (2.2.6)

These coecients allow to evaluate the relative importance of all the functions in the decomposition and give a sense of the sensitivity of the function to each of its parameters and to the interactions between any subset of the parameters.

By construction,

2 d -1 n=1 S n = 1.
(2.2.7)

Anchored ANOVA Computing all the terms in the ANOVA decomposition implies evaluating many integrals in a high-dimensional spaces, which becomes extremely expensive whenever the number of variables is higher than 2 or 3. To alleviate those costly computations, a common approximation is performed by introducing the so-called anchor point c ∈ D and evaluating all the expectations with the Dirac measure δ c . As a conquence, computing integrals is equivalent to particularizing the function in the anchor point, except for the coordinates which are not integrated upon. For instance: The point c can be any point of the domain. It is typically chosen as the center of the domain, but it is also very convenient to dene it as the nominal point of the problem when it exists.

     E(u) = u(c) E(u|p k ) = u(c

ANOVA-PGD

The principle of ANOVA-PGD is to combine the ANOVA decomposition with the PGD-based techniques described in section 2.1. First, the standard anchored ANOVA is used to evaluate the constant term and the one-dimensional functions depending on each parameter u k (p k ), k = 1, . . . , d, by using an adequate sampling, a sort of multidimensional cross centered at the anchor point c. Since these problems are one-dimensional, the sampling covers the space quite well, and interpolation methods that fail in high-dimensional problems are well-suited here. Hence, for each parameter, u k (p k ) can be approximated with piecewise polynomials, splines, or gaussian processes [START_REF] Williams | Gaussian processes for regression[END_REF][START_REF] Bernardo | Regression and classication using gaussian process priors[END_REF] for instance. Then, the residual is evaluated in sampling points which fall outside of the multidimensional cross, and the sPGD, rs-PGD, or s 2 -PGD is applied to that residual, which contains the dierent interactions between the parameters. In that case, an enhanced sparse-sampling can be considered, trying to sample as well as possible the area near the borders of the parametric domain, since the center is already well handled by the ANOVA terms.

Advances in non-intrusive Model Order Reduction (2.2.10)

The multidimensional cross is represented in Fig. 2.8. The design of experiments is then constructed by combining the anchor point, a sampling of each branch and a Latin Hypercube Sampling. The combination of the anchor point and the samplings from all the branches, which is a discretization of D c , is referred to as the ANOVA sampling. The n k points which sample the branch related to parameter p k are written p k,j , j = 1, . . . , n k and are "sorted" such that ∀(j 1 , j 2 ),

j 1 < j 2 =⇒ p k k,j 1 < p k k,j 2 .
An example of such a design of experiments is shown in Fig. 2.9.

ANOVA-PGD for scalar quantities of interest

The ANOVA-PGD strategy is rst detailed in its version aiming at predicting a scalar quantity of interest u which depends on d variables:

u(p 1 , . . . , p d ) : D ⊂ R d → R.
It should be noted that the parametric space is required to be a hypercube: D = × d k=1 D k . Data for u is rst collected according to the design of experiments strategy described above. The parameter values in the design of experiments are p i , i = 1, . . . , n s , among which are the adequate p k,j and the anchor point, p 0 = c.

The solution is searched in the form:

               ũ(p 1 , . . . , p d ) = u ANOVA (p 1 , . . . , p d ) + u PGD (p 1 , . . . , p d ) u ANOVA (p 1 , . . . , p d ) = u 0 + d k=1 u k (p k ) u PGD (p 1 , . . . , p d ) = M m=1 d k=1 ψ k m (p k ) (2.2.11)
The dierent functions are computed according to:

1. u 0 is obtained from u(p 0 ) = u(c).

2. Each function u k is constructed by setting u k (p k k,j ) = u(p k,j ) -u 0 for j = 1, . . . , n k and interpolating between the points using a spline or a gaussian process regression.

3. The residual is computed by removing the "ANOVA" term from all the available data: ∀i = 1, . . . , n s , r i = u(p i ) -u ANOVA (p i ). It can be noted that ∀i = 1, . . . , n s , p i ∈ D c =⇒ r i = 0.

4. The function u PGD is then obtained using a method described in 2.1 to predict the residual r.

Once these steps are completed, ũ can be evaluated for any p ∈ D to approximate u.

ANOVA-PGD for vector and full-field outputs

The function of interest is again a function dened in a hypercube of dimension d, the only dierence being that it is now vector-valued.

u(p 1 , . . . , p d ) : D ⊂ R d → R N .
The same design of experiments is built and data is collected accordingly. The solution now takes the form:

               ũ(p 1 , . . . , p d ) = u ANOVA (p 1 , . . . , p d ) + u PGD (p 1 , . . . , p d ) u ANOVA (p 1 , . . . , p d ) = u 0 + d k=1 u k (p k ) u PGD (p 1 , . . . , p d ) = M m=1 ϕ m d k=1 ψ k m (p k ) (2.2.12)
The workow to construct the dierent functions is similar to the one introduced above, but adapted to the vector form of the solution. Moreover, two additional changes take into account some considerations which appear when dealing with vector outputs: the problem may be more complex than for a scalar output because of the presence of features in the vectors, which calls for caution: a selection criterion is introduced to decide whether or not the dierent terms dened in the ANOVA decomposition are relevant to predict the output; the increased volume of data available allows to extract information which should improve the performance of the algorithm: the one-dimensional setting empowered by the design of experiments yields N one-dimensional functions of each parameter, from which an adequate basis can be extracted and then used for the PGD method.

The construction of the functions thus follows the following steps:

1. u 0 is obtained from u(p 0 ) = u(c).

2. One-dimensional functions ûk are constructed by setting ûk (p k k,j ) = u(p k,j )-u 0 for j = 1, . . . , n k and interpolating between the points using a spline or a gaussian process regression.

3. u k is constructed by checking whether it has a relevant role in the prediction of the output:

u k =      ûk if ns i=1 u(p i ) -ûk (p k i ) ≤ ns i=1 ∥u(p i )∥ 0 otherwise. (2.2.13)
4. The residual is computed by removing the "ANOVA" term from all the available data: ∀i = 1, . . . , n s , r i = u(p i ) -u ANOVA (p i ). It can be noted that ∀i = 1, . . . , n s , p i ∈ D c =⇒ r i = 0. m refers to the fact that, when using the Modal Adaptivity Strategy, some of the basis functions are set to 0 in early iterations.

6. The function u PGD is then obtained using a method described in 2.1 to predict the residual r using the one-dimensional function bases N k m obtained in the previous step.

Using a basis extracted from the data in the PGD-based regression improves significantly the performance of the method when the searched functions can not be well described with a polynomial representation. A very usual example of this is when one of the parameters is a physical angle in the system; its associated responses are usually periodic functions that are very poorly decscribed with low-dimensional polynomials, especially if they are high-frequency functions. This behavior can be successfully extracted using the strategy detailed above, under the condition that enough points are used in the discretization of D k c . Therefore if a parameter p k is expected to yield rich responses, its associated branch D k c should be sampled with a ner discretization than for the other parameters. For instance geometry parameters tend to require richer functions than material properties such as conductivity, Young's modulus, ...

ANOVA-PGD: A novel strategy 2.2.3 Numerical experiments

This section presents a few examples which aim at testing the performance of a few variants of the ANOVA-PGD algorithm, and comparing them with other regression methods.

Simple example

For the rst example, a very simple case aiming at showcasing the power of ANOVAbased regression, a 2D function is considered: f (x, y) = -2 cos(3x 1.75 ) + 10 log(y -0.6) 4 + 6 cos(x)(y -0.3y 2 ), (2.2.14) that perfectly ts the ANOVA structure, because it contains complex terms which depend on each coordinate independently, 2 cos(3x 1.75 ) and 10 log(y -0.6) 4 respectively, and then a term which is less complex coupling both coordinates, 6 cos(x)(y -0.3y 2 ).

When considering the ANOVA-based sampling consisting of the center point of the parametric domain acting as the anchor c = (x c , y c ), 10 additional points on the rst branch (of the form (x, y c )) and 10 additional points on the second branch (of the form (x c , y)), functions f x (x) and f y (y) are constructed using a cubic spline interpolation. Then, a standard 2D nonlinear regression using basis functions of the form (xx c ) m (y -y c ) n , m, n ≥ 1 (due to the low dimensionality of the treated problem the employ of separated representations is not needed) is employed to calculate the coupled term f x,y (x, y) using 4 extra sample points.

The constructed solution is depicted in Fig. 2.10 where it is compared with the exact solution as well as with the solution obtained by using the standard sPGD (with a Latin Hypercube Sampling containing 25 points), while Figs. 2.12 and 2.11 compare the predictions and the reference values. These results stress an excellent performance of the ANOVA-based regression. Among other things, it can be noted that the ANOVA is very helpful to capture the fast changing behavior near lower values of y due to the presence of the logarithm.

f (x) = d k=1 cos(ω c k x k + ϕ c k ) + sin(ω s k x k + ϕ s k ) + d-1 k=1 cos(ω c k x k + ϕ c k ) sin(ω s k+1 x k+1 + ϕ s k+1 ) 2 (2.2.15) where ω c k , ω s k ∈ [0, 3π] and ϕ c k , ϕ s k ∈ [0, 2π] are chosen randomly.
The design of experiments is constructed by setting the center point as the anchor point and by selecting 6 sample points on each branch (not including the anchor point).

A Latin Hypercube Sampling containing 399 points is used for the PGD term, which is solved with the standard sPGD. The total number of samples is therefore 1000. For comparison, the regression is also performed using a SVR with RBF kernel

[156158] and a fully-connected neural network with 4 layers of 50 neurons and ELU activation function [START_REF] Clevert | Fast and accurate deep network learning by exponential linear units (elus)[END_REF].

The results obtained using these three methods are shown in Figs. This problem can be considered to be extremely challenging because of the extreme dimension of the parametric space combined with nonlinear functions, which is why none of these methods produce results which are close to perfect. However, by exploiting the idea of dealing with each parameter one by one, the ANOVA strategy is able to extract valuable information from the data and clearly outperforms the other algorithms.

Cylindrical indentation testing of a lithium-ion cell

This example adresses the construction of a parametric surrogate model for the 3D displacement eld of a lithium-ion cell in a cylindrical indentation test. The cell and impactor are illustrated in Fig. 2.16. The results are shown in Fig. 2.17.

The model is able to predict quite well the solutions. However due to the localization of the impact related to the orientation of the impactor, a bit of noise is present in Advances in non-intrusive Model Order Reduction 49 the model. The next sections address this type of issue.

Figure 2.17: Comparison between the solutions (displacement norm in mm) obtained with nite elements using LS-DYNA (left) and the surrogate model (right)

Reduced Modelling of curves

In the context of parametric surrogates, several nontrivial issues arise when a whole curve needs be predicted from given input features. For instance, dierent sampling or ending points lead to non-aligned curves. This also happens when the curves exhibit a common pattern characterized by critical points at shifted locations (e.g., in mechanics, the elastic-plastic transition or the rupture point for a material). In such cases, classical interpolation methods fail to give physics-consistent results and appropriate pre-processing steps are required. Moreover, when bifurcations occur into the parametric space, to enhance the accuracy of the surrogate, a coupling with clustering and 50 2.3 Reduced Modelling of curves classication algorithms is needed. In this section, several methodologies are proposed to overcome these issues. The surrogates thus created are then exploited to quantify and propagate uncertainty, furnishing parametric stastistical bounds for the predicted curves. The procedures are exemplied over problems in Computational Mechanics.

Data alignment and uncertainty propagation

This section presents the curve parameterization based on data alignment to obtain an accurate physics-informed interpolation. The procedure is applied to an example to study the mechanical response of parametric materials loaded in tension.

The problem considered is a parametric study over dog bone tensile test samples, as sketched in Fig. 2.18. It aims at studying the inuence of the 3 parameters (n, K, ε 0 ) characterizing the Krupkowski hardening law (also known as Swift hardening law ), widely used in FEM software

σ = K(ε + ε 0 ) n ,
linking True Strength and True Strain. ε denotes the eective plastic strain, ε 0 the oset strain, n the strain hardening exponent and K the material constant. Force-Displacement curves. However, for the sake of generality, such curves are referred to as generic functions g(x), presenting two characteristic behaviors in the so-called primary and secondary zones. In the specic case of Force-Displacement, the primary zone is the elastic response of the material, up to the yield point x E . The secondary zone is the post yield behaviour up to the failure point x F , as illustrated in Fig. 2.20. x E is referred to as the transition point" and x F as the end point", related to the specimen fracture. It is assumed that the behaviors in the primary and secondary zone, g 1 (x) and g 2 (x) respectively, and the transition and end points, x E and x F respectively, depend on a series of parameters grouped in vector p, i.e. g 1 (x; p) ≡ g(x ∈ [0, x E ]; p), g 2 (x; p) ≡ g(x ∈ [x E , x F ]; p), x E (p) and x F (p). Indeed, when considering dierent choices of the model parameter p i = (K i , n i , ε 0,i ), i = 1, . . . , n s , one obtains a set of curves, as the ones shown in Fig. 2.21, for instance. Such curves correspond to a sparse DoE (Latin Hypercube) of 20 points in the 3-dimensional parametric space D = I K × I n × I ε 0 , considering the parameters bounds specied in Table 2.2. Numerical simulations have been carried out with VPS simulation software from ESI Group. The variable x corresponds to the displacement in mm, while the function g(x) to the force in kN.

p i = (K i , n i , ε 0,i ), i = 1, . . . , n s .
Once the transition and end points of each curve have been determined, the curves can be rediscretized over the same number of points (through a standard piecewise linear interpolation, for instance). To align them, a dimensionless coordinate is dened in each zone, y in the primary zone, x ∈ [0, x E ], and z in the secondary zone, x ∈ [x E , x F ],

dened through the change of variable

y = x x E , y ∈ [0, 1] and x ∈ [0, x E ], (2.3.1)
and

z = x -x E x F -x E , z ∈ [0, 1] and x ∈ [x E , x F ], (2.3.2)
expressions that hold for each curve g(x; p i ), i = 1, . . . , n s , with

y = x x i E , y ∈ [0, 1] and x ∈ [0, x i E ],
(2.3.3)

and z = x -x i E x i F -x i E , z ∈ [0, 1] and x ∈ [x i E , x i F ].
(2.3.4) Fig. 2.22 depicts functions g 1 i (y) ≡ g 1 (y; p i ) and g 2 i (z) ≡ g 2 (z; p i ). Once the curves have been aligned, the nonlinear regressor presented in section 1.1.4 can be invoked to build the parametric metamodel of the curve. This can be done separately in each zone or over the whole newly dened coordinate x. However, before proceeding with the regression, an ulterior parametrization via the Proper Orthogonal Decomposition is addressed to achieve a further Model Reduction.

POD modes extraction

In order to extract the most signicant modes able to describe these functions, the POD can be applied in each group of curves in Fig. 2.22. This amounts to build the snapshot matrix within each group and perform a truncated SVD. In the case that serves here to illustrate the procedure, a single mode suces to describe the almost linear functions in the primary zone, that will be noted by ξ 1 (y), whereas in the secondary zone two functions are needed, ϕ 1 (z) and ϕ 2 (z).

Thus, any function g 1 i (y) can be expressed ∀i as

g 1 i (y) = α i 1 ξ 1 (y), (2.3.5)
whereas functions g 2 i (z), ∀i, read

g 2 i (z) = β i 1 ϕ 1 (z) + β i 2 ϕ 2 (z). (2.3.6)
The α and β coecients can be easily computed by simple projection, i.e.

1 0 g 1 i (y)ξ 1 (y) dy = α i 1 , (2.3.7)
where the normality of ξ 1 (y) was used. In the same way, and taking into account the orthonormality of functions ϕ 1 (z) and ϕ 2 (z),

1 0 g 2 i (z)ϕ 1 (z) dz = β i 1 , (
g 2 i (z)ϕ 2 (z) dz = β i 2 .
(2.3.9)

Thus, for each curve g i (x) can be extracted its ve main descriptors: x i E , x i F , α i 1 , β i 1 and β i 2 , all of them related to the features grouped in vector p i . Now, each of these descriptors can be expressed parametrically, x E (p), x F (p), α 1 (p), β 1 (p) and β 2 (p), by using any regression technique. In this work, the regression used is the sPGD as described in 1.1.4 in the special case of scalar output.

Curves reconstruction

When considering a choice of the parameters p, the curves descriptors are extracted from the regressions x E (p), x F (p), α 1 (p), β 1 (p) and β 2 (p), the dimensionless coordinates dening both zones are calculated from

y = x x E (p) → x = y x E (p), (2.3.10)
and

z = x -x E (p) x F (p) -x E (p) → x = x E (p) + z (x F (p) -x E (p)), (2.3.11)
and, nally, the curve in each zone reconstructed according to g 1 (y; p) = α 1 (p)ξ 1 (y),

(2.3.12) and

g 2 (z; p) = β 1 (p)ϕ 1 (z) + β 2 (p)ϕ 2 (z), (2.3.13)
from which the curve g(x; p) can be straightforward obtained via

g(x; p) =        α 1 (p)ξ 1 x x E (p) , x ∈ [0, x E (p)] β 1 (p)ϕ 1 x -x E (p) x F (p) -x E (p) + β 2 (p)ϕ 2 x -x E (p) x F (p) -x E (p) , x ∈ [x E (p), x F (p)].
(

2.3.14)

To build the parametric metamodel, 17 curves have been used to train the sPGD regressor, while the remaining 3 for testing. Fig. 2.24 shows the resulting predictions over 3 training points and test points.

Advances in non-intrusive Model Order Reduction

Real-time calibration

Now, given an experimental curve g(x), its parameters are extracted according to

x E from the point at which the change of behavior occurs (for instance, computing the function derivatives by means of nite dierences);

x F is the terminal point;

α 1 follows from y = x x E and 1 0 g(y)ξ 1 (y) dy = α 1 ;

β 1 follows from z = x-x E x F -x E and 1 0 g(z)ϕ 1 (z) dz = β 1 ; β 2 follows from z = x-x E x F -x E and 1 0 g(z)ϕ 2 (z) dz = β 2 .
Then, from the regression models x E (p), x F (p), x 1 (p), β 1 (p) and β 2 (p), the inverse problem is solved to extract the associated parameters, p.

Statistical model derived by parametric curves

With the previously built surrogate model, the curve related to any possible value of p can be computed in real-time, i.e. g(x; p). In this section, this surrogate will be employed for uncertainty quantication.

Each feature p k in vector p is assumed characterized by a Gaussian distribution dened by its mean value µ k and its variance σ 2 k , that is p k ∼ N (µ k , σ 2 k ). Assuming all p k being independent, we get

p ∼ N (µ, Σ), µ = (µ k ) d k=1 , Σ = diag(σ), σ = (σ 2 k ) d k=1 ,
where diag(•) is the diagonal matrix of diagonal •.

The aim is to link the sensitivity over the input features with the one over the output curve. This means computing some estimators of the average M and the variance Σ of the curve descriptors for dierent choices of µ and σ, and from them, by using the sparse PGD presented in Subsection 1.1.4, build the set of statistical surrogates: where O(p) denotes any quantity of interest (QoI) involved in the curves parametrization (i.e., an output depending on the input parameters; e.g., x E , x F , α 1 , β 1 and β 2 in the example presented before) and M and Σ the corresponding estimators for mean and variance, respectively. This allows calculating the envelopes, for a given condence, of the curves, as sketched in Fig. 2.25.

To build the surrogate (2.3.15), for instance for the curve descriptor O(p), a training dataset of N s points shall be generated:

{(µ j , σ j ), ( MO(p j ) , ΣO(p j ) )} Ns j=1 .

This can be achieved by means of a Monte Carlo sampling, which gives the estimators of mean and variance for the curves g(x; p j (µ j , σ j )), and of any descriptor O(p j ), for j = 1, . . . , N s .

The whole procedure is summarized in Algorithm 1.

Data alignment and data clustering

This section focuses on the study of crack propagation in notched specimens loaded in tension, whose geometry is sketched in Fig. 2.27. The test piece has a V-shaped notch defect which is always in the same location (bottom-middle). On the other side of the test piece there is a half-circle groove. The goal is to predict the crack propagation from the defect in dierent congurations (dierent location S and radius R of the groove and test piece thickness h). Depending on the location of the groove, the crack will propagate dierently from the defect, sometimes towards the groove and in other conditions straight towards the other side of the specimen. Randomly sample (e.g., LHS) the model features means and variances

(µ j , σ j ), µ j = (µ j,k ) d k=1 , σ j = (σ 2 j,k ) d k=1 .
3:

Perform a Monte Carlo sampling of the curves statistical descriptor O(p):

4:
1. generate a population of N MC vectors of features p j = (p k j ) d k=1 , by sampling

N MC points from p k j ∼ N (µ j,k , σ 2 j,k ), k = 1, . . . , d; 5:
2. generate the population of the corresponding N MC curves (and of any QoI involved in their parametrization), by using the curves surrogate f X (p), that is, g(x; p j,l ), O(p j,l ) = f X (p j,l ), l = 1, . . . , N MC ; 6:

compute the population mean and variance to obtain the corresponding

Monte Carlo estimators for the curve g(x; p j ) and its descriptor O(p j ):

( Mg(x;p j ) , Σg(x;p j ) ), ( MO(p j ) , ΣO(p j ) ). Same procedure holds for the whole curve g(x; p). The focus is placed on the prediction of the Force-Displacement curves plotted in Fig. 2.28, which are considered as the generic functions g(x), following the same notation as in section 2.3.1.

It can be observed that all the curves present a similar pattern in the rst zone, monotonically increasing, while the response appears much dierent in the secondary zone. A rst pre-processing step consists in splitting the zones as illustrated in Fig.

2.29

, where x M denotes the point where the curve reaches its maximum value, while x F is its endpoint.

Cutting the curves yields the two groups of functions plotted in Fig. 2.30, which are of course not aligned. However, they can be expressed as functions of normalized coordinates y and z, respectively, and aligned following the dilatation procedure discussed in section 2.3.1. Figure 2.30: Functions g 1 i (x) ≡ g 1 (x; p i ) (left) and g 2 i (x) ≡ g 2 (x; p i ) (right), with p i = (R i , S i , h i ), for i = 1, . . . , n s .

Advances in non-intrusive

Once the alignment has been performed, using the usual nonlinear regression techniques and same notations as in section 2.3.1, two regression models, one for each 60 2.3 Reduced Modelling of curves group, can be established:

g 1 (x; p) := g(x ∈ [0, x M (p)]) = f X 1 (p) g 2 (x; p) := g(x ∈ (x M (p), x F (p)]) = f X 2 (p).
(2.3.16)

In Eq. (2.3.16), for the sake of clarity, x M and x F are specied since these points are involved into the parametrization of the functions g 1 (x) and g 2 (x), respectively, and thus expressed parametrically.

As previously pointed out, the second group of functions g 2 i (x), for i = 1, . . . , n s , presents really dierent shapes depending on the features p i . When bifurcations occur in the parametric space, the system responses related to two choices of the model parameters can be completely dierent. In such cases, a standard nonlinear regression over the full space can lead to inaccurate and nonphysical solutions. To enhance the accuracy of the model f X 2 (p), a more valuable route consists in exploring the parametric space prior to interpolation. This can be done via a clustering of the system responses. Once the clusters have been established, several regression submodels can be built, minimizing the risk of mixing spurious eects coming from other clusters.

Clustering

To exemplify the bifurcation problem in the parametric space, two dierent congurations of the model parameters are considered, resulting into the specimens shown in

     g 1 (x; p) = f X 1 (p) g 2 (x; p) = f X 2,1 (p) for C 1 f X 2,2 (p) for C 2 .
(2.3.17) Fig. 2.34 shows the functions in the secondary zone after the clustering.

Reduced Modelling of curves

Curves reconstruction and classification

For a newly dened choice of model features p * , the curve g(x; p * ) is obtained via

g(x; p * ) = g 1 (x; p * ), 0 ≤ x ≤ x M (p * ) g 2 (x; p * ), x M (p * ) < x ≤ x F (p * ), (2.3.18)
where g 1 and g 2 are obtained through Eq. (2.3.17).

The training of the regression models has been performed using 40 points of the DoE, while the remaining 10 have been used for testing. Moreover, a Support Vector Machine classier [START_REF] Mammone | Support vector machines[END_REF] (a Random Forest classier [START_REF] Biau | A random forest guided tour[END_REF] could also be used, for instance) has been trained to select the best regression submodel to predict g 2 (x; p * ). This classier has shown perfect accuracy, as shown by the Confusion Matrices in Fig.

Reduced Modelling based on Optimal Transport

The type of problems faced in section 2.3 can be very challenging when trying to make reduced models of one-dimensional entities such as curves, but become extremely troublesome when the model concerns full-eld solutions in two, three or more dimensions.

One of the usual diculty which hinders the application of Reduced Order Modelling is the presence of features in the solution which move with time or as a result of a change in the parameters. The most popular and ecient ROM algorithms employ linear reduced bases, which are very powerful in many cases, but are not able to capture a continuous displacement of solution features resulting of phenomena such as advection or localized parametric constraints (e.g. loadings, impacts, sources, cracks, localized damage, . . . ). The adequacy of a reduced basis with a problem is studied within the theory of Kolmogorov n-width [START_REF] Kolmogoroff | Uber die beste annaherung von funktionen einer gegebenen funktionenklasse[END_REF][START_REF] Pinkus | N-widths in Approximation Theory[END_REF], which quanties the ability of a basis to represent faithfully the solutions depending on its dimension. Moving features and advection result in a slow evolution of the n-width with respect to the dimension of the reduced basis, a phenomenon known in MOR as the Kolmogorov barrier [131133].

The ineeciency of reduced bases in the aforementioned situations is related to their use of a Euclidean metric. This metric is not adequate for problems regarding misalignement or feature localization, because it only considers point by point dierence which never allows to detect information such as the distance between noteworthy features. However, there exists a collection of metrics which are designed for this precise objective, the ones derived from the theory of Optimal Transport [16,17].

Optimal Transport is witnessing a renewed interest because of its benets for data science and machine learning which can rely on a strong mathematical background [START_REF] Villani | Topics in optimal transportation[END_REF],

wide array of formulations [164166], as well as solid foundations for computational strategies [167169] and ecient implementations [START_REF] Flamary | Pot: Python optimal transport[END_REF].

Advances in non-intrusive Model Order Reduction 65

The theory of Optimal Transport introduces the notion of spaces of probability measures, the Wasserstein spaces [START_REF] Panaretos | An invitation to statistics in Wasserstein space[END_REF]. Although these spaces are often too restrictive to work with solutions of PDEs, which live in much larger spaces, they oer an interpolation framework [START_REF] Kell | On interpolation and curvature via wasserstein geodesics[END_REF][START_REF] Chewi | Fast and smooth interpolation on wasserstein space[END_REF] that is an attractive alternative to standard interpolation in Euclidean spaces. Besides the interpolation between two objects, the concept of barycenter in the Wasserstein space enables interpolation from multiple objects simultaneously, which could be a valuable route for regression and model reduction.

However, the high computational cost and inadequacy of barycenters in high dimension limit their applicability in many cases.

To apply Optimal Transport theory to regression, a very appealing strategy consists in using Linearized Optimal Transport [START_REF] Mérigot | Quantitative stability of optimal transport maps and linearization of the 2-wasserstein space[END_REF][START_REF] Sarrazin | Linearized optimal transport on manifolds[END_REF] to translate objects from the Wasserstein space to the tangent space, and perform the usual operations of machine learning in this Euclidian space, before mapping the results back to the Wassertsein space. This workow introduces two problems, the challenges and cost related to computing the Exponential and Logarithmic maps [START_REF] Sarrazin | Linearized optimal transport on manifolds[END_REF], and the risks of operating away from the initial point [START_REF] Petersen | Riemannian geometry[END_REF].

This section presents a methodology combining the framework of ANOVA described in section 2.2 and Optimal Transport to create parametric reduced models. The theory of Optimal Transport is rst briey introduced, then the ANOVA-empowered OT regression strategy is described before it is applied to a numerical example.

Short introduction to Optimal Transport

This section introduces the basics of Optimal Transport and a few technical aspects relevant to describe the methodology which will be described later. For further insight and extensive information on the theory and applications of Optimal Transport, the interested reader can refer to [START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Peyré | Computational optimal transport[END_REF][START_REF] Peyré | Computational optimal transport: With applications to data science[END_REF].

Measure couplings and transport maps

Considering two measures µ and ν on two metric spaces X and Y, coupling of µ and ν means to construct a measure π on X × Y such that π admits µ and ν as marginals on X and Y respectively.

The most trivial coupling is the measure µ ⊗ ν which is the probability law of (X, Y ), where X and Y are two independent random variables such that µ is the law of X and ν is the law of Y . The opposite extreme case is when there exists a measurable function T : X → Y such that Y = T (X). In that case, π is called a deterministic coupling of µ and ν, and T is called a transport map. Deterministic couplings between µ and ν do not always exist, especially in discrete problems.

The objective of Optimal Transport is to nd a coupling between µ and ν that is optimal in the sense that it minimizes a functional J µ,ν : Π µ,ν → R + (where Π µ,ν is the set of couplings between µ and ν) which is dened as the cost of transportation of a coupling π between µ and ν.

In many cases, the optimal coupling is deterministic, but this not always the case which is not surprising since deterministic couplings do not always exist. In the methodology proposed in section 2.4.2, when the optimal coupling is not deterministic, the transport map will be replaced by an approximation thanks to the notion of 66 2.4 Reduced Modelling based on Optimal Transport barycentric projection: T = E(π|x). This can also be expressed in the following way:

∀z ∈ X , T (z) = X ×Y yδ z (x)dπ(x, y) (2.4.1)

Kantorovitch problem

Introducing a cost function c : X × Y → R + , the formulation introduced by Kantorovitch to dene the optimal coupling is the following:

π * = arg min π∈Πµ,ν X ×Y c(x, y)dπ(x, y) (2.4.2)
In the problems tackled in section 2.4.2, X and Y are both subsets of R 2 or R 3 and c is the Euclidean distance such that c(x, y) = ∥x -y∥ 2 .

Furthermore, this formulation denes a metric, the p-Wasserstein distance W p [START_REF] Rüschendorf | The wasserstein distance and approximation theorems[END_REF], which in turn allows to construct the p-Wassertein spaces of measures when µ and ν operate on the same space with the right properties:

W p (µ, ν) = inf π∈Πµ,ν X ×Y d (x, y) p dπ(x, y) 1 p (2.4.3)
In order to solve Eq. 2.4.2, µ and ν and π are discretized and represented in the following manner:

µ = m i=1 α i δ x i (2.4.4) ν = n j=1 β j δ y j (2.4.5) π = m i=1 n j=1 γ ij δ (x i ,y j ) (2.4.6)
Now, by dening C = (c(x i , y j )) i,j and Π = (γ ij ) i,j , the optimal coupling π * can be obtained by solving the linear program:

min Π ⟨Π, C⟩ s.t. Π1 = α Π T 1 = β Π ≥ 0 (2.4.7)
where ⟨Π, C⟩ = m i=1 n j=1 γ ij c(x i , y j ).

An approximation of the optimal coupling can be obtained faster with the use of Entropic Optimal Transport [START_REF] Nutz | Introduction to entropic optimal transport[END_REF] but it was not necessary in this work.

Advances in non-intrusive Model Order Reduction 67

Linearized Optimal Transport

In order to combine the theory of Optimal Transport with linear methods and interpolation or regression in vector spaces, the Linearized Optimal Transport (LOT) [START_REF] Wang | A linear optimal transportation framework for quantifying and visualizing variations in sets of images[END_REF][START_REF] Lott | Some geometric calculations on wasserstein space[END_REF] is a very valuable tool.

Thanks to the theory of dierential geometry, it is possible to locally linearize the 2-Wasserstein space and thus to work on the tangent space at a certain point µ 0 . Any measure can be embedded into the tangent space at µ 0 using the logarithmic map Log µ 0 . Computing it requires nding the optimal coupling between µ 0 and µ and its associated transport map T . The logarithmic map is then computed as:

Log µ 0 (µ) = T -Id (2.4.8)
The inverse operation, the exponential map Exp µ 0 computes a measure from an em- bedding v of the tangent space in the following manner:

Exp µ 0 (v) = (Id + v) ♯ µ 0 (2.4.9)
where T ♯ µ = µ • T -1 is the pushforward measure dened such that, considering a suitable σ-algebra Σ:

∀A ∈ Σ, T ♯ µ(A) = µ(T -1 (A)). (2.4.10)
In practice, it is convenient to represent a discretized measure as a weighted point cloud, upon which the pushforward operation results in moving each point from x to T (x). It can be noted that the continuous interpolation between two measures dened by operating in the tangent space consists in moving each point along the trajectory

(1 -t)x + tT (x) for t ∈ [0, 1].

ANOVA-based Optimal Transport methodology

This section presents a strategy which combines the ANOVA-PGD described in section 2.2 with the aforementioned tools to build non-intrusive parametric reduced models.

To this aim, in what follows, the PDE solutions u will be transformed into measures by applying a few transformations. At rst, the methodology will be described for positive solutions with localized support, and then it will be adapted to a more general setting.

The motivation for this work is that the natural interpolation in the tangent space allows to move localized features of the solutions in the parametric domain. This is a very interesting ability in contrast with the interpolations performed in the vector space in which the solutions originally live, which do not move features around the space and are only able to increase or decrease their magnitude. The dierence between the two approaches is illustrated in Fig. 2.39.

2.4 Reduced Modelling based on Optimal Transport 3. The matrix Π 1→2 coupling µ 1 and µ 2 is obtained by solving the linear program 2.4.7 with α = µ 0 , β = µ 2 and matrix C computed from T 1 (X 0 ) and X 2 .

4. The transport map is computed by barycentric projection: ∀j = 1, . . . , M 0 , T 2 (X j 0 ) = T 1→2 (T 0→1 (X j 0 )) = (Π 1→2 X 2 ) j /µ j 0 . Measure µ 2 is now represented by (T 2 (X 0 ), µ 0 ).

5.

Steps 3 and 4 are repeated until reaching the end of the branch.

Remark: in order to have a more general reference representation, it is interesting to discretize µ 0 with a decomposition similar to the ones performed in Smoothed Particle Hydrodynamics (SPH) [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF]. Instead of computing µ 0 from Eq. 2.4.11, this consists in choosing a value of M 0 , setting µ j 0 = 1 M 0 for j = 1, . . . , M 0 and then constructing matrix X 0 by solving:

X 0 = arg min X Ω 1 M 0 M 0 j=1 K(∥X j -x∥) - u(x; p 0 ) Ω u(x; p 0 )dx dx (2.4.13)
where K is a suitable kernel function. △

T i (X 0 ) is stored in matrix T i ∈ R M 0 ×2
for all the i related to the ANOVA sampling (note that T 0 = Id). Then, the function T ANOVA : D → R M 0 ×2 (similar to u ANOVA from Eq. 2.2.12) is constructed.

Note: In section 2.2.2, the methodology is dened for vector-valued outputs. For a matrix, the same methodology can be applied by reshaping it into a vector.

△

For any p ∈ D, T ANOVA (p) is an approximation of the mapping T p , which transports µ 0 into µ p , evaluated in the values of X 0 . Now, for all the p i in the design of experiments which do not belong to the ANOVA sampling:

1. T ANOVA (p i ) is evaluated, allowing to dene µ ANOVA i discretized by (T ANOVA (p i ), µ 0 ) as a rst approximation of µ i .

2. The matrix Π i→i coupling µ ANOVA i and µ i is obtained by solving the linear program 2.4.7 with α = µ 0 , β = µ i and matrix C computed from T ANOVA (p i ) and X i .

The transport map is computed by barycentric projection

: ∀j = 1, . . . , M 0 , T i (X j 0 ) = T i→i (T 0→1 (X j 0 )) = (Π i→i X i ) j /µ j 0 . Measure µ i is now represented by (T i (X 0 ), µ 0 ).
Once the (T i (X 0 ) have been computed for i = 1, . . . , n s , the function T PGD (similar to u PGD from Eq. 2.2.12) can be built, and therefore the model T :

D → R M 0 ×2 is complete.
Similarly, a scalar ANOVA-PGD model ũ :

D → R relating p i to Ω u(x; p i )dx is constructed.
Now, for any choice of the parameters p ∈ D, u(•, p) can be approximated by:

∀x ∈ Ω, u(x, p) = ũ(p) M 0 M 0 j=1 K(∥ T (p) j -x∥).
(2.4.14) To recap, the complete workow is the following:

1. For i = 1, . . . , n s , the solutions u(•, p) are discretized and normalized, in order to obtain coordinates X i and weights µ i , which together form a discrete representation of measure µ i .

2. Solution u(•, p 0 ) is normalized, and discretized according to a decomposition similar as the ones performed for SPH, to obtain (X 0 , µ 0 )

3. Now, following the coupling route described in Fig. 2.40, the transport maps which transform µ 0 into µ i are obtained for all i such that p i ∈ D c , using optimal couplings and barycentric projection.

4. A rst model of the parametric transport map is created using the so-called "ANOVA terms" of the ANOVA-PGD.

5. This model is used to evaluate an approximation of the measure for all the i such that p i / ∈ D c . This estimated measure is coupled to the discrete representation of µ i , thus allowing to compute the transport map T i thanks to barycentric projection.

6. These new transport maps allow to complete the construction of the ANOVA-PGD model.

7.

A model of the magnitude of the solution is also constructed with the ANOVA-PGD.

8. These two models can be combined via Eq. 2.4.14 to approximate the parametric solution u.

This workow is also illustrated in Fig. 2.41.

ANOVA-based Optimal Transport for more general solutions

In the previous part, it was assumed that the solutions u had a localized support. This was not explicitly necessary, but it is important for the quality of the solution for two reasons: Approximating a constant in the whole domain with a decomposition in the spirit of SPH is extremely expensive and inecient.

If the localized features are small with respect to the global part, they will be erased by the approximations and error of the model, therefore using Optimal Transport becomes useless.

It was also assumed that the solutions were positive, because the conservation of mass is at the heart of Optimal Transport, which in the general case does not allow negativevalued measures. However, some works in the literature study such an extension [START_REF] Ambrosio | Gradient ow of the chapman rubinsteinschatzman model for signed vortices[END_REF][START_REF] Mainini | A description of transport cost for signed measures[END_REF], and although they are more interested in calculating the distance between objects, they can inspire ways to obtain transport maps.

The following proposes a revised version of the method to enable applying to more general solutions, by attempting to remove the two aforementioned constraints.

In order to deal with global solutions, such that supp(u) = Ω, the idea is to separate the solutions into two contributions, one that is well represented by global modes and which should be addressed directly using a method such as the ANOVA-PGD, and one which has a localized support and which can be dealt with using Optimal Transport.

Therefore, the solutions read: ∀(x, p) ∈ Ω × D, u(x; p) = u glob (x; p) + u loc (x; p). In the simple case, u glob (x; p) = u glob (x) does not depend on the parameters. In this case, u glob can be extracted as the median of all the solutions. Indeed, the part of the solutions corresponding to u loc are outliers in this context, and can be ltered out very eciently by the median. This idea can be extended to solve the problem when u glob has parametric dependen- cies. u glob is searched in the separated form

u glob (x; p) = m k=1 Ξ k (x)Ψ k (p) (2.4.16)
using the following methodology, derived from the PGD and based on the ltering of sparse outliers:

1. The terms of the sum are constructed successively in a greedy manner until convergence.

2. Ψ k is initialized randomly, and a rst approximation of Ξ k is obtained by solving:

min Ξ ∥ΞΨ k -u∥ 1 .
(2.4.17)

Advances in non-intrusive Model Order Reduction (2.4.18) 4. From Ξ k , Ψ k is updated by solving problem 2.4.17.

5.

Steps 3 and 4 are repeated until convergence.

6. If Ψ k or Ξ k is too sparse (usually due to unlucky initialization of Ψ k ), the term is rejected and another initialization is attempted at step 2.

7. After a few successive rejections, the algorithm stops.

Numerical experiments have shown that this algorithm yields a good results when the solutions fullll the requirements of the decomposition, which is nevertheless not guaranteed.

u loc = u -u glob is therefore a good candidate for the Optimal Transport strategy previously described because it has a localized support. However, in many cases, and even if u itself is positive, u loc has negative values. In order to be able to apply the ANOVA-based Optimal Transport strategy to negative-valued functions, the following describes how to perform couplings and compute transport plans between signed measures.

Once again, a solution u(•, p) is normalized by dividing it by its sum Ω u(•, p), which allows to keep the idea of mass conservation, even if negative mass is now introduced.

For this reason, the integral of u needs to have always the same sign, which is a special requirement for u loc without which divisions by 0 could occur.

Once the normalization is done, and the zero values are removed, the signed measure µ is divided into two parts, the positive part µ + > 0 and the negative part µ -> 0 such that µ = µ + -µ -. There are innitely many ways to obtain such a decomposition.

The one chosen is the trivial one, the so-called Jordan decomposition [START_REF] Fischer | Existence, uniqueness, and minimality of the jordan measure decomposition[END_REF], which is the only one in which µ + and µ -have disjoint supports.

When coupling two signed measures, mass can be coupled according to three mechanisms:

Positive mass of the rst measure can be coupled to positive mass of the second measure.

Negative mass of the rst measure can be coupled to negative mass of the second measure.

Positive mass of the "largest" measure, that is the one such that |u| u is the largest, can be coupled to negative mass of itself.

The coupling is thus divided in three steps which correspond to these three mechanisms:

1. Partial Optimal Transport [START_REF] Figalli | The optimal partial transport problem[END_REF] is used to cancel out the excess of positive and negative mass of the "largest" measure. By applying the barycentric projection, two transport maps are obtained, T + which transports the positive mass and T -which transports the negative mass.

The remaining negative and positive mass from the largest measure is added to the other measure in a way that it cancels itself out: for a unit of positive mass in position

x + coupled with a negative unit of mass in position x -, T + and T -(or their inverse depending on if the largest measure is the source or target) are extrapolated to estimate T + (x + ) and T -(x -), and the same unit of mass is added in both the positive and negative parts of the other signed measure, in position

T + (x + )+T -(x -) 2
.

Remark: For the same reason that transport maps do not always exist, the inverse of the transport maps obtained by barycentric projection may not exist, but they can be approximated by using barycentric projection in the opposite way.

△

This strategy, which allows to couple together signed measures and to compute the transport maps between them, can be introduced in the complete workow discussed previously, hence making it possible to apply it for more general types of solutions.

However, many special cases remain dicult to address, and require further investigation.

Numerical experiment

The test case proposed to evaluate the method is the construction of a parametric solution of a parametric wave equation, formulated as follows:

         ∆u(x, t) = ∂ 2 u ∂t 2 (x, t) ∀(x, t) ∈ Ω × [0, T ] ∂u ∂n (x, t) = 0 ∀(x, t) ∈ ∂Ω × [0, T ] u(x, 0) = u 0 (x; p) ∀x ∈ Ω (2.4.19)
where Ω = [0, 4] × [0, 4], T = 0.5 and u 0 (x; p)

= 3 k=1 e -(x 1 -p k ) 2 -(x 2 -k) 2 .
The problem is solved for 15 sets of parameters chosen following the design of experiments prescribed for ANOVA-PGD using Lagrange P1 nite elements. The solutions of this equation have a small support therefore no particular preprocessing is required, however they contain negative values so the coupling procedure described in the previous section is necessary. This problem is quite unsolvable with linear reduced bases because the location of the features depend strongly on the value of the parameters.

The parametric solution built is the 2D eld function of time and of the parameters, which make a total of 4 parameters.

Results for four dierent sets of parameters are shown and compared to the equivalent solution obtained using nite elements in Figs.

Data assimilation and Hybrid Modelling

In the realm of scientic modeling, the application of physics-based models has long been a cornerstone of industrial practices. These models, often formulated using partial dierential equations, provide a valuable framework for understanding and predicting the behavior of complex systems. However, despite their widespread use, these physics-based models are not without their limitations. Discrepancies between model predictions and real-world measurements are frequently observed, indicating the presence of inherent imperfections within the models.

The recognition of such disparities between model predictions and experimental data has given rise to the concept of hybrid modeling, which seeks to rene and enhance physics-based models through the integration of measured data. By incorporating realworld observations, it becomes possible to correct and enrich these models, bridging the gap between theoretical formulations and the actual behavior of the system under investigation.

Yet, this idea of using measurements to modify existing physics-based models presents several important questions. Firstly, can the collected data alone be used to construct a model, whether an entirely new model or the enrichment of an existing one, or is it necessary to acquire additional or improved data ? Secondly, how can an existing, well-posed physics-based formulation be modied to incorporate the assimilation of measurement data, without compromising its fundamental properties and integrity?

They questions have been addressed extensively in the literature, both the diculty of working with partial data [START_REF] Ayed | Learning dynamical systems from partial observations[END_REF][START_REF] Li | Propagation networks for model-based control under partial observation[END_REF] and using data to enrich models [186188]. They are revisited in this chapter, which delves into the realm of hybrid modeling, exploring dierent strategies to build, enrich, and correct physics-based models using measured data. The chapter is structured into three main sections, each addressing a distinct aspect of hybrid modeling.

The rst section addresses the challenge of modeling systems under partial observability. It investigates how incomplete or limited data can be leveraged to construct a reliable model that captures the essential dynamics of the system. Through this exploration of hybrid modeling techniques, this chapter aims to shed light on the challenges, approaches, and opportunities associated with integrating measured data into physics-based models. By harnessing the power of data, we can advance our understanding of complex systems, improve the reliability of predictions, and unlock new possibilities for scientic and industrial applications.

Modeling systems from partial observations

When proceeding with data for modeling purposes, whether in the context of fully data-driven or bybrid modeling, a recurrent issue concerns data accessibility. Sometimes, the considered system is not globally accessible, only a small part of it being reachable to perform measurements.

The present section addresses a conceptual issue, that will be discussed on an example simple enough to be fully understood, and at the same time complex enough to encompass all the modeling issues discussed here.

If there is a part of a system inaccessible for observation in which a loading that is applied cannot be either observed or measured, and that inuences the measures performed in the observable part of the system, dierent questions arise:

Is there a model connecting the observable input(s) to the corresponding output(s), knowing that they are impacted by the hidden dynamics of the system? Is it unique?

Under which conditions such a model could exist? How to nd it?

How to formulate it correctly? Is it well-posed and consistent?

How to learn it?

Which is the impact of these hidden dynamics on the learning process?

This section aims at revisiting the construction of models in the domains exhibiting partial observability, in both the steady and transient cases, while following a double approach: the usual algebraic formulation and the one concerned by Machine Learning approaches.

Data assimilation and Hybrid Modelling 79

On the existence of models relating observable features

Let us rst assume a large system, whose state is described by a number of state variables. The variables involved in the state description are supposed well dened.

However, the model governing the state or its time evolution is assumed unknown and the data describing the state is only observable and measurable on part of the system, remaining unattainable in the rest of the system. Previous analysis on the eld can be found in [7477].

For instance, in the case of the two-masses oscillator depicted in Fig. 3.1, it is assumed that the state is perfectly dened by the position and momentum of each mass, however only the state of the second mass is accessible (and thus, measurable). A natural question concerns the possibility to learn the model that governs the observable state (q 2 , p 2 ) while ignoring the state of the rst mass (q 1 , p 1 ). In what follows, this question is addressed using a quite generic algebraic rationale in two situations: a model that does not depend on time and a transient problem. First, generic settings are considered, then the reasoning will be applied to multiple-mass oscillators.

Time-independent problem

A generic linear time-independent model can be expressed from

KU = F, (3.1.1)
which, considering the observable variables U o and the internal ones U i , can be rewritten as

K oo K oi K io K ii U o U i = F o F i . (3.1.2)
Developing the last equation yields

K io U o + K ii U i = F i → U i = K -1 ii F i -K -1 ii K io U o , (3.1.3)
and introducing the resulting expression of U i into the development of the rst, it results (this is known as static condensation [START_REF] Wilson | The static condensation algorithm[END_REF] or Guyan reduction [START_REF] Chen | Guyan reduction[END_REF])

(K oo -K oi K -1 ii K io )U o = F o -K oi K -1 ii F i , (3.1.4)
which can be rewritten as

Koo U o = F o -Fi , ( 3
= (K oo -K oi K -1 ii K io ) Fi = K oi K -1 ii F i . (3.1.6)
Remark 1:

If F i = 0, a direct relation exists between U o and F o .

In the case of a 1D system in which only the borders of the interval are accessible (observable), U o and F o contains two components. By applying U o T = (1, 0), the resulting F o represents the rst column of Koo , and the solution F o related to U o T = (0, 1) will represent the second column of Koo .

In the same one-dimensional system, when F i ̸ = 0, there are two eective internal variables, the components of Fi . Thus, all the richness of F i boils up to these two components generating some sort of irreversibility: Fi can be obtained from F i , but the last one cannot be obtained from the former. The condensation of the internal degrees of freedom into the observable one produces an entropy increase in the theory of information sense: there are many micro-states F i associated to the macro-state Fi .

Computing these two eective internal variables just described requires an extracalculation. For example, if U o = 0, then F o = Fi .

Time-dependent problem

A general linear second-order dynamical system can be expressed from

M Ü + C U + KU = F, (3.1.7)
which, applying Fourier transform, becomes

-ω 2 MU + jωCU + KU = F, (3.1.8)
with j the imaginary (j 2 = -1) and U and F the Fourier transforms of U and F

respectively. The previous equation can be rewritten as

K * U = F, (3.1.9)
with K * = -ω 2 M + jωC + K, that can be separated in the same way considered in the time-independent case, but now, for each possible frequency (ω) involved in the loading and operating in the complex domain, leading to

K * oo U o = F o -Fi , (3.1.10)
which proves that all the discussion previously addressed in the time-independent case remains valid when the Fourier transform applies.

Thus, one could expect that a model relating observable variables might exist as well (and could be learned from collected data) in the time domain, under certain constraints, as the one referred in Remark 2 below, due to the dependence of Fi on the internal loading F i . This would imply the consideration of the history of the variables, which is naturally implicit in the Fourier transform.

Remark 2:

The just described rationale applies in the forced regime, i.e. far from the transient eects induced by the initial condition. In order to address transient regimes, the Laplace transform could be employed instead of the Fourier one. However, it is well known that the Laplace inverse transform is more challenging from the numerical point of view than Fourier's. It is also important to note that the Fourier transform of the internal loading considered in the training stage should remain invariant to ensure the validity of the learned model.

Neural Network-based modeling

In many cases Articial Intelligence, and more concretely Machine Learning, aims at extracting the model that relates measured inputs to the corresponding outputs [START_REF] Brunton | Data-driven science and engineering: Machine learning, dynamical systems, and control[END_REF][START_REF] Liu | Machine learning conservation laws from trajectories[END_REF]. In general, the measured output depends on the whole internal state.

For instance, in a structural dynamics problem where the loading (evolving in time)

constitutes the problem's input, the corresponding response is the displacement at each point and time; whereas the corresponding output-data is the measured displacement in a certain observable point of the structure.

In physics-based structural mechanics, the internal response (displacement at any location and time instant) is obtained by discretization of the continuum mechanics model, consisting of the momentum balance and the constitutive equations and; from this internal state, the output of interest is directly extracted at each time instant.

Alternatively, Machine Learning looks for the direct relation between observables, the input action and the measured response that, as just mentioned, can depend on the present and past values of a series of non observed internal variables [START_REF] Lee | Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders[END_REF].

Recurrent Neural Networks (RNN) [START_REF] Medsker | Recurrent neural networks[END_REF] and their Long-Short Time Memory counterparts (LSTM) [START_REF] Hochreiter | Long short-term memory[END_REF] address such situations by trying to model the time evolution of the internal state at the same time it constructs the model relating the observable input and output (action and response).

Addressing time-dependent problems in the time domain with K * = C + ∆tK and F * ,n = ∆tF n .

Modeling systems from partial observations

The sequencing of these equations can be written, inspired by the Dynamic Mode Decomposition [START_REF] Tu | Dynamic mode decomposition: Theory and applications[END_REF], in the matrix form K * U n , ..., U 1 = F * ,n , ..., F * ,1 + C U n-1 , ..., U 0 , (3.1.14) and by dening the extended vectors U and F, where the solution U is, in general, computed from the K matrix pseudo-inverse.

U T = [U n T , U n-1 T , ..., U 0 T ] F T = [F n T , F n-1 T , ..., F 1 T ] , ( 3
This algebraic system can be addressed by using the same rationale as the one applied before, but this time, the model explicitly involves the time evolution of the input(s) and output(s), reinforcing the result already obtained when using the Fourier transform.

Another alternative formulation, more aligned with the use of machine learning techniques that will be presented afterwards, consists in writing the explicit integration

CU n = ∆tF n -∆tKU n-1 + CU n-1 , (3.1.18)
that can be reformulated as

U n = AF n + BU n-1 , ( 3

Results for the analytical case in the linear setting

As mentioned, some numerical examples, simple enough to be perfectly understood, but complex enough to underline all the issues and methodological aspects just discussed.

The simplicity of the problem is only apparent. There are forces being applied to the internal masses, unknown and unobserved by the modeler, who, furthermore, totally ignores how many hidden masses are involved in the system. Three masses are considered in the present example, but they could come in any number, from one to the thousands.

When introducing all the system's degrees of freedom -in this case, the state of the three masses -in a model, the last one becomes larger but nally simpler to interpret and to learn, since all the data needed to properly describe the system is there, fully available. On the contrary, when considering only the data associated to one mass, while ignoring all the data related to all the other masses, the model seems simpler from its size, but very intricate nonetheless.

For this reason, and this is the motivation, the simplicity is only apparent, and allows for a more fruitful discussion on the issues and the conceptual questions previously addressed.

Modeling systems from partial observations

This way, the momentum no longer appears in the state variables, since it derives directly from the measurable position. Now, the partition between the internal and the observable degrees of freedom can be enforced:

A oo (ω) A oi (ω) A io (ω) A ii (ω) Q o (ω) Q i (ω) = R o (ω) R i (ω) (3.1.26)
that, following the aforementioned rationale, leads to

Ãoo (ω)Q o (ω) = R o (ω) -Ri (ω), (3.1.27) with Ri (ω) = A oi (ω)A -1 ii (ω)R i (ω) and Ãoo (ω) = A oo (ω) -A oi (ω)A -1 ii (ω)A io (ω),
where the same remarks that were previously discussed apply.

This is illustrated with a system composed of 3 identical masses (m The dynamical model reads:

1 = m 2 = m 3 = m), springs (k 1 = k 2 = k 3 = k)
        q1 ṗ1 q2 ṗ2 q3 ṗ3         =         0 1/m 0 0 0 0 -2k -2c/m k c/m 0 0 0 0 0 1/m 0 0 k c/m -2k -2c/m k c/m 0 0 0 0 0 1 m 0 0 k c/m -k -c/m                 q 1 p 1 q 2 p 2 q 3 p 3         +         0 k 1 l 1 -k 2 l 2 0 k 2 l 2 -k 3 l 3 0 k 3 l 3         +         0 F 1 (t) 0 F 2 (t) 0 F 3 (t)         , (3.1.28)
which, in the linear case and taking into account that k 1 = k 2 = k 3 = k and l 1 = l 2 = l 3 = l, after applying Fourier transform leads to:

        jω -1/m 0 0 0 0 2k jω + 2c/m -k -c/m 0 0 0 0 jω -1/m 0 0 -k -c/m 2k jω + 2c/m -k -c/m 0 0 0 0 jω -1/m 0 0 -k -c/m k jω + c/m                 q1 p1 q2 p2 q3 p3         =         0 0 0 0 0 kl         δ(ω) +          0 F1 0 F2 0 F3          , (3.1.29)
where the hat operator, •, refers to the Fourier transform of the masses positions, momentum and applied forces.

It is important to note that, in the nonlinear case described later on, since the spring stinesses depend on the spring elongation and the latter will obviously be dierent for each node -unlike here in the linear case -, the rst vector of the right hand member will contain three non-vanishing spring contributions:

k 1 l 1 -k 2 l 2 , k 2 l 2 -k 3 l 3 and k 3 l 3 .
By reordering the previous system, the position and momentum degrees of freedom can be grouped:

Tqq (ω) T qp T pq Tpp (ω)         q1 q2 q3 p1 p2 p3         =         0 0 0 0 0 kl         δ(ω) +          0 0 0 F1 F2 F3          , (3.1.30)
and then, the momentum degrees of freedom pi condensed into the ones related to the masses positions qi , as previously discussed:

T pq -Tpp (ω)T -1 qp Tqq (ω) Z q (ω) = S p (ω). which allows making the model involving the observable degree of freedom, q3 , explicit:

A 33 (ω) -A 31 (ω) A 32 (ω) A 11 (ω) A 12 (ω) A 21 (ω) A 22 (ω) -1 A 13 (ω) A 23 (ω) q3 (ω) = klδ(ω) + F3 + A 31 (ω) A 32 (ω) A 11 (ω) A 12 (ω) A 21 (ω) A 22 (ω) -1 F1 (ω) F2 (ω) , (3.1.33)
Rearranging this in a more compact manner yields: .34) which represents the system transfer function.

Ã33 (ω)q 3 (ω) = klδ(ω) + F3 (ω) + Fi3 (ω), (3.1
Now, the nal point concerns the data-driven model identication, that is, how to extract from the given data the dierent model components: Ã33 and Fi3 , for each involved frequency ω. In the last equation, the index i associated with Fi3 reects all the eects coming from the unresolved degrees of freedom (internal unobserved masses).

Conceptually, the system identication could proceed as follows:

1. The free response associated to F 3 = 0 (only the loads on the internal masses apply), q f 3 (t), is obtained (measured), the superscript • f refers to the fact that the observable mass remains free of loading.

2. Then, for a non-null (and measurable) applied loading on the observable mass, F 3 ̸ = 0, the system response q 3 (t) is recorded, which is now a consequence of all the loading terms involved in the right hand member of the previous equation.

3. The dierence between the forced and free displacement can be obtained from ∆q 3 (t) = q 3 (t) -q f 3 (t), allowing for the computation of its Fourier transform ∆q 3 (ω). In this example, the loading applied to the system reads: The 3-masses dynamical system is again considered for the time response scenario.

       F 1 (t) = 2 cos(2πt) F 2 (t) = 2 cos π 4 t F 3 (t) = 2 cos π 2 t , ( 3
The dynamical problem is integrated numerically to obtain the ground truth, that is, the reference solution. The computed data will be used to train the dierent neural networks, the RNN and the LSTM.

In both cases, the input data consists of the force F 3 and position q 3 in the previous time steps, that results in the surrogate H:

qi 3 = H           F i 3 F i-1 3 . . . F i-n+1 3      ,      q i-1 3 q i-2 3 . . . q i-n 3           , (3.1.37)
where qi 3 is the prediction of q 3 at time step i. As Eq. (3.1.37) reects, dierent memory lengths from the use of positive integer n, (n ≥ 0), are taken into account. For n ̸ = 0, an initialization issue occurs.

In the case considered here, the larger memory is the price paid ignore the internal forces, whose consequences on the observed variables are learned from the time evolution of the last.

The initialization can be carried out following two routes:

If one is interested in the forced regime, the long-time solution does not depend on the initialization.

To obtain the transient solution, one could consider a coarser model that updates the state from the just previous state until completing the rst n values. Then the LSTM can take over.

In the present case, as previously indicated, the focus is on proving under which conditions a model relating observable inputs and outputs exists, despite the existence of hidden dynamics, resulting in a noticeable larger memory. For that reason, in the simulations considered here, the rst n values are assumed known.

Using a simple recurrent neural network

Firstly, a RNN surrogate model is used, with the memory dened in Eq. (3.1.37) n = 2. The considered data for training comes from the integration of the dynamical system, in both the linear and nonlinear cases.

The data consists of 10000 states of the observable variables (coming as indicated from the standard integration of the dynamical system). This data are divided into two sets, the training and the testing ones, the former containing 80% of the points and the latter the remaining 20%. When addressing the linear case, the computed results are given in Fig. 3.9, with a MAPE of 0.84% in the training set and 1.33% in the testing set. The results in the nonlinear case are reported in Fig. 3.10, and again, for the sake of clarity, the associated absolute error is presented in Fig. 3

Hybrid modeling by learning a source term

After addressing the issue of collecting the right data to create models, arises the question of how to combine the data at hand with a physics-based model to create a hybrid model.

As mentioned in section 1.2.2, the PINN is a tempting route to answer this question in a simple manner: the equations coming from the physics can be embedded in the loss function to dene the solution of the well-posed physical problem, and then by adding constraints coming from data, the solution becomes a compromise between all the constraints, trying to model the data while retaining physical sense.

Although it is sucient for many applications, the downsides of such an approach is that it produces a hybrid solution but not a hybrid model. This may seem like a minor detail, but creating a hybrid model has two advantages compared to computing a hybrid solution:

The model can be used to obtain solutions, not only in the setting in which it was built, but also other settings (being cautious about extrapolation).

The nature of the correction applied to the model can in certain conditions be interpreted in a physical sense (this point is emphasized in section 3.3).

The present section proposes a hybridization strategy applied to classical discretization methods then improved within the PINN rationale, which corrects a physical model thanks to data measured on a real system. It is rst described, and then applied to an industrial case in the context of Hybrid Twins for power transformers.

Proposed methodology

Consider a system described by a state u(x, t) governed by a partial dierential equation in the form:

       L(u)(x, t) = f (p(x, t)) ∀(x, t) ∈ Ω × [0, T ] u(x, t) = f d (x, t) ∀(x, t) ∈ Γ d × [0, T ] ∂ n u(x, t) = f n (x, t) ∀(x, t) ∈ Γ n × [0, T ] u(x, 0) = u 0 (x) ∀x ∈ Ω (3.2.1)
This problem is assumed to have a unique solution u which perfectly ts within all the constraints, which could be computed methods such as nite elements, nite dierences or PINN. Once measurement data becomes available and a discrepancy between the simulation and reality is observed, the hybrid model comes into play.

Since the model dened by Eq. 3.2.1 is not able to represent reality, it must be wrong, or at least imprecise, and should therefore be modied to t the measurements obtained on the real system.

Therefore a dierent strategy than the one described in section 1.2.2 is proposed, where the enrichment is applied directly to the model rather than to the solver. Eq.

3.

2.1 is modied in the following fashion:

       L(u)(x, t) = f (p(x, t)) + g(u(x, t), p(x, t)) ∀(x, t) ∈ Ω × [0, T ] u(x, t) = f d (x, t) ∀(x, t) ∈ Γ d × [0, T ] ∂ n u(x, t) = f n (x, t) ∀(x, t) ∈ Γ n × [0, T ] u(x, 0) = u 0 (x) ∀x ∈ Ω (3.2.2)
where g is the enrichment term which will be used to learn the data-based correction.

Hybrid model using a classical solver

If enough data is available and the problem is low-dimensional, g could be learned from the solution generated by a classical solver. From measured data sampled in various locations and at dierent times according to a mesh that is ne enough:

u m (x i m , t i m ), i = 1, . . . , N m , an estimation of the residual of Eq. 3.2.1 is computed by approximating L with a discretization method on the mesh related to the sampling:

r i = L(u)(x i m , t i m ) -f (p(x i m , t i m )) (3.2.3)
Now, a regression model ĝ is created to link u and p to the residual:

∀i = 1, . . . , N m , ĝ(u(x i m , t i m ), p(x i m , t i m )) = r i (3.2.4)
Since a lot of data is available, the regression method of choice is usually the neural network, which can be trained extensively without real risk of overtting, unless the dimension of p is very high.

Hybrid modeling by learning a source term

Physics-based model

The IEC and IEEE dynamic thermal models [START_REF] Standards | Power transformer-part 7: Loading guide for mineraloil-immersed power transformers[END_REF][START_REF]Standards Association, Ieee guide for loading mineral-oilimmersed transformers and step-voltage regulators[END_REF] pertain to the standardization of the physical modeling of oil temperature in transformers. Despite being aware of their limitations, such as the lack of information regarding cooling conditions and changes in material properties, an IEC/IEEE thermal model will serve as a benchmark physicsbased model for the construction of the hybrid model and to compare the results.

Interested readers can refer to [START_REF] Alvarez | Transformer thermal capacity estimation and prediction using dynamic rating monitoring[END_REF] for a comprehensive overview of the physics-based thermal modeling process in a power transformer. The physics-based model utilized for this application is derived from the IEC and IEEE standards.

Based on the heat transfer equations, the proposed physics-based thermal model is composed of the following rst-order system:

P (K load , Θ W ) = C th (Θ T O ) d∆Θ T O dt + ∆Θ T O R th (Θ T O ) (3.2.6) Θ T O = T amb + ∆Θ T O (3.2.7) Θ W = Θ T O + ∆Θ OW (3.2.8)
where ∆Θ T O represents the predicted temperature rise of the top oil, Θ W corresponds to the estimated average temperature of the windings, P denotes the total estimated losses generated by the transformer, T amb represents the ambient temperature, ∆Θ OW represents the constant oil-to-winding temperature rise (provided by commissioning data), and C th and R th refer to the overall estimated thermal capacitance and resistance, respectively.

In order to enhance the accuracy of the model, it is important to consider that each parameter is temperature-dependent. The following observations can be made:

Parameter P depends on the transformer load and winding temperature due to the variation in material resistivity as a function of temperature.

The thermal capacitance C th is inuenced by the temperature of the top oil, as the heat capacity of the transformer's mineral oil changes with temperature.

The thermal resistance R th is updated because the viscosity of the oil changes with temperature.

The primary advantage of this physics-based model is its ability to obtain parameters solely based on commissioning data. The model only requires two inputs: the load factor K load and the ambient temperature T amb . Figure 3.12 provides a representation of the given inputs and compares the measured and estimated top-oil temperature for one of RTE's transformers. For the selected 2000 data points, the mean absolute percentage error (MAPE) is 10.4%. As anticipated, given the limited knowledge of the transformer's parameters and the model's simplicity, there is a signicant dierence between the measured and estimated topoil temperature, rendering it an inadequate estimation. Notably, when the top-oil temperature reaches 45°C (around the 1600th data point), a second group of fans is automatically activated, which is not accounted for in the proposed model. Although it is possible, as demonstrated in the example presented by [START_REF] Paulhiac | Dynamic thermal model for oil directed air forced power transformers with cooling stage representation[END_REF], to enhance the model's accuracy by considering variable cooling conditions, it often necessitates heuristic or regression-based parameter tuning, which cannot be universally applied to every transformer.

In conclusion, despite its inaccuracies, the proposed physics-based model serves as the reference for this application of the hybrid modeling concept, as it captures a signicant portion of the system's behavior and retains its physical signicance. The quantities measured are ambient temperature, load factor and top-oil temperature, at an hourly rate. The hourly load factor being a rapidly uctuating variable, it is calculated as the average of the values sampled every 5 minutes. In this work, the presence of a time delay in measurements is not considered, as it is negligible compared to the one-hour sampling period used.

As suggested by the method, the source term in Eq. 3.2.6 is enriched with an unknown function g:

P E = P (K load , Θ W ) + g K load , T amb , Θ T O , d∆Θ T O dt .
(3.2.9) Eq. 3.2.6 therefore becomes :

P E = C th (Θ T O ) d∆Θ T O dt + ∆Θ T O R th (Θ T O
) . that the discrepancy between the model and reality can not be entirely described by the extra source term g included in the model. Another approach may have produced better results, but it is also possible that some of the real behavior is due to unknown outside factors such as wheather conditions.

Hybrid modeling by identifying properties

In the pursuit of rening physics-based models through hybrid modeling, an alternative approach emerges: one that assumes the model itself is accurate but acknowledges the uncertainty surrounding the physical properties associated with the system under study. This uncertainty may arise from modeling errors or approximations, the use of coarse representations, or the presence of damage or aging within the system. This section delves into the latter scenario, focusing on the specic case of localized damage. Works aiming at correcting models which suered damage have been reviewed

in [START_REF] Alkayem | Structural damage detection using nite element model updating with evolutionary algorithms: a survey[END_REF].

When a system exhibits localized damage, accurately predicting its behavior becomes a complex task. However, by leveraging the framework of hybrid modeling, it becomes possible to correct the physical parameters associated with the system, enabling the development of a model capable of accurately capturing the behavior of the damaged system. Additionally, this approach oers the opportunity to not only correct the model but also to identify the nature and location of the damage itself.

The model correction methodology discussed in this section holds particular relevance in the context of Structural Health Monitoring (SHM) [START_REF] Balageas | Structural health monitoring[END_REF]. Extensive research eorts Data assimilation and Hybrid Modelling 97 have been dedicated to SHM, with a focus on investigating and addressing similar objectives. These studies often involve the analysis of the system in either the time domain [197202] or the frequency domain [203207] to detect and assess the presence of damage. Furthermore, the incorporation of machine learning techniques in SHM has gained considerable popularity, and notable applications of such techniques can be found in the literature [START_REF] Azimi | Data-driven structural health monitoring and damage detection through deep learning: State-of-the-art review[END_REF].

This section introduces a novel methodology centered around the idea of parameter correction. By updating the physical parameters within the existing model, it becomes possible to rene its predictive capabilities, aligning it more closely with the observed behavior of the damaged system. The aim is to obtain a model that not only accounts for the presence of damage but also accurately predicts its impact on the system's overall behavior.

Moreover, this parameter correction methodology oers the means to conduct an indepth analysis of the updated model. By comparing the updated model predictions with measured data, it becomes possible to identify the nature and location of the damage. This information is invaluable for understanding the underlying causes of the observed discrepancies and provides crucial insights for maintenance and repair strategies.

Key to the success of this approach is the utilization of sparse regularization techniques, which makes this methodology original. Sparse regularization allows for the selection of the most probable damage scenario from a range of possibilities that could explain the measurement data. By employing regularization techniques that favor sparse solutions, the methodology can eectively pinpoint the most likely damage locations and characterize the severity and extent of the localized damage within the system.

This approach is developed for two dierent types of numerical methods: the Finite Element Method in the context of linear equations and Physics-Informed Neural Networks.

Finite elements approach

The methodology is described on a linear elasticity problem. Let K be the original stiness matrix representing the structural behavior. K is an N ×N matrix, where N is the number of degrees of freedom in the structure. Let U be the nominal displacement vector obtained from the original model. U is an N × 1 vector. Let F be the applied loading vector. F is an N × 1 vector. The linear model describing the sytem reads KU = F.

(3.3.1) Eq. 3.3.1 is assumed to perfectly model the nominal system at the beginning of its lifecycle. However, after some time, a discrepancy is observed between the model's predictions and the strain or displacement measured on the system. This discrepancy is assumed to be due to a localized impact on the structure which deteriorates its mechanical properties. The objective is to nd a corrected stiness matrix ∆K and a corrected displacement vector ∆U that minimize the dierence between the predicted and measured displacements. This can be achieved by solving the following system of equations:

(K + ∆K)(U + ∆U) = F. where E is the number of elements involved in the mesh that covers the domain Ω, α e is the correction coecient for element e, and K e is the nominal stiness matrix related to element Ω e expressed in the global nodal numbering. This can be interpreted as a element-wise constant correction of the system's Young Modulus.

Alongside Eq. 3.3.4, the correction should also satisfy that the corrected displacement matches the measured displacement in the sensor locations:

U + ∆U = U m (3.3.6)
where the tilde • refers to the extraction of the part of the vector that correspond to the sensor locations and U m contains the displacement measurements in the sensor locations.

Remark: when strain measurements are performed instead of displacement measurements, Eq. 3.3.6 is adapted using the discrete (linear) strain operator E:

EU + E∆U = ε m (3.3.7)

△

The system of equations thus obtained enforces both the physics and the measurements while the constraints on the value of the Young's modulus are removed. However, in most cases, this creates an under-determined system which therefore has an ininite number of solutions. To circumvent this issue, the constraints of the Young's modulus can be introduced back in a relaxed form. Moreover, the assumption that the degradation of properties is localized motivates the use of the ℓ 1 -norm which promotes sparsity in the solution.

The optimization problem can thus be formulated as follows: where α = α + -α -similarly as in [START_REF] Tillmann | 2.17 Comparison between the solutions (displacement norm in mm) ob[END_REF].

The solution of this problem is not always completely suitable because there are sometimes still too many non-zero elements in α. When this happens, a relevant solution is to partition the domain into a certain number of parts, and to remove the degrees of freedom in α related to the parts in which α is smaller than a certain threshold, before solving the problem again. This can be iterated a few times until convergence.

Numerical experiment

To test the method, a numerical experiment is carried out by considering a homogeneous plate depicted in Fig. 3.14 and its damaged counterpart described in Fig. 3.15, from which simulations are performed to generate the "measured" data. The plate is discretized with a tetrahedral mesh, depicted in Figure 3.15, upon which the dierent mechanical elds (displacement, strain and stress) are computed using nite elements. The damaged elements are represented in red in Fig. 3.15; in these elements, the Young's modulus is reduced to 10% of its healthy state. For the sake of simplicity, the sensors are placed on nodes of the mesh, representd in blue in Fig.

PINN-based approach

A similar approach can be formulated in the PINN setting. In this case, the problem is not discretized and directly uses the continuous equations of solid mechanics. The following problem is considered:

               ∇ • σ + b = 0 x ∈ Ω u(x) = ϕ(x) x ∈ Γ d σ • n = ψ(x)
x ∈ Γ n ε = 1 2 (∇u + ∇u T ) σ = C : ε (3.3.11) where x is the spatial coordinates, Ω ⊂ R 2 the computational domain and Γ = Γ d ∪ Γ n its boundary. σ is the Cauchy stress tensor, ε the innitesimal strain tensor and C the fourth-order elasticity tensor. b is the force density, and ϕ and ψ are respectively the Dirichlet and Neumann boudary conditions applied to the system. C is dened from the Young's modulus E and Poisson's ratio ν such that: The solution is approximated using two fully-connected neural networks to represent u and σ, respectively denoted as û and σ.

To solve Eq. 3.3.11 with the PINN, the loss function L Φ is dened as:

L Φ = λ Ω L Ω + λ d L d + λ n L n + λ f L f (3.
L n = 1 N n Nn i=1 σ(x i n ) • n -ψ(x i n ) 2 (3.3.16) L f = 1 N Ω N Ω i=1 σ(x i Ω ) -C : ε(x i Ω ) 2 (3.3.17)
and λ Ω , λ d , λ n , λ f ∈ R * + are the weights associated to each term. ε is obtained by applying the adequate derivatives with respect to the inputs of neural network û.

Displacement measurements u m (x i m , t i m ), i = 1, . . . , N m (or strain measurements ε m (x i m , t i m ), i = 1, . . . , N m ) are collected on the real structure. To make the model coincide with this data, the Young's modulus E is once again assumed to have been degraded in a localized area and the displacement (or strain) is modied to match the value of the measurements in the sensor locations. This is performed by adding a third neural network Ê to represent the Young's modulus and replacing the loss function with the following hybrid one, in which two terms are added and one is modied: 3.18) where:

L H = λ Ω L Ω + λ d L d + λ n L n + λ f L f + λ m L m + λ r L r (3.
L f = 1 N Ω N Ω i=1 σ(x i Ω ) -Ĉ : ε(x i Ω ) 2 (3.3.19) L m = 1 N m Nm i=1 û(x i m ) -u m (x i m ) 2 (3.3.20) L r = 1 N Ω N Ω i=1 Ê(x i Ω ) -E(x i Ω ) (3.3.21)
and λ m , λ r ∈ R * + are the weights associated to each term. Ĉ is obtained from Eq. 3.3.12 applied to Ê. The term L r corresponds to the regularization discussed in the section 3.3.1. In this case, it aims at applying a L 1 -norm regularization to Ê -E in order to reduce its support, and therefore to have the nominal value of the Young's modulus in most of the domain.

Remark: when the strain is measured and not the displacement, L m becomes:

L m = 1 N m Nm i=1 ε(x i m ) -ε m (x i m ) 2 .
(3. This approach performs quite well, it reconstructs the full displacement eld from the sparse sensor measurements, and it is able to locate the damage even when the number of sensors is extremely small.

Data assimilation and Hybrid Modelling for representing data in order to build better reduced models. The technologies based on curve alignment and optimal transport proved to be eective in dealing with a priori irreducible data. While the methodologies dealing with curves are simple and easy to apply to any kind of problem, they are not automatic, therefore some work could be performed to decrease the number of tasks required to be performed by the modeler. Meanwhile, Optimal Transport allows an automatic processing of the data but it is very complex and the techniques proposed to allow to extend its domain of applicability are even more intricate.

3.1

Oscillator composed of two masses, two linear springs of stiness k 1 and k 2 , reference lengths l 1 and l 2 , and whose state is dened by the position and momentum of each mass (q 1 , p 1 , q 2 , p 2 ). . . . . . . . . . .

3.2

N-mass dynamical system . . . . . . . . . . . . . . . . . . . . . . . . .

3.3

Free response (F 3 (t) = 0): (left) q f 3 (t); and (right) qf 3 (ω). . . . . . . . .

3.4

Response: (left) q 3 (t); and (right) q3 (ω). . . . . . . . . . . . . . . . . .

3.5

Response dierence: (left) ∆q 3 (t); and (right) ∆q 3 (ω). . . . . . . . . .

3.6

Prediction of the observable position q3 (t) computed from a trained RNN with n = 2 (colors green and red mark the training and testing sets, respectively). It can be noted that the blue curve is not visible because it is almost exactly under the green and red curves. . . . . .

3.7

Prediction of the observable position q3 (t) in the nonlinear case, computed by a trained RNN with n = 3 (again, colors refer to the training and testing sets). It can be noted that the blue curve is not visible because it is almost exactly under the green and red curves. . . . . .

3.8

Error in the prediction of the observable position q3 (t)

0. 3

 3 Content 0.3 ContentThis thesis is divided in three chapters: rst, some tools and methodologies which are key components in the construction of Hybrid Twins are describes, then personal contributions in the eld of non-intrusive Model Order Reduction are presented and nally, personal contributions to hybrid modelling are described. 0.3.1 Chapter 1: Review of methods enabling Hybrid Twins Section 1: Intrusive and non-intrusive Model Order Reduction Engineering problems often requires fast and accurate simulations using complex models and large amounts of data. Although the development of hardware has allowed to improve the quality of simulations, there are still challenges in solving highdimensional problems due to numerical complexity and real-time processing limitations. Model Order Reduction (MOR) techniques, categorized as intrusive and nonintrusive methods, have gained attention to address these issues. Proper Orthogonal Decomposition (POD)[START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF][START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent ows[END_REF] extracts essential system characteristics and creates a reduced order model, signicantly reducing computation time. Proper Generalized Decomposition (PGD)[START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF][START_REF] Ladevèze | The latin multiscale computational method and the proper generalized decomposition[END_REF] is an intrusive method which employs the separation of variables to reduce the complexity of the resolution algorithm. Proper Orthogonal Decomposition with Interpolation (PODI)[START_REF] Tezzele | Shape optimization through proper orthogonal decomposition with interpolation and dynamic mode decomposition enhanced by active subspaces[END_REF][START_REF] Rama | Proper orthogonal decomposition with interpolation-based real-time modelling of the heart[END_REF] is a non-intrusive methods that handle parametric without modifying the solver. The sparse Proper Generalized Decomposition (sPGD)[START_REF] Ibáñez | A multidimensional datadriven sparse identication technique: the sparse proper generalized decomposition[END_REF] constructs a parametric solution using a sparse set of solutions computed by a PDE solver. Overall, MOR techniques provide ecient solutions for complex engineering problems in real-time applications. Section 2: Physics-Informed Neural Networks

Section 2 :

 2 ANOVA-PGD: A novel strategyThe previous section's workow presents new challenges and questions. In the context of s 2 -PGD, the choice of basis functions signicantly impacts the solution's quality, unlike in sPGD. Additionally, PGD-based regression methods currently use sparse design of experiments without exploiting smart sampling techniques, leaving untapped potential. These observations motivate the development of a new technique called ANOVA-PGD, which combines the approaches from the previous section with the ANOVA (Analysis of variance) decomposition [6466] strategy and its low-cost counterpart the anchored ANOVA [6769]. ANOVA-PGD utilizes a hierarchical orthogonal basis and suitable sampling, scaling linearly with the number of parameters involved. Section 3: Reduced Modelling of curves

Section 4 :

 4 Reduced Modelling based on Optimal TransportReduced Order Modelling in contexts such as the one described in the previous section become even more challenging when dealing with full-eld solutions in multiple dimensions. Linear reduced bases, commonly used in ROM algorithms, struggle to capture continuous displacements of solution features caused by phenomena like advection or localized parametric constraints. The limitations of reduced bases stem from their use of a Euclidean metric, which fails to consider misalignment or feature localization. To address this, Optimal Transport theory[16, 17] provides metrics designed for detecting distances between noteworthy features. This section introduces a methodology that combines the ANOVA framework and Optimal Transport to create parametric reduced models. 0.3.3 Chapter 3: Data assimilation and Hybrid Modelling Section 1: Modeling systems from partial observations This section addresses the issue of data accessibility in modelling. When only a portion of a system can be measured or observed, questions arise regarding the existence, formulation, and learning of a model that connects the observable inputs to the corresponding outputs, considering the inuence of hidden dynamics [7477]. The section revisits model construction in domains with partial observability, exploring both traditional algebraic formulations and Machine Learning approaches. 0.4 Scientic contributions Section 2: Hybrid modeling by learning a source term This section explores the combination of data and a physics-based model to create a hybrid model. While the PINN approach oers a simple solution by embedding physics equations in the loss function, it only provides a hybrid solution, not a hybrid model. Creating a hybrid model oers advantages such as versatility in dierent settings and the potential for physical interpretation of the model corrections. The section proposes a hybridization strategy established in the context of classical discretization methods and then improved within the PINN framework.

4 .

 4 S. Torregrosa, V. Champaney, A. Ammar, V. Herbert and F. Chinesta. Predicting high-delity data from coarse-mesh Computational Fluid Dynamics corrected using Hybrid Twins based on Optimal Transport. Mechanics & Industry. 5. A. Schmid, A. Pasquale, C. Ellersdorfer, V. Champaney, M. Raffler, S. Guevelou, S. Kizio, M. Ziane, F. Feist and F. Chinesta. PGD based Meta Modelling of a Lithium-Ion Battery for Real Time Prediction. Frontiers in Materials. 6. C. Ghnatios, X. Kestelyn, G. Denis, V. Champaney, and F. Chinesta.

1 Review of methods enabling Hybrid Twins 1 . 1

 111 7. M. Jacot, V. Champaney, S. Torregrosa, F. Chinesta and J. Cortial. Empowering Optimal Transport matching algorithm for the construction of surrogate parametric metamodel. Mechanics & Industry. 0.4.3 Conference proceedings 1. V. Champaney, A. Sancarlos, F. Chinesta, E. Cueto, D. González, I. Alfaro, S. Guevelou, J. L. Duval, A. Chambard, and P. Mourgue, Hybrid twins-a highway towards a performance-based engineering. part i: Advanced model order reduction enabling real-time physics, in ESAFORM 2021, Intrusive and non-intrusive Model Order ReductionEngineering often involves intricate models that demand swift and precise solutions, requiring processing of vast amounts of data. Meanwhile, from the requirements of their applications, these models must frequently be computed in real-time on deployed platforms. Despite the signicant advancements made in modelling, numerical analysis, discretization techniques, and computer science over the past decades, there are still numerous problems in science and engineering that cannot be solved due to their numerical complexity or limitations imposed by specic needs such as real-time processing on current technological platforms. Very often, these limitations come from the necessity to solve high-dimensional problems on very ne discretizations of complex domains. Recently, Model Order Reduction has attracted a lot of interest because of its ability to address previously intractable problems in numerical simulation. Model Order Reduction techniques are often classied in two groups: intrusive and non-intrusive methods. In this chapter, a review of a selection of intrusive and non-intrusive Model Order Reduction techniques is presented. First, Proper Orthogonal Decomposition (POD) is a general method used to extract the most essential characteristics of a system's behavior and express them with a low-dimensional representation through a set of basis vectors. One of the most common uses of POD is to project the governing equations of the system onto the reduced-order subspace spanned by the extracted basis. This produces a reduced order model that can be resolved instead of the initial system, in a usually much lower amount of time. Section 1.1.2 presents the Proper Generalized Decomposition (PGD), a method making use of separated representations to circumvent the curse of dimensionality. A greedy enrichment procedure coupled with a xed-point algorithm enable replacing a high-dimensional problem with a set of lower-dimensional problems solved iteratively. This strategy is especially useful when dealing with parametric problems or problems that can be separated in a multi-scale setting. Section 1.1.3 describes the Proper Orthogonal Decomposition with Interpolation (PODI), a non-intrusive technique based on a similar approach as POD, which 20 1.1 Intrusive and non-intrusive Model Order Reduction consists in extracting a reduced basis from collected data and performing an interpolation or regression to express the reduced coordinates as a function of the parameters of the problem. Unlike POD, this method does not require to modify the solver, hence it can be combined with any commercial software without needing access to the source code. Finally, section 1.1.4 introduces the sparse Proper Generalized Decomposition (sPGD), the non-intrusive counterpart of the PGD, which constructs a parametric solution from a set of solutions computed by a PDE solver based on a design of exper-

  Functions ψ k m , m = 1, . . . , M and k = 1, . . . , d are expressed from a standard approximation basis N k m , via coecients a k m :

32 2. 1

 321 Regularization of separated representations

Figure 2 .

 2 1 shows a plot of the original function f (x 1 , x 2 , x 3 = 0, x 4 = 0, x 5 = 0.7071). It can be noticed that the rs-PGD corrects the shape of the function in the indicated areas in Fig. 2.1, improving the performance of the regression. This improvement occurs over the whole ve-dimensional domain. Another part of the domain is shown in Figure 2.2 that depicts f (x 1 , x 2 , x 3 = -0.17069, x 4 = -0.17069, x 5 = -0.015517).

Figure 2 . 1 :

 21 Figure 2.1: Comparing the reference Eq.2.1.8 and its associated sPGD and rs-PGD regressions, at points (x 1 , x 2 , x 3 = 0, x 4 = 0, x 5 = 0.7071)

Figure 2 . 2 :

 22 Figure 2.2: Comparing the reference Eq.2.1.8 and its associated sPGD and rs-PGD regressions, at points (x 1 , x 2 , x 3 = -0.17069, x 4 = -0.17069, x 5 = -0.015517)

.1. 10 )

 10 by using a Chebyshev basis for the one-dimensional functions of the PGD. The above function is intended to be reconstructed in the domain D = [-1, 1] 3 . Furthermore, the training set is sampled using a sparse grid based on the Smolyak quadrature rule[START_REF] Smolyak | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF][START_REF] Kaarnioja | Smolyak quadrature[END_REF] with a level 3, employing the Clenshaw-Curtis univariate quadrature rule. Consequently, only these points are used for constructing the model, employing either the sPGD or the s 2 -PGD methodology. Figure 2.3 illustrates the mesh utilized for the training set.

Figure 2 .

 2 Figure 2.5 showcases the results of the s2-PGD methodology. As observed, the predictions are nearly perfect. Upon examining the solution, it becomes apparent that the model correctly identies four modes, representing four sums of the PGD decomposition.

Figure 2 . 3 :

 23 Figure 2.3: Plot of the original function and the training set (circles) used to construct the PGD models.

( 2 . 1 . 11 )

 2111 by using a Chebyshev approximation basis for the one-dimensional functions involved in the PGD constructor. The objective is to reconstruct the given function within the domain D = [-1, 1] 5 . The sampling for the training set contains 200 points. In addition, the Latin hypercube sampling is used to generate this random set of data. The training set comprises 290 points, and the Latin hypercube sampling is utilized to generate this randomized dataset. Additionally, a testing set consisting of 2000

Figure 2 . 4 :

 24 Figure 2.4: Problem dened in Eq.(2.1.10): Comparison of predicted sPGD values with the reference ones in the testing set (the black line represents a perfect prediction)

Figure 2 .

 2 Figure 2.7 presents the results of the s 2 -PGD methodology, showcasing excellent agreement between the real function and the proposed approach. Upon examining the s 2 -PGD solution, it becomes apparent that the model has accurately identied the non-zero elements. Additionally, this PGD solution requires 104 modes, representing 104 sums of the PGD decomposition, which can be further compressed by invoking the PGD again [129].

Figure 2 . 5 :Figure 2 . 6 :

 2526 Figure 2.5: Problem dened in Eq.(2.1.10): Comparison of predicted s 2 -PGD values with the reference ones in the testing set (the black line represents a perfect prediction)

Figure 2 . 7 :

 27 Figure 2.7: Problem dened in Eq. (2.1.11): Comparison of predicted s 2 -PGD values with the reference ones in the testing set (the black line represents a perfect prediction)

  -PGD: A novel strategy that allows dening the so-called Sobol sensitivity coecients S n S n = Var(u n ) Var(u) .

41 Design( 2

 412 of experiments (DoE) Since ANOVA-PGD requires calculating the expressions dened by the anchored ANOVA, a specic design of experiments is needed. Indeed, the solution needs to be evaluated in points which belong to a multidimensional cross D c centered in c. This cross is constructed from a set of one-dimensional domains D k c which are referred to as branches in the following manner: D k c = p ∈ D, ∀l ̸ = k, p l = c l

Figure 2 . 8 :

 28 Figure 2.8: Multidimensional cross and its branches for d = 3

Figure 2 . 9 :

 29 Figure 2.9: Example of a design of experiments for d = 3

5 .

 5 For each k = 1, . . . , d, the snapshot matrix U k ∈ R N ×n k is built by collecting the snapshots u(p k,j ) for j = 1, . . . , n k . The Singular Value Decomposition is performed as described in 1.1.1. The D k right singular vectors (columns of matrix R) which correspond to the largest singular values are extracted. The number of vectors D k is chosen according to an "energy" criterion on the singular values, as usual in POD. These vectors are interpolated in D k in a similar fashion as the one-dimensional functions of the ANOVA decomposition, and form the basis N k m as presented in 1.1.4. Index m in N k

Figure 2 . 10 :

 210 Figure 2.10: Comparing sPGD and ANOVA-PGD regressions

Figure 2 . 11 :Figure 2 . 12 :

 211212 Figure 2.11: Problem dened in Eq. (2.2.14): Comparison of predicted sPGD values with the reference ones in the testing set (the black line represents a perfect prediction)

Figure 2 . 13 :

 213 Figure 2.13: Problem dened in Eq. (2.2.15): Comparison of predicted SVR values with the reference ones in the training and testing sets

Figure 2 . 14 :

 214 Figure 2.14: Problem dened in Eq. (2.2.15): Comparison of predicted Neural Network values with the reference ones in the training and testing sets

Figure 2 . 15 :

 215 Figure 2.15: Problem dened in Eq. (2.2.15): Comparison of predicted ANOVA-PGD values with the reference ones in the training and testing sets

Figure 2 . 16 :

 216 Figure 2.16: Cylindrical indentation test of a lithium-ion cell

Figure 2 . 18 :

 218 Figure 2.18: Parametric dog bone specimen loaded in tension.

Fig. 2 .

 2 Fig. 2.19 shows two patterns of the Force-Displacement curve, obtained for two different choices of the Krupkowski parameters (blue and orange lines). A classical "euclidian" interpolation of these two patterns would result in the non-physical black dashed pattern.

Figure 2 . 19 : 51 K

 21951 Figure 2.19: Main issue encountered when using standard interpolations on non-aligned curves (the black dashed line represents the interpolation between the two colored lines).

Figure 2 . 20 :

 220 Figure 2.20: Behavior zones, transition and end points, for one function g(x).

Figure 2 . 21 :

 221 Figure 2.21: Curves g(x; p i ) related to dierent choices of the model features p i =(K i , n i , ε 0,i ), i = 1, . . . , n s .

Figure 2 . 22 :

 222 Figure 2.22: Functions g 1 i (y) ≡ g 1 (y; p i ) (left) and g 2 i (z) ≡ g 2 (z; p i ) (right), for i = 1, . . . , n s .

Figure 2 . 23 :

 223 Figure 2.23: Functions gi (x), for i = 1, . . . , n s , obtained after dilatation.

55

 55

Figure 2 .

 2 Figure 2.24: sPGD predictions (green line for training, red for testing) versus true curve (blue line).

S

  g(x;p) : (µ, σ) → Mg(x;p) , Σg(x;p) , S O(p) : (µ, σ) → MO(p) , ΣO(p) .

Figure 2 . 25 :

 225 Figure 2.25: Sketch of curve envelopes.

Fig. 2 .

 2 Fig. 2.26 shows the parametric curve and its statistical sensing, for a given choice of the input features distribution parameters. Condence Intervals have been computed using Algorithm 1, for the curve and the rupture point.

7: end for 8 :

 8 Using the previously built population {( MO(p j ) , ΣO(p j ) )} Ns j=1 , train a regression model for the statistical sensing of O(p) involved in the curve parametrization: S O(p) : (µ, σ) → MO(p) , ΣO(p) .

9 :

 9 Given a new couple (µ * , σ * ) of model features means and variances, corresponding to the features p * , one can obtain a Condence Interval CI at a given condence level for the output. For instance, at level 0.95, one can build a CI for the curve g(x; p * ): g(x; p * ) ∈ Mg(x;p * ) -2 Σg(x;p * ) , Mg(x;p * ) + 2 Σg(x;p * ) .

Figure 2 . 26 :

 226 Figure 2.26: Condence Interval of level 0.95 for the parametric Force-Displacement curve and for the rupture point, for a given choice of µ and σ.

Figure 2 . 27 :

 227 Figure 2.27: Parametric notched dog bone specimen loaded in tension (top and side views).

Figure 2 . 28 :

 228 Figure 2.28: Curves g i (x) = g(x; p i ) related to dierent choices of the model features p i =(R i , S i , h i ), i = 1, . . . , n s .

Figure 2 . 29 :

 229 Figure 2.29: Behavior zones, transition and end points, for one function g(x).

Fig. 2

 2 Fig. 2.31.

Figure 2 . 31 :

 231 Figure 2.31: Two dierent parameters congurations. Top: R = 7.59, S = 18.23, h = 0.84; bottom: R = 3.75, S = 5.58, h = 1.51 (all dimensions are provided in mm).The red zone is the part subject to rigid body constraints.

Fig. 2 . 61 Figure 2 . 32 :

 261232 Fig. 2.32 shows four snapshots of the displacement eld related to the specimens in Fig. 2.31, under axial tensile loading. The crack propagation follows two completely dierent patterns, drastically inuencing the Force-Displacement curve, as shown in Fig. 2.33.

Figure 2 . 33 :

 233 Figure 2.33: Force-Displacement curves corresponding to the two parameters congurations in Fig. 2.31.

Figure 2 . 34 :

 234 Figure 2.34: Functions g 2 i (x) of Fig. 2.30 (right) after clustering, for i = 1, . . . , n s .

2. 35 .

 35 Moreover, Fig. 2.36 shows the separating surface and classied points in the 3-dimensional parametric space.

Figure 2 . 35 :

 235 Figure 2.35: Confusion Matrices for the SVM classier (left: training data, right: test data).

Figure 2 . 36 :

 236 Figure 2.36: Parametric space and classied points (marker + + + is used for test points). The red plane is the separation surface.

Figs. 2 .

 2 Figs. 2.37 and 2.38 represent the plots of predictions for train and test, respectively, for 4 data points.

Figure 2 .

 2 Figure 2.37: sPGD predictions (green line) versus true curve (blue line) for training data.

64 2. 4

 644 Reduced Modelling based on Optimal Transport

Figure 2 .

 2 Figure 2.38: sPGD predictions (red line) versus true curve (blue line) for test data.

Figure 2 . 41 :

 241 Figure 2.41: The worklow follows these four steps: (A) The solutions are normalized and the reference solution is decomposed in SPH fashion (B) Couplings and transport maps are obtained following the arrows to visit the entire ANOVA sampling (C) The measures outside of the ANOVA sampling are approximated using the ANOVA model and corrected using the solutions, thus allowing to enrich the model (D) The model can be evaluated in real-time to approximate the solution for any values of the parameters.

73 3 .

 3 From Ξ k , Ψ k is updated by solving: min Ψ ∥Ξ k Ψ -u∥ 1 .

2 .

 2 Regular Optimal Transport is used to couple the remaining positive mass with the positive mass of the other distribution.

2. 4

 4 Reduced Modelling based on Optimal Transport 3. Regular Optimal Transport is used to couple the remaining negative mass with the negative mass of the other distribution.

Figure 2 . 42 :

 242 Figure 2.42: Solution of the problem for t = 0.2 and p = (2.9, 2.3, 1.2 using nite elements (left) and the surrogate model based on Optimal Transport (right).

Figure 2 . 43 :

 243 Figure 2.43: Solution of the problem for t = 0.22 and p = (2.7, 1, 1.7) using nite elements (left) and the surrogate model based on Optimal Transport (right).

Figure 2 . 44 :

 244 Figure 2.44: Solution of the problem for t = 0.32 and p = (2.4, 2.7, 2.4) using nite elements (left) and the surrogate model based on Optimal Transport (right).

Figure 2 . 45 :

 245 Figure 2.45: Solution of the problem for t = 0.36 and p = (1.3, 1.5, 1.7) using nite elements (left) and the surrogate model based on Optimal Transport (right).

78 3. 1

 781 Modeling systems from partial observationsThe second section of this chapter explores the process of enriching a physics-based model from data. It explores the learning of source terms or loading terms that can be incorporated into the existing equations, eectively augmenting the model to better match experimental measurements.Finally, the third section focuses on the correction of a physics-based model from data by identifying and rening the physical properties of the system. By comparing model predictions with measured data, it becomes possible to identify areas of discrepancy and adjust the model parameters to align more closely with reality. This section examines parameter calibration by promoting sparsity, showcasing its potential to enhance the accuracy and predictive power of physics-based models.

Figure 3 . 1 :

 31 Figure 3.1: Oscillator composed of two masses, two linear springs of stiness k 1 and k 2 ,reference lengths l 1 and l 2 , and whose state is dened by the position and momentum of each mass (q 1 , p 1 , q 2 , p 2 ).

.1.5) 80 3. 1

 801 Modeling systems from partial observations with Koo

Finally, it is

  possible to address time-dependent problems modeling but directly operating in the time domain, instead of operating in the Fourier domain as was considered before. For the sake of simplicity, this is illustrated with the rst order dynamical system C U + KU = F, (3.1.11) whose implicit time discretisation reads CU n + ∆tKU n = ∆tF n + CU n-1 , (3.1.12) with ∆t being the considered time step. This equation can be rewritten in the more compact form K * U n = F * ,n + CU n-1 , (3.1.13)

  and dampers (c 1 = c 2 = c 3 = c), with the springs having a reference length also identical (l 1 = l 2 = l 3 = l). Forces can be applied on both the internal masses (the rst two) as well as on the observable one, the third. The following values are considered: m = 0.5 kg, c = 0.8 Ns/m, k = 1 N/m, l = 1 m.

( 3 . 1 . 31 )

 3131 After separating the internal and observable degrees of freedom, it reads:  A 11 (ω) A 12 (ω) A 13 (ω) A 21 (ω) A 22 (ω) A 23 (ω) A 31 (ω) A 32 (ω) A 33 (ω)

4 .

 4 Finally, by means of the just calculated ∆q 3 and the Fourier transform of the measurable force F3 (ω), the model coecient Ã33 (ω) is learned from

.1. 36 ) 86 3. 1

 36861 Modeling systems from partial observationsThe free and forced responses and their Fourier transforms are depicted respectively in Figs.3.3 and 3.4. This loading is used to generate the synthetic data that will serve to identify the model's output q 3 (t) later on as a function of the observed load F 3 (t). During the training process of that model, F 1 (t) and F 2 (t) are fully ignored.

Figure 3 . 3 :

 33 Figure 3.3: Free response (F 3 (t) = 0): (left) q f 3 (t); and (right) qf 3 (ω).

Figure 3 . 4 :

 34 Figure 3.4: Response: (left) q 3 (t); and (right) q3 (ω).

Figure 3 .

 3 Figure3.5 shows the response dierence ∆q 3 (t) = q 3 (t) -q f 3 (t) and its Fourier transform ∆q 3 (ω) on the domain in which the dierence ∆q 3 (t) becomes almost stabilized, meaning the transient component almost vanishes.

Figure 3 . 5 :

 35 Figure 3.5: Response dierence: (left) ∆q 3 (t); and (right) ∆q 3 (ω).

3. 2

 2 Hybrid modeling by learning a source term Using a LSTM recurrent neural networkThe same linear and nonlinear dynamical systems are now processed by LSTM cells, with the same network parameters and initializations used for the RNN.

  .11, presenting a MAPE of 0.15% in the training set and 0.14% in the testing set. The error is again slightly larger in the training set for the same reasons given before. As expected, LSTM outperforms RNN for a large number of epochs. It was noticed that, by reducing the number of epochs, RNN outperforms LSTMs because convergence is more easily achieved using a lower number of parameters. The error in the linear case is larger, possibly due to the fact that it involves close to zero values which negatively impact the error calculation.

Figure 3 . 9 :

 39 Figure 3.9: Prediction of the observable position q3 (t) computed by a trained LSTM with n = 2 (the same color code is employed). It can be noted that the blue curve is not visible because it is almost exactly under the green and red curves.

Figure 3 . 10 :

 310 Figure 3.10: Prediction of the observable position q3 (t) in the nonlinear case, computed by a trained LSTM with n = 3 (same color code). It can be noted that the blue curve is not visible because it is almost exactly under the green and red curves.

Figure 3 . 11 :

 311 Figure 3.11: Error in the prediction of the observable position q3 (t) in the nonlinear case, computed by a trained LSTM with n = 3 (same color code).

Figure 3 . 12 :

 312 Figure 3.12: Physics-based top-oil temperature prediction

3

 3 Hybrid modeling by identifying propertiesA fully-connected neural network is used to model g and is trained by discretizing the equations with nite dierences and then minimizing the dierence between g and the residual of Eq. 3.2.6.

Figure 3 .

 3 Figure 3.13: Hybrid-based top-oil estimation using the source term correction approach

( 3 . 3 .2) 98 3. 3

 33983 Hybrid modeling by identifying propertiesExpanding the product and introducing Eq. 3.3.1 yields: K∆U + ∆KU + ∆K∆U = 0.

( 3 . 3 . 3 )( 3 . 3 . 4 )

 333334 To linearize the problem, the following approximation is considered:K∆U + ∆KU ≈ 0In order to maintain the physical mathematical properties and physical meaning of the stiness matrix, its correction term ∆K is approximated as:

min α,∆U ∥α∥ 1 s 1 T

 11 .t. K∆U + E e=1 α e K e U = 0 U + ∆U = U m (3.3.8)This problem can be transformed in the following linear program: (α + + α -)s.t. α + ≥ 0 α -≥ 0 K∆U + E e=1 (α + e -α - e )K e U = 0 U + ∆U = U m (3.3.10)

Figure 3 . 14 :

 314 Figure 3.14: Structural model considered for illustrating the local model correction methodology.

3. 15 .

 15

100 3. 3 Figure 3 . 15 :

 1003315 Figure 3.15: Damaged elements (red) and location of the 30 sensors (blue).

Figure 3 . 16 :

 316 Figure 3.16: Identied elements most aected by the damage (in red). (A) Top view; (B) bottom view.

Figure 3 . 17 : 2 • 100 = 6 . 2 • 100 = 2 . 2 • 100 = 3 .

 317210062100221003 Figure 3.17: Top view of x (left), y (center) and z (right) components of the displacement eld obtained with the nominal model (top); reference model that takes into account the real damaged region (middle) and corrected model (bottom).

Figure 3 . 18 :Fig. 3 .

 3183 Figure 3.18: Top view of the dierence in absolute value between the corrected and the reference displacement eld, for x (left), y (middle) and z (right) components.

Figure 3 . 19 :

 319 Figure 3.19: Top view of sensor locations for dierent congurations.

Figure 3 . 20 :

 320 Figure 3.20: Relative error U Corr x,y,z -U Ref x,y,z 2 / U Ref x,y,z 2 with respect to number of sensors in the three directions (x in blue, y in red and z in yellow). The 6 congurations are the ones shown in Fig. 3.19.

  method, a numerical experiment is carried out by considering a homogeneous plate with a Young's modulus of 1M P a depicted in Fig.3.21a and its damaged counterpart described in Fig.3.21b, from which direct simulations are performed to generate the "measured" data.

Figure 3 . 22 :

 322 Figure 3.21: Problem conguration

Figure 3 . 24 :

 324 Figure 3.24: Comparison between the real (left) and computed (right) Young's modulus

Figure 3 . 25 :

 325 Figure 3.25: Sensors location (left) for dierent congurations, and computed Young's modulus (right)

  Studies should be performed to understand how to apply Optimal Transport based methodologies to general data.Partial observability: This part of the thesis involved a theoretical study on creating models from partial observability. This investigation explored the creation of machine learning models for dynamical systems in scenarios where the system cannot be completely observed. While a simple case study was presented, future research could expand upon this work to explore more complex scenarios and test the limitations of modeling under partial observability.Hybrid models: This part of the thesis focused on constructing hybrid models from physics-based models and measurement data. The developed strategies aimed at creating hybrid models rather than hybrid solutions, emphasizing explicability and generalization of the model. Two methodologies were proposed 110 Conclusion to enrich models based on partial dierential equations. One approach introduced an unknown as a source term in the problem, while the other aimed to identify properties of the physical system. Both methods successfully improved the model's t to the data. Furthermore, the second methodology demonstrated the utility of model explicability in Structural Health Monitoring and damage detection. These techniques were formulated using classical discretization methods as well as Physics-Informed Neural Networks. Leveraging the exibility of PINNs could yield original and powerful extensions to the techniques presented in the thesis. In conclusion, this thesis has laid the foundation for advancing the eld of hybrid twinning by introducing novel methodologies and exploring various aspects of model construction, representation, data-driven considerations, and hybridization. The contributions made in this work provide a valuable framework for further research and applications in the development of hybrid twins for addressing industrial challenges. List of Figures 2.40 Example in two dimensions showing the coupling route followed in the ANOVA sampling. Each red arrow represents a step of measure coupling, from which the transport map is extracted and then composed with the previous ones to obtain a global transport map from point c.2.41 The worklow follows these four steps: (A) The solutions are normalized and the reference solution is decomposed in SPH fashion (B) Couplings and transport maps are obtained following the arrows to visit the entire ANOVA sampling (C) The measures outside of the ANOVA sampling are approximated using the ANOVA model and corrected using the solutions, thus allowing to enrich the model (D) The model can be evaluated in real-time to approximate the solution for any values of the parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.42 Solution of the problem for t = 0.2 and p = (2.9, 2.3, 1.2 using nite elements (left) and the surrogate model based on Optimal Transport (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.43 Solution of the problem for t = 0.22 and p = (2.7, 1, 1.7) using nite elements (left) and the surrogate model based on Optimal Transport (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.44 Solution of the problem for t = 0.32 and p = (2.4, 2.7, 2.4) using nite elements (left) and the surrogate model based on Optimal Transport (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.45 Solution of the problem for t = 0.36 and p = (1.3, 1.5, 1.7) using nite elements (left) and the surrogate model based on Optimal Transport (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Top view of x (left), y (center) and z (right) components of the displacement eld obtained with the nominal model (top); reference model that takes into account the real damaged region (middle) and corrected model (bottom). . . . . . . . . . . . . . . . . . . . . . . . . 101 3.18 Top view of the dierence in absolute value between the corrected and the reference displacement eld, for x (left), y (middle) and z (right) components. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 3.19 Top view of sensor locations for dierent congurations. . . . . . . . . 102 3.20 Relative error U Corr x,y,z -U Ref x,y,z 2 / U Ref x,y,z 2 with respect to number of sensors in the three directions (x in blue, y in red and z in yellow). The 6 congurations are the ones shown in Fig. 3.19. . . . . . . . . . . . . 103 3.21 Problem conguration . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 3.22 Comparison between the real (left) and computed (right) x-displacement105 3.23 Comparison between the real (left) and computed (right) y-displacement105 3.24 Comparison between the real (left) and computed (right) Young's modulus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 3.25 Sensors location (left) for dierent congurations, and computed Young's modulus (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

  Réduction de modèle avancée et approches basées sur les données pour la construction de jumeaux numériques augmentés par la physique xv favorisation de la parcimonie dans la régression[START_REF] Brunton | Discovering governing equations from data by sparse identication of nonlinear dynamical systems[END_REF][START_REF] Ibañez | Some applications of compressed sensing in computational mechanics: model order reduction, manifold learning, data-driven applications and nonlinear dimensionality reduction[END_REF]. Cette section propose des méthodologies de régression robustes et ecaces pour les représentations séparées. Deux techniques sont présentées : la rs-PGD, qui combine la régularisation des normes L2 et L1, et la s 2 -PGD, une méthode de régression doublement parcimonieuse utilisant la régularisation Lasso. Context and Motivations . . . . . . . . . . . . . . . . . . . . . . . . . .

	xx		Contents
	Section 3 : Modèles réduits de courbes 2 Advances in non-intrusive Model Order Reduction	27
	2.1	Regularization of separated representations . . . . . . . . . . . . . . .	28
	La prédiction de courbes entières (notamment des fonctions scalaires du temps ou d'un Section 1 : Modélisation à partir d'observations partielles
	paramètre) par des modèles paramétriques de substitution présente des dicultés. Le 2.1.1
	non-alignement des courbes, qui peut résulter de fréquences ou fenêtres d'échantillon-Cette section aborde la question de l'accessibilité des données dans la modélisation. . . 31
	Ces deux techniques reposent initialement sur un mécanisme commun : l'ajout d'un terme de pénalisation lors de la dénition des problèmes de minimisation. Suivant nages diérentes, et le décalage de points critiques dans les courbes présentant des caractéristiques communes (par exemple, la transition élastique-plastique en méca-nique), font que les méthodes d'interpolation classiques ne fournissent pas des résul-Lorsque seule une partie d'un système peut être mesurée ou observée, des questions se posent concernant l'existence, la formulation et l'apprentissage d'un modèle qui relie Contents 2.1.3 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . 32 2.2 ANOVA-PGD: A novel strategy . . . . . . . . . . . . . . . . . . . . . . Introduction 38 les entrées observables aux sorties correspondantes, en tenant compte de l'inuence l'algorithme de la PGD et plus particulièrement de la sPGD, la solution recherchée tats cohérents du point de vue de la physique ; cela nécessite d'appliquer aux courbes 2.2.1 ANOVA decomposition and anchored ANOVA . . . . . . . . . 39 des dynamiques cachées [7477]. est exprimée sous une forme séparée, une somme de produits de fonctions unidimen-un prétraitement approprié. Par ailleurs, lorsque des bifurcations se produisent dans Dans un premier temps, la question est étudiée dans le cadre de formulations al-2.2.2 ANOVA-PGD . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 sionnelles, puis chacune de ces fonctions est calculée en résolvant successivement un certain nombre de problèmes unidimensionnels. Dans la s 2 -PGD (resp. rs-PGD), un terme pénalisant la norme L1 (resp. une combinaison des normes L1 et L2) est in-l'espace paramétrique, le couplage avec des algorithmes de partitionnement (cluste-ring) [70, 71] et de classication [72] améliore la précision du modèle de substitution. gébriques classiques. En utilisant la condensation statique (ou réduction de Guyan) Introduction 2.3 Reduced Modelling of curves . . . . . . . . . . . . . . . . . . . . . . . 49 [78, 79], les conditions selon lesquelles un modèle peut être reconstitué à partir d'ob-2.2.3 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . 44 0.1 Context and Motivations
	troduit an de réduire le surapprentissage qui intervient lorsque peu de données sont En combinant des traitements des courbes basés sur l'alignement de caractéristiques servations partielles sont analysées. La même étude est conduite en régime dynamique
	disponibles. remarquables avec une stratégie basée sur un clustering en amont et une classication en appliquant également cette condensation dans l'espace de Fourier. 0.1 2.3.1 Data alignment and uncertainty propagation . . . . . . . . . . 50
		2.3.2
			.	91
		3.2.2	Application to an industrial problem . . . . . . . . . . . . . . .	93
	3.3	1.1.1 Hybrid modeling by identifying properties . . . . . . . . . . . . . . . . Proper Orthogonal Decomposition . . . . . . . . . . . . . . . . 96
		3.3.1	1.1.2 Finite elements approach . . . . . . . . . . . . . . . . . . . . . Proper Generalized Decomposition . . . . . . . . . . . . . . . . 97
		3.3.2	1.1.3 PINN-based approach . . . . . . . . . . . . . . . . . . . . . . . 103 Proper Orthogonal Decomposition with interpolation . . . . . .
			xix

Par ailleurs, la s 2 -PGD ajoute à cette régularisation une étape de sélection des fonctions de base qui permettent une bonne représentation de la solution. Une fois les fonctions adéquates sélectionnées, la solution est recalculée sans pénalisation, mais dans un espace de dimension réduit mais adapté aux données qui permet donc d'éviter le surapprentissage tout en représentant bien la solution.

Les exemples proposés montrent que ces stratégies permettent une amélioration de la capacité prédictive de la sPGD.

Section 2 : ANOVA-PGD : Une stratégie originale Les techniques présentées dans la section précédente soulèvent de nouveaux dés et de nouvelles questions. Dans le contexte de la s 2 -PGD, le choix des fonctions de base a un impact signicatif sur la qualité de la solution, contrairement à la sPGD, qui donne des résultats dépendants de l'espace engendré par la base mais indépendant du choix de la base en elle-même. En outre, les méthodes de régression basées sur la PGD utilisent actuellement des plans d'expériences parcimonieux qui n'exploitent pas de technique d'échantillonnage intelligent ; cela laisse de côté un potentiel inexploité. Ces observations motivent le développement d'une nouvelle technique appelée ANOVA-PGD, qui combine les approches de la section précédente avec la décomposition ANOVA (Analyse de la variance) [6466] et sa version "low-cost" dite "ancrée"

[6769]

.

La décomposition ANOVA permet d'approximer une fonction de manière hiérarchique grâce à une base orthogonale qui contient des fonctions d'un nombre de variables croissant. A l'ordre 0, la fonction est approchée par une constante puis enrichie à l'ordre 1 par une somme de fonctions unidimensionnelles. Grâce à la stratégie ANOVA "ancrée", le calcul de ces fonctions nécessite un échantillonnage approprié, qui sépare la contribution de chacune des variables. Cet échantillonnage a la particularité de croître linéairement avec le nombre de paramètres impliqués, ce qui en fait une technique très intéressante pour des problèmes en grande dimension.

Par ailleurs, l'étude de ces contributions séparées des diérents paramètres permet d'établir des bases construites à partir des données et adaptées à chaque paramètre.

Ainsi, l'approximation à l'ordre 1 est combinée avec une technique basée sur la représentation séparée munie de ces bases adaptées an de construire une solution paramétrique la plus dèle possible, même dans des cas assez fortement non-linéaires et en très grande dimension. xvi Réduction de modèle avancée et approches basées sur les données pour la construction de jumeaux numériques augmentés par la physique en aval, il est montré que les modèles réduits permettant de prédire des courbes à partir de paramètres donnés sont largement améliorés. Ces méthodes ne sont pas totalement automatiques mais simples à mettre en ÷uvre pour améliorer drastiquement la qualité des résultats, notamment en s'assurant de la bonne représentation des caractéristiques de courbes qui sont chargées de sens physique.

Ensuite, la propagation d'incertitudes

[START_REF] Lee | A comparative study of uncertainty propagation methods for black-box-type problems[END_REF]

dans ces modèles est étudiée ; supposant que les paramètres d'entrée du modèle sont connus de manière incertaine et modélisés par une loi de probabilité, un intervalle de conance est calculé pour la courbe de sortie.

Section 4 : Réduction de modèle basée sur le transport optimal La réduction de modèle dans des contextes tels que celui décrit dans la section précédente devient encore plus dicile lorsque les solutions sont des champs complets en plusieurs dimensions. Les bases réduites linéaires, couramment utilisées dans les algorithmes de réduction, peinent à capturer les déplacements continus des caractéristiques de la solution causés par des phénomènes tels que l'advection ou des chargements paramétriques localisés. Les limites des bases réduites proviennent de leur utilisation d'une métrique euclidienne, qui ne prend pas en compte le désalignement ou la localisation des caractéristiques. Pour y remédier, la théorie du transport optimal fournit des métriques conçues pour détecter les distances entre les caractéristiques remarquables.

Le transport optimal

[16, 17]

construit notamment une manière d'interpoler entre des fonctions en prenant en compte cette problématique de caractéristiques localisées. Cette interpolation est ecace quand seulement deux fonctions sont concernées, mais plus dicile quand un grand échantillonnage est concerné, particulièrement dans un vaste domaine paramétrique. Pour répondre à cette problématique, une stratégie reprenant l'échantillonnage de l'ANOVA "ancrée" et permettant de visiter l'espace paramétrique de proche en proche est proposée. Cette combinaison entre l'ANOVA et le transport optimal permet d'agrandir le domaine d'application de la réduction de modèle. Cependant, les outils liés au transport optimal nécessitent des conditions strictes sur les fonctions étudiées, notamment qu'elles soient positives. Pour remédier à ce problème, une deuxième stratégie plus complète est proposée pour étendre encore un peu davantage ce domaine d'application. Réduction de modèle avancée et approches basées sur les données pour la construction de jumeaux numériques augmentés par la physique xvii Chapitre 3 : Assimilation de données et modélisation hybride Les idées développées dans le cadre algébrique simple sont ensuite utilisées pour résoudre un problème similaire mais non-linéaire à l'aide de méthodes d'apprentissage machine adaptées. En eet, les RNN et LSTM sont capables de modéliser la contribution des parties du système qui n'ont pas été observées, dans des conditions similaires à celles établies plus tôt. Section 2 : Modélisation hybride par apprentissage de terme source Cette section explore la combinaison de données et d'un modèle basé sur la physique pour créer un modèle hybride. Bien que l'approche PINN ore une solution simple en intégrant des équations physiques dans la fonction de perte, elle ne fournit qu'une solution hybride, et non un modèle hybride. La création d'un modèle hybride ore des avantages tels que la généralisation à d'autres conditions et potentiellement une interprétation physique des corrections du modèle. Pour enrichir le modèle à partir des données, une inconnue est ajoutée au modèle sous la forme d'un terme source. Dans le cadre des méthodes classiques de discrétisation d'équations aux dérivées partielles, ce terme est appris à partir de l'approximation des opérateurs diérentiels appliqués aux données. En revanche, cela n'est possible que lorsque les données sont abondantes puisqu'elles doivent être associées aux n÷uds d'un maillage susamment n pour approximer correctement les opérateurs. Pour atténuer cette contrainte très forte, le problème est reformulé à l'aide d'un PINN. Dans cette nouvelle formulation, les dérivées sont calculées de manière exacte sans contrainte particulière, et la solution peut donc être obtenue à partir de données parcimonieuses. Section 3 : Modélisation hybride par identification de propriétés Dans le but d'enrichir les modèles physiques par le biais de la modélisation hybride, une approche alternative émerge : une approche qui suppose que le modèle lui-même est précis mais qui reconnaît l'incertitude entourant les propriétés physiques associées au système étudié. Cette incertitude peut résulter d'erreurs de modélisation ou d'approximations, de l'utilisation de représentations grossières ou de la présence de dommages ou de vieillissement dans le système. Cette section se penche sur ce dernier scénario, en xviii Réduction de modèle avancée et approches basées sur les données pour la construction de jumeaux numériques augmentés par la physique 0.1.1 From virtual to digital twins . . . . . . . . . . . . . . . . . . . 0.1.2 In between models and data . . . . . . . . . . . . . . . . . . . . 0.1.3 Physics in real-time: Model Order Reduction . . . . . . . . . . 0.1.4 Data and learning . . . . . . . . . . . . . . . . . . . . . . . . . 0.1.5 Hybridization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.2 Objectives of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . 0.3 Content . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.3.1 Chapter 1: Review of methods enabling Hybrid Twins . . . . . 0.3.2 Chapter 2: Advances in non-intrusive Model Order Reduction . 0.3.3 Chapter 3: Data assimilation and Hybrid Modelling . . . . . . 0.4 Scientic contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.4.1 Publications in international journals . . . . . . . . . . . . . . . 0.4.2 Submitted articles . . . . . . . . . . . . . . . . . . . . . . . . . 0.4.3 Conference proceedings . . . . . . . . . . . . . . . . . . . . . . 0.4.4 International conferences . . . . . . . . . . . . . . . . . . . . . 1 Review of methods enabling Hybrid Twins 1.1 Intrusive and non-intrusive Model Order Reduction . . . . . . . . . . . 1.1.4 Sparse Proper Generalized Decomposition . . . . . . . . . . . . 1.2 Physics-Informed Neural Networks . . . . . . . . . . . . . . . . . . . . 1.2.1 Neural Networks as PDE solvers . . . . . . . . . . . . . . . . . 1.2.2 Neural Networks for hybrid solutions . . . . . . . . . . . . . . . rs-PGD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.1.2 s 2 -PGD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . POD modes extraction . . . . . . . . . . . . . . . . . . . . . . . 53 2.3.3 Curves reconstruction . . . . . . . . . . . . . . . . . . . . . . . 54 2.3.4 Real-time calibration . . . . . . . . . . . . . . . . . . . . . . . . 55 2.3.5 Statistical model derived by parametric curves . . . . . . . . . 55 2.3.6 Data alignment and data clustering . . . . . . . . . . . . . . . 56 2.3.7 Clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 2.3.8 Curves reconstruction and classication . . . . . . . . . . . . . 62 2.4 Reduced Modelling based on Optimal Transport . . . . . . . . . . . . 64 2.4.1 Short introduction to Optimal Transport . . . . . . . . . . . . 65 2.4.2 ANOVA-based Optimal Transport methodology . . . . . . . . . 67 2.4.3 Numerical experiment . . . . . . . . . . . . . . . . . . . . . . . 74 3 Data assimilation and Hybrid Modelling 77 3.1 Modeling systems from partial observations . . . . . . . . . . . . . . . 78 3.1.1 On the existence of models relating observable features . . . . . 79 3.1.2 Results for the analytical case in the linear setting . . . . . . . 82 3.1.3 Results of RNN and LSTM time simulations in both the linear and the nonlinear settings . . . . . . . . . . . . . . . . . . . . . 87 3.2 Hybrid modeling by learning a source term . . . . . . . . . . . . . . . 90 3.2.1 Proposed methodology . . . . . . . . . . . . . . . . . . . . . .

  4 Scientific contributions 0.4.1 Publications in international journals 7. T. Loreau, V. Champaney, N. Hascoet, J. Lambarri, M. Madarieta, I. Garmendia, and F. Chinesta, Parametric analysis and machine learningbased parametric modeling of wire laser metal deposition induced porosity, International Journal of Material Forming, 15 (2022), p. 33 [105] 8. V. Champaney, A. Pasquale, A. Ammar, and F. Chinesta, Parametric curves metamodelling based on data clustering, data alignment, pod-based modes extraction and pgd-based nonlinear regressions, Frontiers in Materials, 9 (2022), p. 363 [96] 9. D. Di Lorenzo, V. Champaney, C. Germoso, E. Cueto, and F. Chinesta, Data completion, model correction and enrichment based on sparse identication and data assimilation, Applied Sciences, 12 (2022), p. 7458 [110] 10. V. Champaney, V. J. Amores, S. Garois, L. Irastorza-Valera, C. Gh-Scientic contributions 19. D. Di Lorenzo, V. Champaney, J. Marzin, C. Farhat, and F. Chinesta, Physics informed and data-based augmented learning in structural health diagnosis, Computer Methods in Applied Mechanics and Engineering, 414 (2023),

	p. 116186 [111]
	0.4.2 Submitted articles
	1.

natios, F. J. Montáns, E. Cueto, and F. Chinesta, Modeling systems from partial observations, Frontiers in Materials, 9 (2022)

[START_REF] Champaney | Modeling systems from partial observations[END_REF]

11. M. El Fallaki Idrissi, F. Praud, V. Champaney, F. Chinesta, and F. Meraghni, Multiparametric modeling of composite materials based on nonintrusive pgd informed by multiscale analyses: Application for real-time stiness prediction of woven composites, Composite Structures, 302 (2022), p. 116228

[START_REF] Fallaki Idrissi | Multiparametric modeling of composite materials based on nonintrusive pgd informed by multiscale analyses: Application for real-time stiness prediction of woven composites[END_REF]

12. S. Torregrosa, V. Champaney, A. Ammar, V. Herbert, and F. Chinesta, Hybrid twins based on optimal transport, Computers & Mathematics with Applications, 127 (2022), pp. 1224 [109] 13. C. Ghnatios, V. Champaney, A. Pasquale, and F. Chinesta, A regularized real-time integrator for data-driven control of heating channels, Computation, 10 (2022), p. 176 [94] 14. A. Pasquale, V. Champaney, Y. Kim, N. Hascoët, A. Ammar, and F. Chinesta, A parametric metamodel of the vehicle frontal structure accounting for material properties and strain-rate eect: application to full frontal rigid barrier crash test, Heliyon, 8 (2022), p. e12397 [100] 15. M. Jacot, V. Champaney, F. Chinesta, and J. Cortial, Parametric damage mechanics empowering structural health monitoring of 3d woven composites, Sensors, 23 (2023), p. 1946 [102] 16. A. Sancarlos, V. Champaney, E. Cueto, and F. Chinesta, Regularized regressions for parametric models based on separated representations, Advanced Modeling and Simulation in Engineering Sciences, 10 (2023), p. 4 [93] 17. A. Runacher, M.-J. Kazemzadeh-Parsi, D. Di Lorenzo, V. Champaney, N. Hascoet, A. Ammar, and F. Chinesta, Describing and modeling rough composites surfaces by using topological data analysis and fractional brownian motion, Polymers, 15 (2023), p. 1449 [106] 18. C. Ghnatios, D. di Lorenzo, V. Champaney, E. Cueto, and F. Chinesta, Optimal velocity planning based on the solution of the euler-lagrange equations with a neural network based velocity regression, Discrete and Continuous Dynamical Systems-S, (2023), pp. 00 [95] 0.4 S. Torregrosa, V. Champaney, A. Ammar, V. Herbert and F. Chinesta. Physics-based Active Learning for design space exploration and surrogate construction for multiparametric optimization. Communications on Applied Mathematics and Computation. 2. M. Kazemzadeh-Parsi, A. Pasquale, D. Di Lorenzo, V. Champaney, A. Ammar and F. Chinesta. NURBS-based shape parametrization enabling PGD-based space separability: methodology and application. Finite Elements in Analysis & Design. 3. A. Schmid, A. Pasquale, C. Ellersdorfer, M. Ziane, M. Raffler, V. Champaney, F. Feist and F. Chinesta. Application of PGD separation of space to create a reduced order model of a lithium-ion cell structure. Frontiers in Materials.

  1.2 Physics-Informed Neural Networks learning techniques heavily rely on vast amounts of labeled data, which may not always be available in domains governed by physical laws. PINNs bridge this gap by incorporating prior knowledge of physical laws into the learning process, allowing for accurate predictions and capturing underlying physics even with limited data. This

	section provides a concise introduction to PINN, describing how the governing equa-
	tions and constraints are combined with a neural network architecture to leverage
	data while ensuring physical consistency.
	1.2.1 Neural Networks as PDE solvers
	Consider a system described by a state u(x, t) governed by a partial dierential equa-
	tion in the form:

  The objective of this section is to propose regression methodologies that are robust, general, ecient, and accurate, specically designed for separated representation settings. To achieve this goal, two techniques are presented and analyzed. The rst technique is called rs-PGD, and is based on an Elastic Net regularized formulation combining Ridge (L2 norm regularization) and Lasso (L1 norm regularization) regressions. The rs-PGD technique unites these regularizations with a rich approximation basis, and ensures specic solutions with smaller coecients to prevent overtting.

Subsequently, the section introduces the s 2 -PGD technique, which is a doubly sparse regression method. This technique utilizes Lasso regularization which enforces the sparsest possible solution by selecting the most signicant contributors in each dimension to the solution approximation.

Table 2 .

 2

1: Parametric ranges

Table 2 . 2 :

 22 Parametric ranges

  X (p), p = (p 1 , . . . , p d ): curves regression model; 2. N s : number of training points for the statistical surrogate model S O(p) ; 3. N MC : number of Monte Carlo sampling points. , Σ O(p) ): regression model for mean and variance of curve descriptor O(p).

	Advances in non-intrusive Model Order Reduction	57
	Algorithm 1 Statistical sensing based on parametric curves
	Input:
	1. f Output:
	(M O(p) 1: for j = 1, . . . , N s do
	2:

Table 2 .

 2

3: Parametric ranges.

  2.42, 2.43, 2.44 and 2.45. The results that are shown correspond to the test set, they were not used to build the model.

  in the nonlinear case, computed by a trained RNN with n = 2 (the same color code applies). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Prediction of the observable position q3 (t) computed by a trained LSTM with n = 2 (the same color code is employed). It can be noted that the blue curve is not visible because it is almost exactly under the green and red curves. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 3.10 Prediction of the observable position q3 (t) in the nonlinear case, computed by a trained LSTM with n = 3 (same color code). It can be noted that the blue curve is not visible because it is almost exactly under the green and red curves. . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 3.11 Error in the prediction of the observable position q3 (t) in the nonlinear case, computed by a trained LSTM with n = 3 (same color code). . . 92 3.12 Physics-based top-oil temperature prediction . . . . . . . . . . . . . . 95 3.13 Hybrid-based top-oil estimation using the source term correction approach 96 3.14 Structural model considered for illustrating the local model correction methodology. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 3.15 Damaged elements (red) and location of the 30 sensors (blue). . . . . . 100 3.16 Identied elements most aected by the damage (in red). (A) Top view; (B) bottom view. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 3.17

	List of Figures
	3.9

gourmandes en calculs et ont du mal à gérer les systèmes à grande échelle et l'incertitude. Elles nécessitent souvent l'utilisation de techniques de réduction de modèle[13] pour devenir adéquates avec les contraintes de fonctionnement en temps réel. En parallèle, les approches basées sur les données sont prometteuses lorsque des modèles précis ne sont pas disponibles, mais elles sont confrontées à des problèmes liés au coût des données[START_REF] Mauroy | Koopman operator in systems and control[END_REF], aux risques d'extrapolation[8, 9], et au manque d'explicabilité et de certications[10, 11].Par conséquent, la combinaison des deux approches semble être le choix optimal, car elle permet de trouver un équilibre entre leurs avantages et leurs inconvénients. En intégrant les perspectives basées sur la physique et les données, nous pouvons tirer parti des forces de chaque approche pour obtenir de meilleurs résultats en matière d'ingénierie.Un avantage signicatif de cette alliance est la réduction des besoins en données pour la construction de modèles. Cette réduction est obtenue en exploitant des lois connues ou en se concentrant sur l'approximation de l'écart entre le modèle et la réalité. En outre, l'utilisation de modèles basés sur la physique permet d'expliquer les aspects fondamentaux et les prédictions qui en résultent, ce qui permet de certier les résultats et d'atténuer les problèmes d'extrapolation.En outre, la réduction des besoins en données n'est pas seulement une conséquence vii viii Réduction de modèle avancée et approches basées sur les données pour la construction de jumeaux numériques augmentés par la physique de l'utilisation des lois physiques, mais aussi de la manière dont la physique peut guider la sélection des lieux et des instants optimaux pour la collecte des données.Ceci est particulièrement évident dans l'apprentissage actif[12], où l'incorporation de connaissances existantes basées sur la physique peut améliorer son ecacité.En combinant la compréhension basée sur la physique avec des techniques basées sur les données, les ingénieurs peuvent exploiter la puissance des deux approches, ce qui conduit à des pratiques d'ingénierie plus ecaces et plus perspicaces.Dans cette thèse, divisée en trois chapitres, nous nous penchons sur les outils, les méthodologies et les contributions personnelles qui sont des éléments cruciaux dans la construction de jumeaux hybrides. Le chapitre 1 passe en revue les méthodes permettant d'obtenir des jumeaux hybrides, notamment les techniques intrusives et non intrusives de réduction de modèle (MOR) et les réseaux neuronaux informés par la physique (PINN)[13]. Le chapitre 2 explore des évolutions en réduction de modèle non intrusive, en se concentrant sur la régularisation des représentations séparées, une nouvelle stratégie ANOVA-PGD basée sur l'ANOVA (analyse de la variance)[14, 15], la modélisation réduite des courbes, et l'utilisation du transport optimal[16, 17] dans la construction de modèles réduits. Enn, le chapitre 3 traite de l'assimilation de données et de la modélisation hybride, en discutant de la modélisation de systèmes à partir d'observations partielles et en présentant des approches de modélisation hybride basées sur l'apprentissage d'un terme source et l'identication de propriétés.En combinant ces chapitres, cette thèse contribue au développement, à l'amélioration et à l'application de méthodologies qui permettent la construction de jumeaux hybrides. Ces méthodologies comblent le fossé entre la science des données et la simulation numérique, en répondant aux dés industriels actuels. Elles orent des modèles paramétriques de haute délité et de haute dimension, accélèrent les modèles basés sur la physique grâce à la réduction de modèle et l'apprentissage machine, s'attaquent aux dés de la construction de modèles sous observabilité partielle, et créent des modèles hybrides qui combinent des principes basés sur la physique avec des données. Grâce à ces contributions, cette thèse vise à faire progresser le domaine de l'ingénierie en tirant parti de la puissance des données et de la simulation pour résoudre des problèmes complexes du monde réel.

se concentrant sur le cas spécique de l'endommagement localisé. Cet aspect permet notamment d'appliquer la méthode au contrôle de santé des structures[START_REF] Balageas | Structural health monitoring[END_REF].En mettant à jour les paramètres physiques du modèle existant, il devient possible d'aner ses capacités prédictives, en l'alignant plus étroitement sur le comportement observé du système endommagé. L'objectif est d'obtenir un modèle qui non seulement tient compte de la présence d'un dommage, mais qui prédit également avec précision son impact sur le comportement global du système.En outre, cette méthodologie de correction des paramètres ore les moyens de mener une analyse approfondie du modèle mis à jour. En comparant les prédictions du modèle mis à jour avec les données mesurées, il devient possible d'identier la nature et la localisation des dommages. Ces informations sont précieuses pour comprendre les causes sous-jacentes des écarts observés et fournissent des indications cruciales pour les stratégies de maintenance et de réparation.La clé du succès de cette approche est l'utilisation de techniques de régularisation parcimonieuses, ce qui rend cette méthodologie originale. La régularisation parcimonieuse permet de sélectionner le scénario d'endommagement le plus probable à partir d'un éventail de possibilités pouvant expliquer les données de mesure. La méthodologie peut ainsi eectivement localiser les endroits où les dommages sont situés et caractériser la gravité et l'étendue des dommages dans le système.Cette approche est développée pour deux types diérents de méthodes numériques : la méthode des éléments nis dans le contexte des équations linéaires et les réseaux neuronaux informés par la physique.

Remerciements Acknowledgements

warmly welcoming me to Singapore, providing administrative assistance, and collaborating on joint work. I extend my heartfelt thanks to Jean-Yves Marzin

Reduced Modelling based on

The problem considered in this section is the construction of the parametric solution u : (x; p) ∈ Ω × D → R of a parametric partial dierential equation, where Ω is a compact convex subset of R 2 (for the sake of simplicity in this example but it could also be R 3 ) and D ⊂ R d . Some of the notations are those dened in section 2.2. The solutions corresponding to a design of experiments such as the one described in section 2.2.2 are computed using an existing solver. The discretized counterpart of u(•; p i ), i = 1 . . . n s is denoted by u(p i ) ∈ R N , a vector which contains the nodal values of u and the nodal coordinates are stored in matrix X ∈ R N ×2 . To transform them into measures suitable for the theory described in section 2.4.1, only the M i points inside the support of u are considered and the solutions are normalized so that their integral over Ω is 1. X i ∈ R M i ×2 contains the points in the support of u(•; p i ) and v i the corresponding nodal values. Then, the discretized measure is represented by vector µ i ∈ R M i given by:

and corresponding coordinates X i .

The idea is to compute the transport maps from the center point c and then create two regressions using the ANOVA-PGD method: one for the transport maps and the other one for the values Ω u(x; p i )dx in order to recover the unnormalized eld.

The main issue with this strategy is related to the risks of operating far away from the tangent point. The transport map between solutions which are really dissimilar can produce unwanted behaviors. However, the particular structure of the ANOVA Advances in non-intrusive Model Order Reduction 69 sampling oers an interesting way to circumvent this issue based on the following example.

Suppose a problem with one parameter with D = [0, 2] and a 3 sample points and associated measures (p 1 = 0, µ 1 ), (p 2 = 1, µ 2 ), (p 3 = 2, µ 3 ). For the sake of the example, the chosen reference point is µ 1 . Then, to take into account that p 2 lies between p 1 and p 3 , the transport map from µ 1 to µ 3 , denoted here as T 1→3 , is computed as:

(2.4.12)

This is done by solving the transport problem between µ 1 and µ 2 and then the one between µ 2 and µ 3 , but never the one between µ 1 and µ 3 .

The same principle is applied in multiple dimensions in the ANOVA sampling. The On the half branch associated with the positive values of p 1 (as shown in Fig. 2.40 with points p 1 and p 2 ), the coupling algorithm is the following:

1. Starting from p 0 = c, the matrix Π 0→1 coupling µ 0 and µ 1 is obtained by solving the linear program 2.4.7 with α = µ 0 , β = µ 1 and matrix C computed from X 0 and X 1 .

2. The transport map is computed by barycentric projection: ∀j = 1, . . . , M 0 , T 1 (X j 0 ) = T 0→1 (X j 0 ) = (Π 0→1 X 1 ) j /µ j 0 where exponent j refers to the extraction of the j-th line of the matrix. The discretized version of measure µ 1 can now both be represented by (X 1 , µ 1 ) and (T 1 (X 0 ), µ 0 ) where T 1 (X 0 ) = (T 1 (X j 0 )) j=1,...,M 0 .

The latter is preferred because it allows to compute easily the composition of transport maps.

Learning in the Fourier space

The N -mass dynamical system considered here is illustrated in Fig. 3.2. The modeling includes inertia, elastic and damping behaviors.

The state of each mass is represented by z i = (q i , p i ), q i and p i being the i-th mass position and momentum respectively. The system state is described by the extended

The usual model, coming from Newton's equation, can be expressed by

where matrix T includes the system properties, masses, spring stiness and viscosity of the dampers. On the other hand, J is a constant vector (in the linear case addressed below) and F contains the external forces applied on the dierent masses, appearing at the odd positions in vector F (an explicit form of that matrix and those vectors will be given later). By dening the eective loading S = J + F, and T = -T + jωI, the matrix form that separates the degrees of freedom related to the measurable position (noted by q) and the derived momentum (p) is now written:

Since S q (ω) = 0, Z p (ω) can be expressed in terms of Z q (ω):

that, introduced into the second equation, leads to:

which can be reshaped into the more compact form:

with A(ω) = T pq -Tpp (ω)T -1 qp Tqq (ω), Q(ω) ≡ Z q (ω) and R(ω) ≡ S p (ω).

Modeling systems from partial observations

The RNN consists of a single layer with one output qi 3 , in reference to Eq. (3.1.37).

The network parameters and the initialization choices are the ones reported in [START_REF] Glorot | Understanding the diculty of training deep feedforward neural networks[END_REF].

The algorithm is trained during 1500 epochs, although the use of fewer epochs also leads to similar results.

The linear problem considers, once more:

, and l 1 = l 2 = l 3 = 1m, expressing the applied loading the following way:

with t max = 500s. This loading is used to generate the synthetic data that will serve afterwards to identify the model q 3 (t) as a function of the observed load F 3 (t). During the training process of that model, F 2 (t) and F 3 (t) are again completely neglected.

The computed results from the trained network are given in Fig. 3.6, the mean absolute percentage errors (MAPE) are 1.38% on the training set and 2.18% on the testing set. It can be noted that the blue curve is not visible because it is almost exactly under the green and red curves.

The same RNN (now with n = 3 in reference to Eq. (3.1.37)) is now employed to tackle a nonlinear dynamical system, with similar parameters to the ones considered in the linear case, except in what concerns the springs stinesses, now given by:

with k 01 = k 02 = k 03 = 10 N/kg, α = 10 -4 m -1 (arbitrary, although carefully tuned to maintain the stability of the simulation) and where ∆l • is the elongation of the corresponding spring, i.e. ∆l 2 = q 2 -q 1 -f 2 , ∆l 3 = q 3 -q 2 -f 3 and ∆l 1 = q 1 -l 1 .

Data assimilation and Hybrid Modelling 89

The considered loading reads:

The results concerning the nonlinear dynamical system are reported in Fig. 3.7, and, for the sake of clarity, the associated absolute error in reported in Fig. It can be noted that the blue curve is not visible because it is almost exactly under the green and red curves. Once the model ĝ is built, it can be introduced in the solver according to Eq. 3.2.2 to produce the enriched solution.

Although appealing, this route is rarely applicable in practice, except to tackle ODE problems, which are one-dimensional and do not require full-ed measurements in space.

Hybrid model strategy based on the PINN

In most cases, the data availability is quite reduced. In this case, it is impossible to approximate L(u m ). However, the PINN formulation allows to circumvent this issue.

Following the PINN rationale, the solution is approximated by a neural network û(x, t).

An additional neural network ĝ : R × R d → R is dened to approximate function g. Both neural networks are trained simultaneously by minimizing the loss function L H dened in Eq. 3.3.18. However the rst term L Ω is adapted from Eq. 1.2.3 to include the change in the model and allow to train ĝ:

When the training is complete, the hybrid solution u is obtained, but also function g and therefore a hybrid model. In certain cases, g could be analyzed by the modeler to understand the nature of the enrichment which was necessary for the model to match the data. It could also be used to compute the solution for a dierent loading or source by changing p, or for dierent initial and boundary conditions.

Application to an industrial problem

Power transformers play a vital role in transmitting electrical energy and demand special focus. Monitoring and predicting the oil temperature of power transformers is crucial because it signicantly impacts the expected lifespan of the insulation and therefore of the transformer. In particular, excessively high or low oil temperatures can cause irreversible harm to the transformer.

Since the usual lifespan of power transformers is between 20 and 60 years, even if a transformer initially has a well calibrated model, it is not expected to match reality forever. The hybrid modeling route is therefore very appealing, to enrich and correct the model whenever the measurements performed on the real system diverge from the model's predictions. This could have two major advantages:

Quantifying the correction necessary to match the measured data is a tool for diagnosis, because the discrepancy between the model and reality is often due to aging or degradation.

Having a model which matches reality enables prognosis and testing the impact of various scenarios of operation on oil temperature and therefore aging of the system.

List of Figures

2.1

Comparing the reference Eq.2.1.8 and its associated sPGD and rs-PGD regressions, at points (x 1 , x 2 , x 3 = 0, x 4 = 0, x 5 = 0.7071) . . . . . . .

2.2

Comparing the reference Eq.2.1.8 and its associated sPGD and rs-PGD regressions, at points (x 1 , x