Lie Transform Based Polynomial
Neural Networks for Dynamical Systems
Simulation and Identification
Andrei Ivanov, Sergei Andrianov, and Alena Sholokhova

arXiv:1802.01353v1 [] 5 Feb 2018

Saint Petersburg State University,
Faculty of Applied Mathematics and Control Processes, Russia
05x.andrey@gmail.com, sandrianov@yandex.ru, al.sholokhova@gmail.com

Abstract. In the article, we discuss the architecture of the polynomial
neural network that corresponds to the matrix representation of Lie transform. The matrix form of Lie transform is an approximation of general
solution of the nonlinear system of ordinary differential equations. Thus,
it can be used for simulation and modeling task. On the other hand, one
can identify dynamical system from time series data simply by optimization of the coefficient matrices of the Lie transform. Representation of
the approach by polynomial neural networks integrates the strength of
both neural networks and traditional model-based methods for dynamical
systems investigation. We provide a theoretical explanation of learning
dynamical systems from time series for the proposed method, as well as
demonstrate it in several applications. Namely, we show results of modeling and identification for both well-known systems like Lotka–Volterra
equation and more complicated examples from retail, biochemistry, and
accelerator physics.
Keywords: matrix Lie transform, deep polynomial neural network, dynamical systems simulation, time series based system identification

1

Introduction

Modeling and control of complex dynamical systems require techniques for nonlinearities and uncertainties consideration. On the face of it, artificial neural
networks could provide a suitable approach for dynamical systems ’learning’.
The applications cover different problems like solving of ordinary differential
equations [9,12], signal processing [10], feedback control systems [11], dynamical
systems modeling and identification [3] and others.
One group of researchers suggests solving a known differential equation with
predefined mathematical formulas by neural networks. Another one tries to utilize benefits of neural networks for system identification and inverse problem
solving when the equations of the system’s dynamics are not defined. In both
cases, the neural network is considered as a black-box model and its applicability
for dynamical systems modeling is explained only by the universal approximation
theorem.

2

Dynamical Systems and Polynomial Neural Networks

In the article, we describe a polynomial neural network, that is suitable for
simulation and identification of dynamical systems in a mathematically rigorous
style. The proposed architecture is a neural network representation of Lie propagator for dynamical systems integration that is introduced in [4] and commonly
is used in charged particles beam dynamics simulation. The neural network representation allows incorporating parallel architecture, performance increasing
and uncertainties accounting with the rigorous theory of dynamical systems investigation. We consider dynamical systems that can be described by nonlinear
ordinary differential equations:

d
X = F(t, X),
dt

(1)

where t is independent variable, X ∈ Rn is state vector. There is an assumption
that function F can be expanded in Taylor series with respect to the components
of X. Note that independent variable t can arise in the equation as an arbitrary
nonlinear function. Such nonlinear systems arise in different fields, namely automated control, robotics, mechanical and biological systems, chemical reactions,
drug development, molecular dynamics and so on.

X

+

+

+

+

Y

X[2]

X
X[3]
Fig. 1. Neural network representation of matrix Lie transform up to the 3rd order of
nonlinearities.

Dynamical Systems and Polynomial Neural Networks

2

3

Proposed Neural Network

