# On the construction of non-intrusive multifidelity models for computer codes with time-series output: Comparison of three paradigms on a transient thermal problem

H. Nasser, G. Perrin, R. Chakir, S. Demeyer, Julien Waeytens

► To cite this version:

H. Nasser, G. Perrin, R. Chakir, S. Demeyer, Julien Waeytens. On the construction of non-intrusive multifidelity models for computer codes with time-series output: Comparison of three paradigms on a transient thermal problem. Journal of Computational Physics, 2025, 543, pp.114411. 10.1016/j.jcp.2025.114411. hal-05327760

HAL Id: hal-05327760

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

Submitted on 23 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# On the construction of non-intrusive multifidelity models for computer codes with time-series output: Comparison of three paradigms on a transient thermal problem

H. Nasserª, G. Perrinª, R. Chakirª, S. Demeyerᵇ, J. Waeytensª

## Abstract

This study evaluates three multifidelity meta-modeling approaches for time-varying systems: the Non-Intrusive Reduced Basis method (NIRB), the Gaussian Process Regression (GPR), and the Long Short-Term Memory Recurrent Neural Networks (LSTM). The focus is on predicting, in a small data context, the output of an expensive numerical code (high-fidelity version) that takes a parameter vector as input and produces a time-dependent function as output, leveraging one (or more) less expensive code(s) referred as low-fidelity versions. The goal is to analyze the strengths and limitations of these approaches in this context, and to make these three approaches accessible, even to non-specialists. To this end, we compare these methods based on prediction performance (accuracy), ease of hyper-parameter tuning, offline versus online computational costs, ability to quantify uncertainty, but also based on their complexity of implementation. The analysis also explores the advantages of using multiple codes and strategies for distributing computational costs across each code. Numerical experiments are conducted using a thermal model designed to predict surface heat fluxes and temperatures of a multi-layer wall based on its geometry and the thermal properties of its materials.

**Keywords:** multifidelity approach, Meta-model, Functional outputs, Reduced Basis, Gaussian Process Regression, Recurrent Neural Networks

## 1. Introduction

Numerical simulations can be computationally expensive, for example the resolution of a 3D-hydrothermal model can take more than two days for one evaluation. In such a configuration, it becomes extremely challenging, to exploit such costly simulators to carry out uncertainty quantification [1], sensitivity analyses [2, 3] or solve optimization problems [4, 5]. A common approach to reduce the computational cost is to use a meta-model (or surrogate model), that is to say a fast-to-evaluate mathematical function designed to replicate the behavior of the complex model of interest. However, building a meta-model typically requires multiple evaluations of the costly simulator, which can also represent a significant cost, depending on the complexity of the code's input-output relationship. To address this, multifidelity (MF) techniques combine simulations at different fidelity levels

<footnote>ªUniversité Gustave Eiffel, COSYS, IMSE, 77420 Champs-sur-Marne, France</footnote>

<footnote>ᵇLaboratoire National de Métrologie et d'Essais, 78190 Trappes, France</footnote>

<footnote>Email addresses: hadi.nasser@univ-eiffel.fr (H. Nasser), guillaume.perrin@univ-eiffel.fr (G. Perrin), rachida.chakir@univ-eiffel.fr (R. Chakir), severine.demeyer@lne.fr (S. Demeyer), julien.waeytens@univ-eiffel.fr (J. Waeytens)</footnote>
during the learning phase to reduce the number of evaluations of the most accurate model, which is called high-fidelity (HF) model, while maintaining equivalent prediction quality. These less expensive models are known as low-fidelity (LF) models, and are generally divided into several categories [6] such as: 1) simplified physical model, which are generally obtained by modifying the governing equations to be solved, as realised in [7, 8], 2) simplified discrete solution, for example by using a coarser discretization as lower fidelity approximation, as realised in [9, 10].

A review of MF models is presented in [6], exploring methods for integrating different fidelity levels into surrogate models. The paper classifies publications based on their application domains, surrogate model techniques, and fidelity management strategies, focusing on the trade-offs between computational cost and accuracy. It also provides recommendations for standardizing the reporting of cost and accuracy metrics to enhance comparability across studies. The authors conclude that MF models can achieve significant computational time savings, up to 90%, depending on the application, while maintaining or even improving predictive accuracy. In [11], a Bayesian method using a multi-level auto-regressive model with Gaussian processes is proposed for fast scalar output prediction. Applied to oil reservoir simulation with varying mesh qualities, this approach enhances prediction accuracy for complex models while significantly reducing computational time and cost. More recent study [12] introduce a Bayesian approach for parameter estimation across multiple fidelity levels is introduced, modeling relationships with Gaussian processes and a scaling factor. This method simplifies covariance matrix inversion by providing closed-form formulations, reducing computational complexity. Applied to thermodynamics, it improves accuracy and lowers computational cost but relies on high correlation between fidelity levels. The authors later extend this work with a recursive approach [13], integrating simulations hierarchically to distribute complexity. Applied to a hydrodynamic simulator, it outperforms traditional co-Kriging methods, achieving greater accuracy and efficiency. Recent improvements in the implementation of neural networks have also motivated research on multifidelity neural networks [14] to optimize engineering functions, particularly in scenarios requiring costly numerical simulations. For instance, in [15], the authors combine a fully connected neural network (NN) and a linear system for the interactions between codes to approximate scalar function and in inverse problems. The low-fidelity surrogate model can also be built using a fully connected NN, see [16] for a direct application, while the method in [17] provides options for single-fidelity active learning and is more general for multifidelity modeling.

To extend these approaches to simulations that produce time-dependent functions, particularly relevant to our work, adaptations of the formalism are required. In the study [18], the temporal nature of the code output is directly integrated into the structure of the neural network under consideration, using Long Short-Term Memory Recurrent Neural Networks (LSTM). To exploit several levels of fidelity, LSTMs for each level of fidelity are constructed in a sequential manner, by giving as input to each LSTM of higher fidelity the outputs of the LSTMs of lower fidelity. A common strategy for transitioning from single-output to multi-output predictions involves projecting the simulation output onto a finite set of carefully chosen deterministic functions. The relationships between the input parameters and each projection coefficient are then approximated, often independently. When these approximations leverage multiple models with varying fidelity, they form a MF model suited for multi-output simulations, as demonstrated in [19]. In this context, the choice of the projection basis plays a crucial role in the
accuracy of the prediction model. The projection functions must be well adapted to the highest fidelity code, as its output is the primary target for prediction. Simultaneously, to maximize the utility of lower-fidelity simulations, these functions should also maintain a reasonable level of compatibility with the lower-fidelity codes. For more details on these approaches combining projection and multifidelity, the interested reader can refer to [20]. In this work, various forms of MF meta-modeling and dimension reduction were applied to different Computational Fluid Dynamics (CFD) benchmark problems. However, the study is limited to cases where lower fidelity simulations are defined by coarser meshes only, to reduce the overall computational cost. This strategy was applied in the context of Gaussian Process Regression (GPR) in [21, 22, 23] working instead on tensorised forms of the covariance matrix of the Gaussian process under consideration. Additionally, Reduced Basis Methods (RBM) can be used to replace high fidelity models [24, 25, 26]. These methods project high-dimensional simulation outputs onto a lower-dimensional subspace, enabling efficient approximations while preserving key system dynamics. Efficient implementation of the traditional RBM is an intrusive procedure that requires manipulations of the HF code [27, 28]. In contrast, non-intrusive RBMs involve simulations on a coarser mesh [29, 30] or with a simplified code [31], allowing for faster computations while still capturing essential features.

In this work, we propose to study and compare in details three frameworks for MF meta-modeling with increasing complexity, which seem to be particularly interesting with respect to the existing literature for predicting the response of costly numerical codes with time-series output: 1) the Non-Intrusive Reduced Basis method (NIRB) 2) the Gaussian Process Regression (GPR), and 3) the Recurrent Neural Networks with Long Short-Term Memory (LSTM). These three approaches are presented in sufficient detail to make them accessible to non-specialists in various applications. The first approach is an extension of the two-grid NIRB method introduced in [32], and is chosen for its simplicity and for its ease of interpretation, as it only relies on an assumption of linearity between the outputs of different fidelity. The second approach is an extension of GPR [8] and [22], and also relies on a certain linearity between the different versions of the code of interest. But the GPR formalism allows for non-linear compensation in terms of projection coefficients. In this sense, it can be adapted to more complex code structures, but relies on a larger number of hyperparameters to be calibrated. Lastly, the LSTM we are considering, which is strongly inspired by the work of [18], makes no *a priori* assumption of linearity between the different versions of the code, and in this sense proposes a much more general formalism. However, this is achieved by introducing a much larger number of hyperparameters, the optimization of which being much more computationally demanding.

By focusing on these only three formalisms, we limit ourselves to non-intrusive MF modeling approaches that can easily be applied to existing physical models without requiring major modifications to the source codes. These approaches are expected not only to provide high accuracy, but also fast evaluation for real-time applications or for analysis requiring numerous simulations. The three families of methods are applied to a transient heat problem with a high dimensional input vector and three levels of fidelity. A comparative analysis of the different meta-models is made to address their respective performance in terms of accuracy and computational cost, but also to emphasize the varying degrees of ease with which these methods can be implemented and refined, without necessarily being specialists in them.

The articulation of this paper is thus divided into two main parts. Section 2 first presents the description of the
general problem as well as the theoretical bases of each approach used in the proposed study. The physical problem of interest is then presented in Section 3, and the numerical comparison of the three approaches on this case is presented.

## 2. Theoretical framework

### 2.1. Problem definition

Let $d$ and $N_t$ be two positive constants, and $\mathbb{X}$ be a connected subspace of $\mathbb{R}^d$. In addition, let $\mathcal{H}$ be the space of functions with values in $\mathbb{R}^{N_t}$ which are square integrable on $\mathbb{X}$, in the sense that for each $\mathbf{h}$ in $\mathcal{H}$,

$$ \|\mathbf{h}\|_{\mathcal{H}}^2 := \int_{\mathbb{X}} \|\mathbf{h}(\mathbf{x})\|^2 d\mathbf{x} < +\infty, \quad (1) $$

with $\|\cdot\|$ the Euclidean norm in $\mathbb{R}^{N_t}$. In this work, we are interested in constructing a predictor $\hat{\mathbf{y}}$ of a specific function $\mathbf{y} \in \mathcal{H}$ so that $\|\mathbf{y} - \hat{\mathbf{y}}\|_{\mathcal{H}}^2$ is as small as possible. We restrict ourselves here to the case where $\mathbf{y}(\mathbf{x})$ is a time series, i.e. $\mathbf{y}(\mathbf{x})$ can be seen as the projection onto a fixed grid of a function (a priori regular) indexed by a single index $t \in \mathbb{R}$ (typically the time). This amounts to assuming that there exist $N_t$ distinct instants $t_1 < \dots < t_{N_t}$ ($N_t$ being potentially large) and a real-valued function $\mathbf{y}$ defined on $\mathbb{X} \times \mathbb{R}$ so that:

$$ \mathbf{y}(\mathbf{x}) := (\mathbf{y}(\mathbf{x}; t_1), \mathbf{y}(\mathbf{x}; t_2), \dots, \mathbf{y}(\mathbf{x}; t_{N_t})). $$

The order of the components of $\mathbf{y}(\mathbf{x})$ thus makes sense, and it is reasonable to think that the closer $t$ and $t'$ are, the closer $\mathbf{y}(\mathbf{x}; t)$ and $\mathbf{y}(\mathbf{x}; t')$ will be. However, depending on the methods presented in the rest of this work, we will see that this assumption is not always required.

We also assume that each evaluation of $\mathbf{y}$ is computationally expensive, so that this construction of $\hat{\mathbf{y}}$ can only be based on a finite (and relatively small) number of calls to $\mathbf{y}$, which we denote by $n^*$. In this work, one of the central points is that we assume to have access to $s-1 \ge 1$ other functions of $\mathcal{H}$, which can be seen as sorted low-fidelity versions of $\mathbf{y}$, in the sense that the larger $1 \le \ell \le s-1$ is, the closer $\mathbf{y}_\ell$ is to $\mathbf{y}$:

$$ \|\mathbf{y} - \mathbf{y}_1\|_{\mathcal{H}}^2 \ge \|\mathbf{y} - \mathbf{y}_2\|_{\mathcal{H}}^2 \ge \dots \ge \|\mathbf{y} - \mathbf{y}_{s-1}\|_{\mathcal{H}}^2. \quad (2) $$

It is also assumed that the closer a code is to $\mathbf{y}$, the longer its numerical cost (execution time). In order to simplify the notations to follow, let us denote $\mathbf{y}_s := \mathbf{y}$, and for all $1 \le \ell \le s$, let us denote by $c_\ell$ the numerical cost of evaluating $\mathbf{y}_\ell$ (these costs are assumed to be known and to be independent of the value of $\mathbf{x}$). Thus, given these different codes, we can decide to construct the predictor $\hat{\mathbf{y}}$ of $\mathbf{y}$ not only on the basis of evaluations of $\mathbf{y}_\ell$ for a given $\ell$, but by coupling them to evaluations of $\mathbf{y}_1, \dots, \mathbf{y}_{\ell-1}$. And we can say that the learning set

$$ \mathcal{S}(n_1, \dots, n_s) := \{ \mathbf{y}_1(\mathbf{x}_1^{(1)}), \dots, \mathbf{y}_1(\mathbf{x}_1^{(n_1)}), \mathbf{y}_2(\mathbf{x}_2^{(1)}), \dots, \mathbf{y}_s(\mathbf{x}_s^{(n_s)}) \} $$

containing $n_\ell$ calls to $\mathbf{y}_\ell$ for all $1 \le \ell \le s$ is admissible for this construction (in terms of calculation time) if:

$$ \sum_{\ell=1}^{s} c_{\ell} \times n_{\ell} \le c_s \times n^* := c_{\text{tot}}. \quad (3) $$
By restricting ourselves to strategies with a total numerical cost lower than $c^{tot}$, several directions can be explored to propose the best (in the $L^2$ sense) predictor of $y$. For example, the following options could be considered:

* optimize the allocation of calculations according to the different codes, i.e. optimize the values of $n_1, \dots, n_s$,

* optimize the positions in $\mathbb{X}$ of all the input vectors to be evaluated for all the codes, i.e. optimize $S(n_1, \dots, n_s)$ for given values of $n_1, \dots, n_s$,

* optimize the choice of the class of functions in which the predictor will be looked for.

It is on this third option that the essential part of this work will focus, by studying the advantages and disadvantages of three classes of approximation of increasing complexity, namely the Non-intrusive Reduced Basis (NIRB) method, the Gaussian process regression (GPR), and the Long Short-Term Memory Recurrent Neural Networks (LSTM). In the following, we therefore assume that we have access to an admissible set $S(n_1, \dots, n_s)$ (with respect to $c^{tot}$) aggregating $n_\ell$ evaluations of the code $y_\ell$ for each $1 \le \ell \le s$, from which we aim to construct the best predictor of $y(x)$ at any $x$ of $\mathbb{X}$. In general, the evaluation points for these different codes can be chosen anywhere in $\mathbb{X}$. In this case, the better these points tend to fill the space $\mathbb{X}$ (for example using Latin hypercube sampling [33]), the more accurate we can expect the predictions of $y = y_s$ to be. However, we will see that working with nested designs of experiments, i.e. such that each point evaluated by $y_\ell$ is also part of the points evaluated by the code $y_{\ell-1}$ (implying $n_\ell \le n_{\ell-1}$) may often simplify the construction and improve the accuracy of the predictors. To simplify the reading, if $S(n_1, \dots, n_s)$ has this nested structure, we will impose a perfect match between the first $n_\ell$ points of each experimental design:

$$x_{\ell}^{(1)} = x_{\ell-1}^{(1)}, \dots, x_{\ell}^{(n_{\ell})} = x_{\ell-1}^{(n_{\ell})}, \quad 2 \le \ell \le s.$$

As indicated in the introduction, the aim of this work is not specifically to establish a hierarchy between these three methods, but

* to see to what extent the mathematical assumptions on which they are based can (more or less easily) be extended to functional outputs and different versions of $y$ are present and where $y(x)$,

* to highlight the main difficult points of such a construction (rather like a tutorial),

* and finally to discuss the conditions under which each method is likely to be particularly relevant.

*Remark.* The three methods described in the following, Non-Intrusive Reduced Basis (NIRB), Gaussian Process Regression (GPR), and Long Short-Term Memory networks (LSTM), are not combined in a hybrid fashion. Each is used independently as a surrogate modeling technique to predict high-fidelity outputs. In our study, the concept of multi-fidelity does not refer to a combination of different surrogate model types (e.g., NIRB + GPR), but rather to the use of different levels of physical model fidelity, that vary in both computational cost and physical complexity. The goal is to accurately predict specific quantities of interest (QoIs), from the high-fidelity model, not the full solution field. For each of the three surrogate models (NIRB, GPR, LSTM), we implement a multi-fidelity
learning strategy using simulations generated at different fidelity levels. These simulations serve as training data, allowing each method to learn a mapping from low- and medium-fidelity outputs to high-fidelity QoIs. Therefore, the multi-fidelity approach is applied at the output level, and focuses on efficiently approximating the high-fidelity model responses.

## 2.2. Non-intrusive Reduced Basis Methods

### 2.2.1. Definition of the approximation class

The first considered approach is the NIRB method [29, 34]. This method relies on the fact that for each $1 \le l \le s$, the manifold

$$ \mathcal{M}_l := \{ \mathbf{y}_l(\mathbf{x}) \in \mathbb{R}^{N_l} : \mathbf{x} \in \mathbb{X} \} \subset \mathcal{H} $$

can be well-approximated by a linear space spanned by a $m_l$-dimensional set of well-chosen functions. In its non-intrusive form and in the presence of several levels of code fidelity, its implementation is separated into two stages: the offline and the online stage. The offline phase assumes that the training set $S(n_1, \dots, n_s)$ has a nested structure (points evaluated by a higher-fidelity code are also evaluated by lower-fidelity codes), and is itself divided into two phases: an estimation phase for the different manifolds $\mathcal{M}_l$, followed by a rectification phase to maximise compatibility between two successive approximated manifolds. During the online phase, in order to predict the value of $\mathbf{y}_s$ in a new value of $\mathbf{x}$, the code with the lowest fidelity is evaluated, and $\tilde{\mathbf{y}}_s(\mathbf{x})$ is then deduced from $\mathbf{y}_1(\mathbf{x})$ using various projections and rectifications.

To go into more detail, for any $2 \le l \le s$, let $\mathbf{M}_l$ be the $(N_l \times n_l)$-dimensional matrix grouping together the different evaluations of $\mathbf{y}_l$:

$$ \mathbf{M}_l := [\mathbf{y}_l(\mathbf{x}_1) \cdots \mathbf{y}_l(\mathbf{x}_{n_l})], \qquad (4) $$

and $\mathbf{X}_l^{(m_l)}$ be the $(N_l \times m_l)$-dimensional matrix that gathers the $m_l$ eigenvectors of $\mathbf{M}_l \mathbf{M}_l^T$ associated with its $m_l$ highest eigenvalues (similar to what is obtained from Principal Component Analysis (PCA) [35]):

$$ \mathbf{M}_l \mathbf{M}_l^T = \mathbf{X}_l^{(N_l)} \mathbf{L}_l^{(N_l)} \mathbf{X}_l^{(N_l)^T}, \quad \mathbf{X}_l^{(N_l)} := [\mathbf{X}_l^{(m_l)} \mathbf{X}_l^{(N_l-m_l)}], \qquad (5) $$

