# Free-Decay Nonlinear System Identification via Mass-Change Scheme

Authors: Dario Anastasio; Stefano Marchesiello
Journal: Shock and Vibration
 (2019)

---

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


---

## Body

## 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):(1)My¨t+Cvy˙t+Kyt+fnlt=0,yt=0=y0,y˙t=0=y˙0.The adopted nomenclature is reported in Nomenclature. The system matricesM, Cv, and K have dimensions ℝN×N in this case, while yt∈ℝN. The term fnlt∈ℝN represents the nonlinear restoring force, i.e., the nonlinear part of the equation, and it is generally a function of both displacements yt and velocities y˙t. As in [19, 20], it is assumed that fnl can be decomposed into J distinct nonlinear contributions using a linear-in-the-parameters model, thus yielding(2)fnlt=∑j=1JμjLjξjt,where μj is the coefficient of the jth nonlinearity and ξj the jth nonlinear basis function, which defines the shape of the considered nonlinearity. The vector Lj∈BN 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 fnl is shifted to the right-hand side of equation (1), becoming a forcing term to the underlying linear system:(3)My¨t+Cvy˙t+Kyt=−fnlt,yt=0=y0,y˙t=0=y˙0.A forcing vector of nonlinear basis functions is defined as(4)ξnlt=−ξ1t,…,−ξJtT,so as to rewrite the equations of motion in a state space formulation. The state vector x=yT,y˙TT can be introduced, yielding(5)xτ+1=Axτ+Bnlξnlτ,yτ=Cxτ+Dnlξnlτ.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 Gω hereafter, in contrast to the original formulation of NSI for input-output data. Instead, the FRF matrix of the nonlinear feedbacks Gnlω can be defined as(6)Gnlω=Dnl+CzI−A−1Bnl,z=eiωΔτ,where I is the identity matrix and i is the imaginary unit. Gnlω has the same structure as the vector of nonlinear basis functions ξnl:(7)Gnlω=Gωμ1L1,…,GωμJLJ.The matrixGω is the abovementioned FRF matrix of the underlying linear system, and since it is unknown, it is not possible to compute the coefficients μj at the current state.Nevertheless, the modal parameters of the underlying linear system can be estimated from the state matrixA by classical eigenvalue decomposition [4], obtaining the natural frequencies ωr, the damping ratios ζr, and the mode shapes ψr, for each identified mode r=1,…,N. The unit-scaled normalization will be assumed hereafter for the mode shapes ψr, while the mass-normalized mode shapes will be called ϕr. The two set of mode shapes are related by the following expression:(8)ϕr=αrψr,r=1,…,N,where αr is the (unknown) modal scaling factor of the rth mode, and it is related to the modal mass mr by the following relation [21]:(9)αr=1mr.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 yields(10)Gpqω=∑r=1NrApqωr2−ω2+2iωωrζr,where rApq=αr2ψprψqr is the residue of the rth mode. Thus, the FRF matrix Gω can be estimated if the scaling factors αr 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 matrixM by a quantity ΔM. Therefore, a new equation of motion of the same form of equation (1) can be written for the modified structure, indicated with subscript ·I:(11)M+ΔMy¨It+Cvy˙It+KyIt+fnlt=0,yIt=0=yI0,y˙It=0=y˙I0.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 AI,BInl,CI,DInl. The set of modal parameters of the modified underlying linear system can be obtained by performing the eigenvalue decomposition of AI, obtaining the natural frequencies ωIr, the damping ratios ζIr, and the unit-scaled mode shapes ψIr for r=1,…,N. Following [17], the modal scaling factors of the unmodified structure can be computed by(12)αr2=ωr2−ωIr2ωr2ΓrrψrTΔMψIr,r=1,…,N,where Γrr is the rth diagonal entry of the matrix Γ, estimated from(13)Γ=Ψ†ΨI,where Ψ† is the pseudoinverse of the modal matrix Ψ containing the eigenvectors ψr of the unmodified structure and ΨI 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 matrixG can be estimated from equation (10), and thus the coefficients of the nonlinearities μj 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(14)μjGLj=μjG11G12⋯G1NG21G22⋮⋱⋮GN1⋯GNNLj=G1jnlG2jnl⋮GNjnl.The coefficientμj can be computed from equation (14) starting from any row of G. 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 G 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 μj, 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 μjid=μjidω, as it is the ratio between complex frequency response functions. Recalling that the true coefficient μj is a real number, it follows that μjidω 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 μjidω 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 Wjω the weighting vector of the jth nonlinearity, it is possible to write the following minimization problem:(15)argminμjidωrjHdiagWjrj,rjω=G1jnl⋮GNjnl−μjidGLj,with rjω being the residue at the frequency ω and superscript H indicating the Hermitian transpose. The LS solution of equation (15) gives the complex-valued quantity μjidω.Moreover, it is possible to split the single modes’ contributions to thejth nonlinearity. Indeed, G is already built as a sum of modes (equation (10)), and the same can be done for Gnl in the discrete state-space formulation as in [13]. In particular, if displacements are measured, the corresponding receptance would be(16)Gnlω=∑r=1Ncrmbrmz−λr=∑r=1NGnlrω,z=eiωΔτ.where crm is the rth column of Cm=CΨ, brm is the rth column of Bm=Ψ−1Bnl, and λr is the rth discrete eigenvalue of A. The contribution of the rth mode to Gnl is called rGnl. If accelerations are measured, the corresponding inertance can be written as(17)Gnlω=∑r=1Nz−1crmbrmλr−1z−λr=∑r=1NGnlrω,z=eiωΔτ.Note that the FRF matricesG and Gnl should be consistent, so if G is a matrix of receptances as in equation (10), then Gnl should be written in a receptance formulation as well. It is therefore possible to rewrite equation (14) for each mode r as(18)μjidGrωLj=rGnlω,r=1,…,N,where rGω is the contribution of the rth mode to Gω. 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 Figure1.Figure 1
Flow diagram of the identification strategy.


