An Optimization-Based Framework for Nonlinear Model Selection and Identification

This paper proposes an optimization-based framework to determine the type of nonlinear model present and identify its parameters. The objective in this optimization problem is to identify the parameters of a nonlinear model by minimizing the differences between the experimental and analytical responses at the measured coordinates of the nonlinear structure. The application of the method is demonstrated on a clamped beam subjected to a nonlinear electromagnetic force. In the proposed method, the assumption is that the form of nonlinear force is not known. For this reason, one may assume that any nonlinear force can be described using a Taylor series expansion. In this paper, four different possible nonlinear forms are assumed to model the electromagnetic force. The parameters of these four nonlinear models are identified from experimental data obtained from a series of stepped-sine vibration tests with constant acceleration base excitation. It is found that a nonlinear model consisting of linear damping and linear, quadratic, cubic, and fifth order stiffness provides excellent agreement between the predicted responses and the corresponding measured responses. It is also shown that adding a quadratic nonlinear damping does not lead to a significant improvement in the results.


Introduction
Mathematical models have been increasingly used in nonlinear structural dynamics. However, the type of nonlinearity and the parameters of these nonlinear models need to be identified from experimental data. This is mainly due to the lack of knowledge about the mechanism of the nonlinearity in structures while in service. Therefore, there has been considerable interest in nonlinear model identification using vibration test data. Thus the identified mathematical model should not only be capable of predicting the real-life behavior of an unknown structure, but also must be physically meaningful. In other words, it is essential to ensure that any mathematical model developed through the process of identification is a valid physical model. Having a valid physical model in the identification process using experimental data will lead to more reliable and meaningful identified parameters.
One deficiency in using the conventional methods, such as H 1 and H 2 , to estimate the frequency response functions of nonlinear systems is that such methods are not able to eliminate the effect of nonlinearities in the system and subsequently predict the underlying linear system accurately. To solve this problem, Richards and Singh [11] presented a spectral identification method based on the "reverse path" formulation to identify multi-degree of freedom nonlinear systems using multi-input multi-output data. For this purpose, broad-band Gaussian random excitation is applied to obtain the nonlinear response of the system. In this approach, the nonlinear elements are considered as feedback forces on the underlying linear system. The application of this method can be found in References [12][13][14][15]26]. Kerschen and Golinval [14] utilized the aforementioned method to introduce a two-step identification approach to generate accurate finite element models of nonlinear structures. In their proposed approach, frequency response functions of the underlying linear system resulting from the Conditioned Reverse Path (CRP) method are correlated with the analytical FRFs to update the linear parameters of the system. CRP was used in [12] to identify and investigate the behavior of a cantilever beam with a geometric nonlinearity. Lenaerts et al. [10] used the updating technique of proper orthogonal decomposition to study the test system of the companion paper [12].
Noël and Kerschen [16] introduced a frequency-domain method to identify nonlinear systems. In this subspace-based technique, the nonlinearity of the system is considered as a feedback force applied to the underlying linear system and is identified utilizing the frequency domain data. Noël et al. [23] investigated a nonlinear spacecraft using two subspace methods in both time and frequency domains. Using these methods, they exploited periodic random measurements to identify the modal parameters of the underlying linear system and the parameters of nonlinear elements of the system.
Hill et al. [17] studied the behavior of a two-degree-of-freedom system subjected to an external force utilizing the backbone curves of the system. They presented a method to find the backbone curves of the system using the second order normal form technique. Renson et al. [18] developed a method to extract the backbone curve of the underlying conservative system. They adapted control-based continuation to track the locus where the response and the excitation are in quadrature. Londoño et al. [19] introduced a method to identify systems with discrete nonlinearities by exciting the structure at a single resonant frequency making the structure vibrate with large displacement where significant effects of nonlinearities are assured. In this method, assuming a form of nonlinearity, the system backbone curves are used to identify the linear and nonlinear parameters of the system. The backbone curves are estimated using measured decay data.
Having measured input data is an important part of conventional identification methods. In order to ignore the necessity of input data measurement, Haroon et al. [24] took advantage of both time and frequency domain techniques to present an approach for identifying the nonlinear elements of structures in the absence of measured input data. Time domain data is used to find a proper assumption for the form of nonlinearity, and the assumed nonlinear form is used with frequency domain data to find an appropriate model to predict the linear and nonlinear response of the system. This paper proposes a framework using an optimization-based method for nonlinear system identification. The method assumes that the nonlinear force can be expanded using a Taylor series. This mathematical model consists of a linear part (first order term) that can be identified from well-known sensitivity-based model updating [2][3][4] and, for the nonlinear part, the order of the nonlinearity increases sequentially, and optimization-based method can be used to identify the parameters for each case. The proposed method is applied to a cantilever beam subjected to an electromagnetic force. The whole system is base excited, and thus a cantilever beam is clamped to the movable bed of a shaker. Two permanent magnets attached to both sides of the tip of the beam along with two electromagnets located symmetrically on the two sides of the permanent magnets generate the nonlinear restoring force applied on the beam. The structure is excited by constant acceleration base motion of the shaker. First, the underlying linear system is updated using the linear response of the system obtained from low amplitude excitation. Afterwards, the nonlinear frequency response of the system measured from a high amplitude vibration test is used to estimate the coefficients of the possible models for nonlinear electromagnetic force.

