# Nonlinear Dynamics of a Duffing-Like Negative Stiffness Oscillator: Modeling and Experimental Characterization

Authors: D. Anastasio; A. Fasana; L. Garibaldi; S. Marchesiello
Journal: Shock and Vibration
 (2020)

---

## Abstract

In this paper, a negative stiffness oscillator is modelled and tested to exploit its nonlinear dynamical characteristics. The oscillator is part of a device designed to improve the current collection quality in railway overhead contact lines, and it acts like an asymmetric double-well Duffing system. Thus, it exhibits two stable equilibrium positions plus an unstable one, and the oscillations can either be bounded around one stable point (small oscillations) or include all the three positions (large oscillations). Depending on the input amplitude, the oscillator can exhibit linear and nonlinear dynamics and chaotic motion as well. Furthermore, its design is asymmetrical, and this plays a key role in its dynamic response, as the two natural frequencies associated with the two stable positions differ from each other. The first purpose of this study is to understand the dynamical behavior of the system in the case of linear and nonlinear oscillations around the two stable points and in the case of large oscillations associated with a chaotic motion. To accomplish this task, the device is mounted on a shaking table and it is driven with several levels of excitations and with both harmonic and random inputs. Finally, the nonlinear coefficients associated with the nonlinearities of the system are identified from the measured data.


---

## Body

## 1. Introduction

Devices and materials exhibiting a negative stiffness phase are often used as vibration isolators due to their amplified damping properties [1, 2]. In particular, in the case of engineering structures, such devices are usually designed adopting discrete macroscopic elements, such as postbuckled beams, plates, shells, and precompressed springs, arranged in appropriate geometrical configurations. Examples can be found in automotive suspensions [3, 4] or seismic isolation [5, 6].While the practical applications of negative stiffness systems are relatively recent, the theoretical studies about their dynamical behavior have a much longer story. This is because of the wide kind of motions they can exhibit, ranging from linear to highly nonlinear and chaotic [7]. In particular, when the negative stiffness effect is coupled to a nonlinear polynomial stiffness contribution, the so-called double-well Duffing oscillator is retrieved. This oscillator exhibits two stable equilibrium positions plus an unstable one, and the oscillations can either be bounded around one stable point (in-well or intrawell, small oscillations) or include all the three positions (cross-well or infrawell, large oscillations). In both cases, periodic oscillations can evolve to steady in-well or cross-well chaotic motions under external excitation [8, 9]. The occurrence of irregular motion, consisting of random-like crossings from oscillations around the two stable equilibrium positions, was first observed experimentally in 1971 [10] and the motion was called “snap-through”. This behavior is associated with an unstable phase generated by the negative stiffness effect; thus, the stability analysis plays a crucial role in the design process of these kinds of systems [11]. Many studies have been conducted starting from the 1970s to analyze the dynamics of this oscillator and the effects of the unstable paths, including phenomena such as bifurcations, subharmonics, period doublings, and chaos. A comprehensive literature review can be found in [9].In this work, a negative stiffness oscillator is studied and tested. The oscillator is part of a device designed to improve the current collection quality in railway overhead contact lines, attempting to alter their damping distribution and reducing the wave propagation [12].A two-fold objective is pursued: first, the dynamical properties of the oscillator are analyzed, replicating experimentally the possible kind of motions it can exhibit (in-well, cross-well, and chaotic); second, nonlinear system identification is performed to extract the model parameters directly from the measurements. The latter in particular seems to be quite a challenging task, given the rich nonlinear dynamics the device is capable of showing and the strength of the nonlinear response. To accomplish these tasks, the oscillator is mounted on a shaking table and it is driven through several levels of excitations with both harmonic and random inputs. The nonlinear system identification is performed by applying the nonlinear subspace identification (NSI) method [13–15]. Generally, NSI estimates a nonlinear state-space model starting from the input-output data, by separating the nonlinear part of the response from the so-called underlying linear system. This is not a straightforward operation for the case considered here though. The bistable nature of the device implies that the underlying linear system is not stable, as it exhibits a negative linear stiffness. Therefore, NSI cannot be directly applied and some prior data manipulation is required. Nonetheless, the identification strategy proposed in this work is capable of estimating both the coefficients of the nonlinearities and the underlying linear systems associated with the stable equilibria starting just from one cross-well measurement, with satisfying accuracy. The nonexistence of a stable underlying linear system and the strong nonlinear behavior make the test set really interesting for the purposes of nonlinear system identification, confirming the effectiveness of the nonlinear subspace identification method. Eventually, the restoring force surface (RFS) method [16, 17] is also applied as a comparison.


## 2. A Negative Stiffness Oscillator

