# Identification of a Duffing Oscillator under Different Types of Excitation

Authors: E. Gandino; S. Marchesiello
Journal: Mathematical Problems in Engineering
 (2010)

---

## 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 theinverse 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.


---

## Body

## 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 theinverse 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 aMatlab 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)ẋ=Acx+Bcu,y=Cx+Du,
where the output y(t) is a q-dimensional column vector, t is time, the input u(t) is an m-dimensional column vector, and the order of the model, that is, the dimension of the state vector x(t), is n.A dynamical system withh degrees of freedom and with lumped nonlinear springs and dampers can be described by the following equation of motion:(2.2)Mz̈(t)+Cvż(t)+Kz(t)=f(t)-∑j=1pμjLnjgj(t)=f(t)+fnl(t),
where M,Cv, 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 gj(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 Lnj, 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 fnl(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 asx=[zTżT]T∈Rn×1 and to an input vector u=[f(t)T-g1(t)⋯-gp(t)]T∈Rm×1, is(2.3){żz̈}=[0h×hIh×h-M-1K-M-1Cv]{zż}+[0h×h0h×1⋯0h×1M-1M-1μ1Ln1⋯M-1μpLnp]{f(t)-g1(t)⋮-gp(t)},(2.4)y=[Ih×h0h×h]{zż}+[0h×h0h×1⋯0h×1]{f(t)-g1(t)⋮-gp(t)},
and matrices Ac∈Rn×n,Bc∈Rn×m,C∈Rl×n, and D∈Rl×m of (2.1) are consequently defined.Then the continuous model of (2.1) may be converted [11] into the following discrete state-space model:
(2.5)xr+1=Axr+Bur,yr=Cxr+Dur,
where A=eAcΔt∈Rn×n and B=(eAcΔt-I)Ac-1Bc∈Rn×m.


### 2.2. Subspace Identification

Given a deterministic-stochastic state-space model withs measurements of the input and of the output
(2.6)xr+1=Axr+Bur+wr,yr=Cxr+Dur+νr,
where wr 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 (2.7)U0|2i-1=def[u0u1⋯uj-1u1u2⋯uj⋮⋮⋱⋮u2i-1u2i⋯u2i+j-2]∈R2mi×j,
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 Y0∣2i-1∈R2li×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 matricesA,B,C, and D have been estimated by a subspace method in the time domain. In fact, system parameters (included in M,Cv,K, and μj) are contained in the matrix
(2.8)HE(ω)=D+C(iωI-Ac)-1Bc,
which is invariant under the similarity transformation corresponding to the application of a subspace method [11].


## 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)ẋ=Acx+Bcu,y=Cx+Du,
where the output y(t) is a q-dimensional column vector, t is time, the input u(t) is an m-dimensional column vector, and the order of the model, that is, the dimension of the state vector x(t), is n.A dynamical system withh degrees of freedom and with lumped nonlinear springs and dampers can be described by the following equation of motion:(2.2)Mz̈(t)+Cvż(t)+Kz(t)=f(t)-∑j=1pμjLnjgj(t)=f(t)+fnl(t),
where M,Cv, 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 gj(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 Lnj, 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 fnl(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 asx=[zTżT]T∈Rn×1 and to an input vector u=[f(t)T-g1(t)⋯-gp(t)]T∈Rm×1, is(2.3){żz̈}=[0h×hIh×h-M-1K-M-1Cv]{zż}+[0h×h0h×1⋯0h×1M-1M-1μ1Ln1⋯M-1μpLnp]{f(t)-g1(t)⋮-gp(t)},(2.4)y=[Ih×h0h×h]{zż}+[0h×h0h×1⋯0h×1]{f(t)-g1(t)⋮-gp(t)},
and matrices Ac∈Rn×n,Bc∈Rn×m,C∈Rl×n, and D∈Rl×m of (2.1) are consequently defined.Then the continuous model of (2.1) may be converted [11] into the following discrete state-space model:
(2.5)xr+1=Axr+Bur,yr=Cxr+Dur,
where A=eAcΔt∈Rn×n and B=(eAcΔt-I)Ac-1Bc∈Rn×m.


## 2.2. Subspace Identification

Given a deterministic-stochastic state-space model withs measurements of the input and of the output
(2.6)xr+1=Axr+Bur+wr,yr=Cxr+Dur+νr,
where wr 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 (2.7)U0|2i-1=def[u0u1⋯uj-1u1u2⋯uj⋮⋮⋱⋮u2i-1u2i⋯u2i+j-2]∈R2mi×j,
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 Y0∣2i-1∈R2li×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 matricesA,B,C, and D have been estimated by a subspace method in the time domain. In fact, system parameters (included in M,Cv,K, and μj) are contained in the matrix
(2.8)HE(ω)=D+C(iωI-Ac)-1Bc,
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 Figure1, whose motion is described by the following Duffing equation:
(3.1)mz̈(t)+cż(t)+kz(t)+k3z3(t)=f(t)
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 = k3,g1(t) = -z3(t), and obviously Ln1 = 1. The system is excited by two different types of force.Table 1
System parameters.
m (kg)k (N/m)c (Ns/m)k3 (N/m3)1.38001.31.5×106Figure 1
The nonlinear system described by the Duffing equation.Case 1.
The first one is a linearly varying frequency sweep (of amplitudeA=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Δt=10-3 s) of the equation of motion has been performed and a total number of s=105 samples has been generated (so tfin = 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 matrixHE(ω), defined in (2.8), can be easily computed for ω = 0: (3.2)HE(0)=D-CAc-1Bc=[00]-[10][-ck-mk10][001mk3m]=[1kk3k].
From the eigenvalues of the system matrix Ac 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:
(3.3)ωn=km,ζ=cccrit=c2km.
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 Table2: 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).
mkck3(Case1-up) Upward sweep4.634.014.045.86(Case1-down) Downward sweep1.711.302.643.97(Case2) Random0.130.540.730.73Effect of noise corruption for Case1-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)(b)A slightly better result for Case1 can be obtained by considering k3 as depending on ω: for each ω, matrix HE(ω) defined in (2.8) simply reduces to a vector h̅E with two elements as in (3.2), and it is possible to compute k3=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 k3 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).Figure 3
Real part of the estimated nonlinear coefficientk3, in the frequency range considered.In Figure4(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) 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.
(a)(b)


### 3.2. Output Prediction

The NSI method presented in this paper is also attractive for its predictive capability. In fact, once the system matricesA,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 Case2 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 tfin=1000 s, in order to have a slower frequency sweep and to obtain a more accurate representation of jump phenomena.In Figure5 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.Downward prediction. (a) Comparison between true and predicted output, near the jump-up. (b) Zoom just after the jump.
(a)(b)


## 3.1. Identification

The invariant matrixHE(ω), defined in (2.8), can be easily computed for ω = 0: (3.2)HE(0)=D-CAc-1Bc=[00]-[10][-ck-mk10][001mk3m]=[1kk3k].
From the eigenvalues of the system matrix Ac 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:
(3.3)ωn=km,ζ=cccrit=c2km.
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 Table2: 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).
mkck3(Case1-up) Upward sweep4.634.014.045.86(Case1-down) Downward sweep1.711.302.643.97(Case2) Random0.130.540.730.73Effect of noise corruption for Case1-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)(b)A slightly better result for Case1 can be obtained by considering k3 as depending on ω: for each ω, matrix HE(ω) defined in (2.8) simply reduces to a vector h̅E with two elements as in (3.2), and it is possible to compute k3=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 k3 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).Figure 3
Real part of the estimated nonlinear coefficientk3, in the frequency range considered.In Figure4(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) 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.
(a)(b)