Experimental Set-Up
A cantilever steel beam of length 30 cm, width 3 cm, thickness 1.5 mm, density 7800 kg/m 3 , and modulus of elasticity 193 GPa, is considered in this study. The beam is clamped to the movable shaker bed and is subject to constant acceleration base excitation using a stepped-sine test. A nonlinear electromagnetic restoring force is applied to the tip of the cantilever beam. This nonlinear force is generated by a pair of permanent magnets attached to the tip of the beam and a pair of electromagnets mounted on a movable shaker bed on both sides of the permanent magnets. The properties of the permanent magnets and the electromagnets are given in Table 1. There is a gap of 6 cm between the permanent and electro magnets at the rest position of the beam. A voltage level of 20 V is applied to the two electromagnets to create attractive magnetic forces. One accelerometer is used to measure the movement of the shaker table. Three accelerometers, labelled AM1, AM2, and AM3, with masses of 8 g each are attached to the beam, respectively, at distances 5, 15, and 25 cm from the clamped end, to measure the translational acceleration of the beam. Figure 1 shows the experimental set-up and its schematic. w b (t) in Figure 1 denotes the base motion of the shaker bed in the direction of the transverse vibration of the beam. acceleration base motion of the shaker. First, the underlying linear system is updated using the linear response of the system obtained from low amplitude excitation. Afterwards, the nonlinear frequency response of the system measured from a high amplitude vibration test is used to estimate the coefficients of the possible models for nonlinear electromagnetic force.

