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

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

The lower surface of the frame is attached to a shaking table so that a displacement 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 , 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 , 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 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 .

The equilibrium equation along the vertical coordinate is

Calling , one obtains

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 variable is taken as an independent variable; therefore, both and should be written as a function of . 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.

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 of 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 . Instead, a qualitative representation of its total vertical component is depicted in Figure 4 as a function of .

It can be seen that has three roots and crosses the origin with a negative slope. Also, it can be approximated as a polynomial expansion of degree 3:

The equation of motion can eventually be written as

Equation (4) has the form of a negative stiffness Duffing equation, whose coefficients , and 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 and the potential can be defined as

A qualitative representation of the potential is shown in Figure 5, 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 . Two out of three positions represent a stable equilibrium, namely, and , while the central position is an unstable equilibrium point.

When the oscillations are bounded around one equilibrium position , the system acts like a stable nonlinear oscillator, and the associated linear natural frequency can be computed bywith being the second derivative of computed in or . 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 . A new variable can be defined as

Substituting equation (8) in equation (4) yieldswith

Finally, equation (9) becomes

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 Figure 6. The moving mass is instrumented with an accelerometer to measure its absolute acceleration and a laser vibrometer to measure its absolute displacement . The zero position of corresponds to the horizontal configuration of the rods . The acceleration of the base is also recorded through a second accelerometer, and the displacement is computed as the difference between the laser measure and the displacement of the base , obtained by integrating twice its measured acceleration.

**(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 is , and the duration of the acquisition is . Several forcing levels are applied, expressed as RMS values of the acceleration of the base , and the starting position is . The results are depicted in Figure 7 in terms of time series and transmissibility , defined as the ratio between the output acceleration and the input . 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 . 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].

**(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 positions and . This situation is represented in Figure 8, where the displacement clearly shows repeated crossings between negative and positive values. Also, the statistical distribution of 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.

**(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 is , and the length of the up-down sweep is . Two forcing levels are taken into account, expressed as the amplitude of the base displacement, and the starting position is the negative equilibrium . The up and down sweeps are shown in Figure 9.

**(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 Figure 9(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 (Figure 9, 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.

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 , and three different amplitudes are considered. The results are presented in the phase diagrams in Figure 11.

**(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 is (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 , 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 and of its subharmonic .

As for the highest excitation level in Figure 11(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.

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.

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

**(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 force and the potential 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:where is the forcing term and is the nonlinear restoring surface, generally a function of the displacement and the velocity . 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 , the obtained slice of the restoring surface approximates the restoring force . On the contrary, when small displacements around the equilibrium positions are considered, such that , an approximation of the damping force, called , is retrieved.

For the case considered here, the cross-well measurements are used to build . The velocity is obtained by integrating the accelerations and and subtracting them. The obtained restoring surface is reported in Figure 16, together its sections and , computed setting .

The restoring force is eventually fitted to a polynomial expansion of degree 3 as in equation (5), obtainingwhere is the static component due to gravitational effects. Once the expression of is obtained, the potential can be computed from equation (6). Both and 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.

**(a)**

**(b)**

As for the damping plot in Figure 16, it is rather difficult to estimate a proper damping model due to the very high dispersion of . 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 Section 2, and it consists of a simple shift of the reference axis with respect to one stable equilibrium position , 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 variable can be defined, as in equation (11). The coefficients of the two nonlinear functions become and , and nonlinear damping is also considered in the form . Equation (11) can then be rewritten according to the feedback formulation:where a linear viscous damping contribution is also considered and is the assumed nonlinear restoring surface, equal to

Referring to Section 3.1, the random test with excitation amplitude (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 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.

NSI is performed considering first the negative equilibrium position as a reference value, meaning . 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.

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, meaning . 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.

The FRFs of the two underlying linear systems are depicted in Figure 20 together with the measured (nonlinear) one, already shown in Figure 18.

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 coefficients and 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, depends on the choice of . The identified coefficients are depicted in Figure 21 in their real parts for the two reference positions.

**(a)**

**(b)**

**(c)**

A summary of the identified coefficients is also reported in Table 2 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 is also shown.

It 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 of and 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 Figure 22, where the damping force estimation is also shown.

**(a)**

**(b)**

The agreement between RFS and NSI estimation of 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 , 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 ), 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 output , 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.

**(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.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.