Subspace-based identification of nonlinear structures

Abstract. Conventional linear estimators give results contaminated in presence of nonlinearities and the extraction of underlying linear system properties is thus difficult. To overcome this problem, the implementation of a recently developed method, called Nonlinear Subspace Identification (NSI), is considered in this paper, by using the perspective of nonlinearities as unmeasured internal feedback forces. Although its formulation is very simple, particular care has to be taken to reduce the ill-conditioning of the problem, in order to find numerically stable solutions. To this purpose, the robustness and the high numerical performances of the subspace algorithms are successfully exploited, as shown by the implementation of the proposed method on simulated multidegree-of-freedom systems with typical nonlinear characteristics as well as on an experimental case. These examples demonstrate that the application of subspace algorithms to nonlinear system identification gives better conditioning and computational efficiency with respect to the most recent nonlinear techniques. Moreover, the capability of the NSI method of simultaneously dealing with several nonlinear terms, with a light computational effort, may be also exploited in those situations where no a priori knowledge of the location and the type of nonlinearities is given, being this method well capable of detecting the contribution of the dominant nonlinearities. On the basis of the results discussed in this paper, and compared with those of other well-assessed nonlinear techniques, the proposed method appears having the chances to become a robust procedure to be widely exploited in many industrial fields, being its capability of separating linear and nonlinear contribution terms widely requested in mechanical and civil engineering field.


Introduction
As demonstrated by recent papers [1,2], an increasing attention has been directed towards nonlinear systems in recent years and many techniques have been developed, seeking to indicate the presence of nonlinearity and to estimate the amount of nonlinear contributions, for a better understanding of nonlinear systems response.In particular, much of this research has been developed on simple single-degree-of-freedom (SDOF) systems with relative simple nonlinearities.Recently techniques have focused the attention on multi-degree-of-freedom (MDOF) systems with multiple nonlinearities.A very synthetic list should include the coherence inspection [3], the causality check via the Hilbert transform [4], the Wiener and Volterra series approaches [5], the time series models such as the NARMAX [6], the force-state mapping technique [7].Recently, the "Reverse Path" formulation (RP) [8] was proposed with its "conditioned" version (CRP), based on the construction of a hierarchy of uncorrelated response components in the frequency domain.More recently, Adams and Allemang [9] proposed a frequency domain method called Nonlinear Identification through Feedback of the Outputs (NIFO), which has demonstrated some advantages with respect to the CRP, mainly due to the lighter conceptual and computing effort.The basic idea of NIFO consists in the interpretation of nonlinearities as unmeasured internal forces, i.e. considering nonlinearities as cause of distortion in the linear Frequency Response Function matrix.Starting from this idea, the authors have developed [10] an efficient time domain method for identifying nonlinear vibrating structures, which exploits the robustness and the high numerical performances of the subspace algorithms.

Nonlinear subspace identification
The equation of motion of a dynamical discrete system with h degrees of freedom, carrying lumped nonlinear springs and dampers can be expressed as: 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) the generalised force vector, both of dimension h, at time t.The nonlinear term is expressed as the sum of p components, each of them depending on the scalar nonlinear function g j (t), indicating the class of the nonlinearity, through a vector L nj , which indicates the location of the nonlinear component and whose elements may assume the values 1, −1 or 0. By moving the nonlinear term of Eq. ( 1) to the right hand side the original system may be viewed as subjected to the external forces f (t) and the internal feedback forces caused by nonlinearities f nl (t).This point of view, already chosen in [9] to develop the NIFO frequency domain method, is also on the basis of the present time domain identification method, referred to as Nonlinear Subspace Identification (NSI) [10].Assuming that only displacements are measured, a state-space formulation of the equation of motion, corresponding to the state vector x = z ż T and to the input vector or, in compact form, The proposed nonlinear identification procedure is based on the estimation of the state space matrices Âc , Bc , Ĉ and D, obtained within a similarity transformation by a subspace method in the time domain, and on the subsequent computation of system parameters from the matrix which is invariant under the similarity transformation [10] and contains all system parameters (included in M, C v , K, µ j and L nj ).It is possible to write the matrix H E (ω) as follows where is the underlying linear system receptance matrix.Hence H E (ω) is called "extended" Frequency Response Function (FRF) matrix, because it also includes nonlinear terms in the Multiple Input Multiple Output (MIMO) model.Moreover, in the particular case ω = 0, one obtains the latter expression will be utilised in the numerical examples.

