Sparse identification of the dynamics of a nonlinear multistable oscillator

Recently, data‐driven modeling approaches are getting increasingly examined regarding their applicability for nonlinear mechanical or mechatronic systems. With a high data availability and often insufficiently accurate descriptions of complex behavior of real systems using established physical models, statistical models provide promising alternatives. Alongside machine learning techniques like deep neural networks, sparse regression is increasingly used to obtain models from measurement data. With sparse regression, governing equations are estimated from a given function space so that the data are explained with as few terms as possible while maintaining a low model error. This method is implemented in a framework called sparse identification of nonlinear dynamics. This paper demonstrates the application of this method on free and forced vibrations of a two degree of freedom nonlinear rotary oscillator with two stable equilibrium positions. The setup and data acquisition as well as the application of sparse identification are described. The selected function space containing the candidate functions is essential for an accurate representation of the system at hand. Monomials form a commonly used function space because they can approximate a wide variety of nonlinear characteristics. However, in this work, it is shown that monomials alone are insufficient when dry friction appears. Therefore, to account for Coulomb friction, a sign‐function is added as a candidate to the function space of monomials up to the fifth order. Adaptions of the common optimization algorithms turned out to be necessary for the inclusion of Coulomb friction. As a result, it is found that the addition of the sign‐function for Coulomb friction increases the model quality significantly.

sometimes accused of being inaccurate and slow [1].With a simultaneously increasing data availability, data-driven methods for system identification are getting increasingly more important.Recently, sparse regression is applied as a method to obtain models by their equations of motion from data by the sparse identification of nonlinear dynamics (SINDy) method [2].
In this work, the SINDy method is applied to an academic nonlinear bistable oscillator.The aim is to understand the capabilities and limitations of the method and to assess the effect of additional physics-based candidate functions for the sparse regression.Three methods of sparse regression are compared.For the model identification and selection, the root mean square error  is used as the model error.As exemplary results, time series of free oscillations of measurements and simulations are shown.
The present work bases on the unpublished Master thesis of the first author [3], following two other Master theses, where in [4] a first version of the test rig was designed and built up, and [5], where the setup was improved and first tests of the SINDy method were performed.In [3], the setup was further improved and significant improvements of the application of the SINDy method implemented.All these student works were performed in the framework of current activities of the chair MMD at TU Berlin on applications of the SINDy and other data analysis methods, for example, on bistable oscillators [6].The setup itself follows classical ideas of easy-to-built-up rotatory oscillators with translatory (not rotatory) springs, for example, described in [7].

CONSIDERED SYSTEM
In the following, the chosen system and the corresponding experimental setup are explained in detail.The considered system is an academic oscillator consisting of two rotating disks connected via two springs, as shown in Figure 1.The two disks are set up in a horizontal plane and can rotate around their vertical center axis.In general, the springs can be connected at different points on the disk marked with black dots.Also, the springs and the distance of the center axes  12 can be varied.This way, a wide variety of different nonlinearities can be realized.In the considered setup, the springs are connected on opposite points on the disks so that a bistable system is obtained, that is, the system has two stable equilibrium positions as depicted in Figure 2. The two DoFs are described by the angles  1 and  2 of the right and left disk, respectively.An excitation torque () acts upon the right disk.Thereby, based on Newton's second law, the general form of the equation of motion reads The full setup used in this paper is shown in Figure 3.The disks are mounted with a double-row angular contact ball bearing, each.For the excitation, a third disk, actuated by a stepper motor, is mounted on top of the right disk.Due to its inertia, a torque (), proportional to the relative angular acceleration ψ() between the excitation disk and the right disk, acts upon the oscillator.
To measure the positions of the three disks, three laser displacement sensors are used to measure the distance to the red and black measuring arcs fastened on the rim of each disk (see Figure 3).The arcs have a linearly increasing radius, so that the angle of rotation of each disk can be derived from the measured distance with an affine linear transformation.

MEASUREMENT SERIES AND SIGNAL PROCESSING
For data collection, 130 single measurements were made with a harmonic excitation at 13 frequencies ranging from 1.0 to 2.2 Hz.Each measurement consisted of four phases.First, the excitation was started with a linear ramp in the amplitude, so that the jerk of the motor at the beginning is minimized.Then, for the longest part of the measurement, the forced oscillation was recorded.Depending on the excitation frequency, this oscillation was either chaotic or-following a transient phase-an interwell or an intrawell oscillation, meaning an oscillation around both or one stable equilibrium, respectively.At the end of each measurement, the excitation was stopped by a linear ramp, again, and the decaying free oscillation was recorded.This way, a wide range of solution types could be observed.Afterward, the collected data are transformed into the angles of the disks and filtered to reduce the noise.Due to vibrations of the stepper motor, a rather small cut-off frequency of 5 Hz had to be selected.To maintain the phase of the signals, a zero-phase digital filtering is used with an infinite impulse response (IIR) lowpass filter.The filtered signals are differentiated using central second-order finite differences, to obtain the angular velocities and accelerations of the disks.DATA-BASED MODELING In the following, the SINDy method [2] is used to discover a model for the oscillator.