Experimental Set-Up
A cantilever steel beam of length 30 cm, width 3 cm, thickness 1.5 mm, density 7800 kg/m , and modulus of elasticity 193 GPa, is considered in this study. The beam is clamped to the movable shaker bed and is subject to constant acceleration base excitation using a stepped-sine test. A nonlinear electromagnetic restoring force is applied to the tip of the cantilever beam. This nonlinear force is generated by a pair of permanent magnets attached to the tip of the beam and a pair of electromagnets mounted on a movable shaker bed on both sides of the permanent magnets. The properties of the permanent magnets and the electromagnets are given in Table 1. There is a gap of 6 cm between the permanent and electro magnets at the rest position of the beam. A voltage level of 20 V is applied to the two electromagnets to create attractive magnetic forces. One accelerometer is used to measure the movement of the shaker table. Three accelerometers, labelled AM1, AM2, and AM3, with masses of 8 g each are attached to the beam, respectively, at distances 5, 15, and 25 cm from the clamped end, to measure the translational acceleration of the beam. Figure 1 shows the experimental set-up and its schematic. ( ) in Figure 1 denotes the base motion of the shaker bed in the direction of the transverse vibration of the beam. As explained later in this paper, two types of tests are required in the identification process of the present study; first a low amplitude test to measure the linear response of the system is utilized to update the underlying linear model, and second a high amplitude test to measure the nonlinear response of the system is used in the identification process. For the low amplitude test, the structure  As explained later in this paper, two types of tests are required in the identification process of the present study; first a low amplitude test to measure the linear response of the system is utilized to update the underlying linear model, and second a high amplitude test to measure the nonlinear response of the system is used in the identification process. For the low amplitude test, the structure is excited by a constant acceleration excitation of amplitude 0.02 g, where g is the gravitational acceleration. Figure 2 illustrates the linear response of the structure excited by low amplitude base motion with a constant acceleration of 0.02 g. A 1 , A 2 , and A 3 denote the accelerations measured by accelerometers AM1, AM2, and AM3, respectively. In order to obtain a nonlinear response from the high amplitude test, the amplitude of the excitation has to be large enough to guarantee that the nonlinearity of the system is excited. In this study, the base of the structure (shaker bed) is excited by a constant acceleration of amplitude 0.08 g to assure the nonlinear response of the system. 315 response, the electromagnets provide softening nonlinearities. Therefore, a sweep-down steppedsine test is carried out to measure the upper branch of the response. The excitation amplitude should be large enough to guarantee that the nonlinearity of the system is excited. The amplitudes and phases of the experimentally measured nonlinear response of the system at three points along the beam, 5, 15, and 25 cm away from the clamped end of the beam (which are equivalent to the locations of AM1, AM2, and AM3), are shown in Figure 4. The measured nonlinear responses are used to identify the unknown nonlinear restoring force of the system.    [9][10][11] Hz. It is notable that the feedback of the response of the structure, particularly in the neighborhood of the resonant frequency, may change the amplitude of the excitation. Therefore, in order to keep the amplitude of the desired excitation constant, a sine test is carried out at each frequency point to measure the steady state response of the system. As a result, all of the experimental data measured from the sine test at each frequency step are utilized to find the frequency response of the system. To obtain all stable branches of the nonlinear response, both sweep-up and sweep-down tests were carried out. However, as the upper branch (with larger amplitude) is used in the optimization, only measuring the upper stable solution would be satisfactory. For this study, as is explained in Section 3 and also shown in the measured nonlinear response, the electromagnets provide softening nonlinearities. Therefore, a sweep-down stepped-sine test is carried out to measure the upper branch of the response. The excitation amplitude should be large enough to guarantee that the nonlinearity of the system is excited. The amplitudes and phases of the experimentally measured nonlinear response of the system at three points along the beam, 5, 15, and 25 cm away from the clamped end of the beam (which are equivalent to the locations of AM1, AM2, and AM3), are shown in Figure 4. The measured nonlinear responses are used to identify the unknown nonlinear restoring force of the system. satisfactory. For this study, as is explained in Section 3 and also shown in the measured nonlinear response, the electromagnets provide softening nonlinearities. Therefore, a sweep-down steppedsine test is carried out to measure the upper branch of the response. The excitation amplitude should be large enough to guarantee that the nonlinearity of the system is excited. The amplitudes and phases of the experimentally measured nonlinear response of the system at three points along the beam, 5, 15, and 25 cm away from the clamped end of the beam (which are equivalent to the locations of AM1, AM2, and AM3), are shown in Figure 4. The measured nonlinear responses are used to identify the unknown nonlinear restoring force of the system.     In the following section, a mathematical model is derived for the nonlinear system corresponding to the experimental rig.

Mathematical Model
In this section, the nonlinear force obtained from the electromagnets explained in Section 2 is In the following section, a mathematical model is derived for the nonlinear system corresponding to the experimental rig.

