Harmonic-Balance-Based parameter estimation of nonlinear structures in the presence of Multi-Harmonic response and force

Abstract Testing nonlinear structures to characterise their internal nonlinear forces is challenging. Often nonlinear structures are excited by harmonic forces and yield a multi-harmonic response. In many systems, particularly ones with strong nonlinearities, the effect of higher harmonics in the force and responses cannot be ignored. Even if the intended excitation is a single frequency sinusoidal force, the interaction of the shaker and the nonlinear structure can lead to harmonics in the applied force. The effects of these higher harmonics of the input force on nonlinear model identification in structural dynamics are often neglected. The objective of this study is to introduce an identification method, motivated by the alternating frequency/time approach using harmonic balance (AFTHB), which is able to consider both multi-harmonic forces and multi-harmonic responses of the system. The proposed AFTHB method can include all significant harmonics by selecting an appropriate time step and sampling frequency to guarantee the accuracy of the results. An analytical harmonic-balance-based (AHB) approach is also considered for comparison. However, the inclusion of all significant harmonics of the response in the analytical expansion of the nonlinear functions is often cumbersome. Furthermore, the AFTHB method can easily cope with complex nonlinearities such as Coulomb friction and with multi-degree of freedom nonlinear systems. Including the effect of higher harmonics in the identification process reduces the approximation error due to truncation and very accurate approximation of the balanced equations of each harmonic is obtained. The proposed identification method requires prior knowledge or an appropriate estimation of the type of system nonlinearities. However, the method of model selection may be used for a set of candidate models, and avoiding a dictionary of arbitrary candidate basis functions significantly reduces the computational costs. This paper highlights the important features of the AFTHB method to ensure accurate estimation using four simulated and two experimental examples. The effects of the number of harmonics considered, the modelling error, measurement noise and the frequency range on the quality of the estimated model are demonstrated.

Testing nonlinear structures to characterise their internal nonlinear forces is challenging. Often nonlinear structures are excited by harmonic forces and yield a multi-harmonic response. In many systems, particularly ones with strong nonlinearities, the effect of higher harmonics in the force and responses cannot be ignored. Even if the intended excitation is a single frequency sinusoidal force, the interaction of the shaker and the nonlinear structure can lead to harmonics in the applied force. The effects of these higher harmonics of the input force on nonlinear model identification in structural dynamics are often neglected. The objective of this study is to introduce an identification method, motivated by the alternating frequency/time approach using harmonic balance (AFTHB), which is able to consider both multi-harmonic forces and multiharmonic responses of the system. The proposed AFTHB method can include all significant harmonics by selecting an appropriate time step and sampling frequency to guarantee the accuracy of the results. An analytical harmonic-balance-based (AHB) approach is also considered for comparison. However, the inclusion of all significant harmonics of the response in the analytical expansion of the nonlinear functions is often cumbersome. Furthermore, the AFTHB method can easily cope with complex nonlinearities such as Coulomb friction and with multi-degree of freedom nonlinear systems. Including the effect of higher harmonics in the identification process reduces the approximation error due to truncation and very accurate approximation of the balanced equations of each harmonic is obtained. The proposed identification method requires prior knowledge or an appropriate estimation of the type of system nonlinearities. However, the method of model selection may be used for a set of candidate models, and avoiding a dictionary of arbitrary candidate basis functions significantly reduces the computational costs. This paper highlights the important features of the AFTHB method to ensure accurate estimation using four simulated and two experimental examples. The effects of the number of harmonics considered, the modelling error, measurement noise and the frequency range on the quality of the estimated model are demonstrated.

