# Polynomial Neural Networks and Taylor maps for Dynamical Systems Simulation and Learning

The connection of Taylor maps and polynomial neural networks (PNN) to solve ordinary differential equations (ODEs) numerically is considered. Having the system of ODEs, it is possible to calculate weights of PNN that simulates the dynamics of these equations. It is shown that proposed PNN architecture can provide better accuracy with less computational time in comparison with traditional numerical solvers. Moreover, neural network derived from the ODEs can be used for simulation of system dynamics with different initial conditions, but without training procedure. On the other hand, if the equations are unknown, the weights of the PNN can be fitted in a data-driven way. In the paper we describe the connection of PNN with differential equations in a theoretical way along with the examples for both dynamics simulation and learning with data.

## Authors

• 8 publications
• 3 publications
• 3 publications
• ### Matrix Lie Maps and Neural Networks for Solving Differential Equations

The coincidence between polynomial neural networks and matrix Lie maps i...
08/16/2019 ∙ by Andrei Ivanov, et al. ∙ 0

• ### Inverse modified differential equations for discovery of dynamics

We introduce inverse modified differential equations (IMDEs) to contribu...
09/02/2020 ∙ by Aiqing Zhu, et al. ∙ 0

• ### Neural Delay Differential Equations

Neural Ordinary Differential Equations (NODEs), a framework of continuou...
02/22/2021 ∙ by Qunxi Zhu, et al. ∙ 16

• ### Unsupervised Reservoir Computing for Solving Ordinary Differential Equations

There is a wave of interest in using unsupervised neural networks for so...
08/25/2021 ∙ by Marios Mattheakis, et al. ∙ 0

• ### Learning Differential Equations that are Easy to Solve

Differential equations parameterized by neural networks become expensive...
07/09/2020 ∙ by Jacob Kelly, et al. ∙ 60

• ### Composition Methods for Dynamical Systems Separable into Three Parts

New families of fourth-order composition methods for the numerical integ...
06/11/2020 ∙ by Fernando Casas, et al. ∙ 0

• ### Accelerating Neural ODEs with Spectral Elements

This paper proposes the use of spectral element methods canuto_spectral_...
06/17/2019 ∙ by Alessio Quaglino, et al. ∙ 2

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction and Related Works

Traditional methods for solving systems of differential equations imply a numerical step-by-step integration of the system. For some problems such as stiff ODEs, this integration procedure leads to time-consuming algorithms because of the limitations on the time step that is used to achieve the necessary accuracy of the solution. From this perspective, neural network as a universal function approximator can be applied for the construction of the solution in a more efficient way. This property induces wide research and consequently large number of papers devoted to the neural networks design for ODEs and PDEs solutions. We go into some of the publications in more detail. Further, we give a short review and comparison of existing approaches as well as highlight the advantages of proposed in this paper methodology.

In [16]

, the method to solve initial and boundary value problems using feed-forward neural networks is proposed. The solution of the differential equation is written as a sum of two parts. The first part satisfies the initial/boundary conditions. The second part corresponds to a neural network output. The same technique is applied for solving Stokes problem in

[7, 10, 27].

In [28]

, the neural network is trained to satisfy the differential operator, initial condition, and boundary conditions for the partial differential equation (PDE). The authors in

[30]

translate a PDE to a stochastic control problem and use deep reinforcement learning for an approximation of derivative of the solution with respect to the space coordinate.

Other approaches rely on the implementation of a traditional step-by-step integrating method in a neural network framework, cf. [2, 29]. In [29], the author proposes such an architecture. After fitting, the neural network (NN) produces an optimal finite difference scheme for a specific system. The back-propagation technique through an ordinary differential equation (ODE) solver is proposed in [8]. The authors construct a certain type of neural network that is analogous to a discretized differential equation. This group of methods requires a traditional numerical method to simulate dynamics.

Polynomial neural networks are also widely presented in the literature, cf. [32, 31, 24]. In [32], the polynomial architecture that approximates differential equations is proposed. The Legendre polynomial is chosen as a basis in [31]. But it should be noted, that in all these paper, the polynomial architectures are used as black box models, and the authors do not indicate its connection to the ODEs. Thus, the advantages of NN can only be partially exploited.

In all the described approaches, NN are trained to consider the initial conditions of the differential equations. This means that a NN has to be trained each time when the initial conditions are changed. The above-described techniques are applicable to the general form of differential equations but are able to provide only a particular solution of the system, that is a strict limitation of applicability.

