Identification of a Duffing Oscillator under Different Types of Excitation

In many engineering applications the dynamics may significantly be affected by nonlinear effects, which must be accounted for in order to accurately understand and robustly model the dynamics. From a practical point of view, it is very important to solve the inverse problem related to system identification and output prediction. In this paper the recently developed Nonlinear Subspace Identification NSI method is presented and applied to an oscillator described by the Duffing equation, with different types of excitation including random forces, which are demonstrated to be very suitable for the identification process. The estimates of system parameters are excellent and, as a consequence, the behaviour of the system, including the jump phenomena, is reconstructed to a high level of fidelity. In addition, the possible memory limitations affecting the method are overcome by the development of a novel algorithm, based on a specific computation of the QR factorisation.


Introduction
In many applications nonlinear effects may affect significantly the dynamics, even when the amplitude of the motion is sufficiently small.These dynamical effects must be accounted for in order to accurately understand and robustly model the dynamics.
In general, bifurcations of equilibrium positions or periodic orbits of nonlinear systems are the source of additional nonlinear features in the dynamics 1 , which result in a qualitative change in the response and also in a substantial quantitative variation in oscillatory behaviour of the system.For example 2 , in the externally excited pendulum a relatively small amplitude periodic attractor, under the variation of a control parameter such as the frequency , may lose its stability at a saddle-node bifurcation in which the system may then start to oscillate with a relatively large amplitude.
Among essentially nonlinear dynamics caused by bifurcations 1 , such as the possibility of multiple coexisting stable equilibrium positions each with its own separate domain of attraction , this paper focuses on sudden nonlinear transitions between stable attractors jumps caused by nonlinear hysteresis phenomena.
Moving to the inverse problem of nonlinear systems, many studies have been recently conducted: in this case, system parameters are unknown and have to be estimated through an identification procedure, consisting in the development of mathematical models from input and output measurements performed on the real system.
Nonlinear system identification has been thoroughly investigated in recent years and many efforts have been spent leading to a large number of methods.An exhaustive list of the techniques elaborated to identify the behaviour of nonlinear dynamical systems is hard to write and, moreover, there is no general analysis method that can be applied to all systems in all circumstances.A comprehensive list describing the past and recent developments is given in 1 .
One of the established techniques is the Restoring Force Surface RFS method, firstly introduced by Masri and Caughey 3 : this simple procedure allows a direct identification for single-degree-of-freedom SDOF nonlinear systems.There exist in the literature several applications of RFS method to experimental systems: in a recent paper 4 , it is applied for the analysis of a nonlinear automotive damper.A similar approach is the Direct Parameter Estimation DPE method, which may be applied to multidegree-of-freedom MDOF nonlinear systems: a practical implementation of the procedure was made by Mohammad et al. 5 .
Recent methods are suitable for identification of more complex nonlinear systems, in particular MDOF systems.One of them is the Conditioned Reverse Path CRP method, developed by Richards and Singh 6, 7 : this technique is based on the construction of a hierarchy of uncorrelated response components in the frequency domain, allowing the estimation of the coefficients of the nonlinearities away from the location of the applied excitation.One of the examples of experimental application is given by Kerschen et al. 8 .
More recently, Adams and Allemang 9 proposed a frequency-domain method called Nonlinear Identification through Feedback of the Outputs NIFO , which has demonstrated 10 some advantages with respect to the CRP, mainly due to the lighter conceptual and computing effort.This method exploits the spatial information and interprets nonlinear forces as unmeasured internal feedback forces.
Starting from the basic idea of NIFO, the Nonlinear Subspace Identification NSI method has been developed by Marchesiello and Garibaldi 11 , showing a higher level of accuracy with respect to NIFO.NSI is a time-domain method which exploits the robustness and the high numerical performances of the subspace algorithms.
In this paper the NSI method is applied to a Duffing oscillator, which has been studied for many years as representative of many nonlinear systems 12 .This system can be considered in order to simply describe the sudden transitions between coexisting stable branches of solutions.For this type of system there are frequencies at which the vibration suddenly jumpsup or down, when it is excited harmonically with slowly changing frequency.
One of the main topics about the study of the Duffing oscillator consists in searching for analytical expressions of the jump frequencies and the amplitudes of vibration at these frequencies.For example, Worden 13 and Friswell and Penny 14 computed these points by using the harmonic balance method, while Malatkar and Nayfeh 15 determined the minimum excitation force required for the jump phenomenon to appear, by using a method based on the elimination theory of polynomials.A recent paper by Brennan et al. 16 provides a full set of expressions determined by using the harmonic balance approach, as a link between the earlier analytical work and the later numerical studies.
In this paper the NSI estimates of system parameters are excellent and, as a consequence, the behaviour of the Duffing oscillator, including the jump phenomena, is reconstructed to a high level of fidelity.
In addition, the NSI method is enforced by the development of a new algorithm to compute the QR factorisation in a Matlab environment, in those cases in which the data matrix is too large to be stored or factorised.This new algorithm, which exploits some useful features of the Householder transformations, allows the NSI method to reach more accurate results in the parameter estimation.