Mathematical Model
In this section, the nonlinear force obtained from the electromagnets explained in Section 2 is described by a mathematical model. A cantilever beam of length L, width d, thickness h, density ρ, and modulus of elasticity E is considered. The cantilever is subjected to an unknown nonlinear electromagnetic restoring force f NL at the tip of the beam. The beam is mounted on a shaker base and is subject to base excitation.
The equation of motion of the system is derived using Lagrange's equation. Since the type of nonlinearity is unknown, a mathematical form is assumed for the nonlinear force to be considered in the model. This assumption is based on the physics of the problem and also the measured nonlinear response of the system from experimental vibration tests. The attractive nonlinear force resulting from two symmetric electromagnets is generally given as [27,28], where χ denotes the displacement of the tip of the beam, d 0 is the gap between each electromagnet and the tip at the equilibrium position, and c s is a positive constant coefficient dependant on the voltage of the electromagnets. Using the Taylor expansion about the equilibrium position (χ = 0), the nonlinear electromagnetic force of Equation (1) is rewritten as The negative sign of the force indicates that the attractive forces of two symmetric electromagnets generate softening nonlinearities in the system. Equation (2) may be used as an appropriate assumption for the nonlinear force with coefficients only dependent on the constants c s and d 0 . However, it is often very difficult to have a purely symmetric electromagnetic force in an experiment. Moreover, the model of Equation (1) presented in [27] does not consider the damping effect of the electromagnetic force. Hence, in this study, a different framework will be considered to model the electromagnetic force. As shown in Equation (2), the Taylor series expansion may be used to represent the nonlinear stiffness force created by the electromagnets. Care has been taken to ensure that the configuration is symmetric and therefore we may assume only the odd terms are included. The damping force can be also represented by a Taylor series expansion. Based on this, four possible models to represent the nonlinearity of electromagnetic force are shown in Table 2. As can be seen in the table, the order of stiffness and damping nonlinearity is increased gradually. The level of nonlinearity is increased until it does not improve the identification results, and then the model may be considered good enough to represent the nonlinear electromagnetic force. Table 2. Possible models to represent the electromagnetic force.
Model I: χ + k l χ + k n 2 χ|χ| + k n 3 χ 3 + k n 5 χ 5 (6) In Table 2 λ, λ n denote the linear and nonlinear damping coefficients, and k l , k n 2 , k n 3 , k n 5 indicate the linear stiffness, symmetric quadratic stiffness, cubic and fifth order stiffnesses, respectively. Since the electromagnets are mounted on the movable bed of the shaker, as well as the beam, the displacement χ in Equations (3)-(6) should be the relative displacement of the tip of the beam.
Another important issue is to consider the effect of the aluminum L-shaped clamp to which the cantilever beam is attached. In fact, the beam under investigation is clamped between two short thick aluminum beams. These two aluminum support beams act like two very stiff springs attached to the base of the steel beam. Although the clamp is very stiff and will not undergo large deflections, it is not completely rigid. The high stiffness of the clamp may cause significant changes in the frequency response of the system, particularly in the frequencies of the resonances and anti-resonances of the system. Therefore, the effect of the very high stiffness of the clamp is simulated by a linear spring between the shaker bed and the base of cantilever. This linear spring represents the bending stiffness of the L-shaped clamp. It is worth noting that the deflection of the L-shaped clamp is in the same direction as the motion of the beam and the rotational motion of the clamped support acts in the torsional direction of the beam (as seen in Figure 5). Because the torsional motion of the beam is not considered in this study, there is no need for a rotational spring representing the torsional degree of freedom of the beam. All of the bolted joints are assumed to be rigidly connected in this model. The equivalent masses of two aluminum support beams are also taken into account as an additional mass. Figure 5 illustrates the equivalent system of the beam subject to the motion of shaker bed. The torsional deflection of the L-shaped clamp is assumed to be negligible, and hence the beam does not have any rotation at the base. Hence, the assumed mode shapes of the cantilever are used to simulate the dynamics of the beam. The equivalent mass m b and linear stiffness k b of the aluminum beams can be estimated by the following equations, [29].
where ρ Al and E Al are density and modulus of elasticity of aluminium, respectively. V Alb , I Alb , and L Alb denote, respectively, the volume, the second moment of area, and the length of each aluminium support beam.
Vibration 2019, 2 FOR PEER REVIEW 7 7 thick aluminum beams. These two aluminum support beams act like two very stiff springs attached to the base of the steel beam. Although the clamp is very stiff and will not undergo large deflections, it is not completely rigid. The high stiffness of the clamp may cause significant changes in the frequency response of the system, particularly in the frequencies of the resonances and antiresonances of the system. Therefore, the effect of the very high stiffness of the clamp is simulated by a linear spring between the shaker bed and the base of cantilever. This linear spring represents the bending stiffness of the L-shaped clamp. It is worth noting that the deflection of the L-shaped clamp is in the same direction as the motion of the beam and the rotational motion of the clamped support acts in the torsional direction of the beam (as seen in Figure 5). Because the torsional motion of the beam is not considered in this study, there is no need for a rotational spring representing the torsional degree of freedom of the beam. All of the bolted joints are assumed to be rigidly connected in this model. The equivalent masses of two aluminum support beams are also taken into account as an additional mass. Figure 5 illustrates the equivalent system of the beam subject to the motion of shaker bed. The torsional deflection of the L-shaped clamp is assumed to be negligible, and hence the beam does not have any rotation at the base. Hence, the assumed mode shapes of the cantilever are used to simulate the dynamics of the beam. The equivalent mass and linear stiffness of the aluminum beams can be estimated by the following equations, [29].
where and are density and modulus of elasticity of aluminium, respectively. , , and denote, respectively, the volume, the second moment of area, and the length of each aluminium support beam. Considering the assumed type of the nonlinearity and the effect of the clamp as a linear stiffness attached to the base of the beam, the governing equations of the system are derived using the Euler-Lagrange method. In this case, the equation of motion will be derived for case IV but the obtained solution for the nonlinear frequency response functions may be used for the three simpler Considering the assumed type of the nonlinearity and the effect of the clamp as a linear stiffness attached to the base of the beam, the governing equations of the system are derived using the Euler-Lagrange method. In this case, the equation of motion will be derived for case IV but the obtained solution for the nonlinear frequency response functions may be used for the three simpler models, i.e. Model I (k n 5 = k n 2 = λ n = 0), Model II (k n 2 = λ n = 0) and Model III (λ n = 0). The kinetic and potential energy are given as follows where w(x, t) is the transverse vibration of the beam, x is the distance from the clamped end of the beam, w b (t) is the base motion of the shaker in the transverse direction of the beam, M n denotes all additional masses including accelerometers, permanent magnets attached to the tip, and the equivalent mass of the clamp support, and N p is the total number of additional masses. Defining z(x, t) as the relative displacement of the beam with respect to the base motion w b (t) of the shaker bed, then and the kinetic and potential energy can be rewritten as, Considering the proportional damping for the beam, Rayleigh's dissipation function of the system is also derived as, where α, β are the proportional damping coefficients per length of the beam. Using the Euler-Lagrange equation, along with the assumed mode method, the equation of motion of the system is derived. For this purpose, a solution is assumed for the response of the beam as where N m represents the total number of modes used in the analysis, q i denotes the i−th time variable modal coordinate, and φ i is the i−th mode shape of the beam with the following boundary conditions In this study, a single-mode approximation of the system is utilized for the identification process. Therefore, the equations of motion are derived for this single-mode approximation. For the sake of brevity, the subscripts related to the first mode are neglected. In order to apply the assumed mode method, the quadratic stiffness and damping energy terms of Equations (12) and (13) that include absolute functions can be approximated using following series, where c 1 , c 3 , c 5 , . . . and e 1 , e 3 , e 5 , . . . are the coefficients of the approximations of Equations (16), and T n (z) is the Chebyshev polynomial of the first kind, for example Given the approximations of Equation (16) for the quadratic stiffness and damping terms, substituting the assumed solution of Equation (14) into Equations (11)- (13), and applying the Euler-Lagrange method, the equation of motion is obtained as where q denotes the time variable modal coordinate of first mode, and There are many methods that are able to solve the nonlinear differential equations of motion [28,30]. In this study, the Harmonic Balance Method (HBM) is used to investigate the steady state dynamics of the system.
As the primary harmonic is dominant in the nonlinear response of the system, the higher harmonics are neglected in this paper. Hence, in order to apply the HBM to the equations of motion (Equation (18)) of the system, the periodic time variable response q(t) of the system is assumed to be Vibration 2019, 2, 20 320 of 330 The slowly varying sine and cosine coefficients a(t) and b(t) are considered time-variant to study the transient response of the system. However, in this paper only the steady-state behavior is examined.
Substituting Equations (21) and (22) into Equation (18) and balancing the coefficients of sin(Ωt) and cos(Ωt), two differential equations are obtained as follows. where P 4 = λ n , P 5 = k n 2 , P 6 = k n 3 R 7 , P 7 = k n 5 R 8 , Vibration 2019, 2, 20 321 of 330 and S q3 = 3 4 a 3 + 3 2 ab 2 , where W b is the amplitude of the base acceleration .. w b . In order to find the steady-state response of the system, the time-variant terms in Equation (23) are set to zero. Consequently, the nonlinear differential equation of motion is transformed into two nonlinear algebraic equations as follows where Different numerical methods may be used to solve the nonlinear algebraic equations. Here, in this study, the resultant nonlinear algebraic equations of Equation (26) are solved using the arc-length continuation method. The amplitude-frequency response of the system obtained from the numerical method is used in the optimization process to identify the nonlinear system. The following section explains the optimization process.

