Identification of a Duffing Oscillator under Different Types of Excitation

Abstract

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.

1. 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.

2. Nonlinear Subspace Identification

2.1. 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 time-invariant system is first considered, as described by the following continuous state-space model: 𝑥̇=𝐴_𝑐𝑥+𝐵_𝑐𝑢,𝑦=𝐶𝑥+𝐷𝑢, (2.1) where the output 𝑦(𝑡) is a 𝑞-dimensional column vector, 𝑡 is time, the input 𝑢(𝑡) is an 𝑚-dimensional column vector, and the order of the model, that is, the dimension of the state vector 𝑥(𝑡), is 𝑛.

A dynamical system with ℎ degrees of freedom and with lumped nonlinear springs and dampers can be described by the following equation of motion:𝑀𝑧(̈𝑡)+𝐶_𝑣𝑧(̇𝑡)+𝐾𝑧(𝑡)=𝑓(𝑡)−^𝑝_𝑗=1𝜇_𝑗𝐿_𝑛𝑗𝑔_𝑗(𝑡)=𝑓(𝑡)+𝑓_𝑛𝑙(𝑡), (2.2) where 𝑀,𝐶_𝑣, and 𝐾 are the mass, viscous damping, and stiffness matrices, respectively, 𝑧(𝑡) is the generalised displacement vector, and 𝑓(𝑡) is the generalised force vector, both of dimension ℎ, at time 𝑡. Each of the 𝑝 nonlinear components depends on the scalar nonlinear function 𝑔_𝑗(𝑡), which specifies the class of the nonlinearity (e.g., Coulomb friction, clearance, quadratic damping, etc.), and on a scalar nonlinear coefficient 𝜇_𝑗. The vector 𝐿_𝑛𝑗, 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 𝑗th 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 𝑓(𝑡) and the internal feedback forces due to nonlinearities 𝑓_𝑛𝑙(𝑡), expressed as the sum of the 𝑝 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 𝑥=[𝑧^𝑇𝑧̇^𝑇]^𝑇∈𝐑^𝑛×1 and to an input vector 𝑢=[𝑓(𝑡)^𝑇−𝑔₁(𝑡)⋯−𝑔_𝑝(𝑡)]^𝑇∈𝐑^𝑚×1, is𝑧̇𝑧̈=0_ℎ×ℎ𝐼_ℎ×ℎ−𝑀⁻¹𝐾−𝑀⁻¹𝐶_𝑣𝑧𝑧̇+0_ℎ×ℎ0_ℎ×1⋯0_ℎ×1𝑀⁻¹𝑀⁻¹𝜇₁𝐿_𝑛1⋯𝑀⁻¹𝜇_𝑝𝐿_𝑛𝑝𝑓(𝑡)−𝑔₁⋮(𝑡)−𝑔_𝑝𝐼(𝑡),𝑦=_ℎ×ℎ0_ℎ×ℎ𝑧𝑧̇+0_ℎ×ℎ0_ℎ×1⋯0_ℎ×1𝑓(𝑡)−𝑔₁⋮(𝑡)−𝑔_𝑝(𝑡), (2.3)(2.4) and matrices 𝐴_𝑐∈𝐑^𝑛×𝑛,𝐵_𝑐∈𝐑^𝑛×𝑚,𝐶∈𝐑^𝑙×𝑛, and 𝐷∈𝐑^𝑙×𝑚 of (2.1) are consequently defined.

Then the continuous model of (2.1) may be converted [11] into the following discrete state-space model: 𝑥_𝑟+1=𝐴𝑥_𝑟+𝐵𝑢_𝑟,𝑦_𝑟=𝐶𝑥_𝑟+𝐷𝑢_𝑟, (2.5) where 𝐴=e^{𝐴𝑐Δ𝑡}∈𝐑^𝑛×𝑛 and 𝐵=(e^{𝐴𝑐Δ𝑡}−𝐼)𝐴_𝑐⁻¹𝐵_𝑐∈𝐑^𝑛×𝑚.

2.2. Subspace Identification

