A Simple Model of the Energy Harvester within a Linear and Hysteresis Approach

In this article, a model of an energy harvester, the mechanical part of which is an inverted pendulum, is proposed. We investigated the stability of a linearized system. It was proven that the stabilizing control of the pendulum, based on the feedback principle, enables the stabilization of the system. We have identified the zones of stability and the amplitude–frequency characteristics. In the second part of this article, a generalization of the dynamic system for the case of the hysteresis friction in the mechanical joint is considered. The role of nonlinear effects within the design Preisach model and the phenomenological Bouc–Wen model is shown.


Introduction
To this date, improvements in the design and performance of energy harvesting systems are the subject of intense research. Modern advances in materials science, electronics, and control theory have reduced the size, reliability, and cost of such devices; nevertheless, achieving a performance comparable to traditional electric batteries continues to be an urgent scientific and technological challenge.
The process of energy saving, as it is known, consists of a redistribution of kinetic energy of an oscillating massive body into the form of electric energy. The non-linearity of the characteristics of the mechanical and/or electrical subsystems is of crucial significance. The importance of non-linear links in energy converters is emphasized in the works [1][2][3][4]. After all, linear systems have a sufficient coefficient of efficiency only under the condition of approximate equality of proper frequencies of mechanical and electrical subsystems, i.e., under the condition of closeness to the resonance mode of functioning [5,6]. Systematic changes in mechanical parameters are inevitably observed in industrial conditions, arising due to mechanical wear and ageing of materials and leading to degeneration of transfer characteristics. Let us also note the important influence of external noises, which can be correctly assessed only by probabilistic and statistical methods, for which we have to introduce probabilistic characteristics for individual system components [7].
Various mechanisms for maintaining resonance modes by means of vibrating oscillating systems with tunable characteristics were considered in the works [8][9][10][11][12][13][14]. However, the systems considered in these works require precision parameter tuning, which is not always feasible in real technical systems. In particular, much attention has been paid to the values of electric voltage [15], temperature conditions [16], and noise exposure [17,18].
In particular, a number of works are devoted to electric circuits that include piezoelectric materials, i.e., substances in which electric polarization arises under elastic deformation. They are crystalline substances without a center of symmetry and are characterized by a complex, non-linear relationship between the applied mechanical stress and the generated electric field. Among piezoelectrics there is the class of ferroelectric materials, which have a non-zero polarization in a certain temperature range, changing due to external excitation [19,20]. In fact, ferroelectrics are the material used to build microelectromechanical systems (MEMS), sound emitters, accelerometers, and precision microdisplacement sensors [16,[21][22][23][24][25][26][27].
The control problem was considered in the articles [28][29][30]. The characteristics of chaotic dynamic modes arising from certain parameters of the system were investigated in [31][32][33][34].
Relatively small mechanical voltages lead to a linear response of the electrical characteristics of the ferroelectric material. On the other hand, in the practical issue of energy harvester design, as a rule, the values of external excitation have to exceed the linear response threshold. Thus, there is a hysteresis dependence of voltage, charge, and other electrical characteristics on the dynamic parameters of the mechanical subsystem [35,36].
A hysteresis loop occurs during consideration of the dependence of the polarization of a ferroelectric material sample on the applied electric field, whereas in an alternating field, the loop parameters depend considerably on the frequency of field change. Furthermore, a hysteresis relationship appears in the volt-farad characteristics of some ferroelectric films, i.e., the capacity depends on the voltage applied to the sample [37,38]. It is quite interesting that the samples of ferroelectric films sprayed on silicon plates of pand ntype conductivity differ in the direction of loop bypass (clockwise and counter-clockwise, respectively). Hysteresis volt-farad characteristics are also demonstrated by metal-oxidesemiconductors (MOS) [39] and crystalline copolymer thin films) [40].
Among the structural models of hysteresis non-linearities, the Preisach operator [50][51][52] plays an important role. This model was originally formulated to describe the properties of ferromagnetic materials [53]; later its applicability to a wide range of phenomena was proved from various scientific and practical problems [54,55]. The Preisach model is based on a non-linear operator, which is a continuum system of relays connected in parallel. The output of the Preisach operator Γ[u(t)] (where u(t) is a continuous function of time) is a function x(t), whose value at each time moment is determined, as for all hysteresis converters, not only by the input u(t), but also by the background.
The Bouc-Wen model is one of the most commonly used phenomenological hysteresis models [48,56]. This model is based on the consideration of a one-dimensional mechanical system in which the return force, Φ BW (x, t), applied to a material point is represented as the sum of two components: where x(t) is the displacement of the body along the coordinate axis, αkx(t) is the elastic return force, k is the stiffness factor, (1 − α)Dkz(t) is the hysteresis term, α is the relative stiffness, i.e., the ratio of the strain above the elastic threshold to the strain below the elasticity, and D is the plastic displacement strain. The auxiliary dimensionless function, z(t), included in the right side of Equation (1), satisfies a differential equation of first order: which includes the time derivative,ẋ(t), of the displacement on the right hand side. The mapping of x(t) → Φ BW (x, t) thus constructed is called the hysteresis model of Bouc-Wen [48].
Note that the integer parameter n determines the smoothness of the transition from elastic to plastic deformation. The parameters A, β, and γ are also dimensionless; they determine the shape of the hysteresis loop.

