dynamicslab / pysindy

A package for the sparse identification of nonlinear dynamical systems from data
https://pysindy.readthedocs.io/en/latest/
Other
1.36k stars 304 forks source link

Documentation for Handling Trigonometric Functions in SINDy-PI Library Functions #498

Open ibrahimw1 opened 1 month ago

ibrahimw1 commented 1 month ago

Firstly, thank you for making this package. I'm currently working on a research project involving the application of SINDy-PI from the PySINDy library. My specific focus is on modeling a double pendulum system, which involves equations with trigonometric functions such as sine and cosine.

The current documentation for PySINDy, especially higher order ODEs, lacks clear guidance for handling library functions involving trigonometric operations like sine and cosine(also lacking general documentation for trig equations like double pendulum). While the provided examples and documentation cover high-order ODEs, they do not address scenarios where ODEs involve trigonometric terms(eg pendulum on cart and double pendulum).

I am trying to replicate the results of the double pendulum problem outlined in the SINDy-PI paper "SINDy-PI: A Robust Algorithm for Parallel Implicit Sparse Identification of Nonlinear Dynamics", where the equation that was found by SINDY-PI was:

image

So it would be very helpful if you could provide library functions for this problem. Please check out the code I have so far and provide any feedback:

import numpy as np
import math
from matplotlib import pyplot as plt
import pysindy as ps
from pysindy.differentiation import FiniteDifference
import sympy as sp
from scipy.integrate import solve_ivp
from pysindy.utils.odes import double_pendulum

dt = 0.001
x0_train = [np.pi + 1.2, np.pi - 0.6, 0, 0]
t_train = np.arange(0,10,dt)
t_train_span = (t_train[0], t_train[-1])
x_train = solve_ivp(double_pendulum, t_train_span, x0_train, t_eval=t_train).y.T

library_functions = []
library_function_names = []

# Create the PDELibrary
sindy_library = ps.PDELibrary(
    library_functions=library_functions,
    temporal_grid=t_train,
    function_names=library_function_names,
    include_bias=True,
    implicit_terms=True,
    derivative_order=2
)

lib = sindy_library.fit(x_train)
lib.transform(x_train)

print("With function names: ")
print(lib.get_feature_names(), "\n")

sindy_opt = ps.SINDyPI(
    threshold=0.0001,
    tol=1e-5,
    thresholder="l1",
    max_iter=6000,
)

model = ps.SINDy(
    optimizer=sindy_opt,
    feature_library=sindy_library,
)

model.fit(x_train, t=t_train)
model.print()

Additional info - when i include library functions such Polynomial tensored with Fourier, the model ends up outputting 100-200 equations, none of them giving me the desired output from the original SINDy-PI paper. If i don't include any library functions, most coefficients become nan(0).

Jacob-Stevens-Haas commented 1 month ago

Hey, thanks for using the package! One of the challenges with SINDy-PI is that the number of models it searches grows combinatorially. The other issue is that there's often multiple ways of representing complex problems in pysindy. My guess is that the result in this paper required a lot of engineering.

That said, it's best if you can aggressively limit the number of features first to make sure you data generation is correct, then work more generally. Here's some suggestions:

  1. Add a control variable equal to $\phi_1 - \phi_2$
  2. Put them in a generalized library. Use inputs_per_library to make sure that they only act on the control variable.
  3. Try solving it without the polynomial library by explicitly calculating $\dot \phi_1$ and $\dot\phi_2$. You can then treat $\phi_1$ and $\phi_2$ as control variables.
  4. print lib.get_feature_names() to know what features to expect.

Admittedly, this is the part of the library that I understand the least.