Nonlinear Identification
There are different methods introduced to identify the nonlinearities in the system. One of the biggest weaknesses of most of these methods is the accumulative error in the results of the identification process. Errors from various sources are accumulated in the identification process leading to a considerable error in the identified values of the parameters. Using the measured frequency domain data directly in the identification of parameters reduces such errors to the lowest possible. One significant error is when the low amplitude response of a nonlinear mechanical structure is used to update the underlying linear system. Although the effect of nonlinearity on the response of the system is at its lowest level for very low amplitude excitation, this effect cannot be eliminated from the measured response it can affect the results of the updating procedure. Therefore, some modelling error appears in the underlying linear system.
In many practical engineering structures, it is not possible to have spatially complete measurements. However, some of the finite element-based identification methods require all coordinates to be measured [21,22]. Hence, various expansion methods, such as the System Equivalent Reduction Expansion Process (SEREP), may be used to estimate the unmeasured coordinates, which may cause error in the final results of the identification process. Avoiding expansion methods would eliminate the possible error in this process. The approach presented in this study aims to eliminate the aforementioned errors in the process of nonlinear identification.
The proposed method is composed of two main steps. In the first step, a linear model updating method is applied to the linear response of a low amplitude test. The resultant updated parameters are used as an appropriate initial guess for the underlying linear system. The updated parameters of the underlying linear system are then corrected in step two. By correcting the updated parameters obtained from the first step, the above-discussed error due to the effect of nonlinearity in the low amplitude response of the system is avoided.
The second step is to identify the underlying linear system and the nonlinear restoring force utilizing the updated linear model and nonlinear experimental response in an optimization process. To this end, a nonlinear form is assumed for the nonlinearity corresponding to the experimental response. Here, good engineering insight is required to have an appropriate assumed model. Consequently, a numerical model is developed for the nonlinear system using the assumed nonlinear force with estimated linear parameters to be corrected and unknown nonlinear parameters to be estimated. Then, the optimized values of the system parameters are targeted considering an objective function defined by minimizing the difference between the experimental and analytical responses at the measured coordinates, [31], as where X m and X a denote respectively the experimental and numerical responses of the system in the frequency domain at the measured points. Using the logarithmic scale in the objective function penalizes large changes between optimized and experimental responses. Considering only measured coordinates in the identification procedure eliminates the need for a spatially complete measurement and the need to use expansion methods. Therefore, the associated error is eliminated from the calculations. N f is the number of measured frequency lines considered in the optimization. The optimized parameters are subjected to different constraints limiting the solution range of the objective function. For instance, the optimized parameters must have values within the allowable range. Some data selection criteria may be applied to the frequency response before being used in the optimization process as follows: (a) The first data selection is to narrow the frequency response to a range in which the effect of the nonlinearity is considerable. By applying this criterion, the calculation cost is reduced significantly, because there is no significant difference between the linear and the nonlinear responses in the regions which are insensitive to the nonlinearity.
(b) The second is to neglect the unstable branches of the numerical solution which are not measured in the experimental tests.
(c) The third data selection criterion is to choose between the solutions in the multi-response regions. Having more than one solution at a single frequency line may cause the optimization process to slow down or, in some cases, not to converge to a reasonable solution. In this case, usually the upper solution branch (having the larger amplitude) is selected, because the effect of the nonlinearity is more significant compared to the lower branch. Of course, the selected branch must correspond to the experimentally measured response used in the optimization. The lower branch of the response or the responses of different excitation levels can be used to make sure that the results of the optimization are correct.
Applying all the constraints and criteria described above, the parameters of the assumed model of the nonlinear system are optimized. These values are used to estimate the linear (low amplitude test) and nonlinear response (high amplitude test) of the system. If an acceptable match is not observed between the experimentally measured response and the estimated response, the assumed form of the nonlinear force is revised, and the optimization process is repeated. This process is carried out until a good correlation between the experimental results and the analytical model is achieved.

