search

andreikitaitsev / SymINDy

Symbolic Identification of Non-linear Dynamics. The method generalizes the SINDy algorithm by combining sparse and genetic-programming-based symbolic regression.
MIT License
72 stars 6 forks source link

readme

SymINDy - Symbolic Identification of Nonlinear Dynamics

This library is a generalization of SINDy, to be used for the reconstruction of dynamical systems with strong nonlinearities, which require the introduction of a combinatorial search in the elementary functions associated with the linear regression part of SINDy.

Video Tutorial

Click here to see the Video Tutorial:

IMAGE ALT TEXT HERE

About

Let's simulate some well-known dynamical systems and try to reconstruct them with SymINDy. We will use the following dynamical systems:

Project architecture

We follow scikit learn conventions and best practices for the development of python packages. The source code is located in the src directory. It consists of several modules:

Analogously to scikit learn estimators, SymINDy main class implements the methods fit, predict, score and _plottrees.

Linear Damped SHO

_To reconstruct the figure below run the script reconstruct_linear_damped_sho.py from simindy.validation package._

We simulate Linear Damped SHO on a specific time range, separate the simulated data into the train and test set (train-test ratio 0.7: 0.3). Then we fit SymINDy instance on the train set and call the predict method on the test set. The predicted data is plotted along with the original data. Running the script prints the reconstructed equation of the dynamical system to std output.

Estimated library functions: f0: x0
Estimated library functions: f1: x1
Estimated library functions: f2: mul(x0, x0)
Estimated library functions: f3: mul(x0, sin(x0))
Estimated library functions: f4: cos(x0)
Estimated dynamical system:

(y0)' = -0.100 f0(x0,x1) + 2.000 f1(x0,x1)
(y1)' = -2.000 f0(x0,x1) + -0.100 f1(x0,x1)

The original equation:

y0 = -0.1 * x0 + 2 * x1
y1 = -2 * x0 - 0.1 * x1

<div style="margin:2%";>
<img src="src\validation\figures\linear_damped_SHO.svg"; alt="linear_damped_SHO"; width=50%;/>

We correctly reconstruct the numeric equation and predict the test data. The coefficient of determination between original and reconstructed data is 1.0.

Cubic Damped SHO

_To reconstruct the figure below run the script reconstruct_cubic_damped_sho.py from simindy.validation package._

Same as before, we simulate Cubic Damped SHO on a specific time range, separate the simulated data into the train and test set (train-test ratio 0.7: 0.3). Then we fit SymINDy instance on the train set and call the predict method on the test set. The predicted data is plotted along with the original data. Running the script prints the reconstructed equation of the dynamical system to std output.

Estimated library functions: f0: mul(mul(x1, x1), x1)
Estimated library functions: f1: add(x0, x1)
Estimated library functions: f2: add(x1, add(x0, x0))
Estimated library functions: f3: mul(x0, mul(x0, x0))
Estimated library functions: f4: mul(x0, add(x0, x1))
Estimated dynamical system:

(y0)' = 2.000 f0(x0,x1) + -0.100 f3(x0,x1)
(y1)' = -0.100 f0(x0,x1) + -2.000 f3(x0,x1)

The original equation:

y0 = -0.1 * x0**3 + 2 * x1**3

y1 = -2 * x0**3 - 0.1 * x1**3

<div style="margin:2%";>
<img src="src\validation\figures\cubic_damped_SHO.svg"; alt="cubic_damped_SHO"; width=50%;/>

Again, we correctly reconstruct the numeric equation and predict the test data. The coefficient of determination between original and reconstructed data is 1.0.

Lorenz attractor

_To reconstruct the figure below run the script reconstruct_lorenz.py from symindy.validation package._

Same as before, we simulate Lorenz Attractor on specific a time range, separate the simulated data into the train and test set (train-test ratio 0.7: 0.3). Then we fit SymINDy instance on the train set and call the predict method on the test set. The predicted data is plotted along with the original data. Running the script prints the reconstructed equation of the dynamical system to std output.