Several types of Convolutional NN and Recurment NN have been developed in the last three years, cf. [9]. These deep networks corresponds to associated ODE discretization schemes. There is some benefit using these special network architecture, but still no clear explanation of it.

In this paper, we consider nonlinear systems of ODEs with polynomial right-hand side,

 ddtX=F(t,X)=∞∑k=0P1k(t)X[k], (1)

where is an independent variable,

is a state vector, and

means -th Kronecker power of vector . For example, for we have , , after reduction of the same terms.

Such nonlinear systems arise in different fields such as automated control, robotics, mechanical and biological systems, chemical reactions, drug development, molecular dynamics, and so on. Moreover, often it is possible to transform a nonlinear equation (either ODE or PDE) to a polynomial form.

Based on the Taylor mapping technique, it is possible to build a polynomial neural network (PNN) that approximates the general solution of (1). The weights of the PNN should be calculated at once directly from the system of ODEs. By a Taylor map we mean transformation in form of

 X(t1)=W0+W1X0+W2X[2]0+…+WkX[k]0, (2)

where , and matrices are weights. The ransformation (2) is linear in weights and nonlinear with respect to the .

In literature, this map can be referred to as Taylor maps and models [19]

, tensor decomposition

[11], matrix Lie transform [3], exponential machines [21], and others. In fact, the transformation (2) is just a polynomial regression with respect to the components of .

Though many numerical solvers for (1) can be considered as maps, they commonly based on small time steps and weight matrices that depend on . In this paper, by mapping approach we mean transformation (2) with the weight matrices

that are estimated for a large enough time interval

and do not depend on . The greatest advantage of the mapping approach is computational performance. Instead of step-by-step integrating with numerical solvers, one can apply a map (2) that estimates dynamics in a large time interval for different initial conditions at the same accuracy.

In the paper, we briefly consider an algorithm for calculating weight matrices in (2) for an arbitrary time step . Considering map (2

) as a neuron (Fig.

1), it is possible to design a PNN that simulates the dynamics of the given differential equation. While the PNN is connected to differential equations, it is also possible to use this architecture for data-driven identification of physical systems.

The rest of the paper is organised as follows. Sec. 2 describes the general approach for building Taylor maps for the systems of ODEs. Sec. 3 corresponds to the PNNs weights calculating for dynamics simulation of applied physical systems. Sec. 4 is devoted to training PNN. The regularization method along with the examples of data-driven dynamics learning without equation consideration are described.

## 2 Taylor maps for the system of ODEs

The solution of (1) with the initial condition during the time interval in its convergence region can be presented in the power series, cf. [4, 3],

 X(t)=M(t−t0)∘X0=∞∑k=0M1k(t)X[k]0, (3)

Theoretical estimations of accuracy and convergence of the truncated series in solving of ODEs can be found in [6]. In [5], it is shown how to calculate matrices by introducing new matrices . The main idea is replacing (1) by the equation

 ddtMik(t)=k∑j=iPij(t)Mjk(t),1≤i

The last equation does not depend on and should be solved at once with initial condition , where

is the identity matrix. The truncated solution of (

4) for desirable time interval yields Taylor map (2) with .

One of the advantages of the described approach is simulation performance. Indeed, instead of step-by-step integrating of the equation (1) with a small time step every time when the initial condition is changed, one should integrate map at once for the unique initial condition and the whole desirable time interval . Then the same map can be applied for different initial conditions.

## 3 Simulation of dynamical systems

Since the Taylor mapping approach is commonly used in accelerator physics [14, 25], we demonstrate the proposed method with the simplified example of charged particle motion. We also introduce deep PNN for simulation and control of charged particle beam dynamics, and discuss a shallow PNN architecture for simulation of a stiff ODE.

### 3.1 Charged particle dynamics

The particle dynamics in the electromagnetic fields can be described by a system of ODEs that has a complex nonlinear form. For simplicity, let’s consider an approximation [26] of particle motion in cylindrical deflector written in form of (1)

 x′=y,y′=−2x+x2/R, (5)

where is equilibrium radius of particle bending, is deviation from this radius, and is derivative on bending angle.

For an example, let’s consider a deflector with that rotates a reference particle with initial conditions on angle . For simulation, we investigate dynamics of particle with nontrivial initial conditions that lead to particle oscillation.

Traditional approach to solve system (5) is step-by-step integrating with numerical solvers. For this purpose, we use Runge–Kutta method of 4th order with fixed time step. To control the numerical error in this example, we use integrating angle in 30 times smaller than bending angle. Bigger steps introduce nonphysical dissipation in particle motion.

