An Overview of the Euler-Type Universal Numerical Integrator (E-TUNI): Applications in Non-Linear Dynamics and Predictive Control

1. Introduction

One of the first works on artificial neural networks is credited to McCulloch and Pitts in [1] and is a work from 1943. In these eighty years of development of artificial intelligence and artificial neural networks, many interesting things have emerged in this area of human knowledge. In [2], a 1957 work credited to the great Russian mathematician Kolmogorov, it is demonstrated that any n-dimensional function can be represented by linearly combining one-dimensional non-linear functions. In [3], Kolmogorov’s work from 1957 was extensively recognized, and the term Kolmogorov neural networks was introduced in the literature. In [4], the classic back-propagation algorithm is proposed for training Multi-Layer Perceptron (MLP) networks with at least one inner layer. In [5,6], it is demonstrated, independently, that MLP networks are universal approximators of functions. A very detailed summary of the main scientific knowledge achieved by artificial neural networks in the 20th century can be found at [7].

One of the great applications of artificial neural networks is in the modeling of non-linear dynamic systems in the real world. Thus, next, a bibliographical review is reported, only on this theme, considering only neural networks with feed-forward architecture trained through supervised learning.

In [8], the concept of Universal Numerical Integrator (UNI) is defined. A UNI is nothing more than the coupling of a conventional numerical integrator (e.g., Euler, Runge-Kutta, predictive-corrector, among others) with a universal approximator of functions (e.g., MLP, SVM, RBF, and wavelet networks, fuzzy inference system, paraconsistent inference system, among others). However, it is essential to note that this definition only took into account neural networks with feed-forward architecture trained exclusively through supervised learning using input/output training patterns.

Also in [8], an interesting classification for UNIs is given. According to this reference, the UNIs can be divided into three large classes, thus giving rise to three distinct methodologies, namely: (i) the NARMAX model, (ii) the mean derivative methodology (e.g., E-TUNI), and ( iii) the instantaneous derivatives methodology (e.g., Runge-Kutta neural network, Adams-Bashforth neural network, predictive-corrector neural network, among others). In [8], these methodologies are also presented in some depth, and a detailed table comparing these three methodologies is presented. It is important to note that in the classification of UNIs, the NARMAX model was included. The reason for this is that, although in the NARMAX model, the neural network is not included in the numerical integrator structure, the neural network behaves as if it were a numerical integrator.

In [9,10,11], they are all classic 20th-century books on conventional numerical integrators. The term "conventional" is used here to inform the reader that such structures do not have a universal approximator of functions coupled to them. In [12,13], the NARMAX methodology is discussed in detail using artificial neural networks. In [14] is Leonard Euler’s original work on the first-order integrator that bears his name. This is a 1768 work and initially uses instantaneous derivatives rather than mean derivatives. Everyone knows that the first-order Euler integrator, although mathematically quite simple, is relatively imprecise. However, when the instantaneous derivative functions are exchanged for the mean derivative functions, the Euler integrator becomes as accurate as any other higher-order integrator. In this case, such structures are given the special name of Euler-Type Universal Numerical Integrator (E-TUNI).

Furthermore, a correction must be given here. In [8], the E-TUNI is called Euler Neural Networks and is a 2019 reference. However, in the 2009 reference [15], the term Euler Neural Network is also used, but in a different context. Therefore, to avoid any confusion, we decided to name the first-order UNI of E-TUNI. That said, the references [16,17,18] refer directly or indirectly to E-TUNI. However, since reference [16] is relatively old, the term E-TUNI is referred to there only as Euler Neural Integrator or Mean Derivative Methodology.

In [16], the E-TUNI was coupled to a neural network with an MLP architecture and applied in a predictive control structure. In [17], a qualitative analysis on the design of a Universal Numerical Integrator (UNI) is carried out in detail. In this reference, the E-TUNI is compared with the NARMAX Model and the Runge-Kutta Neural Network (RKNN). Reference [18] is a continuation of reference [17]. In reference [18], a quantitative analysis of the design of a universal numerical integrator is presented in some mathematical depth. In [19], an important work is done involving E-TUNI in a backward integration process. It is worth noticing that E-TUNI, as well as the NARMAX model, are dynamic models that are naturally discrete and not continuous.

The first work we find in the UNI literature is credited to Wang and Lin in [20]. This reference is from 1998, and nothing before that has been found in the literature. This work is about the Runge-Kutta Neural Network (RKNN). In [21], the Runge-Kutta neural network is used in control. It is interesting to note that there is a gap of more than 20 years between the [20] and [21] references. It is also essential to note that the Predictive-Corrector Neural Network (PCNN) still does not exist in the literature. In [22], it is an essential and current work on neural differential equations. In [23], it is another work on the Runge-Kutta neural network applied in control.

It is vital to know that both the NARMAX model and the mean derivatives methodology (E-TUNI) are inherently of fixed integration step $\Delta t$. This property means that if it is desired to change the $\Delta t$ integration step in these two methodologies, then it is necessary to do a new neural training. However, in the instantaneous derivatives methodology (RKNN, ABNN, PCNN, among others) the integration step $\Delta t$ can be changed, within certain limits imposed by the order of the referred integrator used, without having to carry out a new training [8,17,18,20]. Finally, some small theoretical models of non-linear ordinary differential equations can be taken from [24] to test any UNI.

The original contributions of this paper are as follows: (i) it presents an overview of how E-TUNI works to solve dynamics and control problems, (ii) it presents three studies of practical technological applications of the proposed model (e.g., including the orbit transfer problem and the non-linear simple pendulum problem), (iii) it presents a practical example of how to generate an approximate continuous solution for the E-TUNI model (see the nonlinear simple pendulum problem), which is naturally discrete in solution, and (iv) it presents a correct proof of the general expression for E-TUNI, since in [16] it is not correct. Furthermore, the E-TUNI has its origins in [16].

Accordingly, the remainder of this article is organized as follows. Section 2 describes in detail the meaning of the main discrete and continuous variables, which are used to describe the proposed E-TUNI model. Section 3 presents a concise yet comprehensive theoretical overview of the Euler-Type Universal Numerical Integrator (E-TUNI), outlining its foundational principles and key mathematical formulations. Section 4 details three numerical simulations that serve to validate the proposed methodology and demonstrate its practical effectiveness. Finally, Section 5 summarizes the main conclusions, highlighting the contributions and potential future directions of this research.

2. Preliminaries and Symbols Used

To facilitate a deeper understanding of the theoretical framework developed in this paper, we present a complete definition of all symbols and variables used throughout the work. The proposed methodology is named Euler Type Universal Numerical Integrator (E TUNI), and it is based on a forward integration scheme that operates in a discrete-time setting. This approach is designed to approximate the behavior of continuous-time autonomous ordinary differential equations by means of a discrete formulation. For clarity, all notation used in this context is listed below and classified into two main categories: (i) variables defined in continuous time and (ii) variables defined in discrete time.

(i) Continuous Variables

• $\dot{y}=f\left(y\right)$… System of continuous differential equations.
• $y={[y_1{y}_2\dots y_n]}^T$… State Variables.
• $f\left(y\right)={[f_1\left(y\right)f_2\left(y\right)\dots f_n\left(y\right)]}^T$… Instantaneous derivative functions.
• $y_j^i\left(t\right)=g_j^i\left(t\right)$… Particular continuous and differentiable curve of a family of solution curves of the dynamical system $\dot{y}=f\left(y\right)$.
• ${\dot{y}}_j^i\left(t\right)={\dot{g}}_j^i\left(t\right)$… First derivative of $y_j^i\left(t\right)$.

(ii) Discrete Variables