Energy Converter Model in Linear Approximation
As stated above, the most important parts of the energy converter are the mechanical oscillating system and the associated electromagnetic system, perceiving the oscillation energy. In this paper, it is proposed to choose the mechanical part of the energy converter in the form of an inverse pendulum.
The model of the investigated system consists of an inverted mathematical pendulum fixed on a light horizontal platform and mechanically connected to a ferroelectric condenser, which is included in a closed electrical circuit. The platform P can move horizontally. We denote the angle of deviation of the pendulum with respect to the vertical by ϕ(t), the coordinate of the platform by u(t) (see Figure 1), the length of the pendulum as l, and its mass as m. The connection between the mechanical and electrical subsystems contains a joint obeying the law of viscous friction with a dissipation factor c F fr = cv(t) = c d dt (l ϕ/2 + u)) . Piezoelectric material forming a condenser with capacitance C is included in an electric circuit with an external ohmic load of R. We denote the voltage at the load by V(t).
The dynamic system under consideration is described by a system of ordinary differential equations: The first of the equations in (3) is the equation of motion of weight m under the effect of mechanical and inertial forces, and the second is the balance of currents in an electric circuit. The dot above the symbol indicates the time derivative of t.
Assuming small deviations of the pendulum from the equilibrium position (sin ϕ ∼ ϕ), we obtain (where we have further introduced the notations γ 0 = c/(2m), ω 2 0 = g/l): In this system, A and B are the coupling parameters of the mechanical and electrical subsystems, respectively, which are derived on the basis of the following considerations. The equation with respect to the unknown function ϕ(t) is the equation of motion of the exited oscillator (additionally to mechanical forces) by external force of nonmechanical nature (F ext (t) = A V(t)) and the inertial force F iner (t) = −mü(t).
The elastic strain law (Hooke's law) relates tension σ and the relative strain (for moderate absolute values of elastic strain): where E is the Young's modulus, depending only on the piezoelectric material and its physical state (in the non-linear case σ = Eε + E 2 ε 2 + O(ε 2 )). If there are no tangential stresses, the polarization of the piezoelectric specimen under tension or compression is determined by the expression: where τ x and τ y are the mechanical tensions acting parallel to the axes Ox and Oy, respectively, and d 11 is the constant called the piezoelectric modulus. Express the value of the charge Q = CV(t) generated on the surfaces of the sample (taking S as the area of one surface and C as the electrical capacitance of the condenser formed by the piezoelectric plate): There is a coupling between the elastic force acting from the piezoelectric on the weight m and the displacement δ(t) = l ϕ/2 + u in the form of linear relation: where A is a constant which depends only on the material and its dielectric thermodynamic properties. In other words, the voltage depends linearly on the displacement by the formula where L is the length of the piezoelectric plate (as Hooke's law can be represented as F/S = E × (δ/L)).
To take advantage of the polarization charges occurring on the opposite surfaces of the quartz plate during its deformation, the surfaces are provided with metallic covers. Charges equal and opposite in sign to polarizing charges are induced on such covers, and electric current is generated in the external wires connecting the covers.
In an electric circuit consisting of an external current source (due to polarization of the dielectric plates) and connected condenser C and resistor R in parallel (which acts as an active load), Ohm's law applies: where Note that the literature on energy storage models describes electrical subsystems built on the same principles, but satisfying other differential equations: For example, in ref. [4]: and in refs. [57,58]: where d(x) is the coupling coefficient between the mechanical and electrical subsystems. Among the many options for mechanical pendulum-based oscillation systems proposed for use in energy harvesters we note the inverted pendulum beam [59], the inverted pendulum with amplitude limiters [60], and the inverted cantilever beam with a tip mass [61].
In work [59] it was proposed to construct an energy storage device in the form of an inverted piezoelectric beam and a pendulum. A layer of piezoelectric material was attached to the base of the beam, which is subject to external disturbance in the horizontal direction. Due to the instability of the inverted beam, even weak excitation leads to oscillations with significant amplitude.Additionally, the existence of two oscillation modes, between which there is a mechanical connection, allows the transmission of high-frequency excitation energy to a low-frequency mode, which generates electrical power.
The model used in ref. [60] was based on a vertical solid flexible beam with a massive tip. Instead of the standard clamp, a joint with an additional gap has been introduced.Due to the presence of contacts with limiters and movement around the hinge point, the pendulum was a system with two modes of oscillation: rotation and bending. Calculations showed that the potential energy of the system was characterized by two potential wells, and the responses are subharmonic solutions with multiples of odd frequencies. The authors demonstrated that this mechanical system could be used as a resonator in a broadband energy collector.
A similar device was discussed earlier in ref. [61], where the linear case was analyzed in detail. It was shown that if the beam has a strong bend, then the system shows nonlinear characteristics, including a chaotic regime.
The inverted pendulum is also a key element in the storage of biomechanical energy in a backpack device [62]. It generates energy due to medial-lateral oscillations of carried weight. In general, the concept of an inverted pendulum arises naturally when considering human walking modelling to harvest electrical energy [63].

Converting the System to a Dimensionless Form
Let us switch from variables t and V to dimensionless variables by the following rule: where the characteristic values of time, t c , and voltage, V c , will be determined below. (The deflection angle of the pendulum, ϕ(t), is a dimensionless value.) The system (3) takes the following form: Hereinafter, the point above the symbol denotes the derivative of the dimensionless time τ (two points denotes the second derivative of τ). Let us perform algebraic transformations: Let the function defining the motion of the platform P be u(t) = u(t c τ), and the values of dimensional coefficients be t c = RC and V c = ml A(RC) 2 . Then, we get the basic system in a dimensionless form: To simplify and clarify further calculations, we introduce notations for dimensionless values of damping coefficient, characteristic frequency of the pendulum, and external influence as γ = γ 0 · (RC), ω = ω 0 · (RC), and w = u(τ)/l, respectively. It is also convenient to denote the coupling coefficient of the electrical and mechanical subsystems or, in a more convenient matrix form: where ψ(τ) =φ(τ) is the angular velocity of the pendulum. Note that according to the physical formulation of the problem, in inequalities ω > 0 and γ 0 are satisfied.

Investigation of the Stability of a Linearized System
Investigate the stability of the solutions of the system (18). To do this, we write out the characteristic equation for (18) (we consider the autonomous case): When converted to a third-degree algebraic equation with constant coefficients it takes the form: As is known, a necessary and sufficient condition for stability of a system of differential equations is the positivity of all three main diagonal minors ∆ 1 , ∆ 2 , and ∆ 3 of the Hurwitz determinant where a i , i ∈ {0, 1, 2, 3} are coefficients in (20); moreover, a 0 = 1 (see, for example, ref. [64]). Since a 3 = −ω 2 < 0, the conditions in Hurwitz criterion ∆ 1 > 0, ∆ 2 > 0, and ∆ 3 > 0, as it can be easily seen, are incompatible for any admissible values of the parameters ω, γ, and σ. Herewith, system (18) is unstable.Nevertheless, the introduction of a control with external force w(τ) makes it possible in some cases to bring the dynamical system into a stable oscillation mode. In the next section, we will consider at what values of the parameters of the problem this is possible.

Amplitude Control of the Pendulum Oscillation
In this section, we will show that the introduction of a control by external force w(τ) makes it possible to stabilize the motion of the pendulum and provides an opportunity to control amplitude ϕ(τ) of its oscillations.
In this section, we use the representation of the values varying according to the harmonic law in the form of the real part of the exponential function of the imaginary argument ∼ exp(iΩτ). It is known, for linear systems of algebraic or differential equations, this approach is correct and allows to perform intermediate calculations in a more compact and clear form.
Let us write the system of Equation (18) under control conditions in the formẇ(t) = αϕ(t) + d e iΩt , where d ∈ C: As can be seen from (24), the coefficients of the matrix have changed as a result of the control excitation. Check the conditions of the Hurwitz criterion in this case. Write out the characteristic equation and the Hurwitz matrix which leads to the following conditions: Therefore, system (24) is stable if and only if its parameters satisfy conditions (27). Its solution with respect to control parameter α has a simple form: where α b is the larger of two roots of (27) and D is its discriminant: Write down the solution to system (24), represented as a sum of the general solution of the homogeneous system and the partial solution of the inhomogeneous one: In (30), λ 1 , λ 2 , and λ 3 are the eigenvalues of matrix (24) and f 1 , f 2 , and f 3 are corresponding eigenvectors.
In a steady-state, i.e., at sufficiently large values of the time, τ max(1/|Re As can be seen from (31), the solution in the steady-state condition contains terms varying according to the harmonic law ∼ e iΩτ . Their amplitudes A ϕ , A ψ , and A v , as it can be easily shown by direct substitution of ϕ(τ) = A ϕ e iΩτ , ψ(τ) = A ψ e iΩτ , v(τ) = A v e iΩτ in (24), satisfy the system of algebraic equations: Therefore, for the amplitude of oscillations of the pendulum as a function of the amplitude of the harmonic term with frequency Ω in the controlẇ(τ) we get the expression: where F (Ω) = ω 2 − 2α(γ + σ) − iΩ (1 + 2γ)α + γ + σ − ω 2 + Ω 2 (α + γ + 1) + iΩ 3 is the determinant of (32) (note that this determinant can also be obtained from the right side of (25) by performing the formal substitution λ → iΩ). Steady-state oscillations are possible under condition F (Ω) = 0. The amplitudes of oscillations of the angular velocity of the pendulum and the voltage are calculated in the same way: The results of numerical simulations of the stability zones, roots of the characteristic equation, and solutions to the system at various parameter values are shown in Figures 2 and 3.
Dependencies of the amplitude |A ϕ | of the forced oscillations in the parameter space {(d, Ω)} are shown in Figures 4 and 5. These figures show the maximum amplitude decreasing with increasing γ, as well as the shift in the resonance region towards higher frequencies of Ω in this case. Due to the linearity of system (32), the value of |A ϕ | is proportional to the parameter d in the control functionẇ(t).

Hysteresis Dependencies in the Energy Harvester Model
In the practical implementation of the design of energy harvesting devices, as a rule, there are hysteresis dependencies of various kinds, both in the mechanical subsystem and of an electromagnetic nature. In this section, we investigate an electromechanical energy harvesting system with hysteresis damping, in other words, it is assumed that the friction in the mechanical subsystem obeys the hysteresis law (see Figure 6). Below, we will look at two approaches to the formalization of hysteresis friction. The first approach is based on the design model of the hysteresis operator (a Preisach operator), while the second one is described within the phenomenological Bouc-Wen model.

System with Hysteresis Friction within the Preisach Model
In Figure 6, a model of an energy harvester based on a classical oscillator m, connected to the electrical subsystem by means of a hysteresis damping link, is schematically presented. The external excitation was determined by periodic force F ext = K sin(Ωτ) applied to m.
In this section, we introduce a hysteresis operator using the approach developed in a book by M. A. Krasnoselskii and A. V. Pokrovskii. Within this approach, the hysteresis operator is interpreted as an operator defined on the space of continuous functions which dynamics are described by the relations: "input-state" and "state-output".
Denote the hysteresis converter corresponding to the non-ideal relay with threshold numbers α and β by R[α, β, x 0 , t 0 ], where x 0 ∈ {0, 1} is the initial state of the converter and t 0 is the initial time instant. The space of states of the non-ideal relay is a two-element set {0, 1}. The input of the system is continuous at t t 0 function u(t). The output is a a step function, x(t), defined by the relation: Note that initial state x 0 must satisfy the condition: In the case where inequalities α u(0) β are satisfied, value x 0 can achieve any value from set {0, 1}. The output values of x(t) with continuous input u(t) for t ∈ (t 0 , ∞) at each t = τ are determined following to the rule: Let us say that the relay is on if the output is equal to one, and that the relay is off otherwise.
The Preisach operator is a continuum analogue of a family of non-ideal relays connected in parallel (a block diagram of such a converter is shown in Figure 7). The state space of this operator consists of pairs {u(t), z(α, β, t)}, where u(t) is the input value at time instant t and z(α, β, t) is the characteristic function of the subset halfplane α < β, taking values 0, 1. The input-output relationships of the Preisach operator are defined by the relations: input-state where z 0 = z(α, β, t 0 ), and state-output The specified operator is widely used to formalize various hysteresis dependencies, its properties as well as its various applications can be found, for example, in [65].
Assume that the support of the measure of the Preisach operator is bounded, and therefore the space of states consists of characteristic functions, the support of which is contained in bounded sets.
The results of the numerical simulations of the system with hysteresis friction within the Preisach model where H is the coefficient characterizing the effect of hysteresis friction, are shown in Figure 6.      Figure 12 shows the hysteresis curve in coordinates (ϕ, Γϕ). As can be seen, after several oscillation cycles the established mode is reached. The values of transmitted power and the electrical power are presented in Table 1.

Bouc-Wen Model
Within the Bouc-Wen model the dynamics of the energy harvester with hysteresis friction will be described by (42): The numerical results for (42) are shown in Figures 13 and 14. The hysteresis curve in coordinates (ϕ, z) is shown in Figure 15. In Figure 16 phase portraits are given. Figure 17 shows a family of resonance curves for the mechanical part of system (42). Each of them exhibits the dependence of the amplitude of the solution on the frequency of the external periodic excitation Ω. We also show plots for various values of parameter H, which is responsible for the effect of the hysteresis block in the system. The results show that increasing this parameter significantly shifts the resonance frequency and also changes the shape and nature of the resonance curves depending on the amplitude of the forcing effect K. It is also worth noting a significant decrease in the amplitude of natural oscillations is observed on increasing the parameter H before the hysteresis block. The energy of the external excitation u(τ) = K sin(Ωτ) is equal to E ext =u 2 2 . Thus, the external excitation power can be calculated as: where T = 2π/Ω is the period of steady-state oscillations. Next, estimate the energy absorbed by the hysteresis block [66,67]. The equation connecting variables z and ϕ, has a single root, z 0 . Here, u max is the amplitude of the oscillatory process, u y = 1, The hysteresis loop area is expressed through the root of Equation (44) as follows: The power absorbed by the hysteresis link is equal to S T , where S is the area of the hysteresis loop in the space {(ϕ, z)}. As is mentioned above, the electrical power is equal to p avr = 1 T T 0 v 2 (τ)dτ. Tables 2 and 3 show the dependence of the power on the system parameters. The last column shows the relative power generated by the electrical part of the system with respect to the power of oscillations of the mechanical subsystem. These tables clearly show that the relative power depends both on the amplitude of the external excitation and on the value of the hysteresis term.

Conclusions
In the first part of this paper, we proposed a model for an energy harvester, the mechanical part of which is supposed to be an inverted pendulum. Within the proposed model, we investigated the stability of the linearized system. Particulary, it was proven that the stabilizing control of the pendulum, based on the principles of the feedback, enables the stabilization of the system. An important characteristic of the energy harvester is that the electrical power under the load is presented for the linearized system in a closed analytical form. During numerical experiments, we identified zones in the parameter space corresponding to the maximum power. Additionally, we identified the stability zones in the parameter space, the amplitude-frequency characteristics, as well as the dependence of the transmitted power on the coupling parameter.
The dynamics of the system in the case of the hysteresis friction in the transmission mechanical joint was given in the second part of this paper. The role of non-linear effects was shown within the design Preisach model and the phenomenological Bouc-Wen model.