The device here studied consists in a U-shaped steel frame connected through rods to a central moving mass. The frame has the purpose to keep the rods under compression during their movement so that a bistable mechanism is achieved. A schematic representation of the device can be seen in Figure1.Figure 1
Schematic representation of the negative stiffness oscillator.The lower surface of the frame is attached to a shaking table so that a displacementbt can be imposed to the structure. It is also assumed that the inertia of the moving parts can be concentrated into one central point with mass m, comprising the mass of the central bushing and the equivalent inertia of the rods. The vertical movement of this point is described by the coordinate yt, while the rotation of the rods is called θ. The elasticity of the frame is also taken into account, and it acts like a compression force pθ to the rods. Since the flexural elasticity of the frame is much higher than the axial elasticity of the rods, the latter is considered infinitely rigid. A schematic representation of the functional model here described is reported in Figure 2, together with the free-body diagram of the mass m.Figure 2
Model of the negative stiffness oscillator and free-body diagram ofm.The equilibrium equation along the vertical coordinate is(1)my¨t+2pθtsinθt+mg=0.Callingzt=yt−bt, one obtains(2)mz¨t+2pθtsinθt+mg=−mb¨t.The system is a single-degree-of-freedom system; thus, just one variable is needed to describe the motion. In particular, since vertical displacements and accelerations are measured in the experimental setup, the variablezt is taken as an independent variable; therefore, both θ and pθ should be written as a function of z. The elasticity of the frame can be analyzed to obtain the force that the frame transmits to the rods. Just half of the frame is considered in the following, as depicted in Figure 3.Figure 3
Model of the half-frame.When the mass moves along the vertical axis, the frame bends and deforms until a new equilibrium position is reached, corresponding to the resting (undeformed) position of the frame. The analytic expression ofpθ can be found by studying the flexibility of the half-frame, considered as two connected cantilever beams under bending stress. It is not the purpose of this work to analytically derive the expression of pθ. Instead, a qualitative representation of its total vertical component pζθ=2pθsinθ is depicted in Figure 4 as a function of zθ.Figure 4
Qualitative graph of the forcepζθ.It can be seen thatpζz has three roots and crosses the origin with a negative slope. Also, it can be approximated as a polynomial expansion of degree 3:(3)pζz=k3z3+k2z2−k1z.The equation of motion can eventually be written as(4)mz¨+k3z3+k2z2−k1z+mg=−mb¨.Equation (4) has the form of a negative stiffness Duffing equation, whose coefficients k3,k2, and k1 have to be estimated. In particular, the final objective of this work is to identify them directly from the experimental measurements, in order to assure a good correspondence between the numerical model and the real measured behavior. Also, the restoring force Kz and the potential Uz can be defined as(5)Kz=k3z3+k2z2−k1z+mg,(6)Uz=14k3z4+13k2z3−12k1z2+mgz.A qualitative representation of the potential is shown in Figure5, where its double-well nature can be clearly identified. The potential is not symmetric because of the gravitational contribution and the asymmetry of the frame. Also, the three equilibrium positions are depicted, obtained by setting Kz∗=0. Two out of three positions represent a stable equilibrium, namely, z−∗ and z+∗, while the central position z0∗ is an unstable equilibrium point.Figure 5
PotentialUz. Orange dots: equilibrium positions.When the oscillations are bounded around one equilibrium positionz±∗, the system acts like a stable nonlinear oscillator, and the associated linear natural frequency ω± can be computed by(7)ω±=U″z±∗m,with U″z±∗ being the second derivative of Uz computed in z−∗ or z+∗. It is worth writing the equation of motion considering the oscillations around one of the possible equilibrium points. To ease the notation, the generic equilibrium position is called z∗. A new variable xt can be defined as(8)xt=zt−z∗.Substituting equation (8) in equation (4) yields(9)mx¨+Kx+z∗=−mb¨,with(10)Kx+z∗=k3x3+k2+3k3z∗x2+3k3z∗2+2k2z∗−k1x=k3x3+k˜2x2+k˜1x.Finally, equation (9) becomes(11)mx¨+k3x3+k˜2x2+k˜1x=−mb¨.The advantage of using this formulation instead of equation (4) is that it allows the definition of a stable underlying linear system. This is a crucial requirement for the nonlinear subspace identification method, adopted in Section 4 to perform the nonlinear system identification of the structure under test.


## 3. Experimental Characterization

Two photos of the experimental setup corresponding to the two stable equilibrium positions are reported in Figure6. The moving mass is instrumented with an accelerometer to measure its absolute acceleration y¨t and a laser vibrometer to measure its absolute displacement yt. The zero position of yt corresponds to the horizontal configuration of the rods θ=0. The acceleration of the base b¨t is also recorded through a second accelerometer, and the displacement zt is computed as the difference between the laser measure yt and the displacement of the base bt, obtained by integrating twice its measured acceleration.Figure 6
Photos of the experimental setup: (a) negative equilibrium positionz−∗; (b) positive equilibrium position z+∗.
(a)(b)The system is excited with both harmonic and random inputs to characterize its linear and nonlinear behavior.


### 3.1. Random Tests

A first set of measurements is performed with a random input over the frequency range 7–50 Hz. The sampling frequency isfs=512Hz, and the duration of the acquisition is t=300s. Several forcing levels are applied, expressed as RMS values of the acceleration of the base b¨, and the starting position is z−∗. The results are depicted in Figure 7 in terms of time series and transmissibility T, defined as the ratio between the output acceleration z¨ and the input b¨. It can be noted that the mean value of the displacements is not null, as in-well oscillations take place in the neighborhood of the equilibrium position, which in this case is approximately at −0.03 m from the horizontal configuration of the rods θ=0. Also, the asymmetry of the displacements increases as the forcing level, as the system tries to cross the negative stiffness region and to reach the positive equilibrium position. This results in a clear change in the transmissibility, where a softening effect can be seen when increasing the excitation level, in accordance with the theoretical studies that show a similar behavior in the case of in-well motion [9].Figure 7
Random tests. Black line:b¨=7m/s2RMS; orange line: b¨=9m/s2RMS; blue line: b¨=26m/s2RMS. (a) Time history of the displacement (first 60 seconds); (b) transmissibility T.
(a)(b)If the energy given to the system is high enough, cross-well oscillations are retrieved and the moving mass oscillates in a wider region, including the two stable equilibrium positionsz−∗ and z+∗. This situation is represented in Figure 8, where the displacement clearly shows repeated crossings between negative and positive values. Also, the statistical distribution of zt depicted in Figure 8(b) highlights the asymmetry of the structure so that oscillations around the negative position are roughly the 66% of the total acquisition length. This result agrees with the shape of the potential retrieved from the model described in the previous section, which shows two wells with different heights. The final confirmation will be given in the following section where the experimental data are processed so as to identify the actual potential of the system.Figure 8
Random test at the highest amplitude,b¨=38m/s2RMS: (a) time history of the displacement (first 60 seconds); (b) statistical distribution of the displacement.
(a)(b)