In step-by-step integrating, to track particle inside the deflector, one has to do 30 steps. In contrast to this, mapping technique allows to calculate output state in single step (see Fig. 2).

Let’s fine a 3rd order Taylor map (2) that transform initial particle state at the entrance of the deflector to the resulting state at the end of it,

 (xy)=W1(x0y0)+W2⎛⎜ ⎜⎝x20x0y0y20⎞⎟ ⎟⎠+W3⎛⎜ ⎜ ⎜ ⎜ ⎜⎝x30x20y0x0y20y30⎞⎟ ⎟ ⎟ ⎟ ⎟⎠ (6)

Combining this map with (5), one can write

 (xy)′=W′1(x0y0)+W′2⎛⎜ ⎜⎝x20x0y0y20⎞⎟ ⎟⎠+W′3⎛⎜ ⎜ ⎜ ⎜ ⎜⎝x30x20y0x0y20y30⎞⎟ ⎟ ⎟ ⎟ ⎟⎠,(xy)′=(01−20)(xy)+(0001/R00)⎛⎜⎝x2xyy2⎞⎟⎠. (7)

Grouping the like terms of powers of and in the last system, one can obtain a system of ODEs that do not depends on and represents dynamics of weight matrices

 W′1=f1(W1,W2,W3),W1(0)=I,W′2=f2(W1,W2,W3),W2(0)=0,W′3=f3(W1,W2,W3),W3(0)=0, (8)

where is functions arising after like terms grouping. By integrating this system during the interval , we can receive the desired map. For example, the map up to two digits is

Using this polynomial transformation, one can calculate state of the particle at the end of the deflector for arbitrary initial conditions . The simulation in this case is 160 times faster than Runge–Kutta based simulation.

For long-term charged particle dynamics investigation, the traditional step-by-step numerical methods are not suitable because of the performance limitation. Instead of solving differential equations directly, one can estimate Taylor maps of high orders of nonlinearity for each control element in an accelerator. By combining such maps consequently, one can obtain a deep polynomial neural network that represents the whole accelerator lattice (see Fig. 3). In such PNN each layer is represented by a Taylor map and corresponds to the physical element in charge particle accelerator.

The described approach was used for simulation of nonlinear spin-orbit dynamics in EDM (electric dipole moment) search project

[25, 14]. The dynamics of particles was described by nine-dimensional state vector. The deep PNN architecture had more then 100 polynomial layers that represent lattice of the accelerator. The computational performance was increased in 1500 times in comparison with traditional step-by-step integrating.

Along with the computational performance, another advantage of the proposed model is uncertainties accounting. It is easy to represent misalignments and field errors in physical accelerator by introducing additional polynomial layers inside the PNN. This feature may be essential for control of accelerators, while the PNN architecture provides a possibility to fit layers with measured data.

### 3.2 The Rayleight-Plesset equation

The Rayleight-Plesset equation governs the dynamics of a spherical gas bubble in an infinite body of in-compressible liquid. It is derived using the conservation of mass and momentum which leads to a highly non-linear second order ODE for the bubble radius . The pressure inside the bubble is denoted by and the pressure outside far a way from the bubble is denoted by , is the density of the surrounding liquid, assumed to be constant. The equation of motion has the form

 Rd2Rdt2+32R(dRdt)2=1ρ(pB−p). (9)

Viscous terms as well as surface tension are neglected. Provided that is known and is given, the Rayleigh–Plesset equation can be used to solve for the time-varying gas bubble radius which includes collapses and rebounds.

Because we study the numerical solvers for the Rayleigh-Plesset equation, we assume that the pressure difference is constant, i.e. and for simplicity. The solution contains very strong gradients during the collapse phase and the following rebound phase, i.e it is a stiff ODE. In order to guarantee the accuracy of the numerical solution, very small time steps as well as a high order of the solver must be used. The Rayleigh-Plesset equation (9) is reformulated in a system of first order ODE’s with polynomial right-hand side

 dy1dt=y2,dy2dt=−(pB−p)ρy3−32y22,dy3dt=−y23y2, (10)

where , , , and the initial conditions are , , and .