$$ \mathbf{L}_l^{(N_l)} := \begin{bmatrix} L_{l,1} & 0 & \cdots & 0 \\ 0 & L_{l,2} & \cdots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & L_{l,N_l} \end{bmatrix}, \quad L_{l,1} \ge \dots \ge L_{l,N_l}, \quad \mathbf{X}_l^{(N_l)^T} \mathbf{X}_l^{(N_l)} = \mathbf{X}_l^{(N_l)} \mathbf{X}_l^{(N_l)^T} = \mathbf{I}_{N_l}, \qquad (6) $$

where $\mathbf{I}_{N_l}$ is the $(N_l \times N_l)$ identity matrix. For all $2 \le l \le s$, the manifold $\mathcal{M}_l$ is then approximated by:

$$ \tilde{\mathcal{M}}_l := \{ \mathbf{X}_l^{(m_l)} \mathbf{R}_l^{(m_l)} \mathbf{X}_l^{(m_l)^T} \mathbf{y}_{l-1}(\mathbf{x}) \in \mathbb{R}^{N_l} : \mathbf{x} \in \mathbb{X} \}. $$

The matrix $\mathbf{R}_l^{(m_l)}$ is called the rectification matrix and it is computed offline by solving the following linear system:
$$ \mathbf{X}_{\ell}^{(m_{\ell})T} \mathbf{M}_{\ell} = \mathbf{R}_{\ell}^{(m_{\ell})} \mathbf{X}_{\ell}^{(m_{\ell})T} \tilde{\mathbf{M}}_{\ell-1}, \qquad (7) $$

where $\tilde{\mathbf{M}}_{\ell-1}$ is the matrix gathering the $n_{\ell}$ first columns of $\mathbf{M}_{\ell-1}$ (so that the columns of $\mathbf{M}_{\ell}$ and $\tilde{\mathbf{M}}_{\ell-1}$ are associated with the same values of $\mathbf{x}$). As it is defined, the rectification process is essentially a simple linear regression between the higher and lower fidelity projection coefficients. In the particular case where $\mathbf{y}_{\ell-1}$ is exactly equal to $\mathbf{y}_{\ell}$ at any point of $\mathbb{X}$, so that $\mathbf{M}_{\ell} = \tilde{\mathbf{M}}_{\ell-1}$, let us notice that the rectification matrix is simply equal to the identity matrix, leading to the standard PCA-based approximation of $\mathbf{y}_{\ell}$:

$$ \tilde{\mathbf{M}}_{\ell} := \left\{ \mathbf{X}_{\ell}^{(m_{\ell})} \mathbf{X}_{\ell}^{(m_{\ell})T} \mathbf{y}_{\ell-1}(\mathbf{x}) \in \mathbb{R}^{N_{\ell}} : \mathbf{x} \in \mathbb{X} \right\} = \left\{ \mathbf{X}_{\ell}^{(m_{\ell})} \mathbf{X}_{\ell}^{(m_{\ell})T} \mathbf{y}_{\ell}(\mathbf{x}) \in \mathbb{R}^{N_{\ell}} : \mathbf{x} \in \mathbb{X} \right\}. $$

In practice, we recommend to regularize this rectification matrix to avoid over-fitting, this time by searching it in the form:

$$ \mathbf{R}_{\ell}^{(m_{\ell})} = \mathbf{X}_{\ell}^{(m_{\ell})T} \mathbf{M}_{\ell} \tilde{\mathbf{M}}_{\ell-1}^{T} \mathbf{X}_{\ell}^{(m_{\ell})} \left( \mathbf{X}_{\ell}^{(m_{\ell})T} \tilde{\mathbf{M}}_{\ell-1} \tilde{\mathbf{M}}_{\ell-1}^{T} \mathbf{X}_{\ell}^{(m_{\ell})} + \lambda_{\ell} \mathbf{I}_{m_{\ell}} \right)^{-1}, \qquad (8) $$

with $\lambda_{\ell} > 0$ a relaxation term, whose value needs to be optimized (see section 2.2.2 for more details). Putting all these projections/rectifications together, the multi-fidelity NIRB predictor of $\mathbf{y}(\mathbf{x})$ in any $\mathbf{x} \in \mathbb{X}$ is finally written:

$$ \tilde{\mathbf{y}}^{\text{NIRB}}(\mathbf{x}) := \Pi^{\text{NIRB}} \mathbf{y}_1(\mathbf{x}), \quad \Pi^{\text{NIRB}} := \prod_{\ell=2}^{s} \mathbf{X}_{\ell}^{(m_{\ell})} \mathbf{R}_{\ell}^{(m_{\ell})} \mathbf{X}_{\ell}^{(m_{\ell})T}. \qquad (9) $$

*Remark.* Since we only need the $m_{\ell}$ first eigenvectors of $\mathbf{M}_{\ell} \mathbf{M}_{\ell}^T$, note that it is generally more efficient numerically, when $N_t \gg n_{\ell}$, to compute $\mathbf{X}_{\ell}^{(m_{\ell})}$ directly from a singular value decomposition (SVD) of $\mathbf{M}_{\ell}$, which does not require the manipulation and diagonalization of a matrix of size $(N_t \times N_t)$ [35].

### 2.2.2. Specific features of the method

By relying exclusively on explicit expressions, the NIRB approach is particularly simple to implement. It should also be noted that this method only depends on a small number of hyperparameters, which also makes it particularly attractive: the values of the truncation parameters $\{m_{\ell}\}_{\ell=2}^s$, and potentially the relaxation terms $\{\lambda_{\ell}\}_{\ell=2}^s$. Of course, the efficiency of the method depends heavily on the fact that the manifolds $\mathcal{M}_{\ell}$ can be approximated by reduced-dimensional linear spaces, and on the fact that once the eigenspaces have been estimated, the dependence of $\mathbf{y} = \mathbf{y}_{\ell}$ on the new value of $\mathbf{x}$ can be carried exclusively by $\mathbf{y}_1(\mathbf{x})$ (matrix $\Pi^{\text{NIRB}}$ does not depend on $\mathbf{x}$). Without surprise, the robustness of the method can therefore be challenged if the biases between the low and the high fidelity models are parameter dependent. In this case, the risk of overfitting is central, hence the importance of the relaxation term in this step.

In addition, the choice of the hyperparameters $m_{\ell}$ and $\lambda_{\ell}$ is generally made by cross-validation: the set $\mathcal{S}(n_1, \dots, n_s)$ is decomposed into a training set and a test set, and the values of the relaxation and truncation parameters that give the best approximations are searched. While the parameters associated with two pairs of codes can be estimated independently, it is important to note that, for a given pair of codes $(\ell, \ell-1)$, the values of $\lambda_{\ell}$ and $m_{\ell}$ have to be chosen in the same procedure. Indeed, the larger $m_{\ell}$ is, the more likely it is that the linear system associated with the estimation of $\mathbf{R}_{\ell}^{(m_{\ell})}$ is likely to be ill-conditioned.
It should also be noted that the NIRB approach was presented in a framework where $S(n_1, \dots, n_s)$ has a nested structure. But it can be generalized to cases where this is not true (the results being often much worse in this case). In this case, calculating the rectification matrices as proposed no longer makes sense, but we can propose to define them such that

$$ \mathbf{R}_{\ell}^{(m_{\ell})} = r_{\ell}^2 \mathbf{I}_{m_{\ell}}, \quad r_{\ell}^2 := \frac{n_{\ell-1}}{n_{\ell}} \frac{\sum_{i=1}^{n_{\ell}} \|\mathbf{y}_{\ell}(\mathbf{x}_{\ell}^{(i)})\|^2}{\sum_{i=1}^{n_{\ell-1}} \|\mathbf{y}_{\ell-1}(\mathbf{x}_{\ell-1}^{(i)})\|^2} \quad (10) $$

in order to ensure at least a reasonable reconstructed $L^2$ norm of $\mathbf{y}_s$.

In terms of computation costs, since the predictor is expressed in terms of the lowest-fidelity code, we can see that the offline cost, on the one hand, breaks down into four groups:

* the $n_{\ell}$ evaluations of each code $\mathbf{y}_{\ell}$,

* the $s-1$ SVDs for the definition of the projection matrices $\mathbf{X}_{\ell}^{(m_{\ell})}$ whose complexities are in $O(N_{\ell}m_{\ell}^2)$ (note that in this approach, we do not need to compute the reduced basis for $M_1$),

* the inversion of $s-1$ linear systems for the definition of the rectification matrices $\mathbf{R}_{\ell}^{(m_{\ell})}$, whose complexities are in $O(m_{\ell}^3)$,

* the aggregation of the global projection matrix $\Pi_{\ell}^{\text{NIRB}}$.

Concerning the online cost, it is limited to the cost of evaluating the lowest-fidelity code $\mathbf{y}_1$, then a matrix vector product. Finally, if the cost $c_1$ of evaluating $\mathbf{y}_1$ is non-negligible, it is important to note that this online cost can be problematic in configurations where the predictor has to be evaluated a very large number of times (when solving optimization or inversion problems, for example).

## 2.3. Gaussian process regression

### 2.3.1. Definition of the approximation class

The second method we are considering is Gaussian process regression (GPR). This method is based on a relatively different paradigm. The functions $\mathbf{y}_1, \dots, \mathbf{y}_s$ are viewed as particular realizations of $s$ dependent Gaussian processes $\mathbf{Y}_1, \dots, \mathbf{Y}_s$. Under this hypothesis, on the basis of the evaluations contained in $S(n_1, \dots, n_s)$, the best predictor (in the $L^2$ sense) of $\mathbf{y}_s(\mathbf{x})$ in a new value $\mathbf{x}$ of $\mathcal{X}$ is then given by:

$$ \hat{\mathbf{y}}(\mathbf{x}) := \mathbb{E}[\mathbf{Y}_s(\mathbf{x}) | S(n_1, \dots, n_s)], \quad (11) $$

where $\mathbb{E}[\cdot]$ is the mathematical expectation, and the symbol $|$ indicates the statistical conditioning. The precise construction of $\hat{\mathbf{y}}$ is therefore based on the knowledge of the statistical dependencies between two evaluations of the same code at two different values of $\mathbf{x}$, but also between the evaluations of two codes with different fidelity.

In its most standard form (which is also the one that best exploits its strengths), the multifidelity GPR method assumes a particular nested and affine relationship between all the code levels [13], in the sense that for each $\mathbf{x} \in \mathcal{X}$,
$$ Y_{\ell}(\mathbf{x}) = \begin{cases} W_1(\mathbf{x}) & \text{if } \ell = 1, \\ B_{\ell}(\mathbf{x})Y_{\ell-1}(\mathbf{x}) + W_{\ell}(\mathbf{x}) & \text{otherwise,} \end{cases} \qquad (12) $$

where $\mathbf{x} \mapsto B_2(\mathbf{x}), \dots, \mathbf{x} \mapsto B_s(\mathbf{x})$ are $s-1$ deterministic matrix-valued functions (including constant functions to model constant correlations between codes of close fidelity) as used in [13], and where $\mathbf{x} \mapsto W_1(\mathbf{x}), \dots, \mathbf{x} \mapsto W_s(\mathbf{x})$ are $s$ statistically independent Gaussian processes with values in $\mathbb{R}^{N_t}$, whose mean and covariance functions are denoted by $\mu_\ell$ and $C_\ell$ so that:

$$ \mu_{\ell}(\mathbf{x}) := \mathbb{E}[W_{\ell}(\mathbf{x})], \quad C_{\ell}(\mathbf{x}, \mathbf{x}') := \mathrm{Cov}(W_{\ell}(\mathbf{x}), W_{\ell}(\mathbf{x}')), \quad \mathbf{x}, \mathbf{x}' \in \mathbb{X}, \ 1 \le \ell \le s. \qquad (13) $$

By construction, the matrices $B_\ell$ characterize the correlations between two codes of close fidelity, so that the larger $B_\ell(\mathbf{x})$ is, the more information the observation of $Y_{\ell-1}(\mathbf{x})$ provides about $Y_\ell(\mathbf{x})$. As a consequence, they play a similar role than the rectification matrices introduced in Section 2.2, in the sense that the closer $Y_{\ell-1}(\mathbf{x})$ is to $Y_\ell(\mathbf{x})$, the more they tend towards the identity matrix. But in this GPR formalism, the compensation role is also performed by $W_\ell$. And unlike the NIRB method, which assumes the rectification matrices constant, both $B_\ell$ and $W_\ell$ a priori depend on $\mathbf{x}$, for greater modeling flexibility. In this formalism, it is important to notice that the evaluation $C_\ell(\mathbf{x}, \mathbf{x}')$ is a $(N_t \times N_t)$ matrix and therefore characterizes a double source of dependence between: the correlation in space associated with its dependence on $\mathbf{x}$ and $\mathbf{x}'$, and the correlation in time when looking at its $N_t \times (N_t - 1)$ off-diagonal terms for fixed values of $\mathbf{x}$ and $\mathbf{x}'$ (for more explanations see [36] and [21]).

Coupling Gaussian processes with the affine modeling presented in Eq. (12) is particularly convenient, as it allows us to derive an explicit expression for $\mathbb{E}[Y_s(\mathbf{x}) | S(n_1, \dots, n_s)]$ (this would not be the case if non-linear transforms or non-Gaussian processes were considered). Indeed, assuming that the functions $B_\ell$, $\mu_\ell$ and $C_\ell$ (presented explicitly in the following) are known for all $\ell$, if we write

$$ \bar{\mathbf{y}} := (\mathbf{y}_1(\mathbf{x}_1^{(1)}), \dots, \mathbf{y}_1(\mathbf{x}_1^{(n_1)}), \mathbf{y}_2(\mathbf{x}_2^{(1)}), \dots, \mathbf{y}_2(\mathbf{x}_2^{(n_2)}), \dots, \mathbf{y}_s(\mathbf{x}_s^{(1)}), \dots, \mathbf{y}_s(\mathbf{x}_s^{(n_s)})), \qquad (14) $$

$$ \bar{\mathbf{Y}} := (\mathbf{Y}_1(\mathbf{x}_1^{(1)}), \dots, \mathbf{Y}_1(\mathbf{x}_1^{(n_1)}), \mathbf{Y}_2(\mathbf{x}_2^{(1)}), \dots, \mathbf{Y}_2(\mathbf{x}_2^{(n_2)}), \dots, \mathbf{Y}_s(\mathbf{x}_s^{(1)}), \dots, \mathbf{Y}_s(\mathbf{x}_s^{(n_s)})), \qquad (15) $$

we obtain an explicit form (the detailed expressions of $\mathbf{m}(\mathbf{x}), \bar{\mathbf{m}}, \mathbf{V}(\mathbf{x}), \mathbf{U}(\mathbf{x})$ and $\bar{\mathbf{V}}$ are provided in Appendix A) for the joint distribution of $(Y_s(\mathbf{x}), \bar{Y})$,

$$ (Y_s(\mathbf{x}), \bar{Y}) \sim \mathcal{N}\left(\left(\begin{array}{c} \mathbf{m}(\mathbf{x}) \\ \bar{\mathbf{m}} \end{array}\right), \left[\begin{array}{cc} \mathbf{V}(\mathbf{x}) & \mathbf{U}(\mathbf{x}) \\ \mathbf{U}^{\mathrm{T}}(\mathbf{x}) & \bar{\mathbf{V}} \end{array}\right]\right), \qquad (16) $$

and the meta-model is obtained by Gaussian conditioning, :

$$ \begin{aligned} \mathbb{E}[Y_s(\mathbf{x}) | S(n_1, \dots, n_s)] &= \mathbb{E}[Y_s(\mathbf{x}) | \bar{Y} = \bar{\mathbf{y}}] \\ &= \mathbf{m}(\mathbf{x}) + \mathbf{U}(\mathbf{x})\bar{\mathbf{V}}^{-1}(\bar{\mathbf{y}} - \bar{\mathbf{m}}) :=: \hat{\mathbf{y}}(\mathbf{x}), \quad \mathbf{x} \in \mathbb{X}. \end{aligned} \qquad (17) $$
Furthermore, the conditional covariance, which the Gaussian conditioning formula provides in the same way as it does for the conditional expectation,

$$ \mathrm{Cov}(\mathbf{Y}_s(\mathbf{x}), \mathbf{Y}_s(\mathbf{x}) | \bar{\mathbf{Y}} = \bar{\mathbf{y}}) = \mathbf{V}(\mathbf{x}) - \mathbf{U}(\mathbf{x})\bar{\mathbf{V}}^{-1}\mathbf{U}(\mathbf{x})^T, \quad (18) $$

makes it possible to quantify the uncertainty (and therefore, by contrast, the accuracy) of the prediction of $\mathbf{y}_s$ in $\mathbf{x}$, which is one of the great advantages of the GPR formalism: the higher this covariance, the less reliable the prediction.

Although the problem of constructing the predictor of $\mathbf{y}_s$ has now been solved from a theoretical point of view, it is still far from easy to construct from a numerical and practical point of view. Indeed, looking at Eq. (2.3.1), we notice that its definition relies on the inversion of the matrix $\bar{\mathbf{V}}$, whose dimension is $(N_t \sum_{\ell=1}^s n_\ell) \times N_t \sum_{\ell=1}^s n_\ell)$, which clearly becomes extremely challenging when $N_t$ increases, or when we have a large number of evaluations of all the codes. Several options were then proposed to reduce this numerical cost and construct this predictor in a more pragmatic way. To begin, instead of trying to condition the $s$ numerical codes by all the code evaluations contained in $S(n_1, \dots, n_s)$, we can try to proceed sequentially, as Eq. (12) invites us to do (see Appendix B for the detailed expressions of the obtained mean and covariance functions). The fact that the processes $\mathbf{W}_1, \dots, \mathbf{W}_s$ are statistically independent assures us that this sequential approach is equivalent to the global approach. By proceeding in this way, instead of inverting a matrix of size $(N_t \sum_{\ell=1}^s n_\ell) \times N_t \sum_{\ell=1}^s n_\ell)$, we now have to invert $s$ matrices of size $(N_t n_1 \times N_t n_1), \dots, (N_t n_s \times N_t n_s)$ [13].

But if $N_t$ is large, this cost may still be too high. In order to greatly reduce these matrix sizes, similar to what is proposed in the NIRB approach, we can then propose to no longer focus on the prediction of the complete functions, but to limit ourselves to predicting their projection into a smaller space which can reduce the risk of over-fitting, which is crucial for models such as Gaussian processes [37]. To do this, for $1 \le m \ll N_t$, let $\mathbf{X}_m^* = \mathbf{I}_m$. Each numerical code can then be projected into the space generated by the columns of $\mathbf{X}_m^*$, and for all $1 \le \ell \le s$ and all $\mathbf{x} \in \mathbb{X}$, we can denote these projection coefficients by

$$ z_{\ell}(\mathbf{x}) := \mathbf{X}_{m}^{*T} y_{\ell}(\mathbf{x}). \quad (19) $$

We can then reinterpret all of the above by no longer directly modeling the functions $y_\ell$ by Gaussian processes, but by now assuming that it is these projection coefficients $z_1, \dots, z_s$ which can be seen as particular realizations of the Gaussian processes $\mathbf{Y}_1, \dots, \mathbf{Y}_s$. To simplify reading, note that in the following, we use the same notations for the Gaussian processes associated with the reduced vectors $z_1, \dots, z_s$ than the ones initially associated with the full output vectors $\mathbf{y}_1, \dots, \mathbf{y}_s$. As a consequence, these different Gaussian processes, whose dimensions are now all equal to $m$ and no longer to $N_t$, can again be written in terms of processes $\mathbf{W}_1, \dots, \mathbf{W}_s$ and matrices $\mathbf{B}_1, \dots, \mathbf{B}_s$ as shown in Eq. (12), but the dimensions $N_t$ of these matrices and vectors would this time be also reduced to $m$. Using matrix $\mathbf{X}_m^*$ to project all the evaluations of $\mathbf{y}_1, \dots, \mathbf{y}_s$, we obtain a new training set, which can be used to (sequentially) condition all these Gaussian processes $\mathbf{Y}_1, \dots, \mathbf{Y}_s$. If we now denote by $\mathbf{Y}_s^{(n_1, \dots, n_s)}$ the Gaussian process $\mathbf{Y}_s$ (used here to model $z_s$) once conditioned by the projections onto the columns of $\mathbf{X}_m^*$ of all
the evaluations of $S(n_1, \dots, n_s)$, the final predictor of $\mathbf{y} = \mathbf{y}_s$ is then given by:

$$ \hat{\mathbf{y}}^{\text{GPR}}(\mathbf{x}) := \mathbf{X}_{\text{m}}^{*} \mathbb{E}[\mathbf{Y}_{s}^{(n_1, \dots, n_s)}(\mathbf{x})], \quad \mathbf{x} \in \mathbb{X}. \qquad (20) $$

The matrix $\mathbf{X}_{\text{m}}^{*} \text{Cov}(\mathbf{Y}_{s}^{(n_1, \dots, n_s)}(\mathbf{x}), \mathbf{Y}_{s}^{(n_1, \dots, n_s)}(\mathbf{x})) \mathbf{X}_{\text{m}}^{*T}$ can again be used to quantify the accuracy of this predictor in $\mathbf{x}$. But if we wanted to be more complete, since this uncertainty only concerns what happens in the space generated by the columns of $\mathbf{X}_{\text{m}}^{*}$, we would need to add to this uncertainty an approximation of the projection error, characterized by $\|\mathbf{y}_s - \mathbf{X}_{\text{m}}^{*} \mathbf{X}_{\text{m}}^{*T} \mathbf{y}_s\|_{H}^2$.

The numerical cost associated with the construction of this predictor is now reduced to the inversion of $s$ matrices of respective sizes $(mn_1 \times mn_1), \dots, (mn_s \times mn_s)$, which is much more reasonable. Assuming, moreover, that the components of projection vectors $\mathbf{z}_1, \dots, \mathbf{z}_s$ are uncorrelated (and therefore independent in this Gaussian framework), i.e. assuming that the matrices $C_1, \dots, C_s$ of the processes $\mathbf{W}_1, \dots, \mathbf{W}_s$ are diagonal, we can treat the subspaces generated by the $m$ columns of $\mathbf{X}_{\text{m}}^{*}$ independently, and the numerical cost now boils down to the inversion, $m$ times (one time for each column of $\mathbf{X}_{\text{m}}^{*}$), of $s$ matrices of sizes $(n_1 \times n_1), \dots, (n_s \times n_s)$ respectively.

For the sake of completeness, we must now return to the choice of $\mathbf{X}_{\text{m}}^{*}$, and to the definition of $\mu_1, \dots, \mu_s$, $\mathbf{B}_1, \dots, \mathbf{B}_s$ and $C_1, \dots, C_s$, which characterize the statistical moments of the $s$ Gaussian processes $\mathbf{Y}_1, \dots, \mathbf{Y}_s$ used for the prediction. In general, we are quite free with these definitions, in the sense that the only conditions concern the covariance functions, which must be positive definite functions. Nevertheless, in a small data context, where a reduced number of evaluations to the different codes are available, it is not very reasonable to make these different quantities depend on too many unknown hyperparameters (otherwise there is a great risk of overfitting). Regarding the choice of $\mathbf{X}_{\text{m}}^{*}$, it is customary to limit oneself to the first vectors of the PCA basis of $\mathbf{y}_s$, which would amount to taking $\mathbf{X}_{\text{m}}^{*} = \mathbf{X}_{s}^{(m)}$ by reusing the notations of Section 2.2, the value of $m$ being chosen in order to restore a certain percentage of the observed variance of $\mathbf{y}_s$. Note nevertheless that contrary to the NIRB approach, this basis is not used to approximate only the manifold associated with $\mathbf{y}_s$, but also to project the evaluations of $\mathbf{y}_{s-1}, \dots, \mathbf{y}_1$. Also, looking for a projection basis that is adapted to $\mathbf{y}_s$, but also to the lower fidelity codes makes perfect sense here. And for this, we can consider the adaptation of the PCA algorithm that is proposed in Appendix C and based on the study [38]. Then, one can choose the different covariance functions from standard parametric classes. In our case, we can for example choose, for all $1 \le \ell \le s$, to write function $C_{\ell}$ (in the diagonal case associated with the independence assumption for the projection coefficients) so that