### 3.2. Sine-Sweep and Harmonic Tests

Harmonic excitations are a powerful tool to understand the behavior of nonlinear systems. This is particularly true for the case considered here, because of the possibility of chaotic motion. A logarithmic sine-sweep in the frequency range 7–21 Hz is considered first. The sampling frequency isfs=512Hz, and the length of the up-down sweep is t=240s. Two forcing levels are taken into account, expressed as the amplitude b0 of the base displacement, and the starting position is the negative equilibrium z−∗. The up and down sweeps are shown in Figure 9.Figure 9
Sine-sweep tests. Black line:b0=4.5mm; orange line: b0=5mm. (a) Sweep up; (b) sweep down.
(a)(b)Both the excitation levels clearly produce a nonlinear behavior for the in-well motion, as classic jump phenomena can be observed in the up and down sweeps. The softening effect can be seen in Figure9(a), where an increase in the forcing amplitude corresponds to an earlier occurrence of the jump-up. Also, two distinct jumps can be noticed, corresponding to the dominant frequency (10-11 Hz) and its second harmonic (20-21 Hz).If the input is high enough, cross-well motion is retrieved also in this case (Figure9, orange line). It is interesting to look at the harmonic contributions in this case by computing the spectrogram of the relative displacement. The result is reported in Figure 10, where the first two minutes refer to the sweep up, while the second two minutes refer to the sweep down.Figure 10
Spectrogram of the highest level sine-sweep test and corresponding time-domain response.Both even and odd harmonics of the instantaneous frequency are present along the whole acquisition, confirming the asymmetrical behavior of the nonlinear system. Subharmonics are also visible in some regions, in particular around 2 minutes. Generally, they are symptomatic of the possibility of bifurcations and chaotic motion [18, 19], and thus, a series of harmonic tests with constant frequency is performed to analyze these effects. The excitation frequency is ν=9Hz, and three different amplitudes b0 are considered. The results are presented in the phase diagrams in Figure 11.Figure 11
Phase diagrams of the response under harmonic excitation: (a)b0=2mm; (b) b0=4.7mm; (c) b0=5mm.
(a)(b)(c)The different excitation amplitudes result in different kinds of motion of the moving mass, ranging from harmonic oscillations to cross-well motion. In particular, when the amplitude isb0=2mm (Figure 11(a)), the phase plane shows one closed orbit centered around the equilibrium position, i.e., one periodic solution [20]. When the amplitude increases to b0=4.7mm, some nested orbits can be seen in the phase diagram (Figure 11(b)). Two paths in a closed loop are generally representative of the period doubling effect [7] so that a subharmonic with twice the period of the dominant shows up. Such a behavior is retrieved in the experimental case of Figure 11(b), and this is clearer when looking at the power spectral density Z of the displacement, represented in Figure 12. It can be seen that the system responds at both integer multiples of the dominant frequency ν2ν,3ν,… and of its subharmonic 1/2v3/2v,5/2v,….Figure 12
Power spectral density of the harmonic response,b0=4.7mm.As for the highest excitation level in Figure11(c), a cross-well motion occurs, and the solution is not periodic anymore. This behavior lasts for the whole acquisition (10 minutes), and a portion of the time response can be seen in Figure 13.Figure 13
Time response under harmonic excitation,b0=5mm.This kind of persisting nonperiodic response to a periodic excitation is a symptom of chaotic behavior. It should be recalled that no definition of chaos is universally accepted, and this is particularly true when experimental data are considered. The reason is that uncertainties and noise in the data acquisition may interfere with the extreme sensitivity to the initial conditions that characterize chaotic systems. An useful way to check whether a system is behaving chaotically or not is to look at its largest Lyapunov exponentλ [20]. A positive sign of λ means chaotic motion, while a negative sign is a representative of a periodic orbit. Several methods exist to compute λ from experimental time series, and the one proposed in [21] is adopted here. The results are shown in Figure 14, where a positive λ is retrieved.Figure 14
Estimation of the Lyapunov exponent. Red line: convergence mean value.Eventually, the experimental Poincaré sections are computed for different phase synchronizationsϕ of the data with the forcing term [22]. The typical shape of a strange attractor is retrieved [8] and depicted in Figure 15(a) in a polar plot, while one of its sections is represented in Figure 15(b).Figure 15
Experimental Poincaré sections,b0=5mm: (a) polar representation of the attractor surface; (b) Poincaré section, ϕ=0°; (c) Poincaré section, ϕ=170°; (d) Poincaré section, ϕ=320°.
(a)(b)(c)(d)The considered structure is then capable of exhibiting a variety of motions, as expected from the theory. Among them, the cross-well is surely the richest in terms of dynamics, as it covers all the possible positions of the moving mass. For this reason, cross-well measurements will be used in the following section to identify the nonlinear model parameters.


## 3.1. Random Tests