Proposed neural network implements map M : X → Y using following polynomial transformation
Y = W0 + W1 X + W2 X[2] + . . . + Wk X[k] ,
where X, Y ∈ Rn , Wi are weights matrices, and X[k] means k-th Kroneker
power of vector X. For instance Fig. 1 presents neural network representation
of map M up to the third order of nonlinearities for 2 dimensional state space.
In each layer the input vector X = (x1 , x2 ) is consequently transformed into
X[2] = (x21 , x1 x2 , x22 ) and X[3] = (x31 , x21 x2 , x1 x22 , x32 ) where weighted sum is
applied. The output Y equals to sum of results from every layers. Note in the
figure we reduce Kroneker powers for decreasing of weights matrices dimension
(e.g. X[2] = (x21 , x1 x2 , x2 x1 , x22 ) → (x21 , x1 x2 , x22 ).
The transformation M can be considered as an approximation of evolution
operator of the (1) for predefined initial time t0 and time interval ∆t. This
means that with appropriate weights Wi = Wi (t0 , ∆t) the evolution of initial
state vector X0 = X(t0 ) during time ∆t can be approximately calculated as
Y = X(t0 , X0 , ∆t, ) = M ◦ X0 . If the system (1) is time independent then
weights Wi are constant for a predefined time interval.

3

Mathematical Background

The dynamics of vector X in (1) can be presented in form of Lie transform
Z t

M(t|t0 ) = T exp
LF (τ )dτ ,
t0

where LF (τ ) is Lie operator associated with vector function F in (1). Transformation M is presented in form of the time-ordered exponential operator and
can be identified with the dynamical system itself.
On the assumption that the function F allows its expansion in Taylor series, the required solution of equation (1) in its convergence region can be also
presented in form of series
X(t) = M ◦ X0 =

∞
X

Rk X[k] ,

F=

k=0

∞
X

Pk X[k] .

(2)

k=0

In [2] it is shown how to calculate matrices Rk either analytically or numerically. The main idea is replacing differential equation (1) by the equation
Pk
0
Rik
(t|t0 ) = j=i Pij (t)Rjk (t|t0 ), 1 ≤ i < k, where Pij = P1(j−i+1) P(i−1)(j−1) ,
P1k = Pk , R1k = Rk .
Another way to estimate matrix coefficients (Wk = Rk ) of a truncated map
[2]

[k]

X(t) = W0 (t) + W1 (t)X0 + W2 (t)X0 + . . . + Wk (t)X0 ,

(3)

4

Dynamical Systems and Polynomial Neural Networks

is utilizing an appropriate numerical step-by-step integration method. Taking
derivative of the X(t) with respect to the (2) one can obtain a system of equations
d
d
d
d
d
[2]
[k]
X = W0 (t) + W1 (t)X0 + W2 (t)X0 + . . . + Wk (t)X0 ,
dt
dt
dt
dt
dt
d
X =P0 (t) + P1 (t)X + P2 (t)X[2] + . . . + Pp (t)X[p] =
dt
P0 (t)+


[2]
[k]
P1 (t) W0 (t) + W1 (t)X0 + W2 (t)X0 + . . . + Wk (t)X0 +

[2]
[2]
[k]
P2 (t) W0 (t) + W1 (t)X0 + W2 (t)X0 + . . . + Wk (t)X0
+ ...+

[p]
[2]
[k]
Pp (t) W0 (t) + W1 (t)X0 + W2 (t)X0 + . . . + Wk (t)X0
,
which leads to a new system of ordinary differential equations with respect to
the weight matrices Wk
p

X
∂X[i]
d
Wk (t) =
Pi (t)
,
[i]
dt
∂(X )T
i=1

k = 0, 1, 2, . . .

(4)

0

Since the right-hand sides of the given equations depend only on Wi , the
system can be numerically integrated with initial condition Wk (0) = 0, k 6=
1; W1 (0) = E during necessary time interval.
To clarify this approach let’s consider 1-dim equation
x0 = kx2 ,
with the known general solution
x(t) =

x0
;
1 − ktx0

x0 = x(0).

(5)

Let’s calculate matrix Lie map using described above technique. This means
that we want to find out coefficients of the transformation
x(t) = W0 + W1 x0 + W2 x20 + W3 x30 + . . .
Taking derivative of x(t) yields the equations
x0 = W00 + W10 x0 + W20 x20 + W30 x30 + . . .
x0 = k W0 + W1 x0 + W2 x20 + W3 x30 + . . .

2

,

wich after combining like terms lead to the initial value problem for Wi
W00 = kW02 ,
W10 = 2kW0 W1 ,
W20 = k(2W0 W2 + W12 ),
W30 = k(2W0 W3 + 2W1 W2 ), . . .

Dynamical Systems and Polynomial Neural Networks

5

with initial condition W1 (0) = 1; Wi (0) = 0, i 6= 1. Consequently solving given
equatons one can found map
W0 = 0,
W1 = 1,
W2 = kt,
W3 = k 2 t 2 , . . .
and the estimation of general solution
x = x0 + ktx20 + k 2 t2 x30 + . . .
As we say above this series equals to the general solution inside the region of
convergence. Indeed, taking into account formula (1 − q)−1 = 1 + q + q 2 + q 3 + . . .
the precise solution (5) can be written in the same form
x0
x(t) =
= x0 + ktx20 + k 2 t2 x30 + . . .
1 − ktx0
3.1

Invariants preserving

Any numerical computational process leads to distortion of qualitative properties
(e.g. dynamical and kinematical invariants). These quantities can be evaluated
using, for example, Casimir operators. According to the Lie groups theory, we
can construct invariants using special forms and use these data for computational
process control.
As an example let’s consider Hamiltonian systems that are very popular in
physics problems. The Hamiltonian nature leads us to preserve of the symplecticity of the map M(t|t0)
M T (t|t0 )J0 M (t|t0 ) = J0 ,



0 E
and M (t|t0 ) = ∂ (M(t|t0 ) ◦ X0 ) /∂XT0 .
−E 0
For truncated map (3) one can apply order-by-order symplectification scheme
[1]. This leads to linear algebraic homogeneous equations for matrix elements
k
Wk = {wij
}. For instance for 2 dimensional state vector X = (x, y) and second
order map

where J0 =

[2]

X = W1 X 0 + W2 X 0 =



1
1
w11
w12
1
1
w21 w22

 2 
   2 2 2 
x
x0
w11 w12 w13  0 
x0 y0
+
2
2
2
y0
w21 w22 w23
y02

the describe above symplectic conditions are
1
1
1
1
w11
w22
−w12
w21
= 1,
1
2
1
2
1
2
1
2
w11
w22
−w21
w12
+ 2w22
w11
− 2w12
w21
= 0,
1
2
1
2
1
2
1
2
w22
w12
−w12
w22
+ 2w11
w23
− 2w21
w13
= 0,
2
2
2
2
2
2
2
2
w11
w23
−w13
w21
= 0, w12
w23
− w13
w22
= 0.

6

Dynamical Systems and Polynomial Neural Networks

This means that some of the elements of matrices Wk are coupled with each
other and whatever one computes these elements the described above condition
should be satisfied. To clarify this statement, let’s consider linear system
x0 = −y,
y 0 = x.
One can estimate numerically linear map (3). Solving equations (4) by Runge–
Kutta 4th order method with constant step size 0.02 results the map M(∆t =
0.1) in form of


0.99517016 −0.0998484
W1 =
0.0998484 0.99517016
Since the integration step is quite big the resulting map W1 is not precise
enough. This leads to physical property violation. However, we can use described
above symplectic condition for linear map and correct weights

W̃1 =


0.99517016 −0.09816492
.
0.09816492 0.99517016

Fig. 2 demonstrates the simulation results of this linear system. Note that
both W1 and W̃1 maps correspond to the same time interval ∆t = 0.1. Weights
of W̃1 are slightly corrected to det(W̃1 ) = 1, while det(W1 ) = 1.0003.

Fig. 2. Phase space of a simple linear system. Red lines correspond to the Runge–Kutta
based integration. Left image represents simulation by Lie map W1 with oscillation pattern violation. Right image shows application of the map W̃1 after its symplectification.

Thus, having the system of differential equations one can build matrix Lie
map up to the necessary order of nonlinearity. This map can replace step-by-step
integration method and be used for simulation tasks. On the other hand, matrix
Lie map can be represented as a polynomial neural network and be a model for
learning dynamical systems from data.

Dynamical Systems and Polynomial Neural Networks

4
4.1

7

Learning simple models
Simple Dynamical Systems Simulation

The proposed Lie transform based neural network can be used as a computer
model for dynamical systems and a numerical approach for dynamics simulation.
Having the differential equations one can build matrix Lie map and use it for
dynamics calculation. As an example, we consider three well-known systems. The
parameters of the equation are selected in an arbitrary way just for example.
Lotka–Volterra
equation
x0 = −y − xy,
y 0 = x + xy

Van der Pol
oscillator
x00 = (1 − x2 )x0 − x
(y = x0 , y 0 = x00 )

Saddle Point
x0 = 0.5x + 0.5y,
y 0 = 1.5x − 0.5y

Fig. 3. Phase space simulation by Lie based neural network (blue points) and traditional Runge–Kutta based method (red line) for Lotka–Volterra equation (left), Van
der Pol oscillator (center), and Saddle point.

For all equations, we compare the simulation results (see Fig. 3) of both
Lie transform of 3rd order of nonlinearity and traditional Runge–Kutta based
integration. Lie map is built in according to the described above algorithms in
0.01 time interval. Saddle point represents a linear system


1.005 0.005
W0 = 0; W1 =
; W2 = 0; W3 = 0.
0.015 0.995
The Lotka–Volterra system is presented by nonlinear equations. Therefore
map will consist of higher order weight matrices. Van der Pol oscillator is also a
nonlinear system that produces a limit cycle where trajectories in phase space
spiral into it. The corresponded nonlinear Lie map approximately equals to


0.99995067 0.01004917
W0 = 0; W1 =
; W2 = 0;
−0.01004917 1.00999984


1.6 · 10−7 −4.9 · 10−5 −3.2 · 10−7 −7.9 · 10−10
W3 =
.
4.9 · 10−5 −1.0 · 10−2 −9.9 · 10−5 −3.3 · 10−7

8

Dynamical Systems and Polynomial Neural Networks

Fig. 4. Train and test data in phase space states for Lotka–Volterra system identification.

4.2

Identification of Lotka–Volterra system

The simplest generalization property can be examined with system of nonlinear
oscillations. For instance, we considered differential equation
x0 = −2x + xy,
y 0 = y − xy,
that describes the predator-prey population dynamics with nonlinear oscillations (see Fig. 4). Starting with an initial point X0 we calculated discrete states
{Xi }i=1,n by numerical integrating of the differential equation. Then the differential equation was not utilized in training procedure and identification process.

Fig. 5. Train and test data in the space states for Lotka–Volterra system identification.

Dynamical Systems and Polynomial Neural Networks

9

According to the equations, we generated 4 different particular solutions
that are presented as time series X(ti ) = (xi , yi ). As a training set, we used only
one of them (blue line in the Fig. 4, 5). For testing, three other solutions were
selected, namely oscillation with higher amplitude, near linear oscillation, and
fixed point. We constructed three neural network architectures: polynomial one
based on 3rd order Lie transform, multilayer perceptron with sigmoid activation
functions (MLP), and long short-term memory network (LSTM). As inputs and
outputs for each neural network pairs (Xi , Xi+1 ) were used.

Fig. 6. Results of system identification and simulation in unknown state space by
Lie based neural network (left), MLP (center), and LSTM (right). Red points are
predictions.

After fitting on the training particular solution, both MLP and LSTM networks tend to learn the given track without any kind of generalization (see
Fig. 6). At the same time proposed Lie based neural network can predict dynamics in both nonlinear area and near linear oscillation around stationary point.
Moreover with the fitted matrix Lie map we can even estimate the stationary
point simply solving numerically equation Xf ixed = M ◦ Xf ixed .
The point of this example is the demonstration of generalization property
of Lie based neural network for physical data learning. Lie transform based
polynomial neural network corresponds to a system of differential equations and
recognizes the dynamical pattern from a single example. In the example, a Lie
transform based neural network was fitted on the oscillation pattern, and tends to
predict oscillations in other initial conditions. MLP and LSTM network require
lots of data to correctly predict dynamics and tend to predict known solution
from the training set.

10

5
5.1

Dynamical Systems and Polynomial Neural Networks

Applications
Retail Price Identification

In the article [6] the authors investigate dynamics of iPhone and iPad sales with
differential equations. They suggest analytical formulas for systems of nonlinear
differential equations and fit parameters based on time-series data. This is a traditional approach for system identification. Using the described above technique
one can identify dynamics utilizing matrix Lie transform without knowledge of
appropriate equations at all.
For the example, one can generate data (two time series for iPad and iPhone
sales) similar as in the article [6] and just fit these time series using 3rd order
Lie based neural network. The built neural network can be used for the products
specific market potential evaluation as well as other system parameters investigation. Note that in this case, we did not use any specific assumption on the
possible view of equations as made in the article [6].
70
60
50
40
30
20
10
0
0

5

10

15

20

25

30

35

Fig. 7. Identification of iPhone and iPad sales with Lie transform based neural network
of 5th order of nonlinearity.

After fitting of neural network, one can receive a Lie map
  
 
 
x
3.1838244e − 05
0.6154481 −0.36771438
x0
=
+
+
y
3.8284583e − 03
−0.16011019 0.45261806
y0
 2 
x0



−0.36270088 −0.2199135 −0.2199135 −1.2503444 
x
 0 y0  + . . . ,
−0.09446148 −0.26765344 −0.26765344 −1.0930712 y0 x0 
y02
where x0 , y0 are initial sales for iPhone and iPad respectively. Note that during
fitting coefficients with respect to the terms x0 y0 and y0 x0 are equal to each other
since they have the same impact on optimization criteria. The built polynomial
neural network is a white-box model of the dynamical systems. Moreover, it is
possible to reconstruct differential equation that corresponds to this map using
the inverse problem solving for the algorithm of map calculation.

Dynamical Systems and Polynomial Neural Networks

5.2

11

Biochemical Reaction Simulation

In this example, we provide an analysis of a biochemical system that is described
in [5] and represents the influence of the Raf Kinase Inhibitor Protein (RKIP).
In the article, the influence of RKIP is investigated via numerical analysis of
nonlinear ordinary differential equations using Matlab ode45 function that is
based on step by step integration (see Fig. 8).

Fig. 8. Simulation of RKIP dynamics using Runge–Kutta based method. Figure is
taken from [5].

Instead of using a step-by-step numerical integration method one can build a
polynomial neural network and utilize it for system simulation. Moreover, based
on this model one can further solve identification problem and clarify the form
of equations based on measured data.
The described in the article system of differential equation consists of 11
nonlinear equations that describe the biochemical network.We built a 2nd order
Lie map for this system and used it for simulation with initial condition from [5].
The results of the simulation are slightly different to [5] due to the application
of different integration methods and orders of nonlinearity.

Fig. 9. Simulation of biochemistry reaction utilizing Lie transform based neural network of 2 order of nonlinearity.

12

Dynamical Systems and Polynomial Neural Networks

5.3

Charged Particle Accelerators

Charged particle accelerator consists of a number of physical equipment (e.g.
quadrupoles, bending magnets and others, see Fig. 10). Design of accelerators
and nonlinear dynamics investigation require accurate computer model of such
complicated system.

Fig. 10. Schematic representation of circular accelerator and it’s Lie map based.

Each of the physical equipment can be described by a system of differential
equation that has a complex nonlinear form. For instance, the equation of radial
motion looks like

qH
2
(H(Ex − x0 Ez )/v− (1 + x0 )By + y 0 (x0 Bx + Bs ) ,
x00 =
m0 γv
where electromagnetic fields and particle state vector are incorporated.
For some problems, such as modeling of long-term dynamics, the traditional
step-by-step integration methods are not suitable due to the performance limitation. Instead of solving differential equations directly one can estimate matrix
Lie map for each physical element in an accelerator. Combining such maps consequently one can obtain a deep polynomial neural network representation of
whole accelerator ring (see. Fig. 11).
The proposed Lie transform based mapping approach was used, for example,
for nonlinear dynamics investigation of spin-orbit dynamics simulation in EDM
search project (see for example [14]). The articles [7, 8, 13] describe problem
formulation and simulation results that are achieved with the application of the
matrix Lie maps for simulation of the systems of nonlinear differential equations
in accelerated physics.

6

Conclusion

In the article, we demonstrate how one can build neural network representation
of dynamical systems. The greatest advantage of the proposed Lie transform
based neural network is its equivalence to the differential equations. This neural network is not a black-box model since each weight matrix corresponds to
the certain order of nonlinear effects in the real system and has a physical explanation. Another important issue is the ability of proposed architecture easy

Dynamical Systems and Polynomial Neural Networks

+



X0

X1



+

+



+

+



X2

X3





+



+



X4

Xn-1







+

+



+

+



X[2]



X[2]







X[3]

drift

sextupole

Xn





X[2]

13

X[3]

drift

bend magnet

...

quadrupole

Fig. 11. Lie transform based neural network for circular accelerator representation.

guarantee physical properties of the system by additional constraints on weights.
Note also that if the form of the differential equation is known, one can compute
weight matrices and utilize them as an initial point for optimization procedure
in neural network fitting.
Neural network representation of Lie transform combines the benefit of deep
learning techniques and traditional modeling theory for dynamical systems identification. It allows to avoid numerical iterative time-stepping methods and use
mapping approach instead of them. As soon as the proposed neural architecture
has a good coincidence with traditional modeling methods, it can be useful for
investigation of dynamics in unknown parameter space. Since weight coefficients
in the proposed neural architecture describe dynamical systems, they can be utilized as features for other machine learning models. For example, one can build
matrix representation of Lie transform for time-series data in unsupervised mode
and use weight coefficients for supervised methods.
The described above examples are implemented in Keras\Tensorflow using
Python. The code can be found at GitHub repository https://github.com/
andiva/DeepLieNet. Directory core consists of Keras layer that implements
matrix Lie transform up to the necessary order of nonlinearity, as well as an
algorithm for matrix Lie map estimation based on a predefined system of differential equations. Directory demo contains realizations of described above
examples.

References
1. Andrianov, S.: Order-by-order symplectification of truncated lie maps. Proc. of the
2001 Particle Accelerator Conference, Chicago, USA 6244/2010, 1787–1789 (2001)
2. Andrianov, S.: A role of symbolic computations in beam physics. Computer Algebra
in Sc, Comp., Lecture Notes in Computer Science 6244/2010, 19–30 (2010)
3. Chen, S., Billings, S.A.: Neural networks for nonlinear dynamic system modelling
and identification. International Journal of Control 56(2), 319–346 (1992)
4. Dragt, A.: Lie Methods for Nonlinear Dynamics with Applications to Accelerator
Physics (2011), inspirehep.net/record/955313/files/TOC28Nov2011.pdf
5. Gilbert, D., Heiner, M.: From petri nets to differential equations - an integrative
approach for biochemical network analysis. Lecture Notes in Comp. Sc. pp. 181–200
(01 2006)

14

Dynamical Systems and Polynomial Neural Networks

6. Guidolin, M., Guseo, R.: On product cannibalization. a new lotka-volterra model
for asymmetric competition in the icts. Tech. rep. (2016), http://paduaresearch.
cab.unipd.it/9738/1/GUIDOLINGUSEO_SETT2016.pdf
7. Ivanov, A., Andrianov, S., Senichev, Y.: Simulation of spin-orbit dynamics in storage rings. Journal of Physics: Conference Series 747(1) (2016)
8. Ivanov, A., Senichev, Y.: Matrix integration of odes for spin-orbit dynamics simulation. Proc. of the Particle Accelerator Conference pp. 400–402 (2014)
9. Lagaris, I.E., Likas, A., Fotiadis, D.I.: Artificial neural networks for solving ordinary and partial differential equations. Tech. rep. (1997), https://arxiv.org/
pdf/physics/9705023.pdf
10. Lapedes, A., Farber, R.: Nonlinear signal processing using neural networks: Prediction and system modelling. In: Langley, P. (ed.) IEEE International Conference
on Neural Networks. San Diego, CA, USA (1987)
11. Lewis, F.L., Ge, S.S.: Neural Networks in Feedback Control Systems, in Mechanical
Engineers’ Handbook: Instrumentation, Systems, Controls, and MEMS, vol. 2.
John Wiley and Sons, Hoboken, NJ, USA, third edition (ed m. kutz) edn. (2005)
12. Mall, S., Chakraverty, S.: Comparison of artificial neural network architecture in
solving ordinary differential equations. Advances in Artif. Neural Syst. (2013)
13. Senichev, Y., Ivanov, A., Lehrach, A., Maier, R., Zyuzin, D., Andrianov, S.: Spin
tune parametric resonance investigation. Proc. of the Particle Accelerator Conference pp. 3020–3022 (2014)
14. Senichev, Y., Lehrach, A., Maier, R., Zyuzin, D., Berz, M., Makino, K.,
S.Andrianov, Ivanov, A.: Storage ring edm simulation: methods and results. Proceedings of ICAP2012, Rostock, Germany pp. 99–103 (2012)

