MODELING , CHAOTIC BEHAVIOR , AND CONTROL OF DISSIPATION PROPERTIES OF HYSTERETIC SYSTEMS

In the present work, hysteresis is simulated by means of internal variables. Analytical models of different types of hysteresis loops enable reproduction of major and minor loops and provide good agreement with experimental data. Using an effective algorithm based on the analysis of wandering trajectories, an evolution of chaotic behavior regions of oscillators with hysteresis is presented in various parametric planes. A substantial influence of a hysteretic dissipation value on the form and location of these regions and also restraining and generating effects of the hysteretic dissipation on chaos occurrence are shown. It was demonstrated that for fixed parameters, which govern the shape of a hysteresis loop, the extent of pinch of this loop could be controlled by means of the amplitude and frequency of external periodic excitation.


Introduction
Hysteresis is caused by very different processes in nature and is characterized above all by the presence of an output delay with respect to input in an input-output correlation, energy dissipation, and memory in a system.The problem of hysteretic systems investigation occurs in many fields of science for mechanical, engineering, physical, and biological systems and even for sociological and economic systems.
The present paper is composed of five sections.Section 2 is devoted to the simulation of hysteretic systems.There is a lot of different phenomenological approaches to hysteresis modeling.General mathematical models of hysteresis are presented and discussed in the monographs [21,28].Recent results of the analysis of hysteretic systems including modeling, experiments, dynamic response, and applications can be found in [27].A large number of publications [7, 9-15, 17-19, 22-26] are devoted to hysteresis simulation, since it is known that hysteretic systems are complicated to investigate and various difficulties occur when applying the existing models.The question of multipurpose and generally valid models describing the wide spectrum of hysteretic phenomena is still open.In the present work, hysteresis is simulated by means of additional state variables (internal variables).In particular, the behavior of magnetorheological/electrorheological (MR/ER) fluids (which are known as smart materials and commercially available now) in a damper/absorber [24,25] is simulated.Another example is modeling of stress-strain hysteresis in a steel rope including minor-loop reproduction.The developed models are effective and contain principally less parameters than, for example, Bouc-Wen or Spencer models.
Section 3 is dedicated to the analysis of chaotic behavior occurring in the hysteretic systems in various parametric spaces.The models describing systems with hysteresis are discontinuous and contain high nonlinearities with memory-dependent properties.The output is delayed with respect to the input and for every input there may be more than one equilibrium state.Investigation of these systems within the framework of an approximate analytical approaches as, for instance, slowly varying parameters or harmonic balance methods, results in a conclusion that irrespective of the values of control parameters occurring in the conditions of external periodic excitation, the hysteretic system has a stable symmetric asymptotic response.However, the recent publications and works [8,16,20], based on numerical and combined numerical-analytical techniques, present frequencyresponse curves and bifurcation diagrams which indicate the presence of solutions and bifurcations mostly unexpected for hysteretic oscillators.At the same time, knowledge of the control parameter spaces is still insufficient.In this connection, the prediction of conditions for stable/unstable behavior of such systems is very topical.In the present work, a chaotic behavior occurring in the dynamic hysteretic system governed by a coupled differential set is investigated in various parametric planes using a methodology described in [1][2][3][4][5][6].This methodology had been successfully applied already in particular to predict stick-slip chaos in a weakly forced oscillator with friction [3], in 2-DOF discontinuous systems with friction [2,4], and chaos in other smooth and nonsmooth systems [4].
Section 4 is devoted to the control of dissipation properties of oscillators with hysteresis.