Since the equation is a stiff ODE, we use an adaptive mapping approach. This means that we built Taylor maps of 7th order for different time intervals ranging from to . Then we apply maps based on the relative error during the simulation. For this approach, the polynomial neurons can be organized in a shallow architecture that provides possibility to simulate of a bubble motion with control of accuracy. To compare the computational performance we run simulation for three initial conditions up to the collapsing time. The adaptive mapping approach is 2.5 times faster then the solver for stiff differential equations ode45 from [20] with similar accuracy (see Fig. 4, 5).

### 3.3 The Burgers’ equation

Burgers’ equation is a fundamental partial differential equation that occurs in various areas, such as fluid mechanics, nonlinear acoustics, gas dynamics, and traffic flow. This equation is also often used as a benchmark for numerical methods. In [13], a feed-forward neural network is trained to satisfy Burgers’ equation and certain initial conditions, but the computational performance of the approach is not estimated. In this section, we demonstrate how to build a PNN that solves Burgers’ equation and does not require training to satisfy different initial conditions.

Burgers’ equation has a form

 ∂u(t,x)∂t+u(t,x)∂u(t,x)∂x=ν∂2u(t,x)∂x2. (11)

Following [1] for benchmarking, we use an analytic solution

 u1(t,x)=−2νϕ(t,x)dϕdx+4,ϕ(t,x)=exp−(x−4t)24ν(t+1)+exp−(x−4t−2π)24ν(t+1),

 un+1i−uniΔt+uniuni−uni−1Δx=νuni+1−2uni+uni−1Δx2, (12)

where stands for the time step, and stands for the grid node.

The equation (12) presents a finite difference method (FDM) that consists of an Euler explicit time discretization scheme for the temporal derivatives, an upwind first-order scheme for the nonlinear term, and finally a centered second-order scheme for the diffusion term. The time step for benchmarking is fixed to with the uniform spacing of , . Thus, for the numerical solution for times from to on , the method requires the mesh with 1000 steps on space coordinate and 2000 time steps.

It is indicated in [1] that the FDM introduces a dispersion error in the solution (see Fig. 6). Such error can be reduced by increasing the mesh resolution, but then the time step should be decreased to respect the stability constraints of the numerical scheme.

To build a PNN for solving the Burger’s equation we translate the eq. (11) to a system of ODEs and build a Taylor map in accordance with Sec. 2. Assuming that the right-hand side of the eq. (11) can be approximated by a function and considering this approximated equation as a hyperbolic one, it is possible to derive the system of ODEs

 ddt(XU)=(Uf(x,U)), (13)

where , , and is vector of discrete stamps on space. This transformation from PDE to ODE is well known and can be derived using the method of characteristics and direct method [12]. If is the same discretization as in (12), then the equation (13) leads to the system of 2000 ODEs

 x′i =ui,u′i =f(ν,ui+1,ui,ui−1,xi+1,xi,xi−1),

which can be easily expanded to the power series with respect to the and up to the necessary order of nonlinearity.

Since there are only first and second order approximations in the benchmarking FDM scheme, we built a Taylor map of only the first order in a time interval . This time step is five times larger than that used in the benchmarks. The numerical solution provided by the resulting PNN is presented in Fig. 6, and the accuracy and performance are compared in Tab. 1.

The PNN numerically estimates dynamics in a larger time interval and provides better accuracy with less computational time in comparison with FDM of the same order of derivative approximations. If the FDM scheme is adjusted to a higher accuracy, the computational time will be increased even more. Accuracy is calculated as the MSE metric between the numerical solution and its analytic form at final time .