## 3.2. Output Prediction

The NSI method presented in this paper is also attractive for its predictive capability. In fact, once the system matricesA,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 Case2 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 tfin=1000 s, in order to have a slower frequency sweep and to obtain a more accurate representation of jump phenomena.In Figure5 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.Downward prediction. (a) Comparison between true and predicted output, near the jump-up. (b) Zoom just after the jump.
(a)(b)


## 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ℋ∈Rj×2(m+l)i, constructed from all input and output measurements:
(4.1)ℋ=1j[U0∣2i-1TY0∣2i-1T]=1j[U0∣i-1TUi∣iTUi+1∣2i-1TY0∣i-1TYi∣iTYi+1∣2i-1T]=QR,
where R∈R2(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∈Rj×2(m+l)i is not needed.


### 4.1. Memory Limitations

Assuming to work in aMatlab environment, matrix R in (4.1) should easily be computed through the standard “qr” function, after constructing the block Hankel matrix ℋ∈Rj×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~104 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 sampless as possible (so j≈s should be of the order of 105 or 106) 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 ℋ 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 (increasingl) or systems having many nonlinear terms (increasing m).


### 4.2. New Algorithm

It is then necessary to conceive a new algorithm to compute the QR factorisation. This algorithm is based onMatlab 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 Figure6.(1)
Load measured datay, 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 sampless for the identification procedure and the number of block rows i; this choice determinates the number of rows and columns of matrix ℋ, respectively, j=s-2i+1 and d=2(l+m)i.(3)
Start aCycle 1, k=1,…,d; define δ as the kth column of matrix ℋ.   δ is constructed by using the input (if it is a column of submatrix U0∣2i-1T) or output (if it is a column of submatrix Y0∣2i-1T) data, as defined in (4.1).(4)
Start aCycle 2,   g=1,…,k-1; for each iteration g:
(a)
“load” from the hard disk vectorQg=[vg,…,vj]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” vectorQg 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,…,δ̅k-1]T and ξ=[δ̅k,…,δ̅j]T. Make a copy ψ of vector ξ.(6)
Apply Algorithm1 to vector ψ, which becomes the new Qk=[vk,…,vj]T obtaining also number βk.(7)
Execute, on vectorξ, the transformations defined in Algorithm 2, in order to obtain the new vector ξ=[ξ̅1,0,…,0].(8)
Attain thekth column of matrix R, denoted here as Rk:
(a)
construct vectorR̃=[γξ̅]T∈Rj;(b)
truncate vectorR̃, by eliminating all unnecessary zeros and keeping only the first d elements, in order to obtain Rk∈Rd.(9)
“save” vectorsQk and Rk on the hard disk, and “clear” them from the virtual memory.End of Cycle 1.(10)
Reconstruct matrixR, by loading (load) the d columns Rk from the hard disk.At the end of the algorithm, all saved vectorsQk and Rk (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 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 indexi 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.


### 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 Section3 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=105 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 Table2 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.Table3 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 k3 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 HE(ω=0), as k and knl 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.Table 3
Identification results (noise 1%): percentage error(100·|estimated–actual|/actual).
imkck3600.130.540.730.73900.130.330.570.491200.080.130.150.211800.070.110.330.18Table 4
Identification results (noise 3%): percentage error(100·|estimated–actual|/actual).
imkck3600.681.871.542.98900.761.371.222.321200.570.740.801.371800.510.660.631.20Table 5
Identification results (noise 5%): percentage error(100·|estimated–actual|/actual).
imkck3601.193.080.536.24901.602.620.155.291201.411.840.733.821801.261.611.593.34From Table3 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 ofi in order to attain acceptable accuracy in the estimates, in particular as regards the nonlinear coefficient k3.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 Tables4 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.


## 4.1. Memory Limitations

Assuming to work in aMatlab environment, matrix R in (4.1) should easily be computed through the standard “qr” function, after constructing the block Hankel matrix ℋ∈Rj×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~104 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 sampless as possible (so j≈s should be of the order of 105 or 106) 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 ℋ 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 (increasingl) or systems having many nonlinear terms (increasing m).


## 4.2. New Algorithm

It is then necessary to conceive a new algorithm to compute the QR factorisation. This algorithm is based onMatlab 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 Figure6.(1)
Load measured datay, 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 sampless for the identification procedure and the number of block rows i; this choice determinates the number of rows and columns of matrix ℋ, respectively, j=s-2i+1 and d=2(l+m)i.(3)
Start aCycle 1, k=1,…,d; define δ as the kth column of matrix ℋ.   δ is constructed by using the input (if it is a column of submatrix U0∣2i-1T) or output (if it is a column of submatrix Y0∣2i-1T) data, as defined in (4.1).(4)
Start aCycle 2,   g=1,…,k-1; for each iteration g:
(a)
“load” from the hard disk vectorQg=[vg,…,vj]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” vectorQg 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,…,δ̅k-1]T and ξ=[δ̅k,…,δ̅j]T. Make a copy ψ of vector ξ.(6)
Apply Algorithm1 to vector ψ, which becomes the new Qk=[vk,…,vj]T obtaining also number βk.(7)
Execute, on vectorξ, the transformations defined in Algorithm 2, in order to obtain the new vector ξ=[ξ̅1,0,…,0].(8)
Attain thekth column of matrix R, denoted here as Rk:
(a)
construct vectorR̃=[γξ̅]T∈Rj;(b)
truncate vectorR̃, by eliminating all unnecessary zeros and keeping only the first d elements, in order to obtain Rk∈Rd.(9)
“save” vectorsQk and Rk on the hard disk, and “clear” them from the virtual memory.End of Cycle 1.(10)
Reconstruct matrixR, by loading (load) the d columns Rk from the hard disk.At the end of the algorithm, all saved vectorsQk and Rk (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 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 indexi 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.


## 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 Section3 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=105 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 Table2 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.Table3 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 k3 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 HE(ω=0), as k and knl 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.Table 3
Identification results (noise 1%): percentage error(100·|estimated–actual|/actual).
imkck3600.130.540.730.73900.130.330.570.491200.080.130.150.211800.070.110.330.18Table 4
Identification results (noise 3%): percentage error(100·|estimated–actual|/actual).
imkck3600.681.871.542.98900.761.371.222.321200.570.740.801.371800.510.660.631.20Table 5
Identification results (noise 5%): percentage error(100·|estimated–actual|/actual).
imkck3601.193.080.536.24901.602.620.155.291201.411.840.733.821801.261.611.593.34From Table3 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 ofi in order to attain acceptable accuracy in the estimates, in particular as regards the nonlinear coefficient k3.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 Tables4 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.


## 5. Conclusions

In this paper the NSI method is presented and applied to an oscillator described by the Duffing equation, in order to handle theinverse 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 aMatlab environment, in those cases in which the data matrix is too large to be stored or factorised.



---

*Source: 695025-2010-05-31.xml*