A first set of measurements is performed with a random input over the frequency range 7–50 Hz. The sampling frequency isfs=512Hz, and the duration of the acquisition is t=300s. Several forcing levels are applied, expressed as RMS values of the acceleration of the base b¨, and the starting position is z−∗. The results are depicted in Figure 7 in terms of time series and transmissibility T, defined as the ratio between the output acceleration z¨ and the input b¨. It can be noted that the mean value of the displacements is not null, as in-well oscillations take place in the neighborhood of the equilibrium position, which in this case is approximately at −0.03 m from the horizontal configuration of the rods θ=0. Also, the asymmetry of the displacements increases as the forcing level, as the system tries to cross the negative stiffness region and to reach the positive equilibrium position. This results in a clear change in the transmissibility, where a softening effect can be seen when increasing the excitation level, in accordance with the theoretical studies that show a similar behavior in the case of in-well motion [9].Figure 7
Random tests. Black line:b¨=7m/s2RMS; orange line: b¨=9m/s2RMS; blue line: b¨=26m/s2RMS. (a) Time history of the displacement (first 60 seconds); (b) transmissibility T.
(a)(b)If the energy given to the system is high enough, cross-well oscillations are retrieved and the moving mass oscillates in a wider region, including the two stable equilibrium positionsz−∗ and z+∗. This situation is represented in Figure 8, where the displacement clearly shows repeated crossings between negative and positive values. Also, the statistical distribution of zt depicted in Figure 8(b) highlights the asymmetry of the structure so that oscillations around the negative position are roughly the 66% of the total acquisition length. This result agrees with the shape of the potential retrieved from the model described in the previous section, which shows two wells with different heights. The final confirmation will be given in the following section where the experimental data are processed so as to identify the actual potential of the system.Figure 8
Random test at the highest amplitude,b¨=38m/s2RMS: (a) time history of the displacement (first 60 seconds); (b) statistical distribution of the displacement.
(a)(b)


## 3.2. Sine-Sweep and Harmonic Tests

Harmonic excitations are a powerful tool to understand the behavior of nonlinear systems. This is particularly true for the case considered here, because of the possibility of chaotic motion. A logarithmic sine-sweep in the frequency range 7–21 Hz is considered first. The sampling frequency isfs=512Hz, and the length of the up-down sweep is t=240s. Two forcing levels are taken into account, expressed as the amplitude b0 of the base displacement, and the starting position is the negative equilibrium z−∗. The up and down sweeps are shown in Figure 9.Figure 9
Sine-sweep tests. Black line:b0=4.5mm; orange line: b0=5mm. (a) Sweep up; (b) sweep down.
(a)(b)Both the excitation levels clearly produce a nonlinear behavior for the in-well motion, as classic jump phenomena can be observed in the up and down sweeps. The softening effect can be seen in Figure9(a), where an increase in the forcing amplitude corresponds to an earlier occurrence of the jump-up. Also, two distinct jumps can be noticed, corresponding to the dominant frequency (10-11 Hz) and its second harmonic (20-21 Hz).If the input is high enough, cross-well motion is retrieved also in this case (Figure9, orange line). It is interesting to look at the harmonic contributions in this case by computing the spectrogram of the relative displacement. The result is reported in Figure 10, where the first two minutes refer to the sweep up, while the second two minutes refer to the sweep down.Figure 10
Spectrogram of the highest level sine-sweep test and corresponding time-domain response.Both even and odd harmonics of the instantaneous frequency are present along the whole acquisition, confirming the asymmetrical behavior of the nonlinear system. Subharmonics are also visible in some regions, in particular around 2 minutes. Generally, they are symptomatic of the possibility of bifurcations and chaotic motion [18, 19], and thus, a series of harmonic tests with constant frequency is performed to analyze these effects. The excitation frequency is ν=9Hz, and three different amplitudes b0 are considered. The results are presented in the phase diagrams in Figure 11.Figure 11
Phase diagrams of the response under harmonic excitation: (a)b0=2mm; (b) b0=4.7mm; (c) b0=5mm.
(a)(b)(c)The different excitation amplitudes result in different kinds of motion of the moving mass, ranging from harmonic oscillations to cross-well motion. In particular, when the amplitude isb0=2mm (Figure 11(a)), the phase plane shows one closed orbit centered around the equilibrium position, i.e., one periodic solution [20]. When the amplitude increases to b0=4.7mm, some nested orbits can be seen in the phase diagram (Figure 11(b)). Two paths in a closed loop are generally representative of the period doubling effect [7] so that a subharmonic with twice the period of the dominant shows up. Such a behavior is retrieved in the experimental case of Figure 11(b), and this is clearer when looking at the power spectral density Z of the displacement, represented in Figure 12. It can be seen that the system responds at both integer multiples of the dominant frequency ν2ν,3ν,… and of its subharmonic 1/2v3/2v,5/2v,….Figure 12
Power spectral density of the harmonic response,b0=4.7mm.As for the highest excitation level in Figure11(c), a cross-well motion occurs, and the solution is not periodic anymore. This behavior lasts for the whole acquisition (10 minutes), and a portion of the time response can be seen in Figure 13.Figure 13
Time response under harmonic excitation,b0=5mm.This kind of persisting nonperiodic response to a periodic excitation is a symptom of chaotic behavior. It should be recalled that no definition of chaos is universally accepted, and this is particularly true when experimental data are considered. The reason is that uncertainties and noise in the data acquisition may interfere with the extreme sensitivity to the initial conditions that characterize chaotic systems. An useful way to check whether a system is behaving chaotically or not is to look at its largest Lyapunov exponentλ [20]. A positive sign of λ means chaotic motion, while a negative sign is a representative of a periodic orbit. Several methods exist to compute λ from experimental time series, and the one proposed in [21] is adopted here. The results are shown in Figure 14, where a positive λ is retrieved.Figure 14
Estimation of the Lyapunov exponent. Red line: convergence mean value.Eventually, the experimental Poincaré sections are computed for different phase synchronizationsϕ of the data with the forcing term [22]. The typical shape of a strange attractor is retrieved [8] and depicted in Figure 15(a) in a polar plot, while one of its sections is represented in Figure 15(b).Figure 15
Experimental Poincaré sections,b0=5mm: (a) polar representation of the attractor surface; (b) Poincaré section, ϕ=0°; (c) Poincaré section, ϕ=170°; (d) Poincaré section, ϕ=320°.
(a)(b)(c)(d)The considered structure is then capable of exhibiting a variety of motions, as expected from the theory. Among them, the cross-well is surely the richest in terms of dynamics, as it covers all the possible positions of the moving mass. For this reason, cross-well measurements will be used in the following section to identify the nonlinear model parameters.


