Nonlinear System Identification of an All Movable Fin with Rotational Freeplay by Subspace-Based Method

: To interpret the nonlinear system identiﬁcation of multi-degree-of-freedom vibrating structures with freeplay, an estimation procedure based on subspace method in time domain is developed in this paper. Two types of freeplay, namely, central freeplay and offset freeplay, are estimated. Besides, the frequency response function matrices of the underlying linear system and overlying linear system are identiﬁed. By considering an all movable ﬁn from a portion of a supersonic ﬂutter test platform containing freeplay on the rigid rotating motion, two series of numerical experiments are applied for demonstration. Results show that the proposed approach indicates a high and reliable recognition of freeplay regions and the nominal linear system. The limitations of the method in this study are also discussed.


Introduction
Nonlinearity is generic in nature, and linear behavior is an exception [1]. Nonlinear vibrating behaviors have been reported in numerous pieces of engineering literature. Some self-excited nonlinear vibrations may lead to a catastrophic disaster, such as the nonlinear flutter in aerospace engineering and the chatter vibrations in manufacture. A few series of literature have discussed the chatter prediction and avoidance as earlier by Fofana [2][3][4], then followed by Urbikain [5][6][7][8]. Moreover, the summary could be found in the recent review [9]. In aerospace engineering, as studied by Tian et al. [10], the trapezoidal-plate-like structure may suffer from different types of nonlinear phenomena, such as jump, internal resonance, periodic, quasi-periodic, and chaotic motions. For a real-life structure with non-negligible nonlinearities, it is essential to provide the precise mathematical model to perform accurate and reliable predictions of the structure's dynamic behavior. Therefore, experimental testing and system identification still play a key role because they help the structural dynamicist to reconcile numerical predictions with experimental investigations [1].
The last several decades have witnessed an increasing interest in the system identification of nonlinear mechanical structures for either academic researches or industry applications. A lot of methodologies were well-established through a variety of studies, such as restoring force surface method [11], NARMAX (nonlinear autoregressive moving average with exogeneous inputs) modeling [12,13], the Hilbert transform [14,15], and Volterra series [16]. For more details of developments and advances on the topic of nonlinear system identification in structural dynamics, readers can refer to the reviews proposed by Kerschen et al. [1] and later by Noël [17]. However, these methods were difficult to apply to complicated structures with more than one nonlinearities. On the other hand, it is difficult to apply the identified mathematical model directly in the control field. reduced-order model. Furthermore, the definition of the underlying linear system and overlying linear system, and the description of two typical kinds of freeplay are also detailed. Section 4 details the identification of the all movable fin with either central freeplay or offset freeplay in its pitch degree of freedom. The performance of the TNSI method for estimating freeplay is also discussed. The frequency response function (FRF) matrix of both the underlying and overlying linear system of the nonlinear all movable fin were identified during these two sets of numerical experiments. Finally, Section 5 summarizes the conclusions.

Nonlinear Subspace Identification of Mechanical Systems
The system identification of a nonlinear mechanical system contains three steps: nonlinearity detection, characterization, and parameter estimation as summarized by Kerschen et al. [1]. However, detection and characterization of nonlinearity are not extended in the paper. The TNSI method discussed in this study is associated with the parameter estimation procedure and to estimate the FRF matrix of the nominal linear system together with the corresponding nonlinear coefficients.

State-Space Nonlinear Feedback Description
The governing equation of vibration of the nonlinear mechanical system under generalized coordinates is where M, C, and K ∈ R m×m are the linear mass, damping, and stiffness matrices, respectively. The superscript m means the nonlinear mechanical system is modeled by an m-DOF finite-element (FE) model. p(t) and q(t) are the generalized external force and displacement vectors, respectively. The nonlinear components are modeled by a lumped form that contains n nonlinearities in the system. g j is the mathematical description of the i-th nonlinearity and needs to be known previously. The location of the j-th nonlinearity is specified by a vector L j , the element of which is 0, 1, or −1.
The unknown nonlinear coefficient is µ j . As subspace methods can obtain an identified state-space model, the second-order vibrating equation in Equation (1) is transformed into a first-order state-space description for better understanding. Defining the state vector x = ( q TqT ) T ∈ R N and moving the nonlinear term to the right-hand side by considering the feedback interpretation [26] as shown in Figure 1, the dynamic equations are recast into where the subscript c indicates the continuous-time system; superscript e stands for extended; A c ∈ R N×N , B e c ∈ R N×(m+n) , C c ∈ R m×N , and D e c ∈ R m×(m+n) are the state, extended input, output, and extended direct feed-through matrices, respectively; and N = 2 m. The extended input vector is The transform between state-space and physical-space model is Equation (2) describe a continuous-time state-space system. However, in practice, the excited and observed signals of a vibrating system are measured by discrete-time sampling. A discrete-time translation of Equation (2) is introduced, and the discrete-time state-space system is where k is the sampled time and subscript d represents discrete-time.