Sparse identification of nonlinear dynamics (SINDy-method)
Based on the proceedings in [2] and [6], the following formulation of the SINDy method is used.The equation of motion from Equation ( 1) is approximated as a linear combination of candidate functions The excitation is separated from the other candidate functions, because the term is known, but is also considered as the th candidate function.Candidate functions in general can be arbitrary functions.Common choices are monomials and trigonometric functions.Since a wide variety of functions can be approximated by polynomials, monomials up to the fifth order are used.In addition, the effect of a physics-based candidate function for Coulomb friction is tested and compared to a purely polynomial model.For that, the sign-function applied to the angular velocities sgn( φ1,2 ) is used as an approximation of the torque resulting from Coulomb friction in the bearing.Since the torque on each disk is only a function of the angular velocity of that disk, in Equation ( 2) only sgn( φ1 ) and in Equation ( 3), only sgn( φ2 ) is considered, and the other coefficient is set to zero.For vectors ⃗ φ1,2 , ⃗ φ1,2 , ⃗  1,2 containing observations from times series, a linear system of equations in the form is obtained.The vectors ⃗  1,2 contain the coefficients  1, and  2, , respectively, and the matrix  is called library matrix, containing the candidate functions evaluated for the states in the state vectors.
Using sparse regression with Equations ( 2) and (3) as the regression functions, a set of nonzero coefficients of ⃗  1,2 is sought so that the observations are still explained well.There are different algorithms to perform the sparse regression in the literature.The common LASSO regression and sequential thresholded least squares (STLSQ) proposed in [2] were considered alongside a stepwise sparse regression with a backward elimination, similar to [8].However, instead of the absolute value of the coefficients, the -value is used as a statistical criterion to eliminate the least important coefficient each iteration.Following [9], the -value is calculated for the estimated coefficient ξ, and the estimated standard error   with Since the variance σ2 is constant for all coefficients and the matrix ( ⊤ ) −1 is already calculated for the regression, the simplified -value t, = ξ, is used, to enhance the performance.This way, a standardization of the data can be avoided.Since the LASSO-regression yielded significantly higher model errors than the other two algorithms, it is not examined further.This result is in accordance to the findings in [10], where it is shown, that LASSO-regression tends to identify coefficients with small instead of zero magnitude.Therefore, only STLSQ and stepwise sparse regression are compared in the following model identification.

Model identification and comparison of the regression methods
For the model identification, the data are split from 80% to 20% in training and test data with 104 and 26 measurements, respectively.A random permutation of the measurements is used to split the data.The sparse regression is only performed on the training data, and the test data are then used for model selection and comparison.
With the stepwise sparse regression, a model for each number of coefficients is found, inherently, within one iteration.The STLSQ algorithm, on the other hand had, to be performed at different thresholds to achieve different model complexities.The results of the model identification are shown in Figure 4A and B as the model error on the test data using the root-mean-square error  over the model complexity represented by the number of nonzero coefficients.For higher model complexities, the STLSQ algorithm yields similar models to the stepwise sparse regression.However, with fewer coefficients, the model error rises faster, and the candidate function gets eliminated rather early.The impact of the candidate function can already be seen in the difference between the model error of the purely polynomial models and those with an additional term for Coulomb friction, especially for the second degree of freedom in Figure 4B.

Model selection
To select a model that balances model accuracy and complexity, the difference between the  for the test data and the training data is used.While Δ is larger than zero, overfitting is likely to appear.Therefore, the zero-crossing of Δ as shown in Figure 4C and D is used as a criterion for the selection of the model including Coulomb friction.With that, the model with 18 and 19 coefficients for the first and second degree of freedom is selected.The values of the nonzero coefficients with whom the equations of motion are approximated are shown in Table 1.These equations of motion are here simply given as a result of SINDy identification with polynomial and sgn candidate functions and without any further consideration with respect to physical suitability.For comparison, a purely polynomial model is selected at a similar model error for the test data at 20 and 24 coefficients.The common coefficients between the two models tend to have similar values.The  damping coefficients differ slightly more.This is expected, since the highly nonlinear zero-crossing of the friction is no longer compensated by the linear and higher order damping.

Model comparison
To evaluate the performance of the models and their predictive capabilities, the free oscillations of the measurements from the test data are compared to the numerically integrated solutions of the selected models for the same initial conditions.Two exemplary comparisons of measurement and simulation are shown in Figures 5 and 6 for an intrawell oscillation and an oscillation with a transition from the interwell to the intrawell regime.In both cases, the model with Coulomb friction performs significantly better.The oscillation subsides only shortly after the measurement, whereas the oscillation persists longer for the other model.Also, the purely polynomial model differs already at the transition between interwell and intrawell oscillation, where the model with Coulomb friction still follows the measurement well.

CONCLUSION AND OUTLOOK
In conclusion, the developed academic nonlinear bistable oscillator was successfully modeled using the SINDy-method.
It is shown that a physics-based candidate function for Coulomb friction is beneficial for the accuracy of the model.To achieve this, a stepwise sparse regression with a backward elimination via the -value was used.With the STLSQ algorithm, on the other hand, the impact of the Coulomb friction could not be captured.Reasons for this may lay in the elimination of multiple coefficients in one iteration, and also the absolute value of the coefficients as the elimination criterion.
F I G U R E 5 Comparison of the free intrawell oscillation of a measurement with the selected models.

F I U R E 6
Comparison of the free inter-and intrawell oscillation of a measurement with the selected models.
In future work, broader sets of candidate functions (e.g., following the nonlinear terms from the equations of motion classically derived from Figure 1) can be considered.Additionally, a further improvement of the excitation mechanism in the experimental setup is considered.

A C K N O W L E D G M E N T S
This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) by WA 1427/32-1 -project number 461715882.The authors thank the audience of the corresponding presentation at GAMM conference 2023 for the interesting and stimulating discussion after the talk.
Open access funding enabled and organized by Projekt DEAL.

1
Mechanical system of the considered oscillator with 2 DoFs.

F
I G U R E 2 (A) first and (B) second stable equilibrium position of the oscillator.F I G U R E 3 Full experimental setup of the oscillator.

4
of the models resulting from the sparse regression on the test data for the (A) first and (B) second degree of freedom and their respective differences Δ to the training data in (C) and (D).
Values of the nonzero coefficients of the two selected models.
TA B L E 1