## 4. Nonlinear System Identification

A first guess of the nonlinear restoring forceK and the potential U is estimated from the measured time series using the restoring force surface method (RFS) [16, 17]. This method is fairly simple and allows to visualize the nonlinearity, providing a very useful insight to the user. On the contrary, its capabilities are very limited when trying to identify a nonlinear model structure, as it is essentially based on raw data processing and simple fitting. For this reason, the final model is identified using the nonlinear subspace identification method (NSI) [13–15]. NSI gives as outcomes the fully nonlinear state-space representation of the system, together with the FRF of the so-called underlying linear system and the estimation of the coefficients defining the nonlinearities. It requires the input-output data and the knowledge of the nonlinear basis functions, whose selection can take advantage from the preliminary RFS analysis.


### 4.1. First Step: RFS Method

The equation of motion describing the dynamical system considered here can be written in the following form:(12)mz¨+ℛz,z˙=ft,where ft is the forcing term and ℛz,z˙ is the nonlinear restoring surface, generally a function of the displacement z and the velocity z˙. If the inertial term is shifted to the right-hand side of the equation, the restoring surface can then easily be computed and its features extracted. In particular, if small velocities are taken into account, such that z˙<εs, the obtained slice of the restoring surface approximates the restoring force Kz. On the contrary, when small displacements around the equilibrium positions are considered, such that z−z∗<εd, an approximation of the damping force, called Dz˙, is retrieved.For the case considered here, the cross-well measurements are used to buildℛz,z˙. The velocity z˙ is obtained by integrating the accelerations y¨ and b¨ and subtracting them. The obtained restoring surface is reported in Figure 16, together its sections Kz and Dz˙, computed setting εs=εd=0.1%.Figure 16
Experimental restoring surface. Blue dots: restoring force,εs=0.1%; orange dots: damping force, εd=0.1%.The restoring force is eventually fitted to a polynomial expansion of degree 3 as in equation (5), obtaining(13)Kz=k3z3+k2z2−k1z+k0=7.35⋅105z3+1.56⋅103z2−550z+2.4,where k0 is the static component due to gravitational effects. Once the expression of Kz is obtained, the potential Uz can be computed from equation (6). Both Kz and Uz are shown in Figure 17. As for the theoretical model presented in Section 2, the potential shows two wells with different heights and three equilibrium positions are identified, two stable and one unstable. The final values of the coefficients defining the nonlinearity are estimated using NSI, and a comparison with the ones obtained using RFS is also reported in Section 4.2.Figure 17
RFS estimation of the restoring forceKz (a) and the corresponding potential Uz (b). Red dots: stable and unstable equilibrium positions; gray dots: restoring surface for z˙<εs.
(a)(b)As for the damping plot in Figure16, it is rather difficult to estimate a proper damping model due to the very high dispersion of Dz˙. Nevertheless, a possible and realistic guess is that some friction is present between the bushing of the moving mass and the vertical steel guide. For this reason, a nonlinear damping function is also added to the set of nonlinear basis functions given to NSI in the following.


### 4.2. Second Step: NSI Method