Results and Discussion
In this section, first the results of the linear model updating are presented and compared with the experimental linear response of the system. Afterwards, the results of the identification of the unknown nonlinear force of the system are discussed.
In order to identify the structural model of the cantilever beam, first the underlying linear model is updated using linear updating methods [2][3][4]. For this purpose, the linear response of the nonlinear system is measured by exciting the structure with very low amplitude base excitation. The first three measured natural frequencies are used to update the underlying linear system. Three parameters of the system are updated, namely the flexural stiffness EI, the base stiffness k b , and the base mass m b . After updating these three parameters using the natural frequencies, the proportional damping coefficients α, β are updated. Table 3 gives the values of the fixed and updated parameters of the system. Table 3. Parameters used in the linear model updating.

Parameters
Parameter Values (in SI Units)  Having updated the parameters, the numerical linear frequency response is obtained using the updated linear model of the underlying linear structure. Figure 8 shows the experimental and updated linear frequency response of the system at the location of 1 st accelerometer (AM1) for the frequency range containing the first three modes. A very important observation is that the updated model can estimate the anti-resonance frequencies, which are not included in the updating method. Therefore, the good agreement between both the resonance and anti-resonance frequencies shows the updated linear model is valid. The low-amplitude experimental and updated analytical linear  Having updated the parameters, the numerical linear frequency response is obtained using the updated linear model of the underlying linear structure. Figure 8 shows the experimental and updated linear frequency response of the system at the location of 1 st accelerometer (AM1) for the frequency range containing the first three modes. A very important observation is that the updated model can estimate the anti-resonance frequencies, which are not included in the updating method. Therefore, the good agreement between both the resonance and anti-resonance frequencies shows the updated linear model is valid. The low-amplitude experimental and updated analytical linear Having updated the parameters, the numerical linear frequency response is obtained using the updated linear model of the underlying linear structure. Figure 8 shows the experimental and updated linear frequency response of the system at the location of 1st accelerometer (AM1) for the frequency range containing the first three modes. A very important observation is that the updated model can estimate the anti-resonance frequencies, which are not included in the updating method. Therefore, the good agreement between both the resonance and anti-resonance frequencies shows the updated linear model is valid. The low-amplitude experimental and updated analytical linear responses are compared in Figure 9 in the neighborhood of the first natural frequency in the presence of the nonlinear electromagnetic force.

