Analysis of a mass-spring-relay system with periodic forcing

The dynamics of a hysteretic relay oscillator with harmonic forcing is investigated. Periodic excitation of the system results in periodic, quasi-periodic, chaotic and unbounded behavior. An explicit Poincaré map is constructed with an implicit constraint on the switching time. The stability of the fixed points of the Poincaré map corresponding to period-one solutions is investigated. By varying the forcing parameters, we observed a saddle-center and a pitchfork bifurcation of two centers and a saddle-type fixed point. The global dynamics of the system exhibits discontinuity induced bifurcations of the fixed points.

superconductivity, adsorption, desorption, and materials with shape memory [2][3][4][5]. Hysteresis is a wellknown phenomenon in relay control systems [6]. In general, relays have two output branches, and the output of a relay jumps discontinuously whenever the input exceeds a certain critical value. For an ideal relay, there is a single critical value for which the output is discontinuous, while for a relay with hysteresis, there are two such critical input values as shown in Fig.  1a. Hassani et al. [7] surveyed various mathematical models for hysteresis such as Preisach, Krasnosel'skii-Pokrovskii, Prandtl-Ishlinskii, Maxwell-Slip, Bouc-Wen and Duhem in terms of their applications in modeling, control, and identification of dynamical systems.
Vaiana et al. [8,9] recently presented an enhanced formulation of a class of uniaxial phenomenological hysteresis models, to simulate complex symmetric and asymmetric hysteresis phenomena typical of several rate-independent mechanical systems and materials.
Andronov et al. [10] and Åström [11] investigated the existence of periodic solutions having two relay switches in one period. The periodic response of an ideal relay system with sinusoidal excitation was studied in [12]. Caughley [13] focused on a system with a bilinear hysteresis model and sinusoidal excitation, showing a soft-type resonance with stable, singlevalued response curves. The response of a bilinear hysteretic oscillator with sinusoidal excitation in a hybrid systems framework was investigated in [14]. Zhusubaliyev et al. [15] presented the study of a fourdimensional system with hysteresis exhibiting periodic and chaotic oscillations. Di Bernardo and coworkers [16,17] studied the dynamic behavior of systems with ideal and non-ideal relays. The dynamics of delayed relay systems is investigated extensively in [18,19]. Sivaselvan [20] presented several hysteretic models formulated using an energy approach with specifying scalar-valued functions: a stored energy function and a dissipation potential. Lacarbonara and Vestroni [21] investigated the responses and codimension-one bifurcations in Masing-type and Bouc-Wen hysteretic oscillators. A periodically driven damped harmonic oscillator coupled to a random field Ising model showing complex hysteresis is studied in [22].
Semenov et al. [23] studied the resonance properties of an oscillating system in the case when the energy pumping is made by an external source of hysteretic nature. The nonlinear dynamic analysis of nonstiffening hysteretic mechanical systems was carried out by using a novel rate-independent phenomenological model with an explicit time integration in [24]. Krasnosel'skii and Rachinskii [25] investigated the bifurcations of forced periodic oscillations for systems with Preisach hysteresis. Kalmár-Nagy et al. [26] studied a single-degree-of-freedom forced hysteretic system.
Hysteresis can be defined as the dependence of the state of a system on its history. This phenomenon can be modeled in several ways. Two of the most essential hysteresis models are the Preisach model [27] and the Bouc-Wen model [28,29]. The Preisach model of hysteresis was introduced to model magnetic hysteresis as the relationship between magnetic field and magnetization of ferromagnetic materials, but it can be used perfectly for other scientific areas. This model is based on non-ideal relay hysterons; see Fig. 1a. The discrete form of the Preisach model was studied by using properties of the underlying graph in [30]. In this paper, we will use a single symmetric elementary hysteretic relay operator (hysteron). The hysteretic relay operator This operator depends on the time history of the continuous function x(t); i.e., the discrete variable e is either −1 or 1, depending on whether the function x(t) enters the hysteretic region (1, −1) from the left or right (see Fig. 1a).
The power generated by the hysteresis element is . This type of hysteretic element can thus inject energy into the system or dissipate energy from the system. This paper aims to understand the dynamical behavior of forced oscillations of systems containing hysteretic relay elements, thereby continuing the work of Kalmár-Nagy et al. [26]. We investigate the dynamics of the one-degree-of-freedom forced oscillator by utilizing an explicit Poincaré map together with implicit constraints [31].
The outline of the paper is as follows. Section 2 introduces the equations describing a one-degree-offreedom mass-spring-relay. This simple model can describe, for example, a bang-bang controlled springmass system with external forcing. The different types of behavior of the system is demonstrated in Sect. 3. The explicit solution of the forced system is shown in Sect. 4. Section 5 is devoted to the switching time equations of the forced system. We show that the switching time can be a discontinuous function of the initial conditions and the forcing parameters. The Poincaré map of the system is constructed in Sect. 6. The fixed points of the Poincaré map and their stability were determined analytically in Sect. 7. The local dynamics around the fixed points is illustrated by Poincaré sections. The global dynamics of the system is discussed in Sect. 8. Finally, conclusions are drawn in Sect. 9.