Hysteresis simulation by means of additional state variables (internal variables)
In modeling of systems showing hysteresis based on the parametric approaches, it is intended to use various linear and nonlinear elements simulating memory in the system as well as the energy dissipated in each cycle.In a physical sense, energy losses can be described by including into the model, for example, dashpots, friction elements, springs, and other mechanical elements.Thus, the classical Masing model of hysteresis merges the element of Hooke (elastic body model) and a number of St. Venant's elements (plastic body models) in parallel.The Biot model of hysteresis includes some number of Newton's elements (Newton flow models) instead of St. Venant's elements.Not infrequently successive joints of mechanical (physical) elements are also used as, for example, in the Spencer model which is an extension of the Bouc-Wen model.The parameters of each model can be defined during a parameter identification process with experimental data.The model validity is confirmed also by the experiment.J. Awrejcewicz and L. Dzyubak 3 Let us consider a combination of N friction elements with the maximum friction forces F 1 ,F 2 ,...,F N and springs k 0 ,k 1 ,k 2 ,...,k N in parallel, as it is shown in Figure 2.1.Here, x is the input (input signal) and z is the output (response) of the hysteretic system.The delay in the arrival of the output with respect to the input can be described with the aid of internal variables (forces) y 1 , y 2 ,..., y N .
If the absolute value of the internal variable |y i | (i = 1,2,...,N) is less than maximum friction force F i , or if both this value is equal to F i and the velocity of the input ẋ and the force y i have different signs (including the case when ẋ = 0), then the evolution of force y i in time is governed by the following equation: In all other cases, this evolution is as follows: Equations (2.1) and (2.2) are generalized into the following equation: Using the approximation equation ( 2.3) can be rewritten in the form and then the response of the hysteretic system is A generalization of (2.5) and (2.6) in the case of friction elements and springs with properties of a more complex character (in particular of nonlinear character) yields where k 0 (x),A i (x) ≥0, F i (x) >0, the input signal belongs to admissible codomain x ∈ [x min ,x max ], α i ,β i ∈ R; or it can be written even in a more general form z = p x, y 1 , y 2 ,..., y N , ẏi = q x, ẋ, y i , i = 1,2,...,N. (2.8) Here p and q are the nonlinear (in a general case) functions of their arguments.They are chosen depending on the properties of a hysteretic system and a loop form.Note that, according to (2.7), the structure of differential equations describing the evolution in time of the internal variables is very similar to the structure of the differential equations of the Bouc-Wen hysteretic oscillator model.
The parameters of functions p(x, y 1 , y 2 ,..., y N ) and q(x, ẋ, y i ) (i = 1,2,...,N) are determined via a procedure minimizing the criterion function Φ c 1 ,...,c j ,α 1 ,...,α k ,...,β 1 ,...,β l = i p x c 1 ,...,c j ,t i , y 1 α 1 ,...,α k ,t i ,..., y N β 1 ,...,β l ,t i − z i 2 , ( which characterizes an error between the experimental and calculated curves.Here z i are the responses of a hysteretic system, which are known from the experiment and the values p(x(c 1 ,...,c j ,t i ), y 1 (α 1 ,...,α k ,t i ),..., y N (β 1 ,...,β l ,t i )) are obtained from the integration of the system which is described by model (2.8).To minimize the criterion function (2.9), the method of gradient descent is used.The step-by-step descent to the minimum of the criterion function is performed in the opposite direction to the criterion function gradient the value of the criterion function Φ increases or remains unchanged, it is necessary to decrease the values of the steps h c1 ,...,h cj ,h α1 ,...,h αk ,...,h β1 ,...,h βl .The model presents the fast enough numerical convergence.Applications to different types of hysteresis loops confirmed that models with internal variables were appropriate to simulate hysteresis.
Consider two examples of such applications.As it was mentioned, hysteresis is widely found in nature.Among others, hysteretic behavior is peculiar to "smart" materials that more and more attract the attention of commercial structures.Electrorheological/magnetorheological (ER/MR) fluids belong to such materials.These fluids are used for construction of dampers, which act as interfaces between electronic control systems and mechanical systems.The dynamic characteristics of ER/MR dampers are strongly nonlinear that is reflected in numerous results of the experimental studies of MRF damper behavior [24].Physically accurate and simple models of dampers are highly required by engineering community and they are necessary for producing effective samples.Due to the nonlinearities of hysteresis and jumps in damper behavior, difficulties arise while developing these models [25].
During analysis of seven known parametric models in [25], an attempt was made to find a compromise between simplicity, "economy" of the model, and its physical accuracy.However, it was concluded that the simplest involution model with two parameters could not be applied to simulate the damper behavior.An extension of the Bouc-Wen model proposed by Spencer enables the most accurate prediction of an actual MR damper behavior, but this model contains nine parameters.
We suggest here the developed model with a single internal variable  Further, we consider the stress-strain hysteresis in a steel rope.It happens that when a force acts upon a steel rope, some transient processes are reflected in minor loops reproduction.In this case, there are two (or more) ways to simulate these processes in hysteresis behavior of the steel rope.The first one is to increase the number of internal variables.But it leads to a complication of the model and increase in the number of parameters.
The model with five parameters, which contains a single internal variable and quite well describes the major stress-strain hysteresis loop, is

Evolution of chaotic behavior regions in control parameter planes of the Masing and Bouc-Wen hysteretic oscillators
Consider classical hysteretic models such as Masing oscillator and Bouc-Wen oscillator.
In both cases, an external periodic excitation with an amplitude F and frequency Ω acts on the mass m which oscillates along an inertial base.These oscillators possess hysteretic properties and it is supposed that there is a linear viscous damper with the coefficient 2μ.
The following set of differential equations governs a motion of the Masing oscillator: (3.1) + νz is the total restoring force with nonlinear elastic part (1 − ν)g(x) and with hysteretic part νz.The case ν = 1 corresponds to the maximum hysteretic dissipation and ν = 0 corresponds to the elastic behavior of the oscillator.The parameter δ characterizes a ratio between the post-and preyielding stiffness.The parameter n governs the smoothness of the transitions from the elastic to the plastic range.The couples ±(x i ,z i ) represent the velocity reversal points at ẋ = 0.According to Masing's rule which is extended to the case of steady-state motion of the hysteretic oscillator, the loading/unloading branches of a hysteresis loop are geometrically similar.So, if f (x,z) = 0 is the equation of a virgin loading curve, then the equations f ((x ± x i )/2,(z ± z i )/2) = 0 describe loading/unloading branches of the hysteresis loop.In the case of non-steady-state motion of Masing's oscillator, it is supposed that the equation of any hysteretic response curve can be obtained by applying the original Masing rule to the virgin loading curve using the latest point of velocity reversal.
A motion of the Bouc-Wen oscillator is governed by the following set of differential equations: where R = δx + (1 − δ)z is the total restoring force; the parameters (k z ,β,n) ∈ R + and γ ∈ R govern the shape of the hysteresis loop.The parameters δ and n have the same sense as in the case of the Masing model.
Describe briefly the approach based on the analysis of wandering trajectories in view of the state vector of systems (3.1) and (3.2) for x ∈ R 3 .In the chaotic behavior of nonlinear deterministic systems, wandering of the trajectories of motion around the various equilibrium states is assumed.They are characterized by unpredictability and sensitive dependence on the initial conditions.By analyzing the trajectories of motion of these systems, it is possible to find the chaotic vibration regions in the control parameter space.
The continuous dependence property on the initial conditions x (0) = x(t 0 ) of the solutions of set (3.1) or (3.2) will be used.For every initial condition x (0) , x (0) ∈ R 3 , for every number T >0, no matter how large, and for every preassigned arbitrarily small ε > 0, it is possible to indicate a positive number δ > 0 such that if the distance ρ between x (0) and x (0) , ρ(x (0) , x (0) ), is less than δ, and |t| ≤ T, the following inequality is satisfied: In other words, if the initial points are chosen close enough, then during the preassigned arbitrary large time interval −T ≤ t ≤ T the distance between simultaneous positions of moving points will be less than a given positive number ε.
For the sake of tracing chaotic and regular dynamics, it is assumed that, with an increase of time, all trajectories remain in the closed bounded domain of a phase space, that is, To analyze trajectories of sets (3.1) and (3.2), the characteristic vibration amplitudes A i of components of the motion are introduced A i = (1/2)|max t1≤t≤T x i (t)−min t1≤t≤T x i (t)|.Index number i runs over three values corresponding to three generalized coordinates x, y, z. [t 1 ,T] ⊂ [t 0 ,T] and [t 0 ,T] is the time interval, in which the trajectory is considered.The interval [t 0 ,t 1 ] is the time interval, in which all transient processes are damped.
For the sake of our investigation, it seems the most convenient to use the embedding theorem and to consider a 3-dimensional parallelepiped instead of a hypersphere with the center in point x.Two neighboring initial points x (0) = x(t 0 ) and x (0) = x(t 0 ) (x = (x, y,z) T or x = (x 1 ,x 2 ,x 3 ) T ) are chosen in the 3-dimensional parallelepiped P δx,δy,δz (x (0) ) such that |x (0) i − x (0) i | < δ i , where δ i >0 is small in comparison with A i .In the case of a regular motion, it is expected that ε i >0 used in inequality |x i (t) − x i (t)| < ε i is also small in comparison with A i .The wandering orbits attempt to fill up some bounded domain of the phase space.In instant t 0 , the neighboring trajectories diverge exponentially afterwards.Hence, for some instant t 1 , the absolute values of differences |x i (t) − x i (t)| can take any value in closed interval [0,2A i ].An auxiliary parameter α is introduced, 0 < α < 1. αA i is referred to as a divergence measure of observable trajectories in the directions of generalized coordinates and with the aid of parameter α one has been chosen, which is inadmissible for the case of the motion "regularity."The domains, where a chaotic behavior of considered systems is possible, can be found using the following condition: If this inequality is satisfied in some nodal point of the sampled control parameter space, then such motion is relative to chaotic one (including transient and alternating chaos).The manifold of all such nodal points of the investigated control parameter space sets up domains of chaotic behavior of the considered systems.
Motion stability depends on all parameters of the considered hysteretic models including initial conditions.We succeeded in tracing irregular responses of the Masing and Bouc-Wen hysteretic oscillators in the damping coefficient-amplitude (μ,F) and frequency-amplitude (Ω,F) of external periodic excitation planes after a coordinate sampling.
The Masing oscillator (3.1) is nonlinear both in the case of a pure elastic behavior without hysteretic dissipation (ν = 0) and in the case of motion with hysteretic dissipation (ν > 0).At ν = 0 chaos has been found too.Figures 3.1 and 3.2 display the evolution of chaotic behavior domains with an increasing hysteretic dissipation value in the mentioned planes.The (Ω,F) and (μ,F) planes had been uniformly sampled by 100 × 100 nodal points in the rectangles (0.01 ≤ Ω ≤ 0.61; 0.01 ≤ F ≤ 1.51) and (0.01 ≤ μ ≤ 0.26; 0.01 ≤ F ≤ 1.61), respectively.The time period for the simulation T is 300π/Ω nondimensional time unit.During computations, half of the time period T corresponds to the time interval [t 0 ,t 1 ], where transient processes are damped.The integration step size is π/100Ω.Initial conditions of the closed trajectories are distinguished by 0.5 percent with ratio to characteristic vibration amplitudes, for example, the starting points of these trajectories are in the 3-dimensional parallelepiped (|x The parameter α is chosen to be equal to 1/3.All domains are multiple connected.There are also a number of scattered points here.Such structure is characteristic of the domains, where chaotic vibrations are possible. One can observe the effect of restraining of the chaotic regions with the increasing hysteretic dissipation value in the (Ω,F) plane (Figures 3.  of the restraining decreases when ν increases.So, in the case of a maximum hysteretic dissipation value ν = 1, the chaotic regions in the (Ω,F) plane are distinguished from the regions (c), Figure 3.1, non-principally.
In the (μ,F) plane (Figures 3.2(a), 3.2(b), 3.2(c)), the form and location of the chaotic domains are changed depending on the hysteretic dissipation of the nonlinear terms in set (3.1).
Figure 3.3 characterizes the obtained domains and demonstrates a different character of motion of the Masing oscillator as chaos and hysteresis loss (a), and periodic response (b).
Other situation occurs for the Bouc-Wen oscillator (3.2) which naturally is linear (when δ = 1).So, addition of the hysteretic dissipation leads to chaotic responses that occur in this system.Figures 3.4 and 3.5 present the evolution of chaotic behavior regions with an increasing hysteretic dissipation value in the (Ω,F) and (μ,F) planes, respectively.One can observe a change in the form and location of the chaotic regions.Note that chaotic responses of the Bouc-Wen oscillator are not observed right up till δ = 0.2 when the influence of the nonlinear terms becomes critical.It demonstrates a generating effect of the hysteretic dissipation on chaos occurring in the hysteretic system, which appears after some critical value δ cr .After δ cr , both the form and location of the chaotic behavior regions are changed with the increasing hysteretic dissipation.In the case of a maximum hysteretic dissipation value (when δ = 0), the chaotic behavior regions are practically the same as in case (c) in Figure 3.4.A "friable" form of the chaotic regions (Figures 3.During simulation, the (Ω,F) and (μ,F) planes had been uniformly sampled by 100 × 100 nodal points in the rectangles (0.01 ≤ Ω ≤ 0.36; 0.01 ≤ F ≤ 2.05) and (0.001 ≤ μ ≤ 0.04; 0.01 ≤ F ≤ 1.71), respectively.The time period for the simulation T is 300π/Ω.During computations, half of the time period corresponds to the time interval [t 0 ,t 1 ], where transitional processes are damped.The integration step size equal to π/40Ω is chosen.As in the case of the Masing oscillator, the initial conditions of the closed trajectories are distinguished by 0.5 percent in comparison to characteristic vibration amplitudes.The parameter α is set equal to 1/3.

Control of pinched hysteresis phenomenon by means of an amplitude and frequency of an external periodic excitation
Pinched hysteresis possesses reduce dissipation properties and are a consequence of residual phenomena diminution or "delay" reduction of the output relative to the input in a hysteretic system.In this connection, a control of pinching phenomena is of great importance for various practical applications.It was shown that for the nonzero fixed value of hysteretic dissipation (which, in changing from minimum to maximum, can produce behavior of the system with/without delay) and for the fixed parameters which govern the shape of a hysteresis loop, it is possible to choose an "appropriate" hysteretic process with desirable residual phenomena with the aid of an amplitude and frequency of an external periodic excitation: (i) Masing oscillator ẋ = y, (ii) Bouc-Wen oscillator ẋ = y, i = 1,2,...,n; j = 1,2,...,m.The relation ε r /ε was chosen as a value that characterizes a pinch in the hysteresis loop.As shown in Figure 4.1, ε is the distance between the projection of the velocity reversal points ±(x i ,z i ) to the input-axis; ε r is the distance characterizing the residual phenomena.
It is shown that the thresholds of the regions with a "pinched hysteresis" in the planes frequency versus the amplitude of external periodic excitation in both cases of Masing and Bouc-Wen hysteretic oscillators (Figures 4.

Conclusions
In the present work, hysteresis is simulated by means of additional state variables (internal variables).In particular, the behavior of magnetorheological/electrorheological (MR/ER) fluids in a damper/absorber as well as stress-strain hysteresis with transient processes in a steel rope are simulated.The developed models are effective, enable production of minor loops, present fast numerical convergence, provide good agreement with experimental data, and contain principally less parameters than, for example, Bouc-Wen or Spencer models.
Highly nonlinear Masing and Bouc-Wen hysteretic models with discontinuous righthand sides are investigated using an effective approach based on the analysis of wandering trajectories.This algorithm of quantifying regular and chaotic dynamics is simpler and faster from the computational point of view comparing with standard procedures and allows us to trace accurately enough the regular/irregular responses of the hysteretic systems.The evolution of chaotic behavior regions of the oscillators with hysteresis is presented in various control parameter spaces: in the damping coefficient-amplitude and the frequency-amplitude of the external periodic excitation planes.A substantial influence of the hysteretic dissipation value on the possibility of chaotic behavior occurring in the systems with hysteresis is shown.The restraining and generating effects of the hysteretic dissipation on a chaotic behavior are demonstrated.
It was shown that for fixed parameters, which govern the shape of a hysteresis loop, the extent of pinch of this loop could be controlled by means of an amplitude and frequency of the external periodic excitation.It was found that a linear viscous damper did not influence the pinching.The regions of hysteresis with various dissipation properties for the Masing and Bouc-Wen models are presented.

5
(a) and 3.5(c)) in the (μ,F) planes is conditioned by the fixation of the frequency Ω = 0.24 and changing location of the corresponding domains in the (Ω,F) planes with an increase of the hysteretic dissipation.The straight line Ω = 0.25 only slightly contacts the chaotic behavior regions in cases (a) and (c), in Figure3.4.

Figure 3 .
6 characterizes the obtained regions of irregular motion and depicts various responses of the Bouc-Wen oscillator as chaos and hysteresis loss (a), and periodic response (b).Figures 3.7(a), 3.7(a) and 3.8(a), 3.8(a) with periodic pinched hysteresis also agree well with the obtained regions of regular/irregular behavior of hysteretic oscillators (see Figures 3.1(b), 3.2(b), 3.4(a), and 3.5(a)).

Table 2 .
1. Final values of the parameters used in model (2.13) for identification of the experimental data.

Table 2 .
2. Final values of the parameters used in model (2.15) for the experimental data identification.