• ${}^{k}y^i=y^i(t_0+k·\Delta t)$… Vector of state variables at time $t_k$.
• ${}^{k}y_j^i$… Scalar state variable for $j=1,2,\dots ,n$ at time $t_k$. It is a generic discretization point of the state variables generated by the integers i, j, and k.
• n… Total number of state variables.
• ${}^{k}u=u(t_0+k·\Delta t)$… Vector of control variables at time $t_k$.
• ${}^{k}u_j$… Scalar control variable for $j=1,2,\dots ,m$ at time $t_k$.
• m… Total number of control variables.
• ${}^{k+1}y^i=y^i[t_0+(k+1)·\Delta t]$… Exact vector of state variables at time $t_{k+1}$.
• ${}^{k+1}y_j^i$… Exact scalar state variable for $j=1,2,\dots ,n$ at time $t_{k+1}$.
• ${}^{k+1}\widehat{y}^i={\widehat{y}}^i[t_0+(k+1)·\Delta t]$… Estimated Vector of state variables by UNI or E-TUNI at time $t_{k+1}$.
• ${}^{k+1}\widehat{y}_j^i$… Scalar state variable estimated by UNI or E-TUNI for $j=1,2,\dots ,n$ at time $t_{k+1}$.
• ${}^{k+1}\tilde{y}^i$… Estimated Vector of state variables when using only the integrator and without using the neural network at time $t_{k+1}$.
• $ta{n}_{\Delta t}{}^k\alpha ^i=ta{n}_{\Delta t}{}^k\alpha ^i={[ta{n}_{\Delta t}{}^k\alpha _1^{i}ta{n}_{\Delta t}{}^k\alpha _2^i\dots ta{n}_{\Delta t}{}^k\alpha _n^i]}^T$… Exact vector of positive mean derivative functions at time $t_k$.
• $ta{n}_{\Delta t}{}^k\alpha _j^i=ta{n}_{\Delta t}{}^k\alpha _j^i={\displaystyle \frac{{}^{k+1}y_j^i-{}^{k}y_j^i}{\Delta t}}$… Scalar positive mean derivative functions for $j=1,2,\dots ,n$ at time $t_k$.
• $ta{n}_{\Delta t}{}^k\widehat{\alpha }^i={[ta{n}_{\Delta t}{}^k\widehat{\alpha }_1^{i}ta{n}_{\Delta t}{}^k\widehat{\alpha }_2^i\dots ta{n}_{\Delta t}{}^k\widehat{\alpha }_n^i]}^T$… Estimated vector of positive mean derivative functions by the E-TUNI at time $t_k$.
• $tan{}^k\theta ^i={[ta{n}^k{\theta }_1^{i}ta{n}^k{\theta }_2^i\dots ta{n}^k{\theta }_n^i]}^T$… Vector of positive instantaneous derivatives at time $t_k$.
• $ta{n}^k{\theta }_j^i=\underset{\Delta t\to 0}{lim}{\displaystyle \frac{{}^{k+1}y_j^i{-}^k{y}_j^i}{\Delta t}}$… Scalar positive instantaneous derivative for $j=1,2,\dots ,n$ at instant $t_k$.
• $t_k$… Time instant $t_k=t_0+k·\Delta t$.
• $t_{k+1}$… Time instant $t_{k+1}=t_0+(k+1)·\Delta t$.
• $\Delta t$… Integration step.
• i… Over-index that enumerates a particular curve from the family of curves of the dynamical system to be modelled ($i=1,2,\dots ,q$).
• j… Under-index that enumerates the state and control variables.
• k… Over-index that enumerates the discrete time instants ($k=1,2,\dots ,r$).
• r… Total number of horizons of the time variable.
• q… Total number of curves from the family of curves curves of the dynamic system to be modelled.
• $t_k^{*}$… Instant of time within the interval $[t_k,t_{k+1}]$ as a result of the Differential Mean Value Theorem (see Theorem 1).
• $t_k^x$… Instant of time within the interval $[t_k,t_{k+1}]$ as a result of the Integral Mean Value Theorem (see Theorem 2).

To interpret the notation used in this work, it is essential to distinguish between the variables ${}^{k+1}y^i$, ${}^{k+1}\tilde{y}^i$, and ${}^{k+1}\widehat{y}^i$. The vector variable ${}^{k+1}y^i$ is the exact value of the state variables at time $t_{k+1}=t_0+(k+1)·\Delta t$. The vector variable ${}^{k+1}\tilde{y}^i$ is the estimated value of the state variables using only the numerical integrator. Finally, the vector variable ${}^{k+1}\widehat{y}^i$ is the estimated value of the state variables using the numerical integrator coupled to an artificial neural network. Figure 1 follows this convention and helps to understand better the notations adopted here. Additionally, it is important to understand the meaning of the auxiliary variables i, j, and k. This is accomplished in the next paragraph.

A key point in this notation is the interpretation of the over-indexes k and i and the under-index j in the variables ${}^{k}y_j^i$ and $ta{n}_{\Delta t}{}^k\alpha _j^i$. When the auxiliary variables i, j, and k are used simultaneously, they uniquely identify the secant ($ta{n}_{\Delta t}{}^k\alpha _j^i$) at the point ${}^{k}y_j^i$. Note that the over-index k indicates the time instant $t_k=t_0+k·\Delta t$ of the secant, the under-index j indicates the state variable in question, where $j=1,2,\dots ,n$, and the over-index i$(i=1,2,\dots ,q)$ indicates the particular curve where the respective secant is located, from the family of possible curves of the system of differential equations considered. The value of q can be as large as desired. The larger the value of q, the more different curves the neural network trains on, and the better its generalization. This convention is quite useful to fully understand the formal proof of the general expression of E-TUNI that is performed in section 3.2.

3. Mathematical Development

In this section, we present a concise and complete description of the Euler-Type Universal Numerical Integrator (E-TUNI). To this end, we present a formal mathematical proof for the general expression that governs E-TUNI’s operation to generate discrete solutions for autonomous nonlinear dynamical systems governed by ordinary differential equations. Additionally, we also present the correct way to use E-TUNI in a predictive control framework. We conclude this section by obtaining an approximate mathematical expression for E-TUNI to provide a continuous solution, rather than a discrete solution, for autonomous dynamical systems.

3.1. Basic Mathematical Development of E-TUNI

We provide a brief mathematical description of E-TUNI below. Then, in the following sub-section, we perform a formal mathematical demonstration of the general expression of the first-order Euler Integrator designed with mean derivative functions. So, by definition, the secants or mean derivatives $ta{n}_{\Delta t}{}^k\alpha _j^i$ for $j=1,2,\dots ,n$ between the points ${}^{k+1}y_j^i$ and ${}^{k}y_j^i$ are given by:

$$ta{n}_{\Delta t}{}^k\alpha _j^i={\displaystyle \frac{{}^{k+1}y_j^i-{}^{k}y_j^i}{\Delta t}}$$

(1)

where, ${}^{k+1}y_j^i=y_j^i[t+(k+1)·\Delta t]$ is the forward state of the dynamic system, ${}^{k}y_j^i=y_j^i[t+k·\Delta t]$ is the present state of the dynamical system, the over-index k on the left indicates the instant k, the over-index i on the right suggests a discretization of the continuum, the sub-index j on the right indicates the j-th state variable, n is the total number of state variables and $\Delta t$ is the integration step. We talk more about the over-index i later in this article.

A geometric and intuitive difference between the mean and the instantaneous derivative functions is shown in Figure 1. In line with the Figure 1 and Figure 2, the E-TUNI is entirely based on the concept of mean derivative functions and not on the idea of instantaneous derivative functions. As will be seen throughout this article, this change is quite significant when incorporated adequately into the Euler-type first-order integrator. Furthermore, the mean derivative functions can also be obtained from supervised training with input/output patterns using a neural network with any feed-forward architecture (MLP, RBF, SVM, among others).

It is also essential to note that it is possible to train E-TUNI in two different ways, namely: a) through the direct approach or b) through the indirect or empirical approach. In the direct approach, the neural network is trained decoupled from the Euler-type integrator, while in the indirect or empirical approach, the neural network is trained coupled to the structure of the first-order Euler integrator. For this reason, in the indirect approach, the back-propagation algorithm needs to be modified slightly. In the references [8,17,18] this is explained in more detail.

In this way, be a non-linear dynamic system of n simultaneous first-order equations with dependent variables $y_1,y_2,\dots ,y_n$. If each of these variables satisfies a given initial condition for the same value a of t, then we have an initial value problem for a first-order system, and we can write:

$$\left\{\begin{array}{cc}{\dot{y}}_1=f_1(y_1,y_2,\dots y_n),& y_1\left(a\right)={\eta }_1\\ {\dot{y}}_2=f_2(y_1,y_2,\dots y_n),& y_2\left(a\right)={\eta }_2\\ ⋮& ⋮\\ {\dot{y}}_n=f_n(y_1,y_2,\dots y_n),& y_n\left(a\right)={\eta }_n\end{array}\right.$$

(2)

In [14], the mathematician Leonard Euler himself proposed in 1768 the first numerical integrator, in the history of mathematics, to approximately solve non-linear dynamical systems governed by the equation (2). This solution, which is well-known to everyone, is given by:

$$\left\{\begin{array}{c}{}^{k+1}y_1^i\cong tan{}^k\theta _1^i·\Delta t+{}^{k}y_1^i\\ {}^{k+1}y_2^i\cong tan{}^k\theta _2^i·\Delta t+{}^{k}y_2^i\\ ⋮\\ {}^{k+1}y_n^i\cong tan{}^k\theta _n^i·\Delta t+{}^{k}y_n^i\end{array}\right.$$

(3)

where, $tan{}^k\theta _j^i$ for $j=1,2,\dots ,n$ are instantaneous derivative functions, according to the graphical notation shown in Figure 1. On the other hand, there is an attempt in [18] to prove mathematically that if you exchange the instantaneous derivatives in (3) for the mean derivatives, then the solution proposed by Leonard Euler becomes accurate, i.e.,

$$\left\{\begin{array}{c}{}^{k+1}y_1^i=ta{n}_{\Delta t}{}^k\alpha _1^i·\Delta t+{}^{k}y_1^i\\ {}^{k+1}y_2^i=ta{n}_{\Delta t}{}^k\alpha _2^i·\Delta t+{}^{k}y_2^i\\ ⋮\\ {}^{k+1}y_n^i=ta{n}_{\Delta t}{}^k\alpha _n^i·\Delta t+{}^{k}y_n^i\end{array}\right.$$

(4)

Comparing equations (3) and (4), it is observed that the instantaneous derivative functions do not depend on the integration step, but the mean derivative functions do. Additionally, the equations in (4) can also be expressed in a more compact notation, given by:

$${}^{k+1}y^i=ta{n}_{\Delta t}{}^k\alpha ^i·\Delta t+{}^{k}y^i$$

(5)

where, ${}^{k+1}y^i=$[${}^{k+1}y_1^i$${}^{k+1}y_2^i$…${}^{k+1}y_n^i$${]}^T$, $ta{n}_{\Delta t}{}^k\alpha ^i=$[$ta{n}_{\Delta t}{}^k\alpha _1^i$$ta{n}_{\Delta t}{}^k\alpha _2^i$…$ta{n}_{\Delta t}{}^k\alpha _n^i$${]}^T$ e ${}^{k}y^i=$[${}^{k}y_1^i$${}^{k}y_2^i$…${}^{k}y_n^i$${]}^T$. The generalization of the mean derivatives methodology to multiple backward inputs and/or multiple forward outputs is relatively easy. This expression can be obtained by following the equation:

$${}^{k+p}y_j^i=\sum _{m=0}^{p-1}ta{n}_{\Delta t}{}^{k+m}\alpha _j^i·\Delta t+{}^{k}y_j^i$$

(6)

where p is the number of backward and/or forward instants and $j=1,2,\dots ,n$. For example, if $p=$$p_1+p_2$ then it is possible to design an E-TUNI with $p_1$ inputs in the role of backward mean derivatives and $p_2$ outputs in the role of forward mean derivative. In this way, the reader is free to choose the values of $p_1$ and $p_2$ as long as the previous equality is satisfied. However, it should be noted that the first input of the neural network must be an absolute value, not a relative one.

On the other hand, notice that if the reader tries to compare the NARMAX model with the E-TUNI structure, one can see that the former is very similar to the latter. In the NARMAX model, for example, the output of the universal approximator of functions is the forward instant $y_j(t+\Delta t)$. However, in the E-TUNI model, the output of the universal approximator of functions is the mean derivative function $ta{n}_{\Delta t}{}^k\alpha _j^i$ at the present instant. This description is the only practical difference between these two methodologies. However, the E-TUNI may be a little more computationally efficient than the NARMAX model, as explained in the following paragraph.

Considering the direct approach to training the mean derivative functions required by the E-TUNI structure, it is then possible to perform a theoretical analysis of the local error committed by this first-order universal numerical integrator. Thus, let the exact value ${}^{k+1}\overline{y}_j^i$ and the estimated value ${}^{k+1}\widehat{y}_j^i$ be obtained, respectively, by the equations (7) and (8) of a given solution of a generic dynamical system.

$${}^{k+1}\overline{y}_j^i=ta{n}_{\Delta t}{}^k\alpha _j^i·\Delta t+{}^{k}y_j^i$$

(7)

$${}^{k+1}\widehat{y}_j^i=(ta{n}_{\Delta t}{}^k\alpha _j^i\pm e_m)·\Delta t+{}^{k}y_j^i$$

(8)

where $e_m$ is the mean absolute error of the output variables of the universal approximator of functions used to learn the mean derivative functions. Thus, if the equation (7) is subtracted from the equation (8) and the result of this subtraction is squared, we have:

$${\left(\begin{array}{c}{}^{k+1}\overline{y}_j^i-{}^{k+1}\widehat{y}_j^i\end{array}\right)}^2={\Delta t}^2·e_m^2$$

(9)

The equation (9) states that the local squared error made by E-TUNI can dampen the squared error of training the universal approximator of functions, used in neural training, if $0<\Delta t<1$. If $\Delta t>1$, the local error will be amplified. However, for a more accurate analysis of the global training error, further studies are needed.

The equation (9) is fundamental, as it partially explains why training the E-TUNI with a mean square error $e_m^2$, greater than that obtained by training the NARMAX model can also yield good results from estimation and potentially surpass those obtained in the NARMAX methodology. However, as stated earlier, this is true only if $\Delta t<1$.

Therefore, an E-TUNI neural integrator can be used for supervised training of a generic plant. Additionally, this plant can be used in a predictive control framework. Thus, the plant’s neural model can be used, as a model of the system’s internal response to derive a smooth control policy. This control is then able to track a reference trajectory by minimizing a finite-horizon quadratic functional, given by [16]:

$$J=\{\sum _{j=1}^m{[y_r\left(t_j\right)-\widehat{y}\left(t_j\right)]}^T·r_y^{-1}\left(t\right)·[y_r\left(t_j\right)-\widehat{y}\left(t_j\right)]+$$

(10)

$$\sum _{j=0}^{m-1}{[u\left(t_j\right)-u\left(t_{j-1}\right)]}^T·r_u^{-1}\left(t\right)·[u\left(t_j\right)-u\left(t_{j-1}\right)]\}/2$$

where, $y_r\left(t_j\right)$ is the reference trajectory at the instant $t_j$; m is the number of horizons ahead; $r_y^{-1}\left(t_j\right)$ and $r_u^{-1}\left(t_j\right)$ are positive definite weight matrices; $\widehat{y}\left(t_j\right)$ is the output of the previously trained E-TUNI. Thus, in [16] it is shown that it is necessary to know the partial derivatives ${\displaystyle \frac{\partial {}^{k+q}y^i}{\partial {}^{k}u}}$ to solve the problem of optimization given by the equation iteratively (10) through the Kalman filter. Also, in [16] these partial derivatives can be calculated as follows:

$${\displaystyle \frac{\partial {}^{k+q}y^i}{\partial {}^{k}u}}=\left(\begin{array}{c}{\displaystyle \frac{\partial ta{n}_{\Delta t}{}^{k+q-1}\alpha ^i}{\partial {}^{k+q-1}y^i}}·\Delta t+I\end{array}\right)·{\displaystyle \frac{\partial {}^{k+q-1}y^i}{\partial {}^{k}u}}$$

(11)

$$toq=2,3,...,m$$

where,

$${\left[\begin{array}{c}{\displaystyle \frac{{\partial }^{k+q-1}{y}^i}{{\partial }^{k}u}}\end{array}\right]}_{n_y\times n_u}=$$

$$\Delta t·\left[\begin{array}{cccc}{\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-2}{\alpha }_1^i}{{\partial }^k{u}_1}}& {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-2}{\alpha }_1^i}{{\partial }^k{u}_2}}& \dots & {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-2}{\alpha }_1^i}{{\partial }^k{u}_{n_u}}}\\ {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-2}{\alpha }_2^i}{{\partial }^k{u}_1}}& {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-2}{\alpha }_2^i}{{\partial }^k{u}_2}}& \dots & {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-2}{\alpha }_2^i}{{\partial }^k{u}_{n_u}}}\\ ⋮& ⋮& \ddots & ⋮\\ {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-2}{\alpha }_{n_y}^i}{{\partial }^k{u}_1}}& {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-2}{\alpha }_{n_y}^i}{{\partial }^k{u}_2}}& \dots & {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-2}{\alpha }_{n_y}^i}{{\partial }^k{u}_{n_u}}}\end{array}\right]$$

(12)

$${\left[\begin{array}{c}{\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-1}{\alpha }^i}{{\partial }^{k+q-1}{y}^i}}\end{array}\right]}_{n_y\times n_y}=$$

$$\left[\begin{array}{cccc}{\displaystyle \frac{\partial ta{n}_{\Delta t}{}^{k+q-1}\alpha _1^i}{\partial {}^{k+q-1}y_1^i}}& {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-1}{\alpha }_1^i}{{\partial }^{k+q-1}{y}_2^i}}& \dots & {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-1}{\alpha }_1^i}{{\partial }^{k+q-1}{y}_{n_y}^i}}\\ {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-1}{\alpha }_2^i}{{\partial }^{k+q-1}{y}_1^i}}& {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-1}{\alpha }_2^i}{{\partial }^{k+q-1}{y}_2^i}}& \dots & {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-1}{\alpha }_2^i}{{\partial }^{k+q-1}{y}_{n_y}^i}}\\ ⋮& ⋮& \ddots & ⋮\\ {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-1}{\alpha }_{n_y}^i}{{\partial }^{k+q-1}{y}_1^i}}& {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-1}{\alpha }_{n_y}^i}{{\partial }^{k+q-1}{y}_2^i}}& \dots & {\displaystyle \frac{\partial ta{n}_{\Delta t}^{k+q-1}{\alpha }_{n_y}^i}{{\partial }^{k+q-1}{y}_{n_y}^i}}\end{array}\right]$$

(13)

$${\left[\begin{array}{c}{\displaystyle \frac{\partial {}^{k+1}y^i}{\partial {}^{k}u}}\end{array}\right]}_{n_y\times n_u}=$$

$$\Delta t·\left[\begin{array}{cccc}{\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _1^i}{\partial {}^{k}u_1}}& {\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _1^i}{\partial {}^{k}u_2}}& \dots & {\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _1^i}{\partial {}^{k}u_{n_u}}}\\ {\displaystyle \frac{\partial ta{n}_{\Delta t}^k{\alpha }_2^i}{\partial {}^{k}u_1}}& {\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _2^i}{\partial {}^{k}u_2}}& \dots & {\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _2^i}{\partial {}^{k}u_{n_u}}}\\ ⋮& ⋮& \ddots & ⋮\\ {\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _{n_y}^i}{\partial {}^{k}u_1}}& {\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _{n_y}^i}{\partial {}^{k}u_2}}& \dots & {\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _{n_y}^i}{\partial {}^{k}u_{n_u}}}\end{array}\right]$$

