On stability of linear dynamic systems with hysteresis feedback

The stability of linear dynamic systems with hysteresis in feedback is considered. While the absolute stability for memoryless nonlinearities (known as Lure's problem) can be proved by the well-known circle criterion, the multivalued rate-independent hysteresis poses significant challenges for feedback systems, especially for proof of convergence to an equilibrium state correspondingly set. The dissipative behavior of clockwise input-output hysteresis is considered with two boundary cases of energy losses at reversal cycles. For upper boundary cases of maximal (parallelogram shape) hysteresis loop, an equivalent transformation of the closed-loop system is provided. This allows for the application of the circle criterion of absolute stability. Invariant sets as a consequence of hysteresis are discussed. Several numerical examples are demonstrated, including a feedback-controlled double-mass harmonic oscillator with hysteresis and one stable and one unstable poles configuration.


Introduction
The stability of dynamic systems containing hysteresis has been investigated for some time, but still remains a challenging and appealing topic for systems and control theory. Much the same can be said of hysteresis in other fields such as structural mechanics, tribology, magnetism and biology. It should be noted that most of the existing analyses of dynamic systems with hysteresis have their roots in earlier work [3,38]. A combination of passivity [17], Lyapunov [23], and Popov [28] stability theories have been used [27] for analyzing the asymptotic stability and stationary set of systems with hysteresis based on the LMI (linear matrix inequality). Some strong assumptions, including loop transformation with integral nonlinearity and smooth approximation of discontinuities (for example for Preisach [19,29,36] hysteresis operator), have been requested in [27], while the proposed stability theorem required formulating and solving an LMI with certain additional conditions and constraints. Although more recently the absolute stability of dynamic systems with hysteresis operator has been addressed in [26], here the hysteresis was required to be the so-called Duhem operator [24], and a certain matrix inequality, similar to KYP (Kalman-Yakubovich-Popov) lemma [17], had to be solved. Moreover, only examples of second-order systems have been described. The related passivity property, which states that the Keywords and phrases. hysteresis, nonlinear systems, stability, dissipativity feedback connection of two passive systems is passive, has been used for stability analysis of systems with hysteresis in [15]. However, the passivity of hysteresis in feedback has been established for the output rate only and not for the output itself. Several analyses of dynamic behavior of systems incorporating hysteresis have mainly been guided as domain-specific, for example when studying friction where the dissipative properties of hysteresis have been brought into focus, see [1,10,18,21,31].
The objective of this study is to contribute towards stability analysis of systems with linear forward dynamics (of an arbitrary order) and static hysteresis nonlinearity in the feedback channel. While the dynamic hysteresis loops have also been described and analyzed, e.g., in [8], the rate-independent hysteresis nonlinearity will be our main focus, for details see [9]. It is worth noting that this type of hysteresis, where rate-independent maps have a persistent memory of the past after the transients have died out, is particularly relevant for multiple natural and artificial systems, in other words engineering systems. Furthermore it should be stressed that our study focuses on the clockwise hysteresis, for which the (ẏ, ξ) input-output system is passive, cf. [17] for passivity. This is despite the fact that hysteresis with counter-clockwise input-output loops (sometimes referred to as passive hysteresis [27]), are usually considered like Preisach operators, play-type operators and phenomenological or physics-based magnetic hysteresis models. For details see [4]. For analysis of the systems with counter-clockwise input-output dynamics, we also refer to [2]. While the counter-clockwise hysteresis in systems is often in the forward channels of physical signals and hence energy flow (e.g., classical B-H magnetic hysteresis), the clockwise hysteresis usually appears in the feedback interconnection. Examples include nonlinear friction forces, e.g., [1], restoring forces in elasto-viscoplastic structures, e.g., [7], and vibro-impact in the collision mechanisms [16]. Several models for the latter example can also be found in [25].
Before dealing with the main content of the paper, we first need to formalize the hysteresis systems that we will focus on. These have to fulfill two requirements: to be (i) rate-independent and (ii) clockwise in the input-output (I/O) sense. The first requirement means that the hysteresis output trajectories, and consequently the energy conservation or dissipation, do not depend on the input rateẏ. This means a rate-independent (also called static) hysteresis map is invariant to affine transformations of time, i.e., a + bt ∀ a ∈ R, b ∈ R + . The second requirement implies that for any pair of the output trajectory segments which are connected via an input reversal point, the forward segment for whichẏ > 0 always lies above the backward one for whichẏ < 0.
The rest of the paper is organized as follows. In Section 2 we first discuss a dissipative clockwise hysteresis map. This is to highlight its energy dissipation property upon the input reversals. The linear dynamic system with hysteresis feedback is introduced in Section 3 with the associated equilibrium set. The main results in terms of the proposed loop transformation and circle-based absolute stability are described in Section 4. Section 5 contains several numerical examples which illustrate the stabilizing properties of hysteresis in feedback and invariant sets of equilibria. Two boundary cases of equilibria, lying on the coordinate axes of hysteresis map, are demonstrated. Following this, another relevant example of the feedback controlled double-mass harmonic oscillator with hysteresis spring is described in detail. This includes hysteresis with and without discontinuities and a hysteresis-free case of memoryless nonlinear spring. Some general stabilizing properties of hysteresis in the closed-loop are demonstrated and discussed. Finally, the main conclusions are drawn in Section 6.

Dissipative hysteresis map
Considering a static hysteresis map with the input y and output ξ, which represent "energy pair" of a physical system, one can introduce the supply rate w =ẏ(t) · ξ(t). The energy flowing into the system, I/O hysteresis in our case, is obtained as t2 t1 wdt.
Then, for a non-negative real function V (z) : Z → R + of the system state z, called the storage function, the system is said to be dissipative, cf. [37], when satisfying for all (y, ξ, z) and t 2 ≥ t 1 . In other words, for all causal time trajectories, a dissipative system cannot store more energy than is supplied to it, while the difference represents the energy dissipated by the system. It is worth noting that the scalar function V (·) is often also called the energy function, or more specifically the Lyapunov function, since it reflects the instantaneous energy level of the lumped system supplied through y.
Next, we will consider the energy level of hysteresis system ξ(y), which is indirectly analogous to a hysteresis spring which has a displacement-to-force "energy pair". Any changes in energy level and corresponding energy storage, can be considered for two arbitrary points y 2 > y 1 , assumed as exemplary and depicted in Figure 1 without loss of generality. Here y 1 represents an initial state and y 2 represents the next reversal state in the input sequence. Recall that the total energy state of hysteresis system can then be written as where Y (t) is the entire hysteresis input path with Y (0) = y 1 and E(0) is the initial value; for the sake of simplicity we assume E(0) = 0. Note that no kinetic-type energy should be taken into account since the inputoutput hysteresis map is rate-independent [6,36], i.e., E(t) is independent of the input rate. For a monotonically Figure 1. Hysteresis loop ξ(y) with input reversals.
increasing trajectory Y (t) ∈ [y 1 , y 2 ] it can be seen that the change in energy level is given by First, one can show that if ξ(y) has no hysteresis effect, i.e. no loop, then the forward and backward trajectories coincide with each other. Therefore, for the difference in energy levels at both points, one can write Obviously, the total energy change during one y 1 -y 2 -y 1 cycle is ∆E 2−1 + ∆E 1−2 = 0. That means the system behaves as lossless, i.e., without hysteresis-related damping.
Next, we consider an actual hysteresis branching after an input reversal, as shown in Figure 1, so that the corresponding energy change results in The first summand in the above equation constitutes the potential energy which has been stored when moving from y 1 to y 2 , while the second summand corresponds to the potential energy which has been released on the way back. The positive difference ∆E > 0 represents the energy losses due to the hysteresis. Obviously, the amount of loss corresponds to the area of hysteresis loop, gray-shaded in Figure 1. It is worth noting that the loop does not necessarily need to be closed. As the input proceeds after reversal, the growing (gray-shaded) area reflects the increment in energy losses during each cycle. That is, each flip in the sign of the input rate results in energy dissipation, provided the forward and backward paths diverse from each other. This allows us to introduce the following input-output hysteresis map (equally drawn in Figure 1 by the dotted-line), as the upper boundary of energy dissipation after input reversals. Here h > 0 is a parametrization constant and g(·) is a static (generally nonlinear) map satisfying the sector conditions, cf. (15) below. In the following, we will assume a linear relationship g(y) = γ y with the slope constant γ > 0, this is for the sake of simplicity and without loss of generality. Furthermore, it is worth noting that for h → 0 the input-output map (7) approaches the lower boundary case of energy dissipation ∆E → 0, with a corresponding vanishing of hysteresis. It seems appropriate to assume that all clockwise hysteresis maps which lie between the lower and upper boundary cases, as defined above, will behave as dissipative upon the input cycles. Moreover, this is independent of the shape of the hysteresis loops. For the input-output hysteresis map the following storage function, fully analogous to the generalized energy function of hysteresis spring, can be introduced Note that while the first term of (8) describes the potential energy which is stored and therefore can be recuperated, the second term constitutes the energy dissipation over the path Y traversed within the hysteresis loop. It is worth emphasizing that this class of Lyapunov, corresponding to energy functionals, has already been introduced in [38], for analyzing systems with a hysteresis nonlinearity. Taking time derivative of the storage function (8), for the input-output map (7) with g = γ y, one obtainṡ From the above, one can see that the time derivative of the storage (or Lyapunov) function is always negative definite ifẏ = 0. Only for h = 0 does the system behave as lossless, which means that no hysteresis loops arise after input reversals. Otherwise, any input changes would provoke an energy dissipation at a rate proportional to the input rate. Therefore, a clockwise hysteresis, enclosed by any system dynamics, will act as a damping element on the input cycles. However, no exponential damping rate, as otherwise conventional for linear system dynamics, can be associated with hysteresis. It appears natural since no Laplace transform corresponding to poles allocation is available for analysis of hysteresis behavior. Below, we address the stability of linear dynamic systems incorporating hysteresis in feedback. This is in the spirit of the Lure's type systems [22] which have been thoroughly investigated and formalized for nonlinearities without memory.

Feedback system with hysteresis
The considered feedback system largely consists of a linear dynamics system, given in the state-space notation and hysteresis nonlinearity ξ connected in negative feedback. The linear dynamics system is described by the n × n system matrix A and n × 1 input and output coupling vectors B and C, respectively. The system output of interest is y, while x is a n × 1 vector of the system states. The initial hysteresis state is denoted by ξ 0 . The time argument t is omitted for the sake of simplicity, unless in ξ, and that is for explicitly stressing the memory properties of the hysteresis operator. That means the output at t depends from its previous values on the time interval [0, t). Generally speaking, Φ is what we will understand under the dynamic system with hysteresis feedback, in equivalence to the Lure's systems [22] which have memoryless nonlinearities. Note that although u is often assumed to be a control variable, it can equally appear as some coupling quantity for hysteresis nonlinearity in feedback.
Let us now define the equilibrium state, corresponding to the equilibrium set due to the set-valued hysteresis nonlinearity. It was provided in an elegant way in [3]. Introducing the corresponding transfer function which is a dynamic map from u to y, with s as the Laplace variable, one can write a steady-state equation for hysteresis nonlinearity acting in feedback. This one represents a line in the (y, ξ)-coordinates, which cuts the hysteresis loop at the corresponding points {(y 0 , ξ 0 )}. It is obvious that the steady-state of the dynamic system Φ, cf. (10), is given by while x(t) = x 0 constitutes all steady-state solutions of the closed-loop system (10). It implies that the equilibrium set E = {x 0 } is built up from those (x 0 , ξ 0 ) solutions for which (12) and (13) hold, and that they are independent of the initial (ξ(0), x(0)). Recall that due to a multi-valued hysteresis map, the initial ξ(0) should be explicitly taken into account, along with x(0) of the dynamic system.

Loop transformation and absolute stability
Considering the boundary case, cf. with (7), of hysteresis nonlinearity in feedback, one obtains the output of Φ as a Laplace transformed superposition of both input parts in the loop. Here ξ g [·] maps the static part of feedback nonlinearity, not necessarily limited to the linear g-gain, while ξ h [·] = h sign[·] is the sign-nonlinearity with output rate in argument. By means of the block diagram algebra, with corresponding loop transformations, the closed-loop (14) can be brought into the structure as shown in Figure 2. Both closed loops, coupled in the way of superimposed feedback, are denoted by Φ g and Φ h respectively. One can recognize that while both linear transfer functions are subject to the common input u, the outputs of both nonlinearities act as external disturbances, each one for the opposite loop, i.e., u g = ξ g for Φ h and u h = ξ h for Φ g . This allows us to investigate the input-output stability of Φ g and Φ h separately from each other, as shown in the following.
Transformed loop with both nonlinearities The widely accepted circle criterion [5,32,39], which is closely related to and can be derived from the Popov theorem [28] of absolute stability of the Lure's type systems [22], enables analysis of the dynamical systems Φ. The latter includes ξ which should satisfy the [α, β] sector condition, cf. e.g., [17], Not that the [0, ∞) sector, as a particular case, is also included. The circle criterion says that system matrix A has no eigenvalues on the jω axis and v eigenvalues strictly in the right half-plane. The Nyquist locus of G(jω) does not enter the critical disk D(α, β) and encircles it v times counter-clockwise, that for 0 ≤ α ≤ β, cf. [33,34].
Recall that the open disk D(α, β) lies in the complex half-plane Re[s] ≤ 0 passing through −α −1 and −β −1 as shown in Figure 3. Also note that the critical disk shrinks to the conventional critical point as the cf. [33], that ensures the origin of Φ to be asymptotically stable, while the above stability of A is in effect. Note that for the case when inequality (16) contains ≥ 0, meaning the Nyquist locus G(jω) is admitted to touch the vertical axis −β −1 , only a bounded output stability can be guaranteed, cf. [5].
In the following, we will make a restrictive assumption of a stable linear system dynamics, thus requiring v = 0. Therefore, non-entering of the critical disk remains the single sufficient (but not necessary) condition to be fulfilled by the circle criterion. With regard to the transformed closed-loop of Φ, cf. Figure 2, we request that both Φ g and Φ h satisfy the above circle criterion. This is necessary for the entire system Φ is asymptotically stable. For Φ g , the situation is standard and for the given ξ g ∈ [0, β], it is straightforward to evaluated the fulfillment of circle criterion. For Φ h , which matters in the hysteresis loop, one can show that β → ∞ while α = β due to the sign-operator which switches immediately upon zero crossing of the input of ξ h . Here we stress that the sign nonlinearity lies completely in the first and third quadrants, thus complying with the sector condition of nonlinearity in feedback. The latter is of relevance to avoid violating the formulation of Lure's problem, as in (10) cf. [17,22,33], for which the circle criteria have been originally elaborated. The above sector condition yields that for Φ h , the critical disk D(∞, ∞) = 0 shrinks to the single point, so that the Nyquist locus of jωG(jω) is requested to not enter, or correspondingly encircle, the coordinate origin ∀ |ω| > 0.
It is worth emphasizing that the circle criterion specified above for ξ-type hysteresis nonlinearities is sufficient to prove the global asymptotic stability of Φ. However, based on the above, nothing can be said about the possible stabilizing properties of hysteresis in feedback. This means that for the case when G(s) has unstable poles lying in C + . Notwithstanding, a numerical example of stabilizing hysteresis properties will be provided and commented on in Section 5.2.

Boundary cases with simple dynamics
We look at two boundary cases for which the intersection line of hysteresis nonlinearity (7) coincides with one of the (y, ξ) axes. Recall that (7) captures a limiting shape of the clockwise hysteresis that features the maximal energy dissipation upon changes in the input direction. Here and for the rest of subsection, the unity g, h coefficients are assumed for the sake of simplicity. Recall that hysteresis is in negative feedback with the linear dynamics, as in Φ defined in Section 3.
First, assume a double integrator system, with the state vector x = (x 1 , x 2 ) T , so that It is evident that the corresponding transfer function G 1 (0) → ∞ and the Γ 1 -line rests on the y-axis cf. with [3], while it cuts the hysteresis loop as schematically shown in Figure 4. Since A 1 is singular, i.e., not invertible, no Figure 4. Graphical representation of equilibrium sets within hysteresis nonlinearity equilibrium state can be found by x 0 (t) = A −1 1 B 1 ξ 0 for y 0 + G 1 (0)ξ 0 = 0 belonging to Γ 1 . At the same time, constitutes an invariant (or stationary) set, in which the trajectories are independent of the initial state x(0). Each point of S 1 represents a steady-state solution of the hysteresis closed-loop (10) with (7) and (17), while acts as a global attractor in the state-space, as shown in Figure 5 (a). Note that the double integrator dynamics does not fulfill the circle criterion, cf. Section 4, and that for the static (memoryless) part ξ g of feedback nonlinearity. Next, we assume a second-order system for which the transfer function G 2 (0) = 0. From y 0 + G 2 (0)ξ 0 = 0 one obtains the vertical Γ 2 -line which cuts the hysteresis at y 0 = 0, as equally visualized in Figure 4. Here one should stress that the (y, ξ)-coordinates are different for G 1 and G 2 input-output dynamics, while the hysteresis nonlinearity in both cases is the same, given by (7). According to (7), the ξ 0 ∈ [−1, 1] at y 0 = 0, while the sign-dependency is governed by sign ẏ 0 (t − ) , i.e., instantly before reaching the output steady-state y 0 . From the equilibrium x 0 (t) = A −1 2 B 2 ξ 0 , for ξ 0 = ±1, one obtains the lower and upper boundaries x 0 = {−1, 1}, 0 of the invariant set Here again, each point of S 2 represents a steady-state solution of the hysteresis closed-loop (10) with (7) and (19), while α 2 = {−1 ≤ x 1 ≤ 1, x 2 = 0} acts as a global attractor of the state trajectories. A set of x 1 (t) trajectories, converging to α 2 , are shown in Figure 5 (b) for a large number of random x(0) initializations. Referring to the circle criterion, cf. Section 4, one should notice that the Nyquist locus of G 2 (s)s is touching the critical disk for ω = 0. Recall that the critical disk for the sign-nonlinearity collapses to the single critical point, thus approaching zero from the left for β → ∞. In this case, Φ is stable in the sense that all sets of initial conditions lead to a bounded y as t → ∞ [5]. Moreover, the Nyquist locus enters the critical point only in the stationary case (ω = 0), that cancels stability issues for a Φ system which has free differentiators in the loop.

Feedback controlled double-mass harmonic oscillator
Consider a feedback controlled two-mass harmonic oscillator, for which the connecting spring features the nonlinear stiffness with hysteresis. The free body diagram, shown in Figure 6, visualizes the corresponding closed-loop dynamic system. The setup represents a class of common engineering problems where the inertial elements are interconnected by an energy storage which is coupled with (structural) energy dissipation that takes place during harmonic oscillations. Examples for this can be found in the machine elements [12,30], Figure 6. Closed-loop harmonic oscillator with nonlinear spring piezoelectric materials [11,14], micro-frictional and generally tribological systems [1,35], structures with elastoplasticity [7,13] and others. An open-loop system is inherently passive and therefore remains stable even if a harmonic oscillator without auxiliary damping can reveal nontrivial transient trajectories in the presence of the hysteresis. However, in order to explicitly address the stability of the closed-loop dynamics with hysteresis, a negative feedback of x 2 -state is incorporated. This allows for a direct analogy, for instance to an output feedback controller, where K > 0 is the control gain of an appropriate choice. Note that despite mimicking an output feedback control structure, we deliberately omit to insert any control damping. This is intentional since our goal is to demonstrate the stability and possible stabilizing properties of hysteresis in the loop. For the assumed unity masses, the closed-loop system (10) is given by while for the static map ξ we will consider both a memoryless stiffness and stiffness with hysteresis. These are addressed separately in the following subsections. The state vector is x = (x 1 , x 2 , x 3 , x 4 ) T , while non-zero initial conditions x(0) = (1, 0, 0, 0) T will be assumed for the solutions shown below. Note that the linear stiffness of connecting spring, assigned to be g = 100, is already captured in the system matrix A 3 . Further we note that an auxiliary (low-valued) viscous damping term 0.01 × x 4 is added to the second mass, in order to provide a better separation between two pairs of the conjugate-complex poles of A 3 . This is also to numerically stabilize computation of the solutions since we are addressing two boundary cases in the vicinity of the stability margin of the closed-loop system.

Memoryless stiffness
First, we consider a linear memoryless stiffness of the connection spring so that ξ = 0. In this case an oscillating behavior of the system will depend on assignment of the feedback gain. Analyzing the system matrix A 3 , one can show that the closed-loop system becomes marginally stable once K = g. For a stable and unstable system response, we define a boundary-stable and boundary-unstable feedback gain K = 99 and K = 101, respectively. The associated poles-zeros configuration for both cases are shown in Figure 7 (a) and (b). One can see that a relatively small variation of the feedback gain, i.e., by some 2 %, shifts one boundary-stable conjugate complex pole pair into the right half-plane of the s-domain. This results in an unstable configuration of the system poles and, as implication, the overall closed-loop dynamic system becomes unstable, as can be seen from the phase portrait depicted in Figure 8 (b). Opposite, the phase portrait of the same (x 1 , x 2 )-states is shown in Figure 8 (a) for the boundary-stable feedback gain.

Hysteresis stiffness
Next, we give an example of the stabilizing property of hysteresis in the loop, when keeping the boundaryunstable K = 101 feedback gain, as addressed above. The linear dynamics of the system Φ remains the same as in Section 5.2.1, while the single modification is made in the nonlinear part of (7). We are assuming the signdependent hysteresis term with h = 50, further denoted as sign-hysteresis. In this case the feedback nonlinearity  coincides exactly with (7). Another, instead of h sign(ẏ), we then assume a simple single Prandtl-Ishlinskii (PI) stop-type operator [19,20], which represents one rheological element, further denoted by PI-hysteresis. Recall that PI-hysteresis produces an additional slope in the (y, ξ) coordinates, thus avoiding the sign-dependent step-wise transitions. The single rheological element is parameterized by two constants, one giving the output saturation level (analogous to h for the sign-hysteresis), and one determining the slope after the input reversals. For details see [31]. For the closed-loop response, with the same initial conditions as in Section 5.2.1, the input-output map of both nonlinearities are shown in Figure 9 (a), without incorporating the linear stiffness part (already captured in A 3 ). The sign-hysteresis converges to a final state, while the numerical integration at discontinuities produces a (minor) spurious drift, see the red vertical bar in Figure 9 (a). The same can be said of the PI-hysteresis trajectory while multiple reversal slopes still remain inside the surrounding major loop. The corresponding x 2state trajectories, shown in Figure 9 (b), indicate the appearance of a stable limit cycle in both cases, while the second superimposed harmonic, obviously related to additional stiffness of PI-hysteresis, is also visible during initial transients.
Note that in both the above examples, inclusion of the hysteresis provides a stabilizing effect, while the linear feedback dynamics (21) is inherently unstable for K = 101. However, the stability analysis based on the circle criterion given in Section 4, does not deliver conclusive results since one of the transformed closed loops, cf. Figure 2, is unstable. Nevertheless, an additional stabilizing feature of hysteresis points to an important effect in shaping the closed-loop system dynamics. This does not appear to be controversial and can be interpreted as a local stability property, while the absolute stability cannot be guaranteed in this case. Local stability is supposedly related to the amount of energy dissipated at the hysteresis cycles and is therefore dependent on the hysteresis shape and initial conditions in the otherwise unstable closed-loop dynamics. These additional aspects of hysteresis stabilization are equally interesting, but are beyond the scope of this study. They are certainly worth investigating in the future.

Conclusions
In this paper, we have addressed the stability of linear dynamics with clockwise hysteresis in feedback. Starting by analyzing a dissipative hysteresis map, we have demonstrated the property of energy dissipation once the input changes and for each input reversal. In this context, we have derived an upper boundary case of a clockwise hysteresis, constructed with the help of a simple sign operator applied to the input rate. In contrast, the lower boundary case (where the hysteresis vanishes and the hysteresis loops collapse to the static nonlinearity without memory), is shown to be non-dissipative and therefore lossless from an energy viewpoint. Other hysteresis shapes lying in-between the boundary cases, i.e., having non-zero loop area, are shown to be equally dissipative, cf. Figure 1. A general class of linear dynamic systems with hysteresis in feedback has been considered in the spirit of the seminal work [3]. In addition to existing approaches, a simple loop transformation which allows for stability analysis by means of the circle criterion has been introduced. It has been shown that the standard circle criterion should be satisfied for both parallel loops of the transformed dynamics, upon which the absolute stability can be concluded. The numerical examples have been demonstrated for two boundary cases of the second-order system dynamics. In addition, a more demanding fourth-order dynamics of a feedback controlled double-mass harmonic oscillator has been treated and discussed for one stable and one unstable poles configuration. In this regard, some general stabilizing hysteresis properties have also been discussed.