## 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 Figure2, where the nonlinear link is also shown, while the system parameters are summarized in Table 1. The nonlinear basis function ξ1t and the location vector L1 are, in this case, respectively,(19)ξ1t=y3t−y1t3,L1=−1010T.Figure 2
4-DOF nonlinear system with a cubic stiffness located between DOFs 1 and 3.Table 1
Numerical example parameters.
Mass (kg)
Linear stiffness (N/m)
Damping (Ns/m)
Nonlinear stiffness (N/m3)
m
1
=
25
;
m
2
=
18
;
m
3
=
15
;
m
4
=
25
k
1
=
3
⋅
10
4
;
k
2
=
k
3
=
k
4
=
5
⋅
10
4
;
k
5
=
k
6
=
k
7
=
5
⋅
10
4
c
1
=
c
2
=
c
3
=
10
;
c
4
=
c
5
=
c
6
=
c
7
=
10
k
nl
=
5
⋅
10
10The free-decay response of the system is simulated by applying an impulsive force on DOF 2 at the timet=0s. The sampling frequency is fs=500Hz, 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.Figure 3
Displacements of the 4-DOF system with zoom on the first 2 seconds. Blue line: DOF 1; orange line: DOF 2; green line: DOF 3; purple line: DOF 4.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 Figure4(a), while the percentage frequency shifts of the four modes are reported in Figure 4(b) taking as starting values the instantaneous frequencies at t=0s. 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.Figure 4
Time-frequency analysis of the 4-DOF system. (a) Spectrogram of the output no. 2. (b) Frequency variations of the four modes. Blue line: mode 1; orange line: mode 2; green line: mode 3; purple line: mode 4.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 Figure5. 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 ω,ζ,ψr=1,…,4 is extracted and reported in Table 2 in terms of natural frequencies and damping ratios.Figure 5
Stabilization diagram for the 4-DOF system. Grey dot: new pole; blue plus: pole stable in frequency; red square: extra stable MAC; green circle: extra stable damping. Stabilization thresholds: 0.5%, 10%, and 99.5% for frequencies, damping ratios, and MACs, respectively.Table 2
Identified modal parameters of the 4-DOF system.
Mode
Frequency (Hz)
Damping (%)
Unmodified
Set I
Unmodified
1
2.88
2.75
0.29
2
14.69
14.00
0.97
3
16.03
15.28
1.01
4
17.76
16.93
1.09The 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ωI,ζI,ψIr=1,…,4 is extracted and reported in Table 2 in terms of natural frequencies. The damping ratios ζI 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.Table 3
Identified modal masses of the 4-DOF system.
Mode
Identified modal mass (kg)
Theoretical modal mass (kg)
Difference (%)
1
31.634
31.648
0.046
2
3.977
3.982
0.137
3
1.061
1.066
0.488
4
7.746
7.778
0.412The FRF matrix of the underlying linear system is then built according to equation (10). The driving point FRF G22ω 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.Figure 6
Driving point linear FRFG22 for the 4-DOF system. Continuous black line: identified FRF; dashed-dotted red line: theoretical FRF.As for the nonlinear identification, the coefficient of the nonlinearityμ1=knl is identified solving equation (18) in a LS sense. Since the full FRF matrix is available, N = 4 estimations of μ1 are retrieved considering the 4 DOFs as inputs for each mode r=1,…,4. 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 μ¯1id of the real part of the LS solution μ1idω. The result is μ¯1id=4.98⋅1010N/m3 with a standard deviation of 1.5⋅105N/m3, providing a percentage error of 0.40% from the true value.The single estimations are reported in Figure7 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 μ1 and the final identified value μ¯1id are depicted. As expected, the contributions associated to the last mode, i.e., the ones labeled as 4G∗∗nl, are generally weighted more. In particular, the contribution 4G31nl has the highest weight, and this can be explained by considering that the fourth mode shape ψ4 has its maximum in correspondence of the third DOF. Thus, the FRF 4G31nl is the most reliable in this case, as it corresponds to the best excitation point of the “most nonlinear” mode.Figure 7
Coefficient of the nonlinearity for the 4-DOF system. Green line: theoretical value. Red line: real part of the LS solution. Dots: single estimations. The intensity of the color of the dots is proportional to their weight in the LS solution.Eventually, real and imaginary parts of the LS solutionμ1id 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.Figure 8
Coefficient of the nonlinearity for the 4-DOF system in logarithmic scale. Continuous black line: real part of the LS solution; dashed-dotted black line: imaginary part of the LS solution.