(14)

In the equations (11) to (14), we have that $n_y=n$ is the total number of state variables and $n_u$ is the total number of control variables. Furthermore, to derive these equations, it was assumed that the neural network, which learns the mean derivative functions, was designed with only one late input for the state and control variables, and also only one forward output for the mean derivatives.

So, for example, if a neural network with an MLP architecture is used, the calculation of the partial derivatives, necessary for the use of the gradient training algorithm, can be obtained as follows [16]:

$${\left[\begin{array}{cccc}ta{n}_{\Delta t}{}^k\alpha _1^i,& ta{n}_{\Delta t}{}^k\alpha _2^i,& \dots ,& ta{n}_{\Delta t}{}^k\alpha _{n_y}^i\end{array}\right]}^T=$$

(15)

$$f_{NN}^{MLP}={\left[\begin{array}{ccccc}{y}_1\left(t\right),\dots ,& y_{n_y}\left(t\right),& u_1\left(t\right),& \dots ,& u_{n_u}\left(t\right)\end{array}\right]}^T$$

where,

$${\displaystyle \frac{\partial y^l}{\partial {\overline{y}}^k}}={\displaystyle \frac{\partial y^l}{\partial {\overline{y}}^{k+1}}}·W^{k+1}·I_{f^{\prime }\left({\overline{y}}^k\right)}tok=l-1,l-2,\dots ,1$$

$${\displaystyle \frac{\partial y^l}{\partial {\overline{y}}^l}}=I_{f^{\prime }\left({\overline{y}}^l\right)}$$

$$I_{f^{\prime }\left({\overline{y}}^l\right)}=\left[\begin{array}{cccc}{I}_{f^{\prime }\left({\overline{y}}_1^l\right)}& 0& \dots & 0\\ 0& I_{f^{\prime }\left({\overline{y}}_2^l\right)}& \dots & 0\\ 0& 0& \ddots & 0\\ 0& 0& \dots & I_{f^{\prime }\left({\overline{y}}_{n_k}^l\right)}\end{array}\right]$$

It is essential to note that l is a generic layer of the MLP network. So $y^l$ are the output values of the l layer, and when l is the last layer, then these outputs necessarily are the mean derivative functions. The $f^{\prime }(·)$ functions can be any sigmoid function. Furthermore, if another neural architecture is used, for example, RBF networks or Wavelets, it will only change the equation (15). The equations (11), (12), (13) and (14) remain unchanged. The reason for this is that the equations from (11) to (14) refer exclusively to the type of integrator used, which, in this case, is necessarily the E-TUNI.

On the other hand, the equation (15) refers exclusively to the type of feed-forward neural network used, which, in this case, is necessarily the MLP network. Thus, the equation (15) is nothing more than an iterative version of the back-propagation algorithm, which calculates the partial derivatives from the output of the MLP network, concerning its inputs and not concerning the synaptic weights.

It is also worth noticing that, to fully understand the equations from (11) to (14), it is convenient to consult the reference [16], as it is necessary to make a temporal chain of several first-order Euler integrators, which work exclusively with mean derivative functions. This fact is because the horizon m has, in general, temporal advances in an amount greater than the total number of delayed inputs of the E-TUNI used.

3.2. Correct Mathematical Demonstration of the E-TUNI General Expression

In this section, we formally demonstrate the general mathematical expression of E-TUNI, which is nothing more than a first-order Euler-type integrator working with mean derivative functions. This development is done here because, as previously stated, there is a demonstration error in the general E-TUNI expression presented in reference [16]. So, the starting point is Figure 3.

So, for the reader to understand this figure, we explain it, starting from the point ${}^{k}y_j^i$. The point ${}^{k}y_j^i$ is an instant of the solution of the considered autonomous dynamic system. The variable j means the jth state variable of the considered dynamic system, where $j=1,2,\dots ,n$, that is, there are a total of n state variables. The variable i represents the i-th curve of the family of curves confined in the interval $[y_j^{min}\left(t_0\right),$$y_j^{max}\left(t_0\right)]$ and which is the solution of the considered dynamic system. By the continuum hypothesis, one can have infinite curves ($i=1,2,\dots ,\infty$). The variable k represents the k-th instant of time, that is, $t^k=t_0+k·\Delta t$. In the case of the variable k, if the dynamical system solution does not come out of the region of interest [$y_j^{min}\left(t_0\right),$$y_j^{max}\left(t_0\right)$] and the system is autonomous then the variable k can also go to infinity ($k=1,2,\dots ,\infty$). Also, when we write $y_j^i\left(t_k^x\right)$ or $y_j^i\left(t_k^{*}\right)$ it means that $t_k<t_k^x<t_{k+1}$ and $t_k<t_k^{*}<t_{k+1}$.

This notation for ${}^{k}y_j^i$ is sufficient to uniquely map it to its respective region of interest, where its associated secant characteristic $ta{n}_{\Delta t}{}^k\alpha _j^i$ is confined. So, in this case, for example, we can say that for the state ${}^{k}y_j^i=y_j^i(t_0+k·\Delta t)$ we have its associated characteristic secant given by $ta{n}_{\Delta t}{}^k\alpha _j^i$.

Thus, the secant $ta{n}_{\Delta t}{}^k\alpha _j^i$ is associated with the state variable $y_j$ with a curve i specific to the family of solutions, in this case, by $y_j^i\left(t\right)$ and the instant $t_k$. The fact that the secant is associated with the instant $t_k$ means that it is confined to the closed interval $[t_k,t_{k+1}]$. Thus, the secant $ta{n}_{\Delta t}{}^k\alpha _j^i$ starts at $t_k$ and ends at $t_{k+1}$.

Furthermore, the states $y_j^i\left(t_k^x\right)$ and $y_j^i\left(t_k^{*}\right)$ are special states, where both $t_k^x$ and $t_k^{*}$ are also confined to the closed interval $[t_k,t_{k+1}]$. The instants $t_k^x$ as $t_k^{*}$ are directly associated, respectively, with the integral mean value theorem and the differential mean value theorem [25,26]. So, having made these preliminary considerations, we begin now with this demonstration in a precise manner. Thus, let the following autonomous system of non-linear differential equations be given by,

$$\dot{y}=f\left(y\right)$$

(16)

$$y={\left[\begin{array}{cccc}{y}_1& y_2& \dots & y_n\end{array}\right]}^T$$

(17)

$$f\left(y\right)={\left[\begin{array}{cccc}{f}_1\left(y\right)& f_2\left(y\right)& \dots & f_n\left(y\right)\end{array}\right]}^T$$

(18)

Consider also, by definition, $y_j^i=y_j^i\left(t\right)$ for $j=1,2,$$\dots ,n$ a particular trajectory for a family of solutions to the system of differential equations $\dot{y}=f\left(y\right)$ passing through $y_j^i\left(t_0\right)$ at instant $t_0$, i.e., initializing from a domain of interest $[y_j^{min}\left(t_0\right),$$y_j^{max}\left(t_0\right){]}^n$, where $y_j^{min}\left(t_0\right)$ and $y_j^{max}\left(t_0\right)$ are finite. It is also appropriate to introduce the following vector notation over (16):

$$y_0^i=y^i\left(t_0\right)={\left[\begin{array}{cccc}{y}_1^i\left(t_0\right)& y_2^i\left(t_0\right)& \dots & y_n^i\left(t_0\right)\end{array}\right]}^T$$

(19)

$$y^i=y^i\left(t\right)={\left[\begin{array}{cccc}{y}_1^i\left(t\right)& y_2^i\left(t\right)& \dots & y_n^i\left(t\right)\end{array}\right]}^T$$

(20)

In [25], by definition, the secant curve between two points ${}^{k}y_j^i$ and ${}^{k+1}y_j^i$ belonging to the curve $y_j^i\left(t\right)$ to $j=1,2,\dots ,n$ is the line segment joining these two points. Thus, the tangents of the secants between the points ${}^{k}y_1^i$ and ${}^{k+1}y_1^i$, ${}^{k}y_2^i$ and ${}^{k+1}y_2^i$ , …, ${}^{k}y_n^i$ and ${}^{k+1}y_n^i$ are defined as:

$$ta{n}_{\Delta t}{\alpha }^i(t+k·\Delta t)=ta{n}_{\Delta t}{}^k\alpha ^i=$$

$${\left[\begin{array}{cccc}ta{n}_{\Delta t}{}^k\alpha _1^i& ta{n}_{\Delta t}{}^k\alpha _2^i& \dots & ta{n}_{\Delta t}{}^k\alpha _n^i\end{array}\right]}^T$$

(21)

with,

$$ta{n}_{\Delta t}{}^k\alpha _j^i={\displaystyle \frac{{}^{k+1}y_j^i-{}^{k}y_j^i}{\Delta t}}$$

(22)

for $j=1,2,\dots ,n$. Thus, we can state now the two fundamental theorems to consolidate our proof. These two theorems are the differential mean value theorem and the integral mean value theorem, which are presented below without proof [25,26].