Nonlinear Model
The adopted mathematical approach follows the one used in 11 , in order to derive a mathematical model for a nonlinear dynamical system.The expression for a linear timeinvariant system is first considered, as described by the following continuous state-space model: where the output y t is a q-dimensional column vector, t is time, the input u t is an mdimensional column vector, and the order of the model, that is, the dimension of the state vector x t , is n.
A dynamical system with h degrees of freedom and with lumped nonlinear springs and dampers can be described by the following equation of motion: where M, C v , and K are the mass, viscous damping, and stiffness matrices, respectively, z t is the generalised displacement vector, and f t is the generalised force vector, both of dimension h, at time t.Each of the p nonlinear components depends on the scalar nonlinear function g j t , which specifies the class of the nonlinearity e.g., Coulomb friction, clearance, quadratic damping, etc. , and on a scalar nonlinear coefficient μ j .The vector L nj , whose entries may assume the values 1, −1, or 0, is related to the location of the nonlinear element: it specifies the degrees-of-freedom joint by the jth nonlinear component and the sign of the term appearing in the equation of motion 2.2 .Written as in 2.2 , the original system may be viewed as subjected to the external forces f t and the internal feedback forces due to nonlinearities f nl t , expressed as the sum of the p nonlinear components.This concept, already used in 9 to derive the NIFO frequency-domain method, is also on the basis of the present time-domain identification method.
Assuming that the measurements concern displacements only, the state-space formulation of the equation of motion, corresponding to a state vector chosen as x z T żT T ∈ R n×1 and to an input vector u Then the continuous model of 2.1 may be converted 11 into the following discrete state-space model: where

Subspace Identification
Given a deterministic-stochastic state-space model with s measurements of the input and of the output where w r and ν r are unmeasurable vector signals called process error and measurement error, respectively, the subspace identification problem consists in estimating the model order n and the system matrices A, B, C, and D up to within a similarity transformation, which does not affect the parameter estimation.
In the "data-driven approach" 17 the input data are gathered in a block Hankel matrix The nonlinear system described by the Duffing equation.
where the number of block rows i is a user-defined index.The number of columns j is typically equal to s − 2i 1, which implies that all given data are used.The output block Hankel matrix Y 0|2i−1 ∈ R 2li×j is defined in a similar manner by replacing u with y in 2.3 .Subspace methods take advantage of robust numerical techniques such as QR factorisation and Singular Value Decomposition SVD by using geometric tools such as the oblique projections of the row space of matrices.For a complete description of the estimating procedure see 17 .
The nonlinear identification procedure is based on the computation of system parameters, once the state-space matrices A, B, C, and D have been estimated by a subspace method in the time domain.In fact, system parameters included in M, C v , K, and μ j are contained in the matrix which is invariant under the similarity transformation corresponding to the application of a subspace method 11 .

Application: The Duffing Equation
Consider the SDOF system with cubic hardening stiffness depicted in Figure 1, whose motion is described by the following Duffing equation: with system parameters summarized in Table 1.The strength, the type, and the location of the nonlinearity are defined respectively by the three scalar quantities μ 1 k 3 , g 1 t −z 3 t , and obviously L n1 1.The system is excited by two different types of force.Case 2. The second one is a zero-mean Gaussian random input whose r.m.s. is 20 N, selected so that the r.m.s. of the nonlinear force is equal to 67% of the corresponding linear stiffness force.
A fourth-order Runge-Kutta numerical integration with a time step Δt 10 −3 s of the equation of motion has been performed and a total number of s 10 5 samples has been generated so t fin 100 s and then corrupted by adding a zero-mean Gaussian noise 1% of the r.m.s.value of the output .