$$ C_{\ell}(\mathbf{x}, \mathbf{x}') = \text{diag}(\sigma_{1,\ell}^2 k (\|\text{diag}(\mathbf{l}^{1,\ell})^{-1}(\mathbf{x} - \mathbf{x}')\|), \dots, \sigma_{m,\ell}^2 k (\|\text{diag}(\mathbf{l}^{m,\ell})^{-1}(\mathbf{x} - \mathbf{x}')\|)), \qquad (21) $$

where $h \mapsto k(h) = (1 + \sqrt{5}h + \frac{5}{3}h^2) \exp(-\sqrt{5}h)$ is the Matérn-5/2 function (see [39] for alternative choices), and where for any vector $\mathbf{v}$, $\text{diag}(\mathbf{v})$ defines the diagonal matrix constructed from the components of $\mathbf{v}$. Hence, each covariance matrix $C_{\ell}$ is characterized by $(d+1) \times m$ independent hyperparameters associated with $m$ $d$-dimensional correlation vectors $\mathbf{l}^{1,\ell}, \dots, \mathbf{l}^{m,\ell}$ and $m$ variance terms $\sigma_{1,\ell}^2, \dots, \sigma_{m,\ell}^2$. All that remains is to look at the average functions, which are characterized by two quantities: the matrices $\mathbf{B}_{\ell}$ and the vectors $\mu_{\ell}$. With no real a priori on these functions, we could, for example, write their components in the form of weighted sums of polynomials in
x. In this case, the number of hyperparameters to be estimated for these functions corresponds to the number of weights considered. Note that if we restrict ourselves to polynomials of order 1 with a constant term, i.e. if:

$$ \mu_{\ell}(\mathbf{x}) = \left( a_0^{(1)} + \sum_{i=1}^{d} a_i^{(1)} x_i, \dots, a_0^{(m)} + \sum_{i=1}^{d} a_i^{(m)} x_i \right), \quad \mathbf{B}_{\ell}(\mathbf{x}) = \mathrm{diag} \left( b_0^{(1)} + \sum_{i=1}^{d} b_i^{(1)} x_i, \dots, b_0^{(m)} + \sum_{i=1}^{d} b_i^{(m)} x_i \right), \quad (22) $$

the mean function of each Gaussian process $Y_{\ell}$ is characterized by $2m(d + 1)$ independent hyperparameters. We can then try to estimate all these hyperparameters by cross-validation or likelihood maximization. At first sight, as the different Gaussian processes are by construction nested within each other, it may seem necessary to estimate in the same procedure the hyperparameters of all the Gaussian processes. However, if the training set $S(n_1, \dots, n_s)$ has a nested structure, then we can again decompose the global problem in $s$ smaller problems, and estimate the hyperparameters of each Gaussian process $Y_{\ell}$ independently. And if we have assumed diagonal structures for $C_{\ell}$ and $B_{\ell}$, as is the case in Eqs. (21) and (22), then we can even estimate the hyperparameters associated with each component of $Y_{\ell}$ in a completely parallel way, which amounts to solving $m \times s$ independent (and very standard) optimization problems in (small) dimension $3(d + 1)$ (see [39, 40, 13] for more details).

### 2.3.2. Specific features of the method

Unlike the NIRB approach presented in Section 2.2, the GPR approach introduces a relatively large number of parameters: parameters related to the *a priori* mean and covariance functions of the processes $W_1, \dots, W_{\ell}$, as well as parameters for the correlations between codes with close fidelity. This parameterisation makes the model more flexible, enabling it to adapt to a variety of configurations. However, it also makes the model more difficult to learn and manipulate, and requires a greater number of optimization problems to be solved. Although it is possible to propose relatively standard forms for all these parametric shapes, depending on the case studies encountered, it is often possible to propose more suitable forms. But optimizing these parametric shapes while avoiding overfitting remains a difficult problem, requiring a certain amount of experience in handling these GPR models.

It is also interesting to note that the predictor proposed by GPR has an explicit form, and ultimately relies only on manipulations of matrices and vectors, while the hyperparameters on which it relies all have a relatively clear meaning. These characteristics make the results provided by these GPR models relatively interpretable and expli-cable, which is a definite advantage among the set of approximation classes that could be found in the literature. But one of the great strengths of these GPR approaches really lies in the fact that, in addition to the predictor, they provide a natural characterization of the prediction uncertainty (by post-processing the various covariance matrices), which is rare enough among machine learning models to be worth mentioning.

In terms of computational costs, as with the NIRB approach, a distinction can be made between offline and online costs. The offline costs include defining the projection matrix $X_m^*$, solving the optimization problems associated with estimating the hyperparameters, and then inverting the different covariance matrices associated with the code evaluations. If we restrict ourselves to the case where the matrices $B_{\ell}$ and $C_{\ell}$ are diagonal, it should be remembered that there are $m \times s$ problems in dimension $3(d + 1)$ that need to be solved for the hyperparameters estimation, then $s$ matrices (and this must be done $m$ times) of size $(n_1 \times n_1), \dots, (n_s \times n_s)$ that need to be inverted for the construction of the predictor.
Once the optimization problems are solved, and the inverted matrices have been saved, the online costs are limited, as for the NIRB approach, to matrix-vector products (potentially in much greater number than for the NIRB approach), which allows to ensure very reduced prediction times.

## 2.4. Recurrent Neural Networks (RNN)

### 2.4.1. Long Short-Term Memory (LSTM)

The third type of model we consider is based on recurrent neural networks (RNNs), specifically Long Short-Term Memory (LSTM) networks. This choice represents another fundamental shift in paradigm compared to the NIRB and GPR methods. Unlike these more classical approaches, which rely on explicit statistical assumptions, LSTM exploits the power of deep learning to extract and model temporal relationships within data, and has proven particularly effective for predicting complex time series, especially in contexts where the relationships between past and future observations are neither linear nor trivial.

RNNs were developed to process sequential data, where each new observation potentially depends on previous ones, which well corresponds to the context of time series in which we placed ourselves in Section 2.1. Unlike standard neural networks, which treat each input independently, a RNN maintains a memory of past states through recurrent connections. For this internal memory, we generally speak of input and hidden layers associated with hyperbolic tangent activation functions, as well as a dense output layer associated with a linear activation function for final prediction.

This memory allows the model to capture temporal dependencies, but also introduces a major challenge: when training a RNN on long sequences, the gradients used to update the model's weights tend to either vanish or explode as they propagate backward through time. This phenomenon, known as the vanishing and exploding gradient problem (discussed in [41]), makes it difficult for RNNs to learn long-term dependencies, as the contributions of past information fade or become unstable over many time steps. To overcome this limitation, Hochreiter and Schmidhuber introduced Long Short-Term Memory (LSTM) networks [42], which incorporate a specialized mechanism to regulate how information is stored and updated over time. Instead of relying solely on a single hidden state, LSTM introduces a structured way of controlling the flow of information using gates. These gates act as dynamic filters that decide which parts of the input should be remembered, which should be forgotten, and how much of the past information should influence the current state.

In the standard single-fidelity setting, an LSTM model is trained on a set of input sequences, learning to map past data to future predictions. At each time step $1 \le t \le N_t$, the model takes as input the same vector $\mathbf{x}$, and updates its internal memory before producing a prediction for $(\mathbf{y}(\mathbf{x}))_{1\le t\le N_t}$. This classical LSTM framework can easily be extended to a multifidelity setting by incorporating information from lower-fidelity models as an additional input at each time step. Instead of relying solely on high-fidelity observations, the LSTM also receives evaluations from the lower-cost, lower-accuracy codes. This approach allows the network to exploit the structured knowledge of the low-fidelity model, effectively guiding the learning process and reducing the amount of high-fidelity data required for accurate predictions. Mathematically, the input to the network at each time step $1 \le t \le N_t$ consists of both the input vector $\mathbf{x}$ describing the system, and the corresponding low-fidelity values $(\mathbf{y}_1)_{t}, \dots, (\mathbf{y}_{t-1})_{t}$. The advantage
of this approach is twofold (see [18]). First, it reduces the reliance on expensive high-fidelity simulations, since the model can achieve high accuracy even with a limited number of high-fidelity training points. Second, the low-fidelity codes provide informative priors that can help guide the network's training. This may lead to better convergence and improved generalization. However, the introduction of multifidelity data also adds complexity to the training process. Indeed, the network has to properly weigh the contributions of the different low-fidelity codes while avoiding biases that could degrade accuracy. Ensuring that the model correctly balances the different sources of information requires careful hyperparameter tuning, particularly in terms of network architecture, learning rate, and batch size selection [43].

In order to be able to predict $y(x)$ without having to evaluate one of the available codes at this value of $x$, the previous LSTM formalism needs to be slightly adapted again. And in the same way as for the GPR approach, we propose to proceed sequentially:

* **Initialization:** for $\ell = 1$, we build a first simple-fidelity LSTM model, noted $x \mapsto \hat{y}_1^{\text{LSTM}}(x)$ by only exploiting the $n_1$ evaluations of $y_1$ that are in $S(n_1, \dots, n_\ell)$.

* **Iteration:** For any $2 \le \ell \le s$, we then construct a two-fidelity LSTM model of $y_\ell$ using for the training phase the $n_\ell$ evaluations of $y_\ell$ contained in $S(n_1, \dots, n_\ell)$. And we make this model take as input at each time step $t$ the value of $x$ but also the prediction at $t$ of the $\ell - 1$ LSTM model, so that:

$$ \hat{y}_{\ell}^{\text{LSTM}}(\mathbf{x}) = \mathcal{F}_{\ell}^{\text{LSTM}}(\mathbf{x}; \hat{y}_{\ell-1}^{\text{LSTM}}(\mathbf{x})), \qquad (23) $$

where $\mathcal{F}_{\ell}^{\text{LSTM}}$ represents the learned mapping.

The final prediction model proposed by the LSTM approach, noted $\hat{y}^{\text{LSTM}}$, finally exploits all the evaluations of $S(n_1, \dots, n_\ell)$, and is written in the following nested form:

$$ \hat{y}^{\text{LSTM}}(\mathbf{x}) = \mathcal{F}_{s}^{\text{LSTM}}(\mathbf{x}; \mathcal{F}_{s-1}^{\text{LSTM}}(\mathbf{x}; \mathcal{F}_{s-2}^{\text{LSTM}}(\mathbf{x}; \dots \mathcal{F}_{2}^{\text{LSTM}}(\mathbf{x}; \hat{y}_{1}^{\text{LSTM}}(\mathbf{x})))\dots)), \quad \mathbf{x} \in \mathbb{X}. \qquad (24) $$

A scheme of the LSTM neural network used to predict a high fidelity temporal simulation is presented in Figure 1. In the LSTM layers, with hidden states $h_t$ at each time step $t$, the activation function that is generally used is the hyperbolic tangent function (tanh) producing values in the range $[-1, 1]$. This activation function is well-suited for LSTM networks because it helps to introduce non-linearities and control the values of the hidden states. In the dense output layer, a linear activation function is classically used. For more details about the LSTM model, reader can refer to [41, 44].

### 2.4.2. Specific features of the method

LSTM models have been introduced to capture complex, non-linear dependencies between time steps, making them a powerful alternative to traditional statistical models, especially in contexts where the governing equations are unknown [45] or too complex to be modelled explicitly [46]. Unlike methods such as NIRB and GPR, which rely on explicit mathematical formulations and provide interpretable uncertainty quantification, LSTM operates
<figure><img src="image_2.png" /><figcaption>Figure 1: Scheme of an LSTM network for predicting $(\mathcal{G}_t^{\text{LSTM}})_t$ using $x$ and $(\mathcal{G}_{t-1}^{\text{LSTM}})_t$. The hidden states $h_t$ at each time step $t$ are supposed to capture information about the current input $(x, (\mathcal{G}_{t-1}^{\text{LSTM}})_t)$ as well as previous states $h_{t-1}$ using several hidden layers and a dense output layer.</figcaption></figure>

as a black-box model, making its internal structure difficult to interpret. The weights and activation functions of an LSTM network evolve in a high-dimensional space, making it challenging to understand which features of the input data contribute most to a given prediction (see [47]). This lack of interpretability contrasts with GPR, where the covariance function explicitly describes the correlation structure, and NIRB, which decomposes the solution onto predefined reduced basis functions.

Before starting with LSTM models, it is important to notice that they are extremely non-linear, such that small changes in hyperparameters, such as the number of layers, the number of hidden units, or the choice of activation functions, can have drastic effects on performance. Unlike NIRB, where the dimensionality reduction step provides direct insight into the dominant modes of variation, or GPR, where hyperparameters can often be optimized systematically, tuning an LSTM requires extensive experimentation and empirical validation. Beyond hyperparameter sensitivity, LSTMs are also highly dependent on regularization techniques to prevent overfitting. Since deep neural networks tend to memorize training data rather than generalizing, techniques such as dropout, batch normalization, and $L^2$ regularization are necessary to ensure robust performance. However, these techniques introduce additional complexity. Reader can refer to the study [48] exploring the different architectures and applications of LSTM networks for time series prediction, presenting several variants.

Another fundamental challenge when using LSTMs is their offline computational cost, which is significantly higher than that of NIRB or GPR. This cost is particularly demanding as it involves multiple passes through the dataset with gradient-based optimization, requiring extensive GPU acceleration and substantial memory resources. Unlike GPR, which only requires solving series of relatively small linear systems, or NIRB, where the offline phase consists mainly of computing reduced basis functions and a few low-dimensional projections, training an LSTM involves back-propagation through time, which requires storing and updating long sequences of hidden states. This results in high memory consumption and long computation times, especially when dealing with long time series or large datasets. For neural networks and LSTMs in particular, the dependence of the training dynamics to the choice of optimizer, batch size, sequence length, network depth, or hardware architecture makes therefore difficult the definition of a precise computational complexity in the classical sense for the offline phase, and theoretical guarantees on convergence rates are limited.

The online cost, however, is significantly lower: once trained, the LSTM can make predictions in real time with relatively little computation, as it only needs to perform forward propagation through the network. Nevertheless,
despite the efficiency of LSTMs at inference time, their practical deployment is still constrained by hardware requirements, as deep learning models are optimized for GPU architectures.

Despite these challenges, LSTMs offer a powerful solution for multifidelity time series prediction. By combining the flexibility of deep learning with the structured knowledge of low-fidelity models, they enable accurate and computationally efficient forecasting, reducing the need for extensive high-fidelity data.

Despite the underlying complexity of these models, modern deep learning frameworks such as TensorFlow (with Keras) and PyTorch have made LSTMs widely accessible, even to non-specialists. These libraries provide high-level abstractions that allow users to define and train LSTM networks with only a few lines of code, without needing to manually implement back-propagation, gradient updates, or memory management (see for instance Appendix D for an example of code that will be used for numerical experiments in Section 3.3). Nevertheless, users must remain aware of the potential pitfalls, such as overfitting, sensitivity to hyperparameters, and the need for sufficient training data, as these factors can significantly impact predictive accuracy and reliability.

# 3. Application to a multi-layers thermal problem

## 3.1. *Description of the physical problem*

To assess the thermal and energy performance of buildings, it is necessary to develop non-destructive in situ testing and effective identification methods to determine the thermal properties of opaque walls, e.g. thermal resistance and heat capacity. To achieve this purpose, a strategy integrating physical models and sensor outputs was proposed in [22, 49, 50] to provide an estimate of wall thermal properties in a short time, with a limited number of sensors and while being less sensitive to outdoor weather conditions. In this method, an active thermal excitation is applied on a restricted part of the internal surface of the wall and the measured dynamic response of the wall is used for the identification of wall thermal properties by solving an inverse problem. As shown in [22], Bayesian inference can provide operational credibility intervals for the identified thermal parameters, but it requires a large number of direct model simulations resulting in high computational costs, even when considering tractable physical models. Therefore, metamodeling and multifidelity approaches are essential.

In this Section, we thus focus on the dynamic thermal response of an internal insulation wall (IIW) represented in Figure 2. Such a wall is characterized for each layer i by the thickness ($e_i$), conductivity ($\lambda_i$) and volumetric heat capacity ($\rho C_i$). Four time-varying outputs, also called quantities of interest, are considered. These time series corresponds to the temperature of the internal (resp. external) surface, noted $T_{S_i}$ (resp. $T_{S_e}$) and to the flux of the internal (resp. external) surface, noted $\phi_{S_i}$ (resp. $\phi_{S_e}$). These quantities are notably also impacted by the internal and external transfer coefficient ($h_{int}$ and $h_{ext}$). As can be seen in red in Figure 2, a restricted part of the internal surface of the wall is excited with Heaviside-type heat flux excitation ($\phi_{exc}$). In fact, this excitation allows controlled profiles during the time ($t$) of $T_{S_i}(t)$ and $\phi_{S_i}(t)$ whereas fluctuations over time between day and night occur for $T_{S_e}(t)$ and $\phi_{S_e}(t)$ due to natural weather variations.

The three approaches proposed in Section 2 are tested and compared in the following on four time-varying outputs from the considered thermal problem. This thermal problem is particularly interesting because the studied outputs are of different kinds (temperature and heat flux) and display a variety of different behaviors. On internal surface,
the controlled outputs profiles are influenced by 3D thermal effects which are induced by the localized thermal excitation on the internal surface. While outdoor weather conditions cause uniform excitation on the external surface and more complex time-fluctuating profiles for outputs associated with this surface. Hence, the considered thermal problem will be appropriate for studying the effectiveness and limitations of NIRB, GPR and LSTM approaches when using combinations of different model fidelity levels for predicting various time-varying outputs.

<figure><img src="image_3.png" /><figcaption>Figure 2: Geometry of a 4-layer IIW considered in the study displaying the different parameters of the physical model: internal transfer coefficient $h_{\text{int}}$, external transfer coefficient $h_{\text{ext}}$, thickness $e_i$, conductivity $\lambda_i$, Volumetric heat capacity $\rho C_i$ for $i \in [1, 4]$. Also, the quantities of interest ($T_{S_i}$, $T_{S_e}$, $\varphi_{S_i}$, $\varphi_{S_e}$), the solar ($\Phi_{\text{sun}}$) flux and the excitation ($\Phi_{\text{exc}}$) flux are presented. $z$ is the depth coordinate and $r$ is the radius from the center of the heating zone.</figcaption></figure>

To be more precise, the code we are focusing on is characterized by

* fourteen input parameters, noted $h_{\text{int}}$, $h_{\text{ext}}$, $e_i$, $\rho C_i$, $\lambda_i$, for $i \in \{1, 2, 3, 4\}$ and gathered in **x**, which characterize the wall geometry, its material properties, and external conditions properties. The range of variation of these fourteen inputs is given in Appendix F;

* four time-varying outputs associated with the discretization of $T_{S_i}$, $T_{S_e}$, $\varphi_{S_i}$, $\varphi_{S_e}$ on a grid of time steps spaced by $\Delta t = 2$min, over a total duration of $t_f=72$h (leading to $N_t = 2161$). These different outputs will be processed independently and each corresponds to the vector **y** introduced in Section 2.

In this numerical application, we have access to $s = 3$ codes of varying fidelity to compute the four outputs of interest (the equations on which these models are based are detailed in Appendix E):

1. $\mathbf{y}_3$ is the high fidelity (noted HF in the following) response, and corresponds to the output provided by a two-dimensional axisymmetric Finite Element (FE) thermal model with 30271 degrees of freedom (dof). It is equivalent to a three dimensional thermal model. The time required for each evaluation of $\mathbf{y}_3$ is approximately equal to $c_3 = C_{\text{HF}} = 2.3$ min.

2. $\mathbf{y}_2$ is the medium fidelity response (referred to as $LF_{1D}$), and corresponds to the response of a one-dimensional FE thermal model with 401 dof. Each evaluation of this code lasts around $c_2 = C_{LF_{1D}} = 0.2$ min.

3. $\mathbf{y}_1$ is the low fidelity response (referred to as $LF_{RC}$), and corresponds to the $R_5C_4$ thermal model with five resistors and four capacitances with 4 dof and a computational cost of $c_1 = C_{LF_{RC}} \approx 0.01$ sec per simulation. As $C_{LF_{RC}} \ll C_{LF_{1D}} < C_{\text{HF}}$, we propose to neglect this cost for the sake of simplicity in the following.
It should be noted that these calculation times were obtained using a laptop equipped with an Intel Core i7-1185G7 (11th generation) processor with a maximum frequency of 4.8 GHz and 32 GB of RAM memory.

A comparison of the different outputs obtained from these three models is presented in Figure 3, for a particular value of $x$. For a better understanding of the numerical results to come, it is important to notice that the excitation imposed on the interior surface allows better controlled profiles of $T_{S_i}$ and $\phi_{S_i}$, but this comes at the cost of introducing 3D physical phenomena that only the HF model will be able to capture. We therefore expect non-negligible and potentially model parameter-dependent differences on this surface between the HF model and the models $LF_{1D}$ and $LF_{RC}$. On the contrary, on the outside surface (for $T_{S_e}$ and $\phi_{S_e}$), as there is no local excitation, the entire surface is mostly excited by the solar flux, and we will see that the physics is nearly 1D (the 3D phenomena will arrive after a certain time, but to a lesser degree, via internal excitation once the four wall layers have been crossed), and we expect relatively similar behavior between the codes HF and $LF_{1D}$.

<figure><img src="image_4.png" /><figcaption>Figure 3: Comparison of the different levels of fidelity considered in the numerical application, for four time-series simulations: (a) interior surface temperature $T_{S_i}$, (b) exterior surface temperature $T_{S_e}$, (c) interior surface flux $\phi_{S_i}$ and (d) exterior surface flux $\phi_{S_e}$. The HF is presented in red with empty triangles pointing downwards, the $LF_{1D}$ in green with filled triangles pointing upwards and the $LF_{RC}$ in blue with dots.</figcaption></figure>

In other hand, it is important to study the impact of the input parameters first on the HF model, and then on the relation between the different fidelity levels. To a better understanding of how the outputs evolve with respect to variations in the input parameters, we present in the follow:

1. a visualization of time-dependent surface fluxes, generated using the HF model for three distinct parameter sets (presented in Table 1). These parameter sets were selected to highlight the impact of the input parameters used on simulations generated using the HF thermal model.

2. a visualization of time-dependent surface flux using the HF, $LF_{1D}$, and $LF_{RC}$ models, for two distinct parameter sets (presented in Table 2). These parameter sets were selected to emphasize different relationships between fidelity levels: one that maximizes the gap between the HF model and the $LF_{1D}$, and one that minimizes it.
<table><caption>Table 1: Details of the three input parameters used to generate the high fidelity simulations. Where $h_{int}$ (resp. $h_{ext}$) the internal (resp. external) heat transfer coefficient, $e_i$ the thickness, $\lambda_i$ the thermal conductivity and $\rho C_i$ the volumetric heat capacity of the layer i.</caption><thead><tr><th></th><th>h<sub>int</sub></th><th>h<sub>ext</sub></th><th>e<sub>1</sub></th><th>e<sub>2</sub></th><th>e<sub>3</sub></th><th>e<sub>4</sub></th><th>λ<sub>1</sub></th><th>λ<sub>2</sub></th><th>λ<sub>3</sub></th><th>λ<sub>4</sub></th><th>ρC<sub>4</sub></th><th>ρC<sub>2</sub></th><th>ρC<sub>3</sub></th><th>ρC<sub>4</sub></th></tr></thead><tbody><tr><th>P<sub>1</sub></th><td>6.58</td><td>5.38</td><td>0.01</td><td>0.18</td><td>0.19</td><td>0.01</td><td>0.51</td><td>0.02</td><td>0.197</td><td>1.77</td><td>9 × 10<sup>5</sup></td><td>5.9 × 10<sup>4</sup></td><td>1.5 × 10<sup>6</sup></td><td>7.6 × 10<sup>5</sup></td></tr><tr><th>P<sub>2</sub></th><td>5.34</td><td>36</td><td>0.02</td><td>0.05</td><td>0.16</td><td>0.02</td><td>0.51</td><td>0.05</td><td>2.05</td><td>1.37</td><td>1 × 10<sup>6</sup></td><td>5.4 × 10<sup>4</sup></td><td>1.4 × 10<sup>6</sup></td><td>8.1 × 10<sup>5</sup></td></tr><tr><th>P<sub>3</sub></th><td>8.79</td><td>5.55</td><td>0.01</td><td>0.19</td><td>0.26</td><td>0.01</td><td>0.68</td><td>0.03</td><td>1.02</td><td>0.92</td><td>1.3 × 10<sup>6</sup></td><td>6.8 × 10<sup>4</sup></td><td>1.8 × 10<sup>6</sup></td><td>5 × 10<sup>5</sup></td></tr></tbody></table>

<table><caption>Table 2: Details of the two input parameters used to generate the simulations for models with different fidelity levels. Where $h_{int}$ (resp. $h_{ext}$) the internal (resp. external) heat transfer coefficient, $e_i$ the thickness, $\lambda_i$ the thermal conductivity and $\rho C_i$ the volumetric heat capacity of the layer i.</caption><thead><tr><th></th><th>h<sub>int</sub></th><th>h<sub>ext</sub></th><th>e<sub>1</sub></th><th>e<sub>2</sub></th><th>e<sub>3</sub></th><th>e<sub>4</sub></th><th>λ<sub>1</sub></th><th>λ<sub>2</sub></th><th>λ<sub>3</sub></th><th>λ<sub>4</sub></th><th>ρC<sub>4</sub></th><th>ρC<sub>2</sub></th><th>ρC<sub>3</sub></th><th>ρC<sub>4</sub></th></tr></thead><tbody><tr><th>P<sub>1</sub></th><td>5.13</td><td>7.71</td><td>0.02</td><td>0.17</td><td>0.23</td><td>0.02</td><td>0.77</td><td>0.05</td><td>1.05</td><td>0.51</td><td>8.1 × 10<sup>5</sup></td><td>1.2 × 10<sup>5</sup></td><td>1.9 × 10<sup>6</sup></td><td>1.6 × 10<sup>6</sup></td></tr><tr><th>P<sub>2</sub></th><td>9.06</td><td>30</td><td>0.01</td><td>0.08</td><td>0.18</td><td>0.01</td><td>0.47</td><td>0.02</td><td>1.2</td><td>0.97</td><td>6.9 × 10<sup>5</sup></td><td>1.4 × 10<sup>5</sup></td><td>2.1 × 10<sup>6</sup></td><td>1.2 × 10<sup>6</sup></td></tr></tbody></table>

First, Figure 4 compare three simulations of the interior (resp. exterior) surface flux $\phi_{S_i}$ (resp. $\phi_{S_e}$). These simulations are generated using the HF model, for the three input parameters presented in Table 1. The comparison shows that the behavior of the HF model varies significantly depending on the parameter sets used.

<figure><img src="image_5.png" /><figcaption>Figure 4: Comparison of (a) interior surface flux $\phi_{S_i}$ and (b) the exterior surface flux $\phi_{S_e}$ generated using the high fidelity model for three parameter sets ($P_1 \in [1, 3]$).</figcaption></figure>

Then, Figure 5 shows the exterior surface flux under solar excitation. We observe that:

* The different models respond consistently to changes in the input parameters, although the magnitude and timing of thermal peaks vary.

* The degree of mismatch between fidelity levels depends on the specific parameter configuration, demonstrating the need for models that can adapt to local behavior in the parameter space.

* Despite discrepancies in amplitude or transient behavior, the models remain strongly correlated in shape and dynamics, which validates the use of multi-fidelity modeling techniques.
<figure><img src="image_6.png" /><figcaption>Figure 5: Comparison of exterior surface flux $\phi_{S_e}$ generated using the different fidelity model for the same two parameter sets : (a) case showing larger mismatch between fidelity levels and (b) case with near-perfect agreement between all models. The simulations generated using the HF model are shown in red, those generated using the LF<sub>1D</sub> in green and in blue those generated using the LF<sub>RC</sub>.</figcaption></figure>

## 3.2. Description of the comparisons made

In order to simplify the reading of the results, in this section we present only the results associated with the interior and exterior surfaces flux ($\phi_{S_i}$, $\phi_{S_e}$), and place the results for the interior and exterior surfaces temperatures ($T_{S_i}$, $T_{S_e}$) in Appendix G, (the conclusions on $\phi_{S_i}$ and $\phi_{S_e}$ being relatively close to those on $T_{S_i}$ and $T_{S_e}$).

The different approaches described in Section 2 are compared in terms of accuracy using $n_{Test} = 300$ input vectors gathered in $S^{Test} = \{x_j^{Test}, 1 \le j \le n_{Test}\}$, which have been uniformly chosen in $\mathbb{X}$ and which have not been used during the model training. Three kinds of performance indicators are calculated to quantify the quality of the constructed predictors, which all rely on the post-processing of the mean values over time of the absolute error between the HF response and its predictor,

$$ \delta_i = \frac{1}{N_t} \sum_{j=1}^{N_t} |y_{i,j} - \hat{y}_{i,j}|, \qquad (25) $$

where for each $1 \le i \le n_{Test}$, to simplify notation, we note $y_{i,j}$ and $\hat{y}_{i,j}$ the jth component of $y(x_j^{Test})$ and $\hat{y}(x_j^{Test})$ respectively. The first indicator is noted $MAE^* = \frac{1}{n_{Test}} \sum_{i=1}^{n_{Test}} \delta_i$, and corresponds to the average error over the test set. Then, we introduce $Q_{95\%}$ as the empirical 95% inter-quantile on $\delta_i$. The closer the values of these two indicators are to zero, the better the predictor. At last, we consider the determination coefficient $R^{2*}$ (which must be close to 1 for a good prediction) so that:

$$ R^{2*} = 1 - SS_{res}/SS_{tot}, \qquad (26) $$
$$SS_{\text{res}} = \frac{1}{n_{\text{Test}}} \sum_{i=1}^{n_{\text{Test}}} \frac{1}{N_t} \sum_{j=1}^{N_t} (y_{i,j} - \hat{y}_{i,j})^2, \quad SS_{\text{tot}} = \frac{1}{n_{\text{Test}}} \sum_{i=1}^{n_{\text{Test}}} \frac{1}{N_t} \sum_{j=1}^{N_t} (y_{i,j} - \bar{y}_i)^2, \quad \bar{y}_i = \frac{1}{N_t} \sum_{j=1}^{N_t} y_{i,j}.$$

We now propose to compare three types of predictors of decreasing quality by varying the number and the type of low-fidelity codes that are taken into account.

1. The first class of predictors, which is noted $M_{2a}^{HF}$, relies on **two** levels of fidelity and is based on evaluations of the models HF and $LF_{1D}$ only. It requires a call to $LF_{1D}$ for each prediction. For this reason, it is the most expensive predictor during the online phase (from a numerical point of view), but is expected to be most accurate because it is based on the true evaluation of $LF_{1D}$, and not on its approximation based on a finite number of evaluations only. The results associated with this configuration are compared in Table 3 for $\phi_{S_i}$ and $\phi_{S_c}$, and in Appendix G.8 for $T_{S_i}$ and $T_{S_c}$.

2. The second class of predictors, which is noted $M_3^{HF}$, relies this time on **three** levels of fidelity, and is based on evaluations of the models HF, $LF_{1D}$ and $LF_{RC}$. We remind that for this class, the predictors are constructed in two steps: first, a predictor of $LF_{1D}$ noted $M^{LF}$, is constructed from evaluations of $LF_{1D}$ and $LF_{RC}$ using $n_{LF}$ simulations, which correspond to the number of simulations of the $LF_{1D}$ model, requiring a call to the $LF_{RC}$ for each evaluation. The HF predictor is then constructed using the $n_{HF}$ additional evaluations of HF and requires an evaluation of $M^{LF}$ for each prediction. It is expected to be a little less accurate than $M_{2a}^{HF}$, since an error is propagated from $M^{LF}$ to the final predictor $M_3^{HF}$. In this case, two combinations of simulations are considered presented by $M_{(3,1)}^{HF}$ and $M_{(3,2)}^{HF}$ described in the following. The corresponding results are compared in Tables 4 and 5 for $\phi_{S_i}$ and $\phi_{S_c}$ and in Appendix G.9 and G.10 for $T_{S_i}$ and $T_{S_c}$.

3. The third class of predictors, which is noted $M_{2b}^{HF}$, relies again on **two** levels of fidelity, but combines this time evaluations of models HF and $LF_{RC}$ (no evaluation of the 1D code). In terms of online cost, it only requires for each prediction a call to the $LF_{RC}$, which is the least costly. As it is based on $LF_{RC}$ and not $LF_{1D}$, it is expected to be less accurate compared to the other predictors. The associated results are compared in Table 6 for $\phi_{S_i}$ and $\phi_{S_c}$, and in Appendix G.11 for $T_{S_i}$ and $T_{S_c}$.

In this study, the same training sets are used for the construction of the three predictors. In each Table, the performance of the predictors are compared for $c_{tot} \in \{100, 500, 1000\} \times C_{HF}$. But as mentioned previously and studied numerically in [20], there is no standard rule to choose the number of simulations for each fidelity level while respecting the cost constraint. In the following, for $M_{2a}^{HF}$ and $M_{2b}^{HF}$, i.e. for predictors based on only two levels of fidelity, we limit ourselves to the case where the same number of evaluations is carried out for each code. In the case of the 3-level predictors, $M_3^{HF}$, two combinations are considered: we note $M_{(3,1)}^{HF}$ the configuration when $n_{HF} = n_{LF}$ (see Tables 4 and G.9 for the results), and we note $M_{(3,2)}^{HF}$ the configuration when $n_{HF} < n_{LF}$. In that case, inspired by the study proposed in [20], we take ($n_{HF}C_{HF} = 0.9 \times c_{tot}$, $n_{LF}C_{LF_{1D}} = 0.1 \times c_{tot}$) if $c_{tot} < 500 \times C_{HF}$, and ($n_{HF}C_{HF} = 0.75 \times c_{tot}$, $n_{LF}C_{LF_{1D}} = 0.25 \times c_{tot}$) otherwise (see Tables 5 and G.10 for the results).

Regarding the choice of hyper-parameters, we would like to stress once again that there is no absolute rule for the choice of the parametric classes associated with the GPR and LSTM methods. It is therefore perfectly legitimate to try out different classes and different structures, in a more or less empirical way, or with the help of (numerically
low-cost if possible) sensitivity analyses [2, 3]. So as not to make the conclusions of future comparisons too dependent on these choices, we have limited ourselves in this work to relatively standard forms.

1. For the NIRB method, for each $1 \le l \le s$, the size of the projection bases is chosen equal to $m_l = 50$, while the relaxation term is taken equal to $\lambda_l = 10^{-1}$.

2. For the GPR method:

* A weighted PCA is applied including the three fidelity levels studied to construct the projection basis. The truncation parameters $m$ depend on the configurations and are chosen such that the mean absolute error associated with the outputs of each code is lower than $10^{-3}$ for each quantity of interest (see Appendix C for an illustration of the algorithm used).

* Then, the matrix $B_l(\mathbf{x})$, for each fidelity level $l \in \{2, 3\}$, is modeled as a linear function with respect to the 14 input parameters $\mathbf{x}$: $B_l(\mathbf{x}) = \text{diag}(\beta_{l,0}^{(1)} + \sum_{i=1}^{14} \beta_{l,i}^{(1)} x_i, \dots, \beta_{l,0}^{(m)} + \sum_{i=1}^{14} \beta_{l,i}^{(m)} x_i)$.

* Regarding the statistical properties of the Gaussian processes $W_1, \dots, W_s$, the mean functions $\mu_l$ are modeled as combinations of linear terms of all the inputs parameters, while the covariance functions $C_l$ belong to the Matérn-5/2 kernel class.

In our case, all the unknown hyperparameters are estimated by their maximum likelihood estimators.

3. For the LSTM method (see Appendix D for an example of LSTM code written in python), the following options were chosen:

* the LSTM network is made up of 4 stacked hidden layers with 512 neurons in each layer, and a fully connected dense layer with a linear activation function.

* the dropout rate is chosen equal to 0.3 (used after each LSTM layer for regularization), and an Adam optimizer is used with a learning rate of $10^{-4}$ to minimise the mean square error as loss function.

### 3.3. Analysis of the results

Starting with Table 3, we first see that when using the codes HF and $LF_{1D}$ for the prediction of the inner surface flux ($\phi_{S_i}$), the NIRB method seems to be the less relevant. Indeed, despite a decrease in the lower bound of the $Q_{95\%}$ indicator when the size of the training set increases, we see that the upper bound remains more or less the same, while still being far greater than those of the other two methods. This is not surprising, as the presence of a local excitation on the interior surface introduces parameter-dependent 3D effects, which, despite the rectification matrix, are more difficult to take into account with the NIRB approach.

The LSTM method produces more accurate predictions than those obtained by the NIRB approach, especially when the size of the learning base becomes large, which is a standard specificity of the RNN methods. For this configuration $M_{2a}^{HF}$, it is the GPR method that performs best, which can be explained by the fact that the link between the HF and $LF_{1D}$ codes is relatively linear once adjusted by the matrix $B(\mathbf{x})$ (which may depend on the model parameters $\mathbf{x}$).

The results are significantly different when we now focus on the exterior surface. Indeed, as the phenomena are mainly 1D on this surface, smaller differences between the models HF and $LF_{1D}$ are observed. For that, and
contrary to the results associated with $\phi_{S_i}$, the NIRB approach provides the most accurate prediction of $\phi_{S_e}$ when the training set is very small. However, for larger training sets, the GPR method surpass the NIRB method and once again provides the best prediction.

<table><caption>Table 3: Prediction quality of $\phi_{S_i}$ and $\phi_{S_e}$ in W/m² from the three predictors using the one-dimensional thermal model (LF₁D) as low fidelity model, based on the indicators MAE*, Q95% and R2* calculated over a test set of size 300. c<sub>tot</sub> aggregates the costs associated with the HF and LF₁D evaluations.</caption><thead><tr><th rowspan="2">M<sup>HF</sup><sub>2a</sub></th><th colspan="3">φ<sub>S<sub>i</sub></sub></th><th colspan="3">φ<sub>S<sub>e</sub></sub></th><th rowspan="2">n<sub>HF</sub> = n<sub>LF</sub></th><th rowspan="2">c<sub>tot</sub></th></tr><tr><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th></tr></thead><tbody><tr><td>NIRB</td><td>1.291</td><td>[0.101; 4.640]</td><td>0.996</td><td><b>1.029</b></td><td>[0.091;2.756]</td><td>0.999</td><td rowspan="3">92</td><td rowspan="3">c<sub>tot</sub> = 100C<sub>HF</sub></td></tr><tr><td>GPR</td><td><b>0.822</b></td><td>[0.140;2.918]</td><td>0.998</td><td>1.261</td><td>[0.433;3.039]</td><td>0.999</td></tr><tr><td>LSTM</td><td>1.564</td><td>[1.073;5.443]</td><td>0.992</td><td>3.387</td><td>[1.828;6.208]</td><td>0.995</td></tr><tr><td>NIRB</td><td>1.059</td><td>[0.070;4.199]</td><td>0.997</td><td>0.931</td><td>[0.062;3.011]</td><td>0.999</td><td rowspan="3">460</td><td rowspan="3">c<sub>tot</sub> = 500C<sub>HF</sub></td></tr><tr><td>GPR</td><td><b>0.470</b></td><td>[0.063;1.664]</td><td>0.999</td><td><b>0.579</b></td><td>[0.205;1.466]</td><td>0.999</td></tr><tr><td>LSTM</td><td>0.868</td><td>[0.523;3.392]</td><td>0.997</td><td>1.689</td><td>[0.834;4.144]</td><td>0.998</td></tr><tr><td>NIRB</td><td>1.036</td><td>[0.052;4.136]</td><td>0.997</td><td>0.901</td><td>[0.057;2.796]</td><td>0.999</td><td rowspan="3">920</td><td rowspan="3">c<sub>tot</sub> = 1000C<sub>HF</sub></td></tr><tr><td>GPR</td><td><b>0.299</b></td><td>[0.046;1.253]</td><td>0.999</td><td><b>0.571</b></td><td>[0.160;1.570]</td><td>0.999</td></tr><tr><td>LSTM</td><td>0.556</td><td>[0.340;2.670]</td><td>0.998</td><td>1.207</td><td>[0.614;2.914]</td><td>0.999</td></tr></tbody></table>

A better compromise between online computational cost and accuracy can be achieved using the 3-levels meta-model $M_3^{HF}$. For this predictor, we remind that the model LF₁D is also replaced by a 2-level predictor relying on the available evaluations of LF₁D and LFRC. In this case, the online computational cost of $M_3^{HF}$ is very small, as it is more or less equal to the computational cost of LFRC. The prediction capacity of these predictors is shown in Table 4 when considering $n_{HF} = n_{LF}$, and in Table 5 when $n_{LF} > n_{HF}$.

If we start by looking at what is happening on the interior surface, i.e. the side where the differences between HF and LF₁D are the greatest, the first observation we can make is that the predictions provided by $M_3^{HF}$ are only slightly less accurate when compared to those obtained by $M_{2a}^{HF}$. However for the prediction of $\phi_{S_e}$, there is a significant deterioration, which can once again be explained by the fact that the physics is almost 1D on the exterior surface, and as such, being able to use the real LF₁D outputs is particularly interesting, while the predictions associated with $M_3^{HF}$ are affected by errors due to the replacement of LF₁D by its predictor. It should also be noted that the NIRB approach rapidly converges, which is related to its rectification process. But it is not the case for the two others approaches (GPR and LSTM) as the prediction accuracy still improves when adding more points in the training set.

By comparing the Tables 4 and 5, we also confirm the difficulties of finding an optimal configuration for the distribution of calls to the different codes. The choice depends on the quantity studied and the approach used. However, the results are close for the different tests carried out with $M_{(3,1)}^{HF}$ and $M_{(3,2)}^{HF}$ for a similar total cost, which is reassuring in practice.

The results obtained using only calls to the models HF and LFRC are presented in Table 6. Favoring calls to LFRC rather than LF₁D can be explained by the fact that for certain applications, the online numerical cost corresponding
<table><caption>Table 4: Prediction quality of $\phi_{S_i}$ and $\phi_{S_e}$ in W/m² from the three predictors using the R<sub>5</sub>C<sub>4</sub> model (LF<sub>RC</sub>) and the one-dimensional thermal model (LF<sub>1D</sub>) as LF models, based on the indicators MAE*, Q<sub>95%</sub> and R<sup>2*</sup> calculated over a test set of size 300. c<sub>tot</sub> aggregates the costs associated with the HF and LF<sub>1D</sub> evaluations.</caption><thead><tr><th rowspan="2">M<sup>HF</sup><sub>(3,1)</sub></th><th colspan="3">φ<sub>S<sub>i</sub></sub></th><th colspan="3">φ<sub>S<sub>e</sub></sub></th><th rowspan="2">n<sub>HF</sub> = n<sub>LF</sub></th><th rowspan="2">c<sub>tot</sub></th></tr><tr><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th></tr></thead><tbody><tr><td>NIRB</td><td>1.498</td><td>[0.226;4.211]</td><td>0.995</td><td><b>2.191</b></td><td>[0.578; 5.354]</td><td>0.998</td><td rowspan="3">92</td><td rowspan="3">c<sub>tot</sub> = 100C<sub>HF</sub></td></tr><tr><td>GPR</td><td><b>0.992</b></td><td>[0.255;3.041]</td><td>0.997</td><td>3.200</td><td>[1.038;6.721]</td><td>0.995</td></tr><tr><td>LSTM</td><td>1.789</td><td>[1.287;6.197]</td><td>0.989</td><td>3.731</td><td>[1.987;7.106]</td><td>0.994</td></tr><tr><td>NIRB</td><td>1.344</td><td>[0.173;4.505]</td><td>0.995</td><td>1.979</td><td>[0.453;5.612]</td><td>0.998</td><td rowspan="3">460</td><td rowspan="3">c<sub>tot</sub> = 500C<sub>HF</sub></td></tr><tr><td>GPR</td><td><b>0.636</b></td><td>[0.191;1.864]</td><td>0.998</td><td>2.157</td><td>[0.596;6.444]</td><td>0.997</td></tr><tr><td>LSTM</td><td>0.874</td><td>[0.584;4.245]</td><td>0.997</td><td><b>1.577</b></td><td>[0.840;2.886]</td><td>0.999</td></tr><tr><td>NIRB</td><td>1.295</td><td>[0.206;4.394]</td><td>0.995</td><td>1.934</td><td>[0.479; 5.516]</td><td>0.998</td><td rowspan="3">920</td><td rowspan="3">c<sub>tot</sub> = 1000C<sub>HF</sub></td></tr><tr><td>GPR</td><td><b>0.451</b></td><td>[0.129;1.279]</td><td>0.999</td><td>1.616</td><td>[0.388;5.515]</td><td>0.998</td></tr><tr><td>LSTM</td><td>0.687</td><td>[0.502;2.624]</td><td>0.998</td><td><b>1.090</b></td><td>[0.562;2.015]</td><td>0.999</td></tr></tbody></table>

624

625

626

627

628

629

630

631

632

633

634

635

636

637

638

639

640

641

642

643

644

645

to a call to LF<sub>1D</sub> may be too high (when solving inverse problems, for example). In this table, as expected, we first note that the predictions with M<sup>HF</sup><sub>2b</sub> are less accurate than M<sup>HF</sup><sub>2a</sub> for both GPR and NIRB methods but relatively close for the LSTM method (the models in both configurations being still highly relevant in terms of R<sup>2*</sup> indicators). In terms of error trends, the results in this table are relatively similar to those observed in Table 3, both in terms of the differences between interior and exterior, and in terms of the compared performance between NIRB, GPR and LSTM. This configuration also allows us to quantify the benefits of switching to a three-fidelity predictor rather than a 2-fidelity predictor, while maintaining the same online costs as a call to LF<sub>RC</sub>.

However for the prediction of $\phi_{S_e}$ the accuracy is not improved when using M<sup>HF</sup><sub>(3,1)</sub> compared to M<sup>HF</sup><sub>2b</sub>. As for the accuracy of the prediction when using M<sup>HF</sup><sub>(3,2)</sub> with the GPR method, it is significantly improved compared to the M<sup>HF</sup><sub>2b</sub> only when c<sub>tot</sub> is the largest. On the other hand, when applying the NIRB and LSTM approaches, similar results are obtained using the three predictors M<sup>HF</sup><sub>2b</sub>, M<sup>HF</sup><sub>(3,1)</sub> and M<sup>HF</sup><sub>(3,2)</sub>. It means that, for the same online cost, the prediction quality can be improved using the 3-levels predictors M<sup>HF</sup><sub>(3,1)</sub> and/or M<sup>HF</sup><sub>(3,2)</sub>.

Concerning the quality of predictions of the interior surfaces temperatures (T<sub>S<sub>i</sub></sub>) presented in Appendix G, the conclusions are similar to the one made for $\phi_{S_i}$, presented above. While, for the exterior surfaces temperatures (T<sub>S<sub>e</sub></sub>) the NIRB approach gives the most precise predictions for all c<sub>tot</sub>. We note that for T<sub>S<sub>i</sub></sub> and T<sub>S<sub>e</sub></sub>, the three tested approaches give a high accuracy with a R<sup>2*</sup> > 0.99 for the different c<sub>tot</sub>.

## 3.4. General remarks

After analyzing the results obtained with the different approaches (NIRB, GPR, LSTM) for the different configurations tested (M<sup>HF</sup><sub>2a</sub>, M<sup>HF</sup><sub>3</sub>, M<sup>HF</sup><sub>2b</sub>), it appears that, for the considered application, the best compromise between online computational cost and accuracy is provided by the 3-level predictor. Indeed, in the majority of cases, this predictor achieves performance close to that obtained using the true 1D code, but at a much lower cost.

We can also draw some general conclusions from this application.
<table><caption>Table 5: Prediction quality of $\phi_{S_i}$ and $\phi_{S_e}$ in W/m² from the three predictors using the R<sub>5</sub>C<sub>4</sub> model (LF<sub>RC</sub>) and the one-dimensional thermal model (LF<sub>1D</sub>) as LF models, based on the indicators MAE*, Q<sub>95%</sub> and R<sup>2*</sup> calculated over a test set of size 300. c<sub>tot</sub> aggregates the costs associated with the HF and LF<sub>1D</sub> evaluations. In this case we choose n<sub>LF</sub> > n<sub>HF</sub></caption><thead><tr><th rowspan="2">M<sup>HF</sup><sub>(3,2)</sub></th><th colspan="3">φ<sub>S<sub>i</sub></sub></th><th colspan="3">φ<sub>S<sub>e</sub></sub></th><th rowspan="2">n<sub>HF</sub></th><th rowspan="2">n<sub>LF</sub></th><th rowspan="2">c<sub>tot</sub></th></tr><tr><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th></tr></thead><tbody><tr><td>NIRB</td><td>1.562</td><td>[0.224;4.754]</td><td>0.994</td><td><strong>2.282</strong></td><td>[0.631;5.378]</td><td>0.998</td><td rowspan="3">90</td><td rowspan="3">115</td><td rowspan="3">c<sub>tot</sub> = 100C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>1.230</strong></td><td>[0.316;4.414]</td><td>0.995</td><td>3.362</td><td>[1.039;8.114]</td><td>0.995</td></tr><tr><td>LSTM</td><td>1.737</td><td>[1.283;4.818]</td><td>0.991</td><td>3.458</td><td>[1.973;6.755]</td><td>0.995</td></tr><tr><td>NIRB</td><td>1.374</td><td>[0.223;4.472]</td><td>0.995</td><td>1.999</td><td>[0.515;5.600]</td><td>0.998</td><td rowspan="3">375</td><td rowspan="3">1437</td><td rowspan="3">c<sub>tot</sub> = 500C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>0.669</strong></td><td>[0.128;2.555]</td><td>0.998</td><td>2.009</td><td>[0.570;6.140]</td><td>0.998</td></tr><tr><td>LSTM</td><td>0.936</td><td>[0.664;3.677]</td><td>0.997</td><td><strong>1.639</strong></td><td>[0.816;3.147]</td><td>0.998</td></tr><tr><td>NIRB</td><td>1.318</td><td>[0.187;4.699]</td><td>0.995</td><td>1.962</td><td>[0.476;5.582]</td><td>0.998</td><td rowspan="3">750</td><td rowspan="3">2875</td><td rowspan="3">c<sub>tot</sub> = 1000C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>0.415</strong></td><td>[0.038;1.536]</td><td>0.999</td><td>1.454</td><td>[0.348;4.710]</td><td>0.998</td></tr><tr><td>LSTM</td><td>0.728</td><td>[0.465;2.908]</td><td>0.998</td><td><strong>1.358</strong></td><td>[0.614;3.104]</td><td>0.999</td></tr></tbody></table>

* The NIRB approach is particularly efficient at constructing a low-cost predictor for configurations when the physical phenomena represented with the HF model are close to the ones with the LF model (i.e the gap between the HF and LF output are not very important) or when the discrepancy between the HF and LF models is weakly parameters dependent. This seems to be particularly true when confronted to small database.

* If the difference between the HF and LF models is strongly parameters dependent, it seems preferable to turn to the GPR and LSTM approaches. In the presence of a large number of points, LSTM may be the preferred choice, as it seems more capable of managing large databases, whereas for small data, we recommend focusing on GPR.

* Then, if we need a prediction with a quantified uncertainty, the most appropriate of the three methods tested for these purposes seems to be the GPR.

* Finally, we again insist on the fact that the best combination of HF and LF evaluations to constitute the learning set (i.e. the best choice for n<sub>HF</sub> and n<sub>LF</sub>) is case-dependent, and also depends on the total admissible cost c<sub>tot</sub> and on the considered prediction method (NIRB, LSTM or GPR).

### 3.5. Interest of the multi-fidelity framework

The study proposed aligns with findings reported in recent benchmarking studies such as [20], where multi-fidelity approaches were shown to systematically outperform single-fidelity counterparts across various applications involving time-dependent functional outputs. In particular, the authors compare over a dozen surrogate modeling strategies and demonstrate, through a set of benchmark cases (e.g., viscous flow simulations, airfoil pressure profiles), that multi-fidelity models consistently provide more accurate predictions when low-fidelity data is available and well correlated with the high-fidelity outputs. The advantage is especially clear when the number of high-fidelity samples is limited, a common scenario in computational physics, and the performance gains were observed
<table><caption>Table 6: Prediction quality of $\phi_{S_i}$ and $\phi_{S_e}$ in W/m² from the three predictors using the $R_5C_4$ model (LF<sub>RC</sub>) as LF model, based on the indicators MAE*, Q<sub>95%</sub> and R<sup>2*</sup> calculated over a test set of size 300. c<sub>tot</sub> aggregates the costs associated with the HF evaluations.</caption><thead><tr><th rowspan="2">M<sup>HF</sup><sub>2b</sub></th><th colspan="3">φ<sub>S<sub>i</sub></sub></th><th colspan="3">φ<sub>S<sub>e</sub></sub></th><th rowspan="2">n<sub>HF</sub></th><th rowspan="2">c<sub>tot</sub></th></tr><tr><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th></tr></thead><tbody><tr><td>NIRB</td><td>1.474</td><td>[0.203;4.260]</td><td>0.994</td><td><b>2.225</b></td><td>[0.627; 5.248]</td><td>0.998</td><td rowspan="3">100</td><td rowspan="3">c<sub>tot</sub> = 100C<sub>HF</sub></td></tr><tr><td>GPR</td><td><b>1.164</b></td><td>[0.309;4.420]</td><td>0.995</td><td>3.388</td><td>[1.013;7.975]</td><td>0.995</td></tr><tr><td>LSTM</td><td>1.735</td><td>[1.174;6.471]</td><td>0.991</td><td>3.686</td><td>[2.066;7.128]</td><td>0.994</td></tr><tr><td>NIRB</td><td>1.337</td><td>[0.175;4.655]</td><td>0.995</td><td>1.983</td><td>[0.468;5.522]</td><td>0.998</td><td rowspan="3">500</td><td rowspan="3">c<sub>tot</sub> = 500C<sub>HF</sub></td></tr><tr><td>GPR</td><td><b>0.604</b></td><td>[0.152;2.223]</td><td>0.998</td><td><b>1.783</b></td><td>[0.439;5.660]</td><td>0.998</td></tr><tr><td>LSTM</td><td>0.684</td><td>[0.393;3.268]</td><td>0.998</td><td>1.980</td><td>[0.824;5.581]</td><td>0.998</td></tr><tr><td>NIRB</td><td>1.294</td><td>[0.196;4.425]</td><td>0.995</td><td>1.934</td><td>[0.476;5.464]</td><td>0.998</td><td rowspan="3">1000</td><td rowspan="3">c<sub>tot</sub> = 1000C<sub>HF</sub></td></tr><tr><td>GPR</td><td><b>0.601</b></td><td>[0.135;2.240]</td><td>0.998</td><td>1.651</td><td>[0.421;5.590]</td><td>0.998</td></tr><tr><td>LSTM</td><td>0.629</td><td>[0.427;2.150]</td><td>0.998</td><td><b>1.290</b></td><td>[0.671;2.542]</td><td>0.999</td></tr></tbody></table>

both in terms of precision and in computational cost-efficiency. In the application presented above, the physics-driven nature of the hierarchy between 2D, 1D and RC thermal models naturally lends itself to such multi-fidelity treatment, with a high correlation between the different levels of fidelity (see Figure 3). In addition, to further illustrate the advantages of the multi-fidelity approach in the case study presented in this article, the prediction accuracy is compared using the two-level meta-model $M_{2a}^{HF}$ (considering the one-dimensional thermal model (LF<sub>1D</sub>) as the low fidelity model) and a single-fidelity meta-model of the HF fidelity model (two-dimensional thermal model) in Table 7, both of which use the GPR and LSTM approaches. The NIRB approach is a correction of a low-fidelity calculation, so a 1-fidelity meta-model created by this method corresponds to the high-fidelity model, which is no longer a predictor.

Across all outputs, for different number of high-fidelity simulations used in the training base ($100 \le n_{HF} \le 1000$), and using the GPR and LSTM approaches, the multi-fidelity surrogate model $M_{2a}$ consistently outperforms the single-fidelity meta-model counterpart in terms of accuracy and robustness. Using the GPR method, the prediction accuracy is twice as high when applying the two-level multi-fidelity approach ($M_{2a}$) compared to the single-fidelity meta-model ($M_1$); in the case of the LSTM approach, the prediction accuracy is improved by a factor of four under the multi-fidelity setting. These results confirm the benefit of incorporating low-cost information from simpler models (e.g., RC or 1D thermal approximations) to enhance the learning of complex high-fidelity responses, particularly in the small-data regime. The use of multi-fidelity modeling not only improves predictive performance but also reduces the sensitivity to sample variability, which is crucial for practical applications where high-fidelity simulations are computationally expensive.

## 4. Conclusions

This study presents a comparison in terms of computational cost and prediction performance using different non-intrusive multifidelity models based on three approaches: the Non-Intrusive Reduced Basis method (NIRB), a Gaussian Process Regression (GPR) and a Long Short-Term Memory (LSTM). These methods are applied to the
<table><caption>Table 7: Comparison of the prediction quality of $\phi_{S_i}$ and $\phi_{S_c}$ in W/m² from the two-levels meta-model $M_{2a}^{HF}$ and the single level meta-model $M_1^{HF}$, based on the indicators MAE* and $Q_{95\%}$ calculated over a test set of size 300. $c_{tot}$ aggregates the costs associated with the HF evaluations.</caption><thead><tr><th rowspan="2">Qt</th><th rowspan="2">Model</th><th colspan="2">NIRB</th><th colspan="2">GPR</th><th colspan="2">LSTM</th><th rowspan="2">c<sub>tot</sub></th></tr><tr><th>MAE*</th><th>Q<sub>95%</sub></th><th>MAE*</th><th>Q<sub>95%</sub></th><th>MAE*</th><th>Q<sub>95%</sub></th></tr></thead><tbody><tr><td rowspan="2">φ<sub>S<sub>i</sub></sub></td><td>M<sup>HF</sup><sub>2a</sub></td><td>1.291</td><td>[0.101;4.640]</td><td>0.822</td><td>[0.140;2.918]</td><td>1.564</td><td>[1.073;5.443]</td><td rowspan="2">c<sub>tot</sub> = 100C<sub>HF</sub></td></tr><tr><td>M<sup>HF</sup><sub>1</sub></td><td>-</td><td>-</td><td>1.344</td><td>[0.303;5.784]</td><td>4.019</td><td>[2.472;23.86]</td></tr><tr><td rowspan="2">φ<sub>S<sub>c</sub></sub></td><td>M<sup>HF</sup><sub>2a</sub></td><td>1.029</td><td>[0.091;2.756]</td><td>1.261</td><td>[0.433;3.039]</td><td>3.387</td><td>[1.828;6.208]</td><td rowspan="2">c<sub>tot</sub> = 500C<sub>HF</sub></td></tr><tr><td>M<sup>HF</sup><sub>1</sub></td><td>-</td><td>-</td><td>3.720</td><td>[1.038;10.235]</td><td>16</td><td>[5.184;51.893]</td></tr><tr><td rowspan="2">φ<sub>S<sub>i</sub></sub></td><td>M<sup>HF</sup><sub>2a</sub></td><td>1.059</td><td>[0.070;4.199]</td><td>0.470</td><td>[0.063;1.664]</td><td>0.868</td><td>[0.523;3.392]</td><td rowspan="2">c<sub>tot</sub> = 500C<sub>HF</sub></td></tr><tr><td>M<sup>HF</sup><sub>1</sub></td><td>-</td><td>-</td><td>0.658</td><td>[0.149;2.818]</td><td>3.161</td><td>[1.839;22.058]</td></tr><tr><td rowspan="2">φ<sub>S<sub>c</sub></sub></td><td>M<sup>HF</sup><sub>2a</sub></td><td>0.931</td><td>[0.062;3.011]</td><td>0.579</td><td>[0.205;1.466]</td><td>1.689</td><td>[0.834;4.144]</td><td rowspan="2">c<sub>tot</sub> = 1000C<sub>HF</sub></td></tr><tr><td>M<sup>HF</sup><sub>1</sub></td><td>-</td><td>-</td><td>1.212</td><td>[0.344;3.683]</td><td>8.575</td><td>[4.031;20.994]</td></tr><tr><td rowspan="2">φ<sub>S<sub>i</sub></sub></td><td>M<sup>HF</sup><sub>2a</sub></td><td>1.036</td><td>[0.052;4.136]</td><td>0.299</td><td>[0.046;1.253]</td><td>0.556</td><td>[0.340;2.670]</td><td rowspan="2">c<sub>tot</sub> = 1000C<sub>HF</sub></td></tr><tr><td>M<sup>HF</sup><sub>1</sub></td><td>-</td><td>-</td><td>0.527</td><td>[0.114;2.903]</td><td>2.087</td><td>[1.522;8.901]</td></tr><tr><td rowspan="2">φ<sub>S<sub>c</sub></sub></td><td>M<sup>HF</sup><sub>2a</sub></td><td>0.901</td><td>[0.057;2.796]</td><td>0.571</td><td>[0.160;1.570]</td><td>1.207</td><td>[0.614;2.914]</td><td rowspan="2">c<sub>tot</sub> = 1000C<sub>HF</sub></td></tr><tr><td>M<sup>HF</sup><sub>1</sub></td><td>-</td><td>-</td><td>0.917</td><td>[0.251;2.830]</td><td>6.433</td><td>[3.348;11.354]</td></tr></tbody></table>

prediction of four time-series for a multi-layer transient thermal problem, with different specificity and difficulty. In this work we propose general models that are likely to be accurate, fast and adaptable to contexts requiring a reduced number of simulations, while being able to incorporate several levels of fidelity.

For the numerical case studied, the GPR method demonstrates its ability to quantify uncertainty, while offering the best overall accuracy, but at a non-negligible model training cost (offline phase). The NIRB method, although less accurate in cases where the gap between the HF and LF models is parameter-dependent, has the advantage of being less costly in the offline phase and may offer good performance.

In addition, we observed that the accuracy of the Recurrent Neural networks (LSTM) was highly dependent on the size of the training set. The analysis of the three proposed multifidelity models $M_{2a}^{HF}$, $M_{2b}^{HF}$ and $M_3^{HF}$ showed that the models with three levels of fidelity $M_3^{HF}$ represented the best trade-off between online cost and accuracy. This model, using the three codes, i.e. the HF, the LF<sub>1D</sub> and the LF<sub>RC</sub> models, obtained good accuracy with a reduced prediction cost (online phase).

In conclusion, the choice of the method and of the optimal configuration depends on the intended application, its numerical cost, the accuracy constraints, and on the data availability. It appears that multifidelity models may offer considerable flexibility in balancing these factors, by integrating different levels of complexity to make the best use of the simulation resources available.

To go further and improve prediction accuracy, a number of directions could be explored, such as systematic optimization of the combinations between high-fidelity (HF) and low-fidelity (LF) simulations or the iterative addition of new points in poorly-predicted areas. These lines of research open up opportunities to extend the use and performance of multifidelity models in a variety of scientific and industrial contexts.
711 Acknowledgements

This work, carried out as part of the “RESBIOBAT” project, was funded by the French National Research Agency (ANR) under the reference ANR-21-CE22-0018-01.

714 Credit author statement

715 **Hadi Nasser:** Methodology, Software, Validation, Formal analysis, Investigation, Writing - Original Draft, Visualization; **Guillaume Perrin:** Conceptualization, Methodology, Software, Validation, Formal Analysis, Writing - Original Draft, Supervision; **Rachida Chakir:** Conceptualization, Methodology, Software, Validation, Formal Analysis, Writing - Original Draft; **Séverine Demeyer:** Conceptualization, Writing - Review & Editing; **Julien Waeytens:** Conceptualization, Resources, Writing - Review & Editing, Supervision, Project administration.

720 References

[1] Leo Wai-Tsun Ng and Michael Eldred. Multifidelity uncertainty quantification using non-intrusive polynomial chaos and stochastic collocation. In *53rd aiaa/asme/asce/ahs/asc structures, structural dynamics and materials conference 20th aiaa/asme/ahs adaptive structures conference 14th aiaa*, page 1852, 2012.

[2] Max D Morris. Factorial sampling plans for preliminary computational experiments. *Technometrics*, 33(2):161–174, 1991.

[3] Ilya M Sobol. Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates. *Mathematics and computers in simulation*, 55(1-3):271–280, 2001.

[4] Alexander IJ Forrester, András Sóbester, and Andy J Keane. Multi-fidelity optimization via surrogate modelling. *Proceedings of the royal society a: mathematical, physical and engineering sciences*, 463(2088):3251–3269, 2007.

[5] Felipe AC Viana, Timothy W Simpson, Vladimir Balabanov, and Vasilli Toropov. Special section on multidisciplinary design optimization: metamodeling in multidisciplinary design optimization: how far have we really come? *AIAA journal*, 52(4):670–690, 2014.

[6] M Giselle Fernández-Godino. Review of multi-fidelity models. *arXiv preprint arXiv:1609.07196*, 2016.

[7] Theresa Robinson, Michael Eldred, Karen Willcox, and Robert Haimes. Strategies for multifidelity optimization with variable dimensional hierarchical models. In *47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference 14th AIAA/ASME/AHS Adaptive Structures Conference 7th*, page 1819, 2006.

[8] Loic Le Gratiet. *Multi-fidelity Gaussian process regression for computer experiments*. PhD thesis, Université Paris-Diderot-Paris VII, 2013.

[9] Alexander IJ Forrester, Neil W Bressloff, and Andy J Keane. Optimization using surrogate models and partially converged computational fluid dynamics simulations. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 462(2071):2177–2204, 2006.

[10] Alexander IJ Forrester. Black-box calibration for complex-system simulation. *Philosophical Transactions of The Royal Society A: Mathematical, Physical and Engineering Sciences*, 368(1924):3567–3579, 2010.

[11] Marc C Kennedy and Anthony O'Hagan. Predicting the output from a complex computer code when fast approximations are available. *Biometrika*, 87(1):1–13, 2000.

[12] Loic Le Gratiet. Bayesian analysis of hierarchical multifidelity codes. *SIAM/ASA Journal on Uncertainty Quantification*, 1(1):244–269, 2013.

[13] Loic Le Gratiet and Josselin Garnier. Recursive co-kriging model for design of computer experiments with multiple levels of fidelity. *International Journal for Uncertainty Quantification*, 4(5), 2014.

[14] Shibo Li, Wei Xing, Robert Kirby, and Shandian Zhe. Multi-fidelity bayesian optimization via deep neural networks. *Advances in Neural Information Processing Systems*, 33:8521–8531, 2020.

[15] Xuhui Meng and George Em Karniadakis. A composite neural network that learns from multi-fidelity data: Application to function approximation and inverse pde problems. *Journal of Computational Physics*, 401:109020, 2020.
752

753

[16] Xinshuai Zhang, Fangfang Xie, Tingwei Ji, Zaoxu Zhu, and Yao Zheng. Multi-fidelity deep neural network surrogate model for aerodynamic shape optimization. *Computer Methods in Applied Mechanics and Engineering*, 373:113485, 2021.

754

755

[17] Xuhui Meng, Hessam Babaee, and George Em Karniadakis. Multi-fidelity bayesian neural networks: Algorithms and applications. *Journal of Computational Physics*, 438:110361, 2021.

756

757

[18] Paolo Conti, Mengwu Guo, Andrea Manzoni, and Jan S Hesthaven. Multi-fidelity surrogate modeling using long short-term memory networks. *Computer methods in applied mechanics and engineering*, 404:115811, 2023.

758

759

[19] Baptiste Kerleguer. Multifidelity surrogate modeling for time-series outputs. *SIAM/ASA Journal on Uncertainty Quantification*, 11(2):514–539, 2023.

760

761

762

[20] Lucas Brunel, Mathieu Balesdent, Loïc Brevault, Rodolphe Le Riche, and Bruno Sudret. A survey on multi-fidelity surrogates for simulators with functional outputs: Unified framework and benchmark. *Computer Methods in Applied Mechanics and Engineering*, 435:117577, 2025.

763

764

[21] Guillaume Perrin. Adaptive calibration of a computer code with time-series output. *Reliability engineering & system safety*, 196:106728, 2020.

765

766

767

[22] Hadi Nasser, Guillaume Perrin, Rachida Chakir, Séverine Demeyer, and Julien Waeytens. Incorporating multi-source uncertainties in fast building wall thermal resistance estimation through physics-based and multi-fidelity statistical learning models. *Journal of Building Engineering*, page 111027, 2024.

768

769

[23] Rémi Stroh. *Planification d'expériences numériques en multi-fidélité: Application à un simulateur d'incendies*. PhD thesis, Université Paris Saclay (COmUE), 2018.

770

771

[24] Jan S. Hesthaven, Gianluigi Rozza, and Benjamin Stamm. *Certified Reduced Basis Methods for Parametrized Partial Differential Equations*. 2015.

772

[25] Alfio Quarteroni, Andrea Manzoni, and Federico Negri. *Reduced basis methods for partial differential equations: An introduction*. 2015.

773

774

[26] Guglielmo Padula, Michele Girfoglio, and Gianluigi Rozza. A brief review of reduced order models using intrusive and non-intrusive techniques. *PAMM*, 24(4):e202400210, 2024.

775

776

[27] C. Prud'homme, D.V. Rovas, K. Veroy, L. Machiels, Y. Maday, A.T. Patera, and G. Turinici. Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods. 124(1):70–80, 2002.

777

778

[28] Martin A. Grepl. Certified reduced basis methods for nonaffine linear time-varying and nonlinear parabolic partial differential equations. *Mathematical Models and Methods in Applied Sciences*, 22(3), 2012.

779

780

[29] Rachida Chakir and Janelle K Hammond. A non-intrusive reduced basis method for elastoplasticity problems in geotechnics. *Journal of Computational and Applied Mathematics*, 337:1–17, 2018.

781

782

[30] Elise Grosjean and Yvon Maday. Error estimate of the non-intrusive reduced basis (nirb) two-grid method with parabolic equations. *The SMAI Journal of computational mathematics*, 9:227–256, 2023.

783

784

[31] Mariella Kast, Mengwu Guo, and Jan S. Hesthaven. A non-intrusive multifidelity method for the reduced order modeling of nonlinear problems. *Computer Methods in Applied Mechanics and Engineering*, 364:112947, 2020.

785

786

[32] Rachida Chakir and Yvon Maday. A two-grid finite-element/reduced basis scheme for the approximation of the solution of parameter dependent PDE. In *9e Colloque national en calcul des structures*, Giens, France, May 2009. CSMA.

787

788

[33] Michael D McKay, Richard J Beckman, and William J Conover. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. *Technometrics*, 42(1):55–61, 2000.

789

790

[34] Rachida Chakir, Yvon Maday, and Philippe Parnaudeau. A non-intrusive reduced basis approach for parametrized heat transfer problems. *Journal of Computational Physics*, 376:617–633, 2019.

791

[35] Ian T Jolliffe. *Principal component analysis for special types of data*. Springer, 2002.

792

793

[36] Stefano Conti and Anthony O'Hagan. Bayesian emulation of complex multi-output and dynamic computer models. *Journal of statistical planning and inference*, 140(3):640–651, 2010.

794

795

[37] Svante Wold, Kim Esbensen, and Paul Geladi. Principal component analysis. *Chemometrics and intelligent laboratory systems*, 2(1-3):37–52, 1987.

796

797

[38] Guillaume Perrin, Christian Soize, Denis Duhamel, and C Funfschilling. Karhunen–loève expansion revisited for vector-valued random fields: Scaling, errors and optimal basis. *Journal of Computational Physics*, 242:607–622, 2013.

798

[39] TJ Santner. The design and analysis of computer experiments, 2003.

799

[40] Robert B Gramacy. *Surrogates: Gaussian process modeling, design, and optimization for the applied sciences*. Chapman and Hall/CRC,
[41] Michael A Nielsen. Neural networks and deep learning, volume 25. Determination press San Francisco, CA, USA, 2015.

[42] S Hochreiter. Long short-term memory. Neural Computation MIT-Press, 1997.

[43] Elad Hoffer, Itay Hubara, and Daniel Soudry. Train longer, generalize better: closing the generalization gap in large batch training of neural networks. Advances in neural information processing systems, 30, 2017.

[44] Jason Brownlee. Long short-term memory networks with python: develop sequence prediction models with deep learning. Machine Learning Mastery, 2017.

[45] Hikmet Kerem Cigizoglu and Murat Alp. Generalized regression neural network in modelling river sediment yield. Advances in Engineering Software, 37(2):63–68, 2006.

[46] Rohit K Tripathy and Ilias Bilionis. Deep uq: Learning deep neural network surrogate models for high dimensional uncertainty quantification. Journal of computational physics, 375:565–588, 2018.

[47] Prajit Ramachandran, Barret Zoph, and Quoc V Le. Searching for activation functions. arXiv preprint arXiv:1710.05941, 2017.

[48] Benjamin Lindemann, Timo Müller, Hannes Vietz, Nasser Jazdi, and Michael Weyrich. A survey on long short-term memory networks for time series prediction. Procedia CIRP, 99:650–655, 2021.

[49] Thanh-Tung Ha, Vincent Feuillet, Julien Waeytens, Kamel Zibouche, Laurent Peiffer, Yann Garcia, Véronique Le Sant, Rémi Bouchie, Alain Koenen, Jean-Pierre Monchau, and Laurent Ibos. Measurement prototype for fast estimation of building wall thermal resistance under controlled and natural environmental conditions. Energy and Buildings, 268:112166, 2022.

[50] Thanh-Tung Ha, Vincent Feuillet, Julien Waeytens, Kamel Zibouche, Simon Thebault, Rémi Bouchie, Véronique Le Sant, and Laurent Ibos. Benchmark of identification methods for the estimation of building wall thermal resistance using active method: Numerical study for iwi and single-wall structures. Energy and Buildings, 224:110130, 2020.

## Appendix A. Derivation of the joint distribution of $(Y_s(x), \bar{Y})$

To begin, let us notice that Eq. (12) implies that, by noted $n_d$ the number of sampling used, for each $1 \le l \le s$ and all $x \in \mathcal{X}$,

$$Y_l(x) = A_l(x)\bar{W}(x), \quad \bar{W}(x) := (W_1(x), \dots, W_s(x)), \qquad (A.1)$$

$$\begin{bmatrix} A_1(x) & A_2(x) & \dots & A_s(x) \end{bmatrix}^T = \begin{bmatrix} I_{N_t} & 0 & \dots & 0 \\ B_2(x) & I_{N_t} & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ \prod_{l=2}^s B_l(x) & \prod_{l=3}^s B_l(x) & \dots & I_{N_t} \end{bmatrix}$$

As a consequence,

$$\begin{bmatrix} Y_s(x) \\ \bar{Y} \end{bmatrix} = \bar{A}(x) \begin{bmatrix} \bar{W}(x) \\ \bar{W}(x_1^{(1)}) \\ \vdots \\ \bar{W}(x_d^{(n_d)}) \end{bmatrix}, \quad \bar{A}(x) := \begin{bmatrix} A_s(x) & 0 & \dots & 0 \\ 0 & A_1(x_1^{(1)}) & \dots & \vdots \\ \vdots & \vdots & \ddots & 0 \\ 0 & \dots & 0 & A_s(x_s^{(n_d)}) \end{bmatrix}$$

from which we deduce:

$$[m(x) \quad \bar{m}]^T = \bar{A}(x) [E[\bar{W}(x)] \quad E[\bar{W}(x_1^{(1)})] \quad \dots \quad E[\bar{W}(x_d^{(n_d)})]]^T,$$
$$
\begin{bmatrix} V(\mathbf{x}) & U(\mathbf{x}) \\ U^T(\mathbf{x}) & \bar{V} \end{bmatrix} = \bar{A}(\mathbf{x}) \begin{bmatrix} \operatorname{Cov}(\bar{W}(\mathbf{x}), \bar{W}(\mathbf{x})) & \operatorname{Cov}(\bar{W}(\mathbf{x}), \bar{W}(\mathbf{x}_1^{(1)})) & \cdots & \operatorname{Cov}(\bar{W}(\mathbf{x}), \bar{W}(\mathbf{x}_d^{(n_d)})) \\ \operatorname{Cov}(\bar{W}(\mathbf{x}_1^{(1)}), \bar{W}(\mathbf{x})) & \operatorname{Cov}(\bar{W}(\mathbf{x}_1^{(1)}), \bar{W}(\mathbf{x}_1^{(1)})) & \cdots & \vdots \\ \vdots & \vdots & \ddots & \operatorname{Cov}(\bar{W}(\mathbf{x}_d^{(n_d-1)}), \bar{W}(\mathbf{x}_d^{(n_d)})) \\ \operatorname{Cov}(\bar{W}(\mathbf{x}_d^{(n_d)}), \bar{W}(\mathbf{x})) & \cdots & \cdots & \operatorname{Cov}(\bar{W}(\mathbf{x}_d^{(n_d)}), \bar{W}(\mathbf{x}_d^{(n_d)})) \end{bmatrix} \bar{A}(\mathbf{x})^T.
$$

825 where for each $\mathbf{x}, \mathbf{x}' \in \mathcal{X}$,

$$
\mathbb{E}[\bar{W}(\mathbf{x})] = \begin{bmatrix} \mu_1(\mathbf{x}) \\ \vdots \\ \mu_d(\mathbf{x}) \end{bmatrix}, \quad \operatorname{Cov}(\bar{W}(\mathbf{x}), \bar{W}(\mathbf{x}')) = \operatorname{diag} \begin{bmatrix} C_1(\mathbf{x}, \mathbf{x}') & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & C_2(\mathbf{x}, \mathbf{x}') & \cdots & \vdots \\ \vdots & \vdots & \ddots & \mathbf{0} \\ \mathbf{0} & \cdots & \mathbf{0} & C_d(\mathbf{x}, \mathbf{x}') \end{bmatrix}.
$$

## Appendix B. Derivation of the statistical moments of the s conditioned Gaussian processes

In this section, to simplify the reading of the expressions, we denote by $\mathcal{X}_1, \dots, \mathcal{X}_s$ the numerical designs of experiments gathering respectively all the evaluation points of the codes $y_1, \dots, y_s$. Note that in the case where $S(n_1, \dots, n_s)$ has a nested structure, $\mathcal{X}_s \subset \mathcal{X}_{s-1} \subset \dots \subset \mathcal{X}_1$.

In addition, for all function $f$ defined on $\mathcal{X}$ with values in $\mathbb{R}^p$, all function $h$ defined on $\mathcal{X} \times \mathcal{X}$ with values in $N_t \times N_t$, and all $1 \le l \le s$, we note:

$$
f(\mathcal{X}_l) = (f(\mathbf{x}_l^{(1)}), \dots, f(\mathbf{x}_l^{(n_l)})), \quad h(\mathcal{X}_l, \mathcal{X}_l) = \begin{bmatrix} \mathbf{h}(\mathbf{x}_l^{(1)}, \mathbf{x}_l^{(1)}) & \cdots & \mathbf{h}(\mathbf{x}_l^{(1)}, \mathbf{x}_l^{(n_l)}) \\ \vdots & \ddots & \vdots \\ \mathbf{h}(\mathbf{x}_l^{(n_l)}, \mathbf{x}_l^{(1)}) & \cdots & \mathbf{h}(\mathbf{x}_l^{(n_l)}, \mathbf{x}_l^{(n_l)}) \end{bmatrix}.
$$

832 Given these notations, we can proceed sequentially.

• **Initialisation:** for $l = 1$, we denote by

$$
\mathbf{Y}_1^{(n_1)} \sim \mathcal{GP}(\mu_1^{(n_1)}, C_1^{(n_1)}), \qquad (\text{B.1})
$$

834 the Gaussian process $\mathbf{Y}_1$ once conditioned by the $n_1$ evaluations of $\mathbf{y}_1$. By Gaussian conditioning, for each $\mathbf{x}, \mathbf{x}' \in \mathcal{X}$, we have:

$$
\mu_1^{(n_1)}(\mathbf{x}) := \mu_1(\mathbf{x}) + C_1^{(n_1)}(\mathbf{x}, \mathcal{X}_1) (C_1^{(n_1)}(\mathcal{X}_1, \mathcal{X}_1))^{-1} (\mathbf{y}_1(\mathcal{X}_1) - \mu_1(\mathcal{X}_1)), \quad (\text{B.2})
$$

$$
836 \qquad C_1^{(n_1)}(\mathbf{x}, \mathbf{x}') := C_1(\mathbf{x}, \mathbf{x}') - C_1^{(n_1)}(\mathbf{x}, \mathcal{X}_1) (C_1^{(n_1)}(\mathcal{X}_1, \mathcal{X}_1))^{-1} C_1^{(n_1)}(\mathcal{X}_1, \mathbf{x}'). \qquad (\text{B.3})
$$

837 • **Iteration:** for each $2 \le m \le l \le s$, we then denote by

$$
\mathbf{Y}_{\ell}^{(n_1, \dots, n_m)} \sim \mathcal{GP}(\mu_{\ell}^{(n_1, \dots, n_m)}, C_{\ell}^{(n_1, \dots, n_m)}), \qquad (\text{B.4})
$$

838 the Gaussian processes $\mathbf{Y}_\ell$ once conditioned by all the evaluations of $\mathbf{y}_1, \dots, \mathbf{y}_m$ (but not the ones of $\mathbf{y}_{m+1}, \dots, \mathbf{y}_\ell$).
839 Assuming that the different mean and covariance functions have already been calculated up to iteration $\ell - 1$ from the
840 evaluations of $\mathbf{y}_1, \dots, \mathbf{y}_{\ell-1}$, we then have:

$$
\mu_{\ell}^{(n_1, \dots, n_{\ell-1})}(\mathbf{x}) := \mathbf{B}_{\ell}(\mathbf{x}) \mu_{\ell-1}^{(n_1, \dots, n_{\ell-1})}(\mathbf{x}) + \mu_{\ell}(\mathbf{x}), \qquad (\text{B.5})
$$

$$
841 \qquad C_{\ell}^{(n_1, \dots, n_{\ell-1})}(\mathbf{x}, \mathbf{x}') := \mathbf{B}_{\ell}(\mathbf{x}) C_{\ell-1}^{(n_1, \dots, n_{\ell-1})} \mathbf{B}_{\ell}(\mathbf{x}')^T + C_{\ell}(\mathbf{x}, \mathbf{x}'). \qquad (\text{B.6})
$$

842 We deduce:

$$
\mu_{\ell}^{(n_1, \dots, n_{\ell-1})}(\mathbf{x}) := \mu_{\ell}^{(n_1, \dots, n_{\ell-1})}(\mathbf{x}) + C_{\ell}^{(n_1, \dots, n_{\ell-1})}(\mathbf{x}, \mathcal{X}_{\ell}) (C_{\ell}^{(n_1, \dots, n_{\ell-1})}(\mathcal{X}_{\ell}, \mathcal{X}_{\ell}))^{-1} (\mathbf{y}_{\ell}(\mathcal{X}_{\ell}) - \mu_{\ell}^{(n_1, \dots, n_{\ell-1})}(\mathcal{X}_{\ell})), \quad (\text{B.7})
$$

$$
843 \qquad C_{\ell}^{(n_1, \dots, n_{\ell})}(\mathbf{x}, \mathbf{x}') := C_{\ell}^{(n_1, \dots, n_{\ell-1})}(\mathbf{x}, \mathcal{X}_{\ell}) - C_{\ell}^{(n_1, \dots, n_{\ell-1})}(\mathbf{x}, \mathcal{X}_{\ell}) (C_{\ell}^{(n_1, \dots, n_{\ell-1})}(\mathcal{X}_{\ell}, \mathcal{X}_{\ell}))^{-1} C_{\ell}^{(n_1, \dots, n_{\ell-1})}(\mathcal{X}_{\ell}, \mathbf{x}'). \quad (\text{B.8})
$$
## Appendix C. Adaptation of the PCA algorithm for multifidelity data

A reduction in dimension using a weighted principal component analysis (PCA) is realized as proposed in [38]. The objective is to apply a weighted PCA to a multifidelity data by optimizing the weights assigned to each fidelity level and determining the optimal number of principal components. This process is presented in Algorithm 1. It guarantees minimum reconstruction error within a defined tolerance.

<caption>**Algorithm 1:** Weighted optimization of fidelity level weights and principal components</caption>

**Input:**
- $\mathbf{Y}^{(\ell)} \in \mathbb{R}^{N_i \times n_\ell}$, $\ell \in \{1, \dots, s\}$: outputs for the $s$ models;
- $\tau$ and $\epsilon$: tolerance thresholds for the projection error and the convergence of the weights.

**Initialisation:**
- Compute the initial weights: $w_\ell \leftarrow \frac{1}{\sum_{j=1}^s 1/n_j}$, $\ell \in \{1, \dots, s\}$.
- **Combine the weighted data:** $\mathbf{Y}_{c, w} \leftarrow [\mathbf{Y}^{(1)} \cdot w_1 \quad \mathbf{Y}^{(2)} \cdot w_2 \quad \dots \quad \mathbf{Y}^{(s)} \cdot w_s]$.
- $k \leftarrow 1$, $\text{Error}_{(s)} \leftarrow \infty$.

**while** $\text{Error}_{(s)} > \tau$ **do**
    **repeat**
        Perform PCA with SVD:
        - Compute the SVD of $\mathbf{Y}_{c, w} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^\top$, where $(\mathbf{\Sigma})^\top$ stands for the transpose operator.
        - Extract the projection matrix $\mathbf{\Pi} = \mathbf{U}[:, 1: k]$.
        Compute the projection errors for each fidelity level:
        $$ \text{Error}_{(\ell)} \leftarrow \frac{1}{n_\ell N_\ell} \sum_{i=1}^{n_\ell} \sum_{j=1}^{N_\ell} \left| \mathbf{Y}_{k,i,j}^{(\ell)} - \mathbf{Y}_{i,j}^{(\ell)} \right|, \text{ with } \mathbf{Y}_k^{(\ell)} = \mathbf{\Pi}\mathbf{\Pi}^\top \mathbf{Y}^{(\ell)}, \quad 1 \le \ell \le s. $$
        optimize the weights:
        - Update the weights: $w_\ell^{\text{old}} \leftarrow w_\ell$, $w_\ell^{\text{new}} \leftarrow w_\ell^{\text{old}} \cdot \text{Error}_{(\ell)}$, $1 \le \ell \le s$.
        - Normalize the weights: $w_\ell \leftarrow \frac{w_\ell^{\text{new}}}{\sum_{j=1}^s w_\ell^{\text{new}}}$, $1 \le \ell \le s$.
        - **Update the weighted matrix data:** $\mathbf{Y}_{c, w} \leftarrow [\mathbf{Y}^{(1)} \cdot w_1 \quad \mathbf{Y}^{(2)} \cdot w_2 \quad \dots \quad \mathbf{Y}^{(s)} \cdot w_s]$.
        until convergence ($\text{Error}_{(s)}^{\text{new}} - \text{Error}_{(s)}^{\text{old}} < \epsilon$);
        Increment the number of principal components:
        $k \leftarrow k + 1$; $\text{Error}_{(s)} \leftarrow \text{Error}_{(s)}^{\text{new}}$.
    end
    **return** $w_1, w_2, \dots, w_s$, $\mathbf{\Pi}$.

## Appendix D. Example of LSTM code

In this section, we present one of the LSTM structures that has been used for the numerical applications presented in this work. The idea is to illustrate the simplified nature of such an implementation, even for non-specialists. We are limiting ourselves here to an implementation with the Keras library, but other libraries can work just as well.

853 import pandas as pd

854 from tensorflowkeras.models import Sequential

855 from tensorflowkeras layers import LSTM, Dense, Dropout, Input

856 from tensorflowkeras.optimizers import Adam

857 5
```python
858 6 # Function to create an LSTM model for time series prediction
859 7 def lstm_model(input_size, timesteps, learning_rate, lstm_units, dropout_rates):
860 8 model = Sequential() # Sequential model allowing layers to be stacked linearly
861 9 model.add(Input shape=(timesteps, input_size))) # Input layer defining time steps and features
862 10 for units, dr in zip(lstm_units, dropout_rates):
863 11 model.add(LSTM(units, return_sequences=True)) # Adding LSTM layers with units neurons
864 12 model.add(Dropout(dr)) # Dropout for regularization (reducing overfitting)
865 13 model.add(Dense(timesteps, activation='linear')) # Output layer with linear activation
866 14
867 15 # Compiling the model using Mean Squared Error (MSE) loss
868 16 model.compile(optimizer=Adam(learning_rate), loss='mse', metrics=['mae'])
869 17 return model
870 18
871 19 # Load Low-Fidelity (LF) simulation data and split into training and test sets
872 20 LF_data = pd.read_csv('path_LF')
873 21 LF_train, LF_test = LF_data[:, :Ntrain], LF_data[:Ntrain:Ntrain+Ntest] # Train and test split
874 22
875 23 # Load High-Fidelity (HF) simulation data and split into training and test sets
876 24 HF_data = pd.read_csv('path_HF')
877 25 HF_train, HF_test = HF_data[:, :Ntrain], HF_data[:Ntrain:Ntrain+Ntest]
878 26
879 27 # Load input parameters (features) for the model
880 28 Input = pd.read_csv('pathInputs')
881 29
882 30 # Concatenating input parameters with LF simulation data for training and testing
883 31 Input_train, Input_test = (Input[:Ntrain,:], LF_train), (Input[Ntrain:Ntrain+Ntest,:], LF_test)
884 32
885 33 # Train the LSTM model using the training dataset
886 34 model = train_model(lstm_model(), Input_train, HF_train, Input_test, HF_test, batchSize)
887 35
888 36 # Perform predictions on the test set
889 37 y_pred = model Predict(Input_test)
```

## Appendix E. Thermal models associated with the three fidelity levels

The proposed approaches presented in Section 2 are applied to a transient thermal problem with three levels of fidelity presented in the following. The high fidelity (HF) model corresponds to a two-dimensional axisymmetric thermal model, equivalent to a three dimensional thermal model. By simplifying the physics of the HF model, two low fidelity (LF) thermal models are derived: a one-dimensional thermal model (LF<sub>1D</sub>) and a Resistance-Capacitance thermal model (LF<sub>RC</sub>).

### Appendix E.1. High fidelity model: 2D axisymmetric thermal model

The HF model is a 2D axisymmetric thermal partial differential equation (PDE) in cylindrical coordinate for a wall with four layers of thickness $e_{i,1 \le i \le 4}$ and a thermal excitation zone on the internal wall surface (noted $\Gamma_c = \{r \in [0, 1], z = 0\}$) of size $h_c$. We define by $\Gamma_{\text{ext}} = \{r \in [0, 4], z = 4\}$ the external surface of the wall, $\Gamma_{\text{int}} = \{r \in [0, 4], z = 0\}$ the internal surface of the wall and $\Gamma_s = \{r = 0, z \in [0, 4]\} \cup \{r = 4, z \in [0, 4]\}$ the bottom and top surface of the wall. We thus set $\Omega = \cup \Omega_i$, $\Omega_i = \{r \in [0, 4], z \in [i-1, i]\}$ with $1 \le i \le 4$, $\partial\Omega = \Gamma_{\text{ext}} \cup \Gamma_{\text{int}} \cup \Gamma_s$ and $\Gamma_{i,j} = \Omega_i \cap \Omega_j$ the interface boundary between the layer $i$ and $j$. For all $1 \le i \le 4$, let $T_i = T(r, z, t)|_{\Omega_i}$ be the temperature in layer $i$ at radius $r \in [0, 4]$, at depth $z \in [i-1, i]$ and time $t \in [0, t_f]$, where $t_f$ is the final time of the simulated period. We note $\rho_i$, $c_i$ and $\lambda_i$ are respectively the density, the heat capacity per unit mass, and thermal conductivity of layer $i$. The governing equations for the 2D axisymmetric problem in cylindrical coordinate are
$$
\left\{
\begin{array}{l}
e_i \rho_i c_i \frac{\partial T_i}{\partial t} - \left(\frac{\lambda_i e_i}{h_c^2}\right) \frac{1}{r} \frac{\partial}{\partial r} \left(r \frac{\partial T_i}{\partial r}\right) - \frac{\lambda_i}{e_i} \frac{\partial^2 T_i}{\partial z^2} = 0 \quad \text{in } \Omega_i \times ]0, t_f]
\\[1.5em]
T_i(t=0) = T_i^0 \quad \text{in } \Omega_i,
\\[1.5em]
\frac{\lambda_i}{e_i} \frac{\partial T_i}{\partial z} = \frac{\lambda_j}{e_j} \frac{\partial T_j}{\partial z}, \quad \text{and } T_i = T_j \text{ on } \Gamma_{ij} \times ]0, t_f],
\\[1.5em]
-\frac{\lambda_i}{e_i} \frac{\partial T_i}{\partial r} = 0 \quad \text{on } \Gamma_s \cap \Omega_i \times ]0, t_f]
\\[1.5em]
-\frac{\lambda_4}{e_4} \frac{\partial T_4}{\partial z} = h_{\text{ext}}(T_4 - T_{\text{ext}}(t)) - \alpha \phi_{\text{sun}}(t) \text{ on } \Gamma_{\text{ext}} \times ]0, t_f]
\\[1.5em]
-\frac{\lambda_1}{e_1} \frac{\partial T_1}{\partial z} = h_{\text{int}}(T_{\text{int}}(t) - T_1) \text{ on } \Gamma_{\text{int}} \setminus \Gamma_c \times ]0, t_f]
\\[1.5em]
-\frac{\lambda_1}{e_1} \frac{\partial T_1}{\partial r} = h_{\text{int}}(T_{\text{int}}(t) - T_1) + \phi_{\text{exc}}(t) \text{ on } \Gamma_c \times ]0, t_f].
\end{array}
\right.
\qquad \text{(E.1)}
$$

