Identification and nonlinearity compensation of hysteresis using NARX models
Petrus E. O. G. B. Abreu, Lucas A. Tavares, Bruno O. S. Teixeira, Luis A. Aguirre
IIdentification and nonlinearity compensation ofhysteresis using NARX models
Petrus E. O. G. B. Abreu , Lucas A. Tavares , Bruno O. S.Teixeira and Luis A. Aguirre Graduate Program in Electrical Engineering, Universidade Federal de MinasGerais, Av. Antˆonio Carlos 6627, 31270-901, Belo Horizonte, MG, Brazil Departament of Electronic Engineering, Universidade Federal de Minas Gerais,Av. Antˆonio Carlos 6627, 31270-901, Belo Horizonte, MG, BrazilE-mail: [email protected] , [email protected] , [email protected] and [email protected] December 2019
Abstract.
This paper deals with two problems: the identification andcompensation of hysteresis nonlinearity in dynamical systems using nonlinearpolynomial autoregressive models with exogenous inputs (NARX). First, basedon gray-box identification techniques, some constraints on the structure andparameters of NARX models are proposed to ensure that the identified modelsdisplay a key-feature of hysteresis. In addition, a more general frameworkis developed to explain how hysteresis occurs in such models. Second, twostrategies to design hysteresis compensators are presented. In one strategy thecompensation law is obtained through simple algebraic manipulations performedon the identified models. It has been found that the compensators based ongray-box models outperform the cases with models identified using black-boxtechniques. In the second strategy, the compensation law is directly identifiedfrom the data. Both numerical and experimental results are presented to illustratethe efficiency of the proposed procedures.
Keywords : Hysteresis, gray-box identification, compensation of nonlinearities, NARXmodel
1. Introduction
Hysteresis is a nonlinear behavior that is present in several systems and devices. It iscommonly related to the phenomena of ferromagnetism, plasticity, and friction, amongothers [1]. Some examples include mechanical, electronic and biomedical systems,as well as sensors and actuators such as magneto-rheological dampers, piezoelectricactuators and pneumatic control valves [2–4]. An intrinsic feature of such systemsis the memory effect, meaning that their output depends on the history of thecorresponding input.In addition to the memory effect, the literature provides different definitions andconditions to distinguish such systems and characterize the hysteretic behavior. Insome cases, the occurrence of hysteresis has been associated with the existence of a r X i v : . [ ee ss . S Y ] J a n dentification and nonlinearity compensation of hysteresis H formed by sets of stable equilibria and its implicationon the existence of the hysteresis loop in the identified models has also introducedin [6]. However, for more general cases, this concept and conditions need to beadapted. For instance, the conditions proposed in [6] are not sufficient to ensurethe existence of several equilibrium points at steady-state [5, 30]. Also, the conceptof bounding structure is limited to cases in which the sets of equilibria that form thisstructure are stable. Recently, a NARX model was identified for an experimentalelectronic circuit with hysteresis [29]. In this paper, we will propose ways to overcomesome shortcomings pointed out in the aforementioned references and a more generalframework will be put forward to explain how hysteresis takes place in identifiedmodels.The main contributions of this work are: the proposition of a specific parameterconstraint that ensures reproducing a key-feature of hysteresis through identified dentification and nonlinearity compensation of hysteresis
2. Background
A NARX model can be represented as [22]: y k = ˜ F (cid:0) y k − , · · · , y k − n y , u k − τ d , · · · , u k − n u (cid:1) , (1)where y k ∈ R is the output at instant k ∈ N , u k ∈ R is the input, n y and n u are themaximum lags for the output and input, respectively, τ d ∈ N + is the pure time delay,and ˜ F ( · ) is a nonlinear function of the lagged inputs and outputs.This work considers a linear-in-the-parameters extended model set [31] of theNARX model (1) with the addition of specific functions, such as absolute value,trigonometric, and sign function. The goal is to choose functions that allow themodels to predict systems whose nonlinearities cannot be well approximated using onlyregressors based on monomials of lagged input and output values. For instance, [31]recommends the addition of absolute value and sine functions as candidate regressorsfor the identification of a damped and forced nonlinear oscillator. In the case of theidentification of systems with hysteresis, [6] shows that including the regressor given bysign of the first difference of the input, i.e. sign( u k − u k − ), in addition to polynomialterms is a sufficient condition to reproduce hysteresis. Therefore, in this work themodels are of the general type: y k = F (cid:96) (cid:0) y k − , · · · , y k − n y , u k − , · · · , u k − n u ,φ , k − , φ , k − (cid:1) , (2)where φ , k = u k − u k − , φ , k =sign( φ , k ), and F (cid:96) ( · ) is a polynomial function of theregressor variables up to degree (cid:96) ∈ N + .Evaluating model (2) along a data set of length N , the resulting set of equationscan be expressed in matrix form as: y = Ψ ˆ θ + ξ , (3)where y (cid:44) [ y k y k − · · · y k +1 − N ] T ∈ R N is the vector of output measurements,Ψ (cid:44) [ ψ Tk − ; · · · ; ψ Tk − N ] ∈ R N × n θ is the matrix composed by measurements of theregressors vector ψ k − ∈ R n θ which contains linear and nonlinear combinations of dentification and nonlinearity compensation of hysteresis F (cid:96) ( · ) in (2) weighted by the parameter vector ˆ θ ∈ R n θ , ξ (cid:44) [ ξ k ξ k − · · · ξ k +1 − N ] T ∈ R N is the residual vector and T indicates the transpose.The unconstrained least squares batch estimator is given byˆ θ LS = (Ψ T Ψ) − Ψ T y . Assume the set of equality constraints on the parameter vector written as c = S θ ,where c ∈ R n c and S ∈ R n c × n θ are known constants. Then, the constrained leastsquares estimation problem is given byˆ θ CLS = arg min θ : c = S θ (cid:2) ξ T ξ (cid:3) , (4)whose solution is [32]:ˆ θ CLS = ˆ θ LS − (Ψ T Ψ) − S T [ S (Ψ T Ψ) − S T ] − ( S ˆ θ LS − c ) . (5)For the model structure selection, we use the error reduction ratio (ERR) [33]together with Akaike’s information criterion (AIC) [34]. Other approaches that haveproved to be useful in more demanding contexts are found in [35–39].
3. Identification of Systems with Hysteresis
The scientific community has been investigating which relevant features must bepresent in a model to reproduce hysteresis. Some of these features are: a characteristicloop behavior displayed on the input-output plane [30], several stable equilibriumpoints [5], and multi-valued mapping [26]. However, which and how these features canbe used in the identification procedure remains an open problem.Here, a constraint is proposed to ensure a key-feature of hysteresis. Also, itis shown how the hysteresis loop can be seen as an interplay of attracting andrepelling regions in certain models. Then, the resulting models will be used to designcompensators.First, a property of hysteresis based only on the main aspect discussed in [5,6,30]is presented. In the sequel, this property is used to obtain constraints on the structureand parameters of the model.
Property 1
An identified model of hysteresis, under a constant input, has two ormore real non-diverging equilibria. (cid:3)
In [6], Property 1 was attained by ensuring that the model had at least oneequilibrium point under loading-unloading inputs, with different values for loadingand unloading. Thus, in (2) φ , k = u k − u k − and φ , k =sign( φ , k ), with φ , k =1 forloading, and φ , k = − dentification and nonlinearity compensation of hysteresis By means of static analysis it is possible to determine the fixed points of a model.
Assumption 1
In order to comply with Property 1, considering the recommendationof the literature, the identified models should not have the following regressors:(i) y pk − τ y , y pk − τ y φ m , k − τ u and y pk − τ y φ m , k − τ u for p> , ∀ m [41],(ii) sign( u k − τ u − u k − τ u − ) m = φ m , k − τ u for m > y pk − τ y u mk − τ u , ∀ p, m ,where τ y and τ u are any time lags. (cid:3) The steady-state analysis of a model that complies with Assumption 1 is doneby taking y k = ¯ y, ∀ k , u k = ¯ u, ∀ k and, consequently, φ , k = u k − u k − = 0 φ , k = sign( φ , k ) = 0 , ∀ k , thus yielding ¯ y = Σ y ¯ y , where Σ y is the sum of allparameters of all linear output regressors. For Σ y (cid:54) = 1 the model has a single fixedpoint at ¯ y = 0 with stability domain given by: − < θ < . (6)If | θ | < | θ | > y = 0 is only non-diverging (diverging) equilibrium and,as a result, Property 1 is not satisfied. In order to solve this problem, we start bereviewing the following definition. Definition 1 (Continuum of equilibrium points [5]) . A model has a continuum ofequilibrium points if for any constant value of the input its corresponding output insteady-state is an equilibrium solution. (cid:3)
Based on Definition 1 and the problem aforementioned, the following lemma isstated.
Lemma 1
Given that Assumption 1 holds, if Σ y = 1 is verified, then the identifiedmodel has a continuum of equilibrium points at steady-state. (cid:3) Proof.
The steady-state analysis of a model that satisfies Assumption 1 andLemma 1 yields ¯ y = ¯ y which is trivially true for any value ¯ y . Hence, the model has a continuum of equilibrium points and Property 1 is satisfied. (cid:3) The core idea of the framework proposed in [6] to identify models with a hysteresisloop is to build a bounding structure H made of sets of equilibria and to ensure thatone set is stable during loading and the other one, during unloading. Such a scenario iseffective, but it does not help to understand models with more complicated structuresand with both attracting and repelling regions in the u × y plane. This section aimsat enlarging the scenario developed in [6].In quasi-static analysis, it is assumed that the input u k is a loading-unloadingsignal that is much slower than the system dynamics to the point that, at a given time k , the system will be in a certain attracting region, avoiding any possible repellingdentification and nonlinearity compensation of hysteresis u k , φ , k and φ , k . More specifically, there willbe two sets of regions, one for loading and another for unloading.In quasi-static analysis, we assume that y k ≈ y k − j = ˜ y, j = 1 , , . . . , n y , suchthat (2) is given by˜ y ≈ F (cid:96) (cid:0) ˜ y, u k − , · · · , u k − n u , φ , k − , φ , k − (cid:1) , (7)which can be usually solved for ˜ y , especially if higher powers of the output are notin F (cid:96) ( · ) [41]. This is achieved in practice by removing such group of terms from theset of candidates. If the model has no inputs, then ˜ y coincides with the fixed points.Alternatively, if the inputs are all constant, then ˜ y is a family of fixed points thatdepends on the set of constant inputs.Given the slow input, if ˜ y is in an attractive region, then the model outputmoves towards ˜ y . In what follows, ˜ y aL and ˜ y aU are, respectively, the solutions to (7)in attracting regions under loading and unloading. Likewise, ˜ y rL and ˜ y rU are theircounterparts in repelling regions. The conditions for ˜ y to be attracting is (cid:12)(cid:12)(cid:12)(cid:12) eig (cid:18) ∂F (cid:96) ( y , u k − , φ , k − , φ , k − ) ∂ y (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) < , (8)where y = [ y k − . . . y k − n y ] T . This procedure resembles that of determining thestability of fixed points. Here the Jacobian matrix is not evaluated at fixed points.Hence we do not speak in terms of stable and unstable fixed points.To illustrate how this helps to understand the formation of a hysteresis loop,consider the schematic representation in figure 1. The input is a loading-unloadingsignal such that u min ≤ u k ≤ u max , ∀ k . The sets ˜ y aL , ˜ y aU , ˜ y rL and ˜ y rU are shown.Consider the point A, which takes place under loading. Given that the system isunder the direct influence of ˜ y rL , which is responsible for pushing upwards (see verticalcomponent y A ), and it is the loading regime, there is a horizontal component u A (related to the input) that points to the right. The resulting effect is to pull the systemalong the loop in the NE direction. The same can be said for point B; however, atthat point the vertical component is the result of the attracting action of ˜ y aL . A similaranalysis can be readily done for the unloading regime, given by points D and E. Atthe turning points C and F, φ , k switches from 1 to -1 and from -1 to 1, respectively.Hence the analysis also switches from using ˜ y aL and ˜ y rL , to using ˜ y aU and ˜ y rU . Thisanalysis will be useful in section 5 to understand the formation of hysteresis loops inidentified models.As a final remark, it is important to point out that the assumption that the set˜ y comes in two disjoint parts, either for loading or unloading, is a consequence of thesolution of (7) being rational instead of polynomial. This is useful to analyse modelswith more general model structures.The following example illustrates the application of this analysis and show whichconstraints should be considered to comply with Property 1. Example 1
Consider the following NARX model that complies with Assumption 1: y k = θ y k − + θ φ , k − + θ φ , k − u k − + θ φ , k − φ , k − y k − + θ φ , k − . (9)In this case, the constraint θ = 1 will be used such that, according to Lemma 1, theresulting model will have a continuum of equilibrium points. This can be achievedusing estimator (5) with c = 1 and S = [1 0 0 0 0]. dentification and nonlinearity compensation of hysteresis Smart Mater. Struct. ˜ y aL y uu max u min ˜ y aU ˜ y rU ˜ y rL DE CF A B y A u A u B y B u D y D u E y E U LLU
Figure 1.
Schematic representation of hysteresis loop in the u × y plane. Attracting sets are shown in black continuouslines, whereas the repelling sets are indicated in dash-dot. Thehysteresis loop is indicated by dotted lines. Remark 1 test
Figure 1.
Schematic representation of hysteresis loop in the u × y plane.Attracting sets are shown in black continuous lines, whereas the repelling setsare indicated in red dash-dot. The hysteresis loop is indicated by dotted lines. For a more complicated model structure, the constraint in Lemma 1 is still inthe form 1 = S θ (4) but with S having more that one element equal to one, e.g. asshown in [42] to obtain NARX models able to reproduce dead-zone and in [43] for aquadratic nonlinearity.The quasi-static analysis of model (9) is performed following the steps providedin section 3.2. So rewriting this model as (7), we have˜ y ≈ θ ˜ y + θ φ ,k − + θ φ ,k − u k − + θ φ ,k − φ ,k − ˜ y + θ φ ,k − , which can be described by˜ y ( u, φ , φ )= θ + θ φ u + θ φ − θ − θ φ , for φ = 1; − θ + θ φ u + θ φ − θ + θ φ , for φ = − , (10)where the time indices have been omitted for simplicity. Therefore, the solution givenat the top in (10) represents the set ˜ y L , while the bottom is the set ˜ y U .To define whether the solutions to (10) are in the attracting or repelling regions,(8) should be computed for model (9) as − < θ + θ φ ,k − φ ,k − < , − − θ θ φ ,k − < φ ,k − < − θ θ φ ,k − . (11)Since it is assumed that the known input u k is a non-zero loading-unloading signal,then the conditions (11) to ensure that the solutions (10) are in attracting regionscan be readily verified numerically. In sections 5 and 6, the same analyses will beperformed for the identified models. (cid:3) dentification and nonlinearity compensation of hysteresis
4. Compensator Design
The proposed strategies to design compensators based on NARX models are detailedin this section, starting with some preliminary assumptions. In this paper, a key pointis to investigate if hysteresis in the models estimated according to section 3 have anyimpact on the regulation performance.
Given a nonlinear system S , the first step is to obtain hysteretic models for S ; seefigure 2(a). To achieve that, two procedures will be followed. The first one aimsat identifying a model M based on the direct relationship between u and y , whosesimulation yields ˆ y k , according to section 4.2. The second procedure is based on theidentification of the inverse relationship, in which case a model ˘ M is obtained toyield ˆ u k , as illustrated in figure 2(a) and following section 4.3. In the second step,the identified model is used to design a compensator C that yields the compensationsignal m k for a given reference r k ; see figure 2(b). JOURNAL OF L A TEX CLASS FILES, VOL. 14, NO. 8, FEB 2020 4 ˜ y aL y uu max u min ˜ y aU ˜ y rU ˜ y rL DE CF A B y A u A u B y B u D y D u E y E U LLUFig. 1. Schematic representation of hysteresis loop in the u × y plane.Attracting sets are shown in black continuous lines, whereas the repelling setsare indicated in dash-dot. The hysteresis loop is indicated by dotted lines. As a final remark, it is pointed out that the assumption thatthe set ˜ y (either for loading or unloading) comes in two disjointparts is a consequence of the solution of (7) being rationalinstead of polynomial. This is useful to analyse models withmore general model structures.The following example illustrates the application of thisanalysis and define which constraints should be consideredto comply with Definition 1. Example 1.
Consider a NARX model that complies withRemark 2 y k = θ y k − + θ φ , k − + θ φ , k − u k − + θ φ , k − φ , k − y k − + θ φ , k − . (9)In this case, the extra constraint required by Proposition 1 issimply θ = 1 and it can be imposed using estimator (4) with c = 1 and S = [1 0 0 0 0] .For a more complicated model structure, the constraint inProposition 1 is still in the form S θ (3) but with S havingmore that one element equal to one, e.g. as shown in [39] toobtain NARX models able to reproduce dead-zone.The quasi-static analysis of model (9) is performed follow-ing the steps provided in Sec. III-B. So rewriting this modelas Eq. (7), we have ˜ y ≈ θ ˜ y + θ φ ,k − + θ φ ,k − u k − + θ φ ,k − φ ,k − ˜ y + θ φ ,k − , which can be described by ˜ y ( u, φ , φ )= θ + θ φ u + θ φ − θ − θ φ , for φ = 1; − θ + θ φ u + θ φ − θ + θ φ , for φ = − , (10)where the time indices have been omitted for the sake ofsimplicity. Therefore, the solution given at the top in Eq. (10)represents the set ˜ y L while the bottom is the set ˜ y U . To define if the solutions to (10) are in the attracting orrepelling regions, Eq. (8) should be computed for model (9): − < θ + θ φ ,k − φ ,k − < , − − θ θ φ ,k − < φ ,k − < − θ θ φ ,k − . (11)As for this analysis it is assumed that the known input u k is a non-zero loading-unloading signal, the conditions (11) toensure that the solutions (10) are in attracting regions can bereadily verified numerically. In Sections V and VI, these sameanalyses will be performed for the identified models.IV. D ESIGN OF C OMPENSATORS
The proposed strategies to design compensators based onNARX models are detailed in this section, starting with somepreliminary assumptions. A key point in this paper is to investi-gate if hysteresis in the models estimated following the resultsin Sec.III have any impact on the regulation performance.
A. Preliminaries
Given a nonlinear system S , the first step is to obtainhysteretic models for S (see Fig. 2(a)). In doing this, twoprocedures will be followed. The first (Sec. IV-B) aims toidentify a model M based on the direct relationship between u and y , that is, the model will produce ˆ y k . The second(Sec. IV-C) is based on the identification of the inverserelationship, in which case a model ˘ M is obtained to yield ˆ u k , as illustrated in Fig. 2(a). InputOutputData CollectedIdentificationIdentification Procedure Identified NARX ModelsIdentified ModelDesign C Design of Compensator System Compensation M / ˘ M (a) u k Input ˆ y k u S System y OutputCompensationInputReference (b) M Output y k Input ˆ u k ˘ M Output u Input S System Output y r k C Compensator m k InputOutput M ˘ M Fig. 2. Compensator design based on identified NARX models. (a) Modelidentification; (b) compensator design based on identified models.
In the second step, the identified model is used to design acompensator C that will produce the compensation signal m k for a given reference r k (see Fig. 2(b)).In this paper the NARX models use the extended represen-tation and the following additional assumptions are made. Remark 3.
In the models M and ˘ M , y k can be replaced with r k and u k with m k for compensation design. The motivationbehind this is that y k should ideally be equal to r k undercompensation, that is, when m k is used as an input to thesystem. Figure 2.
Compensator design based on identified NARX models. (a) Modelidentification, and (b) compensator design based on identified models.
In this paper,the following additional assumptions are made for NARX models (2).
Remark 1
For compensation design, y k can be replaced by r k , and u k by m k ,respectively, in the models M and ˘ M . (cid:3) The motivation behind this is that y k should ideally be equal to r k undercompensation, that is, when m k is used as an input to the dynamical system.In what follows, the main idea is to use an identified model to determine thecompensation input m k − τ d +1 . The aim here is to specify a general model structure for M in order to find m k analytically from this model. To achieve that, the following assumptions are needed. dentification and nonlinearity compensation of hysteresis Assumption 2
Assume that: (i) the only regressor involving u k − τ d is linear; (ii) n u >τ d ; (iii) the compensation signal m k is known up to time k − τ d ; and (iv) the referencesignal r k is known up to time k + 1. (cid:3) Assumption 2 imposes conditions on the selection of the model structure. Notethat (i) ensures that u k − τ d can be isolated in the identified models; (ii) allows thatinput terms with a delay longer than τ d to be regressors in the identification procedure;and the other constraints guarantee that the control action can be computed fromknown values. Therefore, the model M is rewritten as A ( q ) y k = B ( q ) u k + f (cid:0) y k − , · · · , y k − n y , u k − τ d − , · · · , u k − n u (cid:1) , (12)where q − is the backward time-shift operator such that q − u k = u k − , and the linearregressors are grouped in A ( q ) y k and B ( q ) u k with A ( q ) = 1 − a q − − a q − − · · · − a n y q − n y , (13) B ( q ) = b τ d q − τ d + b τ d +1 q − τ d − + · · · + b n u q − n u (cid:124) (cid:123)(cid:122) (cid:125) B ∗ ( q ) , (14)and f ( · ) includes all the nonlinear terms and possibly the constant term of the NARXmodel (2). Using (14), model (12) can be rewritten as A ( q ) y k = b τ d u k − τ d + B ∗ ( q ) u k + f (cid:0) y k − , · · · , y k − n y ,u k − τ d − , · · · , u k − n u (cid:1) . (15)From Remark 1, we have A ( q ) r k +1 = b τ d m k − τ d +1 + B ∗ ( q ) m k +1 + f (cid:0) r k , · · · ,r k − n y +1 , m k − τ d , · · · , m k − n u +1 (cid:1) , (16)which, for convenience, has been written an instant of time ahead, i.e. k → k + 1.From the Assumption 2, the compensation input can be obtained from (16) as m k − τ d +1 = 1 b τ d (cid:104) A ( q ) r k +1 − B ∗ ( q ) m k +1 − f (cid:0) r k , · · · ,r k − n y +1 , m k − τ d , · · · , m k − n u +1 (cid:1)(cid:105) . (17)To illustrate the application of this strategy, assume that the constraints discussedin section 3 and Assumption 2 are verified in the identification procedure. Example 2
Consider the NARX model described by y k = θ y k − + θ φ , k − + θ φ , k − u k − + θ φ , k − φ , k − y k − + θ φ , k − . (18)Since φ , k = u k − u k − and φ , k =sign( φ , k ), we have y k = θ y k − + θ sign( u k − − u k − )+ θ [ u k − − u k − ] u k − + θ sign( u k − − u k − )[ u k − − u k − ] y k − + θ [ u k − − u k − ] , which is in the form (12) and, therefore, A ( q ) y k = B ( q ) u k + f (cid:0) y k − , u k − , u k − , sign( u k − − u k − ) (cid:1) , (19) dentification and nonlinearity compensation of hysteresis A ( q ) = 1 − θ q − , (20) B ( q ) = θ q − − θ q − , (21) f ( · ) = θ sign( u k − − u k − )+ θ [ u k − − u k − ] u k − + θ sign( u k − − u k − )[ u k − − u k − ] y k − . (22)From Remark 1, the model (19) is recast as A ( q ) r k +1 = θ m k − θ m k − + f (cid:0) r k , m k − , m k − , sign( m k − − m k − ) (cid:1) , (23)and m k = 1 θ (cid:104) A ( q ) r k +1 + θ m k − − f (cid:0) r k , m k − , m k − , sign( m k − − m k − ) (cid:1)(cid:105) , = 1 θ (cid:104) r k +1 − θ r k + θ m k − − θ sign( m k − − m k − ) − θ [ m k − − m k − ] m k − − θ sign( m k − − m k − )[ m k − − m k − ] r k (cid:105) , (24)which is computed due to Assumption 2. (cid:3) Here, the strategy is to identify NARX models ˘ M that are able to describe the inverserelationship between the input u and output y signals of the system S . The advantageof this strategy is that the compensator C is obtained directly from ˘ M , according tothe variable solutions presented in Remark 1. However, some issues related to theidentification procedure of these models need to be addressed. For nomenclaturesimplicity, in this section, we assume that τ d = 1.For the inverse model ˘ M , the output ˆ u k depends on y k . Hence in order toavoid the lack of causality, y k should be delayed by τ s time steps with respect to u k ,yielding [44]: ˆ u k = ˘ F (cid:0) ˆ u k − , · · · , ˆ u k − n u , y k − τ s , · · · , y k − n y + τ s (cid:1) , (25)where ˘ F ( · ) is the inverse nonlinear function and ˆ u k ∈ R and y k ∈ R are relatedas shown in figure 2(a). It should be noted that τ s ≥ τ d + 1, where usually theequality is preferred. Similar ways to avoid noncausal models can be found in theliterature [29, 45]. Assumption 3
Assume that: (i) there is at least one regressor of the output ( y k ) j for j ≥
1; (ii) the compensation signal m k is known up to time k −
1; and (iii) thereference signal r k is known up to time k − τ s . (cid:3) Assumption 3 should be imposed during the structure selection of the inversemodel ˘ M . Note that (i) ensures that there is at least one input signal y k in theidentified models; (ii) and (iii) ensure that the compensation input m k to be computedat time k is the only unknown variable. Given Assumption 3 and Remark 1, thecompensation signal m k can be obtained directly from ˘ M as m k = ˘ F (cid:0) m k − , · · · , m k − n u , r k − τ s , · · · , r k − n y + τ s (cid:1) . (26) dentification and nonlinearity compensation of hysteresis
5. Numerical Results
This section identifies models to predict the behavior of a hysteretic system fromsimulated data, and evaluates the performance of these models in predicting dynamicsand compensating the nonlinearity of the simulated system.
Consider the piezoelectric actuator with hysteretic nonlinearity modeled by the Bouc-Wen model [12] and whose mathematical model is given by [45] (cid:40) ˙ h ( t ) = A ˙ u ( t ) − β | ˙ u ( t ) | h ( t ) − γ ˙ u ( t ) | h ( t ) | ,y ( t ) = d p u ( t ) − h ( t ) , (27)where y ( t ) is the displacement, u ( t ) is the voltage applied to the actuator, d p = 1 . µ mV is the piezoelectric coefficient, h ( t ) is the hysteretic nonlinear term and A = 0 . µ mV , β = 0 .
008 V − and γ = 0 .
008 V − are parameters that determine the shape and scaleof the hysteresis loop.Model (27) was integrated numerically using a fourth-order Runge-Kutta methodwith integration step δt = 0 .
001 s. The excitation signal was generated by low-passfiltering a white Gaussian noise [6]. In this work, a fifth-order low-pass Butterworthfilter with a cutoff frequency of 1 Hz was used; see figure 3(a). The sampling time isset to T s = δt = 0 .
001 s and a frequency of 1 Hz is chosen to validate the identifiedmodels [45]. The data sets are 50 s long ( N = 50000).The constraints defined in section 3 are here considered to build NARX modelsfor system (27). Input and output signals are shown in figure 3. In additionto the monomial regressors in u k and y k , the following regressors are also used: φ , k = u k − u k − and the sign of this first difference φ , k = sign( φ , k ). Themaximum nonlinear degree for regressors is cubic, (cid:96) = 3, and the maximum delaysare n y = n u = 1. This choice is based on the fact that discrete models of hysteresisthat have only unit delayed regressors typically perform well and result in models withsimple structures, which are advantageous for model-based control [6, 29].The model structure is selected using the ERR criterion to rank the regressorsaccording to importance and the AIC determine the final number of model terms. Inthis case, the standard least squares solution is used for parameter estimation. Also,as proposed in section 3, to obtain models that describe some features of hysteresis,the constrained parameter estimation was used in order to comply with the conditionestablished in Property 1. M . In this example we take the following metaparameters: n u = 2which is the smallest value that complies with Assumption 2-(ii); while n y and (cid:96) maintain the values which were determined above. Using the data shown in figure 3and considering Assumption 1, the following model structure is obtained y k = θ y k − + θ φ , k − + θ φ , k − φ , k − u k − + θ φ , k − φ , k − y k − + θ φ , k − u k − + θ φ , k − u k − y k − , (28)where θ i for i = 1 , · · · , y k = ¯ y, ∀ k , φ , k = 0 , ∀ k ; hence, the resulting expressionis ¯ y = θ ¯ y . Therefore, based on Lemma 1 and Example 1, for model (28) to fulfill dentification and nonlinearity compensation of hysteresis Smart Mater. Struct. y ( µ m ) t (s)(a) u ( V ) t (s)(b) Figure 1.
Signals used to identify system ( ?? ). (a) excitation,and (b) simulated output.-80 -60 -40 -20 0 20 40 60 80 u (V)-250-200-150-100-50050100150200250 y ( µ m ) Figure 2.
Results for model ( ?? ) with input u k =70 sin(2 πk ) V.The hysteresis loop indicated with ( · · · ) is a result of theinteraction of (—) attracting (˜ y aL , ˜ y aU ) and (- · -) repelling (˜ y rL ,˜ y rU ) sets. ( A ) indicates the orientation of the hysteresis loop. 5 5.5 6 6.5 7-40-2002040 5 5.5 6 6.5 7-40-2002040-40 -20 0 20 40-40-2002040 5.5 6 6.5 7-40-20020405.5 6 6.5 7-40-2002040-40 -20 0 20 40-40-2002040 u ( V ) t (s)(a) y ( µ m ) y ( µ m ) t (s)(b) t (s)(c) u (V)(e) t (s)(d) u (V)(f) u ( V ) y ( µ m ) y ( µ m ) Figure 3.
Free-run simulation of model ( ?? ). This figure isarranged in columns, which have: (a) sinusoidal input of voltage u k =40 sin(2 πk ) V and in (b) the case where this input becomesconstant during a loading ( • ) and unloading ( (cid:7) ) regime with thefinal value of 16 . N = 50000 data points. Figure 3.
Signals used to identify system (27). (a) excitation, and (b) simulatedoutput.
Property 1, the constraint Σ y = θ = 1 should be imposed. This can be done using(5) with the constraint written as: c = 1; S =[1 0 0 0 0 0] . (29)Hence, the parameter values estimated by the constrained least squares estimator (5)are shown in Table 1. Table 1.
Model parameters obtained with (5) and (29).
Model
Values(28) θ =1 . θ =0 . θ =1 . × − θ = − . × − θ =3 . × − θ = − . × − (33) θ =1 . θ =1 . θ = − . × − θ =1 . × − θ = − . × − θ =7 . × − Next, a quasi-static analysis of the identified model is performed as discussed insection 3.2 and illustrated in Example 1. First, we write for (28) the corresponding to(7) as ˜ y ≈ θ ˜ y + θ φ ,k − + θ φ ,k − φ ,k − u k − + θ φ ,k − φ ,k − ˜ y + θ φ ,k − u k − + θ φ ,k − u k − ˜ y, yielding ˜ y ( u, φ , φ )= θ φ + θ φ u + θ φ u − θ − θ φ − θ φ u , for φ = 1; θ φ − θ φ u + θ φ u − θ + θ φ − θ φ u , for φ = − , (30)where the time indices have been omitted for brevity.The top expression in (30) gives the set ˜ y L , while the bottom one, ˜ y U . Computingthe derivative of (28) with respect to y k − and using (8), we obtain − < θ + θ φ ,k − φ ,k − + θ φ ,k − u k − < , − − θ − θ φ ,k − φ ,k − θ φ ,k − < u k − < − θ − θ φ ,k − φ ,k − θ φ ,k − . (31) dentification and nonlinearity compensation of hysteresis φ ,k − =1 or φ ,k − = −
1, the conditions for attracting regions under loador unloading, respectively, are obtained. Considering the parameter values presentedin Table 1 and a loading-unloading input signal, the points (30) and their attractionconditions (31) are computed numerically and shown in figure 4. Figure 4 should becompared to figure 1, whose main elements are analogous. Hence, in this way it ispossible to see how model (28) is able to describe the hysteresis nonlinearity.
Smart Mater. Struct. y ( µ m ) t (s)(a) u ( V ) t (s)(b) Figure 1.
Signals used to identify system ( ?? ). (a) excitation,and (b) simulated output.-80 -60 -40 -20 0 20 40 60 80 u (V)-250-200-150-100-50050100150200250 y ( µ m ) Figure 2.
Results for model ( ?? ) with input u k =70 sin(2 πk ) V.The hysteresis loop indicated with ( · · · ) is a result of theinteraction of (—) attracting (˜ y aL , ˜ y aU ) and (- · -) repelling (˜ y rL ,˜ y rU ) sets. ( A ) indicates the orientation of the hysteresis loop. 5 5.5 6 6.5 7-40-2002040 5 5.5 6 6.5 7-40-2002040-40 -20 0 20 40-40-2002040 5.5 6 6.5 7-40-20020405.5 6 6.5 7-40-2002040-40 -20 0 20 40-40-2002040 u ( V ) t (s)(a) y ( µ m ) y ( µ m ) t (s)(b) t (s)(c) u (V)(e) t (s)(d) u (V)(f) u ( V ) y ( µ m ) y ( µ m ) Figure 3.
Free-run simulation of model ( ?? ). This figure isarranged in columns, which have: (a) sinusoidal input of voltage u k =40 sin(2 πk ) V and in (b) the case where this input becomesconstant during a loading ( • ) and unloading ( (cid:7) ) regime with thefinal value of 16 . N = 50000 data points. Figure 4.
Results of quasi-static analysis for model (28) with input u k =70 sin(2 πk ) V. The hysteresis loop indicated with ( · · · ) is a result of theinteraction of (—) attracting (˜ y aL , ˜ y aU ) and (- · -) repelling (˜ y rL , ˜ y rU ) sets. ( (cid:65) )indicates the orientation of the hysteresis loop. Model (28) is simulated with a loading-unloading input (see left side of figure 5)and, in cases where the input becomes constant, either during loading or unloading (seeright side of figure 5), the system remains at the corresponding point of the hysteresisloop.
This is a direct consequence of using Lemma 1. This feature is not generallypresent in identified models found in the literature.The mean absolute percentage error (MAPE)MAPE = 100 (cid:80) Nk =1 | y k − ˆ y k | N | max( y ) − min( y ) | , (32)computed for the case in figure 5(c), is shown in Table 2. ˘ M . The identified model that complies with Assumptions 1 and 3is given by ˆ u k = θ ˆ u k − + θ ˘ φ , k − + θ ˘ φ , k − ˘ φ , k − ˆ u k − + θ ˘ φ , k − ˘ φ , k − y k − + θ ˘ φ , k − y k − ˆ u k − + θ ˘ φ , k − y k − , (33)where ˘ φ , k = y k − y k − , ˘ φ , k = sign( ˘ φ , k ), ˆ u k is the estimated input (model output),and y k is the output of system (27) (model input).Note that the regressors of (28) and of (33) are different. In both cases, theregressors are automatically chosen from the pool of candidates using the ERRcriterion. Nevertheless, also for (33), the steady-state analysis yields ¯ˆ u = θ ¯ˆ u , which issimilar to the result found for model (28). Proceeding as before, the constrained leastsquares estimated parameters are shown in Table 1. dentification and nonlinearity compensation of hysteresis Smart Mater. Struct. y ( µ m ) t (s)(a) u ( V ) t (s)(b) Figure 1.
Signals used to identify system ( ?? ). (a) excitation,and (b) simulated output.-80 -60 -40 -20 0 20 40 60 80 u (V)-250-200-150-100-50050100150200250 y ( µ m ) Figure 2.
Results for model ( ?? ) with input u k =70 sin(2 πk ) V.The hysteresis loop indicated with ( · · · ) is a result of theinteraction of (—) attracting (˜ y aL , ˜ y aU ) and (- · -) repelling (˜ y rL ,˜ y rU ) sets. ( A ) indicates the orientation of the hysteresis loop. 5 5.5 6 6.5 7-40-2002040 5 5.5 6 6.5 7-40-2002040-40 -20 0 20 40-40-2002040 5.5 6 6.5 7-40-20020405.5 6 6.5 7-40-2002040-40 -20 0 20 40-40-2002040 u ( V ) t (s)(a) y ( µ m ) y ( µ m ) t (s)(b) t (s)(c) u (V)(e) t (s)(d) u (V)(f) u ( V ) y ( µ m ) y ( µ m ) Figure 3.
Free-run simulation of model ( ?? ). This figure isarranged in columns, which have: (a) sinusoidal input of voltage u k =40 sin(2 πk ) V and in (b) the case where this input becomesconstant during a loading ( • ) and unloading ( (cid:7) ) regime with thefinal value of 16 . N = 50000 data points. Figure 5.
Free-run simulation of model (28). This figure is arranged in columns,which have: (a) sinusoidal input of voltage u k =40 sin(2 πk ) V and in (b) the casewhere this input becomes constant during a loading ( • ) and unloading ( (cid:7) ) regimewith the final value of 16 . N = 50000 data points. Consider now the quasi-static analysis of model (33). The formation of thehysteresis loop for this model is shown in figure 6. Interestingly, the ability of model(28) to describe hysteresis is also present in model (33). The main difference betweenthem is the orientation of the hysteresis loop, as discussed in [46] and illustrated infigures 4 and 6.
Smart Mater. Struct. -80 -60 -40 -20 0 20 40 60 80 y ( µ m)-250-200-150-100-50050100150200250 u ( V ) Figure 1.
Results for model ( ?? ). For meaning of line patternsrefer to captions of Figure ?? and of Figure ?? .5 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7-500505 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7-50050-40 -30 -20 -10 0 10 20 30 40-40-2002040 m ( V ) y ( µ m ) y ( µ m ) t (s)(a) t (s)(b) r ( µ m)(c)17.5 18 18.51819 5.2 5.24 5.28363840 5.2 5.24 5.28424446 Figure 2.
Hysteresis compensation for the piezoelectricactuator ( ?? ). (a) Compensation inputs, (b) temporal responsesand in (c) hysteresis loops. (- -) results obtained withcompensator ( ?? ) ( · · · ) results with compensator ( ?? ), (- · -) system output without compensation, and (—) displacementreference r = 40 sin(2 πt ) µ m. 0 5 10 15 20 T s (s)012345678 M A PE ( % ) M A PE ( % ) × − (a) T s (s) × − (b) Figure 3.
MAPE index ( ?? ) computed for the models andcompensators described, respectively, by equations (a) ( ?? ) and( ?? ); (b) ( ?? ) and ( ?? ). ( ◦ ) model and ( • ) tracking accuracies.( N ) accuracy of uncompensated system. β = 0 . β = 0 . β = 0 . · · · · · · Figure 4.
Bouc-Wen hysteresis loops within the investigatedrange.50 55 60 65 70 t (s)2.533.5 u ( V )
50 55 60 65 702.833.23.450 55 60 65 702.833.23.4 50 55 60 65 702.533.52.5 3 3.52.833.23.4 2.8 3 3.2 3.42.533.5 y ( V ) y ( V ) t (s) y ( V ) u ( V ) u ( V ) t (s) u (V) t (s) y (V)(a) (b)(c) (d)(e) (f) Figure 5.
Left column refers to model ( ?? ) and right columnto model ( ?? ). (a) input u k =0 .
56 sin(0 . πk ) + 3 V and (c) thecorresponding measured output (—) y and (- -) model ( ?? )free-run simulation; (b) smoothed version of y in (c); (d) thecorresponding output which is u k in (a) and (- -) model ( ?? )free-run simulation. (e) and (f) show the same data as (c) and(d), respectively. Figure 6.
Results of quasi-static analysis for model (33). For meaning of linepatterns refer to captions of figure 1 and of figure 4.
Table 2 shows the prediction performance of models (28) and (33) are similar. Inaddition, Table 2 reports the prediction performance of a black-box NARX polynomial dentification and nonlinearity compensation of hysteresis
Table 2.
Performance of the modeling step. Simulation results.
Design Strategy Model
MAPESection 4.2 (28) 0 . . . Next, the models identified in the previous section are used to design compensatorsusing the procedure illustrated in figure 2(b).
Applying the steps described insection 4.2 to model (28), the following compensation signal is obtained m k = 1 θ (cid:104) r k +1 − θ r k + θ m k − − [ θ m k − + θ r k ]sign( m k − − m k − )[ m k − − m k − ] − [ θ m k − + θ m k − r k ][ m k − − m k − ] (cid:105) . (34)Similarly, using the method described in section 4.3 to identify the inverse model(33), the following compensation signal is obtained˘ m k = θ ˘ m k − + θ [ r k +1 − r k ]+ [ θ ˘ m k − + θ r k +1 ]sign( r k +1 − r k )[ r k +1 − r k ]+ [ θ r k +1 ˘ m k − + θ r k +1 ]sign( r k +1 − r k ) . (35)Since the parameters of compensators (34) and (35) have been estimated(Table 1), and Assumptions 2 and 3 are satisfied, the compensation inputs m k and˘ m k can be computed. The results are summarizedin figure 7(b). From the hysteresis loops shown in figure 7(c), it is clear that thecompensators enforced a practically linear relation between the reference and theoutput. This would greatly facilitate the design and increase the performance of afeedback controller.The accuracy achieved by each compensator was quantified by the MAPE index(32). In order to quantify the compensation effort, the normalized sum of the absolutevariation of the input (NSAVI)NSAVI = N − (cid:88) k =1 (cid:12)(cid:12) m k +1 − m k (cid:12)(cid:12)(cid:12)(cid:12) r k +1 − r k (cid:12)(cid:12) , (36)is calculated. These indices are shown in Table 3.The results shown in figure 7 and Table 3 indicate that the compensators mayprovide a significant improvement in the tracking performance of system (27). Thetracking error was reduced by about 93% at the cost of a 14% increase in the dentification and nonlinearity compensation of hysteresis Smart Mater. Struct. -80 -60 -40 -20 0 20 40 60 80 y ( µ m)-250-200-150-100-50050100150200250 u ( V ) Figure 1.
Results for model ( ?? ). For meaning of line patternsrefer to captions of Figure ?? and of Figure ?? .5 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7-500505 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7-50050-40 -30 -20 -10 0 10 20 30 40-40-2002040 m ( V ) y ( µ m ) y ( µ m ) t (s)(a) t (s)(b) r ( µ m)(c)17.5 18 18.51819 5.2 5.24 5.28363840 5.2 5.24 5.28424446 Figure 2.
Hysteresis compensation for the piezoelectricactuator ( ?? ). (a) Compensation inputs, (b) temporal responsesand in (c) hysteresis loops. (- -) results obtained withcompensator ( ?? ) ( · · · ) results with compensator ( ?? ), (- · -) system output without compensation, and (—) displacementreference r = 40 sin(2 πt ) µ m. 0 5 10 15 20 T s (s)012345678 M A PE ( % ) M A PE ( % ) × − (a) T s (s) × − (b) Figure 3.
MAPE index ( ?? ) computed for the models andcompensators described, respectively, by equations (a) ( ?? ) and( ?? ); (b) ( ?? ) and ( ?? ). ( ◦ ) model and ( • ) tracking accuracies.( N ) accuracy of uncompensated system. β = 0 . β = 0 . β = 0 . · · · · · · Figure 4.
Bouc-Wen hysteresis loops within the investigatedrange.50 55 60 65 70 t (s)2.533.5 u ( V )
50 55 60 65 702.833.23.450 55 60 65 702.833.23.4 50 55 60 65 702.533.52.5 3 3.52.833.23.4 2.8 3 3.2 3.42.533.5 y ( V ) y ( V ) t (s) y ( V ) u ( V ) u ( V ) t (s) u (V) t (s) y (V)(a) (b)(c) (d)(e) (f) Figure 5.
Left column refers to model ( ?? ) and right columnto model ( ?? ). (a) input u k =0 .
56 sin(0 . πk ) + 3 V and (c) thecorresponding measured output (—) y and (- -) model ( ?? )free-run simulation; (b) smoothed version of y in (c); (d) thecorresponding output which is u k in (a) and (- -) model ( ?? )free-run simulation. (e) and (f) show the same data as (c) and(d), respectively. Figure 7.
Hysteresis compensation for the piezoelectric actuator (27).(a) Compensation inputs, (b) temporal responses and in (c) hysteresis loops.(- -) results obtained with compensator (34) ( · · · ) results with compensator (35),(- · -) system output without compensation, and (—) displacement reference r = 40 sin(2 πt ) µ m. Table 3.
Performance of the compensation step. Simulation results.
Design Strategy Compensator
MAPE NSAVISection 4.2 (34) 0 .
322 1 . .
425 1 . .
819 1 . .
536 1 . T s is also investigated. In figure 8, it can be seen that the model dentification and nonlinearity compensation of hysteresis T s is increased. It should be noted that even thelargest values of T s in figure 8 are still comfortably small in terms of the samplingtheorem. However, since one of the regressors is the first difference of the input,then the identification of systems with hysteresis seems to be particularly sensitive tothe sampling time [47]. Another conclusion that can be drawn from figure 8 is that,for both design strategies, the compensation performance is correlated to the modelaccuracy, and that the strategy in section 4.2 (figure 8(a)) is somewhat less sensitiveto such accuracy. Smart Mater. Struct. -80 -60 -40 -20 0 20 40 60 80 y ( µ m)-250-200-150-100-50050100150200250 u ( V ) Figure 1.
Results for model ( ?? ). For meaning of line patternsrefer to captions of Figure ?? and of Figure ?? .5 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7-500505 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7-50050-40 -30 -20 -10 0 10 20 30 40-40-2002040 m ( V ) y ( µ m ) y ( µ m ) t (s)(a) t (s)(b) r ( µ m)(c)17.5 18 18.51819 5.2 5.24 5.28363840 5.2 5.24 5.28424446 Figure 2.
Hysteresis compensation for the piezoelectricactuator ( ?? ). (a) Compensation inputs, (b) temporal responsesand in (c) hysteresis loops. (- -) results obtained withcompensator ( ?? ) ( · · · ) results with compensator ( ?? ), (- · -) system output without compensation, and (—) displacementreference r = 40 sin(2 πt ) µ m. 0 5 10 15 20 T s (s)012345678 M A PE ( % ) M A PE ( % ) × − (a) T s (s) × − (b) Figure 3.
MAPE index ( ?? ) computed for the models andcompensators described, respectively, by equations (a) ( ?? ) and( ?? ); (b) ( ?? ) and ( ?? ). ( ◦ ) model and ( • ) tracking accuracies.( N ) accuracy of uncompensated system. β = 0 . β = 0 . β = 0 . · · · · · · Figure 4.
Bouc-Wen hysteresis loops within the investigatedrange.50 55 60 65 70 t (s)2.533.5 u ( V )
50 55 60 65 702.833.23.450 55 60 65 702.833.23.4 50 55 60 65 702.533.52.5 3 3.52.833.23.4 2.8 3 3.2 3.42.533.5 y ( V ) y ( V ) t (s) y ( V ) u ( V ) u ( V ) t (s) u (V) t (s) y (V)(a) (b)(c) (d)(e) (f) Figure 5.
Left column refers to model ( ?? ) and right columnto model ( ?? ). (a) input u k =0 .
56 sin(0 . πk ) + 3 V and (c) thecorresponding measured output (—) y and (- -) model ( ?? )free-run simulation; (b) smoothed version of y in (c); (d) thecorresponding output which is u k in (a) and (- -) model ( ?? )free-run simulation. (e) and (f) show the same data as (c) and(d), respectively. Figure 8.
MAPE index (32) computed for the models and compensatorsdescribed, respectively, by equations (a) (28) and (34); (b) (33) and (35).( ◦ ) model and ( • ) tracking accuracies. ( (cid:78) ) accuracy of uncompensated system. Finally, the same analysis was carried out for situations with different shapesof the hysteresis loop varying β in the range 0 . ≤ β ≤ . .
002 (see figure 9). The results are quite similar to the ones described so far andare not shown.
Smart Mater. Struct. -80 -60 -40 -20 0 20 40 60 80 y ( µ m)-250-200-150-100-50050100150200250 u ( V ) Figure 1.
Results for model ( ?? ). For meaning of line patternsrefer to captions of Figure ?? and of Figure ?? .5 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7-500505 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7-50050-40 -30 -20 -10 0 10 20 30 40-40-2002040 m ( V ) y ( µ m ) y ( µ m ) t (s)(a) t (s)(b) r ( µ m)(c)17.5 18 18.51819 5.2 5.24 5.28363840 5.2 5.24 5.28424446 Figure 2.
Hysteresis compensation for the piezoelectricactuator ( ?? ). (a) Compensation inputs, (b) temporal responsesand in (c) hysteresis loops. (- -) results obtained withcompensator ( ?? ) ( · · · ) results with compensator ( ?? ), (- · -) system output without compensation, and (—) displacementreference r = 40 sin(2 πt ) µ m. 0 5 10 15 20 T s (s)012345678 M A PE ( % ) M A PE ( % ) × − (a) T s (s) × − (b) Figure 3.
MAPE index ( ?? ) computed for the models andcompensators described, respectively, by equations (a) ( ?? ) and( ?? ); (b) ( ?? ) and ( ?? ). ( ◦ ) model and ( • ) tracking accuracies.( N ) accuracy of uncompensated system. β = 0 . β = 0 . β = 0 . · · · · · · Figure 4.
Bouc-Wen hysteresis loops within the investigatedrange.50 55 60 65 70 t (s)2.533.5 u ( V )
50 55 60 65 702.833.23.450 55 60 65 702.833.23.4 50 55 60 65 702.533.52.5 3 3.52.833.23.4 2.8 3 3.2 3.42.533.5 y ( V ) y ( V ) t (s) y ( V ) u ( V ) u ( V ) t (s) u (V) t (s) y (V)(a) (b)(c) (d)(e) (f) Figure 5.
Left column refers to model ( ?? ) and right columnto model ( ?? ). (a) input u k =0 .
56 sin(0 . πk ) + 3 V and (c) thecorresponding measured output (—) y and (- -) model ( ?? )free-run simulation; (b) smoothed version of y in (c); (d) thecorresponding output which is u k in (a) and (- -) model ( ?? )free-run simulation. (e) and (f) show the same data as (c) and(d), respectively. Figure 9.
Bouc-Wen hysteresis loops within the investigated range.
6. Experimental Results
Both identification and compensation strategies are now applied to an experimentalpneumatic control valve. This type of actuator is widely used in industrial processes,for which control performance can degrade significantly due to valve problems causedby nonlinearities [48] such as friction [49, 50], dead-zone, dead-band and hysteresis[2]. Hence, in this section we aim at compensation hysteresis using the developedtechniques.The measured output is the stem position of the pneumatic valve and the input isa signal that, after passing V/I and I/P conversion, becomes a pressure signal appliedto the valve. The sampling time is T s = 0 .
01 s. For model identification, the inputis set as white noise low-pass filtered at 0 . . N = 20000). Theidentification of the direct M and inverse ˘ M models was performed as in section 5. dentification and nonlinearity compensation of hysteresis (cid:96) = 3, n y = 1 and n u = 2. The modelparameters are estimated using (5) in order to comply with Lemma 1.The estimated model M is y k = y k − − . φ , k − + 19 . φ , k − +9 . φ , k − φ , k − u k − − . φ , k − φ , k − y k − , (37)and the inverse model ˘ M isˆ u k = ˆ u k − + 86 .
67 ˘ φ , k − − .
02 ˘ φ , k − − .
98 ˘ φ , k − y k − + 1 .
72 ˘ φ , k − ˘ φ , k − y k − − .
13 ˘ φ , k − ˘ φ , k − ˆ u k − , (38)which was estimated from a smoothed version of y k obtained by quadratic regression.This is done only to estimate ˘ M to avoid the error-in-the-variables problem, since y k serves as the input for ˘ M . Each model performance is given in figure 10 and Table 4. Smart Mater. Struct. -80 -60 -40 -20 0 20 40 60 80 y ( µ m)-250-200-150-100-50050100150200250 u ( V ) Figure 1.
Results for model ( ?? ). For meaning of line patternsrefer to captions of Figure ?? and of Figure ?? .5 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7-500505 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7-50050-40 -30 -20 -10 0 10 20 30 40-40-2002040 m ( V ) y ( µ m ) y ( µ m ) t (s)(a) t (s)(b) r ( µ m)(c)17.5 18 18.51819 5.2 5.24 5.28363840 5.2 5.24 5.28424446 Figure 2.
Hysteresis compensation for the piezoelectricactuator ( ?? ). (a) Compensation inputs, (b) temporal responsesand in (c) hysteresis loops. (- -) results obtained withcompensator ( ?? ) ( · · · ) results with compensator ( ?? ), (- · -) system output without compensation, and (—) displacementreference r = 40 sin(2 πt ) µ m. 0 5 10 15 20 T s (s)012345678 M A PE ( % ) M A PE ( % ) × − (a) T s (s) × − (b) Figure 3.
MAPE index ( ?? ) computed for the models andcompensators described, respectively, by equations (a) ( ?? ) and( ?? ); (b) ( ?? ) and ( ?? ). ( ◦ ) model and ( • ) tracking accuracies.( N ) accuracy of uncompensated system. β = 0 . β = 0 . β = 0 . · · · · · · Figure 4.
Bouc-Wen hysteresis loops within the investigatedrange.50 55 60 65 70 t (s)2.533.5 u ( V )
50 55 60 65 702.833.23.450 55 60 65 702.833.23.4 50 55 60 65 702.533.52.5 3 3.52.833.23.4 2.8 3 3.2 3.42.533.5 y ( V ) y ( V ) t (s) y ( V ) u ( V ) u ( V ) t (s) u (V) t (s) y (V)(a) (b)(c) (d)(e) (f) Figure 5.
Left column refers to model ( ?? ) and right columnto model ( ?? ). (a) input u k =0 .
56 sin(0 . πk ) + 3 V and (c) thecorresponding measured output (—) y and (- -) model ( ?? )free-run simulation; (b) smoothed version of y in (c); (d) thecorresponding output which is u k in (a) and (- -) model ( ?? )free-run simulation. (e) and (f) show the same data as (c) and(d), respectively. Figure 10.
Left column refers to model (37) and right column to model (38).(a) input u k =0 .
56 sin(0 . πk ) + 3 V and (c) the corresponding measured output(—) y and (- -) model (37) free-run simulation; (b) smoothed version of y in (c);(d) the corresponding output which is u k in (a) and (- -) model (38) free-runsimulation. (e) and (f) show the same data as (c) and (d), respectively. Table 4.
Performance of the modeling step. Experimental results.
Design Strategy Model
MAPESection 4.2 (37) 3 . . m k = 119 . (cid:104) r k +1 − r k +19 . m k − +19 . m k − − m k − ] − . m k − − m k − )[ m k − − m k − ] m k − +12 . m k − − m k − )[ m k − − m k − ] r k (cid:105) , (39) dentification and nonlinearity compensation of hysteresis m k = ˘ m k − +86 . r k +1 − r k ] − . r k − r k − ] − . r k +1 − r k ] r k + 1 . r k − r k − )[ r k − r k − ] r k − . r k − r k − )[ r k − r k − ] ˘ m k − . (40)The results of the experiment are shown in figure 11 and assessed in Table 5. Notethat both approaches significantly reduce the tracking error. Based on both numericaland experimental results, it seems that the performance of the compensators is directlyrelated to the accuracy of the identified model; see Table 4 and Table 5. Hence, asbefore, these results suggest that the compensation effort tends to be lower and moreeffective whenever the identified models are more accurate.The compensation produced by (40) is smoother than the one obtained with(39); see figure 11(a). This occurs because, for the compensator (40), the argumentof the sign function depends on the difference of the reference signal, while, for thecompensator (39), it depends on the difference of the autoregressive variable whichusually produces stronger oscillations and sudden changes; see figure 11(a), e.g. in therange of 51 −
53 s. As a result, larger compensation effort is required as quantified byNSAVI (36) in Table 5.
Smart Mater. Struct.
50 52 54 56 58 60 62 64 66 68 70 t (s)2.533.5 m ( V )
50 55 60 65 70 t (s)2.833.23.4 50 55 60 65 70 t (s)2.833.23.42.8 3 3.2 3.4 r (V)2.833.23.4 2.8 3 3.2 3.4 r (V)2.833.23.4 y ( V ) y ( V ) y ( V ) y ( V ) (a)(d)(b) (e)(c) Figure 1.
Hysteresis compensation for the pneumatic valve.(a) Compensation inputs, (b) and (c) its temporal responsesand in (d) and (e) the hysteresis loops. (- -) illustrates theresults obtained with compensator ( ?? ), ( · · · ) refers to the resultsby using compensator ( ?? ), (- · -) the system output withoutcompensation, and (—) the reference r =0 .
56 sin(0 . πt )+3 V. Figure 11.
Hysteresis compensation for the pneumatic valve. (a) Compensationinputs, (b) and (c) its temporal responses and in (d) and (e) the hysteresis loops.(- -) illustrates the results obtained with compensator (39), ( · · · ) refers to theresults by using compensator (40), (- · -) the system output without compensation,and (—) the reference r =0 .
56 sin(0 . πt )+3 V.
7. Conclusions
This work addressed the problems of identification and compensation of hystereticsystems. In the context of system identification , the contribution is twofold. First, dentification and nonlinearity compensation of hysteresis Table 5.
Performance of the compensation step. Experimental results.
Design Strategy Compensator
MAPE NSAVISection 4.2 (39) 5 .
514 1 . .
939 1 . .
602 1 . continuum of equilibrium points at steady-state, which is an importantingredient for hysteresis [5,30]. To this aim, a particular constraint on the parametersis presented in Lemma 1. As a consequence, the identified models are able to describeboth dynamical and static features of the hysteresis nonlinearity. Second, followinga quasi-static analysis of these models, a schematic framework is proposed to explainhow the hysteresis loop occurs on the input-output plane (see Figure 1).In the context of hysteresis compensation , this paper introduces two strategiesto design compensators. An important aspect of such procedures is that they showhow to restrict the pool of candidate regressors aiming at solving the compensationproblem. Such strategies are not limited to hysteresis and can be extended to othernonlinearities.The effectiveness of the compensation schemes is illustrated by means of numericaland experimental tests. For the strategy described in Section 4.2, the compensationlaw is obtained from the identified model by simple algebraic manipulations. Inthe case of the strategy introduced in Section 4.3, the compensators are identifieddirectly from the data. This, however, leads to models that are noncausal by nature.Guidelines to circumvent this problem are provided. The compensators designed byboth strategies can be readily employed in online compensation schemes.As a general remark, we observed in all our examples that the quality of theachieved compensation is correlated with the accuracy of the identified model (compareTable 2 with Table 3 and Table 4 with Table 5). Finally, we noticed that the identifiedmodels have a discontinuity due to the sign function used in some regressors. Whenthe model has many such terms, it sometimes happens that the compensation signalpresents abrupt transitions. The use of smoother functions in place of the signfunction, in order to alleviate this problem, will be investigated in the future. Acknowledgments
The authors would like to thank Arthur N Montanari for the insightful discussions.PEOGBA, BOST and LAA gratefully acknowledge financial support from CNPq(Grants Nos. 142194/2017-4, 310848/2017-2 and 302079/2011-4) and FAPEMIG(TEC-1217/98).
References [1] A Visintin.
Differential Models of Hysteresis . Springer, Berlin Heidelberg, 1994.[2] M A A S Choudhury, S L Shah, and N F Thornhill.
Diagnosis of Process Nonlinearities andValve Stiction: Data Driven Approaches . Springer, Berlin Heidelberg, 2008.[3] M Rakotondrabe.
Smart Materials-Based Actuators at the Micro/Nano-Scale: Characteriza-tion, Control and Applications . Springer, New York, 2013. dentification and nonlinearity compensation of hysteresis [4] J Peng and X Chen. A Survey of Modeling and Control of Piezoelectric Actuators. ModernMechanical Engineering , 3(1):1–20, 2013.[5] K A Morris. What is Hysteresis?
Applied Mechanics Reviews , 64(5):050801, 2011.[6] S A M Martins and L A Aguirre. Sufficient Conditions for Rate-Independent Hysteresis inAutoregressive Identified Models.
Mechanical Systems and Signal Processing , 75:607–617,2016.[7] G Tao and P V Kokotovic. Adaptive Control of Plants with Unknown Hystereses.
IEEETransactions on Automatic Control , 40(2):200–212, 1995.[8] C Visone. Hysteresis Modelling and Compensation for Smart Sensors and Actuators.
Journalof Physics: Conference Series , 138(1):012028, 2008.[9] H Chaoui and H Gualous. Adaptive Control of Piezoelectric Actuators with Hysteresisand Disturbance Compensation.
Journal of Control, Automation and Electrical Systems ,27(6):579–586, 2016.[10] S Yi, B Yang, and G Meng. Ill-conditioned dynamic hysteresis compensation for a low-frequencymagnetostrictive vibration shaker.
Nonlinear Dynamics , 96(1):535–551, 2019.[11] V Hassani, T Tjahjowidodo, and T N Do. A survey on hysteresis modeling, identification andcontrol.
Mechanical Systems and Signal Processing , 49(1-2):209–233, 2014.[12] Y K Wen. Method for Random Vibration of Hysteretic Systems.
Journal of the EngineeringMechanics Division , 102(2):249–263, 1976.[13] J Oh and D S Bernstein. Semilinear Duhem model for rate-independent and rate-dependenthysteresis.
IEEE Transactions on Automatic Control , 50(5):631–645, 2005.[14] P Ge and M Jouaneh. Tracking Control of a Piezoceramic Actuator.
IEEE Transactions onControl Systems Technology , 4(3):209–216, 1996.[15] M Brokate and J Sprekels.
Hysteresis and Phase Transitions . Springer-Verlag, New York, 1996.[16] A W Smyth, S F Masri, E B Kosmatopoulos, A G Chassiakos, and T K Caughey. Developmentof Adaptive Modeling Techniques for Non-linear Hysteretic Systems.
International Journalof Non-Linear Mechanics , 37(8):1435–1451, 2002.[17] G Quaranta, W Lacarbonara, and S F Masri. A review on computational intelligence foridentification of nonlinear dynamical systems.
Nonlinear Dynamics , 2020.[18] R W K Chan, J K K Yuen, E W M Lee, and M Arashpour. Application of Nonlinear-Autoregressive-Exogenous model to predict the hysteretic behaviour of passive controlsystems.
Engineering Structures , 85:1–10, 2015.[19] H V H Ayala, D Habineza, M Rakotondrabe, C E Klein, and L S Coelho. Nonlinearblack-box system identification through neural networks of a hysteretic piezoelectric roboticmicromanipulator.
IFAC-PapersOnLine , 48(28):409–414, 2015.[20] J Fu, G Liao, M Yu, P Li, and J Lai. NARX neural network modeling and robustness analysisof magnetorheological elastomer isolator.
Smart Materials and Structures , 25(12):125019,2016.[21] L A Aguirre. A Bird‘s Eye View of Nonlinear System Identification. arXiv:1907.06803 [eess.SY] ,2019.[22] I J Leontaritis and S A Billings. Input-Output Parametric Models for Non-Linear Systems PartI: Deterministic Non-Linear Systems.
International Journal of Control , 41(2):303–328, 1985.[23] I J Leontaritis and S A Billings. Input-Output Parametric Models for Non-Linear Systems PartII: Stochastic Non-Linear Systems.
International Journal of Control , 41(2):329–344, 1985.[24] R K Pearson.
Discrete-Time Dynamic Models . Oxford University Press, Oxford, 1999.[25] A Leva and L Piroddi. NARX-based Technique for the Modelling of Magneto-RheologicalDamping Devices.
Smart Materials and Structures , 11(1):79–88, 2002.[26] L Deng and Y Tan. Modeling Hysteresis in Piezoelectric Actuators Using NARMAX Models.
Sensors and Actuators A: Physical , 149(1):106–112, 2009.[27] K Worden and R J Barthorpe. Identification of hysteretic systems using NARX models, Part I:evolutionary identification. In
Topics in Model Validation and Uncertainty Quantification ,volume 4, pages 49–56. Springer, 2012.[28] R Dong and Y Tan. Inverse Hysteresis Modeling and Nonlinear Compensation of Ionic PolymerMetal Composite Sensors. In
Proceeding of the 11th World Congress on Intelligent Controland Automation , pages 2121–2125, Shenyang, China, 2014.[29] W R Lacerda J´unior, S A M Martins, E G Nepomuceno, and M J Lacerda. Control of HystereticSystems Through an Analytical Inverse Compensation based on a NARX model.
IEEEAccess , pages 1–1, 2019.[30] D S Bernstein. Ivory Ghost [Ask The Experts].
IEEE Control Systems Magazine , 27(5):16–17,2007.[31] S A Billings and S Chen. Extended Model Set, Global Data and Threshold Model Identification dentification and nonlinearity compensation of hysteresis of Severely Non-Linear Systems. International Journal of Control , 50(5):1897–1923, 1989.[32] N R Draper and H Smith.
Applied regression analysis . John Wiley & Sons, New York, 3 edition,1998.[33] S Chen, S A Billings, and W Luo. Orthogonal Least Squares Methods and their Application toNon-Linear System Identification.
International Journal of Control , 50(5):1873–1896, 1989.[34] H Akaike. A New Look at the Statistical Model Identification.
IEEE Transactions on AutomaticControl , 19(6):716–723, 1974.[35] L Piroddi. Simulation Error Minimisation Methods for NARX Model Identification.
International Journal of Modelling, Identification and Control , 3(4):392–403, 2008.[36] S A M Martins, E G Nepomuceno, and M F S Barroso. Improved Structure Detection ForPolynomial NARX Models Using a Multiobjective Error Reduction Ratio.
Journal of Control,Automation and Electrical Systems , 24(6):764–772, 2013.[37] A Falsone, L Piroddi, and M Prandini. A randomized algorithm for nonlinear model structureselection.
Automatica , 60:227–238, 2015.[38] P F L Retes and L A Aguirre. NARMAX Model Identification Using a Randomised Approach.
International Journal of Modelling, Identification and Control , 31(3):205–216, 2019.[39] I B Q Ara´ujo, J P F Guimar˜aes, A I R Fontes, L L S Linhares, A M Martins, and F M U Ara´ujo.NARX Model Identification Using Correntropy Criterion in the Presence of Non-GaussianNoise.
Journal of Control, Automation and Electrical Systems , 30(4):453–464, 2019.[40] F Ikhouane and J Rodellar.
Systems with Hysteresis: Analysis, Identification and Control Usingthe Bouc-Wen Model . John Wiley & Sons, 2007.[41] L A Aguirre and E M A M Mendes. Global Nonlinear Polynomial Models: Structure, TermClusters and Fixed Points.
International Journal of Bifurcation and Chaos , 6(2):279–294,1996.[42] L A Aguirre. Identification of smooth nonlinear dynamical systems with non-smooth steady-state features.
Automatica , 50(4):1160–1166, 2014.[43] L A Aguirre, M F S Barroso, R R Saldanha, and E M A M Mendes. Imposing steady-stateperformance on identified nonlinear polynomial models by means of constrained parameterestimation.
IEE Proceedings - Control Theory and Applications , 151(2):174–179, 2004.[44] P Q Xia. An Inverse Model of MR Damper using Optimal Neural Network and SystemIdentification.
Journal of Sound and Vibration , 266(5):1009–1023, 2003.[45] M Rakotondrabe. Bouc-Wen Modeling and Inverse Multiplicative Structure to CompensateHysteresis Nonlinearity in Piezoelectric Actuators.
IEEE Transactions on AutomationScience and Engineering , 8(2):428–431, 2011.[46] G Y Gu, M J Yang, and L M Zhu. Real-Time Inverse Hysteresis Compensation of PiezoelectricActuators with a Modified Prandtl-Ishlinskii Model.
Review of Scientific Instruments ,83(6):065106, 2012.[47] W R Lacerda J´unior, S A M Martins, and E G Nepomuceno. Influence of Sampling Rate andDiscretization Methods in the Parameter Identification of Systems with Hysteresis.
Journalof Applied Nonlinear Dynamics , 6(4):509–520, 2017.[48] R Srinivasan and R Rengaswamy. Stiction Compensation in Process Control Loops: AFramework for Integrating Stiction Measure and Compensation.
Industrial & EngineeringChemistry Research , 44(24):9164–9174, 2005.[49] R A Romano and C Garcia. Valve friction and nonlinear process model closed-loop identification.
Journal of Process Control , 21(4):667–677, 2011.[50] J R Baeza and C Garcia. Friction compensation in pneumatic control valves through feedbacklinearization.