Identification
The invariant matrix H E ω , defined in 2.8 , can be easily computed for ω 0: From the eigenvalues of the system matrix A c it is possible to obtain 18 estimates for the angular frequency ω n of the undamped system and for the damping factor ζ, so that all system parameters can be estimated from 3.2 and from the following relationships: It is observed here that in each of the identification procedures performed, the model order n 2 is determined by inspecting a singular value plot with i 60 block rows , as shown in 11 .
The identification results for all system parameters are presented in Table 2: the best estimates are obtained by applying a random input.In fact, for Case 1, it should be observed that the added noise is related to the r.m.s. of the entire time history, which is nonstationary; so, samples corresponding to small displacements are more deeply corrupted by noise and are consequently counterproductive for the identification procedure.This is shown in Figure 2 for Case 1-up, in which this concept is more evident because the system reaches higher values of response amplitudes and then a higher r.m.s. of the time histories .
A slightly better result for Case 1 can be obtained by considering k 3 as depending on ω: for each ω, matrix H E ω defined in 2.8 simply reduces to a vector h E with two elements    as in 3.2 , and it is possible to compute k 3 h E 2 /h E 1 .The estimated coefficient of the nonlinear term is frequency dependent and complex, albeit its imaginary part is some orders of magnitude smaller than the real part.A single value can be obtained by performing a spectral mean in the frequency range from 3 to 6 Hz Figure 3 .In this way, the percentage errors related to the k 3 estimates become 2.74 for Case 1-up and 1.78 for Case 1-down.Note that this procedure is not applicable to get a spectral mean for k, because for ω > 0 vector h E is not defined as in 3.2 .In Figure 4 a the true Frequency Response Functions FRFs of the nonlinear and underlying linear system are shown in comparison with the NSI estimates, computed from the identified system parameters in Case 2. As a consequence of the results reported in Table 2, the curves are almost overlaid: an excellent agreement can be observed, even in estimating the jump-up and jump-down frequencies and responses.The values for the jumpdown and the jump-up Figure 4 b have been obtained from the approximate expressions derived in 16 : the approximation of the true jump is obtained with the real system parameters of Table 1 while the approximation of the estimated jump is obtained with the NSI estimates of Case 2.

Output Prediction
The NSI method presented in this paper is also attractive for its predictive capability.In fact, once the system matrices A, B, C, and D in 2.5 have been estimated, it is possible to predict the system behaviour when it is subject to a different type of excitation.
It is important to remark that recent methods such as CRP 6, 7 and NIFO 9 would require a second step to perform output prediction in a general case of MDOF systems.In fact, these methods only produce estimates of the underlying linear FRFs and of nonlinear coefficients.On the contrary, the NSI capability of predicting the output is intrinsic in its formulation, since a state-space model is used.In other words, system parameter estimation is not strictly necessary and this represents a great advantage of NSI in case of MDOF systems.However, for simplicity's sake, in this paper an SDOF numerical example is considered, so estimating system parameters out of state-space matrices is both possible and easy to perform.
Starting from the best estimates of system parameters, obtained through Case 2 identification procedure, it is possible to generate new time histories considering the system as excited by the frequency sweeps described in Case 1.Now the numerical integration has been performed for t fin 1000 s, in order to have a slower frequency sweep and to obtain a more accurate representation of jump phenomena.
In Figure 5 the results are shown, in terms of a comparison between the true i.e., system parameters as in Table 1 and the predicted i.e., identified system parameters time histories, for Case 1-down.In Figure 5 a it can be observed that the predicted jump-up occurs at a higher frequency at a lower time instant in the downward sweep , as expected from the FRFs zoom shown in Figure 4 b .After the jump-up, this slight shift has no longer effect on the prediction: as shown in Figure 5 b , the true and the predicted output are almost overlaid just a few seconds after the jump.Notice the high global level of accuracy of the prediction results, albeit system parameters have been estimated starting from a time history corrupted by measurement noise.

QR Factorisation
A common feature in the implementation of all algorithms concerning the subspace methods is the following QR factorisation of a block Hankel matrix H ∈ R j×2 m l i , constructed from all input and output measurements: where R ∈ R 2 m l i×2 m l i is an upper triangular matrix; note that, as shown in 17 , the computation of the orthonormal matrix Q ∈ R j×2 m l i is not needed.

