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

HAL
open science

ARegularized Real-Time Integrator for Data-Driven
Control of Heating Channels

Chady Ghnatios, Victor Champaney, Angelo Pasquale, Francisco Chinesta

► To cite this version:

Chady Ghnatios, Victor Champaney, Angelo Pasquale, Francisco Chinesta. ARegularized Real-Time Integrator for Data-Driven Control of Heating Channels. Computation, 2022, 10 (10), pp.176. 10.3390/computation10100176 . hal-05335609

HAL Id: hal-05335609

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

Submitted on 5 Dec 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# A Regularized Real-Time Integrator for Data-Driven Control of Heating Channels

Chady Ghnatios<fnref n="1" />,*, Victor Champaney<fnref n="2" />, Angelo Pasquale<fnref n="2" />,<fnref n="3" /> and Francisco Chinesta<fnref n="4" />,<fnref n="5" />

**Abstract:** In many contexts of scientific computing and engineering science, phenomena are monitored over time and data are collected as time-series. Plenty of algorithms have been proposed in the field of time-series data mining, many of them based on deep learning techniques. High-fidelity simulations of complex scenarios are truly computationally expensive and a real-time monitoring and control could be efficiently achieved by the use of artificial intelligence. In this work we build accurate data-driven models of a two-phase transient flow in a heated channel, as usually encountered in heat exchangers. The proposed methods combine several artificial neural networks architectures, involving standard and transposed deep convolutions. In particular, a very accurate real-time integrator of the system has been developed.

**Keywords:** machine learning; artificial intelligence; integrator; data-driven modeling

## 1. Introduction

Most phenomena in engineering are described by sequential data. In the context of data-driven simulation and control [1], time series classification and forecasting play therefore a fundamental role. A huge variety of applications, such as financial mathematics, language and speech recognition, weather forecasting, physical and chemical unsteady systems, are extensively discussed in the literature [2–5].

In the last decade, data-driven approximations of the Koopman operator have successfully been applied to build reduced-models of complex dynamical systems, allowing their identification and control [6–9]. Several approaches are based on combining algorithms derived from the dynamic mode decomposition–DMD—with deep learning–DL–techniques, as discussed in [10,11] and references therein.

Many algorithms have been applied in artificial intelligence–AI—for the estimation of output time series from input time series. Among them, convolution and recurrent neural networks–CNNs, RNNs– [12–14], time delay neural networks–TDNNs– [15–17], nonlinear autoregressive neural networks–NARNETs– [18,19] and nonlinear autoregressive neural networks with exogenous inputs–NARXNETs– [20–24].

In this work, several surrogates based on standard and transposed convolution neural networks–CNNs– [25–27] are applied to a thermohydraulic problem in the context of heat exchangers involving change of phase [28–30]. The problem will be described in detail later. It basically consists in fixing the input heat power $u(t)$ injected into a two-phase flow and monitoring some quantities $y(t)$ at the outlet of a pipeline in which the gas-liquid mixture is confined and flows.

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

### check for updates

**Citation:** Ghnatios, C.; Champaney V.; Pasquale, A.; Chinesta, F. A
Regularized Real-Time Integrator for Data-Driven Control of Heating Channels. *Computation* **2022**, *10*, 176.
https://doi.org/10.3390/computation10100176

Academic Editor: Ali Cemal Benim

Received: 11 September 2022

Accepted: 27 September 2022

Published: 5 October 2022

Corrected: 11 March 2024

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2022 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ Department of Mechanical Engineering, Notre Dame University-Louaize, Zouk Mosbeh P.O. Box 72, Lebanon</footnote>

<footnote n="2">² PIMM Laboratory, Arts et Métiers Institute of Technology, 151 Boulevard de l'Hôpital, 75013 Paris, France</footnote>

<footnote n="3">³ LAMPA Laboratory, Arts et Métiers Institute of Technology, 2 Boulevard du Ronceray, BP 93525, CEDEX 01, 49035 Angers, France</footnote>

<footnote n="4">⁴ PIMM Laboratory, Arts et Métiers Institute of Technology, CNRS, Cnam, HESAM Université, 151 Boulevard de l'Hôpital, 75013 Paris, France</footnote>

<footnote n="5">⁵ ESI Group, 3 bis Rue Saarien, 94528 Rungis, France</footnote>

<footnote>* Correspondence: cghnatios@ndu.edu.lb</footnote>
Mathematically speaking, let us consider an input and output time series, respectively **u** and **y**, tracking $N_t$ time instants. The problem is the estimation of the non-linear transfer function $h$ such that $\mathbf{y} = h(\mathbf{u})$. Learning the function $h$ allows the real-time integration of a newly defined input to know the impact on some outputs of interest.

A first approach is based on exploiting a parametrization of the input function $u(t)$, through a few values {$U_1, U_2, \dots, U_K$}. In this case, the modeling framework is non-incremental, since the whole output time series is computed from the inputs {$U_1, U_2, \dots, U_K$}, that is:

$$ \mathbf{y} = h(U_1, U_2, \dots, U_K). \qquad (1) $$

The estimation of $h$ will be done through the use of artificial neural networks. However, since the output consists of $N_t$ evaluations, while the input of $K < N_t$ parameters, a transposed convolution (also known as fractionally-strided convolution) is employed to reach the right dimension at the output layer. As a matter of fact, transposed convolutions increase (up sample) the spatial dimensions of intermediate feature maps, reversing down sample operations obtained by classical convolutions [25].

A second approach consists of considering an incremental modelling framework, in which the model (also known as integrator) computes the next value of the output sequence from a given number of previous inputs and outputs. Such approach has some similarities with the NARXNETs [18]. Denoting with $y_i = y(t_i)$, for $i = 1, \dots, N_t$, we can write:

$$ y_i = h(y_{i-1}, y_{i-2}, \dots, y_{i-k}, u_i, u_{i-1}, \dots, u_{i-k+1}), \qquad (2) $$

where the parameter $k \ge 0$ is a delay term, known as the process dead time. We can also note $U = y_{\text{past}} \cup u, y_{\text{past}}$ being the known $y$ outputs used to predict future values of $y$. A neural-network architecture based on a convolution approach can then be defined. Two non-incremental models (for different values of the delay coefficient $k$) will be compared, in terms of speed and accuracy.

The paper aims to build a data-driven model that is capable of reproducing the data with high fidelity, while possibly performing the integrator task over a large time frame. First of all, the model of the addressed problem and the simulation data structure are described in Section 2. Section 3 presents three different surrogates:

* A non-incremental modelling framework is presented in Section 3.1, which predicts the whole time frame in a single shot

* Two incremental alternatives are presented in Section 3.2. These alternatives used different memory size to predict the future variation in the quantities of interest.

* The incremental surrogate model with large memory is then used as an integrator to predict the variation over a large time series

Finally, Section 4 addresses the conclusion summarizing the results and highlighting the main differences among the learned models.

## 2. Problem Statement

High-fidelity numerical simulations of a thermohydraulic problem are performed using the software CATHARE for nuclear safety analysis [31–33]. The simulation consists into a two-phase (steam-water) flow moving upward against gravity into a channel heated by hot walls, as sketched in Figure 1. The model input is the time-varying thermal power $q(t)$ injected into the fluid through the walls, while the outputs are the time histories of several fluid properties at the channel outlet. The simulation covers the interval $[0, 100]$ s, with a time step $\Delta t = 0.1$ s.
Using the bold notation for vectors (sequences of 1000 elements corresponding to the time discretization), the quantities of interest are the steam fraction in the flowing system A, the liquid velocity V, the steam temperature T<sub>G</sub> and the liquid temperature T<sub>L</sub>. This work aims to produce a surrogate model h for inferring the outputs as a function of the input variable q:

$$ \begin{pmatrix} \mathbf{A} \\ \mathbf{V} \\ \mathbf{T}_G \\ \mathbf{T}_L \end{pmatrix} = h(\mathbf{q}). \qquad (3) $$

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Fluid channel with wall surface heating, controlled by the heat power density law $q(t)$. The green point represents the location of the measurement node at the outlet.</figcaption></figure>

The surrogate model is built in a non-intrusive manner, defining a Design of Experiments–DoE–on the power density and combining the outputs of the corresponding offline simulations.

A user-defined power control law allows to assign the value of the provided power density at specific time instants. Therefore, the input function $q(t)$ is defined through a few values $Q_i = q(t_i)$, normalized by the maximum power density $Q_{max} = 8.4 \times 10^8 \text{ W/m}^3$, and then linearly interpolating in-between. In particular, the time instants (in seconds) $\{0, 1, 25, 50, 75, 100\}$ are considered for the power changes. The corresponding values $Q_i$, for $i = 0, \dots, 5$, define a parametrization of the input function $q(t)$.

A DoE of 23 data points is established using different combinations of the values $Q_i$, as listed in Table 1. For the rest of this work, the last 4 datasets are considered in the testing set and therefore not fed into the training loop of the built surrogate model. Note that some of the values are fixed ($Q_0 = 0$, $Q_1 = Q_5 = 1$) thus the model input parameters are only $Q_2, Q_3$ and $Q_4$. Figure 2 shows an example of the heat power density $q(t)$ defined from values $Q_i$.
<table><caption>**Table 1.** Available plan of experiment (DoE) and sequences of control points values, as a fraction of $Q_{max} = 8.4 \times 10^8$ W/m³.</caption><thead><tr><th>Sequence Number</th><th>Q<sub>0</sub></th><th>Q<sub>1</sub></th><th>Q<sub>2</sub></th><th>Q<sub>3</sub></th><th>Q<sub>4</sub></th><th>Q<sub>5</sub></th></tr></thead><tbody><tr><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td></tr><tr><td>2</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td></tr><tr><td>3</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td></tr><tr><td>4</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr><tr><td>5</td><td>0</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td></tr><tr><td>6</td><td>0</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td></tr><tr><td>7</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td></tr><tr><td>8</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>9</td><td>0</td><td>1</td><td>0.19</td><td>0.07</td><td>0.9</td><td>1</td></tr><tr><td>10</td><td>0</td><td>1</td><td>0.06</td><td>0.8</td><td>0.28</td><td>1</td></tr><tr><td>11</td><td>0</td><td>1</td><td>0.91</td><td>0.76</td><td>0.58</td><td>1</td></tr><tr><td>12</td><td>0</td><td>1</td><td>0.69</td><td>0.19</td><td>1</td><td>1</td></tr><tr><td>13</td><td>0</td><td>1</td><td>0.12</td><td>0.57</td><td>0.38</td><td>1</td></tr><tr><td>14</td><td>0</td><td>1</td><td>0.2</td><td>0.67</td><td>0.02</td><td>1</td></tr><tr><td>15</td><td>0</td><td>1</td><td>0.99</td><td>0.92</td><td>0.47</td><td>1</td></tr><tr><td>16</td><td>0</td><td>1</td><td>0.55</td><td>0.94</td><td>0.74</td><td>1</td></tr><tr><td>17</td><td>0</td><td>1</td><td>0.47</td><td>0.24</td><td>0.83</td><td>1</td></tr><tr><td>18</td><td>0</td><td>1</td><td>0.76</td><td>0.36</td><td>0.18</td><td>1</td></tr><tr><td>19</td><td>0</td><td>1</td><td>0.82</td><td>0.49</td><td>0.4</td><td>1</td></tr><tr><td>20</td><td>0</td><td>1</td><td>0.66</td><td>0.66</td><td>0.72</td><td>1</td></tr><tr><td>21</td><td>0</td><td>1</td><td>0.31</td><td>0.43</td><td>0.6</td><td>1</td></tr><tr><td>22</td><td>0</td><td>1</td><td>0.37</td><td>0.31</td><td>0.11</td><td>1</td></tr><tr><td>23</td><td>0</td><td>1</td><td>0.44</td><td>0.05</td><td>0.2</td><td>1</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Example of provided heat power density, as a fraction of $Q_{max} = 8.4 \times 10^8$ W/m³ (this corresponds to the run 9 of the DoE).</figcaption></figure>

