Evolutionary-Based Sparse Regression for the Experimental Identification of Duffing Oscillator

In this paper, an evolutionary-based sparse regression algorithm is proposed and applied onto experimental data collected from a Duffing oscillator setup and numerical simulation data. Our purpose is to identify the Coulomb friction terms as part of the ordinary differential equation of the system. Correct identification of this nonlinear system using sparse identification is hugely dependent on selecting the correct form of nonlinearity included in the function library. Consequently, in this work, the evolutionary-based sparse identification is replacing the need for user knowledge when constructing the library in sparse identification. Constructing the library based on the data-driven evolutionary approach is an effective way to extend the space of nonlinear functions, allowing for the sparse regression to be applied on an extensive space of functions. ,e results show that the method provides an effective algorithm for the purpose of unveiling the physical nature of the Duffing oscillator. In addition, the robustness of the identification algorithm is investigated for various levels of noise in simulation. ,e proposed method has possible applications to other nonlinear dynamic systems in mechatronics, robotics, and electronics.


Introduction
e Duffing oscillator is a nonlinear dynamic system with a considerable number of engineering applications and presents a key benchmark in nonlinear system analysis. e ordinary differential equation of this system consists of a cubic nonlinear term which can result in chaotic behavior and bifurcation. Suppression strategies are required to accommodate for this behavior in flexible robotic manipulators and high-precision mechatronic systems to increase their efficacy [1]. e control performance is however drastically affected by modeling errors in the system parameters [2,3]. Also, the design of flexible manipulators and high-precision systems depend on the characteristics of the Duffing oscillator parameter variations [4]. Additionally, the design of harvesting devices from vibrations relies on the Duffing-type dynamic equations and when characterized well can be used as a tool for further analysis [5].
is work focuses on identifying the nonlinear ordinary differential equation of the Duffing system that consists of difficult to discover friction terms. Nonlinear system identification is a vast research field. e progress of this research area can be followed via several surveys including earlier works by Billings [6] and Mehra [7] as well as more recent studies [8][9][10].
In cases that the nonlinear model structure can be obtained from the first principles and is a priori known, the identification problem boils down to parameter estimation. Many works have been done in this area such as [11] where the physical parameter values are directly estimated using measured data. In many works such as [12][13][14], the least squares method is used in order to estimate the parameter values. Others report the usage of genetic programming for the same purpose [15,16].
Parameter estimation using a fixed model structure based on captured data has been previously applied on Duffing oscillator-type systems. In [17], the parameters of a numerical fractional-order Duffing system have been identified using the sequential differential evolution method. Other algorithms such as nonlinear subspace identification method, particle swarm optimization, Volterra-Wiener-based model, and Wiener-type cascade model were used to numerically estimate the parameters of Duffing-type systems [18][19][20][21]. In a more recent attempt, authors in [22] have used a tailored sequential Monte Carlo algorithm within a Markov Chain Monte Carlo (MCMC) scheme to identify the parameters of Duffing in a Bayesian manner.
Alternatively, when the model structure is not a priori known, the form of the model needs to be discovered. Different black-box model structures can be considered to form the system equations. In [23], a modeling method for nonlinear systems using polynomial nonlinear state space equations was introduced. Furthermore, NARMAX models have been used in [24,25] to represent nonlinear systems. Genetic algorithm and genetic programming have been also introduced in this field. In [26], genetic programming is used in a multiobjective fashion to generate global nonlinear models. Authors in [27] apply genetic programming to discover nonlinear differential equations. More examples of genetic algorithm application for system identification are [28][29][30]. Other modeling methods include but are not limited to neuro-fuzzy methods [31] and high-order neural network structures [32].
Black-box identification of the Duffing equation has also been a matter of investigation. In [33], explorative genetic programming is used to identify the model of a noisy Duffing-van der Pol oscillator using numerical simulation data. Artificial neural networks have been used to determine the mathematical model of the damped Duffing in [34]. A similar approach was proposed based on a set of basis functions and applying least-squares in [35].
A more exploitation-based nonlinear system identification approach was recently proposed in [36]. In this approach, a fixed matrix of candidate terms is first built upon prior expert knowledge. Subsequently, a linear system of equations is formulated using this matrix. e dominant terms in the constructed matrix later form the identified equation of the system. e sequential threshold leastsquares algorithm is applied to find the true model of the system, depending on choosing the accurate value of the regularization parameter. A revised version of this method using the alternating direction method of multipliers (ADMM) has been successfully implemented on captured data from an experimental Duffing setup [37]. e biggest criticism towards sparse identification method lies in selecting ad hoc the appropriate library functions.
is problem can be observed in [37] as the identification fails to discover the friction terms existing in the experimental data as these complex nonpolynomial terms are lacking in the library of functions. When identifying an experimental dataset, the friction forces within Duffing oscillators form an important model uncertainty that also arises in many other mechatronic applications such as in hydraulic actuators [38].
Nonlinear friction model parameters are reconstructed mostly based on a priori given friction model structures such as Coulomb and Stribeck friction models. Once these friction models are correctly identified, they can be used in control algorithms [39]. Consequently, in this paper, we aim at implementing an evolutionary-based sparse identification algorithm on the numerical and experimental Duffing system. e combination of genetic programming and sparse identification algorithm has been previously suggested in [40]; however, no methodology has been proposed so far.
In this paper, a revised version of sparse identification using the evolutionary-based sparse identification algorithm is studied for the first time to the authors' knowledge applied on a set of real-world experimental data.
is paper is organized as follows. Section 2 describes the Duffing oscillator and the collected data from the setup and simulation. Section 3 provides details on the sparse regression algorithm. Section 4 briefly introduces the genetic programming method as the base for the evolutionary construction of the library and presents the evolutionary-based sparse identification algorithm to identify the model structure and parameters of the Duffing oscillator. Results and discussions are provided in Section 5, applying the identification method on experimental and numerical Duffing oscillator data. Conclusions are drawn in Section 6.