Mathematical Problems in Engineering
Hard disk Step 6 Step 7 Step 3 Step 4 Step 5 Step 8

Memory Limitations
Assuming to work in a Matlab environment, matrix R in 4.1 should easily be computed through the standard "qr" function, after constructing the block Hankel matrix H ∈ R j×2 m l i .This procedure is certainly valid and efficient for linear systems, because an accurate identification does not require the values of i and j to be so large to fall into the problem described below typically j ∼ 10 4 and i do not exceed some tens .
In order to apply subspace methods to nonlinear systems with satisfactory results, it is necessary to consider as many samples s as possible so j ≈ s should be of the order of 10 5 or 10 6 and in particular to extend the index i to some hundreds, especially in presence of noisy measurements.The consequent problem consists in dealing with a matrix H which results in being too large to be stored nor factorised.
Therefore, it is clear that the NSI method undergoes severe limitations in its applicability, in particular as regards MDOF systems increasing l or systems having many nonlinear terms increasing m .

New Algorithm
It is then necessary to conceive a new algorithm to compute the QR factorisation.This algorithm is based on Matlab commands "save" and "load", which allow to save and load variables directly from the hard disk, and the command "clear", useful to clean virtual memory.
Moreover, it is observed that the development of this new procedure exploits the particular structure of the matrix H to be factorised and the useful features of Householder transformations: in particular, from now on, Algorithms 1 and 2 reported in the appendix will be considered.
The new algorithm is described in the following and a flow chart representation is given in Figure 6.
1 Load measured data y, representing the l system outputs, and the values of the external force f; compute from these data the vector u of the m system inputs.
2 Choose the number of samples s for the identification procedure and the number of block rows i; this choice determinates the number of rows and columns of matrix H, respectively, j s − 2i 1 and d 2 l m i.
3 Start a Cycle 1, k 1, . . ., d; define δ as the kth column of matrix H. δ is constructed by using the input if it is a column of submatrix U T 0|2i−1 or output if it is a column of submatrix Y T 0|2i−1 data, as defined in 4.1 .4 Start a Cycle 2, g 1, . . ., k − 1; for each iteration g: a "load" from the hard disk vector Q g v g , . . ., v j T ; b execute, on part δ δ g , . . ., δ j T of vector δ, the transformations defined in Algorithm 2, also using number β g ; vector δ is obtained; c "clear" vector Q g from virtual memory.
6 Apply Algorithm 1 to vector ψ, which becomes the new Q k v k , . . ., v j T obtaining also number β k .
8 Attain the kth column of matrix R, denoted here as R k : b truncate vector R, by eliminating all unnecessary zeros and keeping only the first d elements, in order to obtain R k ∈ R d .
9 "save" vectors Q k and R k on the hard disk, and "clear" them from the virtual memory.
End of Cycle 1.
10 Reconstruct matrix R, by loading load the d columns R k from the hard disk.
At the end of the algorithm, all saved vectors Q k and R k and β also will be deleted from the hard disk.
Note referring in particular to step 3 of the above algorithm that in this way it is not necessary to store the entire matrix H, and the already discussed memory problems can be avoided.It is indeed sufficient to construct and factorise a new column for each iteration k of Cycle 1.
As a final consideration, it should be observed that this new algorithm does not present any limitations about the choice of index i and the number of samples s to be considered in the NSI procedure.The only limitation may be represented by a larger depending on the system considered and on the choice of i and s amount of time requested for the computation of matrix R.

