Global stability using SINDy

[waleedaero]

Hello,

I'm trying SINDy for system identification of oscillating cylinder in Re100 fluid. However, When I change the polynomial library from 2 to 3 order, I get an error (ValueError: Input contains NaN, infinity or a value too large for dtype('float64')) for the given threshold value.
Moreover, the first 10 eigenvalues of the system identified using the DMD method gave eigenvalues [0.99871893 0.99871893 0.99999925 0.99999925 1.00241263 1.00241263 0.99940356 0.99940356 1.00349622 1.00349622].

Later I applied trappingSINDy as the optimizer. The new eigenvalues of the system are -3.64818603e-01 -3.23647315e-01 -5.18377506e-03 -3.83677227e-03
-1.77580647e-05 6.60492746e-04 1.14485518e-03 1.72192342e-01
1.18335711e+00 1.29927824e+00] . For constraints, I used the same methodology as in example 8.

Shouldn't all values go negative, as the original eigenvalues are close to 1?

I might be missing out some concepts and understanding. Looking forward to suggestions to make the system globally stable.

[akaptano]

Hi Waleed,

Thanks for opening a GitHub discussion on this. There are a few things to address here.

1. I am assuming that your error occurs when you call model.simulate(...)? Is this correct? If that is the case, it means SINDy has identified a model that is unstable for the initial condition you are using with model.simulate. This is quite a common issue and in general is hard to address without training on additional data, or trying the TrappingSR3 optimizer as you have done here. This issue is in fact the motivation for the Trapping SINDy method.

2. I can tell you that your DMD eigenvalues are incorrect (or perhaps your DMD eigenvalues should be the logarithm of the ones you reported here), since the oscillating cylinder is known to have modes that are roughly harmonics that go like the Strouhal number. You should check these results and verify that the DMD is giving you modes you can recognize.

3. In the trapping method, the eigenvalues are not guaranteed to be all negative, they are merely "promoted" to be so. You will need to play with the amount of data and the algorithm hyper parameters to achieve the all negative eigenvalues. Moreover, I recommend switching to the trapping_extended branch of the PySINDy code for now, where we have made some bug fixes (in particular there is a small bug for systems with size > 3 so important for your problem). Let me know if you have any issue switching over.

[waleedaero]

Hello,

Thank you for your response,

1. Yes, I get the error, When I do model.simulate
2. I was using the following method to get the Eigen values and looking at only the real part of the discrete system for stability.

D= Fields[:, :439] #Doing SVD again to get dimensions of D and Dprime accurate for DMD method
Dprime= Fields[:, 1:]
u, s, vh = np.linalg.svd(D, full_matrices=False)
sdiag= np.diag(s)
u1= u[:, :439]
u1T= u1.transpose()
vhT= vh.transpose()
sdiag_inv = np.linalg.inv(sdiag)
Atil1= np.dot(u1T,Dprime)
Atil2= np.dot(Atil1,vhT)
Atilde= np.dot(Atil2,sdiag_inv)
eigen= np.linalg.eigvals(Atilde)

3)Moreover, do you think I need to modify the constraints matrix? I'm using the same as in example 8. Von Karman vortex

[akaptano]

With regards to 3), no you shouldn't have to modify the constraints matrix unless you want to include the possibility of a constant term in the model. For that case, I will have a version of the constraints working on the trapping_extended branch of the PySINDy code by next week, which you can find on that branch in the example 17 jupyter notebook on trapping theorem extensions.

[waleedaero]

Thanks