The nonlinear subspace identification method is based on a nonlinear state-space representation of the system, obtained by considering the nonlinear terms as feedbacks to the underlying linear system and processing the measured input-output data using the subspace formulation. The existence of the underlying linear system is an essential requirement for the method to work. In the present study, this requirement is not fulfilled when the system goes through cross-well oscillations because of the bistable nature of the device. A possible overcome to this issue has been already shown in Section2, and it consists of a simple shift of the reference axis with respect to one stable equilibrium position z∗, called “reference value” hereafter. It should be noted that the two stable equilibrium positions of the system correspond to two underlying linear systems mutually exclusive, each one existing only when the moving mass oscillates around that equilibrium position. Thus, the two underlying linear systems can be identified with NSI from a single cross-well measurement by selecting the corresponding equilibrium position as the reference value.Whatever reference position is used, a new displacement variablext=zt−z∗ can be defined, as in equation (11). The coefficients of the two nonlinear functions become k3 and k˜2, and nonlinear damping is also considered in the form x˙x˙. Equation (11) can then be rewritten according to the feedback formulation:(14)mx¨+cx˙+k˜1x=−ℛnlx,x˙−mb¨,where a linear viscous damping contribution is also considered and ℛnlx,x˙ is the assumed nonlinear restoring surface, equal to(15)ℛnlx,x˙=k3x3+k˜2x2+cnlx˙|x˙|.Referring to Section3.1, the random test with excitation amplitude 38m/s2RMS (Figure 8) is considered for the nonlinear system identification with NSI. In particular, the last 60 seconds are used as a validation set for the evaluation of the residuals with the measured output, while the identification is performed on the remaining part of the acquisition. The inertance of the system G is depicted in Figure 18, together with the coherence estimation. As expected, the FRF is very distorted due to the strong nonlinear contributions, and the coherence drops to very low values in the frequency region up to the resonance peak.Figure 18
Inertance of the system, random input with amplitude38m/s2RMS. Blue dotted line: coherence.NSI is performed considering first the negative equilibrium position as a reference value, meaningz∗=z−∗. Stability is checked varying the model order for frequencies, damping ratios, MACs, and modal masses [15], and the stabilization diagram is reported in Figure 19. Since the modal parameters show stability from a model of order 2, this is the selected value.Figure 19
Stabilization diagram for NSI with reference positionz∗=z−∗. Stabilization thresholds for natural frequency, damping ratio, MAC, and modal mass are 1%, 20%, 99.5%, and 20%, respectively. Gray dot: new (not stable) pole; blue plus: pole stable in frequency; red square: pole stable in frequency and MAC (modal assurance criterion); orange circle: pole stable in frequency, MAC, and damping; green cross: pole stable in frequency, MAC, damping, and modal mass.Once the model order has been selected, the fully nonlinear state-space model is retrieved, together with the underlying linear system and the coefficients of the nonlinear basis functions. This is repeated also when the reference value is equal to the positive equilibrium position, meaningz∗=z+∗. The stabilization diagram for this case is not reported since it is similar to the one in Figure 19, and a model order equal to 2 is accomplished also in this case. The identified modal parameters of the two underlying linear systems related to the two reference values are reported in Table 1 in terms of natural frequencies and damping ratios.Table 1
Modal parameters of the two underlying linear systems identified with NSI.
Reference positionNatural frequency (Hz)Damping ratio (%)z−∗11.4111.2z+∗9.1920.3The FRFs of the two underlying linear systems are depicted in Figure20 together with the measured (nonlinear) one, already shown in Figure 18.Figure 20
Inertances of the two underlying linear systems. Gray dots: measured inertance of the nonlinear test; blue line: NSI estimation of the linear inertance associated with the negative equilibriumz−∗; orange line: NSI estimation of the linear inertance associated with the positive equilibrium z+∗.The coefficients of the nonlinear basis functions can also be computed from the nonlinear state-space model, recalling that they are treated by the method as frequency-dependent and complex-valued quantities [13]. Assuming that the true coefficients are real numbers, the imaginary part of the identified counterparts should then be zero and the real part should not depend on the frequency. Since this happens only in complete absence of noise and nonlinear modeling errors, the ratio between real and imaginary parts can be taken as an indicator of the goodness of the nonlinear basis functions choice.Running NSI for the two different reference values means that not only two independent underlying linear systems are retrieved but also two nonlinear state-space models. This result in a double estimation of the coefficients of each nonlinearity. In particular, the estimation of the coefficientsk3 and cnl should be the same when NSI is applied to the two reference positions, as they are both invariant in the equation of motion (equations (14) and (15)). On the contrary, k˜2 depends on the choice of z∗. The identified coefficients are depicted in Figure 21 in their real parts for the two reference positions.Figure 21
Real parts of the identified coefficients as frequency-dependent quantities. Blue line: NSI with reference positionz∗=z−∗; orange line: NSI with reference position z∗=z+∗. (a) Coefficient k3N/m3; (b) coefficient k˜2N/m2; (c) coefficient cnlNs2/m2.
(a)(b)(c)A summary of the identified coefficients is also reported in Table2 for the two reference values. The final value of each coefficient is computed as the spectral mean of its real part over the selected frequency range, and the average ratio between real and imaginary parts Eℜ/ℑ is also shown.Table 2
Identified coefficients with NSI for the two reference positions.
Reference positionCoefficientk3Coefficientk˜2CoefficientcnlValueN/m3Eℜ/ℑValueN/m2Eℜ/ℑValueNs2/m2Eℜ/ℑz−∗7.35⋅10532−6.45⋅10434−4.7649z+∗7.86⋅105295.81⋅10429−4.3937It is worth highlighting that the imaginary part is always much lower than the absolute value of the real part in the selected frequency range, which assesses the goodness of the identification. Also, the real parts ofk3 and cnl computed with respect to the two reference values are very close to each other, as expected. Since the system has a “favorite” equilibrium position, which is the negative one, the number of samples associated with oscillations around this position are higher than the other. For this reason, the estimation of the “negative” underlying linear system is considered more reliable and so is the identification of the nonlinearity. The chosen final values of the coefficients are then the ones related to the negative reference position.These coefficients can then be compared with the ones estimated by the RFS method, by looking at the identified restoring force. The comparison is reported in Figure22, where the damping force estimation is also shown.Figure 22
Estimation of the restoring forceKz (a) and the damping force Dz˙ (b). Gray dots: sliced restoring surface; orange line: RFS estimation (only in a); blue line: NSI estimation.
(a)(b)The agreement between RFS and NSI estimation ofKz is very high so that the two curves in Figure 22(a) appear to be almost overlapped, and the percentage deviation is around 1%. As for the damping force Dz˙, it is difficult to evaluate the goodness of the estimation due to the high dispersion of the points. However, at least for the linear part (for small z˙), it is possible to see that the slope of the NSI estimation and the data points agree.Finally, the identified nonlinear state-space model can be used for reproducing the output time series, given the input force. Calling the simulated outputxNSI, a comparison with the measured one can be carried out to validate once more the identification. The validation set corresponds to the last 60 seconds of the acquisition, and the results are reported in Figure 23 as simulated time history against its residual with the measured one in time and frequency domains.Figure 23
Validation of the nonlinear identification in the time domain (a) and frequency domain (b). Black line: measured output; blue line: simulated output; red dots: residual with the measured output spectrum in dB scales (ref. 1 m2/Hz).
(a)(b)The identified state-space model is capable of catching the cross-well motion with a very good accuracy, providing a percentage RMS deviation from the measurement of approximately 8%.