Mathematical model
The model of the system is shown in Fig. 1b. The equation of motion of the system is where F [x(t)] is the hysteretic relay operator defined in Eq. (1). The parameters A ≥ 0, ω > 0, and φ ∈ S 1 are the amplitude, frequency, and phase of the forcing, respectively. Without loss of generality, the initial conditions of system (2) can be specified as subsystem I I :ẍ I I (t) + x I I (t) + 1 = A cos(ωt + φ I I ).
The dynamics of system (2) switches from subsystem I to subsystem I I when the solution trajectory reaches x I (t) = 1 from the left, i.e.,ẋ I (t) ≥ 0. The switch/transition from subsystems I I to subsystem I occurs when the solution reaches x I I (t) = −1 from the right, i.e.,ẋ I I (t) ≤ 0.
To simplify the analysis, time is reset to zero when the transition occurs between subsystems. To account for this artificial time-shift, the phase of the forcing is updated at the transitions. With this convention, the evolution of the system is fully described by I I :ẍ I I (t) + x I I (t) + 1 = A cos(ωt + φ I I ), where the switching times t I and t I I are defined as The velocity and the phase at the transition is the new initial velocity of the subsequent subsystem, i.e.

System behavior
Here we demonstrate the different types of behavior of system (2) by showing time series and x −ẋ portraits. The unforced system (A = 0) exhibits unbounded solution, see Fig. 2. For the forcing amplitude A = 2 and frequency ω = 0.605 the system (2) exhibits unbounded, periodic, quasi-periodic and chaotic solutions depending on the phase φ and on the initial conditions, see Figs. 3, 4, 5, 6. Our goal in this paper is to better understand the behavior of the system through analytical means.

The evolution of the forced system
Here we solve the forced differential Equations (6) and (7) explicitly for ω ∈ R + \ {1} and A > 0. The case of the unforced system (A = 0) and the case of resonant forcing (ω = 1 and A > 0) is relegated to Appendices A and B.
x I (t) = Aω where t I is the switching time defined in Eq. (8). The final velocity and phase for the solution of subsystem I at the transition t = t I will serve as initial velocity and phase for the solution of subsystem I I , (see Eq. (10)), i.e.
v I I =ẋ I (t I ), φ I I = ωt I + φ I . The position and velocity of subsystem I I (Eq. (7)) satisfying initial conditions where t I I is the switching time defined in Eq. (9). Similarly, the final values of velocity and phase of the solution of subsystem I I at the transition t = t I I will provide the next initial conditions for the solution of subsystem I (see Eq. (11)), i.e.
v I =ẋ I I (t I I ), φ I = ωt I I + φ I I . For a detailed analysis of the response of the forced system, we will introduce a Poincaré map [32,33] and study this discrete map instead of the continuous evolution of the system. In the next section, we investigate the properties of the switching times that are essential to obtain the Poincaré map.

