Simulating SINDy-PI generated model of a higher dimensional ODE

[ainsleylaiuw]

Hello

I am trying to model a 2D ODE system using SINDy-PI with a PDE library to get models in the form of $\dot{x_0} = f_{0}(x_0, x_1), \quad \dot{x_1} = f_1(x_0, x_1)$ (like how Kadierdan does in his video) for simulation purposes. However, following the process in the example notebook, my generated model eqns tend to have both $\dot{x_0}$ and $\dot{x_1}$ in them, making them simplify to $\dot{x_0}=f(x_0, x_1, \dot{x_1})$ and unable to integrate.

Calling model.predict() on my training data and comparing corresponding $\dot{x}$ eqns to model.differentiate() tends to have similar results so I suspect that I'm missing something in my code or my model/library creation logic is flawed.
Any ideas on how to 'filter' out other derivative terms and get both derivatives rather than just one?

The system I am working on:

$$ \dot{u} = -u + \frac{2}{1+v^2} \longrightarrow \dot{u} + v^{2}\dot{u} + u + uv^{2} = 2$$

$$ \dot{v} = -v + \frac{1}{1+u^2} \longrightarrow \dot{v} + u^{2}\dot{v} + v + vu^{2} = 1 $$

And some of my beginning code:

import numpy as np
from scipy.integrate import solve_ivp
import pysindy as ps
import sympy as sp

# system definition
a, b = 2, 1
beta, gamma = 2, 2
dudt = lambda u, v: -u + (a / (1+v**beta))
dvdt = lambda u, v: -v + (b / (1+u**gamma))
netswitch = lambda t, x: np.array([dudt(x[0], x[1]), dvdt(x[0], x[1])])
var = ['u', 'v']
integrator_keywords = {}
integrator_keywords['rtol'] = 1e-12
integrator_keywords['method'] = 'LSODA'
integrator_keywords['atol'] = 1e-12

# training data creation
r = len(var) # appears to be # of var in example (x0, x1, x2,...)
dt = 0.01
T = 20
t = np.arange(0, T + dt, dt)
t_span = (t[0], t[-1])
netswitch0_train = np.random.random(r) * 5
netswitch_train = solve_ivp(netswitch, t_span, netswitch0_train, t_eval=t, **integrator_keywords).y.T

library_functions = [
    lambda x: x,
    lambda x, y: x * y,
    lambda x: x ** 2,
    lambda x, y: x * y ** 2,
    lambda x: x ** 3,
    lambda x, y: x * y ** 3,
    lambda x: x ** 4
]
x_dot_library_functions = [lambda x: x]
library_function_names = [
    lambda x: x,
    lambda x, y: x + y,
    lambda x: x + x,
    lambda x, y: x + y + y,
    lambda x: x + x + x,
    lambda x, y: x + y + y + y,
    lambda x: x + x + x + x
]
pde_library = ps.PDELibrary(
    library_functions=library_functions,
    temporal_grid=t,
    function_names=library_function_names,
    include_bias=True,
    implicit_terms=True,
    derivative_order=1,
    include_interaction=False#,
    #interaction_only=False
)
sindy_opt = ps.SINDyPI(
    threshold=1e-6,
    tol=1e-8,
    thresholder="l1",
    max_iter=20000,
)
model = ps.SINDy(
    optimizer=sindy_opt,
    feature_library=pde_library,
    differentiation_method=ps.FiniteDifference(drop_endpoints=True)
)
model.fit(netswitch_train, t=dt)
model.print()

[akaptano]