Theorem 1 (The Differential Mean Value Theorem): If a function $g_j^i\left(t\right)$ for $j=1,2,\dots ,n$ is a continuous function and defined over the closed interval $[t_k,t_{k+1}]$ is differentiable over the open interval $(t_k,t_{k+1})$, then there is at least one number $t_k^{*}$ with $t_k<t_k^{*}<t_k+\Delta t=t_{k+1}$ such that,

$${\dot{g}}_j^i\left(t_k^{*}\right)={\displaystyle \frac{{}^{k+1}g_j^i-{}^{k}g_j^i}{\Delta t}}$$

(23)

Theorem 2 (The Integral Mean Value Theorem): If a function $g_j^i\left(t\right)$ for $j=1,2,\dots ,n$ is a continuous function and defined over the closed interval $[t_k,t_{k+1}]$, then there is at least one inner point $t_k^x$ in $[t_k,t_{k+1}]$ such that,

$$g_j^i\left(t_k^x\right)·\Delta t={\int }_{t_k}^{t_{k+1}}{g}_j^i\left(t\right)dt$$

(24)

It is important to note that generally $t_k^{*}$ is different from $t_k^x$. Also, the mean value theorems say nothing about how to determine the value of $t_k^{*}$ and $t_k^x$. These two theorems simply state that $t_k^{*}$ and $t_k^x$ are confined to the closed interval $[t_k,t_{k+1}]$.

Property 1: Applying Theorem 2 on the curve ${\dot{g}}_j^i\left(t\right)$ is equivalent to applying Theorem 1 on the curve $g_j^i\left(t\right)$ both on the same interval closed $[t_k,t_{k+1}]$, that is, ${\dot{g}}_j^i\left(t_k^x\right)={\dot{g}}_j^i\left(t_k^{*}\right)$.

Proof: from Theorem 2 applied to the continuous curve ${\dot{g}}_j^i\left(t\right)$ results in a ${\dot{g}}_j^i\left(t_k^x\right)$ for $t_k<t_k^x<t_{k+1}$ such that ${\dot{g}}_j^i\left(t_k^x\right)·\Delta t=$${\int }_{t_k}^{t_{k+1}}{\dot{g}}_j^i\left(t\right)dt=$${}^{k+1}g_j^i-{}^{k}g_j^i$ as a result of the fundamental theorem of calculus. Thus, ${\dot{g}}_j^i\left(t_k^x\right)=$${\displaystyle \frac{{}^{k+1}g_j^i-{}^{k}g_j^i}{\Delta t}}\stackrel{\mathrm{def}}{=}$$ta{n}_{\Delta t}{}^k\alpha _j^i$. On the other hand, the application of Theorem 1 about the continuous and differentiable curve $g_j^i\left(t\right)$ implies the existence of a ${\dot{g}}_j^i\left(t_k^{*}\right)$ for $t_k<t_k^x<t_{k+1}$, such that ${\dot{g}}_j^i\left(t_k^{*}\right)=$${\displaystyle \frac{{}^{k+1}g_j^i-{}^{k}g_j^i}{\Delta t}}\stackrel{\mathrm{def}}{=}$$ta{n}_{\Delta t}{}^k\alpha _j^i$. Thus, ${\dot{g}}_j^i\left(t_k^x\right)=$${\dot{g}}_j^i\left(t_k^{*}\right)$. However, this proof does not prove that $t_k^x=$$t_k^{*}$. □

That said, we prove now the main theorem of this section, namely Theorem 3. But before that, it is worth commenting on some preliminary considerations. In Theorem 3 when we write ${\dot{y}}^i=f\left(y^i\right)$ we are specifying that the general dynamical system $\dot{y}=f\left(y\right)$ is harnessed to the particular solution of the curve i belonging to the family of curves that is a solution of $\dot{y}=f\left(y\right)$. Thus, $\dot{y^i}=f\left(y^i\right)$ is an alternative way of saying that the general dynamical system $\dot{y}=f\left(y\right)$ is harnessed to some initial condition $y\left(t_0\right)={\eta }_0$. Thus, the notation with the index i is more didactic to accept the uniqueness of the mapping explained at the beginning of this sub-section in Figure 3. However, in this article, we did not prove the uniqueness of the mapping discussed in Figure 3, but we assume it is true. To better understand this type of mapping, see a practical example of the same that appears in [19] (see Fig. 7 in this same reference).

Another essential consideration for understanding the proof of Theorem 3 is that it is necessary to know the solution of the dynamic system beforehand, which can then be solved using the Euler integrator that works with mean derivative functions. Also, notice that, since we are working with supervised learning with input/output training patterns, this is not absurd. This fact is easy to understand, as it is necessary to know the references ${}^{k}y_j^i$ and ${}^{k+1}y_j^i$ to accurately estimate the mean derivative $ta{n}_{\Delta t}{}^k\alpha _j^i$, between these two points, through a universal approximator of functions with exclusively feed-forward architecture.

Thus, if the solution of the non-linear dynamical system that is wanted to estimate, then it is necessary to treat this problem as n independent equations with only one time variable t each (black-box approach), rather than n vector coupled equations (white-box approach). This last statement can be justified with the help of Principle 1, which is proved, by reduction to absurdity, right after the proof of Theorem 3.

Therefore, if the solution of the considered dynamic system is previously known and given by $y_j^i\left(t\right)=g_j^i\left(t\right)$ for $j=1,2,\dots ,n$ and $i=1,2,\dots ,\infty$ then the function ${\dot{y}}_j^i\left(t\right)=$${\dot{g}}_j^i\left(t\right)$ are easily obtained by direct differentiation (numerically or analytically). In this way, the original dynamic system ${\dot{y}}_j^i=$$f_j(y_1^i,y_2^i,\dots ,y_n^i)$$=f_j\left(y^i\right)$ can be replaced by ${\dot{y}}_j^i\left(t\right)=$${\dot{g}}_j^i\left(t\right)$ for $j=1.2,\dots ,n$ and $i=1.2,\dots ,\infty$, for a time interval confined in $[t_o,t_f]$ with the help of Principle 1.

However, it should be noted that the math function ${\dot{g}}_j^i\left(t\right)$, in general, has nothing similar to the math function $f_j\left(y^i\right)$. Note also that this simplification makes it possible to reduce the proposed proof to a differential calculus of only one variable. Thus, it allows the use of the differential and integral mean value theorems. Note that ${\dot{g}}_j^i\left(t\right)$ must also have the over-index i, as this function must necessarily also depend on the initial condition as well as its equivalent instantaneous derivative function $f_j\left(y^i\right)$.

Theorem 3: The discrete and exact general solution for the autonomous system of non-linear ordinary differential equations of the type ${\dot{y}}^i=f\left(y^i\right)$ can be established through the first-order Euler relation of the type ${}^{k+1}y_j^i=$$ta{n}_{\Delta t}{}^k\alpha _j^i·\Delta t$$+{}^{k}y_j^i$, for ${}^{k}y_j^i$ and $\Delta t$ fixed; since that the general solution of this dynamical system, given by, $y_j^i\left(t\right)=g_j^i\left(t\right)$ and ${\dot{y}}_j^i\left(t\right)={\dot{g}}_j^i\left(t\right)$ are, previously, known for $j=1,2,\dots ,n$; $i=1,2,\dots ,\infty$ and t$ϵ$$[t_o,t_f]$. Furthermore, the solutions $g_j^i\left(t\right)$ for $j=1,2,\dots ,n$ must all be continuous and differentiable. However, note that ${\dot{g}}_j^i\left(t\right)$ for $j=1,2,\dots ,n$ is suffice to be continuous.

Proof: Let the autonomous non-linear dynamical system of first-order be given by ${\dot{y}}_j^i\left(t\right)=$$f_j(y_1^i,y_2^i,\dots ,y_n^i)$$=f_j\left(y^i\right)$ for $j=1,2,\dots ,n$ and $i=1,2,\dots ,\infty$. If the dynamical system solution is known and equals to $y_j^i\left(t\right)=$$g_j^i\left(t\right)$ then the function $f_j\left(y^i\right)$ can be replaced by ${\dot{g}}_j^i\left(t\right)=$$h_j^i\left(t\right)$, that is, ${\dot{y}}_j^i\left(t\right)=h_j^i\left(t\right)$. In this way, we can write that ${\int }_{{}^{k}y_j^i}^{{}^{k+1}y_j^i}d{y}_j^i\left(t\right)=$${\int }_{t_k}^{t_{k+1}}{h}_j^i\left(t\right)dt$. Note that the indices i, j, and k uniquely map the secant that interests us according to Figure 3. This last integral can still be simplified as ${}^{k+1}y_j^i-{}^{k}y_j^i=$${\int }_{t_k}^{t_{k+1}}{h}_j^i\left(t\right)dt$ as a result of the fundamental theorem of calculus. So, applying the integral mean value theorem, in this last expression, we get ${}^{k+1}y_j^i-{}^{k}y_i^j=$$h_j^i\left(t_k^x\right)·\Delta t=$${\dot{g}}_j^i\left(t_k^x\right)·\Delta t$ where $t_k^x$$ϵ$$[t_k,t_{k+1}]$ and for $h_j^i\left(t\right)=$${\dot{g}}_j^i\left(t\right)$. On the other hand, applying the differential mean value theorem on the function $y_j^i\left(t\right)=$$g_j^i\left(t\right)$ for $j=1,2,\dots ,n$ in the interval $[t_k,t_{k+1}]$ we get that ${\dot{g}}_j^i\left(t_k^{*}\right)=$${\displaystyle \frac{{}^{k+1}g_j^i-{}^{k}g_j^i}{\Delta t}}$$\stackrel{\mathrm{def}}{=}ta{n}_{\Delta t}{}^k\alpha _j^i$. So, by Property 1 we have that ${\dot{g}}_j^i\left(t_k^{*}\right)=$${\dot{g}}_j^i\left(t_k^x\right)=$$ta{n}_{\Delta t}{}^k\alpha _j^i$ to $j=1,2,\dots ,n$. So, ${}^{k+1}y_j^i=$$ta{n}_{\Delta t}{}^k\alpha _i^j·\Delta t$$+{}^{k}y_j^i$ to $j=1,2,\dots ,n$ or, in vector form, it turns out that ${}^{k+1}y^i=$$ta{n}_{\Delta t}{}^k\alpha ^i$$·\Delta t$$+{}^{k}y^i$${.}_{\square }$

