Integrate-and-Differentiate Approach to Nonlinear System Identiﬁcation

: In this paper, we consider a problem of parametric identiﬁcation of a piece-wise linear mechanical system described by ordinary differential equations. We reconstruct the phase space of the investigated system from accelerometer data and perform parameter identiﬁcation using iteratively reweighted least squares. Two key features of our study are as follows. First, we use a differentiated governing equation containing acceleration and velocity as the main independent variables instead of the conventional governing equation in velocity and position. Second, we modify the iteratively reweighted least squares method by including an auxiliary reclassiﬁcation step into it. The application of this method allows us to improve the identiﬁcation accuracy through the elimination of classiﬁcation errors needed for parameter estimation of piece-wise linear differential equations. Simulation of the Dufﬁng-like chaotic mechanical system and experimental study of an aluminum beam with asymmetric joint show that the proposed approach is more accurate than state-of-the-art solutions.


Introduction
A relevant simulation of nonlinear systems requires well-identified computer models. Numerous approaches for nonlinear system identification have been proposed recently: differential equations [1], NARMAX models [2,3], neural networks [4], spline adaptive filters [5], deep state-space models [6], etc. Differential equations attract researchers and engineers due to the possibility of developing accurate white-box models, often having a well-established physical explanation. In many practical cases, the structure of governing differential equations is known, but it is not possible to directly measure parameters of the system such as stiffness, damping coefficients, etc. Nevertheless, it is possible to reconstruct them from data obtained during the system functioning or from a certain test [7]. Even if a single state variable can be measured, it is possible to reconstruct the entire dynamics of the system using an observer and then perform the parameter estimation [8].
An accelerometer is probably the most versatile tool for recording the dynamical response of an arbitrary mechanical system, as it is a small, light-weight device that can be mounted on almost any part of the system. The obtained acceleration time-domain series can then be differentiated or integrated to reconstruct other phase coordinates for further processing [9,10]. Additional information extracted from the signal, such as frequency, can also be used for identification purposes, providing more accurate results [11].
Even a small amount of data recorded with the accelerometer may give sufficient information about the system. For example, the paper [12] describes parametric identification of a linear electro-mechanical positioning system using a short accelerometer signal. A broad spread of accelerometer-based measurement toolkits, including digital amount of phase space variables is not observed. This method is especially applicable to the identification problems utilizing the accelerometer data series.
The outline of the paper is as follows. In Section 2 the problem of system identification via least squares is briefly discussed, the theory of the IRLS method is given, and the differentiated equations of motion are derived for systems with piece-wise linear stiffness. In Section 3 we apply the proposed technique to a mechanical system based on a piece-wise linear modification of the Duffing equations and estimate stiffness in a chaotic oscillations mode. Section 4 addresses the use of the proposed technique for estimating the nonlinear stiffness of an aluminum beam obtained on the experimental stand. Finally, Section 5 gives some conclusions and discussion.

Identification of Mechanical Systems Using Least Squares
In this section, the mathematical background of the least squares method and its modifications in application to system identification is given. Special attention is paid to a particular case of mechanical systems with piece-wise linear stiffness.

Least Squares for ODE Reconstruction
Let us consider the dynamical system described by the ordinary differential equation (ODE), where a state vector derivative is a function of time t, the state vector x, and a vector of unknown parameters h:ẋ (t) = f(t, x, h). (1) Let the system have m dimensions, and let one make n observations, after which a set of system states X = {x 1 , x 2 , . . . , x n } and corresponding derivatives Y = {ẋ 1 ,ẋ 2 , . . . ,ẋ n } related to time moments t = {t 1 , t 2 , . . . t n } can be obtained. The system identification problem has an algebraic solution if Here τ i (t, x) are arbitrary real functions of vector x entries, e.g., monomials such as x 1 x 2 x 4 3 or trigonometric functions such as sin(x 2 ), and L is their quantity. This allows representing each dimension of the function f(t, x, h) as a weighted sum of some known entries where weights are unknown. From a set X = {x 1 , x 2 , . . . , x n } a matrix of entries τ i (t, x) can be composed: Write down Y as a matrix: . . . . . . . . . . . .
and let its k-th column be denoted as . . .
From the equation of dynamics (1) the identity follows: In real world measurements, all observations are contaminated with noise z k : Nevertheless, a large number of observations allow finding the coefficients h k via minimizing the squared error: or, in another notation, This problem is referred to as ordinary least squares (OLS) and the solution of (7) is well-known: If it is required to set as many entries of h k to zero as possible (obtain a sparse solution of the identification problem) we need using 1 -regularization: where |h| = ∑ L i=1 |h i | denotes 1 norm of h and λ is the regularization parameter. This problem has no solution in algebraic operations and needs numerical optimization.
The OLS and 1 -regularized least squares are applicable if the noise is Gaussian and homoscedastic, i.e., its variance is independent of the observation. When the noise z k is heteroscedastic, weighted least squares should be used. In this method weights w i = 1/σ 2 i are multiplied by each squared error, where σ 2 i is variance of the i-th measurement: The weighted least squares can also be used with 1 -regularization. In practice, variance might not be known exactly. Moreover, the real noise may be non-Gaussian, and the observations may contain many outliers. In this case, the iteratively reweighted least squares (IRLS) method should be applied. Algorithm 1 outlines this method.
Inputs of the IRLS method are tolerance tol needed to stop the search when improvements become negligible and a regularization parameter δ needed to avoid division by zero in points where the error |y ki − E i h j | becomes zero. The parameter δ is used as the minimal absolute value of the denominator, and limits the highest value of w i with 1/δ.