Subspace methods
Linear identification through subspace methods take advantage from robust numerical techniques and, although an exhaustive review is impossible, the books of Ljung [11] and of Van Overschee and De Moor [12] are two milestones to be surely mentioned.In this subsection, only the principal steps of the procedure proposed by [12] are illustrated for brevity's sake.
In the data-driven approach the input data are gathered in a block Hankel matrix the subscript p denotes the "past" and the subscript f denotes the "future" and the number of block rows i is a user defined index.The output block Hankel matrices Y 0|2i−1 , Y p and Y f are defined in a similar manner by replacing u with y in Eq. ( 10) and both input and output data are used to define an oblique projection matrix O i .Then the Singular Value Decomposition of the following weighted oblique projection is performed where W 1 and W 2 are two weighting matrices, whose definition depends on the chosen subspace algorithm.As previously discussed, the state-space matrices in Eq. ( 5) can be obtained only within a similarity transformation, that is function of the above mentioned weighting matrices.The model order n is determined by inspecting singular values and accordingly U 1 , V 1 and S 1 are defined.Finally an estimate of the state space matrices Âc , Bc , Ĉ and D can be obtained, as shown in [12].

Nonlinear identification through feedback of the outputs
This procedure [9] treats nonlinearities like unmeasured internal nonlinear feedback forces that, together with the measured external input, act on the underlying linear system to produce the measured output of the system, as shown by the equation of motion Eq. ( 2).Applying the Fourier Transform results in the frequency domain version with Premultiplying both sides of Eq. ( 12) by the linear FRF matrix which represents the basic formulation of the Nonlinear Identification through Feedback of the Outputs.By comparing Eq. ( 13) with the expression of H E (ω) in Eq. ( 7) it is clear that the NIFO formulation may be viewed as the frequency domain counterpart of the NSI formulation; this is obvious, since both methods start from the same equation of motion and consider the nonlinear terms as internal feedback forces.Moreover, both methods estimate nonlinear structural parameters as well as the underlying linear parameters in a single computational step.
Since nonlinearities convert uncorrelated random noise into correlated noise, the matrix involved in Eq. ( 13) is ill-conditioned, this determining the need of numerical techniques to reduce the error propagation.An approach similar to that discussed in [13] may be used to solve for the nonlinear parameters and the underlying linear FRFs simultaneously.

Numerical example
A numerical example is considered hereafter to show the performances of the NSI technique with simulated data from a MDOF system with a single cubic nonlinear stiffness characteristic.This example demonstrates that the application of subspace algorithms to nonlinear system identification gives good conditioning and computational efficiency even with short input-output time histories.
The four-degree-of-freedom nonlinear system shown in Fig. 1 is excited by a zero-mean gaussian random force (with r.m.s.= 50 N) at DOF 2 only; a nonlinear term is included in the equation of motion Eq. ( 1): with system parameters summarised in Table 1.Time histories are obtained by a Runge-Kutta numerical integration and 1% of zero-mean gaussian noise has been added to each simulated displacement.The first 5 × 10 4 samples with sampling frequency f s = 4 kHz are used for the NSI analysis, choosing i = 60 block rows for the block Hankel matrix in Eq. (10).In this case the extended FRF matrix is where H is the FRF matrix of the underlying linear system.Since the force is only applied at DOF 2, only the second column of this matrix can be estimated.However, the reciprocity relationships H 21 = H 12 , H 23 = H 32 and H 24 = H 42 (which hold since H is related to a linear system) can be applied, and this is sufficient to compute µ 1 from the second element of the following vector (the symbol ?denotes unknown terms) The model order n = 8 is selected by inspecting the singular value plot in Fig. 2, where a jump of two orders of magnitude (more in absence of measurement error) between model order eight and nine is observed.
As shown in Fig. 3, the classic H 1 and H 2 procedures (the mean value is depicted) are unable to estimate the linear FRFs: there is a considerable peak shift for the last two modes.On the contrary, excellent agreement may be observed between the NSI estimate and the true FRF of the underlying linear system (the range between 14.5 and 15.5 Hz has been also magnified to show the differences).
The nonlinear stiffness contribution is estimated by using the invariant matrix H E (ω = 0), which is nonzero since displacement measurements are here considered.From Eqs (15) and ( 16) it is possible to estimate the nonlinear coefficient µ 1 : in Fig. 4 the identified nonlinear stiffness characteristic is compared with the true one and with the linear characteristic, showing a good level of agreement (the case without measurement errors gives excellent results).
Due to its capability of treating many nonlinear terms simultaneously, NSI has been adopted for the identification of the same system (with a single nonlinear element), but seeking for possible nonlinear stiffness elements located at each joint.Clearly, the identified coefficients µ j (2 j 7) are not exactly zero, however, as shown in Fig. 5,  the identified nonlinear contributions are negligible with respect to the linear ones, except between the mass 1 and 3 (where the nonlinear element is located).This example shows the possibility of detecting nonlinear lumped elements (also due to damping) by using the NSI method treating the system as a black box and without much computational effort.NIFO shows similar performances as it will be discussed in the experimental application, while CRP would have required a huge computational effort in constructing a set of uncorrelated response components.

