Abstract
Methods for nonlinear system identification of structures generally require input-output measured data to estimate the nonlinear model, as a consequence of the noninvariance of the FRFs in nonlinear systems. However, providing a continuous forcing input to the structure may be difficult or impracticable in some situations, while it may be much easier to only measure the output. This paper deals with the identification of nonlinear mechanical vibrations using output-only free-decay data. The presented method is based on the nonlinear subspace identification (NSI) technique combined with a mass-change scheme, in order to extract both the nonlinear state-space model and the underlying linear system. The technique is tested first on a numerical nonlinear system and subsequently on experimental measurements of a multi-degree-of-freedom system comprising a localized nonlinearity.
1. Introduction
Understanding the dynamical behavior of mechanical structures exhibiting nonlinearities has become of great interest in the last decades, due to the continual interest in improving design and performances. In particular, when experimental data are considered, a common way to gather a reliable model is via system identification. This task can be particularly challenging when the system is behaving nonlinearly due to the breaking of the basic principles of linear modal analysis. A variety of methods have been developed to perform nonlinear system identification, and an exhaustive literature review can be found in [1, 2]. The different methods are generally based on different assumptions, but they all share the same need of providing and measuring a persistent forcing input to the considered structure, for instance with an electromagnetic shaker. This is a strict requirement originating from an intrinsic property of nonlinear systems: the breaking of the superposition principle and thus of the invariance of the FRFs [3]. Therefore, it is not possible to quantify the nonlinear behavior if the energy provided to the system is unknown (or unmeasured). Conversely, linear system identification with output-only data is a consolidated practice nowadays, and it is referred to as stochastic identification if the unmeasured input is assumed to be a realization of a stochastic process [4].
In this paper, a method for performing output-only nonlinear system identification from free-decay measurements is presented and validated. The method is based on the nonlinear subspace identification (NSI) technique [5–9], and it will be referred to as output-only NSI. In its original form, NSI is an input-output identification algorithm relying on a nonlinear state-space representation of the system, obtained by considering the nonlinear terms as unmeasured internal feedbacks to the so-called underlying linear system [10]. The measured data are then processed using the subspace formulation, derived from the linear system identification theory [11, 12] and adapted to the nonlinear case. NSI has proved to be very efficient in several occasions, with both localized and distributed nonlinearities [13, 14].
When output-only free-decay measurements are considered, NSI in its original form cannot be used to fully estimate the nonlinear model. Indeed, it is still possible to estimate the modal parameters of the underlying linear system, but a quantification of the nonlinear contribution is not available [15]. For this reason, the classical NSI algorithm is combined in this paper with a mass-change technique to fill the missing information needed to complete the nonlinear model. The mass-change method is generally used in linear operational modal analysis to estimate the so-called scaling factor from the identified modal parameters by adding a known mass to the structure [16–18]. In this work, the same approach is brought to the nonlinear case into the NSI algorithm to estimate the FRFs of the underlying linear system and thus the full nonlinear state-space model.
The output-only NSI method is intended to be used in situations where no forcing input can be provided to the system, either for practical reasons or because it may alter the properties of the structure itself. Indeed, the need of having free-decay data does not allow to consider, for instance, structures excited with ambient excitation. Therefore, the method cannot be addressed to as “stochastic identification”. Nevertheless, the identification of nonlinear structures with unmeasured (stochastic) inputs still remains an unreached goal for the research community [2], and the method proposed in this study can be considered as a first step towards this direction. The methodology is first tested on a numerical dataset of a multi-degree-of-freedom system comprising a polynomial nonlinearity and finally on an experimental test bench of a scaled nonlinear multistory building.
2. Output-Only Nonlinear Subspace Identification via Mass-Change Method
Let us consider the free response of a discrete nonlinear vibrating system with N degrees of freedom (DOFs):
The adopted nomenclature is reported in Nomenclature. The system matrices , , and have dimensions in this case, while . The term represents the nonlinear restoring force, i.e., the nonlinear part of the equation, and it is generally a function of both displacements and velocities . As in [19, 20], it is assumed that can be decomposed into distinct nonlinear contributions using a linear-in-the-parameters model, thus yieldingwhere is the coefficient of the jth nonlinearity and the jth nonlinear basis function, which defines the shape of the considered nonlinearity. The vector is the Boolean location vector of the jth nonlinearity, whose entries can be 1, −1, or 0. The a priori knowledge of the nonlinear basis functions is required; therefore, this identification method can be classified as grey box. Adopting the nonlinear feedback interpretation, the term is shifted to the right-hand side of equation (1), becoming a forcing term to the underlying linear system:
A forcing vector of nonlinear basis functions is defined asso as to rewrite the equations of motion in a state space formulation. The state vector can be introduced, yielding
The reader can refer to Nomenclature. Subspace identification can be performed to identify the state-space matrices in equation (5), rearranging the measured output data into Hankel-type block matrices. The idea is borrowed from the linear subspace identification theory (SI) [11], and detailed steps in the nonlinear case can be found in [5].
It is worth noticing that the state-space model retrieved in equation (5) does not allow to compute the FRFs of the underlying linear system, called hereafter, in contrast to the original formulation of NSI for input-output data. Instead, the FRF matrix of the nonlinear feedbacks can be defined aswhere I is the identity matrix and is the imaginary unit. has the same structure as the vector of nonlinear basis functions :
The matrix is the abovementioned FRF matrix of the underlying linear system, and since it is unknown, it is not possible to compute the coefficients at the current state.
Nevertheless, the modal parameters of the underlying linear system can be estimated from the state matrix by classical eigenvalue decomposition [4], obtaining the natural frequencies , the damping ratios , and the mode shapes , for each identified mode . The unit-scaled normalization will be assumed hereafter for the mode shapes , while the mass-normalized mode shapes will be called . The two set of mode shapes are related by the following expression:where is the (unknown) modal scaling factor of the rth mode, and it is related to the modal mass by the following relation [21]:
In the case of underdamped modes, the FRF of the underlying linear system can be assembled as a sum of single modes’ contributions. In terms of receptance, it yieldswhere is the residue of the rth mode. Thus, the FRF matrix can be estimated if the scaling factors are known.
A common way to estimate the scaling factors in linear operational modal analysis consists of adding some known lumped masses to the structure in order to exploit the changes in natural frequencies and mode shapes. The same technique can be applied in the nonlinear case when NSI is used for nonlinear system identification. Let us assume that one or more known lumped masses are attached to the structure, modifying the mass matrix by a quantity . Therefore, a new equation of motion of the same form of equation (1) can be written for the modified structure, indicated with subscript :
The passages from equations (2) to (5) can be repeated for the modified structure, leading to a new state-space formulation represented by the matrices . The set of modal parameters of the modified underlying linear system can be obtained by performing the eigenvalue decomposition of , obtaining the natural frequencies , the damping ratios , and the unit-scaled mode shapes for . Following [17], the modal scaling factors of the unmodified structure can be computed bywhere is the rth diagonal entry of the matrix , estimated fromwhere is the pseudoinverse of the modal matrix containing the eigenvectors of the unmodified structure and is the modal matrix of the corresponding eigenvectors of the modified structure. It should be noted that if only a truncated set of mode shapes is retained, equation (13) gives just an approximation of the true matrix [18].
Once the scaling factors are computed, the FRF matrix can be estimated from equation (10), and thus the coefficients of the nonlinearities can be estimated from equation (7). It is worth highlighting that the full FRF matrix can be retrieved, which is supposed to be symmetric since it is linear, thus yielding
The coefficient can be computed from equation (14) starting from any row of . This leads to an intrinsic redundancy of the methodology, as there are in principle N estimations of each coefficient. Conversely, when standard input-output NSI is used, the number of known rows of is equal to the number of physical forcing inputs, that is, generally one. Practically, it is possible to solve equation (14) in a least-square (LS) sense with respect to , also considering an appropriate weighting function. A convenient choice is obtained considering that the result of equation (14) gives a complex-valued and frequency-dependent quantity , as it is the ratio between complex frequency response functions. Recalling that the true coefficient is a real number, it follows that is expected to have a flat dependence on the frequency and a null imaginary part. This is true only if both noise and nonlinear modelling errors are absent, which is not the case of real measurements. Therefore, the ratio between real and imaginary parts of can be taken as an indicator of the goodness of the identification, and it can be used as a weight in the LS solution. Calling the weighting vector of the jth nonlinearity, it is possible to write the following minimization problem:with being the residue at the frequency and superscript H indicating the Hermitian transpose. The LS solution of equation (15) gives the complex-valued quantity
Moreover, it is possible to split the single modes’ contributions to the jth nonlinearity. Indeed, is already built as a sum of modes (equation (10)), and the same can be done for in the discrete state-space formulation as in [13]. In particular, if displacements are measured, the corresponding receptance would bewhere is the rth column of , is the rth column of , and is the rth discrete eigenvalue of . The contribution of the rth mode to is called . If accelerations are measured, the corresponding inertance can be written as
Note that the FRF matrices and should be consistent, so if is a matrix of receptances as in equation (10), then should be written in a receptance formulation as well. It is therefore possible to rewrite equation (14) for each mode aswhere is the contribution of the rth mode to . The redundancy of the FRF matrix exploited in equation (14) still holds in equation (18); thus, the latter can be solved in a least square sense as well by rewriting the minimization problem of equation (15) in terms of single modes’ contributions. This operation can hypothetically be repeated for all the N participating modes. In practice, a convenient choice might be to consider the subset of modes most influenced by the nonlinearity and discard the others. This idea is deepened in the following sections considering both numerical and experimental data.
A flow diagram of the whole identification process is depicted in Figure 1.
3. Numerical Application
A four-degree-of-freedom nonlinear system is considered hereafter, comprising a single cubic nonlinear stiffness between DOFs 1 and 3. A representation of the system is depicted in Figure 2, where the nonlinear link is also shown, while the system parameters are summarized in Table 1. The nonlinear basis function and the location vector are, in this case, respectively,
The free-decay response of the system is simulated by applying an impulsive force on DOF 2 at the time . The sampling frequency is , and the total time of the simulation is 60 seconds. Time histories are obtained with the Newmark time integration scheme, and 1% of zero-mean Gaussian noise is added to each simulated output. The displacements of the four DOFs are reported in Figure 3.
The spectrograms of the outputs are also computed to check whether the nonlinearity has been properly activated. The nonlinear behavior is expected to be enhanced during the first instants of the response, as a consequence of the impulsive force. Afterwards, the linear behavior should be predominant. This is reflected to the instantaneous frequencies of the system, which will generally tend to the linear natural frequencies starting from a shifted value. The spectrogram of the output no. 2 is shown in Figure 4(a), while the percentage frequency shifts of the four modes are reported in Figure 4(b) taking as starting values the instantaneous frequencies at . These curves are directly extracted by the spectrogram and interpolated with a polynomial function. A decent frequency shift is detected, especially for the last mode of the structure.
It is also clear from the figure that the nonlinear part of the response expires after almost 10 seconds, as the rest of the time series shows a linear behavior, i.e., constant instantaneous frequencies. This may be a drawback of the method, since there are only a few samples contributing to the identification of the nonlinearity. On the other hand, the decoupling strategy seen in the previous section might help the identification, considering, for instance, just the modes more influenced by the nonlinearities. For the case considered here, it seems that the fourth mode is the “most nonlinear”; thus, it may be a convenient choice for the identification of the nonlinearity, according to equation (18).
The stabilization diagram of the underlying linear system obtained by increasing the model order from 2 to 20 is depicted in Figure 5. Stabilization is checked for frequencies, damping ratios, and MACs, and a model order equal to 8 is eventually chosen. The set of identified modal parameters of the underlying linear system is extracted and reported in Table 2 in terms of natural frequencies and damping ratios.
The whole process is then repeated with the modified structure, where a change in the mass distribution is accomplished by increasing each lumped mass by 10%. The new set of identified modal parameters of the underlying linear system for the modified structure is extracted and reported in Table 2 in terms of natural frequencies. The damping ratios are not reported as they are not needed.
The modal scaling factors and the modal masses are eventually computed from equations (12) and (9), respectively. The identified modal masses are then compared with the theoretical ones in Table 3, where an excellent agreement can be noted.
The FRF matrix of the underlying linear system is then built according to equation (10). The driving point FRF is depicted in Figure 6 and compared with the theoretical one. A great correspondence is also retrieved in this case between identified and theoretical results.
As for the nonlinear identification, the coefficient of the nonlinearity is identified solving equation (18) in a LS sense. Since the full FRF matrix is available, N = 4 estimations of are retrieved considering the 4 DOFs as inputs for each mode . Thus, a total of 16 possible estimations can be obtained. As pointed out in the previous section, not all the modes equally contribute to the nonlinear part of the response, and for the case considered here, the last one seems to be the most affected by the nonlinearity. Indeed, the LS solution should take into account this information, and this can be done automatically by choosing as weighting function the ratio between real and imaginary parts of the single estimations. With this choice, the final value for the nonlinear coefficient can be obtained by considering the spectral mean of the real part of the LS solution . The result is with a standard deviation of , providing a percentage error of 0.40% from the true value.
The single estimations are reported in Figure 7 in terms of spectral mean of their real parts. The darkness of the dots is proportional to their weight in the LS solution. Also, the true value and the final identified value are depicted. As expected, the contributions associated to the last mode, i.e., the ones labeled as , are generally weighted more. In particular, the contribution has the highest weight, and this can be explained by considering that the fourth mode shape has its maximum in correspondence of the third DOF. Thus, the FRF is the most reliable in this case, as it corresponds to the best excitation point of the “most nonlinear” mode.
Eventually, real and imaginary parts of the LS solution are depicted in Figure 8. It should be noted that the imaginary part is always several orders of magnitude lower than the real part, assessing the goodness of the identification.
4. Experimental Application
The experimental application consists of five aluminum decks linked by thin steel beams (Figure 9) [13]. The rig may be reasonably considered as a 5-DOF system, as it can be assumed that the vertical beams provide just a flexural stiffness contribution. Three photos of the experimental setup are reported in Figure 10, and the characteristics of the structure are reported in Table 4.
A thin pretensioned metallic wire is connected to the fifth floor (Figure 10(c)), and it acts like a nonlinear stiffness when the wire undergoes large amplitude oscillations. The restoring force produced by the wire can be written as a series expansion comprising a linear stiffness term plus a cubic one [22], thus having a nonlinear restoring force .
The free-decay response is recorded with 5 accelerometers positioned at each floor plus one on the ground, with sampling frequency and duration of 40 s. The displacement of the fifth floor is obtained by double integrating its measured acceleration .
The spectrogram of the output is shown in Figure 11(a), while the percentage frequency shifts of the five modes are reported in Figure 11(b) taking as starting values the instantaneous frequencies at . These curves are directly extracted by the spectrogram and interpolated with a polynomial function. A relatively high frequency shift is detected for the first mode of the structure, and it progressively vanishes as the mode number increases. The frequency associated to the fifth mode seems not to be affected by the nonlinearity.
Nonlinear system identification is performed with output-only NSI considering the following nonlinear basis functions and location vectors:
A quadratic nonlinearity is also added to DOF 5 to account for possible asymmetries in the nonlinear restoring force, generally present in real structures.
The stabilization diagram of the underlying linear system obtained by increasing the model order from 2 to 20 is depicted in Figure 12. Stabilization is checked for frequencies, damping ratios, and MACs, and the set of identified modal parameters of the underlying linear system is eventually extracted and reported in Table 5 in terms of natural frequencies and damping ratios. The model order for each mode is selected according to the median-damping criterion [23].
The whole process is then repeated with the modified structure, where a change in the mass distribution is considered. In this case, this is accomplished by adding a mass equal to first on the fifth floor and then on the third floor. The reason of this choice can be found when looking at the identified mode shapes of the unmodified structure, reported in Figure 13. It can be seen that the fifth identified mode has a node on the fifth floor. Thus, adding a mass there does not affect the corresponding mode, making the estimation of the fifth modal mass unreliable. For this reason, the latter is estimated by adding the mass on the third floor.
Output-only NSI is applied with both the modifications, leading to two sets of underlying linear modal parameters: when the mass is on the fifth floor and when the mass is on the third floor.
As for the variations in the mode shapes, the MACs between the unmodified and the modified mode shapes are shown in Figure 14. It can be seen that the MAC between and is equal to 1, as expected. Equation (12) is then applied to the two configurations, and the final modal scaling factors are computed by averaging the two set of estimations, except for the cases of unitary MAC for either one of the two sets. The identified values for the modal masses are listed in Table 6. Also, a comparison with the results obtained performing the (linear) stochastic subspace identification on the last part of the response (linear behavior) is presented in the following subsection.
(a)
(b)
4.1. Comparison with the Linear Identification
Since the mass-change method has been originally developed for linear systems, it is useful to check what happens with linear measurements, also to validate the results obtained with the output-only NSI method. Thus, the same free-decay dataset used so far is considered here, but the first part is cut off to let just the linear response to be present. Referring to Figure 11, the instantaneous frequencies stabilize after almost 15 seconds, so this value is chosen as cutting time. Stochastic subspace identification (SSI) is then performed considering the unmodified structure and the two sets of modifications, to retrieve the linear modal parameters and to estimate the modal masses of the unmodified configuration.
The stabilization diagram of the unmodified structure is depicted in Figure 15, and the model order for each mode is selected according to the median-damping criterion. The identified modal parameters for the three situations are listed in Table 7.
A decent correspondence is retrieved between the results listed in Table 5 (underlying linear systems with output-only NSI) and the results of Table 7 (linear identification with SSI). In particular, the deviation on the identified natural frequencies is generally below 1%, while a higher dispersion is retrieved for the damping ratios. This is very common, as uncertainties in the damping estimation are always quite high. Furthermore, the activation of the stiffness nonlinearity in the complete decay response is likely to trigger some nonlinear dissipation phenomenon as well, possibly related to the contact between the aluminum decks and the vertical slender beams or also to the thin metallic wire undergoing large amplitude oscillations.
As for the computation of the modal masses, equations (9) and (12) are applied again, with the considerations previously made about the position of the added mass still holding. The identified modal masses using the output-only NSI method (on the full nonlinear decay) and the SSI method (on the truncated linear decay) are listed in Table 6.
Generally, there is a good agreement between the estimation of the modal masses. The only exception is the fifth mode, showing a high percentage deviation. This can be explained by considering that the identification with SSI is performed by cutting away the nonlinear part of the response from the decay. Since the fifth mode decays faster, it is possible that this mode is badly identified from the truncated linear decay. This is also confirmed by the spectrogram of the output no. 2 in Figure 11: the power associated to the fifth mode is highly reduced after 15 seconds, in contrast to the other modes. Eventually, the FRFs of the underlying linear system are estimated for both output-only NSI and SSI from equation (10). A comparison is depicted in Figure 16 for the driving point FRFs and .
(a)
(b)
As expected, the agreement is very good except around the fifth mode, for the aforementioned reasons.
4.2. Identification of the Nonlinear Restoring Force
After validating the identification of the underlying linear system, the nonlinear part of the model of equation (1) can be estimated as well. According to equation (20), two nonlinear feedbacks are considered in this case, respectively, cubic and quadratic. Thus, two coefficients should be identified solving twice equation (18) in a LS sense. The weighting function is defined again as the ratio between real and imaginary parts of the single estimations, with N = 5 inputs. As for the single modes’ contributions, the first and second identified modes are considered in this case, as they show a higher frequency shift in Figure 11. With this choice, the final values for the two coefficients can be obtained as the spectral mean of the real parts of the LS solutions, leading to (with a standard deviation of ) and (with a standard deviation of ).
The single estimations are reported in Figure 17 in terms of spectral mean of their real parts. The darkness of the dots is proportional to their weight in the LS solution. Also, the final identified values are depicted.
(a)
(b)
The contributions associated to the second mode are generally weighted more, except for the third DOF, i.e. the dots labeled as and . Interestingly, the third DOF of the mode shape is almost a node (Figure 13), thus making the estimation of the coefficients unreliable for the combination (DOF), (mode). This is correctly caught by the weighting function, which puts almost to zero the corresponding weight. Eventually, the real and imaginary parts of the LS solutions are depicted in Figure 18. It should be noted that the imaginary parts are always several times lower than the real parts.
(a)
(b)
The identified nonlinear restoring force is depicted in Figure 19 as a function of . A zoom around the origin of its shape is also reported in Figure 19(b), where the asymmetry introduced by the quadratic term is visible.
(a)
(b)
The RMS value of cubic component of is approximately 5 times higher then the quadratic one, thus the response in mostly symmetric.
5. Conclusions
The purpose of the proposed study is to perform nonlinear system identification of vibrating structures starting from output-only free-decay measurements. The method is based on the nonlinear subspace identification (NSI) technique combined with a mass-change scheme. The final objectives of the identification are the fully nonlinear state-space model, including the nonlinear restoring force and the FRFs of the underlying linear system. Although generally free-decay measurements are not convenient for nonlinear system identification, as the nonlinearity is poorly excited, the decoupling capability of the presented method allows to maximize the confidence in the identification. This is carried out by weighting the single modes according to their participation to the nonlinear behavior. The technique has been tested first on a numerical system involving stiffness nonlinearity and subsequently on an experimental test bench of a scaled multistory building. The latter comprises a localized nonlinearity on the top floor, consisting of a thin wire undergoing large amplitude oscillations. Results have confirmed the capability of the methodology of identifying the underlying linear and nonlinear parameters of the considered systems with satisfying confidence. Therefore, the presented method is suitable if no forcing input can be provided to a nonlinear structure, relying on a much easier free-decay test.
Nomenclature
General Nomenclature| : | Mass matrix |
| : | Stiffness matrix |
| : | Viscous damping matrix |
| : | Time variable |
| : | Frequency variable |
| : | Displacement vector |
| : | Nonlinear restoring force |
| : | Coefficient of the jth nonlinearity |
| : | Nonlinear basis function of the jth nonlinearity |
| : | Location vector of the jth nonlinearity |
| : | Frequency response function (FRF) matrix |
| : | Subscript for the modified structure |
| : | Superscript for “nonlinear”. |
| : | State vector at sampled time |
| : | Dynamical matrix |
| : | Input matrix (nonlinear feedbacks) |
| : | Output matrix |
| : | Direct feedthrough matrix (nonlinear feedbacks). |
| : | Natural frequency of the rth mode |
| : | Damping ratio of the rth mode |
| : | Unit-scaled mode shape of the rth mode |
| : | Mass-normalized mode shape of the rth mode |
| : | Modal mass of the rth mode |
| : | Scaling factor of the rth mode. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.