Mechanical System Identification
Let a mechanical system state be described by a vector of coordinates x and a vector of velocitiesẋ. When the system is stationary and linear, the equation of dynamics is written as where M stands for a matrix of inertia, C stands for a damping matrix, K is a stiffness matrix, and f(t) is external driving force.
In many practical cases, matrices of inertia and damping are linear, but K depends on x; therefore, we should investigate the following case: A piece-wise linear F(x) can be described as where X j is a certain region in the phase space, K j is a matrix and x 0j is a certain term. Let a function µ j (x) return 1 if x ∈ X j , and 0 otherwise. The functions µ j may be obtained analytically or via support vector machine technique or logistic classification.
Then, the Equation (10) reads: Let us write this equation as the first order ODE, assuming that M is non-degenerate: Speaking of identification, the first equation in (13) is trivial, and the second can be treated as a special case of a general equation of dynamics. Suppose that the term g(t) = M −1 f(t) and the acceleration a =v are known. When we have n observations on the system state, we can use the identity for reconstructing the unknown matrices When we have only time series for {a i } and {t i } and should reconstruct two others, we perform double integration due to relation This approach may cause some problems with state space reconstruction in practical applications: 1. Integration of noisy signal introduces a trend, which is challenging to eliminate, and the more integrations are performed, the more complicated problem is.
2. Integration of noisy signal tends to smooth the signal, masking high-speed, highfrequency processes with the noise, leading to decreasing the signal-to-noise ratio of these components in the entire signal.
The latter circumstance leads to underestimation of the matrices A and, especially, B j .

Proposition 1.
For practical purposes, in the case of a noisy signal, both sides of the Equation (14) may be differentiated, and unknown matrices may be found from the obtained relation more precisely.
The following relation is underlying the Proposition 1: where j denotes jerk. As stated by the Proposition 1, the Equation (14) turns into: Derivative of µ j (x(t)) equals to: where B j denotes the border of X j . In practice, one may treat the probability of x i being on the border infinitesimally small, so (15) can be rewritten as The relation (16) is further used for estimating A and B j . The following advantages over using (14) for identification can be achieved: 1. Only single integration and single differentiation are performed for obtaining independent variables, and double integration is needed for obtaining only values of µ j (x). Therefore, numerical errors due to these operations are minimal.
2. High-frequency signal component suppression after integration and high-frequency noise amplification after differentiation are minimal.
These two advantages lead to more precise parameter estimation. Nevertheless, rigorous theoretical proof of the Proposition 1 meets sufficient difficulties since it involves several unknown parameters such as noise, parameters of motion, type of numerical method used for integration and differentiation, and parameters of numerical noise introduced by the computer. These circumstances force us to study the practical advantage of the proposed approach via numerical simulation and experiment.
Further, we utilize the abbreviation DIA (double integration approach) in relation to the Equation (14) and IDA (integrate-and-differentiate approach) in relation to the Equation (16). Algorithm 2 outlines system identification via DIA.
Algorithm 2: DIA for parametric identification of piece-wise linear mechanical system The IDA is described with Algorithm 3. It is rather similar to Algorithm 2 with another definition of the matrix E and using the modification of the IRLS: IRLS with swapping. The motivation and description of this modification are as follows.

Perform weighted least squares
3. Calculate possible values of y k for each B jk y kj ← D k A k + B jk V, j = [1.
.p] 4. Find errors between actual and possible values of y k , estimate variances:

end
In a realistic scenario, there is some significant errors in x i , so that the values of µ j (x i ) are incorrect in some cases; thus, the vectors I j are determined with errors. In order to correct this error, we include a swapping step needed to correspond the vectors v i multiplied by I j to correct stiffness matrices B j .
Suppose, it was determined initially that x i ∈ X q , but the particular value y ik is better predicted if x i ∈ X s . Then, we should swap the values in the incidence vectors: Suppose σ 2 s is the variance of prediction when x i ∈ X s . Let r 2 kqi be the squared error of estimating y ki obtained when x i ∈ X q and let r 2 ksi be the squared error of estimating y ki obtained when x i ∈ X s . To avoid overtraining of the model, we must perform swapping only if where tol is the tolerance parameter.
The IRLS with swapping is outlined in Algorithm 3 in detail.
The proposed approach can be applied to the general form of the mechanical system, which equation after differentiating of (10) reads: Nevertheless, the special case of K(x) written as (11) is very common and makes the problem formulation easier, so we focused on it directly.