where $\phi_{\text{sun}}$ is the global incident solar radiation, $\alpha$ is the external solar absorption coefficient, $\phi_{\text{exc}}$ is the heat flux excitation on the heated internal wall surface $\Gamma_c$, $h_{\text{int}}$ and $h_{\text{ext}}$ are respectively the indoor and outdoor global heat exchange coefficients. Let us notice that the quantities of interest ($T_{S_i}$, $T_{S_e}$, $\varphi_{S_i}$, $\varphi_{S_e}$) derived with this model, can be written:

$$ T_{S_i}(t) = T_1(0, 0, t), \quad T_{S_e}(t) = T_4(0, 4, t), \quad \varphi_{S_i}(t) = -\frac{\lambda_1}{e_1} \frac{\partial T_1(0, 0, t)}{\partial z}, \quad \varphi_{S_e}(t) = -\frac{\lambda_4}{e_4} \frac{\partial T_4}{\partial z}(0, 4, t) \quad \text{for } t \in [0, t_f] \quad \text{(E.2)} $$

## Appendix E.2. One-dimensional thermal model

For all $1 \le i \le 4$, let $T_i(z_i, t)$ be the temperature in layer $i$ at depth $z_i \in [0, 1]$ and time $t \in [0, t_f]$ and $T_{0i}$ the initial temperature (at $t=0$) in layer $i$ at depth $z_i \in [0, 1]$. The scaled 1D partial differential equation is therefore given by