Introduction
Most practical structures are nonlinear if excited with large forces. Many of these structures have weak nonlinear behaviour due to the nature of the structure or the working condition and can be approximated using conventional linear theories. However, linear theories are not applicable if the nonlinearity is significant in dynamical systems. In order to investigate the exact dynamics of structural systems, a very accurate mathematical model of the system is required. In addition, many practical systems require continuous monitoring to avoid any unexpected failure or damage in the system. Vibration-based structural health monitoring (SHM) is a process to detect any damage by continuously observing any changes in the material or geometry of a structure. To be able to monitor the condition of a structure using periodic measurement, an accurate model of the system is required. Thus, many identification methods have been introduced in the literature. One may refer to review papers [1][2][3][4] for comprehensive summaries of the identification methods.
There are different identification methods which can be categorized in different ways. One way is to categorize them based on the type of measured data used for identification: time-domain or frequency domain. There are various techniques exploiting the time domain data to characterize the nonlinear model of the structure. Masri and Caughy [5] introduced a nonparametric identification method to describe the properties of the system in terms of orthogonal functions. In the so-called force-state mapping method, the measured response and applied force signals in the time domain are exploited to determine the force transmission characteristics of a nonlinear structure or sub-structure in terms of its mechanical state. Then, the transmitted force is plotted in a three-dimensional mapping versus displacement and velocity in order to demonstrate a surface called the force-state map that identifies the nonlinear characteristics of the system. Crawley and Aubert [6] developed the force-state mapping method to identify and quantify the nonlinear characteristics of space truss joints. For the curve-fitting process, they used a least-squares algorithm instead of using orthogonal polynomials. Al-Hadid and Wright [7] compared the use of orthogonal and ordinary polynomials in the curve-fitting process of the force-state mapping technique. They showed multiple advantages of the ordinary polynomial approach, such as being more accurate and using special functions to model non-polynomial types of nonlinearities. They also modified the method to detect the location of discrete nonlinear elements within lumped parameter systems. Kerschen et al. [8] applied the restoring force surface method to a cantilever beam with two different nonlinearities in order to investigate the performance of the method in the identification of nonlinear systems. Firstly, a cantilever beam was considered with a piecewise linear stiffness, and secondly, a bilinear stiffness was taken into account.
Yasuda et al. [9] presented a technique to identify nonlinear multi-degree-of-freedom systems. The. method requires the nonlinear elements of the system to be assumed using polynomial functions and unknown parameters. The measured response and applied force in the time domain are used to express the quantities and nonlinear functions in a Fourier series. The unknown parameters are determined by applying the harmonic balance method. Feldman [10] utilized the Hilbert transform to present a technique for the identification of nonlinear systems. His method is capable of identifying instantaneous modal parameters of the system exploiting the time-domain measured data obtained from free vibration analysis and forced vibration tests with various types of excitation. Feldman [11] later exploited the capability of the modern Hilbert Vibration Decomposition approach to introduce a non-parametric identification method to characterize the nonlinear elastic force functions and backbone curves of the asymmetric nonlinear vibration structures. He also proposed that the accuracy of identifying the time-frequency distribution of a time-varying nonlinear system response signal can significantly be improved by using HT processing along with Hilbert vibration decomposition and congruent functions [12]. Some identification methods have been dually presented in both the time and frequency domains [13]. Adams and Allemang [14] introduced the nonlinear identification through the feedback of the outputs (NIFO method) in the frequency domain, which is the dual of the restoring force surface (RFS) method presented by Masri and Caughy [5] in the time domain. This method determines the parameters of nonlinear models by exploiting internal feedback to be considered as nonlinear elements. Marchesiello and Garibaldi [15] proposed a time-domain subspace identification (TNSI) method of nonlinear oscillating systems considering the nonlinear elements of the structures as internal feedback forces. Despite the very simple formulation of this method, finding numerically stable solutions requires special care to be taken to avoid ill-conditioning in the problem. Hence, they proposed a numerically robust implementation of a time-domain algorithm applicable to nonlinear mechanical systems introduced by Lacy and Bernstein [16]. Likewise, Noël and Kerschen [17] considered the nonlinearities as feedback forces applied to the underlying linear system and introduced a subspace-based identification method (FNSI) utilizing the frequency domain measured data. Noël et al. [13] utilized two subspace nonlinear identification methods formulated in the time and frequency domains, called TNSI and FNSI, in order to identify a strongly nonlinear satellite structure. Another duality can be seen in the conditioned reverse path (CRP) identification method presented by Richards and Singh [18] and its alternative time domain decorrelation technique, referred to as the orthogonalised reverse path (ORP) method, introduced by Muhamad et al. [19]. Kerschen et al. [20] addressed the comparison between the Restoring Force Surface (RFS) method, which is a time-domain approach, and the frequency domain Conditioned Reverse Path (CRP) method for the identification of nonlinear structures.
Wang and Zheng [21] introduced the Equivalent Dynamic Stiffness Mapping (EDSM) technique utilizing the steady-state primary harmonic response in the frequency domain to identify the nonlinear structural elements in dynamical systems. Both parametric and non-parametric nonlinear models were considered. The method requires the response at all coordinates of the structure to be measured or estimated. It is usually impossible to measure all coordinates due to the complicated physics of the structure or a lack of measuring equipment, and hence unmeasured degrees of freedom need to be estimated using expansion methods. This may lead to error in the identification results. Taghipour et al. [22] developed a new optimization based identification framework to minimize the discrepancy between the experimentally measured data and the regenerated response obtained from a mathematical model. Avoiding the requirement of spatially complete measurement, they avoided the inaccuracy resulting from expansion methods to estimate the unmeasured coordinates.
Other than the inaccuracy coming from the expansion methods, there are various factors and different sources of noise or error affecting the identification results. Practically, it is difficult to control the excitation force in experiments, and therefore the measured force and response signals usually include higher harmonics. There are other experimental research works [23,24] showing the high amplitude of second harmonic of the force signal, while this does not result in a considerable higher harmonic response. Assuming that the response is dominated by only the primary harmonic may be acceptable for cases in which higher harmonics have no considerable contribution in the dynamics of the system. However, this assumption may result in significant errors in the identification of strongly nonlinear systems (e.g. systems with frictional contact interfaces). One of the most important sources of inaccuracy in identifying nonlinear systems that have received less attention in the literature is neglecting the effect of higher harmonics, particularly where higher harmonics play a significant role in the response. Some research studies focused on the nonlinear vibration analysis considering the effect of higher harmonics [25], or investigated the effect of neglecting the higher harmonics on the results of identification [26]. Therefore, the objective of this study is to propose an identification method to consider the effect of significant harmonics of the multiharmonic dynamics of the system. To this end, a method is introduced with two approaches: Analytical Harmonic-Balance-based (AHB) approach and the Alternating Frequency/Time approach using Harmonic Balance (AFTHB). The AHB approach is based on expanding the nonlinear functions in frequency domain implementing analytical methods such as Fourier Integral (FI), Complex Averaging (CXA) technique [27,28] or Harmonic Balance Method (HBM) [29]. Indeed, this approach is the extended version of EDSM technique [21] for multi-harmonic identification. This approach exploits the measured data in the frequency domain. In order to identify nonlinear systems using AFTHB approach, the nonlinear functions are expanded in the frequency domain by calculating the Fourier Transform of the measured time response. It is shown that AHB approach is very accurate in theory if all significant harmonics of force and response are taken into account. Nevertheless, it is often cumbersome to include all significant harmonics of the response in the analytical expansion of nonlinear functions. Besides, it would be difficult to use AHB for the structures with complex forms of nonlinearities such as Coulomb friction and also for multi-degree of freedom nonlinear systems. Therefore, the AFTHB approach is developed based on the corresponding method used for harmonic balance [30,31]. It should be noted that the AFTHB approach is different from well-known time-frequency methods such as the Wavelet and Hilbert Transforms. The great advantage of AFTHB is that it can be applied to any type of nonlinear functions, where there are significant harmonics in the response. It is worth noting that selecting the appropriate time step and the sampling frequency is important to guarantee the accuracy of the results obtained by the AFTHB method. However, many well-developed algorithms such as MATLAB ode solvers automatically adjust the time steps during time marching.
This study develops a method with two different approaches to identify nonlinear dynamical systems exploiting multi-harmonic responses and force signals. Section 2 provides a discussion on various types of model identification and discusses different model selection approaches. The proposed identification procedure is then described; the analytical (AHB) and the numerical alternating frequency/time (AFTHB) identification approaches are introduced. In Section 3, four simulated examples of nonlinear systems are considered and used to highlight different features of the method. Section 4 applies the AFTHB approach to two experimental examples. Finally, Section 5 gives brief conclusions of the study.

Identification procedure
Section 2.1 discusses and compares parametric and non-parametric methods for the identification of nonlinear systems. In addition, the advantages and disadvantages of various model selection approaches are discussed. The assumptions made for the proposed identification method are then given before the proposed method is described in Section 2.2.