Switching time
In order to find the switching time t I at which transition takes place from subsystem I to subsystem I I , the switching criterion x I (t I ) = 1 (see Eq. (8)) is substituted into Eq. (12). This results in where the switching time t I is the first positive root of Eq. (18 ) and is a function of v I , φ I , A, and ω. We define the following two quantities to rewrite the switching time equation (18) in the following compact form The initial velocity v I ≤ 0 of subsystem I is nonpositive (cf. Eq. (6)), thus the domain of β I (defined in Eq. (19)) is Similarly, to find the switching time t I I of the transition from subsystem I I to subsystem I , the switching criterion x I I (t I I ) = −1 (see Eq. (9)) is substituted into equation Eq. (15 ) resulting in where the switching time t I I is the first positive root of Eq. (22). The quantities α I I (φ I I ) and β I I are defined as and the domain of β I I is 5.1 Continuity of the switching time Not all forcing parameters A, ω yield a continuous switching time function. Figure 8a shows the switching time t I (β, φ) for the forcing parameters A = 2, ω = 0.4. The discontinuous jump of the switching time by varying φ for fixed β = 0.5 is shown in Fig. 8b.
What is the cause of the discontinuity? Geometrically, the roots of the transcendental Equations (20) and (22) represent points of intersection of two cosine functions with different amplitude, phase and frequency, see Fig. 9. The first intersection of the two cosine functions is shown by a red square. It can be seen from Fig. 9 that a small change in the initial phase causes a large discontinuous change in the first intersection point, which corresponds to the switching time t I . This type of discontinuous jump in the switching time can also be observed for systems containing bilinear hysteresis elements [14].
Relating back to the dynamics of the system, the discontinuous change in the switching time is related to reaching of the transition points x I (t I ) = 1 or x I I (t I I ) = −1 with zero velocity, i.e.ẋ I (t I ) = 0 oṙ x I I (t I I ) = 0. Reaching the transition point of a piecewise dynamical system with zero velocity is the socalled grazing phenomenon. Grazing can lead to complex dynamical behavior and bifurcations of dynamical systems [34][35][36]. The resulting complexity for our system is demonstrated in Sect. 8. In Fig. 10 the qualitative change in the dynamical behavior of the system (6,7) is illustrated on x −ẋ portraits.
The initial conditions were chosen, as shown in Fig.  10a, separated by the curve associated with the discontinuous jump of the switching time function t I (β I , φ I ).
The trajectories until the first and second transition are shown in Fig. 10b and c. The initial conditions β I = 0.5, φ I = 1.527 and β I = 0.5, φ I = 1.588 of the trajectories are indicated by squares in Fig.   10a. The switching times for the initial conditions β I = 0.5, φ I = 1.527 are t I = 2.4, t I I = 1.1, while for the initial conditions β I = 0.5, φ I = 1.588 the switching times are are t I = 7.7, t I I = 3.1.
Let us now investigate the conditions under which the switching time function is continuous. First we recast the switching time equations (20) and (22) in the following general form where τ is the smallest positive root of (25), and the quantities θ and B are The switching time function is continuous if The proof is given in Appendix C.