## 4.1. First Step: RFS Method

The equation of motion describing the dynamical system considered here can be written in the following form:(12)mz¨+ℛz,z˙=ft,where ft is the forcing term and ℛz,z˙ is the nonlinear restoring surface, generally a function of the displacement z and the velocity z˙. If the inertial term is shifted to the right-hand side of the equation, the restoring surface can then easily be computed and its features extracted. In particular, if small velocities are taken into account, such that z˙<εs, the obtained slice of the restoring surface approximates the restoring force Kz. On the contrary, when small displacements around the equilibrium positions are considered, such that z−z∗<εd, an approximation of the damping force, called Dz˙, is retrieved.For the case considered here, the cross-well measurements are used to buildℛz,z˙. The velocity z˙ is obtained by integrating the accelerations y¨ and b¨ and subtracting them. The obtained restoring surface is reported in Figure 16, together its sections Kz and Dz˙, computed setting εs=εd=0.1%.Figure 16
Experimental restoring surface. Blue dots: restoring force,εs=0.1%; orange dots: damping force, εd=0.1%.The restoring force is eventually fitted to a polynomial expansion of degree 3 as in equation (5), obtaining(13)Kz=k3z3+k2z2−k1z+k0=7.35⋅105z3+1.56⋅103z2−550z+2.4,where k0 is the static component due to gravitational effects. Once the expression of Kz is obtained, the potential Uz can be computed from equation (6). Both Kz and Uz are shown in Figure 17. As for the theoretical model presented in Section 2, the potential shows two wells with different heights and three equilibrium positions are identified, two stable and one unstable. The final values of the coefficients defining the nonlinearity are estimated using NSI, and a comparison with the ones obtained using RFS is also reported in Section 4.2.Figure 17
RFS estimation of the restoring forceKz (a) and the corresponding potential Uz (b). Red dots: stable and unstable equilibrium positions; gray dots: restoring surface for z˙<εs.
(a)(b)As for the damping plot in Figure16, it is rather difficult to estimate a proper damping model due to the very high dispersion of Dz˙. Nevertheless, a possible and realistic guess is that some friction is present between the bushing of the moving mass and the vertical steel guide. For this reason, a nonlinear damping function is also added to the set of nonlinear basis functions given to NSI in the following.


## 4.2. Second Step: NSI Method