The reason that the functions $g_j^i\left(t\right)$ for $j=1,2,\dots ,n$ are all continuous and differentiable is that, in the proof of Theorem 3, we used the differential mean value theorem over $g_j^i\left(t\right)$ and this theorem requires these two conditions to be applied. The reason that the functions ${\dot{g}}_j^i\left(t\right)$ for $j=1,2,\dots ,n$ are all continuous is that, in the proof of the same Theorem 3, the integral mean value theorem was applied over the function ${\dot{g}}_j^i\left(t\right)$ and this theorem requires this condition. Note that for the application of the integral mean value theorem, the function ${\dot{g}}_j^i$ is not required to be differentiable.

Principle 1: Given the non-linear autonomous dynamic system ${\dot{y}}_j^i=$$f_j(y_1^i,y_2^i,\dots ,y_n^i)$$=f_j\left(y^i\right)$ and their respective general solutions $y_j^i\left(t\right)=g_j^i\left(t\right)$ for $j=1,2,\dots ,n$ and $i=1,2,\dots ,\infty$ then, the autonomous instantaneous derivative functions $f_j\left(y^i\right)$ can be replaced by the exclusively non-autonomous instantaneous derivative functions given by ${\dot{g}}_j^i\left(t\right)$, that is, $f_j(y_1^i,y_2^i,\dots ,y_n^i)=$$f_j\left(y^i\right)=$${\dot{g}}_j^i\left(t\right)=$$h_j^i\left(t\right)$ for $j=1,2,\dots ,n$ and $i=1,2,\dots ,\infty$.

Proof: Assume, for absurdity, that the referred principle is false. If this is true, then there would be no compromise of veracity between the computational data acquisition system (through sensors) and the real universe, and thus the black-box approach would not work; but this is absurd. So, by exclusion, the principle is true . □

Notice that when a computational data acquisition system is performed, for example, from a particular real-world plant, something interesting happens. Thus, whoever maintains the consistency of the acquired data, through several sensors acting simultaneously and independently of each other, is the omnipresent and immutable manifestation of the natural physical laws, which govern the proper functioning of the universe (e.g., the law of universal gravitation, conservation of energy law, Faraday’s law, among others). In this case, little is allowed to human wisdom and technology.

Finally, the E-TUNI universal numerical integrator is suitable for solving only ordinary differential equations. However, it should be noted that the current use of neural networks in solving Partial Differential Equations (PDEs) is quite extensive [27,28,29,30,31,32,33,34,35]. Additionally, [36] explains how to use the Runge-Kutta numerical integrator to solve partial differential equations (mainly hyperbolic ones). Therefore, using the E-TUNI to solve partial differential equations may be a good option for future work and should be investigated more carefully.

3.3. Mathematical Relationship Between Mean and Instantaneous Derivatives

A very remarkable fact is that it is possible to obtain an equation that relates the instantaneous derivative functions to the mean derivative functions. This equation can be easily obtained using the chain rule for functions on several independent variables [25,26]. In this way, two previously trained neural networks are known to represent the instantaneous and mean derivative functions, respectively, by:

$${\left[\begin{array}{cccc}{\dot{y}}_1\left(t\right),& {\dot{y}}_2\left(t\right),& \dots ,& {\dot{y}}_{n_y}\left(t\right)\end{array}\right]}^T=$$

(25)

$${\left[\begin{array}{cccc}tan{}^k\theta _1^i,& tan{}^k\theta _2^i,& \dots ,& tan{}^k\theta _{n_y}^i\end{array}\right]}^T=$$

$$f_{NN}^{id}{\left[\begin{array}{ccccc}{y}_1\left(t\right),& y_2\left(t\right),& \dots ,& y_{n_y}\left(t\right),& \widehat{w}\end{array}\right]}^T$$

and

$${\left[\begin{array}{cccc}ta{n}_{\Delta t}{}^k\alpha _1^i,& ta{n}_{\Delta t}{}^k\alpha _2^i,& \dots ,& ta{n}_{\Delta t}{}^k\alpha _{n_y}^i\end{array}\right]}^T=$$

(26)

$$f_{NN}^{md}{\left[\begin{array}{ccccc}{y}_1\left(t\right),& y_2\left(t\right),& \dots ,& y_{n_y}\left(t\right),& \widehat{w}\end{array}\right]}^T$$

For reasons of simplification, we do not consider the case using control variables in the equations (25) and (26). So, if we use the chain rule [25,26] in the equation (26) we get:

$${\displaystyle \frac{d}{dt}}ta{n}_{\Delta t}{}^k\alpha _j^i=$$

(27)

$${\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _j^i}{\partial {}^{k}y_1^i}}·{\displaystyle \frac{d{}^{k}y_1^i}{dt}}+\dots +{\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _j^i}{\partial {}^{k}y_{n_y}^i}}·{\displaystyle \frac{d{}^{k}y_{n_y}^i}{dt}}=$$

$$\left[\begin{array}{ccc}{\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _j^i}{\partial {}^{k}y_1^i}},& \dots ,& {\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _j^i}{\partial {}^{k}y_{n_y}^i}}\end{array}\right]·\left[\begin{array}{c}{\displaystyle \frac{d{}^{k}y_1^i}{dt}}\\ ⋮\\ {\displaystyle \frac{d{}^{k}y_{n_y}^i}{dt}}\end{array}\right]=$$

$$\left[\begin{array}{ccc}{\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _j^i}{\partial {}^{k}y_1^i}},& \dots ,& {\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _j^i}{\partial {}^{k}y_{n_y}^i}}\end{array}\right]·\left[\begin{array}{c}tan{}^k\theta _1^i\\ ⋮\\ tan{}^k\theta _{n_y}^i\end{array}\right]=$$

$${\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _j^i}{\partial {}^{k}y^i}}\circ tan{}^k\theta ^{i}toj= 1, 2, \dots , n_y$$

For the sake of simplicity, we can represent the expression (27) as follows:

$$ta{n}_{\Delta t}^k{\Psi }_j^i={\displaystyle \frac{d}{dt}} ta{n}_{\Delta t}^k{\alpha }_j^i={\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha _j^i}{\partial {}^{k}y^i}}\circ tan{}^k\theta ^i$$

(28)

$$toj=1, 2, \dots , n_y$$

The operator ∘ is just the usual dot product applied over a Euclidean space. So, just for completeness, the vector form of the equation (28) can also be expressed by:

$${\displaystyle \frac{d}{dt}} ta{n}_{\Delta t}^k{\alpha }^i={\displaystyle \frac{\partial ta{n}_{\Delta t}{}^k\alpha ^i}{\partial {}^{k}y^i}}·tan{}^k\theta ^i$$

(29)

where,

$${\displaystyle \frac{d}{dt}} ta{n}_{\Delta t}^k{\alpha }^i=\left\{\begin{array}{c}{\displaystyle \frac{d}{dt}} ta{n}_{\Delta t}^k{\alpha }_1^i\\ {\displaystyle \frac{d}{dt}} ta{n}_{\Delta t}^k{\alpha }_2^i\\ ⋮\\ {\displaystyle \frac{d}{dt}} ta{n}_{\Delta t}^k{\alpha }_{n_y}^i\end{array}\right\}$$

(30)

$${\displaystyle \frac{\partial ta{n}_{\Delta t}^k{\alpha }^i}{\partial {}^{k}y^i}}=\left[\begin{array}{ccc}{\displaystyle \frac{\partial ta{n}_{\Delta t}^k{\alpha }_1^i}{\partial {}^{k}y_1^i}}& \dots & {\displaystyle \frac{\partial ta{n}_{\Delta t}^k{\alpha }_1^i}{\partial {}^{k}y_{n_y}^i}}\\ {\displaystyle \frac{\partial ta{n}_{\Delta t}^k{\alpha }_2^i}{\partial {}^{k}y_1^i}}& \dots & {\displaystyle \frac{\partial ta{n}_{\Delta t}^k{\alpha }_2^i}{\partial {}^{k}y_{n_y}^i}}\\ ⋮& \ddots & ⋮\\ {\displaystyle \frac{\partial ta{n}_{\Delta t}^k{\alpha }_{n_y}^i}{\partial {}^{k}y_1^i}}& \dots & {\displaystyle \frac{\partial ta{n}_{\Delta t}^k{\alpha }_{n_y}^i}{\partial {}^{k}y_{n_y}^i}}\end{array}\right]$$