Subspace System Identification
As the nonlinear vibrating equation can be recast into a discrete-time state-space system as presented in Equation (4), the identified state-space model within a similarity transformation can be obtained using the subspace-based method. The procedure of subspace system identification is briefly given. The identification problem statement of subspace identification algorithms is as follows, for a given s measurements of the input u and the output q generated by the unknown system of order N, to determine the order N of the unknown system and the system matrices A d , B e d , C d , and D e d up to within a similarity transformation [33]. In the time domain, the output Hankel matrix is defined as where i is an user-defined index, l is typically equal to s − 2i + 1, and s is the number of data sample.
The subscript p and f stand for past and future, respectively. Similarly, the extended input block Hankel matrix is Then the time-domain output-state-input matrix equation Equations (4) is rewritten as follows by taking account of extended observability matrix Γ i and the lower-block triangular Toeplitz matrix Λ i .
] ∈ R N×l is the state sequence, and Γ i and Λ i are defined as To estimate the observability matrix, the terms in Equation (7), depending on U p , are eliminated using an oblique projection with the symbol as Different left and right user-defined weighting matrices are applied to the oblique projection induce various of subspace algorithms, such as N4SID [34], MOESP [35], and CVA [36]. Then, the singular value decomposition (SVD) of weighted oblique projection, from which the order of estimated model and the estimated observability matrixΓ i are determined. Finally, the estimated matricesÂ d ,B e d ,Ĉ d , andD e d are derived. Readers can refer to [33] for more detailed computation procedures.

Estimation of Underlying Linear System and Nonlinear Coefficients
As mentioned previously, the estimated matricesÂ d ,B e d ,Ĉ d , andD e d are within a similarity transformation. Therefore, the structural parameters can not be extracted directly from these matrices. Fortunately, it is not necessary to compute the similarity transformation matrix that the system parameters can be obtained based on an invariant matrix called extended FRF matrix defined as where the subscript e stands for extended. Here, only in Equation (10), j = √ −1 is the imaginary unit and ω is the angular frequency. The invariance property of the extended FRF matrix up to the similarity transformation has been proofed [27]. However, the identified matrices are discrete-time state-space model and need to transform to the corresponding continuous-time system. In general, the zero-order hold assumption is considered and the transformations are One should point out that, it is difficult to maintain zero-order-hold assumption if the sample frequency is too low. Then, the bilinear transformation should be performed while mapping the discrete-time matrices to continuous-time domain [28].
Assuming the FRF matrix of underlying linear system is H defined likely as Equation (10), the relationship between H e and H is where The underlying linear system and the nonlinear coefficients can be identified from Equation (12).

Freeplay Nonlinearity
Freeplay nonlinearity is a special case of piecewise linear stiffness nonlinearity that the inner stiffness is zero. There are often two typical freeplay in a real structure, namely, central freeplay and offset freeplay. The central freeplay, with the two limits of the freeplay region being symmetrical as shown in Figure 2a, is defined as where q is the generalized displacement and f (q) is the generalized restoring force. k o is the outer stiffness of central freeplay and the subscript o stands for outer. d is the limit of gap, and the freeplay region is δ = 2d.  Because of the existence of freeplay, the nonlinear vibrating system has discontinuous nonlinear stiffness and becomes three piecewise linear systems with different regions. For nonlinear systems with freeplay, the nominal system is the overlying linear system (OLS), that is the system without freeplay (d = 0 as shown in Figure 2a) and full stiffness [37]. Moreover, the underlying linear system (ULS) possessing zero stiffness (k o = 0 in Figure 2a) on the source of freeplay is also important for the nonlinear mechanical system. As for nonlinear aeroelastic analysis, the characteristics of the underlying linear system offer an additional reference to the nonlinear behavior and are paid more attention.
The restoring force of the offset freeplay is depicted in Figure 2b, whereby there is a shift on the generalized displacement, as defined in Equation (15). As freeplay generally occurs due to wear and to loosen of components, it is not necessarily centered [37]. d 1 and d 2 are the left and right limits of the gap, respectively, and the freeplay region is δ = d 2 − d 1 . The underlying linear system and overlying linear system for a system with offset freeplay can also be defined similarly to the one with central freeplay.