A sample set of results, of the available datasets number 1, 2, 3 and 4 of Table 1, is shown in Figure 3 for illustration purposes.
<figure><img src="image_7.png" /><figcaption>Figure 3. Sample results for the datasets number 1 to 4 as shown in Table 1.</figcaption></figure>

### 3. Surrogate Modeling of the Data

Several models are tested and in what follows we explain the results of three of them. Before any training, the variables are normalized by subtracting their corresponding means and dividing by the variation amplitude $\bullet_{max} - \bullet_{min}$.

#### 3.1. A Non-Incremental Approach

In this section, we present a non-incremental modelling framework, aiming to reproduce all the measured 1000 steps time sequence while considering only the three variables $Q_2, Q_3$ and $Q_4$. Therefore, the first model $\mathcal{H}_1$ aims to create an output of the form:

$$ \left( \begin{array}{c} \hat{\mathbf{A}} \\ \hat{\mathbf{V}} \\ \hat{\mathbf{T}}_G \\ \hat{\mathbf{T}}_L \end{array} \right) = \mathcal{H}_1(Q_2, Q_3, Q_4), \qquad (4) $$

with $\bullet$ being the approximation of the data, output of the surrogate model. In this section, since we are using a non-incremental approach, each dataset is nothing but a single data-point, for which we evaluate a single output. We use a neural network based on a transposed convolution approach [25–27]. The network starts with an input of 3 values, for instance $(Q_2, Q_3, Q_4)$, and propagate the variation of the input variables into the required output. The architecture of the used network is illustrated in Table 2 and schematically represented through the diagram in Figure 4. The first layer in the designed neural network transforms the input into a dimension suitable for beginning a convolution transpose. The second layer reshapes the inputs into a simulated 2D picture of dimension $8 \times 4$, with 240 filters. The subsequent layers reduce the number of filtering while increasing the size of the first dimension of the tensor, through multiplying by 5 on each layer, reaching therefore a dimension of $1000 \times 4$ on layer 6, the exact dimension of the required 4 variables output.
<table><caption>**Table 2.** Structure of the deep convolution transpose neural network used for the fitting of $H_1$. tanh stands for the hyperbolic tangent activation function, selu stands for the scaled exponential linear unit [34], while linear stands for linear activation or no activation function. CT stands for “convolution transpose”.</caption><thead><tr><th>Layers</th><th>Shape</th><th>Activation</th></tr></thead><tbody><tr><td>1</td><td>Fully connected dense layer with 7680 neurons</td><td>tanh</td></tr><tr><td>2</td><td>Reshape layer into a 3D Tensor of shape (8 × 4 × 240)</td><td>No activation</td></tr><tr><td>3</td><td>2D CT, 120 filters, kernel (5 × 1), strides (5 × 1), “same” padding</td><td>selu</td></tr><tr><td>4</td><td>2D CT, 60 filters, kernel (5 × 1), strides (5 × 1), “same” padding</td><td>selu</td></tr><tr><td>5</td><td>2D CT, 30 filters, kernel (5 × 1), strides (5 × 1), “same” padding</td><td>selu</td></tr><tr><td>6</td><td>2D CT, 1 filter, kernel (1 × 1), strides (1 × 1), “same” padding</td><td>linear</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Graphical diagram of the deep convolution neural network used for the fitting of $H_1$. All convolutions employ “same” padding (also called “half” padding). If the input and output of a transpose convolution have shapes $I_1 \times I_2 \times I_3$ and $O_1 \times O_2 \times O_3$, respectively, then the dimensions are obtained as $O_i = s_i I_i$, where $s_i$ denotes the stride along the dimension $i$, for $i = 1, 2$. The dimension $O_3$ is determined by the number of filters.</figcaption></figure>

The neural network is trained over 5000 epochs, with a batch size of 10 datasets, knowing that we have only 19 available datasets. The selected loss function is the mean square error, without regularization, and the optimization algorithm is Adam [35] with an initial learning rate $\alpha = 10^{-4}$. The results for a selected training dataset are shown in Figure 5, whereas the results for the selected testing datasets are illustrated for datasets 20 and 22 in Figures 6 and 7, respectively.
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Sample results of $H_1$ surrogate model for the training datasets number 10, as shown in Table 1.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Sample results of $H_1$ surrogate model for the testing datasets number 20, as shown in Table 1.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Sample results of $H_1$ surrogate model for the testing datasets number 22, as shown in Table 1. Note the representation scale in panel (c).</figcaption></figure>

The relative errors generated by the non-incremental approach are illustrated in Table 3. Since the output **A** contains several zero values, the relative error $\mathcal{E}_A$ is normalized with respect to the maximum value of **A**:

$$ \mathcal{E}_A = \frac{\|\mathbf{A} - \hat{\mathbf{A}}\|_2}{\|\mathbf{A}\|_\infty}, \qquad (5) $$