## 4. Experimental Application

The experimental application consists of five aluminum decks linked by thin steel beams (Figure9) [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.Figure 9
Multistory building with nonlinearity produced by a thin wire.Figure 10
Photos of the experimental setup. (a) Overall view with nonlinear link highlighted in red. (b) Overall view with motion directions in green. (c) Particular of the nonlinear link (thin wire).Table 4
Characteristics of the experimental setup.
Element,i
Plate
Vertical beam
Mass (kg)
Width (mm)
Length (mm)
Thickness (mm)
Length (mm)
Section (mm2)
1
4.30
270
250
24
50
60 × 0.3
2
2.15
270
250
12
30
60 × 0.3
3
1.97
270
250
10
60
60 × 0.3
4
1.79
270
250
10
60
60 × 0.3
5
1.99
270
250
10
60
60 × 0.3A thin pretensioned metallic wire is connected to the fifth floor (Figure10(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 klin plus a cubic one [22], thus having a nonlinear restoring force fnl=μ1y53.The free-decay response is recorded with 5 accelerometers positioned at each floor plus one on the ground, with sampling frequencyfs=409.6Hz and duration of 40 s. The displacement of the fifth floor y5t is obtained by double integrating its measured acceleration y¨5t.The spectrogram of the outputy¨2t 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 t=0s. 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.Figure 11
Time-frequency analysis of the multistory building. (a) Spectrogram of the output no. 2. (b) Frequency variations of the five modes. Blue line: mode 1; orange line: mode 2; green line: mode 3; purple line: mode 4; yellow line: mode 5.Nonlinear system identification is performed with output-only NSI considering the following nonlinear basis functions and location vectors:(20)ξ1t=y5t3,L1=00001T,ξ2t=y5t2,L2=00001T.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 Figure12. Stabilization is checked for frequencies, damping ratios, and MACs, and the set of identified modal parameters of the underlying linear system ω,ζ,ψr=1,…,N 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].Figure 12
Stabilization diagram for the multistory building. Grey dot: new pole; blue plus: pole stable in frequency; red square: extra stable MAC; green circle: extra stable damping. Stabilization thresholds: 0.5%, 10%, and 99.5% for frequencies, damping ratios, and MACs, respectively.Table 5
Identified modal parameters of the multistory building.
Mode
Frequency (Hz)
Damping (%)
Unmodified
Set I
Set II
Unmodified
1
3.49
3.02
3.06
0.49
2
6.11
5.46
6.08
0.48
3
10.20
8.93
8.61
0.25
4
14.38
13.34
14.31
0.18
5
27.00
26.61
26.85
0.53The 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 to1.99kg 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.Figure 13
Linear mode shapes of the multistory building.Output-only NSI is applied with both the modifications, leading to two sets of underlying linear modal parameters:ωI,ζI,ψIr=1,…,N when the mass is on the fifth floor and ωII,ζII,ψIIr=1,…,N 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 Figure14. It can be seen that the MAC between ψ5 and ψI5 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.Figure 14
MAC between the mode shapes of the building in the unmodified and modified configurations. (a) Set I, mass added on the fifth floor. (b) Set II, mass added on the third floor.
(a)
(b)Table 6
Identified modal masses of the multistory building.
Mode
Identified modal mass (kg), output-only NSI
Identified modal mass (kg), SSI
Difference (%)
1
5.47
5.69
4.05
2
6.82
6.26
8.20
3
3.94
3.83
2.67
4
1.34
1.37
2.44
5
0.79
0.47
39.48


### 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 Figure11, 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 Figure15, 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.Figure 15
Stabilization diagram for the multistory building, linear system identification. Grey dot: new pole; blue plus: pole stable in frequency; red square: extra stable MAC; green circle: extra stable damping. Stabilization thresholds: 0.5%, 10%, and 99.5% for frequencies, damping ratios, and MACs, respectively.Table 7
Identified modal parameters of the multistory building, linear system identification.
Mode
Frequency (Hz)
Damping (%)
Unmodified
Set I
Set II
Unmodified
1
3.51
3.05
3.08
0.40
2
6.12
5.42
6.10
0.21
3
10.25
8.94
8.62
0.16
4
14.35
13.34
14.29
0.25
5
27.05
26.47
26.80
0.42A decent correspondence is retrieved between the results listed in Table5 (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 Figure11: 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 G11 and G33.Figure 16
Underlying linear FRFs of the multistory building. Continuous black line: output-only NSI estimation; dashed-dotted red line: SSI estimation. (a)G11. (b) G33.
(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 μ¯1id=5.2⋅107N/m3 (with a standard deviation of 2.4⋅106N/m3) and μ¯2id=6.7⋅104N/m2 (with a standard deviation of 5.2⋅103N/m2).The single estimations are reported in Figure17 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 μ¯1,2id are depicted.Figure 17
Coefficients of the nonlinearities of the multistory building. Red line: real part of the LS solution. Dots: single estimations. The intensity of the color of the dots is proportional to their weight in the LS solution. (a) Cubic coefficient. (b) Quadratic coefficient.
(a)
(b)The contributions associated to the second mode are generally weighted more, except for the third DOF, i.e. the dots labeled as2G31nl and 2G32nl. Interestingly, the third DOF of the mode shape ψ2 is almost a node (Figure 13), thus making the estimation of the coefficients unreliable for the combination u=3 (DOF), r=2 (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.Figure 18
Coefficients of the nonlinearities of the multistory building in logarithmic scales. Continuous black line: real part of the LS solution; dashed-dotted black line: imaginary part of the LS solution. (a) Cubic coefficient. (b) Quadratic coefficient.
(a)
(b)The identified nonlinear restoring forcefnl=μ¯1idy53+μ¯2idy52 is depicted in Figure 19 as a function of y5t. 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.Figure 19
Nonlinear restoring force of the multistory building. Dashed black line: cubic term; dashed-dotted black line: quadratic term; red line: total forcefnl. (a) Overall view. (b) Zoom around the origin.
(a)
(b)The RMS value of cubic component offnl is approximately 5 times higher then the quadratic one, thus the response in mostly symmetric.


## 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 Figure11, 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 Figure15, 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.Figure 15
Stabilization diagram for the multistory building, linear system identification. Grey dot: new pole; blue plus: pole stable in frequency; red square: extra stable MAC; green circle: extra stable damping. Stabilization thresholds: 0.5%, 10%, and 99.5% for frequencies, damping ratios, and MACs, respectively.Table 7
Identified modal parameters of the multistory building, linear system identification.
Mode
Frequency (Hz)
Damping (%)
Unmodified
Set I
Set II
Unmodified
1
3.51
3.05
3.08
0.40
2
6.12
5.42
6.10
0.21
3
10.25
8.94
8.62
0.16
4
14.35
13.34
14.29
0.25
5
27.05
26.47
26.80
0.42A decent correspondence is retrieved between the results listed in Table5 (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 Figure11: 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 G11 and G33.Figure 16
Underlying linear FRFs of the multistory building. Continuous black line: output-only NSI estimation; dashed-dotted red line: SSI estimation. (a)G11. (b) G33.
(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 μ¯1id=5.2⋅107N/m3 (with a standard deviation of 2.4⋅106N/m3) and μ¯2id=6.7⋅104N/m2 (with a standard deviation of 5.2⋅103N/m2).The single estimations are reported in Figure17 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 μ¯1,2id are depicted.Figure 17
Coefficients of the nonlinearities of the multistory building. Red line: real part of the LS solution. Dots: single estimations. The intensity of the color of the dots is proportional to their weight in the LS solution. (a) Cubic coefficient. (b) Quadratic coefficient.
(a)
(b)The contributions associated to the second mode are generally weighted more, except for the third DOF, i.e. the dots labeled as2G31nl and 2G32nl. Interestingly, the third DOF of the mode shape ψ2 is almost a node (Figure 13), thus making the estimation of the coefficients unreliable for the combination u=3 (DOF), r=2 (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.Figure 18
Coefficients of the nonlinearities of the multistory building in logarithmic scales. Continuous black line: real part of the LS solution; dashed-dotted black line: imaginary part of the LS solution. (a) Cubic coefficient. (b) Quadratic coefficient.
(a)
(b)The identified nonlinear restoring forcefnl=μ¯1idy53+μ¯2idy52 is depicted in Figure 19 as a function of y5t. 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.Figure 19
Nonlinear restoring force of the multistory building. Dashed black line: cubic term; dashed-dotted black line: quadratic term; red line: total forcefnl. (a) Overall view. (b) Zoom around the origin.
(a)
(b)The RMS value of cubic component offnl 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.



---

*Source: 1759198-2019-07-07.xml*