The nonlinear subspace identification method is based on a nonlinear state-space representation of the system, obtained by considering the nonlinear terms as feedbacks to the underlying linear system and processing the measured input-output data using the subspace formulation. The existence of the underlying linear system is an essential requirement for the method to work. In the present study, this requirement is not fulfilled when the system goes through cross-well oscillations because of the bistable nature of the device. A possible overcome to this issue has been already shown in Section2, and it consists of a simple shift of the reference axis with respect to one stable equilibrium position z∗, called “reference value” hereafter. It should be noted that the two stable equilibrium positions of the system correspond to two underlying linear systems mutually exclusive, each one existing only when the moving mass oscillates around that equilibrium position. Thus, the two underlying linear systems can be identified with NSI from a single cross-well measurement by selecting the corresponding equilibrium position as the reference value.Whatever reference position is used, a new displacement variablext=zt−z∗ can be defined, as in equation (11). The coefficients of the two nonlinear functions become k3 and k˜2, and nonlinear damping is also considered in the form x˙x˙. Equation (11) can then be rewritten according to the feedback formulation:(14)mx¨+cx˙+k˜1x=−ℛnlx,x˙−mb¨,where a linear viscous damping contribution is also considered and ℛnlx,x˙ is the assumed nonlinear restoring surface, equal to(15)ℛnlx,x˙=k3x3+k˜2x2+cnlx˙|x˙|.Referring to Section3.1, the random test with excitation amplitude 38m/s2RMS (Figure 8) is considered for the nonlinear system identification with NSI. In particular, the last 60 seconds are used as a validation set for the evaluation of the residuals with the measured output, while the identification is performed on the remaining part of the acquisition. The inertance of the system G is depicted in Figure 18, together with the coherence estimation. As expected, the FRF is very distorted due to the strong nonlinear contributions, and the coherence drops to very low values in the frequency region up to the resonance peak.Figure 18
Inertance of the system, random input with amplitude38m/s2RMS. Blue dotted line: coherence.NSI is performed considering first the negative equilibrium position as a reference value, meaningz∗=z−∗. Stability is checked varying the model order for frequencies, damping ratios, MACs, and modal masses [15], and the stabilization diagram is reported in Figure 19. Since the modal parameters show stability from a model of order 2, this is the selected value.Figure 19
Stabilization diagram for NSI with reference positionz∗=z−∗. Stabilization thresholds for natural frequency, damping ratio, MAC, and modal mass are 1%, 20%, 99.5%, and 20%, respectively. Gray dot: new (not stable) pole; blue plus: pole stable in frequency; red square: pole stable in frequency and MAC (modal assurance criterion); orange circle: pole stable in frequency, MAC, and damping; green cross: pole stable in frequency, MAC, damping, and modal mass.Once the model order has been selected, the fully nonlinear state-space model is retrieved, together with the underlying linear system and the coefficients of the nonlinear basis functions. This is repeated also when the reference value is equal to the positive equilibrium position, meaningz∗=z+∗. The stabilization diagram for this case is not reported since it is similar to the one in Figure 19, and a model order equal to 2 is accomplished also in this case. The identified modal parameters of the two underlying linear systems related to the two reference values are reported in Table 1 in terms of natural frequencies and damping ratios.Table 1
Modal parameters of the two underlying linear systems identified with NSI.
Reference positionNatural frequency (Hz)Damping ratio (%)z−∗11.4111.2z+∗9.1920.3The FRFs of the two underlying linear systems are depicted in Figure20 together with the measured (nonlinear) one, already shown in Figure 18.Figure 20
Inertances of the two underlying linear systems. Gray dots: measured inertance of the nonlinear test; blue line: NSI estimation of the linear inertance associated with the negative equilibriumz−∗; orange line: NSI estimation of the linear inertance associated with the positive equilibrium z+∗.The coefficients of the nonlinear basis functions can also be computed from the nonlinear state-space model, recalling that they are treated by the method as frequency-dependent and complex-valued quantities [13]. Assuming that the true coefficients are real numbers, the imaginary part of the identified counterparts should then be zero and the real part should not depend on the frequency. Since this happens only in complete absence of noise and nonlinear modeling errors, the ratio between real and imaginary parts can be taken as an indicator of the goodness of the nonlinear basis functions choice.Running NSI for the two different reference values means that not only two independent underlying linear systems are retrieved but also two nonlinear state-space models. This result in a double estimation of the coefficients of each nonlinearity. In particular, the estimation of the coefficientsk3 and cnl should be the same when NSI is applied to the two reference positions, as they are both invariant in the equation of motion (equations (14) and (15)). On the contrary, k˜2 depends on the choice of z∗. The identified coefficients are depicted in Figure 21 in their real parts for the two reference positions.Figure 21
Real parts of the identified coefficients as frequency-dependent quantities. Blue line: NSI with reference positionz∗=z−∗; orange line: NSI with reference position z∗=z+∗. (a) Coefficient k3N/m3; (b) coefficient k˜2N/m2; (c) coefficient cnlNs2/m2.
(a)(b)(c)A summary of the identified coefficients is also reported in Table2 for the two reference values. The final value of each coefficient is computed as the spectral mean of its real part over the selected frequency range, and the average ratio between real and imaginary parts Eℜ/ℑ is also shown.Table 2
Identified coefficients with NSI for the two reference positions.
Reference positionCoefficientk3Coefficientk˜2CoefficientcnlValueN/m3Eℜ/ℑValueN/m2Eℜ/ℑValueNs2/m2Eℜ/ℑz−∗7.35⋅10532−6.45⋅10434−4.7649z+∗7.86⋅105295.81⋅10429−4.3937It is worth highlighting that the imaginary part is always much lower than the absolute value of the real part in the selected frequency range, which assesses the goodness of the identification. Also, the real parts ofk3 and cnl computed with respect to the two reference values are very close to each other, as expected. Since the system has a “favorite” equilibrium position, which is the negative one, the number of samples associated with oscillations around this position are higher than the other. For this reason, the estimation of the “negative” underlying linear system is considered more reliable and so is the identification of the nonlinearity. The chosen final values of the coefficients are then the ones related to the negative reference position.These coefficients can then be compared with the ones estimated by the RFS method, by looking at the identified restoring force. The comparison is reported in Figure22, where the damping force estimation is also shown.Figure 22
Estimation of the restoring forceKz (a) and the damping force Dz˙ (b). Gray dots: sliced restoring surface; orange line: RFS estimation (only in a); blue line: NSI estimation.
(a)(b)The agreement between RFS and NSI estimation ofKz is very high so that the two curves in Figure 22(a) appear to be almost overlapped, and the percentage deviation is around 1%. As for the damping force Dz˙, it is difficult to evaluate the goodness of the estimation due to the high dispersion of the points. However, at least for the linear part (for small z˙), it is possible to see that the slope of the NSI estimation and the data points agree.Finally, the identified nonlinear state-space model can be used for reproducing the output time series, given the input force. Calling the simulated outputxNSI, a comparison with the measured one can be carried out to validate once more the identification. The validation set corresponds to the last 60 seconds of the acquisition, and the results are reported in Figure 23 as simulated time history against its residual with the measured one in time and frequency domains.Figure 23
Validation of the nonlinear identification in the time domain (a) and frequency domain (b). Black line: measured output; blue line: simulated output; red dots: residual with the measured output spectrum in dB scales (ref. 1 m2/Hz).
(a)(b)The identified state-space model is capable of catching the cross-well motion with a very good accuracy, providing a percentage RMS deviation from the measurement of approximately 8%.


## 5. Conclusions

A negative stiffness oscillator has been characterized experimentally to exploit its dynamical behavior. A variety of different kinds of motions can be obtained because of the bistable nature of the device, from in-well to cross-well oscillations, including chaotic motion. These dynamical behaviors have been confirmed by the experimental observations, retrieving linear, nonlinear, and chaotic motions. Eventually, the parameters defining the nonlinearity have been recognized via nonlinear system identification, adopting two different methods. A first guess has been obtained using the restoring force surface method, whose implementation allows to easily visualize the nonlinear behavior and the asymmetric double-well potential of the system. The final identification has been performed using the nonlinear subspace identification (NSI) method with a cross-well random measurement. NSI proved to be a robust technique, capable of handling a very strong nonlinear behavior. The nonlinear coefficients have been eventually estimated, together with the linear modal parameters associated with small oscillations around the two equilibrium positions.



---

*Source: 3593018-2020-05-15.xml*