with **A** being the measured data and $\hat{\mathbf{A}}$ the surrogate model output. The normalization is standard for any other variable **v**:

$$ \mathcal{E}_v = \|\mathbf{v} - \hat{\mathbf{v}}\| \odot \|\mathbf{v}\|_2, \qquad (6) $$

where, similarly, **v** denotes the measured data, $\hat{\mathbf{v}}$ the surrogate model output and the symbol $\odot$ the Hadamard division operator. The results are excellent in either training or testing sets.

<table><caption>**Table 3.** Mean relative errors of the trained non-incremental model $H_1$.</caption><thead><tr><th>Model <i>H</i><sub>1</sub></th><th><i>E</i><sub>A</sub></th><th><i>E</i><sub>V</sub></th><th><i>E</i><sub>T<sub>G</sub></sub></th><th><i>E</i><sub>T<sub>L</sub></sub></th></tr></thead><tbody><tr><td>Training sets</td><td>0.058%</td><td>0.735%</td><td>1.7 × 10<sup>-5</sup>%</td><td>0.037%</td></tr><tr><td>Testing sets</td><td>0.14%</td><td>1.61%</td><td>9.0 × 10<sup>-5</sup>%</td><td>0.265%</td></tr></tbody></table>

### 3.2. Incremental Approaches

Using an incremental approach consists in learning the integrator that allows to update the state of the system at each time step. The model is built by first decomposing the datasets into a sequence chain of length $t_s$. Increasing the chain size $t_s$, by considering therefore a higher number of time steps in the variation history, will definitely improve
the results, by jeopardizing the computation time required for the prediction of the data point. Two incremental models are tested, the first one considers $t_s = 4$, and is an excellent predictor of the next output of the considered sequences but fails as dynamical integrator as will be discussed later, whereas the second model uses $t_s = 20$, with a regularization, for improving stability, that performs correctly as integrator.

### 3.2.1. A Fast Predictor

In this section we explain the possibility of using a fast algorithm based on convolution layers, to predict the next outcome of the sequence at hand. The model will take as an input the variables ($\mathbf{A}$, $\mathbf{V}$, $\mathbf{T}_G$, $\mathbf{T}_L$) at time steps $i-1$, $i-2$, $i-3$ and $i-4$, as well as the input $q$ at time steps $i$, $i-1$, $i-2$ and $i-3$, to predict the quantities of interest at time step $i$. For instance, we can write the second trained surrogate model $\mathcal{H}_2$ as follows:

$$ \begin{pmatrix} \hat{A}_i \\ \hat{V}_i \\ \hat{T}_{Gi} \\ \hat{T}_{Li} \end{pmatrix} = \mathcal{H}_2 \left( \begin{pmatrix} A_{i-1} \\ A_{i-2} \\ A_{i-3} \\ A_{i-4} \end{pmatrix}, \begin{pmatrix} V_{i-1} \\ V_{i-2} \\ V_{i-3} \\ V_{i-4} \end{pmatrix}, \begin{pmatrix} T_{Gi-1} \\ T_{Gi-2} \\ T_{Gi-3} \\ T_{Gi-4} \end{pmatrix}, \begin{pmatrix} T_{Li-1} \\ T_{Li-2} \\ T_{Li-3} \\ T_{Li-4} \end{pmatrix}, \begin{pmatrix} q_i \\ q_{i-1} \\ q_{i-2} \\ q_{i-3} \end{pmatrix} \right). \quad (7) $$

The considered network to approximate the surrogate model $\mathcal{H}_2$ is summarized in Table 4 and schematized in Figure 8. The input of the network is a 3D tensor of dimension $(4 \times 5 \times 1)$, for 4 time steps and 5 variables. The selected optimization algorithm is Adam with the loss function being the mean square errors. The available data after reformulation into a sequence of 4 steps has the size of $19 \times 996$ inputs for the training sets and $4 \times 996$ for the testing set. A batch size of 996 is selected with 1500 epochs for the training. The initial learning rate was set to $\alpha = 10^{-4}$.

<table><caption>**Table 4.** Structure of the deep convolution neural network used for the fitting of $\mathcal{H}_2$. tanh stands for the hyperbolic tangent activation function, *selu* stands for the scaled exponential linear unit [34], while *linear* stands for linear activation or no activation function.</caption><thead><tr><th>Layers</th><th>Shape</th><th>Activation</th></tr></thead><tbody><tr><td>1</td><td>2D convolution, 60 filters, kernel (2 × 1), strides (2 × 1), “same” padding</td><td>selu</td></tr><tr><td>2</td><td>2D convolution, 120 filters, kernel (2 × 1), strides (1 × 1), “same” padding</td><td>selu</td></tr><tr><td>3</td><td>A flatten layer, reshapes all inputs into a single vector</td><td>no activation</td></tr><tr><td>4</td><td>Fully connected dense layer with 250 neurons</td><td>selu</td></tr><tr><td>5</td><td>Fully connected dense layer with 125 neurons</td><td>selu</td></tr><tr><td>6</td><td>Fully connected dense layer with 70 neurons</td><td>selu</td></tr><tr><td>7</td><td>Fully connected dense layer with 25 neurons</td><td>tanh</td></tr><tr><td>8</td><td>Fully connected dense layer with 4 neurons</td><td>linear</td></tr></tbody></table>

The results for a selected training dataset are shown in Figure 9, whereas the results for the selected testing datasets are illustrated for datasets 20 and 22 in Figures 10 and 11, respectively. The relative errors are shown in Table 5, computed using Equations (5) and (6), as described in Section 3.1.

Model $\mathcal{H}_2$ shows an excellent ability to reproduce the datasets when the input is always taken from the physical dataset values and not from the output of the predictor $\mathcal{H}_2$ itself. The errors increase dramatically when the input of $\mathcal{H}_2$ is taken from its previous output, yielding an impossibility to use it as an integrator of the dynamical system.

On the other hand, the performance of $\mathcal{H}_2$ is comparable to the performance of $\mathcal{H}_1$, when comparing the relative errors of both models. From a comparison of the resulting graphical illustrations, $\mathcal{H}_2$ seems to perform better on the testing sets.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Graphical diagram of the deep convolution neural network used for the fitting of $H_2$. All convolutions employ "same" padding (also called "half" padding). If the input and output of a convolution have shapes $I_1 \times I_2 \times I_3$ and $O_1 \times O_2 \times O_3$, respectively, then the dimensions are obtained as $O_i = \lfloor I_i / s_i \rfloor$, where $s_i$ denotes the stride along the dimension $i$, for $i = 1, 2$. The dimension $O_3$ is determined by the number of filters.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Sample results of $H_2$ surrogate model for the training datasets number 10, as shown in Table 1. (a) A; (b) V; (c) $T_G$; (d) $T_L$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 10. Sample results of $H_2$ surrogate model for the testing datasets number 20, as shown in Table 1. (a) A; (b) V; (c) $T_G$; (d) $T_L$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 11. Sample results of $H_2$ surrogate model for the testing datasets number 22, as shown in Table 1.</figcaption></figure>
<table><caption>**Table 5.** Mean relative errors of the trained incremental model $\mathcal{H}_2$ with $t_s = 4$. All convolutions employ “same” padding, meaning that the input is half padded and that the filter is applied to all the elements of the input.</caption><thead><tr><th>Model <i>H</i><sub>2</sub></th><th><i>E</i><sub>A</sub></th><th><i>E</i><sub>V</sub></th><th><i>E</i><sub>T<sub>G</sub></sub></th><th><i>E</i><sub>T<sub>L</sub></sub></th></tr></thead><tbody><tr><td>Training sets</td><td>0.11%</td><td>0.81%</td><td>1.91 × 10<sup>-5%</sup></td><td>0.25%</td></tr><tr><td>Testing sets</td><td>0.039%</td><td>0.898%</td><td>1.94 × 10<sup>-5%</sup></td><td>0.284%</td></tr></tbody></table>

### 3.2.2. A Larger Sequence Time Integrator

In this section we consider again the rationale followed in Section 3.2.1, with an algorithm based on convolution layers, to predict the next outcome of the sequence at hand. The model will take as an input the variables available at the last 20 time steps. According to the previously introduced notation, $t_s = 20$. The built model $\mathcal{H}_3$ will predict the quantities of interest at the subsequent time step. For instance, we can write the trained surrogate model $\mathcal{H}_3$ as shown below:

$$ \begin{pmatrix} \hat{A}_i \\ \hat{V}_i \\ \hat{T}_{Gi} \\ \hat{T}_{Li} \end{pmatrix} = \mathcal{H}_3 \left( \begin{pmatrix} A_{i-1} \\ \vdots \\ A_{i-20} \end{pmatrix}, \begin{pmatrix} V_{i-1} \\ \vdots \\ V_{i-20} \end{pmatrix}, \begin{pmatrix} T_{Gi-1} \\ \vdots \\ T_{Gi-20} \end{pmatrix}, \begin{pmatrix} T_{Li-1} \\ \vdots \\ T_{Li-20} \end{pmatrix}, \begin{pmatrix} q_i \\ \vdots \\ q_{i-19} \end{pmatrix} \right). \quad (8) $$

The surrogate model $\mathcal{H}_3$ is built using the network illustrated in Table 6 and represented through the diagram in Figure 12. The input of the network is a 3D tensor of dimension $(20 \times 5 \times 1)$, for 20 time steps and 5 variables as inputs. The selected optimization algorithm is Adam with the loss function being the mean square errors, modified to add a regularization similar to the DMD regularization [36,37], to perform integration in time without increasing the error significantly. For instance, the loss function $\mathcal{J}$ in this section reads:

$$ \mathcal{J} = \sum_{j=1}^{19} (\|\mathbf{A}_j - \hat{\mathbf{A}}_j\| + \|\mathbf{V}_j - \hat{\mathbf{V}}_j\| + \|\mathbf{T}_{Gi} - \hat{\mathbf{T}}_{Gi}\| + \|\mathbf{T}_{Li} - \hat{\mathbf{T}}_{Li}\|) + \\ + \frac{\lambda_w}{2m} \sum_{j=1}^{m} \theta_j^2 + \frac{\lambda_b}{2p} \sum_{j=1}^{p} \kappa_j^2 \quad (9) $$

where $\theta_j$ are all the weights of the last layer of the network described in Table 6 (the linear layer), and $\kappa_j$ its biases. $\lambda_w$ and $\lambda_b$ are regularization parameters. The available data after reformulation into a sequence of 20 steps has the size of 19 × 980 inputs for the training sets and 4 × 980 for the testing set. A batch size of 980 is selected with 1500 epochs for the training. The initial learning rate was set to $\alpha = 10^{-4}$. The input shape has therefore of a size of 20 × 5 × 1. The results are shown for a $\lambda_w = \lambda_b = 10^{-6}$.

<table><caption>**Table 6.** Structure of the deep convolution neural network used for the fitting of $\mathcal{H}_3$. tanh stands for the hyperbolic tangent activation function, *selu* stands for the scaled exponential linear unit [34], while *linear* stands for linear activation or no activation function.</caption><thead><tr><th>Layers</th><th>Shape</th><th>Activation</th></tr></thead><tbody><tr><td>1</td><td>2D convolution, 30 filters, kernel (4 × 5), strides (2 × 1), “same” padding</td><td>selu</td></tr><tr><td>2</td><td>2D convolution, 60 filters, kernel (2 × 1), strides (2 × 1), “same” padding</td><td>selu</td></tr><tr><td>3</td><td>2D convolution, 120 filters, kernel (5 × 1), strides (1 × 1), “same” padding</td><td>selu</td></tr><tr><td>4</td><td>A flatten layer, reshapes all inputs into a single vector</td><td>no activation</td></tr></tbody></table>
<table><caption>**Table 6.** *Cont.*</caption><thead><tr><td>Layers</td><td>Shape</td><td>Activation</td></tr></thead><tbody><tr><td>5</td><td>Fully connected dense layer with 1200 neurons</td><td>tanh</td></tr><tr><td>6</td><td>Fully connected dense layer with 125 neurons</td><td>selu</td></tr><tr><td>7</td><td>Fully connected dense layer with 70 neurons</td><td>selu</td></tr><tr><td>8</td><td>Fully connected dense layer with 25 neurons</td><td>selu</td></tr><tr><td>9</td><td>Fully connected dense layer with 4 neurons</td><td>linear</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>**Figure 12.** Graphical diagram of the deep convolution neural network used for the fitting of $H_3$. All convolutions employ “same” padding (also called “half” padding). If the input and output of a convolution have shapes $I_1 \times I_2 \times I_3$ and $O_1 \times O_2 \times O_3$, respectively, then the dimensions are obtained as $O_i = \lfloor I_i / s_i \rfloor$, where $s_i$ denotes the stride along the dimension $i$, for $i = 1, 2$. The dimension $O_3$ is determined by the number of filters.</figcaption></figure>

The results using physical inputs for a selected training dataset 10 are shown in Figure 13, whereas the results for the selected testing datasets 20 and 22 in Figures 14 and 15, respectively. The relative errors are shown in Table 7, computed using Equations (5) and (6), as illustrated in Section 3.1.

<figure><img src="image_17.png" /><figcaption>**Figure 13.** *Cont.*</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 13.** Sample results of $\mathcal{H}_3$ surrogate model for the training dataset number 10, as shown in Table 1.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 14.** Sample results of $\mathcal{H}_3$ surrogate model for the testing datasets number 20, as shown in Table 1.</figcaption></figure>

<table><caption>**Table 7.** Mean relative errors of the trained incremental model $\mathcal{H}_3$ with $t_s = 20$.</caption><thead><tr><th>Model <i>H</i><sub>3</sub></th><th><i>ε</i><sub>A</sub></th><th><i>ε</i><sub>V</sub></th><th><i>ε</i><sub>T<sub>G</sub></sub></th><th><i>ε</i><sub>T<sub>L</sub></sub></th></tr></thead><tbody><tr><td>Training sets</td><td>0.058%</td><td>1.21%</td><td>9.18 × 10<sup>-6</sup>%</td><td>0.039%</td></tr><tr><td>Testing sets</td><td>0.012%</td><td>1.4%</td><td>9.11 × 10<sup>-6</sup>%</td><td>0.04%</td></tr></tbody></table>

### 3.2.3. Integrator Results

The neural network described in Table 6, complemented with the regularized loss given in Equation (9), can be used to integrate in time the system behavior. For instance, we can write
$$ \begin{pmatrix} \hat{A}_i \\ \hat{V}_i \\ \hat{T}_{Gi} \\ \hat{T}_{Li} \end{pmatrix} = \mathcal{H}_3 \left( \begin{pmatrix} \hat{A}_{i-1} \\ \vdots \\ \hat{A}_{i-20} \end{pmatrix}, \begin{pmatrix} \hat{V}_{i-1} \\ \vdots \\ \hat{V}_{i-20} \end{pmatrix}, \begin{pmatrix} \hat{T}_{Gi-1} \\ \vdots \\ \hat{T}_{Gi-20} \end{pmatrix}, \begin{pmatrix} \hat{T}_{Li-1} \\ \vdots \\ \hat{T}_{Li-20} \end{pmatrix}, \begin{pmatrix} q_i \\ \vdots \\ q_{i-19} \end{pmatrix} \right). \quad (10) $$

<figure><img src="image_20.png" /><figcaption>Figure 15. Sample results of $\mathcal{H}_3$ surrogate model for the testing datasets number 22, as shown in Table 1.</figcaption></figure>

The integrator is stable enough to compile all the time steps. Figure 16 shows the integrator results in a training case, whereas Figures 17 and 18 show the integrator results on two testing cases. The integrator errors are given in Table 8.

<table><caption>Table 8. Mean relative errors of the trained incremental model $\mathcal{H}_3$ with $t_s = 20$, when performing an integrator operation as per Equation (10).</caption><thead><tr><th>Model <i>H</i><sub>3</sub></th><th><i>E</i><sub>A</sub></th><th><i>E</i><sub>V</sub></th><th><i>E</i><sub>T<sub>G</sub></sub></th><th><i>E</i><sub>T<sub>L</sub></sub></th></tr></thead><tbody><tr><td>Training sets</td><td>1.7%</td><td>1.54%</td><td>2.07 × 10<sup>-5%</sup></td><td>0.091%</td></tr><tr><td>Testing sets</td><td>0.44%</td><td>1.05%</td><td>2.9 × 10<sup>-5%</sup></td><td>0.045%</td></tr></tbody></table>

<figure><img src="image_21.png" /><figcaption>Figure 16. Cont.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 16.** Sample integrator results of $H_3$ surrogate model when integrating the training datasets number 10, as shown in Table 1.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 17.** Sample integrator results of $H_3$ surrogate model when integrating the testing datasets number 20, as shown in Table 1.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 18.** *Cont.*</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure 18.** Sample integrator results of $H_3$ surrogate model when integrating the testing datasets number 22, as shown in Table 1.</figcaption></figure>

It can be clearly noticed that the integrator performs extremely well on the selected datasets.

**4. Conclusions**

In this work several prediction models for output time series from input time series have been compared. Without loss of generality, a coupled problem in thermohydraulics has been considered, but the procedure generalizes to many other phenomena and fields. Moreover, instead of simulation-based data, experimental data can be used.

The models are based on different deep neural networks architectures. If a parametrization of the input function is available, a non-incremental approach involving transposed convolution neural networks allows the prediction of the whole sequence, as described and discussed in Section 3.1. The non-incremental approach is based on convolution transpose layers, which takes a low amount of data as an input and increase the size of the output to reach a desired dimension. These layers appears as a suitable approach in this situation where only few inputs are required to build the whole sequence of time measurements.

The incremental models (dynamical system integrators) presented in Section 3.2 are based on regressing the variable on its own lagged (i.e., past) values together with the ones of the input, through standard deep convolution neural networks. Therefore, the prediction occurs step-by-step. The difference between the two incremental models stands in the delay (or lag) factor $t_s$. A small value of $t_s$ allows really fast predictions. However, the model can be unstable when its input is taken from the previous model outputs (this occurs, for instance, when $t_s = 4$). In such case, the surrogate cannot be used as an integrator of the system. Results are thoroughly improved for $t_s = 20$, where also a regularization term was introduced in the loss function. The incremental approach is based on convolution layers. In fact, these layers will take a large input and reduce the input size to extract relevant information. The filters are designed with a large dimension along the time axis, to allow easier extraction of the inputs' time variation.

The non-incremental approach can be clearly faster to evaluate the output of the solution, however, the integrator clearly outperforms the non-incremental approach. Moreover, the integrator can be used on novel datasets, and potentially outside of the simulated time frame. The initialization requires only $t_s$ time-steps to be simulated, and the remaining dynamics is predicted in real time, allowing huge computational gains. However, the non-incremental approach can't be used on a different time frame, out of the simulated region, as the output is hardly coded to predict only $N_t = 1000$ time steps.

**Author Contributions:** Conceptualization, C.G. and F.C.; methodology, C.G and A.P.; software, C.G.; formal analysis, C.G.; investigation, V.C. and A.P.; data curation, V.C. and A.P.; writing—original draft preparation, C.G.; writing—review and editing, V.C. and A.P.; supervision, F.C.; project
administration, F.C.; funding acquisition, F.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Data are the property of the project's consortium and cannot be made available.

**Acknowledgments:** The authors greatly acknowledge the support of Anoop Ebey Thomas from ESI Group and all the partners of the PSPC Réacteur Numérique project.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Markovsky, I.; Rapisarda, P. Data-driven simulation and control. *Int. J. Control.* **2008**, *81*, 1946–1959. [CrossRef]

2. Box, G.; Jenkins, G.M. *Time Series Analysis: Forecasting and Control*; Holden-Day: San Francisco, CA, USA, 1976.

3. Pollock, D.; Green, R.; Nguyen, T. *A Handbook of Time-Series Analysis, Signal Processing and Dynamics*; Signal Processing and Its Applications; Academic Press: Cambridge, MA, USA, 1999.

4. Kantz, H.; Schreiber, T. *Nonlinear Time Series Analysis*; Cambridge University Press: Cambridge, UK, 2003.

5. Palit, A.K.; Popovic, D. *Computational Intelligence in Time Series Forecasting: Theory and Engineering Applications (Advances in Industrial Control)*; Springer: Berlin/Heidelberg, Germany, 2005.

6. Lian, Y.; Wang, R.; Jones, C.N. Koopman based data-driven predictive control. *arXiv* **2021**, arXiv:2102.05122.

7. Klus, S.; Nüske, F.; Peitz, S.; Niemann, J.H.; Clementi, C.; Schütte, C. Data-driven approximation of the Koopman generator: Model reduction, system identification, and control. *Phys. D Nonlinear Phenom.* **2020**, *406*, 132416. [CrossRef]

8. Otto, S.; Rowley, C. Koopman Operators for Estimation and Control of Dynamical Systems. *Annu. Rev. Control. Robot. Auton. Syst.* **2021**, *4*, 59–87. [CrossRef]

9. Maksakov, A.; Palis, S. Koopman-based data-driven control for continuous fluidized bed spray granulation with screen-mill-cycle. *J. Process. Control.* **2021**, *103*, 48–54. [CrossRef]

10. Williams, M.O.; Kevrekidis, I.G.; Rowley, C.W. A Data-Driven Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition. *J. Nonlinear Sci.* **2015**, *25*, 1307–1346. [CrossRef]

11. Zhang, C.; Zuazua, E. A quantitative analysis of Koopman operator methods for system identification and predictions. 2021, in press.

12. Tijskens, A.; Roels, S.; Janssen, H. Neural networks for metamodelling the hygrothermal behaviour of building components. *Build. Environ.* **2019**, *162*, 106282. [CrossRef]

13. Connor, J.; Martin, R.; Atlas, L. Recurrent neural networks and robust time series prediction. *IEEE Trans. Neural Netw.* **1994**, *5*, 240–254. [CrossRef] [PubMed]

14. Mandic, D.; Mandic, D.; Chambers, J. *Recurrent Neural Networks for Prediction: Learning Algorithms, Architectures and Stability*; Adaptive and Cognitive Dynamic Systems: Signal Processing; Wiley: Hoboken, NJ, USA, 2001.