Problem Statement and Data Acquisition
In this work, the proposed algorithm is applied on the Duffing oscillator both numerically and experimentally. e cubic Duffing equation as a differential equation with a third-power nonlinear term is an example of a dynamic system that exhibits chaotic behavior and bifurcations. Experimental datasets are extracted from this setup to identify the underlying equation. Simulations from the same setup (using its characteristic parameters) implemented in the MATLAB environment [41], furthermore, allow to examine the efficiency of the algorithm and validating the presented method. e robustness of the algorithm is investigated on the data set by increasing added noise.

Duffing Oscillator: eory and Experimental Realization
e mechanical Duffing oscillator with an imposed ground motion depicted as in Figure 1 is characterized by the following dynamic equation: with m being Duffing's mass, c its linear damping, k the linear stiffness, and k 3 the cubic stiffness. A change of coordinates to the relative ground displacement q≜x − z yields or the dynamics expressed in state space:

Design Principle of Mechanical Duffing Oscillator.
To realize the mechanical Duffing oscillator, a mass-spring system, see Figure 2, is constrained to move along a designed track y � f(x). e track's shape determines the linear and nonlinear stiffness, k and k 3 , (3). If the mass is subject to a static force in the x-direction, the mass moves along the track until equilibrium is reached. It is now shown that the spring characteristic is nonlinear. e track exerts a reaction force on the followers attached to the spring, R, perpendicular to the track's curvature. e linear spring is compressed according to the track, imposing a force on the mass in the y-direction, F y � k l y(x). e static applied force, the reaction force on the follower, and the reaction force on the mass are related by static equilibrium: with θ being the tracks curvature's angle, related to the force profile by tan θ with k � 4k l ab and k 3 � 4k l a 2 . By machining a parabolic track, f(x) � ax 2 , the linear coefficient can be simply tuned by shifting the profiles over a distance b.

Experimental Setup.
e realized mechanical Duffing oscillator with the abovementioned design principle is shown in Figure 3(a). A mass with linear springs was fitted on a linear guide rail. Tracks with the shape with a � 4 m − 1 were made from machined steel. e followers on the springs are SKF ball bearings. e linear springs have a stiffness of k l � 16.7 kN, according to the manufacturer, with the cubic stiffness then being k 3 � 1.07 MN/m 3 . e profiles can be shifted for adjusting the b term in equation (5).
To impose the ground motion, the oscillator is put on a shaking table, here a Beckhoff linear permanent magnet motor. To measure Duffing's mass and shaking table displacement, accelerometer signals are integrated with the algorithm in [42]. For this algorithm to perform well, the signals should stay in a certain frequency band. e ground displacement imposed by the shaking table is limited in bandwidth by choosing a sine sweep and a random phase multisine.
e material contact between the followers and the track causes dry friction. e force in the y-direction F y will cause a perpendicular opposing friction force, μF y , with μ the friction coefficient. e total opposing friction is with v x the speed in the x-direction. Including the friction forces in the state space representation of the Duffing oscillator dynamics is in which the viscous damping c, linear stiffness k, and dry friction coefficients μ 1 and μ 2 have to be experimentally identified.

Experimental Data.
e input of the dynamical equation (7), acceleration of the shaking table € z, and the relative acceleration between the mass and shaking table € q, captured from the described setup, are shown in Figures 4(a) and 4(b), respectively. e excitation signal of the experiment is a sine sweep from 2 to 20 Hz. e sampling time equals 0.488 ms.

Duffing
Oscillator: Numerical Data. In order to validate the performance of the algorithm on numerical data, the described Duffing setup has been simulated in MATLAB. e state space model used for the purpose of simulation is the same as (3). e control input of the system is a linear swept-frequency cosine presented in Figure 5(a). e noisy output  acceleration obtained from numerical simulation is presented in Figure 5(b). e sampling time is 0.488 ms. e amplitude is selected such that the bifurcation is observable in the output. Considering both sets of experimental and numerical data in Figures 4 and 5, it can be noted that the bifurcation occurs sooner in simulation. e bifurcation of a Duffing oscillator occurs at a certain frequency of the input sweep.
is frequency depends on the amplitude of the sweep, [43], which is different for the experiment and the simulation, explaining why the bifurcation happens at a different instant. e data for both (experimental and numerical) cases are divided in identification and validation parts.

Sparse Regression
e aim of sparse regression in the field of system identification is to extract a low-dimension (sparse) representation   of the system from a high-dimensional space of candidate representations using input and output data of the system. Considering q ∈ R n×p as the data matrix with p state variables, each presented as a column of the matrix over n time instants, sparse regression determines the state space equation as a general nonlinear function g: where u ∈ R n×1 is the input of the system, q . ∈ R n×p is the time derivative of the states which can be measured or numerically calculated, and the q matrix (with derivatives q . ) is assumed to be fully observable.
By introducing a library of terms as functions of the states and input of the system, the identification problem can be presented as finding the sparse matrix ξ ∈ R m×p [36]: where A is the library of (non)linear terms. Choosing the right form of nonlinearity in the construction of the dictionary is essential in this approach which requires user knowledge. Equation (10) illustrates such a library. Each column, m, corresponds to a linear/nonlinear term as a function of the states or the input: By solving equation (9), the dominant linear and nonlinear elements of the library A(q,u) will be chosen to combine linearly and form the equation of the system g in equation (8). e ξ matrix is determined by minimizing a defined optimization problem. In this paper, we define the optimization problem as the elastic net regulator [44]: where λ 1 and λ 2 are the hyperparameters that are changed discretely. e order of magnitude for these hyperparameters is defined through the parameter sweep.

Evolutionary-Based Sparse
Regression Methodology 4.1. Genetic Programming. Genetic programming (GP) is a subclass of genetic algorithms that was first presented by Koza in 1992 [45]. e basic idea of genetic programming is to evolve populations of equations based on the captured data and the fitness function evaluation of the simulation of each equation, where each equation is presented as a tree.
In the first generation, a population is randomly constructed by combining the numbers, variables, and mathematical operations. Terminal nodes of the trees are occupied by variables and numbers. e operators consisting of basic algebraic operations (+, −, ×, and /), functions (sin, cos, tan, abs, and sgn), or user-defined functions fill in the nonterminal nodes called the primitives. Afterwards, the population can vary in two ways: crossover and mutation. A crossover happens when two parent trees randomly exchange branches to form new offspring ( Figure 6). Mutation Mathematical Problems in Engineering involves random alteration of a parent's subtree (Figure 7). In the next step, the algorithm evaluates the fitness of each tree. e next generation is built based on the fitness evaluation. Following, the algorithm cycles through this loop until it reaches the stopping criteria or its convergence. A typical error metric such as least squares or root mean squared error is used as the fitness measure.

ESparse Algorithm.
Following the description of sparse regression and genetic programming, in this section, the proposed algorithm is described. As presented in Algorithm 1, the identification procedure consists of two main steps: (1) Construction of the library ( A n×m ) using genetic programming (2) Performing the sparse regression In each iteration, an ODE equation is realized by solving a layered optimization problem: the individual trees in the population are used as the functions to build the A n×m library in (10). Next, the sparse regression is performed on the constructed library by solving (11). Based on this approach alternative to a predefined library, the sparse regression is applied on a dynamic set of functions generated from the genetic programming. e advantage of this method is clearly its ability to generate an explorative library consisting of an extensive space of functions derived from the captured data. Moreover, this alternation between the explorative step 1 and the exploitative step 2 allows a reduction in the number of terms for regression.

Results and Discussion
In this section, the ability of identifying the correct form of the Duffing equation using the method from Section 4.2 in case of both numerical and noisy experimental datasets is analyzed. Both sets of data are captured from the Duffing oscillator described in Section 2, as a nonlinear dynamic system benchmark. In case of experimental data, we are specifically looking for the identification of the state space including the friction term as in equation (7). We also investigated the robustness of the algorithm with respect to noise in the data. By changing the level of added noise in simulation and how the accuracy of the identified model is affected by that noise provides a means to assess the robustness.  Table 1. Moreover, q, _ q, and € z are the inputs of the GP denoted as X0, X1, and X2. e theoretical ODE equation together with the identified model for different levels of signal-to-noise ratio (SNR) is given in Table 2. e associated error percentage is calculated using validation data.

Robustness Analysis.
To demonstrate the robustness of the algorithm, various levels of Gaussian white noise with zero mean were added to the data set. Figure 8 presents the tree of the identified equation in case of SNR � 19.5 dB. Moreover, the comparison between actual and identified validation data for SNR � 19.5 dB and SNR � 18.5 dB is presented in Figures 9(a) and 9(b), respectively.
A more general assessment for large ranges of signal-tonoise ratio (SNR) was performed, and the results are presented in Figure 10. Each data point in this figure corresponds to the mean of the accuracy percentage of 20 identification runs. Error bars are as well depicted that relate to the standard deviation. e results suggest that the proposed algorithm possesses the capability to reveal both the structure of the governing equation as well as the parameter values of the Duffing system. For SNR values from 20 to approximately 17 dB that correspond to increasing noise level, the accuracy of the identified parameter values decreases and ultimately the accuracy of the identification procedure itself. However, no additional terms appear in the discovered model, indicating the robustness of the presented algorithm. As can be observed in Figure 10, for low SNR (lower than approximately 17 dB) more terms are added to the equations. is clearly indicates that the data become overfitted by the identified model ultimately resulting in deteriorated accuracies.

Experimental
Duffing. Similar to the numerical Duffing data, ESparse is applied onto experimental Duffing data with the purpose to identify the Duffing equation. We conducted three experiments with the same input acceleration profile under the same conditions, resulting in data presented as in Figure 11. In all cases, the first 90000 data samples of the control input and the output are selected for validation and the rest are used for training. e evolutionary parameters of the genetic programming are given in Table 3. Table 4 summarizes the identified model obtained from the three experiments using the ESparse algorithm. Additionally, Figure 12 presents the tree of the identified equation with 5.6% error. e terms appearing in these equations are well supported by the theoretical model from equation (7) that includes Coulomb friction. e results demonstrate the ability of the algorithm to identify nonpolynomial nonlinearities. Comparison between the actual and identified output acceleration data is illustrated in Figure 13. A clear correlation between the two sets of data can be observed from Figure 13(b).

Comparison with Other Available Methods.
To substantiate the advantages of the proposed ESparse algorithm, a comparison with other available methods with respect to performance measures run time and % error are drawn in Table 5. Sparse regression (see Section 3) and genetic programming (see Section 4.1) are applied on the same dataset. For the purpose of having a fair comparison, crossover and mutation probabilities as well as the employed basis functions applied for genetic programming are the same as those employed in the ESparse algorithm (Table 3). However, to achieve the correct model of the system using genetic programming, the population size and number of generations have to increase to 250 and 80, respectively. As suggested by the results, the ESparse algorithm is capable of converging to the model with the same level of accuracy with much less computational effort. As for the sparse regression method, we had to manually include the sign function as being part of the library of (non)linear terms (coming from knowledge gained with the ESparse algorithm) since otherwise the model structure cannot be discovered, whereas the ESparse algorithm automatically builds the proper library using genetic programming. Generate new population using crossover and mutation (8) End for (9) End procedure ALGORITHM 1: ESparse algorithm. Advantages. e proposed methodology has major benefits in comparison with sparse regression and genetic programming-based methods for nonparametric identification.
e evolutionary-based sparse regression requires lower computational effort relative to genetic programmingbased algorithms. For GP-based algorithms to converge to the true solution, large populations with high number of generations are typically required. Nonetheless, the presented ESparse algorithm has the ability to converge to the correct model with less computational effort and having a balanced model complexity since ESparse alternates between exploration (genetic programming) and exploitation (sparse regression). erefore, the algorithm can discover the system equation with fewer generations and smaller populations.
As for the sparse regression method, the strict model assumptions prior to identification can limit the model complexity while the dynamic library of functions in    evolutionary-based sparse regression allows for discovery of more complex models by extending the search space and replacing the need for user knowledge for the construction of the library with data-driven GP step.

Limitations.
Although the proposed algorithm allows to identify more complex nonpolynomial terms in the equation such as friction terms, the basic building blocks are required to be included in the pool of the basic functions of the GP algorithm. Otherwise, the identified system will only be composed of available blocks which may not represent the nature of the system accurately.

Conclusion
In this paper, an evolutionary-based sparse regression algorithm for discovering both the structure and the parameter values of the system has been proposed. e methodology is used for the purpose of identifying the Duffing oscillator system using both numerical and noisy experimental data. In case of numerical Duffing, the data are    polluted with different levels of noise to study the robustness of the algorithm. Furthermore, the approach is challenged to discover governing dynamics that include nonpolynomial nonlinear Coulomb friction terms, from noisy experimental Duffing data. As shown by the percentage of the identification error, the algorithm is effective in unveiling the  e proposed method has possible applications to other nonlinear systems such as in mechatronics, robotics, and electronics.

Data Availability
All the data are included in the article. If there is a further demand for data, the author can provide depending on their availability.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.