# Adaptive parameters identification for nonlinear dynamics using deep permutation invariant networks

Mouad Elaarabi, Domenico Borzacchiello, Yves Le Guennec, Philippe Le Bot, Sebastien Comas-Cardona

► To cite this version:

Mouad Elaarabi, Domenico Borzacchiello, Yves Le Guennec, Philippe Le Bot, Sebastien Comas-Cardona. Adaptive parameters identification for nonlinear dynamics using deep permutation invariant networks. Machine Learning, 2025, 114 (1), pp.22. 10.1007/s10994-024-06732-7. hal-04892926

HAL Id: hal-04892926

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

Submitted on 20 May 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
# Adaptive parameters identification for nonlinear dynamics using deep permutation invariant networks

Mouad Elaarabi¹,² · Domenico Borzacchiello¹ · Philippe Le Bot² · Yves L. E. Guennec² · Sebastien Comas-Cardona¹

Received: 5 May 2024 / Revised: 15 August 2024 / Accepted: 19 September 2024 / Published online: 16 January 2025
© The Author(s) 2025

## Abstract

The promising outcomes of dynamical system identification techniques, such as SINDy (Brunton et al. in Proc Natl Acad Sci 113(15):3932–3937, 2016), highlight their advantages in providing qualitative interpretability and extrapolation compared to non-interpretable deep neural networks (Rudin in Nat Mach Intell 1(5):206–215, 2019). These techniques suffer from parameter updating in real-time use cases, especially when the system parameters are likely to change during or between processes. Recently, the OASIS (Bhadiraju et al. in AIChE J 66(11):16980, 2020) framework introduced a data-driven technique to address the limitations of real-time dynamical system parameters updating, yielding interesting results. Nevertheless, we show in this work that superior performance can be achieved using more advanced model architectures. We present an innovative encoding approach, based mainly on the use of Set Encoding methods of sequence data, which give accurate adaptive model identification for complex dynamic systems, with variable input time series length. Two Set Encoding methods are used: the first is Deep Set (Zaheer et al. in Adv Neural Inf Process Syst 30, 2017), and the second is Set Transformer (Lee et al. in: International conference on machine learning, PMLR, pp 3744–3753 2019). Comparing Set Transformer to OASIS framework on Lotka–Volterra for real-time local dynamical system identification and time series forecasting, we find that the Set Transformer architecture is well adapted to learning relationships within data sets. We then compare the two Set Encoding methods based on the Lorenz system for online global dynamical system identification. Finally, we trained a Deep Set model to perform identification and characterization of abnormalities for 1D heat-transfer problem.

**Keywords** Data-driven dynamical system identification · Neural networks · Deep permutation invariant networks · Online parameter updating · Sparse regression

Editor: Michael Kamp.

Extended author information available on the last page of the article
# 1 Introduction

Exploring dynamic systems through differential equations is intriguing, providing valuable insights into variable interactions and system behavior. These insights can be useful for many tasks, such as monitoring and controlling industrial processes. However, this could be a challenging undertaking due to their inherent complexity (such as variability of parameters and representativeness of the model). The direct construction of these equations from theoretical principles is a significant challenge and may be sensitive to changes in the parameters or modeling hypotheses, leading researchers to rely on alternative methods such as analyzing empirical data through simulations or experiments. This process involves generating datasets and conducting meticulous analyses using appropriate methods, offering valuable insights about the intrinsic interactions between variables and enhancing the understanding of physical system dynamics.

While machine learning methods offer a potential solution, they suffer from issues such as a lack of interpretability as the trained model is a black box and generally requires extensive training data for high accuracy, which is not always readily available. These methods may also struggle to generalize for unseen data (Rudin, 2019). On the other hand, model discovery methods (Makke & Chawla, 2024; Lejeune et al., 2023) aims to identify the true governing equation through empirical data, providing insights into the interactions between the data, which could be applied to a wide range of engineering tasks.

## 1.1 Background and motivations

In this work, we focus on dynamic system identification methods (North et al., 2023) that have demonstrated excellent results. SINDy (Brunton et al., 2016) (Please refer to Table 1 for abbreviations), based on sparse regression, aims to extract the dynamical links within data. The main idea is to select a few optimal candidate functions from pre-defined over-complete dictionaries, which could be defined based on prior knowledge about the physical laws governing the process and effectively describe the dynamics. This is achieved by using a sparse regression method (Tibshirani, 1996; Manucci et al., 2022; Zheng et al., 2018). On the other hand, PDE Net, proposed by Long et al. (2019), aims to identify PDE equations by discretizing the PDE using forward Euler in time and finite difference in space. This method introduces a different approach mainly based on a succession of symbolic neural networks that aim to link multiple candidates using the Equation Learning (EQL) architecture (Martius & Lampert, 2016) and select the best ones using learnable convolutional kernels. Other methods use a different approach, SGA-PDE (Chen et al., 2022) proposes the use of symbolic mathematics tree representation with classical and differential operators to represent the PDE and a genetic algorithm to find the best representation tree. In the same context, PySR (Cranmer, 2023) introduces a Python library for symbolic regression based on multi-population evolutionary algorithm. This library is particularly well-suited for identifying ODEs, but has limited capabilities in identifying PDEs. This limitation arises from the fact that PySR's operators are restricted to basic algebraic operators. In contrast, SGA-PDE (Chen et al., 2022) extends this capability by incorporating differential operators. These methods offer high potential as the identified parameters can be transferred and used within traditional CAE simulation software finite elements. Additionally, they can be coupled with DNN such as PINN (Raissi et al., 2019) to identify and simulate the dynamics simultaneously as proposed by Chen et al. (2021).
<table><caption>**Table 1** Abbreviation table</caption><thead><tr><td>Abbreviation</td><td>Definition</td></tr></thead><tbody><tr><td>CAE</td><td>Computer-aided engineering</td></tr><tr><td>DMD</td><td>Dynamic mode decomposition</td></tr><tr><td>DMDc</td><td>Dynamic mode decomposition with control</td></tr><tr><td>DNN</td><td>Deep neural network</td></tr><tr><td>EQL</td><td>Equation learning</td></tr><tr><td>ISAB</td><td>Induced self-attention block</td></tr><tr><td>LASSO</td><td>Least absolute shrinkage and selection operator</td></tr><tr><td>MAB</td><td>Multi-head attention block</td></tr><tr><td>MAPE</td><td>Mean absolute percentage error</td></tr><tr><td>MLP</td><td>MultiLayer Perceptron</td></tr><tr><td>MPC</td><td>Model predictive control</td></tr><tr><td>NSODE</td><td>Neural symbolic ordinary differential equation</td></tr><tr><td>OASIS</td><td>Operable adaptive sparse identification of systems</td></tr><tr><td>OASIS-P</td><td>Online adaptive sparse identification of systems for fault prognosis</td></tr><tr><td>ODE</td><td>Ordinary differential equation</td></tr><tr><td>PINN</td><td>Physics informed neural network</td></tr><tr><td>PINN-SR</td><td>Physics informed neural network with sparse regression</td></tr><tr><td>PDE</td><td>Partial differential equation</td></tr><tr><td>PMA</td><td>Pooling by multihead attention</td></tr><tr><td>PySR</td><td>Python symbolic regression</td></tr><tr><td>SAB</td><td>Self-attention block</td></tr><tr><td>SENDI</td><td>Set encoding for nonlinear dynamics identification</td></tr><tr><td>SGA-PDE</td><td>Symbolic genetic algorithm for discovering open-form partial differential equations</td></tr><tr><td>SINDYc</td><td>Sparse identification of nonlinear dynamics with control</td></tr><tr><td>SINDy</td><td>Sparse identification of nonlinear dynamics</td></tr><tr><td>SR3</td><td>Sparse relaxed regularized regression</td></tr><tr><td>rFF</td><td>Row-wise FeedForward</td></tr><tr><td>sMAPE</td><td>Symmetric mean absolute percentage error</td></tr></tbody></table>

Nevertheless, these methods are not pertinent (inference time, adaptability, and leveraging new data to enhance accuracy...) in cases of online identification, where abrupt changes may occur during the process. In such scenarios, the model needs to rapidly adjust to these changes. SINDy-Model Predictive Control (MPC) (Kaiser et al., 2018) demonstrates that a combination of SINDy with control (Fasel et al., 2021; Brunton et al., 2016), built on classical SINDy and DMDc (Narasingam & Kwon, 2017; Proctor, 2016), allows for better control of the system dynamics after a change. This method compared SINDYc, DMDc, and a DNN, revealing that the combination of SINDYc with DMDc yields interesting results. Moreover Quade et al. (2018) suggests an adaptive model for abrupt changes based on SINDy. This method involves two processes: the first detects the changes, while the second identifies and then modifies, adds, or removes parameters from the dynamic system.

These methods present significant potential but have limitations when the identification time is crucial. A slight delay in identification can lead to substantial changes, such as in process systems prognosis (Bhadriraju et al., 2021). Developing an accurate dynamics
identification model using limited data is not viable and may take a considerably long time, hindering its use for interactive applications. Therefore, using SINDYc or abrupt-SINDy online for real applications is challenging.