What equations do you find with your model? The analytic equations look like they can be written x0_dot = f0(x0, x1). You can also choose library terms such that x1_dot is not in the library for the x0_dot equation. Since you are doing SINDy-PI with PDEs, the newest main branch of PySINDy actually uses the PDELibrary for this (https://github.com/dynamicslab/pysindy/blob/master/examples/9_sindypi_with_sympy.ipynb) and this probably a better way to go.

[ainsleylaiuw]

Tweaking the code to use multiple trajectories (singular traj. gives less model eqns), full eqn output with that PDELibrary is:

1 = 0.715 x0 + 0.822 x1 + -0.016 x0x1 + -0.186 x0x0 + -0.069 x1x1 + -0.011 x0x1x1 + 0.022 x0x0x0 + 0.014 x1x1x1 + 0.002 x0x1x1x1 + -0.001 x1x1x1x1 + 0.111 x0_t + 0.600 x1_t
x0 = 1.257 1 + -0.988 x1 + -0.024 x0x1 + 0.323 x0x0 + 0.017 x1x1 + 0.058 x0x1x1 + -0.049 x0x0x0 + -0.007 x1x1x1 + -0.012 x0x1x1x1 + 0.002 x0x0x0x0 + 0.001 x1x1x1x1 + -0.076 x0_t + -0.860 x1_t
x1 = 1.124 1 + -0.769 x0 + 0.058 x0x1 + 0.176 x0x0 + 0.166 x1x1 + -0.018 x0x1x1 + -0.014 x0x0x0 + -0.046 x1x1x1 + 0.003 x0x1x1x1 + -0.001 x0x0x0x0 + 0.005 x1x1x1x1 + -0.130 x0_t + -0.652 x1_t
x0x1 = -0.456 1 + -0.397 x0 + 1.213 x1 + 0.540 x0x0 + -1.417 x1x1 + 0.701 x0x1x1 + -0.106 x0x0x0 + 0.318 x1x1x1 + -0.124 x0x1x1x1 + 0.010 x0x0x0x0 + -0.025 x1x1x1x1 + 0.778 x0_t + -1.464 x1_t
x0x0 = -2.335 1 + 2.304 x0 + 1.610 x1 + 0.235 x0x1 + 0.288 x1x1 + -0.296 x0x1x1 + 0.192 x0x0x0 + -0.046 x1x1x1 + 0.058 x0x1x1x1 + -0.014 x0x0x0x0 + -0.108 x0_t + 2.015 x1_t
x1x1 = -0.735 1 + 0.101 x0 + 1.299 x1 + -0.527 x0x1 + 0.246 x0x0 + 0.375 x0x1x1 + -0.075 x0x0x0 + 0.288 x1x1x1 + -0.068 x0x1x1x1 + 0.009 x0x0x0x0 + -0.030 x1x1x1x1 + 0.338 x0_t + -0.219 x1_t
x0x1x1 = -0.535 1 + 1.609 x0 + -0.646 x1 + 1.190 x0x1 + -1.156 x0x0 + 1.713 x1x1 + 0.228 x0x0x0 + -0.387 x1x1x1 + 0.189 x0x1x1x1 + -0.019 x0x0x0x0 + 0.029 x1x1x1x1 + -0.953 x0_t + 2.412 x1_t
x0x0x0 = 6.000 1 + -7.652 x0 + -2.923 x1 + -1.020 x0x1 + 4.230 x0x0 + -1.939 x1x1 + 1.287 x0x1x1 + 0.448 x1x1x1 + -0.255 x0x1x1x1 + 0.088 x0x0x0x0 + -0.031 x1x1x1x1 + 0.663 x0_t + -5.580 x1_t
x1x1x1 = 1.501 1 + -0.396 x0 + -3.498 x1 + 1.153 x0x1 + -0.384 x0x0 + 2.811 x1x1 + -0.828 x0x1x1 + 0.169 x0x0x0 + 0.154 x0x1x1x1 + -0.022 x0x0x0x0 + 0.120 x1x1x1x1 + -0.158 x0_t + -0.457 x1_t
x0x1x1x1 = 3.250 1 + -8.720 x0 + 2.745 x1 + -5.702 x0x1 + 6.110 x0x0 + -8.393 x1x1 + 5.125 x0x1x1 + -1.224 x0x0x0 + 1.954 x1x1x1 + 0.099 x0x0x0x0 + -0.150 x1x1x1x1 + 4.606 x0_t + -11.726 x1_t
x0x0x0x0 = -10.364 1 + 41.450 x0 + -15.612 x1 + 10.852 x0x1 + -34.406 x0x0 + 25.216 x1x1 + -11.868 x0x1x1 + 9.883 x0x0x0 + -6.511 x1x1x1 + 2.320 x0x1x1x1 + 0.566 x1x1x1x1 + -8.944 x0_t + 20.940 x1_t
x1x1x1x1 = -9.550 1 + 5.474 x0 + 24.340 x1 + -5.761 x0x1 + 0.074 x0x0 + -18.687 x1x1 + 3.988 x0x1x1 + -0.751 x0x0x0 + 7.701 x1x1x1 + -0.760 x0x1x1x1 + 0.123 x0x0x0x0 + -3.029 x0_t + 10.421 x1_t
x0_t = 2.115 1 + -0.824 x0 + -1.813 x1 + 0.516 x0x1 + -0.165 x0x0 + 0.602 x1x1 + -0.372 x0x1x1 + 0.046 x0x0x0 + -0.029 x1x1x1 + 0.066 x0x1x1x1 + -0.006 x0x0x0x0 + -0.009 x1x1x1x1 + 0.190 x1_t
x1_t = 1.147 1 + -0.935 x0 + -0.910 x1 + -0.097 x0x1 + 0.307 x0x0 + -0.039 x1x1 + 0.094 x0x1x1 + -0.039 x0x0x0 + -0.008 x1x1x1 + -0.017 x0x1x1x1 + 0.001 x0x0x0x0 + 0.003 x1x1x1x1 + 0.019 x0_t

As shown, the eqns have both $\dot{x_0}$ and $\dot{x_1}$ (I should probably also raise the optimizer threshold). Using Sympy to solve for either $\dot{x}$ leads to the other $\dot{x}$ being on the other side, causing trouble when integrating (at least I think this is the issue).

You can also choose library terms such that x1_dot is not in the library for the x0_dot equation

Constraining models to only have either $\dot{x_0}$ or $\dot{x_1}$ in a model eqn, as you suggest, would fix my problem. However, I'm not sure what to change to achieve this as lambda x, y allows interaction of variables $u$ and $v$

[akaptano]

you will need to use the GeneralizedLibrary class to achieve this (see Example 1 notebook towards the end).

[ainsleylaiuw]

Does the GeneralizedLibrary compare between two different libraries? I'm trying to understand the documentation on it and it seems like it helps omit unwanted library functions so if my basic PDELibrary is ['1', 'x0', 'x1', 'x0x1', 'x0x0', 'x1x1', 'x0x1x1', 'x0x0x0', 'x1x1x1', 'x0x1x1x1', 'x0x0x0x0', 'x1x1x1x1', 'x0_t', 'x1_t'], I can either choose to omit x0_t or x1_t in all of my models, which isn't what I'm aiming for.
Instead, in the SINDy-PI output with multiple model eqns, I want them to only have one derivative but test for both, if that makes sense.
i.e.

Model 0, 1 = 1 * x0 + 1 * x0x1 + 1 * x0_t
Model 1, x = 1 * x1 + x1_t
...
Model 12, x0_t = f0(x0, x1)
Model 13, x1_t = f1(x0, x1)

Alternatively, creating two models, one omitting x0_t, the other omitting x1_t might be better?

Sorry if this question seems ill-prepared, I'm getting an error from the Ex. 1 notebook GeneralizedLibrary section NotFittedError: This GeneralizedLibrary instance is not fitted yet. Call 'fit' with appropriate arguments before using this estimator, running the second cell.

[ainsleylaiuw]

Update on this: After reviewing Kaheman's SINDy-PI repo, solving for one derivative at a time seems to be the correct course of action. The NotFittedError was from using the pre-release version (to get temporal_grid for the PDELibrary as per Ex. 9 notebook). This was fixed by using the 1.7.1 release. PDELibrary loses temporal_grid in 1.7.1 so I switched to using custom libraries ('regular' functions 1, x0, x1, x0x0, ...) and SINDyPILibrary (x0_dot, x1_dot) and constraining to only using one xn_dot.

[akaptano]

Also for clarification,
model.differentiate -> produce x_dot from your training data
model.predict -> produce x_dot from your model (hence why it looks similar but not the same as model.differentiate)
model.simulate -> produce x from your model

[ainsleylaiuw]

Achieving the behavior of model.simulate is what I'm trying to get since I can't directly call model.simulate when using SINDy-PI due to it giving multiple models. Since I've been having trouble with the $\dot{x}$ solving, I've been using differentiate() and predict() as small check to verify the models.