Adaptive parameters identification for nonlinear dynamics using deep permutation invariant networks

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 $\alpha=0.5$, $\beta=0.025$, $\delta=0.5$, and $\gamma=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.

The dataset is generated through various combinations of control variable ($c$) and initial conditions $x_{0}$ and $y_{0}$, using the Sobol [31, 32] sampling method as given in Tab. 2. The data generation is performed using the LSODA [32] 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 [33]. We use a base of polynomials of degree 3 defined as :

where $\Theta(\widetilde{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.

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 Tab. 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 $\Xi$ and predicted $\hat{\Xi}$ parameters, the second part represents local dynamic loss by calculating the difference between the known derivative of $X$ and $Y$ and $\Theta(\widetilde{X},C)\hat{\Xi}$ where $\Theta(\widetilde{X},C)$ is know, 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{\Xi}$, 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 $\widetilde{X}_{t_{i}}$ and control variable $C_{t_{i}}$ as defined in Eq. 20 :

Where $N^{{}^{\prime}}$ is the time snapshots, $N^{{}^{\prime}}=1$ for OASIS framework and $N^{{}^{\prime}}=\Delta t=10$. Which will help in achieving better results for extrapolation.

To ensure fairness and balance in the comparison, we used hyper-parameters tunning 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.

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 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}^{j}_{t+\Delta t}$ which corresponds to the local dynamics of the following sub domain for the $j$ test dataset. It is then used with initial values $\widetilde{X}^{j}_{t_{i}+\Delta t}$ and control $C_{t_{i}+\Delta t}$ value to extrapolate using ’LSODA’ as ODE-Solver, the output of the solver is a times series $\widetilde{X}^{j}_{prediction}=[x_{j}(t)_{p},~y_{j}(t)_{p}]~~\textup{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 (Eq. 21 and Eq. 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 [34] 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 sections 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 (Eq. 23) and Eq. 24).

This metric is then used to compare the two methods. Based on [34] 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 $\textup{MAPE}^{x_{j}}$ and $\textup{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 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 Tabs. 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 domain dynamics. While the average of MAPE remains good for large extrapolation windows ($2\times\Delta t,~3\times\Delta t~~\text{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\widetilde{X}}{dt}-\Theta(\widetilde{X},C)\hat{\Xi}\right\|_{1}$ as given in Tab. 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 [35]. 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 challenging. Furthermore, even slight noise in the initial conditions may lead to inaccurate extrapolation or complete divergence.

We generate a dataset of 1024 Lorenz systems with variables parameters $\left[\sigma,~\rho,~\beta\right]$ , and initial conditions $\left[x_{0},~y_{0},~z_{0}\right]$ as described in Tab. 7 using LSODA [32] methods implemented in SciPy Python package. The parameters and initial conditions are sampled using the Sobol [31, 32] 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\left[0,0.02,0.05\right]$ as described in Eq. 25 :

Fig. 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 $\widetilde{\frac{\mathrm{d}X}{\mathrm{d}t}}~=~\left[\widetilde{\frac{\mathrm{d}x}{\mathrm{d}t}},~\widetilde{\frac{\mathrm{d}y}{\mathrm{d}t}},~\widetilde{\frac{\mathrm{d}z}{\mathrm{d}t}}\right]$ using total variation method [27] for noisy data and central difference (Eq. 26) for clean data.

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

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 section 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 $\left[\sigma_{k},~\rho_{k},~\beta_{k}\right]$. The identification process is defined as in Eq. 28 .

We implement a consistent training regime for both models (Deeps 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 $\left[\sigma_{k},~\rho_{k},~\beta_{k}\right]$. The training process is summarized in Tab. 8. We employ the Adam [36] 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 $\widetilde{\Xi}$ and $\hat{\Xi}$ denote the fitted parameters using SINDy (Eq. 28) based on denoised data and the predicted parameters. Note here that $\widetilde{\Xi}$ may differ from $\Xi$, 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.

The instances comparison are based on three features :

1. 1.

   Accuracy : the accuracy of selected instances is assessed using $R^{2}$-score metric [37] defined as :

   $$R^{2}=1-\frac{\sum_{i}(y_{i}-\hat{y}_{i})^{2}}{\sum_{i}(y_{i}-\bar{y}_{i})^{2}},$$

   (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{\Xi}$. 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}}~~~\textup{where}~~F_{T}=\sum_{i}1/n_{i},$$

   (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 to different noise levels. We present the $R^{2}_{n_{i}}$ with respect to the quantity of data in appendix A.2 (Fig. 18, 19 and 20).

   

2. 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.6Ghz 48-Core Processor, making these models more suitable for applications like the extrapolation of state variables [22] or online control [3]. 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. 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.

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) :

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, $\mathrm{G}=\mathrm{center_{abnormality}}$ which presents the center of abnormality and $\mathrm{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}$:

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 (Tab. 9) :

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 $i^{th}$ dataset is defined as (Eq. 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 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):

, which is composed of four losses: $\mathscr{L}_{\text{center}}$, $\mathscr{L}_{\text{ratio}}$, $\mathscr{L}_{\alpha}$ and $\mathscr{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]$:

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 4.3 :

  $$R^{2}_{w.a}=\sum_{i}\frac{1/n_{i}}{F_{T}}R^{2}_{n_{i}}~~~\textup{where}~~F_{T}=\sum_{i}1/n_{i},$$

  (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 (Fig. 17, and. 21), we can see that the model can easily predict the $\mathrm{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 $\mathrm{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 $\mathscr{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.