Recently, Kamienny et al. (2022) presented a sequence to sequence Transformer (Vaswani et al., 2017) architecture to train a general model to perform symbolic regression identification, with faster inference time. This technique is then generalized to ODE identifications, as presented in ODEFormer (d'Ascoli et al., 2023) and NSODE (Becker et al., 2023). They introduced a transformer-based method to identify ODEs structure and parameters, utilizing Symbolic Regression techniques. While the first method targets univariate ODEs, the second extends to multivariate ODEs. These methods focus on training a general model to conduct forward identification of symbolic representations of ODEs based on time series observations. ODEFormer d'Ascoli et al. (2023) exhibits promising results compared to alternative methods. However, this approach is limited when a partial or full knowledge of the dynamics terms are available. Also these methods necessitate an embedding step and an additional fine-tuning step to refine the identified ODE parameters, posing challenges for real-world applications. Hence, the high complexity of the architecture may pose limitations, especially with long time series. Finally, since the model output uses natural language, it would be challenging to couple these methods with other deep learning methods, such as PINNSR (Chen et al., 2021).

In the same context, OASIS (Bhadriraju et al., 2020) has presented a framework that integrates both the strengths of SINDy for identifying the dynamics parameters and a DNN trained offline with parameters identified by SINDy. The latter is subsequently used online to provide parameter predictions for predictive control. OASIS has also been applied for fault real-time prognosis and demonstrated interesting and accurate results (Bhadriraju et al., 2021). The trained neural network uses a single time frame input which cannot handle information about sequence data. Additionally, classical neural network architectures may fail to encode the interactions between data, and therefore struggle to capture the dynamical system, especially with complex systems.

## 1.2 Paper contributions

We introduce a deep learning architecture for evolving online dynamics identification based on Set Encoding. Our goal is to enhance the encoding of interactions within data, leading to improved prediction accuracy of system parameters. We employ two distinct architectures Deep Set (Zaheer et al., 2017) and Set Transformer (Lee et al., 2019). These methods use two different encoding techniques, that can learn both simple pairwise relationships and higher-order links between set data, which can be applied to diverse datasets, making it suitable for a wide range of applications including supervised (our work) and unsupervised learning. In this work, we are particularly interested in applying this technique to sequential data, notably time series. Our main focus is to combine these methods with SINDy to provide an identification framework of the dynamics parameters of a fixed ODE for online use cases such as time series forecasting.

The problem addressed in this article is defined in Eq. 1:

$$ \left\{ \begin{array}{l} \frac{dX}{dt} = \Theta(X)\Xi \quad \text{where } X = [x_0, \dots, x_n] \\ \Theta(X) = [1, X, X^2, \dots, X^2 \otimes \sin(X), \cos(X), \dots], \\ G(\Theta(X)) = \Xi \end{array} \right. \qquad (1) $$
where $\Theta(X)$ is a pre-defined search library based on knowledge about the dynamical system, $X$ is a vector containing the variables of the dynamical system, $\Xi$ represents the parameters of the dynamical system to be calculated. The main goal is to train the Set Encoding methods represented as $G$ to predict the $\Xi$ matrix, thereby identifying the dynamics system. The main contributions of the study are to:

* Propose an appropriate approach to encode and decode to effectively learn the interactions between data, which provide accurate predictions of dynamics parameters. This architecture can handle sequential data with variable sizes, offering enhanced predictions compared to single-time point predictions.

* Illustrate that integrating physics loss in offline learning contributes to a better convergence and generalization of the model.

* Demonstrate that these methods can provide robust predictions with data containing noise.

* Provide accurate predictions (after offline training) with an inference time of milliseconds, making it suitable for real-time applications.

This article is structured as follows: the Sect. 2 will be dedicated to a brief presentation of related works such as SINDYc which will be used in the following sections, and the OASIS framework along with its training method. Then, in Sect. 3, we introduce the two set encoding architectures and discuss their use for online dynamics identification. Section 4 will focus on comparing OASIS and set encoding for local dynamics identification using the Lotka–Volterra system, followed by a comparison of the two set encoding architectures for global dynamics identification using the Lorenz system, then using Deep Set for 1D heat-transfer identification and characterization. Finally, a conclusion of this work and areas for improvement are given.

## 2 Related works

In this section, we briefly introduce SINDYc (Sect. 2.1) and OASIS (Sect. 2.2). SINDYc is utilized to generate training data for comparing OASIS and our methods in the first application with the Lotka–Volterra system (Sect. 4.2), as well as for comparing the two Set Encoding methods presented in Sect. 3 for the second application with the Lorenz system (Sect. 4.3).

### 2.1 SINDYc

SINDYc (Brunton et al., 2016) method allows the identification of parameters of a dynamical system with control variables. Drawing upon prior knowledge of the system’s dynamic behavior, a set of candidate functions, referred to as the search space or library denoted by $\Theta$, is defined. These functions can collectively represent the dynamical system. The underlying idea is primarily based on the notion that an identified system should encompass a limited and concise set of functions. This approach can be applied to various dynamic systems (ODE and PDE), as well as to dynamic systems with control variables (Eq. 2), which are external variables.
$$ \dot{X} = f(X, U), \quad X(t_0) = X_0 \qquad (2) $$

Where $X = [x_1 \ x_2 \ \dots \ x_l] \in \mathbb{R}^{n \times l}$ is a matrix with $l$ state variables with $n$ time snapshots, $\dot{X} = [\dot{x}_1 \ \dot{x}_2 \ \dots \ \dot{x}_l] \in \mathbb{R}^{n \times l}$ is the time derivative of $X$, usually approximated using total variation regularized derivatives (Chartrand, 2011), and $U = [u_1 \ u_2 \ \dots \ u_p] \in \mathbb{R}^{n \times p}$ are the $p$ control variables with $n$ time snapshots, and $X(t_0)$ are the initial conditions. The equation identified by SINDYc is defined as:

$$ \dot{X} = \Theta(X, U)\Xi, \qquad (3) $$

whereas $\Xi$ is a coefficient matrix representing the selected functions from the search library $\Theta(X, U)$, which is defined by multiple functions (Eq. 4) usually using partial knowledge of the physics, for example, the library function could be defined as:

$$ \Theta(X, U) = [1 \ X^T \ U^T \ (X \otimes U)^T \ \dots \ \sin(X^T) \ \dots \ U^T \otimes \exp(X^T) \ \dots] \qquad (4) $$

To identify the dynamical system while avoiding overfitting, the $\Xi$ matrix must be sparse, implying that the less significant coefficients need to be equal to zero. This is achieved by solving the LASSO optimization problem (Eq. 5) defined as:

$$ \Xi = \arg \min_{\Xi} \left[ \| \dot{X} - \Theta(X, U) \Xi \|_2 + \lambda \| \Xi \|_1 \right], \qquad (5) $$

where the first term represents the residual error and the second term is a regularizer applied to the coefficients matrix $\Xi$ or the parameters of the ODE/PDE, and $\lambda$ is a weight parameter that imposes regularization. The larger it is, the stronger the regularization becomes, leading to the elimination of more terms from the search library while potentially reducing the accuracy. Other methods could be used such as Sparse Relaxed Regularized Regression (SR3) (Zheng et al., 2018), or thresholded least squares that use the norm 0 to eliminate less significant coefficients (Brunton et al., 2016; Rudy et al., 2017; Manucci et al., 2022).

Despite all the advantages of SINDYc, this method is restricted in cases where system parameters tend to vary during the process, or whenever the system must be identified online based on sensor data. While abrupt-SINDy (Quade et al., 2018) presents a method to address this issue, it has drawbacks, notably the quantity of data required for identification and the computational cost, which can be significant for real-time cases. In this study we will focus on using SINDy and SINDYc as offline identification methods to generate some datasets, which are then used to train Deep Neural Networks to perform the approximation of parameters online.

## 2.2 OASIS

The OASIS framework is primarily based on SINDYc but unlike the latter, it favors the use of offline trained DNN through results from SINDYc, making it more suitable for real-time parameter identification cases.

The dynamical system (Eq. 2) is a set of measured or generated time series. Then, the set is divided into multiple sub domains (Fig. 1a), which are further used for offline identification through SINDYc (Eq. 3) using an identical search library (Fig. 1b).

This dataset is then utilized for training the DNN (Fig. 2) model (Fig. 1c), where the input comprises state and control variables and the output corresponds to the predicted $\Xi$.
<figure><img src="image_2.png" /><figcaption>Fig. 1 OASIS training process. $X_{jk} \in \mathbb{R}^{N' \times l}$ and $U_{jk} \in \mathbb{R}^{N' \times p}$ represent respectively state variables and control variables for the j-th set and the k-th sub domain, where $j \in [1, S]$ and $k \in [1, m]$. Here, S is the number of sets and m is the total number of sub domains. l, p, and N' represent respectively the number of state variables, control variables, and time windows</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2 Deep neural network with two hidden layers with two inputs denoted as X and U, and five outputs denoted as $\Xi$. The function F, which is trainable, depends on a set of variables denoted as $\Theta$ that governs the connections between these hidden layers</figcaption></figure>

The training process involves leveraging the information contained in the datasets, enabling the model to learn and capture the underlying patterns and relationships within the data, allowing accurate and online identification (Fig. 1d). This is achieved by minimizing the error between the predictions of the DNN and the true coefficients (Eq. 6).

$$ \arg \min_{\Theta} [F(X, U; \Theta) - \Xi + \beta \|\Theta\|_2] \quad \text{where} \quad F(X, U; \Theta) = \hat{\Xi} \qquad (6) $$

Here, F is the output prediction of the DNN (Fig. 2), and $\Theta$ denotes its trainable parameters. Additionally, it is important to include regularization terms $\beta \|\Theta\|_2$ to mitigate overfitting problems, where $\Theta$ represents the trainable variables of the DNN.

The trained model is then employed in real-time for online identification of the dynamical system (Eq. 2) coefficients $\Xi$, enabling accurate extrapolation or prediction of optimal control parameters.
Because the DNN inputs are the current states of the system X and controller U, the OASIS framework does not aim to identify the global dynamics; this means the primary goal is not to identify the true and global coefficients of the governing ODE, but rather a local dynamic system limited to a small time range, and that provides limited extrapolation or prediction of optimal control parameters. This approach is motivated by the fact that accurate identification of the global dynamics is not straightforward and not always feasible (Bhadriraju et al., 2021).

In seeking to identify a valid local dynamics within a limited time interval, frequent model updates are necessary, which can lead to computational issues. To address this concern, the strategy of OASIS-P (Bhadriraju et al., 2021) aims at choosing an optimal frequency for model updating, minimizing the identification error and consequently the extrapolation error.

The fact that the networks use only a single state to predict the parameters defining the local system behaviors poses a potential restriction whenever the dynamics cannot be captured from a single time frame but require a longer observation window. This approach is expected to yield accurate predictions (example in Fig. 3 top left). However, the limitation is attributed to the non-uniqueness of the solution with a single state, as demonstrated by the example in Fig. 3 top right where the intersection of the two time series may introduce confusion during the model training process. The Fig. 3 bottom presents two examples of the Lokta–Volterra dynamics with similar initial conditions and control variables but with small changes in the dynamics parameters. The black arrows show cases where the model inputs are similar ($x(t_i)$ and $y(t_i)$) for the two examples (example 1 in blue and example 2 in red), but the dynamics differ. This aligns with the model’s inadequate capacity to capture temporal information effectively, resulting in less accurate predictions. The following section will present a more suitable

<figure><img src="image_4.png" /><figcaption>**Fig. 3** In top: two examples showcase the dynamics where the model can easily capture the local behavior with only one point in time (left image) and another one where the model will fail to capture the local behavior due to confusion (image at the right). In the bottom: example of dynamics with two state variables system (Lotka–Volterra) where the model could fail to capture the local dynamics at certain points</figcaption></figure>

<figure><img src="image_5.png" /></figure>
activation functions (Eq. 8). In our cases the input is a set of $X_{jk} \in \mathbb{R}^{N' \times l}$ and $U_{jk} \in \mathbb{U}^{N' \times p}$, where $l, p$, and $N'$ represent respectively the number of state variables, control variables, and the time series length. This transformation can be made permutation-equivariant by employing suitable layers. This implies that the output remains equivariant to permutations of set elements. This objective can be accomplished by defining the neural network layers as in Eq. 7, as demonstrated by Zaheer et al. (2017).

$$ \left\{ \begin{array}{l} \mathbf{f}(\mathbf{x}) \doteq \sigma(\lambda \mathbf{I} \mathbf{x} + \gamma \mathbf{pool}(\mathbf{x})) \\ \sigma: \text{ nonlinear activation function} \\ \lambda, \gamma: \text{ learnable variables} \\ \mathbf{pool}(\mathbf{x}): \text{ max, sum, mean} \dots \end{array} \right. \qquad (7) $$

Following this encoding, a pooling transformation is applied (such as max, mean, sum, etc.) with respect to the time dimension which yields a set of latent variables $Z \in \mathbb{R}^{F \times L}$, as defined in Eq. 8:

$$ \left\{ \begin{array}{l} \mathbf{Encoder}(X_{jk}, U_{jk}) = \mathbf{Pool}(\mathbf{DNN}(\{X_{jk}(t_0), \dots, X_{jk}(t_{N'}), U_{jk}(t_0), \dots, U_{jk}(t_{N'})\})) = Z \\ \mathbf{Decoder}(Z) = \mathbf{DNN}(Z) = \Xi_{jk}, \end{array} \right. \qquad (8) $$

where $F$ in $Z \in \mathbb{R}^{F \times L}$ represents the number of encoded features. This operation enables the model to manage inputs of variable shapes, thereby achieving an invariant-type (Eq. 9) encoding.

$$ f(\{x_1, \dots, x_n\}) = f(\{x_{\pi(x_1)}, \dots, x_{\pi(x_n)}\}) \quad \text{for any permutation } \pi \qquad (9) $$

The final step involves decoding the latent vector to obtain the desired output. In this application, we will only utilize the invariant property.

## 3.2 Set transformer

The architecture of set transformer (Fig. 6), is based on transformer model (Vaswani et al., 2017). The motivation behind the use of this architecture lies in its ability to encode links and pairwise interactions among the elements of sets through the attention mechanism. The architecture is primarily built on layers called self attention layers, constructed to be permutation-equivariant.

The attention mechanism can be defined by the Eq. 10:

$$ \mathrm{Att}(Q, K, V; \sigma) = \sigma(QK^T)V. \qquad (10) $$

<figure><img src="image_6.png" /><figcaption>**Fig. 6** Set Transformer Architecture: the encoder is a succession of self-attention layers followed by a pooling layer, which is then decoded with self-attention layers and a DNN</figcaption></figure>
Where $Q$ is a matrix of $n$ sets with the same dimension $d_q$, and $K$ and $V$, referred to as Keys and Values, are two matrices of sets with the same number of sets. The dimension of $K$ is the same as that of $Q$, while $V$ has a different dimension $d_v$. The initial part of the equation, $\sigma(QK^T)$, measures the similarity between the two sets of matrices, where $\sigma$ usually denotes the softmax function applied to the product to convert it into a probability distribution. The result of this product represents a weight matrix, which is used to assign different importance to the elements of the matrix $V$.

Utilizing this mechanism empowers the model to capture relevant patterns and dependencies in the data. To achieve this objective, the mechanism is extended to Multi-head Attention (Eq. 11) defined by three transformations.

$$ \left\{ \begin{array}{l} \textbf{MultiHead}(Q, K, V; \sigma) = \text{concatenation}(O_1, \dots, O_j, \dots, O_L)W^o \\ O_j = \textbf{Att}(QW_j^q, KW_j^k, VW_j^v; \sigma_j) \end{array} \right. \qquad (11) $$

The initial step involves projecting the $Q$, $K$, and $V$ matrices to a higher dimension $d_q^H = d_k^H$ and $d_v^H$ using learnable parameters $W_j^q, W_j^k \in \mathbb{R}^{d_k \times d_k^H}, W_j^v \in \mathbb{R}^{d_v \times d_v^H}$. Subsequently, the attention mechanism is applied to the $L$ resulting projections, referred to as Head Attentions $O_j$. These Head Attentions are then concatenated and multiplied by a trainable variable matrix $W^o \in \mathbb{R}^{L \times d_v^H \times d}$. In this application, it was considered that $d_q = d_v = d$, and $d_q^H = d_v^H = d/L$. The Set Transformer proposes set operation blocks, utilizing a Multi-Head Attention mechanism to encode the data and subsequently aggregate and decode the encoded features. The data is encoded through Multihead Attention Blocks (MAB) and Set-Attention Blocks (SAB). Given two input sets matrices $X, Y \in \mathbb{R}^{n \times d}$, MAB is defined as:

$$ \left\{ \begin{array}{l} \textbf{MAB}(X, Y) = \text{Normalization}(H + \text{rFF}(H)), \\ H = \text{Normalization}(X + \text{MultiHead}(X, Y, Y)). \end{array} \right. \qquad (12) $$

The Normalization is the layer normalization (Ba et al., 2016), and rFF represents a row-wise feed forward neural network without activation functions. The SAB block is defined as:

$$ \textbf{SAB}(X) = \textbf{MAB}(X, X) \qquad (13) $$

Where the Keys and Queries are identical, employing the self-attention mechanism to capture interactions between input features. The encoding is then defined by a succession of SAB blocks (Eq. 14).

$$ \textbf{ENCODER}(X) = \textbf{SAB}(\textbf{SAB}(X)) \qquad (14) $$

The output of the latter is determined by the function $F: X \to Z \in \mathbb{R}^{n \times d}$, where $X \in \mathbb{R}^{n \times d_v}$. This output is subsequently aggregated through the PMA (Pooling Multihead Attention) block defined as:

$$ \left\{ \begin{array}{l} \textbf{DECODER}(Z) = \textbf{rFF}(\textbf{SAB}(\textbf{PMA}(Z))), \\ \textbf{PMA}(Z) = \textbf{MAB}(S, \textbf{rFF}(Z)). \end{array} \right. \qquad (15) $$

Where $S \in \mathbb{R}^{k \times d}$ is a trainable matrix, with $k$ being the desired output features (1 in our case). This allows for an aggregation based on the importance of set elements rather than equal importance (as in the cases of max, mean, etc., aggregation). Following this
transformation, a SAB is applied to better capture the connections between the k output features (Eq. 15).

Considering the quadratic time complexity of the Set Attention Block (SAB) as $O(n^2)$, where n represents the number of sets, Lee et al. (2019) introduced an improved version named Induced Set Attention Block (ISAB). The ISAB offers reduced computational complexity while preserving the same performance.

The architectures of both models are more complex compared to OASIS, which may raise concerns, particularly regarding overfitting problems. We address these concerns in some examples in the next section. The selection of hyperparameters can be managed using classical testing and validation approaches, which we have employed in this work. In the following section, we demonstrate that this complexity offers significant advantages through three different applications. Additionally, we provide solutions and perspectives to mitigate potential overfitting issues.

# 4 Results

## 4.1 Tests cases

To demonstrate the advantages of Set Encoding methods, we present three test cases with three different systems. In the Sect. 4.2 we present a comparison of OASIS framework and Set Transformer for online nonlinear dynamics identification. The Lotka–Volterra (Wangersky, 1978) (Eq. 16) system is used to compare these two methods.

$$ \begin{cases} \dot{x} = \alpha x - \beta xy \\ \dot{y} = \delta xy - \gamma y + c \end{cases} \qquad (16) $$

First, we generate a dataset with various initial conditions and control variables. We then divide the time domain into multiple sub domains, and SINDYc is subsequently employed to identify the local dynamics model for each sub domains. Both methods are then trained to capture the local dynamics using time series and their corresponding identified local SINDYc mode.

While OASIS methods are limited to online local system identification, Set Encoding methods could also be applied in global non-linear dynamics identification. In the Sect. 4.3 we compare the two Set Encoding. These methods can handle variable input sizes, which means capturing the true dynamics through a sequence of data is possible. Furthermore, it can capture the interactions between input data, all of which help enhance the identification process. To study the usefulness of these methods, Set Transformer and Deep Set are compared for the Lorenz system (Eq. 17).

$$ \begin{cases} \dot{x} = \sigma(y(t) - x(t)) \\ \dot{y} = \rho x(t) - y(t) - x(t)z(t) \\ \dot{z} = x(t)y(t) - \beta z(t) \end{cases} \qquad (17) $$

Monitoring and control of industrial processes are generally based on fixing the theoretical parameters of the process and material properties. However, these parameters can change during the process. Therefore, it is important to use sensor data to identify these parameters in real-time. If any deviations are detected, control systems can be used to correct the
process by adjusting the parameters, thereby improving its performance. To address these challenges, Set Encoders can be applied for industrial process monitoring. For example, in the thermo-stamping process, a thin composite plate is heated to make it less rigid before being formed using a specific mold. Using real-time sensor data, the Set Encoder can approximate the material properties during the heating step, which helps, for example, to ensure the temperature is homogeneous throughout the depth of the plate to avoid defects during forming. In this application, we begin with a preliminary study of 1D heating using synthetic data before applying it to real data of the heating step in the thermo-stamping process in future work.

In Sect. 4.4, we will study the 1D heat-transfer problem (PDE example) using a Deep Set model. We will train the model to approximate diffusivity and to characterize an abnormality defined as a local decrease in diffusivity, using noisy data. The equation is given by (Eq. 18):

$$ \left\{ \begin{aligned} \dot{T} &= \frac{\partial(\alpha(z)\partial T)}{\partial z^2} \\ -\frac{\partial T}{\partial z}\Big|_{z=0} &= \frac{q_0}{k}(T) \\ \frac{\partial T}{\partial z}\Big|_{z=L} &= 0 \\ T(t=0, z) &= T_0 \end{aligned} \right. \qquad (18) $$

## 4.2 Application 1: local model identifications (OASIS vs. set transformer)

### 4.2.1 Data generation

The Lotka–Volterra system (Eq. 16) is used to compare these two methods. Where, x and y respectively represent the populations of prey and predators. c is a control variable that limits the rate of increase or decrease of the predators. The parameters α = 0.5, β = 0.025, δ = 0.5, and γ = 0.005 represent respectively the rates of increase or decrease of prey and predators, as well as the effect of predators on prey (and vice versa). An example of a solution of Eq. 16 is given in Fig. 7.

<figure><img src="image_7.png" /><figcaption>**Fig. 7** Example of Lotka–Volterra dynamics solution for α = 0.5, β = 0.025, δ = 0.5, and γ = 0.005 with a constant control variable c = 0, x(t₀) = 27.5 and y(t₀) = 10.0</figcaption></figure>
<table><caption>**Table 2** Range of variables</caption><thead><tr><td>Variable</td><td>Range</td></tr></thead><tbody><tr><td>c: Control variable</td><td>[-1; 5]</td></tr><tr><td>x<sub>0</sub>: Initial condition x</td><td>[5; 50]</td></tr><tr><td>y<sub>0</sub>: Initial condition y</td><td>[5; 15]</td></tr><tr><td>t: Time interval</td><td>[0; 30]</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 8** Left: initial value distribution $(x(t_0), y(t_0))$ for training, validation, and test and Right: the control value distribution $c$ for training, validation, and test</figcaption></figure>

The dataset is generated through various combinations of control variable ($c$) and initial conditions $x_0$ and $y_0$, using the Sobol (Sobol, 1967; Virtanen et al., 2020) sampling method as given in Table 2. The data generation is performed using the LSODA (Virtanen et al., 2020) method from the SciPy Python package. In this comparison (Fig. 8), 278 datasets are generated, of which 166 will be used for model training, and 112 will be employed for model validation to mitigate overfitting issues, and a test set to compare the two methods.

Each dataset is divided into several sub domains over time. SINDYc is then used to identify the local model of the different sub domains, which may be different from the global dynamics but can be used within a short time interval. In this case, we assume that there is no initial information about the search library. The identification is performed using the PySINDy Python package (Silva et al., 2020). We use a base of polynomials of degree 3 defined as:

$$ 
\begin{aligned}
\frac{d\tilde{X}}{dt} &= \Theta(\tilde{X}, C)\Xi \\
\Theta(\tilde{X}, C) &= [1, X, Y, \dots, Y^3, Y^2C, YC^2, C^3] \in \mathbb{R}^{N \times 20}, \\
\tilde{X} &= [X, Y]
\end{aligned}
\qquad (19) $$

where $\Theta(\tilde{X}, C) \in \mathbb{R}^{N \times 20}$ is the search library, where $C \in \mathbb{R}^{N \times 1}$ represents the control variables, and $\Xi \in \mathbb{R}^{2 \times 20}$ is a sparse tensor representing the parameters to be fitted.
<table><caption>**Table 3** Training process summary</caption><thead><tr><th></th><th>Set transformer</th><th>OASIS</th></tr></thead><tbody><tr><td>Optimizer</td><td>Adam</td><td></td></tr><tr><td>Loss function</td><td>L<sub>Ξ</sub> + λ<sub>0</sub>L<sub>ode</sub> + λ<sub>1</sub>L<sub>reg</sub> = ||Ξ - ̂Ξ||<sub>2</sub> + λ<sub>0</sub>||d̂X/dt - Θ(X̃, C)̂Ξ||<sub>1</sub> + λ<sub>1</sub>||W||<sub>1</sub></td><td></td></tr><tr><td>Learning rate</td><td>[10<sup>-3</sup>, 4 × 10<sup>-4</sup>, 10<sup>-4</sup>, 10<sup>-5</sup>, 10<sup>-6</sup>, 10<sup>-7</sup>, 10<sup>-8</sup>]</td><td></td></tr><tr><td>Epochs</td><td>[2000, 2000, 4000, 2000, 2000, 2000, 2000]</td><td></td></tr><tr><td>Input shape</td><td>[Batchsize, time, features]: [Batchsize, 10, 4]</td><td>[Batchsize, 1, 4]</td></tr><tr><td>Output shape</td><td>[Batchsize, coefficients]: [Batchsize, 20]</td><td></td></tr></tbody></table>

### 4.2.2 Training process and criteria for comparison

Given the utility and significance of these methods in real-time control and defect detection, it is crucial to evaluate their extrapolation capabilities. For a fair comparison, both methods underwent training on the same dataset and were then assessed and validated using a separate portion designated for validation. The Table 3 summarizes the training process for both methods. The error function is defined as a sum of three parts, the first one represents the difference between the true $Ξ$ and predicted $\hat{Ξ}$ parameters, the second part represents local dynamic loss by calculating the difference between the known derivative of X and Y and $\Theta(\tilde{X}, C)\hat{Ξ}$ where $\Theta(\tilde{X}, C)$ is known, the final term represents the regularization, contributing to the effective generalization of models. Both methods are trained using identical optimizer and a fixed number of epochs, incorporating variable learning rates during the training process. A noteworthy distinction pertains to the input of the two models. The first maps a single time step to $\hat{Ξ}$, whereas the Set Transformer encodes a temporal series of variable length. We also train the models to predict the dynamics of the subsequent sub domain based on the current sub domain's time series $\tilde{X}_{t_i}$ and control variable $C_{t_i}$ as defined in Eq. 20:

$$F(\tilde{X}_{t_i}, C_{t_i}) = \hat{\Xi}_{t_i+\Delta t} \quad \text{where} \quad \tilde{X}_{t_i} \in \mathbb{R}^{N' \times 2}, \quad \tilde{C}_{t_i} \in \mathbb{R}^{N' \times 1}, \quad \hat{\Xi}_{t_i+\Delta t} \in \mathbb{R}^{20 \times 2} \quad (20)$$

where $N'$ is the time snapshots, $N' = 1$ for OASIS framework and $N' = \Delta t = 10$. Which will help in achieving better results for extrapolation.

To ensure fairness and balance in the comparison, we used hyper-parameters tuning for training both methods, varying hyper-parameters like hidden layers, neurons per layer, regularization coefficient $\lambda_1$, and activation functions. We then select the best architecture by assessing performance through validation error.

### 4.2.3 Comparison and conclusion

The test set is not utilized during the training process; instead, the validation portion is employed solely for selecting the best model during training. We also employ two models for each method to predict the coefficients of $\frac{dx}{dt}$ and $\frac{dy}{dt}$, which we believe will facilitate the training process.

The models are then compared based on their ability to capture local dynamics and extrapolate (Fig. 10) for a limited time range. This process is summarized in Fig. 9. The
<figure><img src="image_9.png" /><figcaption>Fig. 9 Comparison method of OASIS and Set Transformer: Extrapolation using local dynamic identification and ODE-Solver</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 10 Example of extrapolation using Set Transformer in blue and OASIS in red compared to the ground truth in black (HF) simulation results (Color figure online)</figcaption></figure>

input of Set Transformer is a sequence, while the input of OASIS is a single point in time, first the models are used to identify the local dynamics $\hat{\Xi}_{t+\Delta t}^j$ which corresponds to the local dynamics of the following sub domain for the j test dataset. It is then used with initial values $\tilde{X}_{t+\Delta t}$ and control $C_{t+\Delta t}$ value to extrapolate using 'LSODA' as ODE-Solver, the output of the solver is a times series $\tilde{X}_{\text{prediction}}^j = [x_j(t)_p, y_j(t)_p]$ where $t \in [t_i + \Delta t : t_i + 2\Delta t]$, and $\Delta t$ corresponds to the time sub domain windows, which is then compared to the ground truth (simulation data) using MAPE (Eqs. 21 and 22).

$$ \mathbf{MAPE}_{t_i+\Delta t}^j = \left\| 100 \times \frac{x_j(t)_p - x_j(t)_{gt}}{x_j(t)_{gt}} \right\|_1 \quad \text{where} \quad t \in [t_i + \Delta t : t_i + 2\Delta t] \text{ and } j \in [1, N] \qquad (21) $$
$$ \mathbf{MAPE}_{t_i+\Delta t}^{y_j} = \left\| 100 \times \frac{y_j(t)_p - y_j(t)_{gt}}{y_j(t)_{gt}} \right\|_1 \quad \text{where} \quad t \in [t_i + \Delta t : t_i + 2\Delta t] \text{ and } j \in [1, N] \qquad (22) $$

This metric is usually used to measure the prediction accuracy. It quantifies the average percentage deviation between the predicted solution and the ground truth, providing a measure of how well the predictions align with the actual outcomes. Based on the analysis by Hewamalage et al. (2023) of time series forecasting evaluation, MAPE is a suitable metric in the absence of outliers and intermittence. We also utilize sMAPE, which is more suitable in the presence of outliers. This metric quantifies the percentage difference between predicted and actual values, with adjustments made to ensure equal treatment of both small and large values in the evaluation process (see Appendix A.1).

The two models are trained using fixed sub domain windows (sequences of 10 time points referred to as $\Delta t$), as explained in Sects. 4.2.2 and 4.2.1, and then compared using the test dataset, with 3 different sub domain windows 10, 15 and 20, to evaluate the ability of the models to generalize and to capture the local dynamics.

Applying the same logic as explained earlier, we compute the Total MAPE for each dataset (Eqs. 23) and 24).

$$ \mathbf{MAPE}_{\text{Total}}^{x_j} = \left\| 100 \times \frac{x_j(t)_p - x_j(t)_{gt}}{x_j(t)_{gt}} \right\|_1 \quad \text{where } t \in [\Delta t : T_f] \text{ and } j \in [1, N] \qquad (23) $$

$$ \mathbf{MAPE}_{\text{Total}}^{y_j} = \left\| 100 \times \frac{y_j(t)_p - y_j(t)_{gt}}{y_j(t)_{gt}} \right\|_1 \quad \text{where } t \in [\Delta t : T_f] \text{ and } j \in [1, N] \qquad (24) $$

This metric is then used to compare the two methods. Based on Hewamalage et al. (2023) MAPE is not a suitable metric in the presence of outliers. To ensure a fair comparison we first, eliminate outlier results using the normal Z-score (or Six Sigma rules). This means that if the MAPE of a dataset performs significantly better or worse compared to others, it will not be included in the comparison of the two methods. Then we calculate the average of MAPE$^{x_j}$ and MAPE$^{y_j}$ over all test datasets for $x(t)$ and $y(t)$ respectively. In addition, we calculate the median and the 90th percentile. This will provide a more comprehensive understanding of the accuracy of the model. We conduct the same comparison without

<table><caption>**Table 4** Comparison between the two methods set transformer and OASIS in terms of MAPE error for extrapolation through local dynamics identification with a time window of 10 steps</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Mean</th><th colspan="2">90%</th><th colspan="2">Divergence (%)</th></tr><tr><th>SETTR</th><th>OASIS</th><th>SETTR</th><th>OASIS</th><th>SETTR</th><th>OASIS</th></tr></thead><tbody><tr><td rowspan="2">1 × Δt<sub>10</sub></td><td>x</td><td><strong>0.43</strong></td><td>1.01</td><td><strong>0.76</strong></td><td>3.75</td><td><strong>5.36</strong></td><td>44.64</td></tr><tr><td>y</td><td><strong>0.67</strong></td><td>1.90</td><td><strong>1.07</strong></td><td>3.27</td><td><strong>5.36</strong></td><td>44.64</td></tr><tr><td rowspan="2">2 × Δt<sub>10</sub></td><td>x</td><td><strong>1.35</strong></td><td>2.35</td><td><strong>3.63</strong></td><td>9.73</td><td><strong>5.36</strong></td><td>42.86</td></tr><tr><td>y</td><td><strong>3.65</strong></td><td>6.31</td><td><strong>17.01</strong></td><td>21.60</td><td><strong>5.36</strong></td><td>42.86</td></tr><tr><td rowspan="2">3 × Δt<sub>10</sub></td><td>x</td><td><strong>1.93</strong></td><td>3.21</td><td><strong>4.21</strong></td><td>10.28</td><td><strong>7.14</strong></td><td>41.07</td></tr><tr><td>y</td><td><strong>3.28</strong></td><td>6.24</td><td><strong>11.06</strong></td><td>19.00</td><td><strong>7.14</strong></td><td>41.07</td></tr><tr><td rowspan="2">4 × Δt<sub>10</sub></td><td>x</td><td><strong>4.09</strong></td><td>5.44</td><td><strong>12.20</strong></td><td>15.40</td><td><strong>3.57</strong></td><td>39.29</td></tr><tr><td>y</td><td><strong>5.34</strong></td><td>8.45</td><td><strong>18.76</strong></td><td>24.83</td><td><strong>3.57</strong></td><td>39.29</td></tr></tbody></table>
<table><caption>**Table 5** Comparison between the two methods set transformer and OASIS in terms of MAPE error for extrapolation through local dynamics identification with a time window of 15 steps</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Mean</th><th colspan="2">90%</th><th colspan="2">Divergence (%)</th></tr><tr><th>SETTR</th><th>OASIS</th><th>SETTR</th><th>OASIS</th><th>SETTR</th><th>OASIS</th></tr></thead><tbody><tr><td rowspan="2">1 × Δt<sub>15</sub></td><td>x</td><td>1.00</td><td><strong>0.83</strong></td><td>1.54</td><td><strong>1.22</strong></td><td><strong>7.14</strong></td><td>51.79</td></tr><tr><td>y</td><td><strong>0.88</strong></td><td>2.92</td><td><strong>1.50</strong></td><td>4.63</td><td><strong>7.14</strong></td><td>51.79</td></tr><tr><td rowspan="2">2 × Δt<sub>15</sub></td><td>x</td><td><strong>2.41</strong></td><td>2.81</td><td><strong>3.86</strong></td><td>5.12</td><td><strong>7.14</strong></td><td>51.79</td></tr><tr><td>y</td><td><strong>2.26</strong></td><td>5.49</td><td><strong>2.75</strong></td><td>8.98</td><td><strong>7.14</strong></td><td>51.79</td></tr><tr><td rowspan="2">3 × Δt<sub>15</sub></td><td>x</td><td><strong>3.49</strong></td><td>5.11</td><td><strong>5.84</strong></td><td>11.79</td><td><strong>5.36</strong></td><td>44.64</td></tr><tr><td>y</td><td><strong>2.88</strong></td><td>7.83</td><td><strong>5.56</strong></td><td>19.44</td><td><strong>5.36</strong></td><td>44.64</td></tr><tr><td rowspan="2">4 × Δt<sub>15</sub></td><td>x</td><td><strong>6.24</strong></td><td>8.09</td><td><strong>14.70</strong></td><td>22.14</td><td><strong>3.57</strong></td><td>44.64</td></tr><tr><td>y</td><td><strong>5.20</strong></td><td>11.80</td><td><strong>16.38</strong></td><td>38.38</td><td><strong>3.57</strong></td><td>44.64</td></tr></tbody></table>

<table><caption>**Table 6** Comparison between the two methods set transformer and OASIS in terms of MAPE error for extrapolation through local dynamics identification with a time window of 20 steps</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Mean</th><th colspan="2">90%</th><th colspan="2">Divergence (%)</th></tr><tr><th>SETTR</th><th>OASIS</th><th>SETTR</th><th>OASIS</th><th>SETTR</th><th>OASIS</th></tr></thead><tbody><tr><td rowspan="2">1 × Δt<sub>20</sub></td><td>x</td><td>2.42</td><td><strong>1.09</strong></td><td>3.99</td><td><strong>1.84</strong></td><td><strong>7.14</strong></td><td>55.36</td></tr><tr><td>y</td><td><strong>1.03</strong></td><td>3.89</td><td><strong>1.78</strong></td><td>6.15</td><td><strong>7.14</strong></td><td>55.36</td></tr><tr><td rowspan="2">2 × Δt<sub>20</sub></td><td>x</td><td>5.10</td><td><strong>4.41</strong></td><td><strong>8.79</strong></td><td>9.15</td><td><strong>7.14</strong></td><td>50.00</td></tr><tr><td>y</td><td><strong>2.08</strong></td><td>7.14</td><td><strong>3.79</strong></td><td>11.86</td><td><strong>7.14</strong></td><td>50.00</td></tr><tr><td rowspan="2">3 × Δt<sub>20</sub></td><td>x</td><td><strong>9.90</strong></td><td>10.17</td><td><strong>15.92</strong></td><td>23.57</td><td><strong>7.14</strong></td><td>44.64</td></tr><tr><td>y</td><td><strong>6.51</strong></td><td>14.75</td><td><strong>16.66</strong></td><td>39.65</td><td><strong>7.14</strong></td><td>44.64</td></tr><tr><td rowspan="2">4 × Δt<sub>20</sub></td><td>x</td><td><strong>9.28</strong></td><td>10.25</td><td><strong>14.16</strong></td><td>24.59</td><td><strong>3.57</strong></td><td>53.57</td></tr><tr><td>y</td><td><strong>4.46</strong></td><td>8.95</td><td><strong>9.60</strong></td><td>23.78</td><td><strong>3.57</strong></td><td>53.57</td></tr></tbody></table>

eliminating outlier results, using the sMAPE metric, which is more suitable in the presence of outliers (see Appendix A.1).

During these tests, we observed that OASIS and Set Transformer may fail to predict the local dynamics, which can result in divergence. In some cases, this divergence is acceptable (Samples 11 and 21 in Fig. 10), but we also found that it is possible for the model to completely fail to capture the dynamics, leading to infinite divergences. In this case, the dataset is not included in the comparison using the MAPE metric. To represent this, we added the divergence metric to highlight the model's ability to avoid complete divergences.

The results are presented in Tables 4, 5 and 6 which correspond to the 3 sub domains windows sizes as previously explained ($\Delta t_{10}$, $\Delta t_{15}$ and $\Delta t_{20}$), with four different extrapolation windows. This means that based on one sub domain window, we try to extrapolate for larger windows. We refer to these as extrapolation windows.

The primary attention lies in the divergence difference between the Set Transformer and OASIS frameworks, which increases with larger sub domain windows with 40–55% of infinite divergences. In contrast, Set Transformer remains stable with significantly lower divergence percentages. This phenomenon is understandable because the set encoding architecture attempts to encode all time series data to generate predictions. Additionally, as the quantity of data increases, so does the accuracy. We can also see that both models give very acceptable results for small extrapolation windows, this could be highly related to the training process of the two models, where we focus on predicting only the next sub
<figure><img src="image_11.png" /><figcaption>**Fig. 11** Loss comparison when using $\mathcal{L}_{ode}$ with $\lambda = 1$ versus when $\mathcal{L}_{ode}$ is not used, with $\lambda = 0$. For the Set Transformer method using the same architecture</figcaption></figure>

domain dynamics. While the average of MAPE remains good for large extrapolation windows ($2 \times \Delta t$, $3 \times \Delta t$ and $4 \times \Delta t$). The inaccuracy, or the 90th percentile represents the maximum MAPE error observed in 90% of the test sets, increases significantly for OASIS, while it increases reasonably for Set Transformer maintaining acceptable extrapolation.

On the other hand, we concluded that adding the ODE loss $\mathcal{L}_{ode} = \left\| \frac{d\tilde{X}}{dt} - \Theta(\tilde{X}, C)\hat{\Xi} \right\|_1$ as given in Table 3 significantly improves the results, especially in cases where the dynamics terms are unknown or only partially known, as well as in cases where the goal is to identify local dynamics valid for short time windows. The Fig. 11 represents the ODE loss for training and validation using the Set Transformer method with the same architecture. When $\mathcal{L}_{ode}$ is used with $\lambda = 1$, and when $\lambda = 0$ indicates that the ODE loss is not used, we can easily observe the difference in convergence. This difference can be explained by the fact that $\mathcal{L}_{ode}$ assigns weighted attention to each coefficient. This means that if the contribution of a coefficient is significant, its loss $\mathcal{L}_{ode}$ will also be significant compared to other coefficients, and vice versa.

Based on this comparison, we can see the importance of using Set Encoding for local dynamics identification, which could be utilized for extrapolation, as demonstrated in the example above, or for other use cases such as predictive control. We demonstrated that OASIS may fail to capture the local dynamic, resulting in inaccurate results and divergence, while Set Transformer remains stable and leverages all available data to predict the local dynamics. These results could be improved by using variable sub domain windows during the training process, which would provide the Set Transformer method the ability to generate more accurate results and enhance extrapolation capabilities. This approach is demonstrated in the next example, where we aim to identify the global dynamics using an expanding time window. Furthermore, using this approach with noisy data is going to be challenging, which is related to the high sensitivity of SINDy to noise (Gelß et al., 2019). Additionally, in this case, we only aim to identify a local dynamic, which means that the identification process is more challenging. Moreover, in extrapolation use cases, it is essential to update the model frequently to obtain accurate results. Using the predictions of Set Transformer (or OASIS) and then an ODE solver to extrapolate poses two challenges. The first one is related to the training process, which involves the generation of training data and using SINDy for the identification of local dynamics with noisy data, which is
<figure><img src="image_12.png" /><figcaption>**Fig. 12** Example of the Lorenz system with different noise levels, left 0% ($\xi$ = 0), middle 2% ($\xi$ = 0.02), and right 5% ($\xi$ = 0.05)</figcaption></figure>

<table><caption>**Table 7** Range of variables</caption><thead><tr><td>Variables</td><td>Range</td></tr></thead><tbody><tr><td>σ variations</td><td>[8; 12]</td></tr><tr><td>ρ variations</td><td>[20; 35]</td></tr><tr><td>β variations</td><td>[0; 4]</td></tr><tr><td>x<sub>0</sub> variations</td><td>[-5; 5]</td></tr><tr><td>y<sub>0</sub> variations</td><td>[4; 50]</td></tr><tr><td>z<sub>0</sub> variations</td><td>[5; 15]</td></tr><tr><td>t train and validation</td><td>[0, 10]</td></tr><tr><td>t test</td><td>[0, 20]</td></tr></tbody></table>

challenging. Furthermore, even slight noise in the initial conditions may lead to inaccurate extrapolation or complete divergence.

## 4.3 Application 2: global model identifications (deep set vs. set transformer)

### 4.3.1 Data generation

We generate a dataset of 1024 Lorenz systems with variables parameters [$\sigma$, $\rho$, $\beta$], and initial conditions [$x_0$, $y_0$, $z_0$] as described in Table 7 using LSODA (Virtanen et al., 2020) methods implemented in SciPy Python package. The parameters and initial conditions are sampled using the Sobol (Sobol, 1967; Virtanen et al., 2020) method implemented in the SciPy Python package. As the number of possible combinations of parameters and initials conditions is large and taking into consideration the chaotic behavior of the Lorenz system, we use 70% data to train the model and the remaining 30% for validation and testing.

To study the robustness of the models concerning data quality, we generate 3 datasets with three levels of Gaussian noise. $\xi \in [0, 0.02, 0.05]$ as described in Eq. 25:
$$ \begin{cases} x(t) = x(t) + \xi \epsilon_x & \text{where } \epsilon_x \sim \mathcal{N}(0, \sigma_x^2) \\ y(t) = y(t) + \xi \epsilon_y & \text{where } \epsilon_y \sim \mathcal{N}(0, \sigma_y^2) \\ z(t) = z(t) + \xi \epsilon_z & \text{where } \epsilon_z \sim \mathcal{N}(0, \sigma_z^2) \end{cases} \qquad (25) $$

Figure 12 shows an example with three levels of noise. Using the same parameters and initial conditions, we then compute the approximation of the derivative with respect to time $\frac{d\tilde{X}}{dt} = [\frac{d\tilde{x}}{dt}, \frac{d\tilde{y}}{dt}, \frac{d\tilde{z}}{dt}]$ using total variation method (Chartrand, 2011) for noisy data and central difference (Eq. 26) for clean data.

$$ f'(x) \approx \frac{f(x+h) - f(x-h)}{2h} \qquad (26) $$

Afterwards, we compute $\tilde{X} = [\tilde{x}, \tilde{y}, \tilde{z}]$, which represents the approximation of the integral of $\frac{\partial \tilde{X}}{\partial t}$ using the cumulative trapezoidal integration method (Eq. 27) (with the SciPy Python package).

$$ \tilde{X} = \int_a^b \frac{\partial \tilde{X}}{\partial t} dt \approx \frac{1}{2} \sum_{n=1}^{N-1} (t_{n+1} - t_n) \left[ \frac{\partial \tilde{X}}{\partial t}_{n+1} + \frac{\partial \tilde{X}}{\partial t}_{n} \right] \qquad (27) $$

This approach contributes to noise reduction in the data (Fig. 13), it will also enable the models to learn effective denoising techniques (further explanations are given in Sect. 4.3.2). PySINDy (Python package) is then used to compute the approximation of the parameters for each dataset. In this application we assume that the ODE search library $\Theta$ is known and well defined, our main focus is on calculating the parameters $[\sigma_k, \rho_k, \beta_k]$. The identification process is defined as in Eq. 28.

$$ \begin{bmatrix} \frac{\tilde{d}x}{dt} \\ \frac{\tilde{d}y}{dt} + \tilde{y} + \tilde{x}z \\ \frac{\tilde{d}z}{dt} - \tilde{x}\tilde{y} \end{bmatrix} = [\tilde{x} - \tilde{y}, \tilde{y}, \tilde{z}]\tilde{\Xi} = \Theta(\tilde{x}, \tilde{y}, \tilde{z})\tilde{\Xi} \qquad (28) $$

<figure><img src="image_13.png" /><figcaption>Fig. 13 Example of the Lorenz system: comparison of noisy data and denoised data using the Total Variation method</figcaption></figure>
<table><caption>**Table 8** Training process summary</caption><thead><tr><th></th><th>Deep Set</th><th>Set Transformer</th></tr></thead><tbody><tr><td>Optimizer</td><td>Adam</td><td></td></tr><tr><td>Loss function</td><td>||Ξ̃ - Ξ||<sub>2</sub> + λ<sub>0</sub> ||d̃X/dt - Θ(Ξ̃, C)Ξ||<sub>1</sub></td><td></td></tr><tr><td>Learning rate</td><td>[4 × 10<sup>-4</sup>, 5 × 10<sup>-5</sup>, 10<sup>-5</sup>, 5 × 10<sup>-6</sup>, 10<sup>-6</sup>]</td><td></td></tr><tr><td>Epochs</td><td>[400, 400, 400, 400, 400]</td><td></td></tr><tr><td>Input shape</td><td>[Batchsize, time, Features]</td><td>[Batchsize, variable, 4]</td></tr><tr><td>Output shape</td><td>[Batchsize, 3]</td><td></td></tr></tbody></table>

### 4.3.2 Train process and criteria of comparison

We implement a consistent training regime for both models (Deep Set and Set Transformer), employing varied architectures for both methods. Following this, we conduct a comparison of the optimal architecture configurations for each method. A core principle underlying the utilization of Set Encoding models is the integration of adaptability, enabling predictions that dynamically adjust to an increasing volume of data over time. We employ a sliding window approach with increasing widths. Through this approach, the models are designed to adaptively incorporate new information for updating predictions of Lorenz parameters.

The model inputs are $[t, X_k(t)] = [t, x_k(t), y_k(t), z_k(t)]$, which represent the noisy data, we added the time dimension, which will be encoded during the pooling phase. This will enable the model to better comprehend the dynamic while increasing input size, and hold good accuracy with variables times steps $\Delta t$. The output is defined as $[\sigma_k, \rho_k, \beta_k]$. The training process is summarized in Table 8. We employ the Adam (Kingma & Ba, 2014) optimizer for training both models, with a step decay learning rate. The loss function comprises two components. The first part represents the loss associated with parameters loss where $Ξ̃$ and $Ξ$ denote the fitted parameters using SINDy (Eq. 28) based on denoised data and the predicted parameters. Note here that $Ξ̃$ may differ from $Ξ$, the true parameters vector. The second part represents the loss associated with the ODE loss utilizing denoised variables, in this study, we fix $\lambda_0$ to 1. Training both methods to predict denoised parameters enables the model to additionally learn denoising techniques while processing noisy input data. Both methods are trained and validated on the same dataset. The best model will be selected based on its accuracy in validation; subsequently, the best model instances are chosen for a comparative analysis between the two models, for more details see Appendix A.4.

The comparison of instances of the two models will be based on their ability to deliver accurate and stable predictions, while considering the variability in both the quality and quantity of data. Note here that we trained individual instances to predict each parameter of the Lorenz systems, which helps the models converge faster.

### 4.3.3 Comparison and discussion

The instances comparison are based on three features:
<figure><img src="image_14.png" /><figcaption>Fig. 14 Accuracy results using $R^2_{w.a}$ metric for Deep Set and Set Transformer instances for $\sigma$, $\beta$ and $\rho$. On the left interpolation results for $t \in [0, 10[$, and on the right extrapolation results for $t \in [10, 20[$</figcaption></figure>

1. Accuracy: the accuracy of selected instances is assessed using $R^2$-score metric (Chicco et al., 2021) defined as:

$$R^2 = 1 - \frac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - \bar{y})^2}, \qquad (29)$$

which measures how well the predicted values match the actual values. In our context, it quantifies the precision of predicting Lorenz parameters. The closer the $R^2$ score is to 1, the better the accuracy. It is important to note that our comparison involves evaluating the model predictions against the true Lorenz parameters, and not the fitted SINDy parameters $\hat{\sigma}$. We assess the instances' capability to provide accurate predictions with variable input sizes, we employ the weighted average of $R^2_{w.a}$ defined as:

$$R^2_{w.a} = \sum_i \frac{1/n_i}{F_T} R^2_{n_i} \quad \text{where} \quad F_T = \sum_i 1/n_i, \qquad (30)$$

which assigns more importance to predictions with fewer data points, where $n_i \in [100, 300, 500, 700, 900]$. This study is based on an expanding window approach, with two windows utilized: one for interpolation study where $t \in [0, 10]$, and another for extrapolation where $t \in [10, 20]$, which is never used during training or validation. Results of $R^2_{w.a}$ are presented in Fig. 14. The overall performance of both methods is good for both interpolation and extrapolation ($R^2_{w.a} > 0.8$). The accuracy drops for extrapolation, which is related to new Lorenz patterns unseen during the training process. It is also related to the time dimension used during training, limited to $t \in [0, 10]$. In interpolation, Method Deep Set yields slightly better results, while in extrapolation, Method Set Transformer tends to be more stable (results do not change much). An exception is the prediction of $\sigma$, for which the accuracy drops for both methods. This can be explained by the data generation process (Eq. 28), where we can observe that the prediction of $\sigma$ is based on $x - y$, whereas for the prediction of $\beta$ and $\rho$, it is based directly on $y$ and $z$, respectively, which are the input of the models, making it easier to learn. Furthermore, in extrapolation, Set Transformer performs better in the prediction of $\sigma$, while in interpolation, Deep Set was slightly better. We can also observe that for both interpolation and extrapolation, both methods demonstrate very good robustness
<figure><img src="image_15.png" /><figcaption>**Table 9** Range of characteristics for train, validation and test</figcaption></figure>

<table><thead><tr><th>Variables</th><th>Range</th><th>Distribution</th></tr></thead><tbody><tr><td>α<sub>ref</sub></td><td>[8 × 10<sup>-7</sup>; 1.2 × 10<sup>-6</sup>]</td><td>N(10<sup>-6</sup>, (5 × 10<sup>-8</sup>)<sup>2</sup>)</td></tr><tr><td>ratio</td><td>[40; 80]</td><td>N(60, (8)<sup>2</sup>)</td></tr><tr><td>G</td><td>[0.001; 0.009]</td><td>Sobol</td></tr><tr><td>Time t</td><td>[0; 200]</td><td>—</td></tr></tbody></table>

to different noise levels. We present the $R_{n_i}^2$ with respect to the quantity of data in Appendix A.2 (Figs. 18, 19 and 20).

2. *Inference time*: Both trained models could be used directly for online applications. This is because the inference time is in the order of a few milliseconds on AuthenticAMD AMD EPYC 9474F@3.6 Ghz 48-Core Processor, making these models more suitable for applications like the extrapolation of state variables (Bhadriraju et al., 2021) or online control (Bhadriraju et al., 2020). These methods also can dynamically use new data to enhance the accuracy of the prediction or to reevaluate the prediction in case of changes in the system dynamics.

3. *Training*: Based on the comparison of both methods, we observed that the implementation of Deep Set is easier compared to Set Transformer. Deep Set relies on simple dense layers, resulting in fewer hyperparameters, thus facilitating hyperparameter tuning. More details about the instance selection are provided in Appendix A.4.

## 4.4 Application 3: characterizing abnormalities in 1D heat-transfer using deep set model

### 4.4.1 Data generation

For this application we aim to train a Deep Set model to approximate online the diffusivity value $\alpha$ from temperature data, we also challenge the model to characterize a local decrease of diffusivity. We start by generating a dataset using finite difference methods (implicit Euler method) with different diffusivity values $\alpha$. The 1D heating equation is
defined as in Eq. 18 in this case $\alpha(z)$ is space-dependent. The initial temperature value is $T(t=0, z) = 20 \,^\circ C$, $L = 0.01$ m, and the heat flux is defined as function as follow (Eq. 31):

$$ \left\{ \begin{array}{l} -\frac{\partial T}{\partial z}\bigg|_{z=0} = \frac{q_0}{k} = 30000 \, \text{K/m} \quad \text{if } T(z=0) < T_{\text{max}} \\ -\frac{\partial T}{\partial z}\bigg|_{z=0} = 0 \, \text{K/m} \quad \text{if } T(z=0) > T_{\text{max}} \end{array} \right. \qquad (31) $$

where $T_{max} = 500 \,^\circ C$, which simply means that we impose a regularization to stop the heating if $T > T_{max}$. The $\alpha(z)$ is defined as in Fig. 15, this function could be simplified as three characteristics, $\alpha_{ref}$ which is the normal diffusivity value, $G = \text{center}_{abnormality}$ which presents the center of abnormality and ratio = $\frac{\alpha_{ref}}{\alpha_{abnormality}}$ which is simply the ratio between normal and abnormal diffusivity. The $\alpha(z)$ function defined as (Eq. 32), where A and B are two values that help define $\alpha_{abnormality}$:

$$ \alpha(z) = \alpha_{ref} - Ae^{-B(G-z)^6} \qquad (32) $$

We generate N datasets with different combinations of the three characteristics defined previously, where N = 300, where 70% of this dataset is used for training and 30% for validation, we also generate separately a $N_{test} = 150$ test datasets to evaluate the performance of the model, the distribution of training, validation and test characteristics is defined as follows (Table 9).

### 4.4.2 Training process

Using these datasets, we train a Deep Set model, where the input $X_i$ combines a time series of temperature at two positions $z=0$ and $z=L = 0.01$ m and the corresponding time $t$ (Example in Fig. 16), where $\Delta t = 1$ s. In this study we limited the model input temperature information to only two positions, the input $X_i$ which corresponds to the ith dataset is defined as (Eq. 33):

<figure><img src="image_16.png" /><figcaption>**Fig. 16** Example of the input temperature with noise of the deep set model</figcaption></figure>
$$X_i = [z_i, t_i, T_i] = \begin{bmatrix} z_{i,1} & t_{i,1} & T(t_{i,1}, z_{i,1}) \\ z_{i,1} & t_{i,2} & T(t_{i,2}, z_{i,1}) \\ \vdots & \vdots & \vdots \\ z_{i,k} & t_{i,j} & T(t_{i,j}, z_{i,k}) \\ z_{i,k} & t_{i,j+1} & T(t_{i,j+1}, z_{i,k}) \\ \vdots & \vdots & \vdots \\ z_{i,L} & t_{i,T-1} & T(t_{i,T-1}, z_{i,L}) \\ z_{i,L} & t_{i,T} & T(t_{i,T}, z_{i,L}) \end{bmatrix} \quad \text{where } k \in \{1, L\}, \quad j \in [1, T] \quad (33)$$

Thus the model input and output could be defined as $F(X_i; \Theta) = [\text{center}_i, \text{ratio}_i, \alpha_i]$, where $\Theta$ are the trainable parameters of Deep Set model. The details about the used architecture are given in Appendix A.5. We used a time-expanding window, which gives the model a flexible approximation of the characteristics using available information, which is improving as new data is available. The loss function is defined as (Eq. 34):

$$\left\{
\begin{aligned}
\mathcal{L}_{\text{Total}} &= \mathcal{L}_{\text{center}} + \mathcal{L}_{\text{ratio}} + \mathcal{L}_{\alpha_{\text{ref}}} + \lambda \mathcal{L}_{\Theta} \\
\mathcal{L}_{\text{center}} &= \|G_{\text{true}} - G_{\text{prediction}}\|_2 \\
\mathcal{L}_{\text{ratio}} &= \|\text{ratio}_{\text{true}} - \text{ratio}_{\text{prediction}}\|_2 \\
\mathcal{L}_{\alpha_{\text{ref}}} &= \|\alpha_{\text{true}} - \alpha_{\text{prediction}}\|_2 \\
\mathcal{L}_{\Theta} &= \|\Theta\|_1
\end{aligned}
\right.
, \quad (34)$$

which is composed of four losses: $\mathcal{L}_{\text{center}}$, $\mathcal{L}_{\text{ratio}}$, $\mathcal{L}_{\alpha}$ and $\mathcal{L}_{\Theta}$, which correspond to the losses of each characteristic and the regularization of the model which helps avoid overfitting problems, in this application we found that $\lambda = 0.01$ gives better generalization.

We also study the model robustness to noise, for that we tested the model for four noise levels defined as (Eq. 35) where $\xi \in [0, 0.02, 0.05, 0.1]$ and $i \in [0, N]$:

$$T_i(t) = T_i(t) + \xi \epsilon_T \quad \text{where} \quad \epsilon_T \sim \mathcal{N}(0, \sigma_{T_i}^2) \quad (35)$$

<figure><img src="image_17.png" /><figcaption>**Fig. 17** Accuracy results using $R_{w.a}^2$ metric for Deep Set for abnormalities characterization: G, $\alpha_{\text{ref}}$ and ratio, for 4 different noise levels from $\xi = 0.0$ to $\xi = 0.1$</figcaption></figure>
<table><caption>**Table 10** Comparison between the two methods set-transformer and OASIS in terms of sMAPE error for extrapolation through local dynamics identification with a time window of 10 steps</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Mean</th><th colspan="2">90%</th><th colspan="2">Divergence (%)</th></tr><tr><th>SETTR</th><th>OASIS</th><th>SETTR</th><th>OASIS</th><th>SETTR</th><th>OASIS</th></tr></thead><tbody><tr><td rowspan="2">1 × Δt<sub>10</sub></td><td>x</td><td><strong>0.47</strong></td><td>1.23</td><td><strong>1.11</strong></td><td>3.58</td><td><strong>5.36</strong></td><td>44.64</td></tr><tr><td>y</td><td><strong>1.46</strong></td><td>3.61</td><td><strong>1.23</strong></td><td>6.06</td><td><strong>5.36</strong></td><td>44.64</td></tr><tr><td rowspan="2">2 × Δt<sub>10</sub></td><td>x</td><td><strong>1.96</strong></td><td>3.40</td><td><strong>6.07</strong></td><td>12.20</td><td><strong>5.36</strong></td><td>42.86</td></tr><tr><td>y</td><td><strong>4.45</strong></td><td>6.44</td><td><strong>18.69</strong></td><td>19.46</td><td><strong>5.36</strong></td><td>42.86</td></tr><tr><td rowspan="2">3 × Δt<sub>10</sub></td><td>x</td><td><strong>2.31</strong></td><td>4.23</td><td><strong>5.78</strong></td><td>13.59</td><td><strong>7.14</strong></td><td>41.07</td></tr><tr><td>y</td><td><strong>4.16</strong></td><td>6.23</td><td><strong>14.40</strong></td><td>18.15</td><td><strong>7.14</strong></td><td>41.07</td></tr><tr><td rowspan="2">4 × Δt<sub>10</sub></td><td>x</td><td><strong>4.51</strong></td><td>5.87</td><td><strong>13.44</strong></td><td>16.14</td><td><strong>3.57</strong></td><td>39.29</td></tr><tr><td>y</td><td><strong>6.85</strong></td><td>9.07</td><td>24.34</td><td><strong>21.02</strong></td><td><strong>3.57</strong></td><td>39.29</td></tr></tbody></table>

<table><caption>**Table 11** Comparison between the two methods set-transformer and OASIS in terms of sMAPE error for extrapolation through local dynamics identification with a time window of 15 steps</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Mean</th><th colspan="2">90%</th><th colspan="2">Divergence (%)</th></tr><tr><th>SETTR</th><th>OASIS</th><th>SETTR</th><th>OASIS</th><th>SETTR</th><th>OASIS</th></tr></thead><tbody><tr><td rowspan="2">1 × Δt<sub>15</sub></td><td>x</td><td><strong>1.14</strong></td><td>1.76</td><td><strong>1.87</strong></td><td>3.01</td><td><strong>7.14</strong></td><td>51.79</td></tr><tr><td>y</td><td><strong>1.84</strong></td><td>4.62</td><td><strong>1.81</strong></td><td>8.41</td><td><strong>7.14</strong></td><td>51.79</td></tr><tr><td rowspan="2">2 × Δt<sub>15</sub></td><td>x</td><td><strong>2.69</strong></td><td>3.78</td><td><strong>5.07</strong></td><td>10.55</td><td><strong>7.14</strong></td><td>51.79</td></tr><tr><td>y</td><td><strong>3.37</strong></td><td>5.76</td><td><strong>7.32</strong></td><td>12.71</td><td><strong>7.14</strong></td><td>51.79</td></tr><tr><td rowspan="2">3 × Δt<sub>15</sub></td><td>x</td><td><strong>4.60</strong></td><td>6.13</td><td><strong>6.03</strong></td><td>13.73</td><td><strong>5.36</strong></td><td>44.64</td></tr><tr><td>y</td><td><strong>3.71</strong></td><td>8.48</td><td><strong>9.53</strong></td><td>20.58</td><td><strong>5.36</strong></td><td>44.64</td></tr><tr><td rowspan="2">4 × Δt<sub>15</sub></td><td>x</td><td><strong>7.40</strong></td><td>9.56</td><td><strong>14.88</strong></td><td>24.25</td><td><strong>3.57</strong></td><td>44.64</td></tr><tr><td>y</td><td><strong>7.05</strong></td><td>10.83</td><td><strong>18.18</strong></td><td>31.84</td><td><strong>3.57</strong></td><td>44.64</td></tr></tbody></table>

### 4.4.3 Results and discussion

We trained different models for the four different noise levels $\xi$, we then select the best model based on the loss of the validation datasets, the selected models are then tested using the test datasets. The selected model architecture is given in Appendix A.5. We assess the model accuracy with respect to the quantity data and robustness to noise using the $R^2$ metrics (Eq. 29), we also use the same temperature position used for the training and validation, the results are:

* **Accuracy and robustness:** We use the same accuracy assessment $R^2_{w.a}$ as in application Sect. 4.3:

$$R^2_{w.a} = \sum_i \frac{1/n_i}{F_T} R^2_{n_i} \quad \text{where} \quad F_T = \sum_i 1/n_i, \qquad (36)$$

for this application we test the models for different data quantities where $t \in [0, 200]$ and data quantities $n_i \in [50, 100, 150, 200]$. From these results (Figs. 17 and 21), we can see that the model can easily predict the G of abnormality and remains stable with respect to noise, also the accuracy of the prediction of $\alpha_{ref}$ is good and decreases reasonably with the increase of noise, we believe that this is related to confusion related to the added noise which makes the precision of diffusivity very difficult, this is also
related to the fact that the range of diffusivity is small while the noise levels represent a higher $\Delta T_{noise}$ related to noise (for $\xi = 0.1 \Delta T_{noise} \in [-40, 40]$) which confuse the training of the model. Finally, the prediction of the ratio is acceptable for low noise levels but decreases as the noise level increases. We believe that this is related to the lack of information (only two temperature positions) and the confusion caused by the noise level. We also believe that increasing the number $N$ of datasets will improve the model accuracy. More details on the accuracy of the model with respect to data quantity are given in Appendix A.3.

* Discussion: In this application we focused on the Deep Set architecture, nevertheless Set Transformer, which gives very good results for Lorenz chaotic systems, could also be applied to this application and is expected to give good results. Also, in this application, we did not add a physics loss $\mathcal{L}_{PDE}$ simply because the approximation of $\frac{\partial(\alpha(z)\partial T)}{\partial z^2}$ in the PDE (Eq. 18) is not an easy task. Specifically, in this case where we use only two separated positions, the approximation of the second derivative with respect to space is challenging. This results in less robustness to noise compared to the second application. This problem will be addressed in a subsequent study where we will attempt to improve this method for application to PDE problems.

# 5 Conclusion and perspectives

In this paper, we presented a combination of two architectures of Set Encoding with the methods of SINDy and SINDYc. The Set Encoding methods are initially trained on a dataset generated by SINDy (or SINDYc) and subsequently employed in real-time to swiftly and accurately predict system dynamics parameters, which could be used for applications such as extrapolating state variables or system control. Our primary contribution lies in encoding variable-sized data sequences using an appropriate architecture, thereby enabling the model to better utilize all temporal information to enhance precision. Additionally, we propose a more suitable loss function. In application 1, we illustrate the need to incorporate the $\mathcal{L}_{ode}$ component in the loss function, which enhances the model training process. Our approach can be applied to local identification and global dynamics, as demonstrated in application 2 and 3. In the second application, we suggest a training approach on noisy data, demonstrating that both methods yield robust and stable results with different noise levels. However, this approach relies on datasets generated by SINDy, implying that

<table><caption>**Table 12** Comparison between the two methods set-transformer and OASIS in terms of sMAPE error for extrapolation through local dynamics identification with a time window of 20 steps</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Mean</th><th colspan="2">90%</th><th colspan="2">Divergence (%)</th></tr><tr><th>SETTR</th><th>OASIS</th><th>SETTR</th><th>OASIS</th><th>SETTR</th><th>OASIS</th></tr></thead><tbody><tr><td rowspan="2">1 × Δt<sub>20</sub></td><td>x</td><td><strong>2.63</strong></td><td>2.66</td><td>3.92</td><td><strong>2.24</strong></td><td><strong>7.14</strong></td><td>55.36</td></tr><tr><td>y</td><td><strong>2.79</strong></td><td>5.81</td><td><strong>2.05</strong></td><td>7.52</td><td><strong>7.14</strong></td><td>55.36</td></tr><tr><td rowspan="2">2 × Δt<sub>20</sub></td><td>x</td><td><strong>6.21</strong></td><td>7.35</td><td><strong>11.32</strong></td><td>24.77</td><td><strong>7.14</strong></td><td>50.00</td></tr><tr><td>y</td><td><strong>4.77</strong></td><td>9.57</td><td><strong>9.95</strong></td><td>21.72</td><td><strong>7.14</strong></td><td>50.00</td></tr><tr><td rowspan="2">3 × Δt<sub>20</sub></td><td>x</td><td><strong>9.82</strong></td><td>11.83</td><td><strong>23.89</strong></td><td>25.61</td><td><strong>7.14</strong></td><td>44.64</td></tr><tr><td>y</td><td><strong>7.08</strong></td><td>12.76</td><td><strong>18.22</strong></td><td>33.02</td><td><strong>7.14</strong></td><td>44.64</td></tr><tr><td rowspan="2">4 × Δt<sub>20</sub></td><td>x</td><td><strong>9.88</strong></td><td>12.51</td><td><strong>16.45</strong></td><td>35.09</td><td><strong>3.57</strong></td><td>53.57</td></tr><tr><td>y</td><td><strong>4.92</strong></td><td>9.72</td><td><strong>12.82</strong></td><td>23.82</td><td><strong>3.57</strong></td><td>53.57</td></tr></tbody></table>
<figure><img src="image_18.png" /><figcaption>**Fig. 18** Accuracy comparison of $\sigma$ prediction between Deep Set and Set Transformer using the $R^2$ metric for the test dataset as a function of data quantity with different noise levels: 0 on the left, 0.02 in the center, and 0.05 on the right for interpolation (in the top) where $t \in [0, 10]$ and for extrapolation (in the bottom) where $t \in [10, 20]$</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 19** Accuracy comparison of $\beta$ prediction between Deep Set and Set Transformer using the $R^2$ metric for the test dataset as a function of data quantity with different noise levels: 0 on the left, 0.02 in the center, and 0.05 on the right for interpolation (in the top) where $t \in [0, 10]$ and for extrapolation (in the bottom) where $t \in [10, 20]$</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. 20** Accuracy comparison of $\rho$ prediction between Deep Set and Set Transformer using the $R^2$ metric for the test dataset as a function of data quantity with different noise levels: 0 on the left, 0.02 in the center, and 0.05 on the right for interpolation (in the top) where $t \in [0, 10]$ and for extrapolation (in the bottom) where $t \in [10, 20]$</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 21** Accuracy results using the $R^2$ metric and noise robustness with respect to data quantity for the Deep Set model for abnormalities characterization: G, $\alpha_{ref}$, and ratio, for four different noise levels ranging from $\xi = 0.0$ to $\xi = 0.1$</figcaption></figure>
<table><caption>**Table 13** Selected architecture for Deep Set method for all noise levels and all the parameters of the Lorenz system</caption><thead><tr><th></th><th>Activation</th><th>Layers</th><th>Pooling</th><th>Trainable parameter count</th></tr></thead><tbody><tr><td>Encode</td><td>ReLU</td><td>320 × 5</td><td>abs, mean</td><td>927,043</td></tr><tr><td>Decode</td><td>ReLU</td><td>320 × 5</td><td></td><td></td></tr></tbody></table>

<table><caption>**Table 14** Selected architecture for set transformer method for all noise levels and all the parameters of the Lorenz system</caption><thead><tr><th></th><th>Attention blocks</th><th>rFF</th><th>Trainable parameter count</th></tr></thead><tbody><tr><td>Encode</td><td>ISAB (d = 45, m = 128, h = 40)</td><td>rFF (45, 2, ReLU)</td><td>1,045,733</td></tr><tr><td>Pooling</td><td>PAM (k = 1, d = 40, h = 40)</td><td>rFF (40, 2, ReLU)</td><td></td></tr></tbody></table>

<table><caption>**Table 15** Selected architecture for the deep set method for all noise levels in the third application</caption><thead><tr><th></th><th>Activation</th><th>Layers</th><th>Pooling</th><th>Trainable parameter count</th></tr></thead><tbody><tr><td>Encode</td><td>GeLU</td><td>256 × 5</td><td>Sum</td><td>1,262,083</td></tr><tr><td>Decode</td><td>GeLU</td><td>3 × 256 × 5</td><td></td><td></td></tr></tbody></table>

accuracy is directly dependent on the precision of SINDy (or SINDYc). In the third application, we addressed a PDE example and showed that the model gives interesting results. Nevertheless, the generation of the training datasets is more challenging, especially when dealing with limited and noisy datasets in the context of partial differential equation (PDE) systems, where calculating derivatives in the presence of noise is particularly challenging. To address this issue, an approach akin to PINN-SR (Chen et al., 2021) can be employed, extending our methodology to handle ODEs and PDEs with noisy and limited data quantities. In addition, the current Set Encoding architecture cannot handle cases where terms in differential equations are prone to deletion or addition in real-time applications. This implies that the models cannot activate or deactivate terms. To address this challenge, one can integrate a new decoding block with an appropriate activation function (such as a sigmoid function), which would be responsible for activating and deactivating dynamic system terms. Finally, we believe that a comparison between Deep Set and Set Transformer in online nonlinear dynamics identification is necessary to better understand the difference between these two methods. While Deep Set encodes data using classical neural network architecture, Set Transformer employs attention mechanisms, which enable better learning of interactions among input data. This has shown interesting results in many application fields. However, these advantages come with drawbacks such as increased training time and architectural complexity.
# Appendix

## Using sMAPE metrics for comparison between OASIS and set transformer

We present the results based on the sMAPE metric (Eqs. 37 and 38). This metric is more suitable for evaluating the forecasting of time series with outliers. In this case, it is not necessary to eliminate outliers; we follow the same approach as that presented in Sect. 4.2.2. According to the results presented in Tables 10, 11 and 12, we can see a close resemblance to the results presented using MAPE with outlier elimination.

$$ \text{sMAPE}_{\text{Total}}^{x_j} = \left\| 200 \times \frac{x_j(t)_p - x_j(t)_{gt}}{|x_j(t)_{gt}| + |x_j(t)_p|} \right\|_1 \quad \text{where} \quad t \in [\Delta t : T_f] \text{ and } j \in [1, N] \qquad (37) $$

$$ \text{sMAPE}_{\text{Total}}^{y_j} = \left\| 200 \times \frac{y_j(t)_p - y_j(t)_{gt}}{|y_j(t)_{gt}| + |y_j(t)_p|} \right\|_1 \quad \text{where} \quad t \in [\Delta t : T_f] \text{ and } j \in [1, N] \qquad (38) $$

## Detailed accuracy results for the comparison between deep set and set transformer

### Detailed model accuracy with respect to data quantity and noise levels

### Instance selection for the comparison between deep set and set transformer

The two Set Encoding methods are trained on different architectures independently, considering three levels of noise and three parameters of the Lorenz system (σ, β, and ρ), with variable time windows. Simultaneously, they are validated on the validation dataset. A check-pointing process is used if the validation error (Table 8) reaches a new minimum to save an instance of the architecture. The same process is conducted with different architectures for both methods, and the best instances are selected for each method, noise level, and each parameter of the Lorenz system. These instances are then used to compare the two methods. The selected architectures are presented in Tables 13 and 14.

### Deep set

The Deep Set model, as presented in Sect. 3, can be summarized into three parts:

* The first part is the encoding part, which takes as input the state variables of the Lorenz system [x(t), y(t), z(t)] and the time dimension t. It projects these state variables into a higher dimension using a simple MLP. The architecture of this MLP is defined in
(Table 13), where we have 5 layers with 320 perceptrons each and a ReLU activation function.

* The second part consists of choosing a pooling function over the time dimension. In this application, we found that applying the absolute value function followed by mean aggregation gives better results. The output of the aggregation function is a timeless feature vector.

* The last part consists of decoding the feature vector into the dynamics parameters. In this application, we used an MLP with the same architecture as for encoding.

## Set transformer

The Set Transformer architecture is more sophisticated compared to Deep Set, as it uses the attention mechanism, which means different encoding, pooling, and decoding approaches. Nevertheless, we can also summarize it as a three-part model:

* As for the Deep Set, this part consists of encoding the Lorenz system variables. For that, we used ISAB (Sect. 3.2), which performs an attention mechanism encoding within the time and the state variables of the Lorenz system. These variables are projected into a higher dimension $d = 45$, and then the attention mechanism is applied with respect to a trainable matrix $I$ called inducing points (Lee et al., 2019), with $m = 128$ rows and $d = 45$ columns. This means that ISAB will start by projecting the $n$ time snapshots into $m$ snapshots and then perform the attention mechanism with respect to the state variables of the Lorenz system. This method decreases the calculation time while providing similar results; more details are given in Lee et al. (2019). For the attention mechanism, we used $h = 40$ attention heads (Table 13). In this part, we apply a feedforward neural network to the output of the attention mechanism, adding more nonlinearity. For this application, the feedforward neural network is simply an MLP with 2 layers, 45 perceptrons each, and a ReLU activation function.

* The second part consists of applying pooling over the encoder feature vector, which is based on PMA, a Multi-Head Attention block that performs aggregation using the multi-head attention mechanism. For this, we used a pooling seed vector $S \in \mathbb{R}^{k \times d}$, where $d = 40$ and $k = 1$.

* Finally, we use an MLP with 2 layers, 40 perceptrons each, and a ReLU activation function, with a final layer of output dimension $d_{out} = 3$ that corresponds to the 3 parameters of the Lorenz system.

## Deep set architecture for the third application

The selected architecture is defined in Table 15. As detailed in Appendix A.4, we selected the best model using the validation datasets and a checkpoint that saves instances when the validation loss reaches a new minimum. The selected architecture consists of three blocks. In the first block, we use an MLP with 5 layers, 256 perceptrons, and a ReLU activation function. For the second part, we use a sum aggregation function. The encoded vector is
then decoded using three different MLPs with the same architecture as the encoder, each
corresponding to one of the characteristics (G, ratio, and α<sub>ref</sub>).

**Author Contributions** Not applicable.

**Funding** Open access funding provided by Centrale Nantes. The Jules Verne Institute of Research and Technology (IRT Jules Verne) (PhD PERFORM Program).

**Declarations**

**Conflict of interest** We affirm that there are no known Conflict of interest that have influenced this work.

**Ethics approval and consent to participate** We confirm that this research complies with authors ethical stand-
ards. All authors consent to participate in the submitted work.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

**References**

Ba, J. L., Kiros, J. R., & Hinton, G. E. (2016). Layer normalization. arXiv:1607.06450

Becker, S., Klein, M., Neitz, A., Parascandolo, G., & Kilbertus, N. (2023). Predicting ordinary differential equations with transformers. In: *International conference on machine learning* (pp. 1978–2002). PMLR.

Bhadriraju, B., Bangi, M. S. F., Narasingam, A., & Kwon, J.S.-I. (2020). Operable adaptive sparse identification of systems: Application to chemical processes. *AIChE Journal*, 66(11), 16980.

Bhadriraju, B., Kwon, J.S.-I., & Khan, F. (2021). OASIS-P: Operable adaptive sparse identification of systems for fault prognosis of chemical processes. *Journal of Process Control*, 107, 114–126.

Brunton, S. L., Proctor, J. L., & Kutz, J. N. (2016). Discovering governing equations from data by sparse identification of nonlinear dynamical systems. *Proceedings of the National Academy of Sciences*, 113(15), 3932–3937.

Brunton, S. L., Proctor, J. L., & Kutz, J. N. (2016). Sparse identification of nonlinear dynamics with control (SINDYc). *IFAC-PapersOnLine*, 49(18), 710–715.

Chartrand, R. (2011). Numerical differentiation of noisy, nonsmooth data. *International Scholarly Research Notices*, 2011.

Chen, Z., Liu, Y., & Sun, H. (2021). Physics-informed learning of governing equations from scarce data. *Nature Communications*, 12(1), 6136.

Chen, Y., Luo, Y., Liu, Q., Xu, H., & Zhang, D. (2022). Symbolic genetic algorithm for discovering open-form partial differential equations (SGA-PDE). *Physical Review Research*, 4(2), 023174.

Chicco, D., Warrens, M. J., & Jurman, G. (2021). The coefficient of determination r-squared is more informative than SMAPE, MAE, MAPE, MSE and RMSE in regression analysis evaluation. *PeerJ Computer Science*, 7, 623.

Cranmer, M. (2023). Interpretable machine learning for science with PySR and Symbolic Regression. jl arXiv:2305.01582

d'Ascoli, S., Becker, S., Mathis, A., Schwaller, P., & Kilbertus, N. (2023). Odeformer: Symbolic regression of dynamical systems with transformers. arXiv:2310.05573

Fasel, U., Kaiser, E., Kutz, J. N., Brunton, B. W., & Brunton, S. L. (2021). Sindy with control: A tuto-
rial. In: *2021 60th IEEE conference on decision and control (CDC)* (pp. 16–21). IEEE.

Gelß, P., Klus, S., Eisert, J., & Schütte, C. (2019). Multidimensional approximation of nonlinear dynamical systems. *Journal of Computational and Nonlinear Dynamics*, 14(6), 061006.
Hewamalage, H., Ackermann, K., & Bergmeir, C. (2023). Forecast evaluation for data scientists: Common pitfalls and best practices. *Data Mining and Knowledge Discovery*, 37(2), 788-832.

Kaiser, E., Kutz, J. N., & Brunton, S. L. (2018). Sparse identification of nonlinear dynamics for model predictive control in the low-data limit. *Proceedings of the Royal Society A*, 474(2219), 20180335.

Kamienny, P.-A., d'Ascoli, S., Lample, G., & Charton, F. (2022). End-to-end symbolic regression with transformers. *Advances in Neural Information Processing Systems*, 35, 10269-10281.

Kingma, D. P., & Ba, J. (2014). Adam: A method for stochastic optimization. arXiv:1412.6980

Lee, J., Lee, Y., Kim, J., Kosiorek, A., Choi, S., & Teh, Y. W. (2019). Set transformer: A framework for attention-based permutation-invariant neural networks. In: *International conference on machine learning* (pp. 3744-3753). PMLR.

Lejeune, C., Mothe, J., Soubki, A., & Teste, O. (2023). Data driven discovery of systems of ordinary differential equations using nonconvex multitask learning. *Machine Learning*, 112(5), 1523-1549.

Long, Z., Lu, Y., & Dong, B. (2019). PDE-Net 2.0: Learning PDEs from data with a numeric-symbolic hybrid deep network. *Journal of Computational Physics*, 399, 108925.

Makke, N., & Chawla, S. (2024). Interpretable scientific discovery with symbolic regression: A review. *Artificial Intelligence Review*, 57(1), 2.

Manucci, M., Aguado, J. V., & Borzacchiello, D. (2022). Sparse data-driven quadrature rules via 1 p-quasi-norm minimization. *Computational Methods in Applied Mathematics*, 22(2), 389-411.

Martius, G., & Lampert, C. H. (2016). Extrapolation and learning equations. arXiv:1610.02995

Narasingam, A., & Kwon, J.S.-I. (2017). Development of local dynamic mode decomposition with control: Application to model predictive control of hydraulic fracturing. *Computers & Chemical Engineering*, 106, 501-511.

North, J. S., Wikle, C. K., & Schliep, E. M. (2023). A review of data-driven discovery for dynamic systems. *International Statistical Review*, 91(3), 464-492.

Proctor, J. L., Brunton, S. L., & Kutz, J. N. (2016). Dynamic mode decomposition with control. *SIAM Journal on Applied Dynamical Systems*, 15(1), 142-161.

Quade, M., Abel, M., Nathan Kutz, J., & Brunton, S. L. (2018). Sparse identification of nonlinear dynamics for rapid model recovery. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 28(6).

Raissi, M., Perdikaris, P., & Karniadakis, G. E. (2019). Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. *Journal of Computational Physics*, 378, 686-707.

Rudin, C. (2019). Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead. *Nature Machine Intelligence*, 1(5), 206-215.

Rudy, S. H., Brunton, S. L., Proctor, J. L., & Kutz, J. N. (2017). Data-driven discovery of partial differential equations. *Science Advances*, 3(4), 1602614.

Silva, B. M., Champion, K., Quade, M., Loiseau, J.-C., Kutz, J. N., & Brunton, S. L. (2020). Pysindy: A python package for the sparse identification of nonlinear dynamics from data. arXiv:2004.08424

Sobol, I. (1967). The distribution of points in a cube and the accurate evaluation of integrals (in Russian). Zhurnal Vychislitel'noi Matematiki i Matematicheskoi Fiziki, 7, 784-802.

Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. *Journal of the Royal Statistical Society Series B: Statistical Methodology*, 58(1), 267-288.

Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., & Polosukhin, I. (2017). Attention is all you need. *Advances in Neural Information Processing Systems*, 30.

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., et al. (2020). SciPy 1.0: Fundamental algorithms for scientific computing in Python. *Nature Methods*, 17(3), 261-272.

Wangersky, P. J. (1978). Lotka-Volterra population models. *Annual Review of Ecology and Systematics*, 9(1), 189-218.

Zaheer, M., Kottur, S., Ravanbakhsh, S., Poczos, B., Salakhutdinov, R. R., & Smola, A. J. (2017). Deep sets. *Advances in Neural Information Processing Systems*, 30.

Zheng, P., Askham, T., Brunton, S. L., Kutz, J. N., & Aravkin, A. Y. (2018). A unified framework for sparse relaxed regularized regression: SR3. *IEEE Access*, 7, 1404-1423.

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
## Authors and Affiliations

Mouad Elaarabi¹,² · Domenico Borzacchiello¹ · Philippe Le Bot² · Yves L. E. Guennec² · Sebastien Comas-Cardona¹

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

Mouad Elaarabi
mouad.elaarabi@ec-nantes.fr; mouad.el-aarabi@irt-jules-verne.fr

Domenico Borzacchiello
domenico.borzacchiello@ec-nantes.fr

Philippe Le Bot
philippe.le-bot@irt-jules-verne.fr

Yves L. E. Guennec
yves.le-guennec@irt-jules-verne.fr

Sebastien Comas-Cardona
sebastien.comas@ec-nantes.fr

¹ Ecole Centrale Nantes, CNRS, GeM, UMR 6183, Nantes Université, 44321 Nantes, France

² IRT Jules Verne, Nantes Université, 44340 Bouguenais, France