15. Waibel, A.; Hanazawa, T.; Hinton, G.; Shikano, K.; Lang, K. Phoneme recognition using time-delay neural networks. *IEEE Trans. Acoust. Speech, Signal Process.* **1989**, *37*, 328–339. [CrossRef]

16. Peddinti, V.; Povey, D.; Khudanpur, S. A time delay neural network architecture for efficient modeling of long temporal contexts. *Proc. Interspeech* **2015**, 3214–3218. [CrossRef]

17. Wöhler, C.; Anlauf, J. Real-time object recognition on image sequences with the adaptable time delay neural network algorithm—Applications for autonomous vehicles. *Image Vis. Comput.* **2001**, *19*, 593–618. [CrossRef]

18. Pinheiro, E.; Bosco, J.; Lima, L.; Geraldo, F. Nonlinear autoregressive neural networks for forecasting wind speed time series. *Int. J. Dev. Res.* **2021**, *10*, 40336–40343. [CrossRef]

19. Ruiz, L.; Cuéllar, M.; Calvo-Flores, M.; Pegalajar Jiménez, M.d.C. An Application of Non-Linear Autoregressive Neural Networks to Predict Energy Consumption in Public Buildings. *Energies* **2016**, *9*, 684. [CrossRef]

20. Xie, H.; Tang, H.; Liao, Y.H. Time series prediction based on NARX neural networks: An advanced approach. In *Proceedings of the 2009 International Conference on Machine Learning and Cybernetics*, Hebei, China, 12–15 July 2009; Volume 3, pp. 1275–1279. [CrossRef]

21. Diaconescu, E. The use of NARX neural networks to predict chaotic time series. *WSEAS Trans. Comput. Res.* **2008**, *3*, 182–191.

22. Guzman, S.M.; Paz, J.O.; Tagert, M.L.M. The Use of NARX Neural Networks to Forecast Daily Groundwater Levels. *Water Resour. Manag.* **2017**, *31*, 1591–1603. [CrossRef]

23. Menezes, J.M.P.; Barreto, G.A. Long-term time series prediction with the NARX network: An empirical evaluation. *Neurocomputing* **2008**, *71*, 3335–3343. [CrossRef]

24. Boussaada, Z.; Curea, O.; Remaci, A.; Camblong, H.; Mrabet Bellaaj, N. A Nonlinear Autoregressive Exogenous (NARX) Neural Network Model for the Prediction of the Daily Direct Solar Radiation. *Energies* **2018**, *11*, 620. [CrossRef]

25. Dumoulin, V.; Visin, F. A guide to convolution arithmetic for deep learning. *arXiv* **2016**, arXiv:1603.07285.

26. Zeiler, M.D.; Fergus, R. Visualizing and Understanding Convolutional Networks. *arXiv* **2013**, arXiv:1311.2901.
27. Zeiler, M.D.; Taylor, G.W.; Fergus, R. Adaptive deconvolutional networks for mid and high level feature learning. In Proceedings of the 2011 International Conference on Computer Vision, Barcelona, Spain, 6–13 November 2011; pp. 2018–2025. [CrossRef]

28. Bernard, M.; Dellacherie, S.; Faccanoni, G.; Grec, B.; Penel, Y. Study of a low Mach nuclear core model for two-phase flows with phase transition I: Stiffened gas law. *ESAIM Math. Model. Numer. Anal.* **2014**, *48*, 1639–1679. [CrossRef]

29. Brown, J.; Maldonado, G. History of PWR and BWR development. *Encycl. Nucl. Energy* **2021**, 157–171. [CrossRef]

30. Brown, J.A. Pressurized Water Reactors. In *Encyclopedia of Nuclear Energy*; Greenspan, E., Ed.; Elsevier: Oxford, UK, 2021; pp. 196–213. [CrossRef]

31. Liu, Z.T.; Qin, B.K.; Xie, H.; Wang, B.H. Status and trends of thermal-hydraulic system codes for nuclear power plants with pressurized water reactors. *Yuanzineng Kexue Jishu/At. Energy Sci. Technol.* **2009**, *43*, 966–972.

32. Zhang, K. The multiscale thermal-hydraulic simulation for nuclear reactors: A classification of the coupling approaches and a review of the coupled codes. *Int. J. Energy Res.* **2020**, *44*, 3295–3315. [CrossRef]

33. Emonot, P.; Souyri, A.; Gandrille, J.; Barré, F. CATHARE-3: A new system code for thermal-hydraulics in the context of the NEPTUNE project. *Nucl. Eng. Des.* **2011**, *241*, 4476–4481.

34. Klambauer, G.; Unterthiner, T.; Mayr, A. Self-Normalizing Neural Networks. In Proceedings of the 31st Conference on Neural Information Processing Systems (NIPS 2017), Long Beach, CA, USA, 4–9 December 2017; pp. 1–102.

35. Kingma, D.; Ba, J. Adam: A Method for Stochastic Optimization. In Proceedings of the ACM SIGMOD International Conference on Management of Data, San Diego, CA, USA, 9–12 June 2003.

36. Sancarlos, A.; Cameron, M.; Le Peuvedic, J.M.; Groulier, J.; Duval, J.L.; Cueto, E.; Chinesta, F. Learning stable reduced-order models for hybrid twins. *Data-Centric Eng.* **2021**, *2*, e10. [CrossRef]

37. Erichson, N.; L, L.M.; Kutz, J.; Brunton, S. Randomized dynamic mode decomposition. *SIAM J. Appl. Dyn. Syst.* **2019**, *18*, 1867–1891.