Estimated library functions: f0: mul(x1, x0)
Estimated library functions: f1: x0
Estimated library functions: f2: x1
Estimated library functions: f3: x2
Estimated library functions: f4: mul(x0, x2)
Estimated dynamical system:

(y0)' = -10.000 f1(x0,x1,x2) + 10.000 f2(x0,x1,x2)
(y1)' = 28.000 f1(x0,x1,x2) + -1.000 f2(x0,x1,x2) + -1.000 f4(x0,x1,x2)
(y2)' = 1.000 f0(x0,x1,x2) + -2.667 f3(x0,x1,x2)

The original equation:

(y0)' =  10*(z1 - z0),
(y1)' =  z0*(28 - z2) - z1,
(y2)' =  z0*z1 - 8/3*z2

<div style="margin:2%";>
<img src="src\validation\figures\lorenz.svg"; alt="Lorenz"; width=50%;/>

Again, we correctly reconstruct the numeric equation and predict the test data. The coefficient of determination between original and reconstructed data is 1.0.

Thus, we accurately reconstruct dynamical linear and non-linear dynamical systems using symbolic regression to look for basis (library) functions. However, the main advantage of SymINDy is revealed with the reconstruction of highly non-linear systems.

Myspring (nonlinearly perturbed oscillator)

_To reconstruct the figure below run the script reconstruct_myspring.py from simindy.validation package._ Again, the fitting procedure is same as above. However, this time we focus on the original system first.

(x0)' = x1
(x1)' = - -4.518 x0 - 0.372 x1 + 9.123*sin(x0**2)

As we can see, the argument passed to sine is also non-linear. This makes it impossible for SINDy to reconstruct the system. But SymINDy can easily do it. Let's run the code reconstruct_myspring.py and see how it performs!

Reconstructed equation:

Estimated library functions: f0: x1
Estimated library functions: f1: x0
Estimated library functions: f2: cos(x1)
Estimated library functions: f3: sin(mul(x0, x0))
Estimated library functions: f4: sin(mul(x1, x0))
Estimated dynamical system:

(x0)' = 1.000 f0(x0,x1)
(x1)' = -0.372 f0(x0,x1) + -4.518 f1(x0,x1) + 9.123 f3(x0,x1)

Bingo, we did it! Let's see the graph.

<div style="margin:2%";>
<img src="src\validation\figures\myspring.svg"; alt="myspring"; width=50%;/>

We have correctly reconstructed non-linearly perturbed oscillator.

But let us not take anything for granted: maybe SINDy would do the same job. Let's reconstruct mysping with both SINDy and SymINDy and compare the results.

SymINDy and SINDy

_To reconstruct the figure below run the script SINDy_vs_SymINDy.py from simindy.validation package._

<div style="margin:2%";>
<img src="src\validation\figures\symindy_vs_sindy.svg"; alt="myspring"; width=50%;/>

Well, probably our statement above still holds: SINDy cannot estimate non-linearly perturbed oscillator, while SymINDy, as we have seen above, accurately recovers the underlying equation.

Summary

SymINDy is a new algorithm for the reconstruction of non-linear dynamics. It uses symbolic regression and SINDy algorithm to recover the systems of equations from time-series observations. It is free from the linearity assumption and thus is able to reconstruct systems unreachable for SINDy.

SymINDy can be applied to multiple theoretical and applied problems from blood dynamics to financial forecasting.



Installing the package

Option 1

Running the shell script

bash install symindy

Option 2

  1. Create a new python virtual environment
python -m venv env
  1. Activate virtual environment
source env/Scripts/activate
  1. Downgrade setuptools to version 57.0.0 (required by DEAP when using python > 3.7)
pip install setuptools==57.0.0
  1. Install the requirements using pip
pip install -r requirements.txt
  1. Install the SymINDy package
pip install -e .

Option 3

Run the package inside of the docker container built from Dockerfile.

docker build -t SymINDy .

Note, you shall have the docker client and daemon installed on your machine.

Commands

Dear developer, before pushing, remember to run

 pre-commit run --all-files

Notes

Relevant works