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

Unsure on the correct use of pySINDy #505

Closed oskarmue closed 1 month ago

oskarmue commented 1 month ago

Hello

as part of a project, I would like to compare different modelling / model identification options. The system to be analysed is a proportional solenoid. More precisely, I would like to derive the magnetic force F_m of the solenoid as a function of the current i, the voltage u and the position x of the actuator. The magnetic force should then be simulated in an OpenLoop process. The current, voltage and position are known in each simulation step. For this purpose, we have recorded a large data set that varies in the position x of the actuator, the current strength i and the rate of current change. Except for the Position x, the recorded data is low noise, here are some snippets: image

Ideally, we would use SINDy to find a sparse model for F_m = f(u, i, x) Up to now I tried to build a Library with polynomial terms up to the 3. degree, sin and cos and the e-function for possible satturation effects.

the important bits of the script are as follows:

features_train = np.concatenate((data_F_1, data_i_1), axis = 1)
control_train = np.concatenate((data_u_1, data_position_1, data_F_2, data_i_2, data_u_2, data_position_2), axis = 1)
feature_names = ['F_1', 'i_1', 'u_1', 'x_1', 'F_2', 'i_2', 'u_2', 'x_2']

# Polynome
poly_lib_1 = ps.PolynomialLibrary(
    degree=3,
    include_interaction =True,
    interaction_only=False,
    include_bias=False
)

poly_lib_2 = ps.PolynomialLibrary(
    degree=3,
    include_interaction =True,
    interaction_only=False,
    include_bias=False
)

#sin, cos Terme
four_lib_1 = ps.FourierLibrary(
    n_frequencies=2
)
four_lib_2 = ps.FourierLibrary(
    n_frequencies=2
)

#weitere Terme
library_functions = [
    lambda x: np.exp(x),
    #lambda x, y: x**y
]

library_function_names = [
    lambda x: 'exp(' + x + ')',
    #lambda x, y: x + '^(' + y + ')'
]

pde_lib_1 = ps.PDELibrary(library_functions=library_functions, 
                        function_names=library_function_names,  
                        include_bias=False, 
                        is_uniform=True,
                        derivative_order=0, 
                        temporal_grid=time_train,
                        implicit_terms =True,
                        interaction_only=True, #wenndas True, dann müssen Polynome etc selber nochmal in die Library als "feature" eingefügt werden
                        differentiation_method =ps.FiniteDifference
)

pde_lib_2 = ps.PDELibrary(library_functions=library_functions, 
                        function_names=library_function_names,  
                        include_bias=False, 
                        is_uniform=True,
                        derivative_order=0, 
                        temporal_grid=time_train,
                        implicit_terms =True,
                        interaction_only=True, #wenndas True, dann müssen Polynome etc selber nochmal in die Library als "feature" eingefügt werden
                        differentiation_method =ps.FiniteDifference
)

inputs_per_library = np.array([
[0, 1, 2, 3],
[4, 5, 6, 7],
[0, 1, 2, 3],
[4, 5, 6, 7],
[0, 1, 2, 3],
[4, 5, 6, 7]
])

lib_generalized = ps.GeneralizedLibrary(
    [poly_lib_1, poly_lib_2, four_lib_1, four_lib_2, pde_lib_1, pde_lib_1],
    inputs_per_library = inputs_per_library
)

#this are alls the library terms which include the Force F_m
indices_feat = [i for i, x in enumerate(lib_feat_namen) if 'x0' in x]

print('SINDyPI Optimizer')

sindy_opt = ps.SINDyPI(threshold=0.1,
                       max_iter=10000, 
                       tol=1e-3,
                       thresholder='l1', 
                       normalize_columns=True,
                       model_subset=indices_feat,
                       verbose_cvxpy = True
)
model = ps.SINDy(
    feature_library=lib_generalized, 
    optimizer=sindy_opt, 
    feature_names=feature_names,
    discrete_time = True,
    t_default = time_train[1] - time_train[0]
)
model.fit(
    features_train, 
    u = control_train, 
    t=time_train[0:19995],
)
model.print()

x = features_train: that includes F_m_1 and i_1 from the current timestep 1 while u includes u = u_1, x_1, F_2, i_2, u_2, x_2

the variables with _2 are from timesetpts k-2. I created two libraries of each, beacause I didnt want to mix up measurements from two different time steps u and x are control inputs and in some issue, I read that data we dont want to include in the left hand size of the discovered equation could be included as u as well, thats why F_2 and i_2 are considered as u.

I chose the SINDy-PI optimizer because I expect rational terms in this dynamic.

Trying to discover only a model for F_m, i dont follow this typical state space formulation and thus I am unsure if this is a viable approach for SINDy. So far I have had no success in discovering an accurate model, so I wanted to ask for some tips/obvios mistakes I am making.