(31)

$$tan{}^k\theta ^i=\left\{\begin{array}{c}tan{}^k\theta _1^i\\ tan{}^k\theta _2^i\\ ⋮\\ tan{}^k\theta _{n_y}^i\end{array}\right\}$$

(32)

There is an essential fact between instantaneous and mean derivative functions. Using these two kinds of derivative functions, that is, $tan{}^k\theta ^i$ and $ta{n}_{\Delta t}{}^k\alpha ^i$ it is possible to interpolate a parabola between the discrete interval $[t_k,t_{k+1}]$, which passes at least at two precise points of the exact solution of the considered non-linear dynamical system, in this same interval. Thus, Table 1 considers the points used to perform this parabolic interpolation, and the equations (34), (35), and (36) are the coefficients of the parabola (33).

$${}^{k}y_j^i\left(t\right)={\alpha }_k·t^2+{\beta }_k·t+{\gamma }_k$$

(33)

$${\alpha }_k={\displaystyle \frac{1}{2}}·ta{n}_{\Delta t}{}^k\Psi _j^i$$

(34)

$${\beta }_k=ta{n}_{\Delta t}{}^k\alpha _j^i-{\displaystyle \frac{1}{2}}·(t_{k+1}+t_k)·ta{n}_{\Delta t}{}^k\Psi _j^i$$

(35)

$${\gamma }_k={\displaystyle \frac{1}{2}}·(t_k·t_{k+1}·ta{n}_{\Delta t}{}^k\Psi _j^i-2·t_k·ta{n}_{\Delta t}{}^k\alpha _j^i+2{·}^k{y}_j^i)$$

(36)

where, $t ϵ [t_k,t_{k+1}].$

4. Results and Analysis

In this article, we perform a complete computational numerical analysis, studying the training of the neural integration structures discussed in the previous sections, i.e., the E-TUNI. Thus, for all the experiments that are presented below, it was standardized to use MLP neural networks with the traditional back-propagation [4] and Kalman filter [37,38,39] training algorithms in Example 2; and the Levenberg–Marquardt [40] training algorithm in Example 3. All experiments were trained using only one inner layer in the MLP network.

Example 1. Check the validity of Principle 1 for the dynamical system $\dot{y}=f\left(y\right)=y^2$.

Proof: from $\dot{y}={\displaystyle \frac{dy}{dt}}=y^2$ then, ${\int }_{y\left(t_0\right)}^{y\left(t\right)}{y}^{-2}dy={\int }_{t_o}^{t}dt$, as a result of the fundamental theorem of calculus. Solving this integral, analytically, results in $y\left(t\right)={\displaystyle \frac{y\left(t_0\right)}{-y\left(t_0\right)t+y\left(t_0\right)t_0+1}}$. Differentiating the function $y\left(t\right)$, with respect to time t, we have that $\dot{y}=\dot{g}\left(t\right)$ = $y\left(t_0\right){\displaystyle \frac{d}{dt}}{[-y\left(t_0\right)t+y\left(t_0\right)t_0+1]}^{-1}$ = ${\displaystyle \frac{y^2\left(t_0\right)}{{[-y\left(t_0\right)t+y\left(t_0\right)t_0+1]}^2}}$. Thus, replacing $y\left(t\right)$ in $\dot{y}=y^2$ results in $\dot{y}=f\left(y\right)=y^2=$${\left[\begin{array}{c}{\displaystyle \frac{y\left(t_0\right)}{-y\left(t_0\right)t+y\left(t_0\right)t_0+1}}\end{array}\right]}^2$$=\dot{g}\left(t\right)$.

The graph in Figure 4 illustrates the drawing of the function $y\left(t\right)$ (upper part) and the drawing of the function $\dot{y}\left(t\right)$ (lower part). For a proper understanding of the graphs in Figure 4, it should be noted that $t_0=0$ and there are seven curves in the graphs of Figure 4 representing the solution family of Example 1. For example, in $y\left(t_0\right)=y\left(0\right)=1$ part a particular solution of the dynamical system $\dot{y}=y^2$ that have a vertical asymptote at $t=1$ (upper part of Figure 4). Note that to the left of this asymptote, the solution of the dynamical system has been plotted in blue. On the other hand, to the right of this asymptote, the solution has been plotted in red. This convention is followed to plot the remaining solutions that appear in Figure 4. It is also important to note that the blue curves have distinct vertical asymptotes for each of the initial conditions. On the other hand, all red curves have horizontal asymptotes at $y=0$, when t tends to infinity.

Note that Principle 1 is valid for discontinuous solutions in $y\left(t\right)$. However, the use of the general expression of E-TUNI to learn this discontinuous solution is not possible according to the theory presented in this article. The reasons for this impossibility are as follows: (i) the integral and differential mean value theorems require that the E-TUNI solution be continuous and differentiable, and (ii) E-TUNI would have to learn a numerical value equal to infinity, at the points of the vertical asymptotes shown in Figure 4 (upper part), and this is impossible, as conventional computers have finite memory. For the use of E-TUNI, in examples like this, further studies are required.

Example 2. Solve the predictive control problem, necessarily using E-TUNI as a model of the non-linear plant, for the system of non-linear ordinary differential equations, which describes the Earth/Mars orbit transfer dynamics for a rocket of mass m, as shown in Figure 5. Solve this problem for two different horizons, that is, $m=1$ and $m=15$, according to equations (10) to (15). The set of ordinary differential equations describing this problem is given below [30]:

$$\begin{array}{c}\dot{m}=-0.0749\\ \dot{r}=w\\ \dot{w}={\displaystyle \frac{v^2}{r}}-{\displaystyle \frac{\mu }{r^2}}+{\displaystyle \frac{T·sin\theta }{m}}\\ \dot{v}={\displaystyle \frac{-w·v}{r}}+{\displaystyle \frac{T·cos\theta }{m}}\end{array}$$

(37)

An illustrative scheme of this dynamic system can be seen in Figure 5. So, the state variables for this predictive control problem are m (rocket mass), r (orbit radius), w (radial velocity), and v (transverse velocity). The only controlling variable is the thrust angle $\theta$ of the rocket. The normalized constants of this problem are $\mu =1.0$ (gravitational constant), $T=0.1405$ (rocket thrust), $t_o=0$ (initial instant), and $t_f=3.3$ (final instant). In this problem, each unit of time equals 58.2 days.

We trained the E-TUNI using an MLP network with only one inner layer. This inner layer contained 41 neurons and one bias neuron. We applied the hyperbolic tangent activation function in the inner layer, while we used a purely linear activation function at the network output. The training and testing yielded mean square errors (MSE) of $9,021x{10}^{-6}$ and $1,067x{10}^{-5}$, respectively.

In the neural training of the E-TUNI structure, two standard training algorithms were used: the classic back-propagation [4] and the Kalman filter with recursive and parallel processing [39]. They have been casually changed. In addition, this training took a considerable amount of time, as a shallow architecture network was used instead of a deep architecture network.

Figure 6 presents the trajectory traced by the Kalman filter, for a predictive control structure whose plant model was obtained through E-TUNI. The solid line is the reference trajectory, and the dotted line is the trajectory tracked by the Kalman filter. Figure 7 presents the control estimated by the Kalman filter in a predictive control structure designed according to equations (10) to (15) and associated with the trajectory traced in Figure 6. However, the Kalman filter equations employed were taken from [16]. The integration step used was $\Delta t=0.1$ and the horizon equal to $m=1$. As can be seen, the estimated control was not very good for just using one horizon in equation (10).

On the other hand, Figure 8 and Figure 9 present the same E-TUNI structure as in the previous case, but now in a new control estimating with the horizon used equally to $m=20$. Note that this small change significantly improved the result obtained. However, while the training time for a horizon of $m=1$ may only take between an hour or two, the training time for a horizon of $m=20$ can be as long as ten days, if the integration step is changed from $\Delta t=0.1$ to $\Delta t=0.01$. For a better understanding of E-TUNI coupled to a predictive control structure see reference [16].

Figure 6, Figure 7, Figure 8 and Figure 9 are an original extension of what can be found in [16]. In [16], the predictive control structure with E-TUNI for reference trajectories displaced from the initial condition was not explored, as presented here. Similar behavior for the other state variables m, w, and v was observed. Therefore, we have omitted the other graphics due to space constraints.

Figure 10 and Figure 11 were taken from [16]. They represent the application of a predictive control structure on the E-TUNI, for the case where a deviation between the reference trajectories and the initial condition of the dynamic system in question was not imposed. However, note that an integration step of $\Delta t=0.01$ and a horizon of $m=10$ was used. Note also that the 10-fold reduction in the integration step, although it makes the algorithm much slower, manages to considerably refine the control policy obtained in Figure 11, compared to Figure 7 and Figure 9.

Example 3. The non-linear simple pendulum is a second-order autonomous system given by $l·\ddot{\theta }+g·sin\theta =0$. Note that l and g are, respectively, the length of the pendulum and the local acceleration due to gravity. Here, we solve this problem by using E-TUNI for several different integration steps. Furthermore, for each of these integration steps, we obtain the interpolating parabolas governed by the equations from (33) to (36).