Description of Parameters Initial Values Updated Values
Vibration 2019, 2 FOR PEER REVIEW 15 15 responses are compared in Figure 9 in the neighborhood of the first natural frequency in the presence of the nonlinear electromagnetic force.  responses are compared in Figure 9 in the neighborhood of the first natural frequency in the presence of the nonlinear electromagnetic force.  Given the updated model of the underlying linear system, the nonlinear system is now identified. For this purpose, the nonlinear amplitude-frequency response of the system at the location of 3rd accelerometer (AM3) is used in optimization process to identify the nonlinear force. As expected from Equation (2), the response of the system demonstrates a softening nonlinearity in the system. Having the updated underlying linear system, the measured nonlinear response of the system is utilized in the optimization process to estimate the unknown parameters of the nonlinear force model I given by Equation (3) as Negative values of k l and k n 3 indicate the softening nonlinearity due to the attractive electromagnetic force. The optimized parameters of the system are used to estimate the nonlinear response of the system, which is then compared to the experimentally measured response. Figure 10 shows the estimated response nearest the beam tip using model I of Equation (3) for the nonlinear force, and demonstrates that the estimated nonlinear model is able to provide a fairly accurate prediction for the response of the nonlinear system. 16 Given the updated model of the underlying linear system, the nonlinear system is now identified. For this purpose, the nonlinear amplitude-frequency response of the system at the location of 3 rd accelerometer (AM3) is used in optimization process to identify the nonlinear force. As expected from Equation (2), the response of the system demonstrates a softening nonlinearity in the system. Having the updated underlying linear system, the measured nonlinear response of the system is utilized in the optimization process to estimate the unknown parameters of the nonlinear force model I given by Equation (3) Negative values of and indicate the softening nonlinearity due to the attractive electromagnetic force. The optimized parameters of the system are used to estimate the nonlinear response of the system, which is then compared to the experimentally measured response. Figure  10 shows the estimated response nearest the beam tip using model I of Equation (3) for the nonlinear force, and demonstrates that the estimated nonlinear model is able to provide a fairly accurate prediction for the response of the nonlinear system. Figure 10. The experimentally measured nonlinear response at the location of AM3 versus the estimated response of the nonlinear system using the first assumed model for the nonlinear force. Now, in order to increase the accuracy of the nonlinear model, Models II, III, and IV, introduced in Table 2, are used in the identification process. The optimization process is repeated for these models of the nonlinear force; the optimized parameter values have been obtained and are presented in Table 4. Subsequently, the optimized parameters are used to estimate the nonlinear response of the system with the same excitation level. Figure 11 compares the measured nonlinear response, including amplitude and phase, at AM3 (25-cm away from the clamped end) with the estimated responses using all four assumed models of the nonlinear force. Figure 11 shows that the predictions from these four models are generally in good agreement with the experimental data. However, it can be seen that there is some improvement in the predictions from model III when compared to model I, while there is very little difference in the predictions between models III and IV. Based on this observation, one may conclude that because no noticeable change in the level of accuracy was achieved when the order of nonlinearity increased from model III to model IV, model III should be considered as the identified model to represent the system. Now, in order to increase the accuracy of the nonlinear model, Models II, III, and IV, introduced in Table 2, are used in the identification process. The optimization process is repeated for these models of the nonlinear force; the optimized parameter values have been obtained and are presented in Table 4. Subsequently, the optimized parameters are used to estimate the nonlinear response of the system with the same excitation level. Figure 11 compares the measured nonlinear response, including amplitude and phase, at AM3 (25-cm away from the clamped end) with the estimated responses using all four assumed models of the nonlinear force. Figure 11 shows that the predictions from these four models are generally in good agreement with the experimental data. However, it can be seen that there is some improvement in the predictions from model III when compared to model I, while there is very little difference in the predictions between models III and IV. Based on this observation, one may conclude that because no noticeable change in the level of accuracy was achieved when the order of nonlinearity increased from model III to model IV, model III should be considered as the identified model to represent the system.    Another important problem in nonlinear system identification is assessing the validity of the model. In this example, the optimization process was performed using the response of only AM3. Therefore, in order to verify the validity of the results, the optimized parameters of model IV of the assumed nonlinear force are used to estimate the response of the nonlinear system at all three measured accelerometers, including both amplitude and phase of the acceleration. The estimated response is compared with the experimental data in Figure 12. The results illustrate a good agreement between the measured and reconstructed response and confirms the validity of the identified model. Another important problem in nonlinear system identification is assessing the validity of the model. In this example, the optimization process was performed using the response of only AM3. Therefore, in order to verify the validity of the results, the optimized parameters of model IV of the assumed nonlinear force are used to estimate the response of the nonlinear system at all three measured accelerometers, including both amplitude and phase of the acceleration. The estimated response is compared with the experimental data in Figure 12. The results illustrate a good agreement between the measured and reconstructed response and confirms the validity of the identified model. Vibration 2019, 2 FOR PEER REVIEW 18 18 Figure 12. Comparison of the experimentally measured nonlinear responses and the numerically estimated responses using nonlinear model IV. , , and denote the accelerations measured by accelerometers AM1, AM2, and AM3, respectively.