Description of the All Movable Fin
The numerical experiments of nonlinear system identification were carried out on a supersonic flutter test platform as shown in Figure 3. The test platform was assembled on a rigid base with an all movable fin connecting by an elastic beam. Based on the convenience of performing the ground vibration tests (GVTs) and the compatibility requirements of mounting on the side walls of the wind tunnel, the rigid base was designed with sufficient weight and stiffness. To simulate the installation environment of the fin, a rotating device was also considered inside the rigid base with rotating fixed. The base was made of steel while the all movable fin, specified as a skeleton-cover structure, was manufactured using aluminum-alloy. With several additional weights inside, the whole mass of the all movable fin is 3.20 kg. The size of the elastic beam was designed to meet the stiffness requirements, which was also made of steel.
The all movable fin was tested and identified by linear system identification techniques to obtain the natural frequency and mode shapes. Figure 4 presents the experimental set-up of GVTs performed using LMS Test.Lab. Burst random excitation was applied on the two selected nodes by a shaker, separately. Furthermore, the displacements of four selected nodes were sampled by four laser sensors. The sampling frequency was set to 2000 Hz to avoid the bias of low-frequency sampling. Table 1 gives the natural frequencies and damping ratios of the first 2 modes of all movable fin up to 150 Hz of interests. The mode shapes of the first two modes are shown in Figure 5.    The finite-element (FE) model was established by MSC Nastran [38] and has been modified based on the test results of GVTs using FEMTools [39], an optimization software developed by Dynamic Design Solutions. Figure 6 is the FE model of the all movable fin and elastic beam. The white-colored component is the all movable fin, whereas the cyan-colored component is the elastic beam. To reduce computational burdens and achieve faster nonlinear calculations, the linear elements of the FE model were condensed using the modal synthesis technique. The whole FE model was reduced in terms of several elastic modes of all movable fin and two assumed modes involve the rigid bending and rotating motion of all movable fin. Freeplay nonlinearity was finally introduced to the rigid rotating DOF. The freeplay regions considered in this paper are listed in Table 2.  A 4-5th-order Runge-Kutta time integration was performed with the same sampling frequency of 2000 Hz as GVTs during the numerical investigation. The mechanical and electrical disturbances in measurements were considered by adding a 0-150 Hz white noise to the synthetic signals, and the noise-signal ratio was set to 2%.

Nonlinear Identification Based on Random Excitation
The nonlinear identification of the all movable fin based on the subspace algorithm was addressed in this section. In the first numerical experiment, a periodic random forcing was applied to the structure with central freeplay. The location of excitation was chosen as same as in GVTs. The responses of the four nodes corresponding to the GVTs were observed by direct time integration, which was performed on the reduced-order model under generalized coordinates description. The nonlinearity in the structure was thoroughly excited by a series of different excitation levels. A second investigation was aimed to freeplay with offset. Both ULS and OLS were identified during the two cases of nonlinear system identification, and the stabilization diagrams were applied to estimate the modal parameters.

Identification of Central Freeplay
As the TNSI method exploited directly Input-Output data in time-domain that will not suffer from leakage and aliasing for no need to transform to frequency-domain, a random noise forcing was applied perpendicularly on the all movable fin at the same exciting location in Figure 4. The excitation signal consisted of 60,000 points over 30 s and band-limited in 0-150 Hz of interests, and repeated 5 times at each excitation level. It also has been discussed in [29] that the transient response was generally found to yield improved results compared with the steady-state regime, especially in the estimation of the linear modal properties. The analysis was conducted using the output data in terms of generalized coordinates established by the modal synthesis method for straightforward identification of freeplay. Therefore, the excitation force has to transform into generalized forces as input data.
However, the accuracy of nonlinear system identification is conditional upon the selection of the mathematical forms g j (t) introduced in Equation (1). In particular, the boundary of different subsystem regions in Figure 2 plays a key rule in freeplay nonlinearity and still remain a challenging task. First, user-defined multiple piecewise-linear functions with varying clearances are used to representing the nonlinearity to be identified. Then, the freeplay can be posterior, determined after the reconstruction of nonlinear restoring force established by the participation of the user-defined functions. In a detailed explanation, we can assume the relative displacement across the nonlinear connection is q r . Then, the displacement range of q r , i.e., [0, max(|q r |)], can be divided into l equally spaced intervals. Given the jth interval [d j−1 , d j ], the nonlinear form is defined as where sign(•) is the sign function. Figure 7 shows the nonlinear coefficients estimated by Equation (12). All estimations of d j indicate that the limit of central freeplay lies at approximately 0.10 • . Then, the reconstruction of the nonlinear restoring force is performed based on the identified nonlinear coefficients µ j and mathematical form g j (t), which is shown in Figure 8. However, it seems that dividing more intervals would not improve or even reduce the accuracy of identification, even though the estimated limits of the freeplay region are within an acceptable range. Therefore, a proper division of the displacement range is essential and the number of the interval l around 10 would lead a better estimation.  An alternative way to improve the accuracy is to narrow the displacement range of varying clearances and repeat a 2nd estimation. A further division of the 0.08 • ∼ 0.12 • range into 10 new sub-intervals based on the 1st estimation of nonlinear coefficient in Figure 7. The identified nonlinear coefficients in the range from 0.08 • to 0.12 • is presented in Figure 9a that a more precise estimate can be inferred, and the identified nonlinear restoring force is shown in Figure 9b, which meets the theoretical nonlinear restoring force with high consistency.  Then, the estimated FRF (H 11 ) of OLS was depicted in Figure 10. TNSI method gives a clear FRF of the OLS; however, FRF estimated by traditional linear H 1 method possesses a significant distortion with no observable peak on the contrary, which may lead to an incorrect estimation of modal parameters. Then, the stabilization diagrams computed by TNSI method indicate the natural frequency and the damping ratio of the overlying linear system in Figure 11 with the max modal order up to 20.
Different from the amplitude independence for the linear system, the dependence on amplitude or frequency of nonlinear behavior is one of the inherent properties of the nonlinear system. Therefore, numerical investigations on varying excitation levels were also considered while applying TNSI method. Figure 12a counted the distributions of nonlinear DOF in different phases as mentioned in Figure 2a. The estimated restoring force is depicted in Figure 12b. The estimations of the freeplay region are excellent; however, the estimation on outer stiffness is still required to be improved by repeating the same identifying procedure in a narrower range.