To precisely understand the resolution of Example 3, we follow the following basic algorithm:

Step 1. Generate the E-TUNI input/output training patterns through the Runge-Kutta 4-5 integrator, applied to the non-linear simple pendulum equation. Note that, as this dynamical system is a non-linear equation, its solution in the phase plane is not a perfect circle.

Step 2. Use the Levenberg-Marquardt algorithm in the direct approach to training two distinct neural networks, namely: (i) a neural network to learn the instantaneous derivative functions and (ii) another neural network to learn the mean derivative functions.

Step 3. Determine the vector ${\displaystyle \frac{d}{dt}}ta{n}_{\Delta t}{}^k\Psi ^i$ using the equations from (29) to (32). Just pay attention to the fact that the partial derivatives ${\displaystyle \frac{\partial }{{\partial }^k{y}^i}}ta{n}_{\Delta t}{}^k\alpha ^i$ were obtained numerically and not analytically. For this, a ${\Delta }^k{y}^i={10}^{-6}$ was used.

Step 4. Determine the interpolating parabolas within a horizon of interest $[t_0,t_f]$ and with a step of $\Delta t$. For this, the equations from (33) to (36) were used recursively.

Step 5. Compare the results obtained for different $\Delta t$ integration steps. Also, compare the solution of the interpolating parabolas with the mean derivative equations.

Step 6. A summary of this algorithm is presented in Figure 12 and Figure 13.

To solve this problem appropriately, all graphical simulation results will be displayed using an MLP neural network with only one inner layer (shallow network). In the inner layer, the hyperbolic tangent sigmoid activation function was used, and in the output of the neural network, a purely linear function was used. The training algorithm used was the traditional Levenberg-Marquardt [40], always using the direct approach of E-TUNI and not the indirect approach. For that, the MATLAB neural networks toolbox was used.

The dynamical system in question was trained within the finite domains $Do{m}^{\theta }=[-\pi /2,+\pi /2]$$rad$ and $Do{m}^{\dot{\theta }}=[-2,+2]$${\displaystyle \frac{rad}{s}}$. For all neural networks trained, 400 training patterns were used. $75\%$ of them were used for training, $15\%$ for validation, and $15\%$ for testing. All training has been standardized to be trained with exactly $10,000$ epochs. Table 2 summarizes the main results achieved in these training sessions.

Analyzing Table 2 it is observed that the greater the integration step $\Delta t$, used in the E-TUNI training, the greater the training error for the validation patterns. This fact experimentally confirms (9), which states that the integration step $\Delta t$ amplifies the Mean Squared Error (MSE) of E-TUNI training when it is greater than $1.00$.

However, the most interesting numerical results, concerning Example 3, are shown in Figure 12 and Figure 13. In Figure 12, the integration steps are decreased from top to bottom. Notice that when the integration step $\Delta t$ is equal to $10.0$, E-TUNI could not accurately learn the angular position prediction $\theta$ (see all blue points in Figure 12). With the integration steps $\Delta t$ equal to or less than $1.00$, E-TUNI was able to learn them perfectly. The reason for this is that several E-TUNI training was carried out for the integration step $\Delta t$ equal to $10.0$, but the best result achieved was the one presented in Table 2, that is, $mse=1.7520\times {10}^{-3}$, confirming, again, the validity of the equation (9), as a limiting factor for this method. The E-TUNI training, which involved integration steps $\Delta t$ equal to or less than $1.00$, was easily learned by the MLP neural networks in the first attempts.

Figure 13 also presents several simulation graphs with the integration steps $\Delta t$, used in the E-TUNI training, also decreasing from top to bottom. The dotted pink lines represent the true values of the dynamical system in question. The solid lines in blue represent the mean derivatives between two consecutive points. Continuous lines in red, interpolating parabolas between two successive points. Note that, for tiny $\Delta t$ integration steps, the parabolic interpolation using mean and instantaneous derivatives functions, simultaneously, is quite accurate and relatively better than using E-TUNI alone.

Analyzing Figure 13 more closely, observe that the E-TUNI equation is equivalent to consecutive uniform rectilinear movements and with constant velocities equal to the mean derivatives (constant from interval to interval). On the other hand, the equations of the interpolating parabolas are consecutive equations of uniformly varied motions and with consecutive accelerations ${tan}_{\Delta t}{}^k\Psi _j^i$ (also constant interval to interval).

Note also that uniformly varied motions (interpolating parabolas) allow at most one change of direction in their motion. However, uniform rectilinear motion (E-TUNI) does not. For this reason, it appears that interpolating parabolas are more accurate than mean derivative curves for tiny integration steps. Because neural networks are universal approximators of functions [2,5,6], the physical variables that can be used as input to E-TUNI can be any.

Additionally, the Euler integrator, designed exclusively with instantaneous derivative functions, is never able to achieve the equivalent precision of the Euler integrator, designed with mean derivative functions. This fact happens mainly if substantial integration steps are used.

Finally, it is only possible to use control variables (in the inputs of E-TUNI) to obtain the parabolic interpolations if the control variable remains with its constant value throughout the entire interval $[t_k,t_{k+1}]$. This property is only practical for control if the E-TUNI integration step is tiny and close to zero.

5. Conclusions

In this paper, a detailed mathematical development of the Euler-Type Universal Numerical Integrator (E-TUNI) was presented. Its main mathematical properties were presented, with a direct and quite objective development. A correct demonstration of the general expression of E-TUNI was also presented, since in previous articles, this demonstration was presented in the wrong way. Three practical case studies using E-TUNI were presented, demonstrating the computational versatility of E-TUNI.

The E-TUNI is an alternative methodology to the traditional NARMAX methodology. The E-TUNI can also be an alternative to the conventional Runge-Kutta Neural Network (RKNN) and Adams-Bashforth Neural Network (ABNN). However, it is essential to note that the NARMAX and mean derivatives (E-TUNI) methodologies are fixed-step, whereas the instantaneous derivatives methodology (RKNN, ABNN, among others) is variable-step, when both are used in the simulation phase. In the training phase, all methods are fixed-step.

To understand how the Euler-Type Universal Numerical Integrator (E-TUNI) works, it is helpful to use Figure 1. This figure illustrates the elementary difference between mean derivative functions and instantaneous derivative functions. As everyone knows, the secant between two points on a curve is the mean derivative, and the tangent line to a point on a curve is the instantaneous derivative function.

Thus, if a Runge-Kutta Neural Network (RKNN) of order four is trained with a tiny training error, then in [20] it is demonstrated that the neural network converges to the instantaneous derivative functions. However, if now this neural network is decoupled from the Runge-Kutta integrator and re-coupled into an E-TUNI, then the output error in the latter integrator will increase because, as illustrated in Figure 1, an Euler integrator designed with an instantaneous derivative function generates a substantial error.

However, if neural training is restarted in E-TUNI, as the neural network is a universal function approximator, then it reduces the training error back to zero. As the solution of autonomous and non-linear dynamical systems is unique, there is only one way to correct this error. So, as Figure 1 illustrates, it is done by making the instantaneous derivative functions converge to the mean derivative functions in the neural architecture. However, as the secant between two points varies with the variation of the horizontal distance between these two points, the E-TUNI works with a fixed integration step. If it is desired to change this step, then new neural training have to be done. Note that the high-order Runge-Kutta Neural Network (RKNN) can vary the integration step within certain limits since the instantaneous derivative function does not depend on the integration step.

Principle 1 can be applied to both continuous and discontinuous solutions of real-world dynamical systems governed by autonomous non-linear ordinary differential equations. However, the general E-TUNI expression using mean derivative functions only applies to continuous solutions, as demonstrated earlier in this article. The integral and differential mean value theorems impose this limitation. A demonstration of the general expression of E-TUNI, for discontinuous solutions, requires further studies.

Note that Principle 1 can provide a rather interesting theoretical framework for explaining why the black-box approach works well in practice. Note also that this principle establishes a precise relationship of equivalence between ${\dot{y}}_j^i=f_j\left(y^i\right)$ (metaphysical property of the universe) and ${\dot{y}}_j^i=g_j^i\left(t\right)$ (the physical property of the universe). Prior knowledge of ${\dot{y}}_j^i=f_j\left(y^i\right)$ is metaphysical, as these equations can be mentally abstracted through mathematics and the laws of physics only. Prior knowledge of ${\dot{y}}_j^i=f_j\left(y^i\right)$ does not require carrying out any laboratory experiments or numerical calculations. On the other hand, knowledge of ${\dot{y}}_j^i=g_j^i\left(t\right)$, in general, requires carrying out a laboratory experiment or analytical (or numerical) resolution of ${\dot{y}}_j^i=f_j\left(y^i\right)$.

Finally, the software used to test the E-TUNI was developed using only shallow neural networks. Furthermore, the predictive control part could not be availed from the Artificial Neural Networks Toolbox of MATLAB, as it required the implementation of the original equations from (10) to (15). In this way, the software was relatively slow in terms of processing time. Thus, modern techniques and algorithms involving Deep Neural Networks (DNNs) or Convolutional Neural Networks (CNNs) can significantly improve the training time and simulation time of E-TUNI in the future.