Parametric vs Non-parametric identification
Nonlinear identification methods can be categorized into two main classes: non-parametric [6,[32][33][34] and parametric [35][36][37][38][39][40][41][42][43][44][45][46][47][48]. The latter group can be divided into two sub-classes: methods with a model selection stage [35][36][37][38][39][40][41][42][43] and methods using an assumed model for the nonlinear system [44][45][46][47][48]. Despite the efforts in developing nonlinear identification methods, there is not yet any identification method that can be generally applied to an extensive range of nonlinear systems. Indeed, the methods from each of the two categories have their own advantages and disadvantages. Most importantly, nonparametric methods do not require prior knowledge or assumptions for the type of nonlinearities, which is an advantage over parametric approaches. However, the high computational cost is a common drawback of nonparametric methods. The time-domain method of force-state mapping (FSM) [32,33] is a commonly used non-parametric nonlinear identification method in structural dynamics. Despite the simple formulation and its non-parametric nature, the high computational costs of FSM limit the applications of this method to the simple practical structures with few degrees of freedom.
To avoid the necessity for surface fitting, which makes FSM computationally expensive, Karaagaçlı and Ö zgüven [34] developed a frequency-domain non-parametric identification method by determining a describing surface of nonlinearity exploiting responsecontrolled stepped-sine measured data. This method is shown to have good performance in identifying systems with relatively simple nonlinearities. The method has the advantage of characterizing the systems with nonlinearity as a function of frequency. However, the method has convergence problems for more complex structures with non-smooth nonlinearities. To overcome this problem, Karaagaçlı and Ö zgüven suggested applying parametric curve fitting to the extracted non-parametric surfaces or using higher order finite differences to obtain the describing surface rather than first-order finite differences. Although this may solve the convergence issue and even increase the accuracy of the result, it will increase the complexity of the method.
Selecting an appropriate model that truly represents the structure is very important. Various model selection algorithms have been introduced in the literature, including forward, backward, and exhaustive search methods [35][36][37]. Although these methods are applicable to linear-in-parameter nonlinear models, they cannot be used for other nonlinear systems. Some identification methods utilize a dictionary of basis functions [21,[38][39][40][41] to select the type of nonlinearity rather than using prior assumptions for the type of nonlinear force. Although this approach is able to consider a wide range of nonlinear functions in the candidate pool, it has high computational costs, and some nonlinearities are difficult to represent with suitable basis functions. Fuentes et al. [39] introduced a Bayesian-based equation discovery method to overcome these problems by setting to zero all parameters of the basis functions not related to the model using a sparsity-inducing prior within a Relevance Vector Machine (RVM) framework. However, the selection of the basis functions is difficult, and various statistical methods have been investigated. Brunton et al. [40] used a manual selection of a threshold that gives an appropriate set of basis vectors for each problem, although the bias induced by the implicit regularization is a limitation. Ben Abdessalem et al. [41] investigated the performance of various approximate Bayesian computations for the purpose of model selection. Their method exploits a predefined dictionary of candidate nonlinear models, but possible combinations of candidate models are missed. In addition, due to the recursive numerical simulation required to select the appropriate model, their method has a high computational cost. Considering the limitations of model selection, an empirical approach (i.e. expert engineer) to make a proper estimate for the type of nonlinearity based on the experimental observation is still a reliable option for nonlinear model identification methods. Exploiting engineering insight avoids the high computational costs of an extensive dictionary of candidate basis functions. This approach forms the basic assumption for a large number of parametric identification methods [44][45][46][47][48]. Although using incorrect nonlinear functions as the assumed type of nonlinear force may lead to inaccurate model identification, using this approach significantly reduces the computational costs. Furthermore, this approach ensures that the selected mathematical model for the nonlinear structure under study is physically meaningful.
Methods of non-parametric model identification and parameter estimation of assumed models are both important in the analysis of nonlinear systems. For many systems, such as nonlinear springs or frictional contact interfaces, it is straightforward to estimate the nature of the model and the important part in characterizing nonlinear model identification. For instance, there are well known models such as the Coulomb friction model, the Lu-Gre model, the Dahl model, the Valanis model, the Iwan model, and the Maxwell slip model, which are usually used in modelling contact interfaces. Using an appropriate identification approach that leads to physically justifiable parameters for these models is also important in structural dynamics.
Accurate knowledge of the properties of the underlying linear system of the nonlinear structure is another important issue. Updating the underlying linear model of nonlinear structures exploiting experimental data obtained from very low-amplitude excitation tests is a practical approach and widely used [46][47][48]. However, this approach may not be applicable to systems with non-smooth nonlinearities such as Coulomb friction. Some identification methods characterize the mathematical model of nonlinear structures without knowledge of the underlying linear system. The physics of the underlying linear system may be modelled and the parameters estimated or updated during the nonlinear identification process.
Guisquet and Amabili [48] presented a procedure to characterize the nonlinear damping and stiffness parameters of a singledegree-of-freedom model using large-amplitude experimental data obtained from a continuous system excited by a harmonic force. They based the parameter estimation of the presented procedure on the harmonic balance method in which they approximate the nonlinear model using a Fourier expansion with different truncation orders. The procedure requires the type and location of the nonlinearities to be known a-priori and the underlying linear system to be updated. In the present study, an analytical harmonicbalance-based (AHB) frequency domain method is introduced. It is shown in the present study that applying the AFD method to systems with complex nonlinearities and significant contributions of higher harmonics in measured force and responses will be cumbersome. Therefore, an alternating frequency-time harmonic-balance-based numerical approach (AFTHB) is developed to overcome the issue of complexity of considering higher harmonics. Hence, the proposed method is based on the following assumptions: • The method requires the nonlinear functions to describe the nonlinear restoring forces in the model to be known a-priori. Engineering insight to know the appropriate type of nonlinearity in the system is essential to guarantee the accuracy of the identified nonlinear model. • The underlying linear model of the structure is updated using the experimental data obtained from very low-amplitude excitation tests prior to the nonlinear identification process. Keeping the excitation amplitude at a very low level assures the nonlinearities of the structure are not excited. However, the proposed method is also able to estimate the parameters of the underlying linear model during the nonlinear identification process.
Thus the proposed method is categorized as a parametric method that requires a-priori assumptions of the type and location of the nonlinearities.

Proposed identification method
The equation of motion of a general vibrating system subject to external forces is given by where M, C, K ∈ R N×N denote the mass, damping, and stiffness matrices, respectively, u is the vector of displacements, N denotes the number of degrees of freedom, f NL (u,u) represents the unknown linear/nonlinear restoring force, and f ex is the external force vector applied to the system. The nonlinear restoring force is assumed to be a combination of linear or nonlinear functions, which is linear in the unknown parameters. Thus, where G ∈ R N×M is the matrix of assumed functions g i,j (u,u) that represent the internal nonlinear restoring forces and a ∈ R M×1 is the vector of M unknown parameters, a p . Determining the form of the nonlinear functions requires detailed engineering insight, as highlighted in Section 2.1, and examples will be given later. Rearranging Eq. (1) and using Eq. (2), the unknown nonlinear restoring force may be written as Often, in practical dynamical systems, the applied force and the corresponding response will be periodic. Define the fundamental frequency of both the periodic force and response to be ω; note this is the fundamental frequency for both the force and response, and if a subharmonic response is excited, then this will be lower than the fundamental frequency of the force. Then the external force may be written as the truncated Fourier series where N f is the number of harmonics considered, F 0 is the constant static force, and F n is the complex force amplitude for the n th harmonic. I denotes the imaginary part of a complex number. The real part of the complex functions may be used rather than the imaginary part in the development of the harmonic balance method; it is the use of complex quantities that makes the development simpler. Similarly, the steady-state periodic response of the system is given by where N u denotes the total number of harmonics of the response, which is not necessarily equal to N f . u 0 denotes the static response and U n is the complex response amplitude for the n th harmonic. As the nonlinear internal forces, G(u,u), are functions of the periodic displacement and velocity, these forces will also be periodic and can be written as the Fourier series where N g denotes the total number of harmonics considered, which is not necessarily equal to N u , V 0 is the constant static part of the nonlinear functions, and V n is the complex coefficient of the nonlinear function for the n th harmonic.
Since we have the measured response, one option is to directly reconstruct G and simply take the FFT of both sides of Eq. (3). Although this works in theory, the accuracy of the reconstruction needs to be considered carefully for experimental data, and therefore a more sophisticated approach, based on the alternating frequency-time harmonic balance method, is required. This is especially important when velocity and displacement signals within the G functions must be obtained from measured acceleration signals. Numerical integration of acceleration signals gives errors in the phase, which greatly affects the accuracy of the identified parameters, especially for the damping terms. Reconstruction from stepped sine tests using measured data for the harmonics of the response can help, although careful sampling is required.
Substituting Eqs. (4), (5) and (6) into Eq. (3), and equating coefficients of each harmonic, gives where E n ∈ C N×1 is the equivalent force, given by where C denotes complex numbers and . Equation (7) gives N linear simultaneous equations for each of the (N H + 1) harmonics considered (including the static response), that is N × (N H +1) equations, for the M unknown parameters in the vector a. The system of linear equations can be written conveniently as where Since the parameters in a are real, Eq. (9) may be split into real and imaginary parts, and hence where R and I denote the real and imaginary parts.
Note that Eq. (9) is a set of equations for the unknown parameter vector, a, at a particular fundamental excitation frequency, ω.
Typically, a number of excitation frequencies will be measured, and the option is to estimate the unknown parameters independently at each excitation frequency. If the assumed nonlinear model is correct, then the estimated parameters would be constant as the frequency varied. The alternative is to combine the measurements for different excitation frequencies into one set of equations for the unknown parameters as ⎡ where ω i is the i th excitation frequency, and N ω is the number of frequencies. If the force amplitude or the spectral composition of the force is also changed between experiments, then these data could be included in Eq. (11). In practical cases, care should be taken in selecting the frequency points to guarantee the accuracy of the results. Nonlinear terms at frequencies away from resonance might not be excited or too small, and accordingly, the contribution of the higher harmonics to the response would be the lowest. On the other hand, the response is contaminated by noise in practical cases, and away from resonant frequencies, noise level makes it difficult to identify the system. Hence, normally points around resonance with high amplitude would be the best choices to minimize the effect of the noise on the identification results. The parameters may now be estimated from the linear equations given in Eq. (9), or equivalently Eq. (11) if V and E are complex, or Eq. (12) if multiple experiments are used. If the number of harmonics in the response is equal to or greater than the number of unknown parameters (N H ≥ M), one can often find the parameters explicitly using a least-squares (pseudo-inverse) method for each excitation frequency. Sometimes this estimation may be ill-conditioned, for example, if the number of unknown parameters is much larger than the number of degrees of freedom and measured harmonics. Or the amplitude of the higher harmonics of the force or response may be too small to accurately estimate the parameters because of noise and ill-conditioning. In this case, a regularised solution may be obtained. Furthermore, the magnitude of the real and imaginary parts in Eq. (11) will depend on the phase of the signal and may cause issues with the pseudo-inverse calculation if one of these components is zero or close to zero. To minimize the effect of this issue, all the real and imaginary parts of the different harmonics are included in the estimation, and the conditioning of the pseudo-inverse is checked.
Furthermore, the assumption was that the internal nonlinear force is linear in the parameters, i.e. represented by Eq. (2). Sometimes this is not possible (for example, the parameter is in an exponent), and in these cases, the system may be identified using an iterative optimization procedure.
The accuracy of the identification process depends on the accuracy of the calculation of the frequency content of the signals. This requires the calculation of the coefficients V n in Eq. (6) that are used to assemble the matrix V. Two approaches are now introduced to estimate these coefficient matrices. This analysis is established based on the harmonic balance approach used to estimate the simulated response of a nonlinear system. Indeed, the present approach utilizes the capability of the harmonic balance approach to consider many higher harmonics of force and response of the structure.

Analytical harmonic balance (AHB) approach
In this approach, an analytical method (e.g. the modified complex averaging technique (MCXA) or the Fourier integral) is used to determine the multi-harmonic analytical expression to approximate the mathematical model of Eq. (6) in the frequency domain. In this study, MCXA [27,28] is applied to determine the harmonics of the nonlinear functions given in Eq. (6) as functions of the assumed response harmonics in Eq. (5). It is worth noting that: • For the case of single harmonic approximation, this approach is the same as the Equivalent Dynamic Stiffness Mapping (EDSM) technique developed by Wang and Zheng [21]. In other words, the AHB approach is the extended version of the EDSM method for multi-harmonic identification. • The AHB approach is applicable to nonlinear systems with smooth nonlinearities. An exact multi-harmonic expression cannot be derived for non-smooth nonlinear functions such as Coulomb friction, and so the AHB approach may not be applicable for the systems with such nonlinearities. • The AHB approach is developed in this paper to show that increasing the number of harmonics used will improve the accuracy of the identification results. Thus, for systems with a multi-harmonic response, neglecting higher harmonics in the identification process may result in considerable error.

Alternating Frequency/Time approach using harmonic balance (AFTHB)
Although the AHB approach is accurate, it is often cumbersome to apply to complicated models, especially in systems with many DOFs and non-smooth nonlinearities. Section 3 demonstrates that applying the AHB approach with higher harmonics becomes very messy even for a cubic nonlinearity, even though the results are accurate. Therefore, the alternating frequency/time (AFTHB) approach used for harmonic balance is proposed in this section [30,31]. The proposed approach is based on a time-frequency transformation of the general vibrating system of Eq. (1). Given the assumed type of nonlinearity and the measured time response u(t) given by Eq. (5), the nonlinear functions G(u,u) of the internal force can be calculated as a function of time. Calculating the Fourier transform of this nonlinear restoring force gives the force in the form of Eq. (6), i.e. the amplitudes V n . The steps required to calculate the coefficient matrix V are: i. Calculate the time period T = 2π ω , and determine how many samples per cycle to use, as a power of 2. This really depends on how many harmonics are required to be taken into account and can be relatively low (e.g. 32 or 64). Sample time points are obtained as t i = i T N , where N is the number of samples. Although increasing the number of samples improves the precision of the time response and hence the results of identification process, a large number of samples may lead to a significant increase in the time required for the experimental and numerical study (computational cost). Therefore, an optimum number of samples should be selected.
ii. Sample the nonlinear function G(t i ) using the periodic response u of the system with time period T. iii. Calculate the Fourier Transform (FT) of the sampled function G(t i ). The FT gives the coefficient V n of the nonlinear function at each frequency. iv. Assemble the coefficient matrix V using the calculated frequency components V n .
It should be noted that the AFTHB approach is different from well-known time-frequency methods such as the Wavelet Transform. The great advantage is that this approach can be used for either smooth or non-smooth nonlinear functions. In addition, almost all of the significant harmonics of the function (depending on sampling frequency) are immediately available; hence this approach makes it feasible to consider the higher harmonics of the force and response.
In the following sections, several simulated and experimental case studies with different types of nonlinear elements are considered. For each case, the nonlinear mathematical model of the system is identified utilizing various methods including the EDSM technique (which is equivalent to the AHB approach considering only primary harmonic), the AHB approach with multi-harmonic approximation, and the AFTHB approach. The results are compared to investigate the applicability of the presented methods.

Simulated case studies
In this section, the proposed identification method with both AHB and AFTHB approaches is applied to four simulated case studies with different types and various numbers of nonlinear elements and excitation force. The objective of using these examples is to show the characteristics of the proposed identification method for simple and known systems.
In the presence of multi-harmonic responses and excitation force, considering only the primary harmonic may result in inaccurate estimates of the system parameters. For instance, even for the simple case of a SDOF Duffing system, it is shown that neglecting the effect of higher harmonics may result in 1-3% error in the estimated parameters. Although 1-3% may not be considered as a high level of inaccuracy, the error due to neglecting higher harmonics in practical systems with more significant higher harmonics would be much higher. The effect of higher harmonics in the identification of non-smooth nonlinearities is illustrated using an SDOF system with Coulomb friction and the results of the proposed method are compared with EDSM, which neglects the higher harmonics. A multi-DOF system is used to assess the performance of the proposed method for various nonlinear elements and multi-harmonic excitation forces. The effect of modelling error and measurement noise is also demonstrated.

Duffing oscillator with harmonic excitation
Consider the equation of motion of a single degree of freedom Duffing oscillator, excited by a harmonic force, given by In this case Now the two approaches described above are applied to the Duffing oscillator system in order to find the unknown parameter k 2 . AHB Approach: The MCXA technique is applied to Eq. (13) in order to obtain the equation of motion in the frequency domain. If three harmonics of the response are considered, then averaging over the 1st, 2nd, and 3rd harmonics of the response gives the three algebraic equations where the overbar denotes the complex conjugate. Assuming the force and response harmonics are measured (i.e. F, U 1 , U 2 , U 3 ), then the unknown parameter k 2 can be identified using each equation of Eq. (15), or from a least-squares solution. However, the small amplitude of the higher harmonics may lead to ill-conditioning and, in this section, only the first equation, which is related to the primary harmonic of the system, is used to find k 2 as Neglecting higher harmonics and utilizing only the primary harmonic of the response, k 2 is obtained as AFTHB Approach: Using this approach, the harmonics generated by the cubic nonlinearity (the terms multiplying by k 2 in Eq. (15)) are generated numerically rather than analytically. For the primary harmonic, this will give V 1 which includes all of the significant harmonics participating in the response of the system.
The parameter values used for the numerical simulation are given in Table 1. For a given excitation frequency, the response is obtained by integrating the differential equation, and the response harmonics are extracted. The accuracy of the identified parameter, k 2 , using different methods as a function of excitation frequency is shown in Fig. 1. As shown, using the EDSM technique and identifying the parameter exploiting only the primary harmonic would lead up to a 3.5% error in the estimated value of the nonlinear stiffness. Considering three harmonics in the calculation will significantly improve the accuracy of the results and reduce the error below 0.01%. The results of AHB approach would be even more accurate if more harmonics are considered. However, the most accurate value of k 2 (with errors of less than 0.002%) has been identified using AFTHB approach in which the significant harmonics (the first five harmonics) of the response are taken into account.

Duffing oscillator with Multi-Harmonic excitation
Consider the Duffing oscillator of Eq. (13), excited by a multi-harmonic force including the primary and third harmonics. Thus where F 1 = 2N, F 2 = 1N, ϕ = π 2 and the other parameters are given in Table 1. The first five harmonics of the response are considered to investigate the steady-state dynamics of the system using MCXA. Fig. 2 shows the amplitude and phase of the primary, third and fifth harmonics of the response. The higher harmonics, particularly the third harmonic, clearly play a significant role in the response of the system. Due to the nature of the nonlinearity of the system (cubic), only odd harmonics are excited. The primary harmonic is dominant after ω = 5.06 rad/s, although the third harmonic is still significant in the response. For lower frequencies, the higher harmonics can dominate the response. Hence, it is vital to account for the higher harmonics in the identification of the nonlinear model.
Suppose that the equation of motion is averaged over the first harmonic to identify the unknown parameter k 2 of the nonlinear force, as for Section 3.1. Applying EDSM technique, the unknown parameter is estimated using Eqs. (16) and (17) considering only the first three and the primary harmonic, respectively. In order to consider the effect of the first five harmonics of the response into account in the AHB approach, the unknown parameter is estimated as where V 1 is the complex amplitude of the cubic nonlinear function expanded utilizing five harmonics of the response. Thus, The unknown parameter may also be estimated using the AFTHB approach, where V 1 is estimated numerically and then used in Eq. (19) to calculate the nonlinear stiffness. Fig. 3 compares the estimated k 2 using various approaches, versus the excitation frequency. The response from the numerical solution was used for the identification, and hence only the higher amplitude branch from the upward frequency sweep was used. The identification via EDSM technique and utilizing a single harmonic of the response (neglecting higher harmonics, yellow circles in   higher harmonics to the identification make the results more accurate, especially for higher excitation frequencies. The unknown parameter k 2 is also estimated using the AFTHB approach and is shown to give accurate results. Overall, it can be concluded that neglecting higher harmonics, especially where they play a significant role in the response, leads to inaccuracy in the parameter identification.

Single DOF system with Coulomb friction
The equation of motion of a single DOF mass-damper-spring oscillator with Coulomb friction is The nonlinear force in this system may be written in terms of the unknown parameter and nonlinear function as a ≡ a = μ, Thus, the equation of motion averaged over the primary harmonic of u(t) and g 1 (u,u) is sufficient to estimate the value of the single unknown parameter μ. Utilizing AHB approach exploiting the single harmonic expansion of Coulomb friction force in the frequency domain, the unknown friction coefficient μ may be calculated as [21,35] Note that there will be significant error in the estimate of the nonlinear friction force when only the primary harmonic is considered in Eq. (23). However, including more harmonics in the MCXA approach is very difficult.
To apply the AFTHB approach, the nonlinear force must be computed over one period of the response, and the fundamental component of g(u,u) obtained, namely V 1 . Then μ is computed as, The steady state response of the system is simulated using the parameter values given in Table 2, and the friction coefficient estimated. Fig. 4 compares the identified values of μ for different excitation frequencies. The results show that using the EDSM technique (AHB approach with only the primary harmonic) and neglecting the higher harmonics has led to considerable error, while the identified μ using AFTHB approach is close to its true value. In the vicinity of the resonance (ω ≈ 4rad/s) the primary harmonic is

Multi-DOF nonlinear system
In this section, the 6-DOF nonlinear discrete system shown in Fig. 5 is considered. The system includes various types of grounded or ungrounded nonlinear elements such as Coulomb friction, nonlinear damping, and cubic stiffness. The equation of motion of the system is derived in matrix form as where the linear mass and stiffness matrices,M and K are easily derived, and proportional damping is assumed so that C = αM + βK.
f NL and f ex denote the nonlinear and external forces, respectively, and are given by  μ denotes the friction coefficient, g is the gravitational acceleration, k N1 ,k N2 , and k N3 are the cubic stiffness coefficients, and c N is the coefficient of nonlinear damping. F and ω denote the amplitude and frequency of the excitation force. The model of the nonlinear forces is given by Eq. (26), which is in the form of Eq. (2), where the parameter vector a and matrix of nonlinear functions, G, are given by and The parameters of the simulated system are To investigate the applicability of the introduced AFTHB identification method, two different cases of the MDOF system of Fig. 5 are considered: (i) the system without the friction force, i.e. μ = 0, and with and without modelling error; and (ii) the system with Coulomb friction force, without and with measurement noise.

Case I: MDOF system without friction
The simulation of the steady state dynamics of the system of Fig. 5 is performed neglecting the friction between the mass m 1 and ground, i.e. μ = 0. The model of the nonlinear forces is given by Eqs. (26) and (27) with μ = 0, which means that the dimension of a and G are reduced by removing μ, g 2 and the corresponding zero row and column. The dynamics of the system is simulated for the frequency range of ω = 20 − 26rad/s with a frequency step of Δω = 0.05 rad/s, using the four levels of excitation F = 2, 10, 20, and 50N. At each frequency point, the dynamics of the system is simulated for 600 cycles with 2 10 data points at each cycle. The first 300 cycles of the response at each frequency are discarded to ensure the steady-state dynamics are measured. The AFTHB method is then applied to identify the parameters of the nonlinear forces using the data from all excitation frequencies and force levels simultaneously. Fig. 6 shows the phase diagram of the steady-state response of x 1 at excitation frequency ω = 24rad/s and force level F = 20N, and shows that the response is periodic with many harmonics excited. This process is repeated for the full range of excitation frequencies and force amplitudes, and the complex response amplitudes for the first nine harmonics are extracted. However, as the amplitude of the even (2, 4, …) harmonics are very small and hence negligible, only the odd harmonics from 1 to 9 are used to identify the parameters of the system. The estimated parameters for Case I are given in Table 3. The true parameters are also given and the comparison shows that AFTHB identification method gives highly accurate results with no model error.
The identification assumes knowledge of the form of the nonlinear restoring force so that the parameters may be identified. In practice, the nonlinear function will not be known with certainty, leading to a model error in the identification process. Here the effect of model error on the performance of the AFTHB approach to identify nonlinear forces is studied by assuming the incorrect forms of nonlinear forces. As Model 1, an incorrect form is assumed for the nonlinear damping between m 4 and m 5 . Hence, a polynomial function with odd terms is considered for the identification, that its Fig. 6. Phase diagram and spectrum of the steady-state dynamic response for x 1 at excitation frequency ω = 24 rad/s and force level F = 20 N.
where f N d is the nonlinear force that replaces the nonlinear damping term c N . c nlq are the unknown coefficients of the equivalent nonlinear force to be identified. Thus, the vector of unknown parameters of the nonlinear force and nonlinear functions used for the identification of Model 1 is now with the corresponding definition of G. In Model 2, parameters of the system are identified considering a modelling error so that the grounded cubic stiffness attached to the oscillator m 4 is neglected. Therefore, the vector of unknown parameters for Model 2 is now The AFTHB approach is now applied using the simulated response of the original nonlinear functions and the models defined by Eqs. (31) and (32) for the identification. Table 4 gives the identified parameters for each model. The results show that in the case of Model 1, the stiffness parameters are estimated accurately. In addition, for the case of Model 2, stiffness parameters k N1 and k N2 are identified accurately, while the identified value of the nonlinear damping c N shows significant error with respect to its true value. To verify the accuracy of the results of the identification process, the parameter set of Table 4 is used to regenerate the simulated response of the system. Fig. 7 compares the phase diagrams and the spectra of the steady state response obtained using the original and equivalent models of the system for x 4 . The results of Fig. 7(a) and 7(b) show very good compatibility between the simulated reconstructed response using the identified parameters of Model 1. On the other hand, Fig. 7(c) and 7(d) illustrate that the modelling error due to neglecting the cubic stiffness k N3 led to significant error in the results of identification.

Case II: System with Coulomb friction
The system of Fig. 5 is now considered with the dry friction between mass m 1 and ground included. Hence the nonlinear force is given by Eq. (26). The same simulation procedure as Case I, with the identical frequency range and step and force levels, is performed to obtain the steady-state dynamics. To examine the performance of the proposed AFTHB identification approach in the presence of measurement noise, the simulated response is polluted with 2% and 5% normally distributed noise. The clean and noisy spectra of the steady state dynamic response of x 1 are shown in Fig. 8 for excitation frequency ω = 23.4rad/s and force level F = 20N. Including Coulomb friction in the dynamics of the system has increased the number of harmonics contributing to the response of the system. The measurement noise has produced a noise floor that limits the number of harmonics that are useful for the identification.
Applying the AFTHB approach to both clean and noisy responses of the system, the unknown parameters are identified utilizing the first five odd harmonics (1,3,5,7,9). Table 5 compares the true values of the parameters and their identified values for both data sets. The performance of the proposed method in identifying the nonlinear systems is excellent even in the presence of noise in the measured/simulated data. Fig. 9(a) shows the experimental set-up of a steel cantilever beam of length 300 mm, width 28.5 mm, thickness 1.45 mm and excited by a shaker. Two permanent magnets are attached to the tip of the beam, and two electromagnets apply a nonlinear electromagnetic force. Three accelerometers measure the response of the beam, and a force transducer measures the excitation from the shaker. The equation of motion of the cantilever beam is derived using Finite Element (FE) analysis with six two-noded beam elements and two degrees of freedom (DOFs) per node, as shown in Fig. 9(b). Point masses are included to represent the accelerometers (nominally 8 g), the force transducer (nominally 26 g) and the permanent magnets attached to the beam (nominally 4 g each), and the shaker attachment is modelled as grounded stiffness k sh and damping c sh . The resulting equations of motion have the form of Eq. (1) with 12 degrees of freedom.

Cantilever beam with electromagnetic actuator
The nonlinear electromagnetic force is assumed to act at the tip of the beam at DOF 11. Five candidate models for the unknown nonlinear electromagnetic force will be identified, as given in Table 6, and include linear stiffness and damping and quadratic, cubic, and fifth-order stiffness. Taghipour et al. [22] showed that for a perfectly symmetric assembly, the nonlinear electromagnetic restoring Table 3 Simulated and estimated parameters for the MDOF system without friction. force only includes odd terms, and hence the first four models are composed only of odd functions. However, it is impossible to enforce exact symmetry in a real experiment, and therefore a nonlinear quadratic stiffness is included in Model 5, although the effect of the quadratic nonlinearity is expected to be low. The vector of unknown parameters a and nonlinear functions G(u,u) for Model 5 are where u 11 and u 11 are transverse deflection and velocity at DOF 11 (the location of electromagnetic force at the tip of the beam). Two vibration tests are performed to identify the nonlinear model of the system. First, a random excitation test is performed to measure the linear response of the system. The underlying linear model is then updated using this measured linear response. Secondly, a stepped-sine vibration test is performed to obtain the steady-state nonlinear response of the structure at various frequencies. The accelerations of the structure are measured using the accelerometers located on the beam according to Fig. 9. The spatially incomplete measured responses of the beam are expanded to the unmeasured coordinates via System Equivalent Reduction Expansion Process (SEREP) method. The displacements of the beam are calculated using the expanded accelerations. The proposed AFTHB approach is then applied to identify the unknown parameters of the nonlinear internal force. To this end, the experimental data within the frequency range 10.6 − 11.4Hz with a frequency step of 0.05Hz is used along with Eq. (12) to estimate the unknown parameters. Fig. 10 shows the experimental results, including linear response obtained from low amplitude random test, the first three harmonics of the nonlinear response obtained from the stepped-sine vibration test, and three harmonics of the applied force used in the stepped-sine test. Only the response at DOF 9 is given, although data from all three accelerometers were used in the identification. In Fig. 10, F 1,1 and W i,H denote, respectively, the amplitude of the primary harmonic of the measured excitation force signal and the amplitude of Table 4 Simulated and estimated parameters for the MDOF system without friction but with modelling error.  H-th harmonic of the response at DOF i. Although the force has significant higher harmonic amplitudes, it is observed that higher harmonics are not significant in the nonlinear response. Therefore, this experimental example addresses the identification of nonlinear system in the presence of higher harmonics of the force. The underlying linear system of the experimental set-up of Fig. 9 is updated using well-known linear model updating techniques [49,50]. The updated linear system is then utilized to identify the nonlinear model of the system using AFTHB approach. The identified parameters of the nonlinear electromagnetic force models are given in Table 7. The optimized parameters of the various models of nonlinear force are then used to regenerate the experimentally measured response of the system. To this end, the measured force signal is applied to the identified models in order to estimate the nonlinear response of the system using various forms of identified nonlinear force given in Table 7.
The nonlinear responses obtained using various identified models are compared with the measured response in Fig. 11. Models 3, 4 and 5 give good predictions of the nonlinear response, although there is still some inaccuracy in the estimated response. This error arises from different sources, such as the expansion method used to estimate the unmeasured degrees of freedom. The quadratic stiffness added to Model 5 to account for the asymmetry of the nonlinear electromagnetic force has a low effect on the steady-state response of the system. This quadratic term does not have a significant effect on the primary harmonic of the response. However, the second and third harmonics of the response given in Fig. 11(b) demonstrates that adding the quadratic term will improve the prediction of the second harmonic of the response.

Cantilever beam with a nonlinear spring
This section considers a stainless-steel cantilever beam with a spring with geometric nonlinearity at the tip, as described by Shaw et al. [51] and shown in Fig. 12. The nonlinearity is implemented as a grounded stiffness attachment at the tip of the beam, as  Table 5 Simulated and estimated parameters for the MDOF system with friction and with and without measurement noise.  Fig. 12(c). A shaker is attached to the beam at x 1 to excite the structure through a force transducer, and three accelerometers are attached to the beam at x 1 ,x 2 , and x 3 along the beam. Further details of the experimental setup and the dimensions and material properties can be found in [51]. The nonlinear restoring force may be modelled based on the geometrical configuration of the nonlinear spring, which may be approximated as a cubic function. Hence, the force applied to the tip of the beam is where w L represents the deflection of the tip of the beam, k denotes the linear stiffness of the springs, l 0 is the free length of the spring,  Table 6 Candidate models of the nonlinear electromagnetic force of the system in Fig. 9. Model 5 fNL 11 = c l1u11 + k l1 u11 + kn 2 u 2 11 + kn 3 u 3 11 + kn 5 u 5 11 and 2d denotes the distance between the fixed ends of the two springs. The coefficients k l and k n represent the coefficients of the cubic approximation. The stinger is modelled as a spring with linear stiffness k s , and the accelerometers and force transducer are modelled as point masses. Proportional damping for the beam is assumed. The beam is modelled using the Galerkin approach based on the lower N m modes of a uniform cantilever beam. The generalised coordinates are the DOFs u in Eq. (1), and w L is a linear function of u. This linear model is updated using the measured natural frequencies. The parameters updated are the stinger stiffness k s , the linear stiffness k l of the nonlinear spring, and the point masses at x 1 , x 2 , and x 3 (denoted m 1 ,m 2 , and m 3 ). Table 8 shows the updated parameters of the underlying linear system of the cantilever beam and compares the experimental and updated natural frequencies. The remaining unknowns are the proportional damping coefficients α, β and the coefficient k n of the cubic stiffness, and these parameters are identified using the AFTHB approach. The residual given by Eq.
(3) is linear in the unknown parameters α, β and k n , where the function g 1 corresponding to k n is g 1 (u) = w 3 L . The experimental steady-state dynamics of the structure was measured in sweep-sine tests in the vicinity of the first natural frequency. Sweep-sine tests were performed using 151 frequency points (df = 0.02Hz) for five different levels of excitation (1.6,2.4,3.2, 3.6, 4.0N) to measure the nonlinear response of the system. Fig. 13 shows the amplitude-frequency diagram of the first and third harmonics of the force signal and the displacement of the beam from the accelerometer near the beam tip. The experimental results of Fig. 13 show the significant contribution of the third harmonic, as well as the primary harmonic, in the dynamic behaviour of the system. Two stable solutions of the system with different amplitudes were measured in response to the same force level using forward (upper branch solution) and backward (lower branch solution) sweep-sine vibration tests.
To apply the AFTHB method to the experimental data, data selection is carried out. To avoid the effect of noisy data on small amplitude response, the experimental data obtained from forward sweep test is used in the identification process. The experimental This identified nonlinear model of the structure was used to regenerate the experimental response of the beam. Fig. 14 compares the  experimental and simulated responses of the system for the excitation level of 4N in the neighbourhood of the first natural frequency. An isolated response of the system was measured for three different excitation levels through controlled forward and backward sweep-sine vibration tests. A sudden excitation force with sufficient amplitude was applied to the structure at ω = 13Hz and the force was controlled at the desired level until the system settled to its steady-state response. Forward and backward sweep-sine tests were then performed within frequency ranges of 13 to 16Hz and 13 to 10Hz, respectively. Detailed explanations of the experiments and measured results are given in [51]. The identified nonlinear model of the experimental test rig is now used to regenerate the isolated nonlinear responses of the beam for an excitation level of 1.6N. Fig. 15 compares the amplitude and phase of the experimental and simulated nonlinear response of the cantilever beam. The results show that the identified nonlinear model of the structure is able to regenerate not only the primary harmonic of the experimental response, but the third harmonic response. The identified model is able to regenerate the isolated response, as well as the non-isolated response, for both primary and third harmonics with good accuracy. A very low level of discrepancy is observed between the estimated and measured data in the isolated response, even though this data was not used in the estimation. Such extrapolation cannot be guaranteed in any estimation method, as it depends on the example and how closely the assumed model represents the real system.

Conclusion
This study proposes an identification method that exploits the multi-harmonic response and force signals of the structure measured from a multi-input/multi-output vibration test in both the time and frequency domains. Two different harmonic balance approaches are used in this method: an Analytical Harmonic-Balance-based (AHB) approach and an Alternating Frequency/Time approach using Harmonic Balance (AFTHB). In both approaches, the underlying linear model of the structure under study is updated using conventional model updating methods prior to starting the nonlinear model identification. In the AHB approach, the proposed method utilizes the expanded nonlinear functions in the frequency domain obtained from analytical methods such as Fourier Integral (FI), Complex Averaging (CXA) technique or Harmonic Balance Method (HBM). In contrast, in the AFTHB approach, the Fourier Transform of the measured time response is calculated to expand the nonlinear functions in the frequency domain. The AHB approach is shown to be very accurate in identifying nonlinear systems if all significant harmonics are considered in the analytical expansion of the nonlinear functions. However, it is often cumbersome to include all significant harmonics of the response in the analytical expansion. Besides, it would be difficult to apply AHB approach to the structures with complex form of nonlinearities like Coulomb friction and also to multidegree of freedom nonlinear systems. The advantages of the proposed AFTHB method are: The AFTHB approach can be applied to all forms of both smooth and non-smooth nonlinear functions. Significant harmonics of the response are included.
Using the AFTHB method results in accurate algebraic equations for each harmonic, including the effect of higher harmonics without truncated error.
The computational cost is relatively low since a dictionary of candidate basis functions is avoided. Although the proposed AFTHB approach is very easy and efficient to use, the following limitations may affect the accuracy of the results of the method: Care should be taken to select the appropriate time step and sampling frequency so that the accuracy of the approach is guaranteed. Indeed, the time step and sampling frequency depend on the frequency range in which the significant harmonics of the response lie. Note that in the AFTHB method the integration time step can be different to the sampling time step for the AFTHB methodusing ode45 in MATLAB will automatically do this.
The method requires prior knowledge or appropriate estimation of the type of system nonlinearities. Hence, engineering insight is necessary to make a proper estimation of the type of nonlinear elements.
The method has been applied to four simulated nonlinear vibration systems with single and multi-degrees of freedom and various types of nonlinearities, including cubic stiffness, Coulomb friction, and quadratic damping, to demonstrate the performance of different approaches. Both methods are able to accurately identify the nonlinear terms in the model if significant higher harmonics of the response are considered as well as the primary harmonic. The results of AHB approach also show that neglecting the higher harmonics leads to inaccurate identification of the nonlinear system, especially when the higher harmonics have significant contribution to the excitation and responses. The application of the proposed method was also demonstrated in two experimental studies of cantilever beam structures with different nonlinear restoring forces at the tip. In the first experimental case study, the successful application of the AFTHB in the identification of different possible nonlinear models when the system is excited by a multi-harmonic force, including the primary, second and third harmonics was shown. In the second experimental case, the AFTHB method has identified the nonlinear stiffness and damping terms of nonlinear model of the system excited by single harmonic force (via controlling the input force) and having significant first and third harmonics in the response. The second case study also indicated the good

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.