Poincaré map
The switched nature of the problem motivates the construction of a Poincaré map to investigate the behavior of the solutions [37]. Poincaré maps are widely used in piecewise-smooth dynamical systems [33,34,[38][39][40][41]. First, consider the mapping (ẋ Aω Using variables introduced in Eq. (19) we recast Eq. as We can now relate initial values of the variables β, φ to their values at the switching by the map I defined as We rewrite Eq. (33) using the quantities defined in (23) to get The relation between the initial values of the variables β, φ to their values at the switching by the map I I defined as where t I I is the smallest positive root of Eq. (22). To specify the range and domain of these maps, we introduce the Poincaré surfaces I = {(β I , φ I )|x(t) = −1} and I I = {(β I I , φ I I )|x(t) = 1}. Clearly, I and I I are maps from I onto I I and from I I onto I , respectively. The Poincaré map (a.k.a. return map) is now defined as the map of the plane I onto itself after a pair of switchings. The map is therefore obtained by composing the two maps I I and I as Substituting Eqs. (32) and (36) Equation (38) defines the Poincaré map to be used in the subsequent analysis. We note that the switching times t I and t I I implicitly depend on the state variables (β I , φ I ). The implicit dependence of the switching time t I on the variables (β I , φ I ) is due to the switching time equation (20), i.e.
where t I is the smallest positive root of Eq. (39). The implicit dependence of t I I on the variables (β I , φ I ) can be expressed by substituting β I I and φ I I from Eqs. (31) and (29) into Eq. (22). This results in where t I I is the smallest positive root of Eq. (40).

Inverse map
The inverse of the Poincaré map (which is used in Sect. 8 to compute backward iterations) can simply be obtained by reversing time and appropriately changing initial conditions. This yields where the switching times t I I and t I are defined by the equations (43)

Fixed points
Fixed points of the Poincaré map (38) correspond to periodic solutions of (2). The fixed points of the Poincaré map are given by From Eq. (44) the fixed points can be expressed as where From Eq. (46) and using that φ ∈ S 1 we get Using Eq. (48) the expression (45) can be written as From Eq. (44) the fixed points can be expressed as Substituting Eqs. (51) and (52) into Eq. (39) we get the following implicit expression for the switching time t I of the period-one solutions Equation (53) can be solved numerically for the possible switching times t I ∈ (0, max(2π, 2π/ω)) for given forcing parameters A and ω. Substituting these switching time values into (51) and (52) we get the potential fixed points of the Poincaré map. Finally, we have to check that these (β, φ) pairs are real fixed points of the Poincaré map (38), since the property that the switching times t I and t I I are the smallest positive roots of Eqs. (39) and (40) must be fulfilled. Figure 11 illustrates the fixed points of the Poincaré map (38) as the function of the forcing frequency ω for forcing amplitude A = 2. For the forcing amplitude A = 2 fixed points exist if ω ∈ [0.375, 1] ∪ (1, 1.721) and k = 1.
Each of the fixed points of the Poincaré map (38) corresponds to a period-one solution of the system (2). Figure 12a illustrates the unique period-one solutions in the x −ẋ plane for A = 2 and ω ∈ {0.4, 0.8, 1.2}.
If the forcing amplitude A = 2 and the frequency ω ∈ (0.575, 0.638) three fixed points coexist (see Fig. 12b). Figures 13 and 14 illustrate the fixed points of the Poincaré map (38) as the function of forcing frequency ω for forcing amplitude A = 8 and A = 15, respectively. For forcing amplitude A = 15 we can also observe 1 : 2 and 1 : 3 subharmonic resonances of the system, thus fixed points for k = 2 or k = 3 coexist for certain ranges of forcing frequency with the k = 1 fixed point, see Fig. 15.
Two of the five different period-one solutions for A = 15 and ω = 2.13 are illustrated in Fig. 16a, where the solid line corresponds to 1 : 1 resonance (k = 1), while the dash-dotted correspond to 1 : 2 resonance   (38) can be written as The partial derivatives ∂t I ∂β , ∂t I I ∂β , ∂t I ∂φ and ∂t I I ∂φ are obtained by differentiating Eqs. (39) and (40). The eigenvalues of the matrix D determine the stability of the fixed points of the Poincaré map (Eq. (38)) or equivalently the stability of the period-one solutions of the system Eq. (2). The eigenvalues of D are given by The trace and the determinant of D are The determinant of D is equal to 1 (for the derivation see Appendix D). Since λ 1 λ 2 = det(D ) = 1, there are three possibilities for the eigenvalues: 1. Both λ 1 and λ 2 are real and distinct. In this case, one has a modulus greater than one (eigenvalue outside the unit circle) and the other smaller than one (eigenvalue inside the unit circle). This fixed point is a saddle. 2. λ 1 and λ 2 are complex conjugate with |λ 1 | = |λ 2 | = 1 (eigenvalues on the unit circle). The fixed point is a center. 3. Either λ 1 = λ 2 = 1 or λ 1 = λ 2 = −1. The fixed point is a non-hyperbolic fixed point and nonlinear analysis is required to determine the behavior of the fixed point. From Eq. (57), we note that the eigenvalues λ 1,2 are real and distinct if |Tr D | > 2, and they are complex conjugate if |Tr D | < 2. Figure 17 illustrates the fixed points and their stability as the function of the forcing frequency ω for forcing amplitude A = 2. The center type fixed point loses its stability at ω = 0.5 and becomes a saddle point. At ω = 0.5 a line of non-isolated fixed points can be observed. Further increasing the forcing frequency ω a supercritical pitchfork bifurcation occurs, which result the coexistence of three fixed points: two centers and a saddle point.
Dynamics in the φ − β plane for A = 2, ω ∈ {0.49, 0.51, 0.6, 0.7} in the vicinity of the fixed points is illustrated in Fig. 18. The invariant curves around the center points can be observed. The attracting and repelling directions of the saddle points are depicted by arrows. Figure 19 illustrates the fixed points and their stability as the function of the forcing frequency ω ∈ [1,4] for forcing amplitude A = 15. Saddle-center bifurcations can be observed at ω = 2 and ω = 4. The dynamics of the system in the φ, β variables on the Poincaré plane for A = 15 and ω ∈ {1.5, 2.13, 2.3, 3.5} are illustrated in Fig. 20.

Higher-period and quasi-periodic solutions
The procedure of finding higher periodic solutions is similar to calculating period-one solutions. A periodn solution of Eq. (2) is a fixed point of the map n . Denoting the fixed point of the map n as (β,φ), from (60) Figure 21a shows the dynamics in the φ−β plane for A = 2, ω = 0.6 in the vicinity of the center type fixed point. Around the center point period-six type invariant curves can be observed. In Fig. 21b we magnified the boxed region of Fig. 21a. In vicinity of the period-six point we observed Period-36, Period-78 points. Similar structures can be observed for the Bogdanov map [42]. Figure 22 illustrates the phase portrait of higher-period and quasi-periodic solutions. Figure 23a shows the dynamics in the φ−β plane for A = 2, ω = 0.7. In Fig. 23b we magnified the boxed region of Fig. 23a to show higher-period solutions and the complex dynamics of the system. In vicinity of the period-three point (empty circle) we observed Period-9 (filled circles), Period-21 (crosses) points. Figure 24 illustrates the x −ẋ portraits of higher-period and quasiperiodic solutions.

Global dynamics
The global dynamics of the Poincaré Map (38) is influenced by the continuity of the switching time function t I (β, φ) + t I I (β, φ). Figure 25a illustrates the sum of switching times t I (β, φ) + t I I (β, φ), and Fig. 25b shows the dynamics in the φ − β plane for A = 2, ω = 0.6. The thick dashed curve is the discontinuity of t I (β, φ) + t I I (β, φ) projected on the φ − β plane. For the forcing amplitude A = 2 three fixed points coexist if ω ∈ (0.575, 0.630). By decreasing the forcing frequency ω the two center points move closer to the discontinuity curve, see Fig. 26. At ω = 0.575 the two center points collide with the discontinuity curve of the switching time and disappear. The presence of centers and saddles in the φ − β plane hints towards possible homoclinic orbits or tangles [26]. We carried out numerical simulations to investigate the complex global dynamical behavior in the presence of saddle-type fixed points. Figure 27 illustrates the dynamics in the φ − β plane for A = 2, ω = 0.605 showing transverse intersection of the stable (blue) and unstable (red) manifolds of the saddle type fixed point. This Figure is the result of forward and backward iterations of 50 × 50 points in a small neighborhood of the saddle. This numerical evidence shows the existence of a Smale horseshoe, which implies the existence of an infinite number of higher periodic and bounded aperiodic/chaotic solutions [32,43].
For the forcing parameters A = 2, ω = 0.605 chaotic solutions were also observed, see Fig. 28.

Conclusions
The dynamics of a spring-mass system with a hysteretic relay operator with simple harmonic forcing were studied in this paper. The energy injected and dissipated by the hysteretic element results in complex dynamical behavior. An explicit Poincaré map with implicit constraints for the switching times has been constructed to facilitate the analysis.
Implicit expression for the existence of periodic solutions have been obtained. On the Poincaré plane center, saddle and non-isolated fixed points were identified. By varying the forcing frequency ω three types of bifurcations occur: saddle-center, supercritical pitchfork bifurcation (two centers and a saddle point) and discontinuity induced bifurcation. At the discontinuity induced bifurcation two fixed points disappear, similarly as in [44].
Higher-period solutions and invariant curves surrounding the center on the Poincaré plane (corresponding to quasi-periodic solutions) have been obtained numerically. The observed homoclinic tangles imply the presence of chaotic solutions.
The limit of the sequence of the switching times is

Appendix B
Resonant case ω = 1 In this section we solve the forced differential equations (6) and (7) explicitly for the case ω = 1. The position and velocity for subsystem I are x I (t) = At 2 cos(t + φ I ) The final velocity and phase for the solution of subsystem I at the transition t = t I will serve as initial velocity and phase for the solution of subsystem I I , i.e.
For the position and velocity of subsystem I I we obtain x I I (t) = At 2 sin(t + φ I I ) + 2 cos(t) x I I (t) = At 2 cos(t + φ I I ) Similarly, the final values of velocity and phase of the solution of subsystem I I at the transition t = t I I will provide the next initial conditions for the solution of subsystem I as v I =ẋ I I (t I I ), In order to find the switching time t I at which transition takes place from subsystem I to subsystem I I , the switching criterion x I (t I ) = 1 is substituted into Eq.
(76) resulting in The switching time t I is the first positive root of Eq.
Substituting the expression of B (defined in Eq. (26)) into condition (88) we get Rearranging condition (92) yields The switching time function τ (β, φ) will be continuous for fixed forcing parameters A, ω if condition (93) is satisfied for ANY φ and β. To establish this criterion we take the maximum of the expression If 0 < ω < 1 and expression (94) is smaller than 0 then switching time τ (β, φ) will be continuous, i.e.
Carrying out a similar calculation for the remaining 3 cases (89)-(91) results in Appendix D