Extracting dynamics and faster simulation

[ricopicone]

Is your feature request related to a problem? Please describe.

This feature request attempts to solve two issues with pySINDy:

1. The lack of a way to extract the dynamics as an evaluable function, usable for instance, in scipy.integrate or control.NonlinearIOSystem.
2. The very slow pysindy.pysindy.SINDy.simulate() method.

Describe the solution you'd like

I'd like there to be a method in pysindy.pysindy.SINDy that returns the dynamics as an evaluable function, either in the form used by scipy.integrate or control.NonlinearIOSystem (an option for either would be nice). This could also be used in the slow pysindy.pysindy.SINDy.simulate() method to significantly improve performance. The current technique has very slow integration performance.

Describe alternatives you've considered

I've written my own parser:

def extract_sindy_dynamics(sindy_model, eps=1e-12):
    """Extract SINDy dynamics"""
    variables = sindy_model.feature_names  # e.g., ["x", "y", "z", "u"]
    coefficients = sindy_model.coefficients()
    features = sindy_model.get_feature_names()  
        # e.g., ["1", "x", "y", "z", "u", "x * y", "x * z", "x * u", "y * z", ...]
    features = [f.replace("^", "**") for f in features]
    features = [f.replace(" ", " * ") for f in features]
    def rhs(coefficients, features):
        rhs = []
        for row in range(coefficients.shape[0]):
            rhs_row = ""
            for col in range(coefficients.shape[1]):
                if np.abs(coefficients[row, col]) > eps:
                    if rhs_row:
                        rhs_row += " + "
                    rhs_row += f"{coefficients[row, col]} * {features[col]}"
            rhs.append(rhs_row)
        return rhs
    rhs_str = rhs(coefficients, features)  # Eager evaluation
    n_equations = len(rhs_str)
    def sindy_dynamics(t, x_, u_, params={}):
        states_inputs = x_.tolist() + np.atleast_1d(u_).tolist()
        variables_dict = dict(zip(variables, states_inputs))
        return [eval(rhs_str[i], variables_dict) for i in range(n_equations)]
    return sindy_dynamics

It returns a function in the form used by control.NonlinearIOSystem, but could be easily adapted for scipy.integrate as well. I think the most fragile part is this:

features = [f.replace("^", "**") for f in features]
features = [f.replace(" ", " * ") for f in features]

Perhaps there is a better way. For a system constructed from Lorenz system data and , I compared the performance to the built-in simulate() method and got:

SINDy built-in simulation took 44.461 s.
SINDy extracted dynamics simulation took 0.113 s.

Clearly, the performance of the extracted dynamics simulation are far superior.
The state trajectories were similar but not identical:

I don't know which one is more correct, actually. Perhaps someone can help with that.

I could do a PR, but I figured it would be best to discuss and work out a good approach first.

[Jacob-Stevens-Haas]

Thanks for your suggestion! We'd love anything that speeds up SINDy. I've got a few clarification questions

The lack of a way to extract the dynamics as an evaluable function, usable for instance, in scipy.integrate or control.NonlinearIOSystem.

SINDy.predict is a callable that takes x and u and returns the derivative of the system. With some reshaping, it can be used Is that what you mean?

For a system constructed from Lorenz system data and , I compared the performance to the built-in simulate() method and got (faster time):

Is this because of the integrator used, or the the function evaluation?

control.NonlinearIOSystem

What is this?

[ricopicone]

Thanks for your suggestion! We'd love anything that speeds up SINDy. I've got a few clarification questions

The lack of a way to extract the dynamics as an evaluable function, usable for instance, in scipy.integrate or control.NonlinearIOSystem.

SINDy.predict is a callable that takes x and u and returns the derivative of the system. With some reshaping, it can be used Is that what you mean?

The SINDy.predict method was fairly slow in my tests. I didn't dig into its code very far, but it looks like there are some validations that are checked every time it's called. I didn't benchmark it closely, but I can do so at some point, if it would help.

For a system constructed from Lorenz system data and , I compared the performance to the built-in simulate() method and got (faster time):

Is this because of the integrator used, or the the function evaluation?

I compared pySINDY's model.simulate() to the control.input_output_response() from the Python Control System Library, which just calls scipy.integrate.solve_ivp(). I think I read in the pySINDY docs that scipy.integrate.solve_ivp() is used there as well, so I think it must be the state time-derivative function (we sometimes call this the "dynamics" function) evaluation time.

control.NonlinearIOSystem

What is this?

This is a nonlinear system dynamics model class from the aforementioned Python Control Systems Library. It encodes the "dynamics" function and is used extensively by controls engineers. It's not really that important that pySINDy be able to return such a model because that can be easily constructed from the "dynamics" function.

What my "parser" function extract_sindy_dynamics() does is essentially extract the "dynamics" and writes them into the returned sindy_dynamics() function. This is akin to how we would manually write the dynamics function to use scipy.integrate.solve_ivp() or another integrator. If SINDy.predict was as fast as the resulting function, that would be ideal. Perhaps there could be a fast_eval option that skips checks for integration purposes, or something to that effect? Thanks!