Given a deterministic-stochastic state-space model with 𝑠 measurements of the input and of the output 𝑥_𝑟+1=𝐴𝑥_𝑟+𝐵𝑢_𝑟+𝑤_𝑟,𝑦_𝑟=𝐶𝑥_𝑟+𝐷𝑢_𝑟+𝜈_𝑟, (2.6) where 𝑤_𝑟 and 𝜈_𝑟 are unmeasurable vector signals called process error and measurement error, respectively, the subspace identification problem consists in estimating the model order 𝑛 and the system matrices 𝐴,𝐵,𝐶, and 𝐷 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 𝑈=_0|2𝑖−1^def𝑢₀𝑢₁⋯𝑢_𝑗−1𝑢₁𝑢₂⋯𝑢_𝑗𝑢⋮⋮⋱⋮_2𝑖−1𝑢_2𝑖⋯𝑢_{2𝑖+𝑗−2}∈𝐑^{2𝑚𝑖×𝑗}, (2.7) where the number of block rows 𝑖 is a user-defined index. The number of columns 𝑗 is typically equal to 𝑠−2𝑖+1, which implies that all given data are used. The output block Hankel matrix 𝑌_{0∣2𝑖−1}∈𝐑^{2𝑙𝑖×𝑗} is defined in a similar manner by replacing 𝑢 with 𝑦 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 𝐴,𝐵,𝐶, and 𝐷 have been estimated by a subspace method in the time domain. In fact, system parameters (included in 𝑀,𝐶_𝑣,𝐾, and 𝜇_𝑗) are contained in the matrix 𝐻_𝐸(𝜔)=𝐷+𝐶(i𝜔𝐼−𝐴_𝑐)⁻¹𝐵_𝑐, (2.8) which is invariant under the similarity transformation corresponding to the application of a subspace method [11].

3. 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: 𝑚𝑧(̈𝑡)+𝑐𝑧(̇𝑡)+𝑘𝑧(𝑡)+𝑘₃𝑧³(𝑡)=𝑓(𝑡) (3.1) 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 𝜇₁ = 𝑘₃,𝑔₁(𝑡) = −𝑧³(𝑡), and obviously 𝐿_𝑛1 = 1. The system is excited by two different types of force.

Table 1 

System parameters.

Figure 1 

The nonlinear system described by the Duffing equation.

Case 1. The first one is a linearly varying frequency sweep (of amplitude 𝐴=1) between 3 and 6 Hz, applied for an upward (up) and a downward (down) frequency sweep.

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 Δ𝑡=10⁻³ s) of the equation of motion has been performed and a total number of 𝑠=10⁵ samples has been generated (so 𝑡_fin = 100 s) and then corrupted by adding a zero-mean Gaussian noise (1% of the r.m.s. value of the output).

3.1. Identification

The invariant matrix 𝐻_𝐸(𝜔), defined in (2.8), can be easily computed for 𝜔 = 0: 𝐻_𝐸(0)=𝐷−𝐶𝐴_𝑐⁻¹𝐵_𝑐=−−𝑐0010𝑘−𝑚𝑘11000𝑚𝑘₃𝑚=1𝑘𝑘₃𝑘. (3.2) From the eigenvalues of the system matrix 𝐴_𝑐 it is possible to obtain [18] estimates for the angular frequency 𝜔_𝑛 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: 𝜔_𝑛=𝑘𝑚𝑐,𝜁=𝑐_crit=𝑐2√𝑘𝑚. (3.3) It is observed here that in each of the identification procedures performed, the model order 𝑛=2 is determined by inspecting a singular value plot (with 𝑖=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).

Table 2 

Identification results: percentage error (100⋅|estimated-actual|/actual).

(a)

(b)

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).