Experimental application
The analysed data were chosen from those proposed by VTT Technical Research Centre of Finland within the framework of the COST action F3 working group WG3 on "Identification of non-linear systems" [14].The structure  shown in Fig. 6 consists of helical wire rope isolators mounted between the load mass and the base mass of an electrodynamic shaker.The acceleration a 2 of the load mass, the acceleration a 1 of the bottom plate, the force F and the relative displacement u 12 between top and bottom plate are measured.Excitation is white noise low pass filtered at 400 Hz, producing records of approximately 5×10 5 samples for each channel, sampled at 4096 Hz, with different levels of excitation and two different load masses.
The choice of the type of nonlinearities still represents a cumbersome matter, since it should be known a priori for the method to be applied to real nonlinear system.However, this information can be obtained by performing a preliminary analysis based on the use of multiple and partial coherence functions between the force vector and any nonlinear vector with attempt form [15].As discussed in the previous section, it is also possible to consider the nonlinearity as sum of several functions and then to remove all of them giving negligible contributions to the  restoring force.In this manner the following type of nonlinearity has been chosen: Both NSI and NIFO analyses have been performed on the time histories corresponding to load mass m 2 = 2.2 kg with excitation level 4 V rms .The first 4 × 10 4 samples have been used for NSI and the first 3 × 10 5 for NIFO, since the latter method requires more samples to give more stable results [10]: this fact represents a drawback of NIFO with respect to NSI.On the other hand, the choice of the model order is a cumbersome matter when applying NSI to experimental vibrating structures, since often it is not possible to find a strong discontinuity in the decreasing order of the singular values, as for the numerical examples (Fig. 2).The use of stabilization diagrams similar to those adopted for a linear analysis is probably useful in order to find the correct model order and to eliminate spurious poles corresponding to computational modes [16].
In the present case, the model order n = 16 has been chosen with i = 60 block rows in Eq. (10), which reduces to a minimum the contribution of computational modes, and results in accordance to the NIFO method have been found.Moreover, in order to quantify the nonlinearity (i.e. to estimate the nonlinear coefficient µ 1 ), Eq. ( 7) should be considered instead of Eq. ( 9), since only the accelerations of the two masses have been measured and the inertance is zero when ω = 0.As a consequence, the estimated coefficient is complex and frequency dependent.However, the imaginary part, without any physical meaning, is much smaller than the real one, and a single value of nonlinear coefficient can be easily obtained by a spectral mean.Figure 7 shows that the peak frequency of the acceleration transmissibility, obtained by linear methods, is shifted toward left, due to softening nonlinear effects.The estimate of the nonlinear coefficient (real part) is also depicted in the same figure, being the spectral mean of the imaginary part (without any physical meaning) much smaller than the spectral mean of the real one and thus negligible.Good agreement between NSI and NIFO can be observed, except near the system resonance, where NIFO underestimates the FRF, since the internal nonlinear forces are highly correlated [9].
Moreover, once the state-space matrices have been estimated by a subspace method, it is easy to reconstruct the output [12].A comparison between measured and reconstructed output for this experimental case is made in Fig. 8, where excellent overlay is achieved.Such capability may be used to predict the experimental response also to other loading conditions, which is often requested in many industrial fields [2].

Conclusions
A time domain method, called Nonlinear Subspace Identification (NSI), is proposed for the identification of nonlinear vibrating structures, by using the perspective of nonlinearities as internal feedback forces.The numerical and experimental examples demonstrate that the application of subspace algorithms to nonlinear system identification in structural dynamics gives good conditioning and computational efficiency.
Furthermore, the capability of the NSI method of simultaneously dealing with many nonlinear functions may be successfully exploited in order to characterize and localize the nonlinearities, being this method well capable of detecting the dominant nonlinearities.Shorter time histories, in comparison with other nonlinear methods such as NIFO, are needed to obtain good results, but the choice of the correct model order may be not immediate, this area requiring further studies.
On the basis of the results discussed in this paper, NSI appears to have the chances to be a robust method for employing in many industrial fields, being its capability of predicting the experimental response widely requested in real structures.

Fig. 8 .
Fig.8.Experimental measurement (solid line) and prediction (dashed line).A zoom is also shown to put in evidence the excellent overlaying.