$$
\left\{
\begin{array}{l}
e_i \rho_i c_i \frac{\partial T_i}{\partial t} - \frac{\lambda_i}{e_i} \frac{\partial^2 T_i}{\partial z_i^2} = 0 \quad \text{in } [0, 1] \times ]0, t_f], \quad \text{for } 1 \le i \le 4 \\
\\
T_i(t=0) = T_{0i} \quad \text{in } [0, 1], \quad \text{for } 1 \le i \le 4 \\
\\
-\frac{\lambda_4}{e_4} \frac{\partial T_4}{\partial z_4} = h_{\text{ext}}(T_4 - T_{\text{ext}}(t)) - \alpha \phi_{\text{sun}}(t) \quad \text{on } \{1\} \times ]0, t_f] \\
\\
-\frac{\lambda_1}{e_1} \frac{\partial T_1}{\partial z_1} = h_{\text{int}}(T_{\text{int}}(t) - T_1) + \phi_{\text{exc}}(t) \quad \text{on } \{0\} \times ]0, t_f]
\end{array}
\right.
\qquad \text{(E.3)}
$$

In addition the continuity conditions must be ensured:

$$ T_i(z_i = 1, t) = T_{i+1}(z_{i+1} = 0, t), \quad \text{&} \quad \frac{\lambda_i}{e_i} \frac{\partial T_i}{\partial z}(z_i = 1, t) = \frac{\lambda_j}{e_{i+1}} \frac{\partial T_{i+1}}{\partial z}(z_{i+1} = 0, t) \quad 1 \le i \le 3 \quad \text{and } t \in [0, t_f] \quad \text{(E.4)} $$