Identification of Offset Freeplay
Another series of numerical experiments of offset freeplay were carried out with similar simulating parameters as in the previous numerical experiment. The excitation was generated a Gaussian white-noise signal of 60,000 points in 30 s and band-limited in 0-150 Hz repeating 5 times. The analysis was also conducted in terms of generalized coordinates. The identification problem is complicated by the existence of offset in freeplay; however, it is still in reach with utilization of the TNSI method considering some adjustments and skills.
The selection of the mathematical forms g j (t) is different from central freeplay because of the offset. The identification of the left limit and the right limit of the freeplay region are operated separately in two different sets of displacement ranges. Similar to the procedure of identifying central freeplay, the 1st two sets of ranges is based on the boundary of relative displacement response, i.e., [min(q r ), 0] and [0, max(q r )]. The displacement ranges can likely be divided into l − and l + equally spaced intervals, respectively, where the subscripts − and + stand for negative and positive. Given the jth negative interval [d j−1,− , d j,− ], the nonlinear form is defined as Similarly, for the jth positive interval [d j−1,+ , d j,+ ], the nonlinear form is and the nonlinear terms is the summary of g j,− (t) and g j,+ (t). Figure 13 gives the magnitude of nonlinear force contributions identified by TNSI method when l = 15 and l = 20. Both of the results indicate that the left limit of the freeplay region is around −0.05 • , whereas the right one is near to 0.15 • . The reconstruction of the nonlinear restoring force is plotted in Figure 14, indicating an accurate estimation of the inflection points in the offset freeplay. A same procedure is also repeated considering division of the narrowed ranges [−0.07 • , −0.03 • ] and [0.13 • , 0.17 • ] to improve the results as the force contribution in Figure 15a and the nonlinear restoring force in Figure 15b. The modal parameters of OLS can still be extracted from the stabilization diagram in Figure 16.    The identified natural frequencies of OLS based on the two numerical experiments previously are summarized in Table 3.

Conclusions
(1) An identification methodology on a nonlinear system with central and offset freeplay was proposed and demonstrated in this paper. With the utilization of test data under random excitation force, both the underlying linear structure and overlying linear structure were estimated, with the company of estimated coefficients of the freeplay by exploiting the TNSI method. Furthermore, the identified state-space model can naturally be used in the control field.
(2) A series of numerical experiments were implemented. With comparison to the theoretical system, the estimated regions of either central or offset freeplay were accurate using the TNSI method. The capability of the present approach to identifying the offset freeplay indicates a more widely potential application to real-life structures. Combined with the modal synthesis technique, the TNSI method can also identify complicated MDOF structures.
(3) However, one limitation is that the outer stiffness of freeplay may shift depending on the division of response range. In the future, an experimental implementation of the TNSI method is in reach. However, a more robust integrated procedure based on TNSI needs to be developed on practical applications.

Conflicts of Interest:
The authors declare no conflicts of interest.

Abbreviations
The following abbreviations are used in this manuscript: