A Local Sparse Screening Identification Algorithm with Applications

Extracting nonlinear governing equations from noisy data is a central challenge in the analysis of complicated nonlinear behaviors. Despite researchers follow the sparse identification nonlinear dynamics algorithm (SINDy) rule to restore nonlinear equations, there also exist obstacles. One is the excessive dependence on empirical parameters, which increases the difficulty of data pre-processing. Another one is the coexistence of multiple coefficient vectors, which causes the optimal solution to be drowned in multiple solutions. The third one is the composition of basic function, which is exclusively applicable to specific equations. In this article, a local sparse screening identification algorithm (LSSI) is proposed to identify nonlinear systems. First, we present the k-neighbor parameter to replace all empirical parameters in data filtering. Second, we combine the mean error screening method with the SINDy algorithm to select the optimal one from multiple solutions. Third, the time variable t is introduced to expand the scope of the SINDy algorithm. Finally, the LSSI algorithm is applied to recover a classic ODE and a bi-stable energy harvester system. The results show that the new algorithm improves the ability of noise immunity and optimal parameters identification provides a desired foundation for nonlinear analyses.


Introduction
For systems analysis, models are generally established using quantitative approaches. However, such quantitative methods are very effective for linear systems modelling not for nonlinear systems [1][2][3]. As most models are nonlinear, researchers have proposed various algorithms to recover nonlinear governing equations from time series. One of the most exciting modelling approaches is the sparse representation. Markus et al. [4] used sparse identification of nonlinear dynamics for rapid model recovery. Fahimeh et al. [5] identified nonlinear dynamical systems using the sparse recovery and dictionary learning. The continuous optimization of algorithm leads to the parameters of nonlinear system models being reconstructed using the neural network algorithm [6][7][8], genetic algorithm [9][10][11][12], and particle swarm algorithm [13][14][15][16]. Despite the promising performance of such algorithms, there are still some defects during the identification procedure. It is difficult for the genetic algorithm to solve the nonlinear constraint problems, which has a strong connection with initial population and empirical parameters [17].
Because of the random generation of the population, the results of the particle swarm algorithm deviate from the optimal solution [18]. Additionally, the neural network algorithm encounters difficulties in dealing with incomplete data. The most important aspect is that the developed strategy is often affected by the well-known overfitting problem [19].
Considering the limitations of general algorithms, Steven et al. [20] leveraged the fact that most physical systems have only a few relevant terms to define the dynamics, which made governing equations sparse in high-dimensional nonlinear function space, and proposed the sparse identification nonlinear dynamics algorithm (SINDy). The algorithm uses symbolic regression to determine the dynamics and conservation laws, and balances the complexity of the model (measured as the number of model terms) with the coherence of data. Moreover, it is a decision-making process for some empirical parameters based on the analysis and evaluation of expert knowledge in terms of dealing with noisy data. Simultaneously, as a result of the interference of nonlinear factors, there are redundant terms in the identified results, so the optimal solution cannot be automatically determined. According to Steven et al. [20], in successful examples, the choice of coordinates and initial conditions is somehow fortunate because they enable sparse representation.
Accordingly, in this article, we propose a local sparse screening identification algorithm (LSSI) that combines the local linear embedding (LLE) [21][22][23], the SINDy algorithm [24,25] and the mean error screening method (MES). The new algorithm replaces all empirical parameters and effectively avoids redundant terms in multiple identified solutions. In addition, the composition of the basic function is enhanced, and it is applied to a class of non-autonomous nonlinear systems. First, the LLE algorithm's kneighbor parameter, which evolved from the hierarchical algorithm [26], is substituted for the selection of all empirical parameters, such as regularization parameters, step size, number of iterations, etc. This reduces the dependence on external experts and improves the precision of parameter selection. Meanwhile, the LLE algorithm has an inherent advantage in terms of data dimensionality reduction and noise filtering, which accordingly enhances noise robustness and accelerates high-dimensional system recovery from scratch. Second, the MES method automatically screens every possible solution to determine the optimal solution that improves the calculation efficiency. Third, the basic function introduces the time variable t to make it more complete and expands the scope of application of the SINDy algorithm.
The rest of the paper is organized as follows. Section 2 introduces the theory of the LSSI algorithm. Some important applications and comparisons between different approaches are presented in Section 3. The algorithm is applied to the data-driven modelling process of a classic ODE and a membrane type electromagnetic vibration energy harvester (EVEH), which shows promising results with respect to system identification, even starting from a strongly nonlinear and noisy reference dataset. The conclusions are drawn in Section 4. substitute for all empirical parameters that are tough to choose, due to they strongly rely upon the analysis and evaluation of expert knowledge. The training dataset filtering process is divided into three steps.
(1) Find the neighbors for each sample point In this study, the Euclidean distance measurement is used: (2) Calculate the reconstruction weight matrix w ij with the constraint condition w ij ¼0; X j = 2NðX i Þ; j¼1; 2; Á Á Á ; K; where NðX i Þ represents the neighbor points. When X j is located in the range of NðX i Þ, w ij ¼1; otherwise, w ij ¼0. Eq. (2) is calculated using the Lagrange coefficient.
Obtain the initial variable: the measurement data x (1) Find the neighbor parameter K for each sample point Using the transition matrix Z i ¼ððX i ÀX j Þ T ðX i ÀX j ÞÞ jl 2 R KÂK ; ðj; l¼1; 2; Á Á Á ; KÞ, Eq. (2) changes to min e i ðW Þ¼ Next, the weight matrix is (3) Linear reconstruction of data using the k-neighbor parameter Substituting X i for the k-neighbor linear regression P K j¼1 w ij X j of sample point X i yields where " X is the filtered training sample, " X i ði 2 ð1; nÞÞ is a column vector of " X .

The SINDy Algorithm
In 2016, Steven et al. [20] proposed the SINDy algorithm to identify nonlinear governing equations. That is _ xðtÞ ¼f ðxðtÞÞ; In this work, we extend the sparse learning framework discussed the nonlinear non-autonomous system. Substituting " X into Eq. (7), we consider dynamics driven by a differential equation with external incentives, where the vector " X ðtÞ 2 R n denotes the state of the system at time t, _ " XðtÞ is the derivative of " XðtÞ and the function f ð " X ðtÞÞ represents the dynamic constraints that define the equations of motion of the system, such as the momentum theorem.
Most physical systems have only a few nonlinear terms in their dynamics, which makes the right-hand side f ð " X ðtÞÞ in Eq. (7) sparse in high-dimensional nonlinear function space. To search those remaining terms, we collect the time series " X ðtÞ from the system structure and measures the derivative _ " X ðtÞ from " X ðtÞ. The dataset is sampled at several times t 1 ; t 2 ; Á Á Á ; t m and the two matrices can be created in Eqs. (9) and (10), Depending on the above analysis, the sparse regression problem is where Θð " X ðtÞ;tÞ is the potential function of the columns " X ðtÞ. It can be observed through simulation or measurement data according to the given initial conditions. Since the research object is a non-autonomous system, the improved basis function is shown in Eq. (12), which commonly consists of constants, polynomials, and trigonometric functions. sinðxtÞ and cosðxtÞ denote the external incentive, where x is the excitation frequency. Each column of Θð " X ðtÞ;tÞ is a candidate function for f ð " X ðtÞÞ. The higher polynomials are denoted as " where " X p j denotes the j th nonlinearities in the state " X ðtÞ, that is The sparse regression problem is set up to determine the sparse coefficients " Ξ¼ " n 1 " n 2 . . . " n n Â Ã , as described in Eq. (11). It starts with a least-squares solution for " Ξ and then filters all coefficients that are smaller than the cut off values. That is, " Ξ becomes the minimizer of In general, the solution" Ξ of Eq. (14) includes multiple solutions, as shown in Tab. 2, particularly for complex nonlinear systems. Consequently, the following MES method is used to determine the optimal one.

The MES Method
To solve that problem, the MES method is introduced to automatically select the optimal one among the identified results, which determines the minimum mean error using Pareto front analysis, where M is the number of the solution. The principle of MES method is that the minimum mean error corresponds to the optimal solution.

Application
The LSSI algorithm reduces the reliance on the selection of empirical parameters, automatically determines the optimal solution and expands the scope of adaptation. We verify the superiority of this new algorithm by modelling a classic ODE and a bi-stable EVEH.

Recovery of a Classic ODE Based on Numerical Simulation Data
The data we obtain from the physical experiments generally contains noise. Therefore, noise contained in a dataset should be considered to simulate a real-sense environment where g is the original data, G d is noisy data, d is the disturbance value, e d represents n random values (n 2 ð0; 1Þ), and C 0 is a constant noise term.
The LSSI algorithm is expanded to consider a general model which is available for a series of nonlinear systems [28,29]. Eq. (17) is a weakly nonlinear system with e ( 1; contrarily, it is a strongly nonlinear system. Additionally, we also assume C 1 ¼x 0 , C 2 ¼b 1 , C 3 ¼a 1 , C 4 ¼a 2 , C 5 ¼b 2 and C 6 ¼F 0 .
The equation is given parameter values, as shown in Tab. 1. We set , and e ¼ 1 as the initial values to numerically obtain the training dataset from Eq. (17), which includes a sequence of time states x and derivatives _ x, where _ x is computed using cubic spline interpolation. The basic function is which determines the equation by calculating the related sparse coefficients" Ξ¼ " n 1 " n 2 . . . " n n Â Ã .
Subsequently, we test the training dataset depending on the given parameter values. Fig. 2 demonstrates the procedure for successful identification from a given simulation dataset. Remarkably, it represents our innovation computing architecture that combines data filtering (LLE), sparse regression (SINDy) and optimal solution selection (MES).

Replace Empirical Parameters Using a k-Neighbor
To demonstrate the effect of empirical parameters in filtering process, different step sizes are considered as a practical example, which directly affects the precision of governing equations. As shown in Fig. 3, contrary to our initial speculation, the larger the step size we increase, the more accurate the identification results are. Consequently, we replace empirical parameters to lessen uncertainty in data filtering.
The value of K relates to the global reconstruction error eðKÞ. Figs. 4a and 4b indicate the obtained global k-neighbor. The error tolerance is limited to obtain a local k-neighbor, which is extracted from the global k-neighbor, as shown in Fig. 4c. The principle is that the minimum residual variance leads to the optimal values, where K xopt ¼ 7 and K _ xopt ¼ 8. Subsequently, data filtering can be completed with K xopt and K _ xopt according to Eq. (6). Fig. 5 denotes the effect of the setting parameters on the identification precision of Eq. (17). In the noise level d ¼ 0:001 and nonlinear disturbance e ¼ 1 case of, the SINDy algorithm fluctuates significantly, while the LSSI algorithm remains relatively stable. It is seen that the LSSI algorithm has high precision and gets rid of the dependence on traditional empirical parameters.

Identify the Optimal Solution Using MES Method
In the SINDy algorithm, we obtain the sparse coefficients with multiple solutions. If the goal is to find the optimal solution that reliably represents the data among the large number of possibilities offered in the function, screening of" Ξ needs to be enforced and the process will be discussed below.
Multiple solutions with some redundant terms exist in Eq. (11), as shown in Tab. 2. A different solution S i leads to a different model structure of the system. However as the model rises, sparse vectors" Ξ consequently x for the k-neighbor, respectively. (c) Residual variance distribution of the optimal solutions K xopt and K _ xopt . The shaded area represents the error tolerance, where the setting range considers sample numbers and nonlinear orders. According to the local k-neighbor, we obtain each residual variance value of K xopt and K _ xopt Figure 5: Calculation error under different parameters for the SINDy and LSSI algorithms. It describes the relationship between the error and each parameter value in Eq. (17) produce mean errors. Hence, the MES method is applied to select the optimal one in the mean error span, and it provides the minimum mean error that corresponds to the optimal solution S 6 in Fig. 6.
Comparison of the two algorithms denotes the LSSI algorithm whose advantage is demonstrated by its high computational precision in Tab. 3. The proposed algorithm with such sturdy noise resistance can determine the coefficients to be within 1% of the true values. We extract from Eq. (17) the time history curve that varies in the range of 5-12, as shown in Fig. 7. The local enlargement plot of the curve inside a period indicates that the LSSI algorithm is close to the original data, which similarly proves its high precision.

Analyse Noise Level and Nonlinear Perturbation
Considering the noise level d and perturbation strength e, the standard for judging the visual quality is presented by the mean error, as shown in Fig. 8a. Clearly, the mean error increases with the growth of d and e. In Fig. 8b, when the nonlinear perturbation e ¼ 1, noise level of different orders of magnitude is applied to the training samples, that is 0.001, 0.01, 0.1. In particular, the mean error of the SINDy algorithm will reach approximately 20% at d ¼ 0:1. The training samples are seriously polluted by noise, which makes it difficult to recover governing equations in testing process. However, in the mentioned above case, the LSSI algorithm has the potential to suppress noise and keeps the error precision around 12%. Note that it is more suitable for noisy data model recovery. Tab. 4 shows the influence of different nonlinear perturbations on  The new algorithm is approximately controlled within 1%, even with the strong perturbation e ¼ 1, so it is applicable to both strongly and weakly nonlinear systems.

Recovery of a Bi-Stable EVEH Based on Experimental Data
For the experimental dataset, we consider a bi-stable harvester. The model of this component is shown in Fig. 9. The working principle is that the concentric permanent magnet driven by the membrane vibrates forwards and backwards in the cavity wall, which changes the magnetic flux of the coil windings with an iron core around the front and back end covers to generate inductive electromotive force. The distance between the magnet and core influences the equilibrium position of the system. When the iron core is adjusted within a certain range, the EVEH performs in the bi-stable oscillation stage.

Theoretical Analysis of the EVEH
Given the physical parameters in Tab. 5, a theoretical equation can be established. In detail, governing equations describe the forced lateral axisymmetric vibration of a circular membrane with a centre magnet. Assuming that the rigid magnet sustains a front-back symmetric transverse vibration, an axial force is exerted as the boundary condition. Regarding the eigenfunction and the boundary conditions, the differential eigenvalue equations can be obtained. Discretising the partial differential equations obtains the ODEs. If no air resistance or random noise effects occur in the EVEH, the equation is given by [30].
where B is the damping coefficient; D and H represent the squared term of the linear frequency and the cubic term coefficient, respectively; and F cosðtÞ is the external excitation, where and F are the frequency and acceleration.
The detailed theoretical analysis of Eq. (19) is derived from [30]. We quote the theoretical results to be compared with experimental measurements and the LSSI algorithm. The following sections describe the procedure.

Experimental Analysis of the EVEH
The layout of the testing system is shown in Fig. 10. The EVEH is implemented using a Shaker (APS 400) that drives a signal generator (Agilent 33250A) with different frequencies and amplitudes. During the experiments, a data acquisition device (B & K3039) is used to record the vibration acceleration, displacement response and voltage outputs. Based on the displacement response, we apply the LSSI algorithm to recover governing equations. The main process is shown in Fig. 11.
In the experimental test, the partial displacement signal which contains noise level is extracted from the time series curve according to the sampling frequency. By analyzing and observing the frequency, the experimental measurement produces a resonant frequency with a value that reaches approximately 3.41 Hz. Fig. 11a shows the displacement signal with burrs smoothed to obtain the noise-free time series  First, data is collected as the time series curve from EVEH. Second, the noise-free data sample at an acceleration level of 0.98 m/s 2 is obtained by filtering the displacement signal, which is extracted from sweep signals. Third, the basis function is constructed from input sample data to solve " Ξ. Finally, the MES method selects the optimal solution S 17 in Tab. 6. Active terms in S 17 are synthesised into an ODE during the data filtering process. Based on noise-free training dataset, the LSSI algorithm to recover governing equations, that is, the active terms (the optimal solution) are synthesizes into an ODE, as shown in Fig. 11b.
During the equation reconstruction process, the sparse coefficients" Ξ with multiple solutions exist in the EVEH system, as shown in Tab. 6. Subsequently, the MES method is applied to select the optimal solution. According to the principle of the minimum mean error, S 17 is the optimal one, as shown in Fig. 12. Generalising the LSSI algorithm makes it possible to recover governing equations with different order nonlinearities. The resonant frequency is obtained from the identified results in Tab. 7, which is compared with theory and experiment in the following sections.

Comparative Analysis in Theory, Experiment and LSSI Identification
For governing equations with different order nonlinearities, the identified resonant frequency fluctuates in the range of 3-3.5 Hz, as shown in Fig. 13. It demonstrates that the equation with 9th-order nonlinear components is close to the original system, which can be regarded as the research foundation for future nonlinear analysis.
The discrepancy of the initial theoretical model in Eq. (19) originates from an insufficient consideration of complicated nonlinear factors in the membrane vibration, in addition to unavoidable measurement errors. According to Williams et al. [31,32], the energy harvester may be further constrained by unwanted damping owing to undesirable effects, such as air resistance. It is note that the deviations are hardness to be effectively exhibited in theoretical modelling process. Hence, data-driven modelling has become an inevitable choice to establish and improve governing equations of nonlinear vibration systems both numerically and experimentally.  Meanwhile, any infinite-dimensional continuous system expressed as partial differential equations can be discretized and expressed as a finite-dimensional matrix equation. Such a finite-dimensional system can often be expressed as a SDOF system using orthogonal transformations. Thus, the main concerned SDOF system in this article provides a good start point to the data-driven modelling process, based on which future research on more complicated systems can be built upon.

Conclusions
The LSSI algorithm displays high precision in recovering governing equations from experimental and numerical data. For the new algorithm, the k-neighbor parameter is substituted for all empirical parameters in data filtering to reduce the reliance on 'fortunate choice'. It also applies the MES method to solve the multisolution problem in the sparse recovery procedure. The basic function has been enhanced to extend the scope of the SINDy algorithm. The innovations have shown promising results in terms of noise immunity and hidden nonlinear factor mining starting from a strongly nonlinear and noisy reference dataset.  The efficiency of the LSSI algorithm is verified in two stages. First, the identified equations are compared with the original SINDy algorithm in a classic ODE. Second, the theoretical model of an EVEH system is considered and then compensated necessary nonlinear components to the mechanical model to simulate the experimental dataset. The results show that there are promising potential applications in data-driven modelling process that arises across the physical, engineering, and biological sciences.
However, other multiscale systems or high-dimensional datasets that involve more complicated nonlinear terms may be encountered in practice. Therefore, how to ensure the physical meaning of nonlinear terms in the identification process is considered as an open question, which will be explored in future work.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.
x 1 Theoretically calculated resonance frequency