Lin-Ewins Mechanical Oscillator
Chaotic vibration in mechanical constructions has gained special attention since the early 1990s, when several simple chaotic systems have been discovered, including a double pendulum [34] and an oscillator with backlash [35], based on the Duffing equation. The latter work is of particular interest since it introduces a system described by the equation of type (10). In this work, R. Lin and D. Ewins proposed a nonlinear mechanical system with backlash, which is shown in Figure 1. x 0 m c c x 0 x 0 The system consists of a mass m suspended between two symmetric springs with stiffness k 1 and dampers with damping coefficients equal to c. There are two additional springs with stiffness k 2 − k 1 with a backlash x 0 between the mass and the springs. The overall force caused by springs for the parameter values x 0 = 5 · 10 −3 m, k 2 = 40, 000 N/m and k 1 = 0 N/m is shown in Figure 2, corresponding to the case when there are only two springs with backlash. The equation of motion of this mechanical system excited by a sinusoidal force is as follows: where F(x) is a piece-wise linear function: where The direct application of Equation (20) may cause malfunction of some ODE solvers because of discontinuous functions. In order to avoid inaccuracies during numerical simulation, we used the following sigmoidal approximations of µ 1 (x) and µ 3 (x): where s = 10 4 is the parameter controlling the slope of the sigmoid near x 0 and −x 0 . The following parameters of the system were used for the study: c = 4 N · s/m, x 0 = 5 · 10 −3 m, k 2 = 40, 000 N/m, k 1 = 0 N/m, A = 100 N, m = 1 kg, ω = 40 rad/s. We use ode45 MATLAB solver with x(0) = (0, 0) , and a constant stepsize h = 10 −3 s. With these parameters, the system demonstrates a chaotic behavior, as shown in Figure 3.

Rewrite the system (18) in a form of the first-order ODE:
The second equation of the system (22) can be used for the system identification.

Identification of Nonlinear Stiffness in Lin-Ewins Oscillator
Let the accelerometer be mounted on a mass m, and its signal {a i } is sampled over time series {t i }. Let the driving force A cos(ωt) be known exactly, as well as the mass m and damping coefficient c, and the series for speed {v i } and displacement {x i } are obtained by integration. Rewrite the second equation of (22) in a form appropriate for parametric identification via least squares, corresponding to DIA: Differentiating both sides of this equation, one can obtain a form appropriate for parametric identification via least squares, corresponding to IDA: where the series for jerk {j i } can be obtained by differentiation, and We perform the comparison of identification approaches in two situations. First, we assume there is no noise in the time series for {a i }, and then, we add a Gaussian noise showing that in both cases, the IDA performs better than DIA.

Noiseless Case
To implement the identification via DIA, we twice integrate time series of {a i } for obtaining {v i } and {x i } and then perform Algorithm 2.
For identification via IDA we integrate time series of {a i } for obtaining {v i } and {x i } and differentiate it for obtaining {j i } and then perform Algorithm 3. The simulation included 2 · 10 4 points corresponding to 20 s of simulation with sampling frequency f = 1 kHz. The comparison of system identification approaches is given in Figure 5.

Case of Additive Gaussian Noise
Accelerometer typical errors include the presence of additive white Gaussian noise and zero point drift, which is especially inherent to widely spread piezoelectric accelerometers [36]. These errors can be modeled with the following equation: where a i is the sample obtained from the accelerometer, a(t i ) is the true acceleration, z i is a standard normal random variable, σ is standard deviation and a 0 is zero point drift. In our experiment, we set the following values: σ = 10 m/s 2 , a 0 = 5 m/s 2 , which is about 5% and 2.5% of the acceleration amplitude, respectively. In order to eliminate possible trends in {v i } and {x i } we split the time series into 20 fragments and integrate them independently from zero and then filter them with an ideal notch filter with band 0-2 Hz, and then eliminate linear trends. The comparison of system identification approaches is given in Figure 6. It is clearly seen that the DIA failed to find the value of k 1 with enough accuracy, while the IDA performed rather well. The overall results are given in Table 1. The relative error in Table 1 was calculated as wherek 1 is the estimated value of k 1 with its actual value k 1 = 40, 000. In all cases, k 1 was found much more accurate via IDA, and the error of estimation of k 2 was below the accepted numerical error.

Experimental Study of Nonlinear Vibration of an Aluminum Beam
Aluminum extruded profile is often used for building frames, racks, enclosures in many industrial applications. It has a large variety of shapes and allows creating sophisticated constructions using various connection elements. A majority of profile manufacturers follow Bosh Rexroth standard [37] and therefore are compatible with each other. In our study, we investigate free vibrations of a vertical 20 × 20 mm beam made of aluminum profile mounted on a 20 × 40 mm fixed aluminum profile with one L-shape inner bracket made of silumin. This joint has sufficiently different stiffness for positive and negative beam tilt direction due to different stress distribution in detail. Figure 7 shows the experimental setup and its schematics. The equation of motion of this mechanical system is as follows: where σ = 2c is a total damping coefficient and F(θ) is a piece-wise linear function: where The parameters of the setup were as follows: beam length L = 0.38 m, mass of the accelerometer with a mounting clamp m 1 = 0.05 kg, distance between top of the beam and the accelerometer h = 0.01 m, mass of the beam m 2 = 0.160 kg, from where the center of inertia is estimated as and moment of inertia is The linear damping parameter of the system can be easily found from the rate of vibration amplitude decay over time and was estimated as σ = 0.045 N · sec/rad. Stiffness identification is much more challenging.
The system (25) in a form appropriate for stiffness identification via DIA reads: Identification via IDA uses the following equation: In our study, we used IMV VP-4200 piezoelectric accelerometer connected to NI PXI system for data acquisition. The sampling frequency was 1000 Hz. We recorded 10 time series up to 1 s length, thus obtaining 5375 points of data. We used the differentiation method of order 4 based on central differences and the integration method of accuracy order 3 based on Simpson's rule for obtaining a derivative and integrals of the acceleration. The following relation was used for calculating angular acceleration from linear one: whereẍ is linear acceleration measured with the accelerometer. The results of stiffness identification are shown in Figure 8. The obtained values for k 1 and k 2 are summarized in Table 2. One can see that in comparison to IDA, DIA underestimates k 1 and overestimated k 2 , which corresponds to dynamics closer to harmonic oscillations rather than nonlinear behavior. The time domain signals for the tilt angle θ shown in Figure 9. This simulation was performed by ode45 MATLAB solver with fixed stepsize h = 10 −3 s.
From Figure 9, the advantage of IDA over DIA is obvious, especially in reconstructing third angle derivative θ (3) , where the shape is reproduced in the model obtained via IDA rather realistically, unlike in the model obtained via DIA.
One can clearly see that there is a notable advantage of IDA over DIA in Figure 10. While some specific features are not reproduced in the model, such as a slight asymmetry of a trajectory shape (compare the left column with two others), the overall shape is reconstructed much more accurate with IDA than with DIA. Comparison of power spectra of acceleration signal obtained in experiment and simulations in models identified with IDA and DIA is given in Figure 11.
From Figure 11, it is also clearly seen that while both approaches allow reproducing the first harmonic rather precisely, the other harmonics are almost absent in the signal from the model identified with DIA.

Conclusions and Discussion
In this study, we considered the problem of nonlinear mechanical system identification in the form of piece-wise linear ODE from accelerometer data using a variant of the iteratively reweighted least squares method. Two aspects of the identification process were highlighted in detail.
First, we demonstrated by simulation and experimental means that using noisy acceleration series may sufficiently decrease the accuracy of the parameter identification. Thus, we propose using a differentiated governing equation containing acceleration and velocity as independent variables, and position as an auxiliary variable needed for data classification. We call an approach based on the conventional governing equation in velocity and position a double integration approach (DIA), and the proposed approach is called the integrate-and-differentiate approach (IDA).
Second, we proposed a modification of the iteratively reweighted least squares method for using with IDA. Our modification includes an auxiliary reclassification step, and the modified method is called iteratively reweighted least squares with swapping. Applying this method allows improving identification accuracy through eliminating errors in classification needed for parameter estimation of piece-wise linear ODEs.
There is a number of practical applications of the proposed approach. Accurate mathematical models given as differential equations are needed in various scientific and engineering areas, having many advantages over the other models. Special attention is paid to them in applications of nonlinear science. For instance, electronic chaotic system design and investigation require high-precision modeling, for which parameter identification can be performed [38,39]. Digital chaotic systems are often used for data encryption, and system identification can be used for a cryptographic attack against a chaotic encryption method [40]. Recent experimental studies of chaos in mechanics confirmed theoretical predictions made several decades ago. For example, a chaotic system resembling the Lin-Ewins oscillator was investigated by R.J. Chang and Y.-C. Wang in 2020 [41]. These are only a few illustrations of why system identification, and particularly, the identification of mechanical systems from accelerometer data, is of great importance.

Conflicts of Interest:
The authors declare no conflict of interest.