Conclusion
This paper has proposed a framework for nonlinear model identification. In this framework, the nonlinear element/force was represented by a Taylor series expansion. Based on the number of terms retained in the expansion, different possible models may be chosen. An optimization-based method can be used to identify the parameters of the chosen models. The objective function of this optimization is the difference between the predicted nonlinear response and the experimental response. The predicted responses from the identified models are then compared with experimental results and one may determine how many terms should be retained in the Taylor series expansion. The proposed method is demonstrated on an experimental test case, where a nonlinear electromagnetic force is applied to the tip of a steel cantilever beam. First, the linear and nonlinear responses of the system were measured from low and high amplitude vibration tests, respectively. Euler Bernoulli beam theory is used to model the cantilever beam. The support of the beam is modelled with a spring and the effect of the mass of the support is also considered in the model. It was found that the linear model of the support is needed to predict the anti-resonances of the linear experimental responses. In the next step, four possible nonlinear models (using Taylor series expansion) were assumed for the electromagnetic force. The order of the nonlinearity in the stiffness and damping terms was increased gradually in these models. The measured nonlinear responses were used in an optimization-based identification process to estimate the unknown parameters of these models. It has been demonstrated that a nonlinear model composed of linear damping and nonlinear stiffness (including linear, quadratic, cubic, and fifth-order terms) provides excellent agreement between the predicted responses and the corresponding measured responses.

Conclusions
This paper has proposed a framework for nonlinear model identification. In this framework, the nonlinear element/force was represented by a Taylor series expansion. Based on the number of terms retained in the expansion, different possible models may be chosen. An optimization-based method can be used to identify the parameters of the chosen models. The objective function of this optimization is the difference between the predicted nonlinear response and the experimental response. The predicted responses from the identified models are then compared with experimental results and one may determine how many terms should be retained in the Taylor series expansion. The proposed method is demonstrated on an experimental test case, where a nonlinear electromagnetic force is applied to the tip of a steel cantilever beam. First, the linear and nonlinear responses of the system were measured from low and high amplitude vibration tests, respectively. Euler Bernoulli beam theory is used to model the cantilever beam. The support of the beam is modelled with a spring and the effect of the mass of the support is also considered in the model. It was found that the linear model of the support is needed to predict the anti-resonances of the linear experimental responses. In the next step, four possible nonlinear models (using Taylor series expansion) were assumed for the electromagnetic force. The order of the nonlinearity in the stiffness and damping terms was increased gradually in these models. The measured nonlinear responses were used in an optimization-based identification process to estimate the unknown parameters of these models. It has been demonstrated that a nonlinear model composed of linear damping and nonlinear stiffness (including linear, quadratic, cubic, and fifth-order terms) provides excellent agreement between the predicted responses and the corresponding measured responses. In addition, it is shown that adding a quadratic nonlinear damping to the assumed nonlinear model may not lead to a significant change in the results. The identified model is capable of predicting the experimental FRFs at points which are not used in the identification. This indicated the identified model is valid.