A slightly better result for Case 1 can be obtained by considering 𝑘₃ as depending on 𝜔: for each 𝜔, matrix 𝐻_𝐸(𝜔) defined in (2.8) simply reduces to a vector ℎ_𝐸 with two elements as in (3.2), and it is possible to compute 𝑘₃=ℎ_𝐸(2)/ℎ_𝐸(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 𝑘₃ 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 𝑘, because for 𝜔>0 vector ℎ_𝐸 is not defined as in (3.2).

Figure 3 

Real part of the estimated nonlinear coefficient $k_3,$ in the frequency range considered.

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 jump-down 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.

(a)

(b)

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.

3.2. Output Prediction

The NSI method presented in this paper is also attractive for its predictive capability. In fact, once the system matrices 𝐴,𝐵,𝐶, and 𝐷 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 𝑡_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.

(a)

(b)

Figure 5 

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

4. 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 ℋ∈𝐑^{𝑗×2(𝑚+𝑙)𝑖}, constructed from all input and output measurements: 1ℋ=𝑗[𝑈^𝑇_{0∣2𝑖−1}𝑌^𝑇_{0∣2𝑖−1}1]=𝑗[𝑈^𝑇_0∣𝑖−1𝑈^𝑇_𝑖∣𝑖𝑈^𝑇_{𝑖+1∣2𝑖−1}𝑌^𝑇_0∣𝑖−1𝑌^𝑇_𝑖∣𝑖𝑌^𝑇_{𝑖+1∣2𝑖−1}]=QR, (4.1) where 𝑅∈𝐑^{2(𝑚+𝑙)𝑖×2(𝑚+𝑙)𝑖} is an upper triangular matrix; note that, as shown in [17], the computation of the orthonormal matrix 𝑄∈𝐑^{𝑗×2(𝑚+𝑙)𝑖} is not needed.

4.1. Memory Limitations

Assuming to work in a Matlab environment, matrix 𝑅 in (4.1) should easily be computed through the standard “qr” function, after constructing the block Hankel matrix ℋ∈𝐑^{𝑗×2(𝑚+𝑙)𝑖}. This procedure is certainly valid and efficient for linear systems, because an accurate identification does not require the values of 𝑖 and 𝑗 to be so large to fall into the problem described below (typically 𝑗∼10⁴ and 𝑖 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 𝑠 as possible (so 𝑗≈𝑠 should be of the order of 10⁵ or 10⁶) and in particular to extend the index 𝑖 to some hundreds, especially in presence of noisy measurements. The consequent problem consists in dealing with a matrix ℋ 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 𝑙) or systems having many nonlinear terms (increasing 𝑚).

4.2. 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 ℋ 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 𝑦, representing the 𝑙 system outputs, and the values of the external force 𝑓; compute from these data the vector 𝑢 of the 𝑚 system inputs.(2)Choose the number of samples 𝑠 for the identification procedure and the number of block rows 𝑖; this choice determinates the number of rows and columns of matrix ℋ, respectively, 𝑗=𝑠−2𝑖+1 and 𝑑=2(𝑙+𝑚)𝑖.(3)Start a Cycle 1, 𝑘=1,…,𝑑; define 𝛿 as the 𝑘th column of matrix ℋ.   𝛿 is constructed by using the input (if it is a column of submatrix 𝑈^𝑇_{0∣2𝑖−1}) or output (if it is a column of submatrix 𝑌^𝑇_{0∣2𝑖−1}) data, as defined in (4.1).(4)Start a Cycle 2,   𝑔=1,…,𝑘−1; for each iteration 𝑔:(a) “load” from the hard disk vector 𝑄_𝑔=[𝑣_𝑔,…,𝑣_𝑗]^𝑇;(b)execute, on part 𝛿=[𝛿_𝑔,…,𝛿_𝑗]^𝑇 of vector 𝛿, the transformations defined in Algorithm 2, also using number 𝛽_𝑔; vector 𝛿 is obtained;(c)“clear” vector 𝑄_𝑔 from virtual memory.

Figure 6 

Flow chart representation of the new algorithm, from step $(3)$ to step $(8).$

End of Cycle 2.(5)Subdivide vector 𝛿 into two vectors 𝛾=[𝛿₁,…,𝛿_𝑘−1]^𝑇 and 𝜉=[𝛿_𝑘,…,𝛿_𝑗]^𝑇. Make a copy 𝜓 of vector 𝜉.(6)Apply Algorithm 1 to vector 𝜓, which becomes the new 𝑄_𝑘=[𝑣_𝑘,…,𝑣_𝑗]^𝑇 obtaining also number 𝛽_𝑘.(7)Execute, on vector 𝜉, the transformations defined in Algorithm 2, in order to obtain the new vector 𝜉=[𝜉₁,0,…,0].(8)Attain the 𝑘th column of matrix 𝑅, denoted here as 𝑅_𝑘:(a)construct vector 𝑅=[𝛾𝜉]^𝑇∈𝐑^𝑗;(b)truncate vector 𝑅, by eliminating all unnecessary zeros and keeping only the first 𝑑 elements, in order to obtain 𝑅_𝑘∈𝐑^𝑑.(9)“save” vectors 𝑄_𝑘 and 𝑅_𝑘 on the hard disk, and “clear” them from the virtual memory.

End of Cycle 1.(10)Reconstruct matrix 𝑅, by loading (load) the 𝑑 columns 𝑅_𝑘 from the hard disk.

At the end of the algorithm, all saved vectors 𝑄_𝑘 and 𝑅_𝑘 (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 ℋ, and the already discussed memory problems can be avoided. It is indeed sufficient to construct and factorise a new column for each iteration 𝑘 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 𝑖 and the number of samples 𝑠 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 𝑖 and 𝑠) amount of time requested for the computation of matrix 𝑅.

4.3. 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 𝑖=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 𝑖,  Matlab goes out of memory and the NSI procedure with the standard “qr” function fails.

The same time histories (𝑠=10⁵ 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 𝑖.

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 𝑖=60 are also reported for comparison purposes. For this value of 𝑖 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 𝑅 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 𝑘 and 𝑘₃ decreases as 𝑖 increases. This trend is not so evident for the estimates of 𝑚 and 𝑐: this is due to the fact that these parameters are not directly estimated from matrix 𝐻_𝐸(𝜔=0), as 𝑘 and 𝑘_𝑛𝑙 in (3.2), but they depend on the estimates of 𝑘,𝜔_𝑛, 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.

Table 3 

Identification results (noise 1%): percentage error (100⋅|estimated-actual|/actual).

Table 4 

Identification results (noise 3%): percentage error (100⋅|estimated-actual|/actual).

Table 5 

Identification results (noise 5%): percentage error (100⋅|estimated-actual|/actual).

From Table 3 it can also be observed that a value of 𝑖=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 𝑖 in order to attain acceptable accuracy in the estimates, in particular as regards the nonlinear coefficient 𝑘₃.

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 𝑖 required in order to obtain the same level of accuracy increases as the noise percentage increases.

5. 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.

Appendix

Householder Transformations

In this appendix some concepts, exploited in Section 4.2 to conceive a new useful algorithm to compute the QR factorisation of a matrix, are presented. For a detailed overview of Householder transformations (also known as elementary reflectors), see [19]. In particular, the algorithms presented below are a revised form of those contained in [19, pages 40-41].

Given a generic vector 𝑥 different from zero, the Householder transformation𝑈=𝐼−𝛽𝑢𝑢^𝑇 (A.1) with 𝑢=𝑥+𝜎⋅𝑒₁, 𝑒₁=[1,0,…,0]^𝑇, 𝜎=±||𝑥||₂ and 𝛽=2/||𝑢||²₂ yields the following relation:𝑈𝑥=−𝜎⋅𝑒₁. (A.2) It can be observed that the couple (𝑢,𝛽), formed of 𝑛+1 real numbers, is sufficient to uniquely determine matrix 𝑈, having 𝑛² elements. Thus, given a vector 𝑥=[𝜉₁,𝜉₂,…,𝜉_𝑛]^𝑇, it is possible to write an efficient algorithm providing the quantities 𝑢 (which is overwritten to 𝑥) and 𝛽 (and also 𝜎).

Algorithm 1. We have the following:(1) 𝜂←max{|𝜉_𝑖|,𝑖=1,…,𝑛}(2) 𝜎←0(3) cycle 1: 𝑖=1,…𝑛(4)    if |𝜉_𝑖|≥𝜂√eps then 𝜎←𝜎+(𝜉_𝑖/𝜂)² (5) end of cycle 1(6) 𝜎=sgn(𝜉₁)𝜂√𝜎(7) 𝜉₁←𝜉₁+𝜎(8) 𝛽←1/(𝜎⋅𝜉₁).

Note that 𝑒𝑝𝑠 stands for the lowest possible machine number, and that this algorithm avoids possible phenomena of overflow, underflow, and numerical cancellation.

The couple (𝑢,𝛽) determined through the above algorithm is sufficient to construct products of the form𝑈𝐴=𝑈[𝑎₁,𝑎₂,…,𝑎_𝑛]=[𝑈𝑎₁,𝑈𝑎₂,…,𝑈𝑎_𝑛]. (A.3) In fact, given the two vectors 𝑢=[𝑣₁,𝑣₂,…,𝑣_𝑛]^𝑇 and 𝑎=[𝛼₁,𝛼₂,…,𝛼_𝑛]^𝑇, and the number 𝛽, the substitution of 𝑎 with vector 𝑈𝑎 can be computed in the following way.

Algorithm 2. We have the following:(1)  𝜏←𝛽∑^𝑛_𝑖=1𝑣_𝑖𝛼_𝑖(2)  𝛼_𝑖←𝛼_𝑖−𝜏⋅𝑣_𝑖,   𝑖=1,…,𝑛.

As an application of the concepts introduced above, it is possible to construct 𝑛−1 elementary reflectors 𝑈₁,𝑈₂,…,𝑈_𝑛−1 such that the new matrix 𝑈_𝑛−1⋯𝑈₂𝑈₁𝐴=𝑄^𝑇𝐴=𝑅 (A.4) is upper triangular; note the orthogonality of 𝑄, which is a product of orthogonal matrices.

As a final observation, the QR factorisation can be computed even if matrix 𝐴 is rectangular 𝑚×𝑛; in this case 𝐴=QR with 𝑄∈𝐑^𝑚×𝑚 and 𝑅∈𝐑^𝑚×𝑛 and the factorisation is attained with 𝑟=min{𝑚−1,𝑛} elementary reflectors 𝑈₁,𝑈₂,…,𝑈_𝑟.