Application
In order to test the new algorithm and to analyse the results of the NSI procedure exploiting it, the numerical application described in Section 3 is considered.Note that the previously adopted i 60 is the maximum index for the calculator used for the computations which allows to avoid the memory limitation problems described in Section 4.1.In fact, for larger values of i, Matlab goes out of memory and the NSI procedure with the standard "qr" function fails.
The same time histories s 10 5 samples as in Section 3 are considered, and the NSI procedure with the novel algorithm is performed for higher values of the number of block rows i.
Since Table 2 shows that the best parameter estimations are obtained in Case 2 Gaussian random input , the results presented in this section refer only to Case 2. Note also that in all the following tables the results obtained by choosing i 60 are also reported for comparison purposes.For this value of i the results are the same as in Table 2, as expected: the novel algorithm does not alter the NSI results, it just proposes a useful way to compute matrix R in those cases in which Matlab produces an "out of memory" message.However it is observed that, when the standard Matlab "qr" function is still applicable, the novel algorithm is about 26 times slower because of its many savings and loadings from the hard disk.
Table 3 shows the identification results relative to an output corrupted by 1% of noise: it is clear that the percentage error in the estimates of k and k 3 decreases as i increases.This trend is not so evident for the estimates of m and c: this is due to the fact that these parameters are not directly estimated from matrix H E ω 0 , as k and k nl in 3.2 , but they depend on the estimates of k, ω n , and ζ through the relationships of 3.3 ; this may cause a sort of error propagation or compensation.This remark is also valid for Tables 4 and 5.
From Table 3 it can also be observed that a value of i 60 is anyway sufficient to obtain an excellent level of accuracy in the estimates, so the application of the new algorithm is not necessary.
The new algorithm appears to be more appealing when the output is corrupted by a higher level of noise: in this case it is necessary to increase the value of i in order to attain acceptable accuracy in the estimates, in particular as regards the nonlinear coefficient k 3 .
For this reason, the previously generated output is corrupted by adding a higher percentage of zero-mean Gaussian random noise, and the results of the identification procedures are shown in Tables 4 and 5 for 3% and 5% noise, respectively.It can be observed that the index i required in order to obtain the same level of accuracy increases as the noise percentage increases.

Conclusions
In this paper the NSI method is presented and applied to an oscillator described by the Duffing equation, in order to handle the inverse problem related to identification and output prediction.
It is shown that the best results in parameter estimation are obtained when the system is excited by a Gaussian random input, in particular in presence of a measurement noise.However, the NSI method is also applicable in case of a linearly varying frequency sweep: with this type of excitation jump phenomena are highlighted, but a reduced level of accuracy is attained.
The best parameter estimates are then exploited in order to predict the system behaviour when it is subject to a frequency sweep excitation: the output reconstruction is excellent, in particular as regards the amplitudes and the frequencies at which jump phenomena occur.
The predictive accuracy depends on the quality of parameter estimates, but their improving implies the need of processing a larger amount of data.To this purpose, the NSI method is enforced by the development of a new algorithm to compute the QR factorisation in a Matlab environment, in those cases in which the data matrix is too large to be stored or factorised.

Mathematical Problems in Engineering
to write an efficient algorithm providing the quantities u which is overwritten to x and β and also σ .Note that eps stands for the lowest possible machine number, and that this algorithm avoids possible phenomena of overflow, underflow, and numerical cancellation.
The couple u, β determined through the above algorithm is sufficient to construct products of the form UA U a 1 , a 2 , . . ., a n Ua 1 , Ua 2 , . . ., Ua n .A.3 In fact, given the two vectors u v 1 , v 2 , . . ., v n T and a α 1 , α 2 , . . ., α n T , and the number β, the substitution of a with vector Ua can be computed in the following way.
Algorithm 2. We have the following: As an application of the concepts introduced above, it is possible to construct n − 1 elementary reflectors U 1 , U 2 , . . ., U n−1 such that the new matrix is upper triangular; note the orthogonality of Q, which is a product of orthogonal matrices.As a final observation, the QR factorisation can be computed even if matrix A is rectangular m × n; in this case A QR with Q ∈ R m×m and R ∈ R m×n and the factorisation is attained with r min{m − 1, n} elementary reflectors U 1 , U 2 , . . ., U r .

Figure 2 :
Figure 2: Effect of noise corruption for Case 1-up.The r.m.s. of the entire time history is 0.0088 m. a Zoom just before the jump-down large amplitudes .b Zoom after the jump small amplitudes .

Figure 3 :
Figure 3: Real part of the estimated nonlinear coefficient k 3 , in the frequency range considered.

Figure 4 :
Figure 4: a Frequency response curves.The crosses and the circles denote the responses at the jump-up and jump-down frequencies, respectively.The dashed lines denote unstable solutions.b Zoom near the jump-up.

Figure 5 :
Figure 5: Downward prediction.a Comparison between true and predicted output, near the jump-up.b Zoom just after the jump.

Figure 6 :
Figure 6: Flow chart representation of the new algorithm, from step 3 to step 8 .

Table 2 :
Identification results: percentage error 100 • |estimated-actual|/actual .The first one is a linearly varying frequency sweep of amplitude A 1 between 3 and 6 Hz, applied for an upward up and a downward down frequency sweep.