Let us notice that the quantities of interest ($T_{S_i}$, $T_{S_e}$, $\varphi_{S_i}$, $\varphi_{S_e}$) derived with the 1D model, can be written:

$$ T_{S_i}(t) = T_1(0, t), \quad T_{S_e}(t) = T_4(1, t), \quad \varphi_{S_i}(t) = -\frac{\lambda_1}{e_1} \frac{\partial T_1(0, t)}{\partial z}, \quad \varphi_{S_e}(t) = -\frac{\lambda_4}{e_4} \frac{\partial T_4}{\partial z}(1, t) \quad \text{for } t \in [0, t_f] \quad \text{(E.5)} $$

## Appendix E.3. RC Thermal model

The considered RC model (denoted as $R_5C_4$) is defined by 5 resistances and 4 capacitances

$$ R_A = \frac{e_1}{2\lambda_1}, \quad R_B = \frac{e_1}{2\lambda_1} + \frac{e_2}{2\lambda_2}, \quad R_C = \frac{e_2}{2\lambda_2} + \frac{e_3}{2\lambda_3}, \quad R_D = \frac{e_3}{2\lambda_3} + \frac{e_4}{2\lambda_4}, \quad R_E = \frac{e_4}{2\lambda_4}, \quad \text{and } C_i = e_i \rho_i C_i \text{ for } i \in [1, 4]. $$