Since the derived from the differential equation PNN is an approximation of the general solution, it can be used for simulation of the dynamics with different initial conditions without weight modification. For example, for another analytic solution

 u2(t,x)=(1+exp(0.5(x−0.5t)/ν)−1,

the same PNN provides a numerical solution in time with MSE error , while FDM yiels . The elapsed times are the same since the mesh sizes were not changed.

## 4 Learning dynamical systems from data

In this section, we discuss training of the PNN. We demonstrate how one can adjust weights of the PNN for certain initial conditions of the Burgers’ equation by additional training of the PNN derived from the equation. We also introduce a regularization of the PNN based on combinatorial binary optimization and demonstrate the proposed methods in a practical application of a data-driven system learning.

### 4.1 Training PNN for certain initial conditions

The PNN derived from differential equation works for any initial condition with the same level of accuracy. The accuracy depends on the size of the PNN, namely the order of non-linearities in the map (2

). But it is still possible to train PNN for certain initial conditions. This fitting increases accuracy only for this particular solution, but probably decrease accuracy for others.

In order to train PNN for the Burgers equations, we consider deep architecture with 400 layers with shared weights. Each layer M is defined by Taylor map (2) for time :

 u(t=0,x)→M→…→M→u(t=0.5,x).

The loss function is defined as the inconsistency of the numerical solution to differential equations. The initial values of weights are calculated from the PDE with corresponding Taylor map.

Initially, the PNN has weights with only of them that are not equal to zero. The point of training is to adjust these weights to achieve better accuracy for the given initial conditions. We implemented a simplified training based on coordinate descent method that allows to increase accuracy in 3 times, but in the same time requires lots of computational time.

### 4.2 Regularization of the PNN

It is common for physical systems, that map is presented by sparse matrices . So it is important to identify non-zero weight in the polynomial neural network. Classical and regularization terms do not work in this case since they tend to decrease weight. Weights in PNN can be large by absolute value but the total amount of non-zero weights should be as small as possible. To solve this issue, one can introduce a binary mask that is applied for weight matrices,

 Xi+1=(B0∗W0)+(B1∗W1)Xi+(B2∗W2)X[2]i+…, (14)

where means element-wise multiplication.

During the fitting of weight matrices , it is necessary also to find an optimal mask . If the weights

are fixed after each training epoch, the map (

14) becomes linear for the components of .

It is possible to formulate regularization as a high-order binary optimization problem and translate it into a quadratic binary minimization problem (QUBO) that can be run on Quantum Annealers. For example,the loss function

translates directly to QUBO, while loss function that introduced for Burgers equation translates to high-order unconstrained binary optimization problem. By introducing additional binary variables it is possible finally to formulate regularization as QUBO.

To demonstrate the advantage of QUBO–based regularization, we implemented a simplified example with the Van der Pole oscillator. The equation is widely used in the physical sciences and engineering and can be used for the description of the pneumatic hammer, steam engine, periodic occurrence of epidemics, economic crises, depressions, and heartbeat. The equation has well-studied dynamics and is widely used for testing of numerical methods, e.g. [23].

The Van der Pol oscillator is defined as the system of ODEs that can be presented in the form of

 x′ =y, (15) y′ =y−x−x2y.

We generate a training data as a particular solution of the system (15) with the initial condition . After this solution was generated, the equation is not used further.

Having this training data set (Fig. 9 ), the proposed neural network can be fitted with the mean squared error (MSE) as a loss function based on the norm

 ||Xi+1−M∘Xi||=||(xi+1yi+1)−W0−…−W3⎛⎜ ⎜ ⎜ ⎜ ⎜⎝x3ix2iyixiy2iy3i⎞⎟ ⎟ ⎟ ⎟ ⎟⎠||.

We implemented the above-described technique in Keras/TensorFlow and fitted a third-order Lie transform–based neural network with an Adamax optimizer. After each training epoch we applied QUBO–based regularization for the weights of the PNN.

Fig. 8 shows that QUBO–based regularization can significantly decrease the number of epochs required to achieve an appropriate level of accuracy. However, the time required to solve QUBO problem is longer then training without regularization. It should be noted that using of special hardware (e.g. quantum or digital annealers) instead of classical QUBO-solvers can potentially provide additional improvement in calculation time.

The generalization property of the network can be investigated by examining prediction not only from the training data set but also for new initial conditions. Fig. 10 demonstrates predictions that are calculated starting also at the new points , and that were not presented to the PNN during fitting. For the prediction starting from the training initial condition, the mean relative error of the predictions is . For the new initial conditions, the mean error is even less and equals to .

### 4.3 Data-driven identification of a zinc-air energy storage system

Zinc-air batteries are one of the most promising energy storage systems. An effective mathematical model is a key part of a battery management system, responsible for its operation control and internal state evaluation. The model-based analysis is an effective tool in the design and manufacture of power supplies, as well as in the analysis of their behavior under different operating conditions, cf. [17].

There are several approaches for building dynamic model of a battery. By now the researches were mostly focused on the zinc-air batteries mathematical models derivation based on electrochemical principles of its functioning. This implies utilization of complicated nonlinear PDEs, describing the battery behavior precisely, but requiring time-consuming numerical procedures for their solution.

Another type is a “black box” model or data-driven model that is formed by learning on experimental data using artificial neural networks. Since this approach does not imply the presence of state equations in the model, it works only on similar data to what the NN was trained on. This requires availability of operational data from real environment of sufficient volume and quality, which in practice cannot always be guaranteed.

Thus, state space models with lumped parameters are more suitable for analysis, control, and optimization of power supplies, since in this case, it includes battery general characteristics, like voltage, current, state of charge and so on. In the most general form without specification of any details, it can be presented as

 x(k+1)=f(x(k),u(k)),y(k)=g(x(k),u(k)). (16)

It should be noted that real power sources exhibits strongly nonlinear effects, that are to be included as a nonlinear part of and functions of mathematical model (16) from [22]. This implies manual derivation of the equations for dynamics description and identification of its parameters. In this work we take an alternative approach by data-driven construction of mathematical model from a set of controlled experiments [18]. We investigate whether it is possible to extract relevant features from current and voltage measurements collected during battery discharge in various fixed loading conditions.

For demonstration of this approach, we consider only voltage dynamics for discharge currents , , and and solve an identification problem only for regions where voltage decreases in time.

Initial data taken from the open source

[18] were filtered with simple moving average window and aligned to constant time stamp of 1ms. To present dynamics of the voltage decreasing, we use polynomial neural network with three layers (see Fig. 12). Each layer is represented by a Taylor map, namely and are Taylor maps of second order, is a Taylor map of fifth order.

To train the PNN, vector of voltage at the given time stamp and the value of discharge current is used as input data, and voltage after three time stamps is used as output data. Also, the described above regularization were used.

After training, the dynamics of the voltage decreasing is calculated as consequent prediction performed by PNN. Starting with initial voltage at time we calculate dynamics by formula

 V0→V1=PNN(V0)→…→Vk=PNN(Vk−1).

In this way, PNN plays role of the model of the dynamical system that can calculate dynamics with different initial conditions .

Fig. 11 shows the results of the voltage dynamics simulation with the trained PNN. While additional features like derivatives of voltage or power can potentially increase the accuracy, the built PNN preserves physical properties of the dynamics. For example, the PNN successfully predicts the lower level of voltage that is the same for all currents in contrast to simple function approximation, that can give aberrant negative values of voltage, cf. [22].

One of the advantage of the polynomial architectures in comparison with other machine learning methods is its coincidence with theory of dynamical systems and differential equations. A PNN trained with data can provide mathematical model with physical properties preservation. For example, in

[15] it is shown, that MLP and LSTM neural networks just memorized data even for simplified physical systems. They can not predict dynamics for data that have not been presented during fitting. PNN architecture allows to extrapolate dynamics on new initial condition and preserve physical properties of the dynamical systems.

## 5 Supplementary code

The program code that contains the considered examples can be found in the following GitHub repositories.

github.com/andiva/DeepLieNet: implementation of PNN in Keras/TensorFlow (sections 1, 3.3, 3.2, and 4.3):

/demo/Deflector: cylindrical deflector simulation
/demo/accelerator.ipynb: storage ring simulation
/demo/Ray_Ples: simulation of the Rayleight-Plesset equation
/demo/Battery: trained PNN for voltage dynamics

github.com/andiva/AQCC: PNN for simulation of the Burgers equation (Sec. 3.3), training PNN for the Burgers equation (Sec. 4.1), QUBO–based regularization and training of PNN for the Van der Pol oscillator (Sec. 4.2).

## 6 Conclusion

Since the Taylor maps and polynomial neural network have strong connection to systems of differential equations, these methods can be used for physical system simulation and data-driven identification. If the differential equations for the system are known, it is possible to calculate weights of the PNN with necessary level of accuracy. In this way, the built PNN can be used for simulation of the dynamics instead of the traditional numerical solvers. It is shown, that Taylor mapping and PNN provides the same or even better accuracy with less computational time in comparison to step-by-step integrating of the equations.

If the equations of the system are not known, the PNN can be trained by scratch with the available data. One can define an architecture of the PNN and fit weights directly with data. Instead of fitting of a parametrized equations, this approach allows to fit weights of the PNN directly.

In the article, we consider simulation of dynamical system arising in practical application with the PNN. For training PNN, we introduces QUBO–based regularization and demonstrated data-driven dynamics learning with both simplified Van der Pol oscillator and applied problem of the batteries development.

We do not consideed questions of accuracy, PNN architecture selection, and optimization methods for training in details. This work should be done in further research. With the provided examples, we rather demonstrated the applicability of Taylor mapping methods that are commonly used in dynamical systems investigation to the theory of neural network and data-driven system learning.

We would like to thank Prof. Dr. Sergei Andrianov for his help with explanation of Taylor mapping techniques for solving of system of differential equations.