To determine the 4 wall nodal temperatures ($T_{m_1}$, $T_{m_2}$, $T_{m_3}$, $T_{m_4}$), we need to solve the system defined below for $t \in [t_0, t_f]$:

$$ C \frac{\partial T_{RC}}{\partial t} + K T_{RC} = F, \qquad \text{(E.6)} $$

$$ \text{where } T_{RC} = (T_{m_1}, T_{m_2}, T_{m_3}, T_{m_4})^T, \quad F = \left( \frac{1}{R_A + R_{S_i}} \text{T}_{\text{int}} + \frac{R_{S_i}}{R_A + R_{S_i}} \phi_{\text{exc}}, 0, 0, \frac{1}{R_E + R_{S_e}} \text{T}_{\text{ext}} + \frac{R_{S_e}}{R_E + R_{S_e}} \phi_{\text{sun}} \right)^T, $$
$$ K = \begin{pmatrix} \frac{1}{R_A + R_{S_i}} + \frac{1}{R_B} & -\frac{1}{R_B} & 0 & 0 \\ -\frac{1}{R_B} & \frac{1}{R_B} + \frac{1}{R_C} & -\frac{1}{R_C} & 0 \\ 0 & -\frac{1}{R_C} & \frac{1}{R_C} + \frac{1}{R_D} & -\frac{1}{R_D} \\ 0 & 0 & -\frac{1}{R_D} & \frac{1}{R_D} + \frac{1}{R_E + R_{S_e}} \end{pmatrix}, \quad C = \text{diag}(C_1, C_2, C_3, C_4). \quad \text{(E.7)} $$
916 where $R_{S_i} = \frac{1}{h_{int}}$, $R_{S_e} = \frac{1}{h_{ext}}$.

917 Based on determined ($T_{m1}$, $T_{m2}$, $T_{m3}$, $T_{m4}$), the quantities of interest ($T_{S_i}$, $T_{S_e}$, $\varphi_{S_i}$, $\varphi_{S_e}$) are written:

$$ egin{cases} T_{S_i}(t) = \left(\frac{R_A}{R_{S_i}+R_A}\right) T_{\text{int}}(t) + \left(\frac{R_{S_i}}{R_{S_i}+R_A}\right) T_{m1}(t) + \left(\frac{R_A R_{S_i}}{R_{S_i}+R_A}\right) \phi_{\text{exc}}(t) & \text{for } t \in [0, t_f] \\ T_{S_e}(t) = \left(\frac{R_E}{R_{S_e}+R_E}\right) T_{\text{ext}}(t) + \left(\frac{R_{S_e}}{R_{S_e}+R_E}\right) T_{m4}(t) + \left(\frac{R_E R_{S_e}}{R_{S_e}+R_E}\right) \phi_{\text{sun}}(t) & \text{for } t \in [0, t_f] \\ \varphi_{S_i}(t) = \frac{T_{S_i}(t)-T_{m2}}{R_A}, \quad \varphi_{S_e}(t) = \frac{T_{m4}-T_{S_e}(t)}{R_E} & \text{for } t \in [0, t_f] \end{cases} \quad (\text{E.8}) $$

## Appendix F. Range of variation of the input parameters

919 The range of variation of the input parameters are presented in the following:

920 1. Transfer coefficients (W/(m².K)) : $5 \le h_{int} \le 9.3$; $4.3 \le h_{ext} \le 37$,

921 2. Thickness of layer i $e_i(m)$ : $0.01 \le e_1 \le 0.03$; $0.04 \le e_2 \le 0.2$; $0.15 \le e_3 \le 0.3$; $0.01 \le e_4 \le 0.03$

922 3. Conductivity of layer i $\lambda_i$ (JK⁻¹m⁻¹) : $0.2 \le \lambda_1 \le 0.8$; $0.02 \le \lambda_2 \le 0.06$; $0.1 \le \lambda_3 \le 2.3$; $0.3 \le \lambda_4 \le 1.8$,

923 4. Volumetric heat capacity of layer i $\rho C_i$ (kJ/Km³) : $600 \le \rho C_1 \le 1500$; $27 \le \rho C_2 \le 210$; $560 \le \rho C_3 \le 2500$; $500 \le \rho C_4 \le 2000$.

## Appendix G. Prediction results for the surfaces temperatures

<table><caption>Table G.8: Prediction quality of $T_{S_i}$ and $T_{S_e}$ in °C from the three predictors using the one-dimensional thermal model (LF<sub>1D</sub>) as low fidelity model, based on the indicators MAE*, Q<sub>95%</sub> and R<sup>2*</sup> calculated over a test set of size 300. c<sub>tot</sub> aggregates the costs associated with the HF and LF<sub>1D</sub> evaluations.</caption><thead><tr><th rowspan="2">M<sub>2a</sub><sup>HF</sup></th><th colspan="3">T<sub>S<sub>i</sub></sub></th><th colspan="3">T<sub>S<sub>e</sub></sub></th><th rowspan="2">n<sub>HF</sub> = n<sub>LF</sub></th><th rowspan="2">c<sub>tot</sub></th></tr><tr><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th></tr></thead><tbody><tr><td>NIRB</td><td>0.174</td><td>[0.014;0.586]</td><td>0.999</td><td><strong>0.061</strong></td><td>[0.004;0.219]</td><td>0.999</td><td rowspan="3">92</td><td rowspan="3">100C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>0.147</strong></td><td>[0.016;0.488]</td><td>0.999</td><td>0.112</td><td>[0.025;0.360]</td><td>0.999</td></tr><tr><td>LSTM</td><td>0.715</td><td>[0.542;1.039]</td><td>0.994</td><td>0.448</td><td>[0.227;0.971]</td><td>0.987</td></tr><tr><td>NIRB</td><td>0.140</td><td>[0.007;0.525]</td><td>0.999</td><td><strong>0.059</strong></td><td>[0.004;0.213]</td><td>0.999</td><td rowspan="3">460</td><td rowspan="3">500C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>0.106</strong></td><td>[0.009;0.471]</td><td>0.999</td><td>0.074</td><td>[0.015;0.228]</td><td>0.999</td></tr><tr><td>LSTM</td><td>0.294</td><td>[0.201;0.803]</td><td>0.998</td><td>0.192</td><td>[0.101;0.354]</td><td>0.997</td></tr><tr><td>NIRB</td><td>0.136</td><td>[0.009;0.524]</td><td>0.999</td><td>0.054</td><td>[0.004;0.197]</td><td>0.999</td><td rowspan="3">920</td><td rowspan="3">1000C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>0.093</strong></td><td>[0.007;0.434]</td><td>0.999</td><td><strong>0.050</strong></td><td>[0.009;0.187]</td><td>0.999</td></tr><tr><td>LSTM</td><td>0.410</td><td>[0.326;0.534]</td><td>0.997</td><td>0.181</td><td>[0.089;0.371]</td><td>0.998</td></tr></tbody></table>
<table><caption>Table G.9: Prediction quality of $T_{S_i}$ and $T_{S_e}$ in °C from the three predictors using the $R_5C_4$ model (LF<sub>RC</sub>) and the one-dimensional thermal model (LF<sub>ID</sub>) as LF models, based on the indicators MAE*, Q<sub>95%</sub> and R<sup>2*</sup> calculated over a test set of size 300. c<sub>tot</sub> aggregates the costs associated with the HF and LF<sub>ID</sub> evaluations.</caption><thead><tr><th rowspan="2">M<sup>HF</sup><sub>(3,1)</sub></th><th colspan="3">T<sub>S<sub>i</sub></sub></th><th colspan="3">T<sub>S<sub>e</sub></sub></th><th rowspan="2">n<sub>HF</sub> = n<sub>LF</sub></th><th rowspan="2">c<sub>tot</sub></th></tr><tr><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th></tr></thead><tbody><tr><td>NIRB</td><td>0.207</td><td>[0.025;0.731]</td><td>0.998</td><td><strong>0.186</strong></td><td>[0.039;<strong>0.691</strong>]</td><td>0.998</td><td rowspan="3">92</td><td rowspan="3">100C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>0.177</strong></td><td>[0.040;<strong>0.531</strong>]</td><td>0.999</td><td>0.328</td><td>[0.078;1.420]</td><td>0.994</td></tr><tr><td>LSTM</td><td>0.621</td><td>[0.522;1.085]</td><td>0.995</td><td>0.448</td><td>[0.212;0.915]</td><td>0.987</td></tr><tr><td>NIRB</td><td>0.175</td><td>[0.025;0.623]</td><td>0.999</td><td><strong>0.139</strong></td><td>[0.026;0.557]</td><td>0.999</td><td rowspan="3">460</td><td rowspan="3">500C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>0.130</strong></td><td>[0.027;<strong>0.530</strong>]</td><td>0.999</td><td>0.257</td><td>[0.051;1.107]</td><td>0.996</td></tr><tr><td>LSTM</td><td>0.386</td><td>[0.288;0.830]</td><td>0.998</td><td>0.262</td><td>[0.146;<strong>0.438</strong>]</td><td>0.996</td></tr><tr><td>NIRB</td><td>0.168</td><td>[0.022;0.623]</td><td>0.999</td><td><strong>0.123</strong></td><td>[0.023;0.459]</td><td>0.999</td><td rowspan="3">920</td><td rowspan="3">1000C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>0.106</strong></td><td>[0.020;<strong>0.447</strong>]</td><td>0.999</td><td>0.137</td><td>[0.032;0.493]</td><td>0.999</td></tr><tr><td>LSTM</td><td>0.309</td><td>[0.246;0.472]</td><td>0.998</td><td>0.128</td><td>[0.068;<strong>0.221</strong>]</td><td>0.999</td></tr></tbody></table>

<table><caption>Table G.10: Prediction quality of $T_{S_i}$ and $T_{S_e}$ in °C from the three predictors using the $R_5C_4$ model (LF<sub>RC</sub>) and the one-dimensional thermal model (LF<sub>ID</sub>) as LF models, based on the indicators MAE*, Q<sub>95%</sub> and R<sup>2*</sup> calculated over a test set of size 300. c<sub>tot</sub> aggregates the costs associated with the HF and LF<sub>ID</sub> evaluations. In this case we choose n<sub>LF</sub> > n<sub>HF</sub></caption><thead><tr><th rowspan="2">M<sup>HF</sup><sub>(3,2)</sub></th><th colspan="3">T<sub>S<sub>i</sub></sub></th><th colspan="3">T<sub>S<sub>e</sub></sub></th><th rowspan="2">n<sub>HF</sub></th><th rowspan="2">n<sub>LF</sub></th><th rowspan="2">c<sub>tot</sub></th></tr><tr><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th></tr></thead><tbody><tr><td>NIRB</td><td>0.207</td><td>[0.025;0.731]</td><td>0.999</td><td><strong>0.18</strong></td><td>[0.034;<strong>0.691</strong>]</td><td>0.998</td><td rowspan="3">90</td><td rowspan="3">115</td><td rowspan="3">100C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>0.204</strong></td><td>[0.053;<strong>0.599</strong>]</td><td>0.999</td><td>0.328</td><td>[0.080;1.443]</td><td>0.994</td></tr><tr><td>LSTM</td><td>0.672</td><td>[0.543;1.120]</td><td>0.994</td><td>0.427</td><td>[0.195;0.929]</td><td>0.988</td></tr><tr><td>NIRB</td><td>0.175</td><td>[0.025;0.623]</td><td>0.999</td><td>0.139</td><td>[0.027;0.557]</td><td>0.999</td><td rowspan="3">375</td><td rowspan="3">1437</td><td rowspan="3">500C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>0.115</strong></td><td>[0.022;<strong>0.462</strong>]</td><td>0.999</td><td><strong>0.130</strong></td><td>[0.031;0.537]</td><td>0.999</td></tr><tr><td>LSTM</td><td>0.468</td><td>[0.364;0.740]</td><td>0.997</td><td>0.212</td><td>[0.115;<strong>0.406</strong>]</td><td>0.997</td></tr><tr><td>NIRB</td><td>0.168</td><td>[0.022;0.623]</td><td>0.999</td><td><strong>0.129</strong></td><td>[0.023;0.459]</td><td>0.999</td><td rowspan="3">750</td><td rowspan="3">2875</td><td rowspan="3">1000C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>0.113</strong></td><td>[0.019;<strong>0.476</strong>]</td><td>0.999</td><td>0.171</td><td>[0.032;0.679]</td><td>0.998</td></tr><tr><td>LSTM</td><td>0.531</td><td>[0.446;0.620]</td><td>0.997</td><td>0.197</td><td>[0.092;<strong>0.425</strong>]</td><td>0.998</td></tr></tbody></table>
<table><caption>Table G.11: Prediction quality of $T_{S_i}$ and $T_{S_c}$ in °C from the three predictors using the R<sub>5</sub>C<sub>4</sub> model (LF<sub>RC</sub>) as LF model, based on the indicators MAE*, Q<sub>95%</sub> and R<sup>2*</sup> calculated over a test set of size 300. c<sub>tot</sub> aggregates the costs associated with the HF evaluations.</caption><thead><tr><th rowspan="2">M<sup>HF</sup><sub>2b</sub></th><th colspan="3">T<sub>S<sub>i</sub></sub></th><th colspan="3">T<sub>S<sub>c</sub></sub></th><th rowspan="2">n<sub>HF</sub></th><th rowspan="2">c<sub>tot</sub></th></tr><tr><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th><th>MAE*</th><th>Q<sub>95%</sub></th><th>R<sup>2*</sup></th></tr></thead><tbody><tr><td>NIRB</td><td>0.194</td><td>[0.031;0.687]</td><td>0.998</td><td><strong>0.175</strong></td><td>[0.038;<strong>0.676</strong>]</td><td>0.998</td><td rowspan="3">100</td><td rowspan="3">100C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>0.182</strong></td><td>[0.045;<strong>0.538</strong>]</td><td>0.999</td><td>0.332</td><td>[0.076;1.363]</td><td>0.994</td></tr><tr><td>LSTM</td><td>0.630</td><td>[0.484;1.131]</td><td>0.994</td><td>0.458</td><td>[0.206;1.033]</td><td>0.985</td></tr><tr><td>NIRB</td><td>0.173</td><td>[0.024;0.619]</td><td>0.999</td><td><strong>0.144</strong></td><td>[0.027;0.536]</td><td>0.998</td><td rowspan="3">500</td><td rowspan="3">500C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>0.129</strong></td><td>[0.025;<strong>0.513</strong>]</td><td>0.999</td><td>0.257</td><td>[0.047;1.099]</td><td>0.996</td></tr><tr><td>LSTM</td><td>0.474</td><td>[0.355;0.897]</td><td>0.997</td><td>0.210</td><td>[0.093;<strong>0.412</strong>]</td><td>0.997</td></tr><tr><td>NIRB</td><td>0.165</td><td>[0.021;0.581]</td><td>0.999</td><td><strong>0.131</strong></td><td>[0.024;<strong>0.449</strong>]</td><td>0.999</td><td rowspan="3">1000</td><td rowspan="3">1000C<sub>HF</sub></td></tr><tr><td>GPR</td><td><strong>0.105</strong></td><td>[0.019;<strong>0.455</strong>]</td><td>0.999</td><td>0.214</td><td>[0.045;0.896]</td><td>0.997</td></tr><tr><td>LSTM</td><td>0.389</td><td>[0.320;0.547]</td><td>0.998</td><td>0.254</td><td>[0.130;0.509]</td><td>0.996</td></tr></tbody></table>