WIN“lWlWllllWlWNWWIWIWINNil

     

   

14
H

-_1m
.5
'(DNO

This is to certify that the
thesis entitled

AN INVESTIGATION OF MODEL REFERENCE ADAPTIVE
CONTROL OF UNKNOWN DYNAMIC HYSTERETIC
SYSTEMS USING SLOW ADAPTATION

 

 

presented by

James J Reynolds

LIBRARY
Michigan State
University

 

has been accepted towards fulfillment

 

of the requirements for the

Master of degree in Electrical and Computer
Science EngineerinL

 

 

Major ProfessoW

DE C e meP I” 752007

Date

MSU is an affinnative-action, equal-opportunity employer

.n-.-.- _..— - ~n---a----n-u-u-.-.-.-._ -<.. _ --

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/07 p:lClRC/DateDue.indd-p.1

 

 

 

AN INVESTIGATION OF MODEL REFERENCE ADAPTIVE CONTROL OF
UNKNOWN DYNAMIC HYSTERETIC SYSTEMS USING SLOW ADAPTATION

’By

James J Reynolds

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE
Department of Electrical and Computer Engineering
2007

ABSTRACT

AN INVESTIGATION OF MODEL REFERENCE ADAPTIVE CONTROL OF
UNKNOWN DYNAMIC HYSTERETIC SYSTEMS USING SLOW ADAPTATION

By

James J Reynolds

This thesis is a first step in analyzing model reference adaptive control of an unknown
linear system preceded by an unknown hysteresis nonlinearity. The parameters of the
hysteresis nonlinearity are adapted much slower than the parameters of the linear system,
which are in turn adapted much slower than the plant dynamics. This approach allows the
system to be considered under an averaging framework, where slower variables are viewed
as constant with respect to faster variables. The problem statement is established and the
first steps in analyzing the system are laid out. Challenges that arise in this problem beyond
a classical averaging or singular perturbation problem are outlined. Simulation is used to
not only show the effectiveness of the control method, but also to gain some additional

insight into properties and behaviors of the closed-loop system.

For my Mother, Earaina Lynn, and my biggest supporter, Jennifer Lynn.

iii

ACKNOWLEDGEMENTS

I would like to express my sincerest thanks to those who made this work possible. With-
out their guidance, support and unwavering patience, I could never have known where or
how to begin this project, or how far I could stretch my own abilities. I would like to
specifically thank Dr. Hassan Khalil and Dr. Xiaobo Tan for all of their time, effort and
invaluable advice. I would also like to thank Dr. Ning Xi for serving on the committee for

this thesis.

iv

TABLE OF CONTENTS

LIST OF TABLES ............................................. vi
LIST OF FIGURES ............................................ vii
1. Introduction ............................................... 1
2. Hysteresis Operator Review .................................... 4
2.1 The Preisach Operator .................................... 4
2.2 Inversion of the Preisach Operator ............................ 10
3. Problem Formulation ........................................ 12
3.1 System Error Dynamics .................................... 12
3.2 Uniqueness of the Solution .................................. 17
3.3 Adaptation Rules ........................................ 19
4. Averaging Analysis .......................................... 22
4.1 An Identification Problem .................................. 22
4.2 The Local Approximate Closed-Loop System ...................... 23
4.3 Local Lyapunov Stability ................................... 29
5. Simulation Results ........................................... 36
5.1 Parameter Convergence .................................... 37
5.2 Region of Attraction ...................................... 40
6. Future Work ............................................... 43
References .................................................. 44

5.1

LIST OF TABLES

Parameter convergence simulation settings.

vi

1.1
2.1
2.2
2.3
2.4
2.5
2.6
3.1
5.1
5.2
5.3
5.4

LIST OF FIGURES

A hysteresis operator preceding linear dynamics ................. 1
Input-output relationship of a hysteron operator. ................ 5
The Preisach plane. .............................. 5
The discretized Preisach plane, with discretization level L = 4. ......... 7
Schematic representation of a Preisach operator. ................ 8
Evolution of the memory curve M with time ................... 9
A Preisach operator and its approximate right inverse. ............. 10
The proposed controller with parameter estimates. ............... 13
Parameter error for case L = 4. ........................ 39
Parameter error for case L = 8. ........................ 40
Output error for case L = 8. .......................... 40
Comparison of regions of attraction ....................... 42

vii

1. Introduction

Model reference adaptive control (MRAC) of systems with hysteresis nonlinearities
has recently become a topic of interest, in particular regarding control of smart materials
[1, 2, 3]. A hysteresis operator preceding linear dynamics will be used to model these
systems. Refer to Fig. 1.1.

There are many different hysteresis models. In [4, 5], Tao and Kokotovic consider
a piecewise linear hysteresis model. This model has a small number of parameters that
describe it. It will lead to computational simplicity at the expense of model accuracy.
Overparameterization by way of a Kronecker product is used to accommodate the bilinear
coupling of parameters. Because of the inaccuracy in this model, Tan and Khalil instead
used a Preisach hysteresis model [6]. The same bilinear coupling of parameters is present
with a Preisach model as the piecewise linear model. For n” and up parameters associated
with the hysteresis operator and linear dynamics, respectively, overparameterization will
give anp parameters. The number of parameters of a Preisach operator can be quite large,
meaning the use of overparameterization will require a cumbersome amount of computa-
tion. For example, in [7], Tan and Baras used a 10-level discretization of the Preisach
plane which results in 56 hysteresis parameters. For a linear plant with denominator of
degree three, overparameterization will give 112 parameters from the Kronecker product.
To alleviate this difficulty, Tan and Khalil instead proposed using slow adaptation to create
a separation of time scales. In this way, the number of parameters remains manageable; 58
parameters in the above example. Slow adaptation also allows the behavior of the plant pa-

rameter adaptation and hysteresis parameter adaptation, as well as the closed-loop stability,

Hysteresis rator

    

 

 

 

Linear Dynamics

Figure 1.1: A hysteresis operator preceding linear dynamics.

to be evaluated using averaging and singular perturbation techniques.

For both hysteresis models (piecewise linear and Preisach) the output is shown to be
the inner product of a regressor vector and a vector of available signals (and possibly an
additive error term that results from inaccuracy in modelling). For this reason, the analysis
of an MRAC with slowly adapted parameters should apply to either model (and any other
hysteresis model whose output has the same structure), though slow adaptation will only
be beneficial when the number of hysteresis parameters is large. The focus of this thesis
will not be on the choice of a hysteresis model, but the performance and robustness of an
MRAC with slow adaptation when the hysteresis model satisfies certain assumptions.

There are two major goals in this thesis: to begin the analysis of an MRAC with
slowly adapted parameters with current averaging techniques and identify the obstacles
in completing the analysis, and to perform simulations that will provide useful insight
into the behavior of the completenonlinear system. The majority of the analysis herein
will rely heavily on “linearization” of the closed-loop system, and will only provide local
information regarding its behavior. Therefore, simulations will be used as motivation for
continuing this work and observing the stability requirements, convergence rates, basins
of attraction, etc.

The thesis is structured as follows. Sections 2.1 and 2.2 derive the Preisach oper-
ator and introduce useful relationships between the operator and its inverse. Section 3.1
formulates the closed-loop system equations, using a similar method to that of Sastry and
Bodson [8]. The uniqueness of the solution is established in Section 3.2. The parameter
adaptation rules are derived in Section 3.3, using a simple gradient law. Section 4.1 reit-
erates an identification problem from [6] to motivate Section 4.2, where the closed-loop
system is analyzed. This is the first known work establishing the closed-loop stability of
MRAC systems with inverse hysteresis compensation and slow parameter adaptation. The
theorems for stability rely on “linearization” of the closed-loop system about the origin,

averaging techniques in [8] and singular perturbation techniques from [9]. Linearization is

used loosely here, because the system is not described by a pure ordinary differential equa-
tion (ODE). Section 4.3 provided local Lyapunov analysis based on a converse Lyapunov
function argument. The Lyapunov stability of the system implies that the approximations
from Section 3.3 are appropriate. Sections 5.1 and 5.2 give simulation results illustrating
stability and parameter convergence and insight on the region of attraction, respectively.

Future work and research goals are highlighted in Section 6.

2. Hysteresis Operator Review

The particular hysteresis operator used in this thesis, both for theoretical and sim-
ulation purposes is the Preisach Operator. For a more detailed discussion of the Preisach
operator one can refer to [10]. In [7] Tan and Baras show that with sufficiently large
numbers of parameters a Preisach operator accurately describes the hysteresis in a magne-
torestrictive actuator, and that continually increasing the number of parameters has limited
benefits. This indicates that for a particular actuator, the dimensionality of the problem can
be changed to optimize the performance, with a tradeoff between computational complex-
ity and model accuracy. Because of this versatility of the Preisach operator it is preferable
to a piecewise continuous Operator as in [4, 5], which can only accurately describe very
simple hysteretic phenomena.

The method for inverting the Preisach operator used in simulations in this thesis is
described in [3]. The algorithm therein results in the exact inverse of a given Preisach
operator after several iterations. The iterative nature of the inversion method requires a
discrete-time implementation, with system sampling period greater than the necessary in-

version time.

2.1 The Preisach Operator

The basic building block of a Preisach operator is a hysteron. A hysteron is a delayed
relay, (pm, [-, -], with output either +1 or --1, where a Z [3 are the switching thresholds.
Fig. 2.1 shows the input-output curve of a generic hysteron. The arrows on the figure
show the direction of the delay behavior. For v E C ([0, T]) and an initial configuration

(po 6 {-1, l}, (o = (p5,a[v,(po] is defined as

+1 a < v(t),
w(t) 3 _1 fl > v0), fort e [0,T] (2.1)

to (t‘) otherwise,

+1

 

A

 

 

 

 

V

 

Figure 2.1: Input-output relationship of a hysteron operator.

Preisach plane -

 

 

Figure 2.2: The Preisach plane.

where (0(0‘) 3 (pa and

- A .
t 2 11m t—E,
e>0,e-+O

as in [3]. From (2.1) it is clear that the current output of the hysteron operator depends on
the output history. The hysteron operator has memory, as is typical of natural hysteretic

phenomena. Define the Preisach plane as
92{(p,a)enznga}. (2.2)

The Preisach plane is clearly one half of R2, as shown in Fig. 2.2. In addition to the

Preisach plane, define the weighting function
#010!) 20V(fi,a) E «9”. (2.3)

The weighting function scales the portion of the Preisach operator output associated with

a ([3, a) pair. The Preisach operator, 1", has output

um = /§H(B,a)¢p,a lv, mm» (t)dBda- (2.4)

In the most general case, the Preisach operator depends on an infinite collection of
hysterons and their respective weights, making it infeasible to implement. Therefore, a
truncated region of the Preisach plane will be used and divided into K cells. The dis-
cretization level, L, is the number of cells in either the a or B direction. This means the

number of resulting cells is K = L(L+ 1) / 2. Define the truncated Preisach plane as
yré{(p,a)en2:pgaandpmgpandagam}.

Truncating the Preisach plane will limit the minimum and maximum output levels of the
Preisach operator, which is acceptable as many real actuators will have a limited range of
motion. Fig. 2.3 shows an example of a truncated Preisach plane with L = 4 and K = 10.
\Vrthin each cell of the truncated Preisach plane, the weighting function is assumed to be
uniform, i.e.

“(13,00 = a,“ v ([3, a) 6 km cell. (2.5)

The signed area of the k’h cell, wHJc, is defined as the lower area, C; , minus the upper

area, Ck- . These areas are defined as

C? = {(p,a) 6 km cell: (pp,a[v,(po](t) = i1}. (2.6)

 

 

 

 

 

 

 

 

Figure 2.3: The discretized Preisach plane, with discretization level L = 4.

The total area of any cell is unity, i.e. if q)”, [v, {pg} (t) = +1 (—1)V(fi,a) 6 km cell, then

with (t) = +1 (— l ). The output of the truncated Preisach Operator can now be written as

K
“(0 = 2 WH,k9H,k + 9H,o, (2.7)
k=l

where 03,0 is the contribution from hysterons outside of the discretized Preisach plane. If
the truncated Preisach plane underestimates the range of the actuator, there will be a set of

hysterons that can never change states. From this point on, (2.7) will be expressed as the

inner product
u(t) .—_— egwfln), (2.8)
where _ _ _ _
1 911,0
WH,1(‘) 911,1
wH(t) = . , 0” = . . (2.9)
.mm . t 9m.

 

 

 

 

Schematically, from this point on, the Preisach Operator will be represented as a block
labeled 1" as in Fig. 2.4.

The regressor in (2.9) is known for all time I, provided the initial memory curve is

 

v( t) “(1)

 

 

 

Figure 2.4: Schematic representation of a Preisach operator.

known. The memory curve indicates the dividing line between all hysterons with output
+1 and with output — 1. Given the initial memory curve, it is straightforward to obtain the
memory curve at time t. Fig. 2.5 illustrates the general idea The initial memory curve
is labeled M (t = O) in Fig. 2.5 (a). The input to the hysteron operator, v, at time t = O is
labeled v0. The input will always be the point where the memory curve intersects the line
[3 = a. The memory curve divides the Preisach plane into two regions. Every hysteron
whose (B, a) pair lies above or to the right (below or to the left) of the memory curve will
have configuration —1 (+1). If the input to the Preisach operator is then monotonically
decreased to V] at time t = 1, then the corresponding memory curve is M (t = 1) in Fig. 2.5
(b). Decreasing the input to any hysteron can only change the output of that hysteron if
the fl-level is crossed. This gives the vertical portion of the memory curve in M (t = 1);
all of the hysterons to the right of the memory curve that were previously (in Fig. 2.5 (a))
configured as +1 have switched to —l with their respective B-level crossings. Figs. 2.5 (c)
and ((1) show the memory curve when the output is increased monotonically to V3 at time
t = 3, then increased monotonically to V4 at time t = 4. Since V4 > v0, the last increase
removes all the comers from the memory curve when the output surpasses its uppermost
portion. Combining any initial memory curve with any discretization, the lower and upper
areas, Cf and C; , of each cell can be obtained by updating the memory curve as previously
illustrated.

One more important property of the Preisach Operators used in this thesis is that
they are rate-independent. This means that the output of each hysteron and therefore the
Preisach Operator itself does not depend on the slope of the input. The parameter vector

0” is invariant with regard to the input. In [2, 3], Tan and Baras show that for a specific

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

M(t = O) / _ M(t = 1)
v0 i
I? .8
.v,
(a) (b)

a a /
V4

M(t = 3) M(t = 4)

V3

19 fl

 

 

 

 

(C) ((0

Figure 2.5: Evolution of the memory curve M with time.

actuator, rate-independence is a valid assumption for low input frequencies. They derive
a dynamical model to approximate the rate-dependence at higher frequencies. In [6], the
linear dynamics in cascade with the Preisach Operator are shown to adequately capture
the rate-dependent behavior of a given piezoelectric positioning scheme, which illustrates
the motivation for both the Preisach operator (memory) and the linear dynamics (rate-
dependence).

TO summarize the previous section, the truncated Preisach operator is a way of ap-
proximating a hysteretic phenomenon. The output of a truncated Preisach Operator can be

expressed as the inner product of a time-varying regressor vector and a constant parameter

 

 

ud(t) f_l v(t)

 

 

 

r u(t) >

 

 

 

 

 

 

Figure 2.6: A Preisach operator and its approximate right inverse.

vector. The regressor vector is known if both the initial memory curve and the input of the
Preisach operator are known. The memory curve can always be set to a known “initial”

value by driving the input to the Preisach operator above am,“ or below Bmin.

2.2 Inversion of the Preisach Operator

The exact inverse of a discrete Preisach operator is a topic Of the work in [3]. The
algorithm for (exactly) inverting the Preisach operator is not of concern in this thesis, so
additional information can be found therein. Certain relationships of a Preisach operator
and its inverse are essential to the development in later sections.

The inversion algorithm in [3] will produce the exact “right” inverse of the Preisach

operator (see Fig. 2.6) when the true parameter vector, OH is known. This gives

“(1) = (For—1) (WW) = “(1(1)-

From [6], if instead of the exact parameters, the inversion is based on a parameter estimate,

9H (t), then the inversion will result in
x 1' ~7-
ud(t) — u(t) = (0H(t) — 0H) WH(t) = GHWHO). (2.10)

The Preisach inverse is a user-defined quantity with discretization level L. When used to
(approximately) invert a Preisach Operator with a continuous weighting function, rather

than discrete, an error term is introduced to (2.10). The new relationship is
ud(t) — u(t) = égwuo) + dH(t), (2.11)

10

where d” (t) is a bounded disturbance with bound established in [3]. The last equation is
the same relationship as that for the piecewise linear hysteresis model in [5]. This means
that the overparameterization method proposed by Tao and Kokotovic could be applied us-
ing a Preisach operator rather than the piecewise linear hysteresis model, if computational

complexity were not an issue.

11

3. Problem Formulation

3.1 System Error Dynamics

The control scheme in this thesis uses multi-time scale slow adaptation to accommo-
date bilinearly coupled parameters without the use of overparameterization. The MRAC
framework will be the same as in [5]. This work was first undertaken by Tan and Khalil
[6] and continued by Reynolds, Tan and Khalil [11].

Let Gp (s) = kpgi—gg, where kp is the high-frequency gain, Pp(s) and Zp (s) are monic,
coprime polynomials Of degree n and m, respectively. The plant has state space represen-

tation (Ap,Bp,Cp,O) so that

xpm = Apxpm +Bp“(‘)

y(t) = Cpxp(t). (3.1)

The goal of the controller design is to make the plant output y(t) track the reference output.
The output of the model is given by ym (t) = G,,. (s) [r] (t), where r(t) is a bounded, piece-
wise continuous reference input and Gm (s) [r] (t) is the time domain output of the transfer
function Gm(s) acting on the signal r(t). The following assumptions, from [5, 6, 11], are

made about the plant and reference model:
— (Al) Zp (s) is a stable polynomial;
— (A2) The degrees n and m are known;
— (A3) kp > O;

— (A4) The model, G,,,(s) = m1}? has stable polynomial Pm(s) of degree n* = n — m.
The proposed controller structure (shown in Fig. 3.1) consists of two parts: a clas-
sical MRAC that can be found in adaptive control texts such as [8, 12], and the Preisach
inverse estimate, which will directly precede the actuator illustrated in Fig. 1.1.

With exact hysteresis cancellation and knowledge of the plant parameters, perfect

12

 

 

 

 

 

 

 

 

 

12,9;

 

 

 

 

 

* 6? lww) 1(5)
AS)

A __1 v(t) u(t) ip=Apxp+BPu y(t) A

E F y=Cpxp '

g\ “d (t)
w2(t) 1(3).
Ell—2o
a

Figure 3.1: The proposed controller with parameter estimates.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

model following can be achieved with

u(t) = 91TW1(‘)+92TW2(‘)+920)’(’)+93’(1), (3.2)
where

A a(s) A a(s) a(s) _ 1 S

W1(I)=m[u](t), w2(t):1(—S)U’I(t)a RES—m r

 

 

and 2(s) is a stable polynomial of degree n — l. The parameters 91 E R""1, 92 E R""1,

920 E R and 93 E R are determined by the matching equation:
911 a(s)Pp(s) + [9{a(s) + 9201(3)] kpr(s) = 1(s) [Pp(s) — 93kpr(s)Pm(s)]. (3.3)
Now define a controllable canonical pair (A, 8;) such that

(sI—A)‘IB)_ = 9E. (3.4)

MS)

13

Using (3.1) and (3.4), the state space form of the system with perfect matching is

Wml

wm2

The closed-loop form of the perfect model following equations are

 

 

Wm1 = B]. OZOCp

 

3,9,7 3,9;
A+BA 9,7“ 3,1927

0 A

This is a non-minimal realization of the model, Gm(s), so

A
A -

 

3,9,? 3,9;
A+BA l9,T 3,19;

0 A

 

are defined with appropriate dimensions, giving

Ym

WT

A
_xT m

wherexm —[ pm, l,w;,2]T.

The signal w] (t) is unavailable because the output of the Preisach operator, u(t), is

unavailable. Since w1(t) and the true plant parameters are both unknown, the implemented

Amxm + er

mem a

14

 

 

I ..
Wm]
Wm2
d L.

||l>

 

 

Apxpm +Bp(01TWml ‘1' 92Tme + OZOCp-xpm + 93")

AWml +Bl(91TWml + QZTWmZ + 620Cpxpm + 93”)

 

3,93
31.93

Bp93
31.93

 

 

control is

ud(t) = éfwn(r)+észz(t)+ézoy(t)+63r(t). (3.5)
where 9,, is an estimate of 9;, and

w“ (t) 2 g—gmm. (3.6)

The true hysteresis parameters are also unknown, so the inversion will be performed based
on estimates of those parameters. The following two assumptions are made about the

hysteresis Operator, 1", and its approximate inverse, 1L]:

— (A5) F is a Preisach operator with piecewise uniform weighting function characterized
by OH 6 RM , as discussed in Sec. 2.2. Furthermore, the weights corresponding to each

cell adjacent to the line a = [3 are nonzero;
- (A6) The inversion error ud — u satisfies (2.11);

— (A7) The control effort ud remains (or is forced to remain) within the saturation limits
corresponding to the output range of the Preisach operator. ‘
The last assumption ensures that the finite-dimensional signal why can accurately model
the internal state Of the Preisach operator. If the control is unrestricted, than an infinite-
dimensional wH would be necessary to capture its effect, which is computationally infea-
sible and unnecessary for many true hysteretic systems. With the control signal ud in (3.5)

and (2.11), the system has the state space representation

Jip = Apxp+Bpu=Apxp+Bp(ud—é§wH—dy)

= Apxp+Bp(élTwn+92Twz+ézonxp+égr—é;wH-dy)
wl = Aw1+Blu = Aw1+B,1(ud — 5,7,.qu —dH)

= Aw1+B)(91Tw11+ észz + ézonxp + é3r— 95w” ——dH)

wz = sz +BACpxp.

15

Now, using (2.11), substituting 9). = 9;. + 9k and adding and subtracting 917‘ w, the closed-

loop state model becomes
JICp xp
W1 3 Am W] (3-7)
W2 W2

 

 

 

 

~

Bm "" ~ A A ~
+ B: (917w, + 627w2+620Cpxp+ 93r+ 91T(W11 —W1)"' 915%! ‘dH)°

The error state is defined as

 

 

 

 

._ 1 .. -
A
el = 812 2" W1 - Wm] a
813 W2 _Wm2

and subtracting (3.5) from (3.7), the error dynamics become
B ~ ~ ~ ~ .. ..
e] = Ame1+ Em (9.7% + 93m + Ozonxp + 93r+ 91T(w11 — w1)-— 915w” —dH). (3.8)

The driving terms on the right-hand side of (3.8) are dependent on the error state, so each

much be separated into internal and exogeneous parts,

 

 

 

 

 

 

r r r
Cpxp _ Cp (x17 _ xpm +1512!!!) _ Cp (e11 +xpm)
W1 W1 - Wml +Wm1 €12 +Wm1
W2 w2 — Wm2 + Wm2 d _ e13 + Wm2 _

The closed-loop dynamics now become

B... .. B... .. .. . ..
é1= [Am+-6:(9£Q)] el+-0-3-(93r+95Wm+91T(W11—Wl)—9liWH-dH), (3-9)

16

where

50-9; 61 ,wm(r)é wm1(t) .92 o 1 o . (3.10)
éZ Win20) 0 0 I

 

 

 

 

 

 

Cp 6 R1 x", I 6 R"“1x"'1 is the identity matrix and the remaining dimensions are defined
appropriately. Notice that the quantity W11 — W1 is left alone, as it depends on the hysteresis

parameter error, not the error state e1. In particular,

WNW—M(t): 781 udl()- Z—is—(lu l()= fiiiléuwwdnll) (3.11)

From this point forward we make the following assumption, so that conclusions on param-

eter convergence can be drawn:

— (A8) The bounded disturbance, d” (t), is identically zero, i.e. dH(t) = 0 V t.

3.2 Uniqueness of the Solution

It is shown by Sastry and Bodson in [8] that the matching equation (3.3) has a unique
set of parameters {91, 02, 020, 93} by which it is satisfied. For this case, with a Preisach
operator and its right inverse preceding the linear dynamics Gp (s), it will now be shown
that there are in fact an infinite number of solutions that produce the same input to the
linear dynamics, and therefore produce the same idea] output y... (t).

Consider the case when the exact Preisach inverse is known, and the problem reduces
to a classical MRAC. The control signal that gives asymptotic tracking of the reference

model then is given in (3.2) as
u(t) = 91TW1(I)+ 92Twz(t) + 920W) + 93r(t)-

This is the case when all parameter errors (plant and hysteresis) are zero. Now instead

17

consider the case

 

 

Q )
2’89’59 .33’9’

 

093

a6” a

for some constant a gé 0. Note that

am) — u(t) = ézwelr) = (a -1)9§wn(t) = (a —1)u(t) => um = au(t).

 

 

 

 

(a-1)02
(0‘1)920
(a—l)93

_ (a—1)9H _

 

(3.12)

Using (2.10) and substituting the values from (3.12) the control signal applied to the plant

is

au(t) =

ud(t) — éngn)

9100)

9110‘")

9W)
9111(8)

[u l(t) + 92 mm + ézoy(t) + ésr(t)- 9mm
[on] (t) +a627wz(t) +a920y(t) +a93r(t)- (a — 1)u(t)

a(91 W10) + 93mm + 020340 + 93r(t)) ,

(3.13)

and eliminating the factor a from both sides of (3.13) gives the same control signal as in

(3.2). Notice that although there are an infinite number of solutions, the value of 91 is the

same in all of them. In order to overcome the issue of ambiguity of multiple equivalent

solutions, and observe convergence Of the parameter error to zero, the reference signal gain

93 will be set to a fixed value. This will ensure that there is again a unique set Of parameter

estimates that produces asymptotic tracking. To simplify the notation, it will be assumed

that 93 = 93 = 21;, as derived in [8]. This way all of the true plant parameter values do not

change from those found by the matching equation (3.3) and the true hysteresis parameter

18

values 9” reflect the true gain of the Preisach operator. Setting the parameter estimate 93

to its true value is merely a convention; if instead 93 were some scaled version of 93 (e.g.

93 E 1), the aforementioned parameter values would be scaled accordingly as in (3.12) (by
_ 1

a factor Of a — 35-).

Having established S3 E 0, (3.9) can now be simplified to

B ~ B - A -
é1= [Am '1' Em(9gQ)] 81+ ‘6':- (95w... + 917‘ (W11 — W1) — 915W”) . (3.14)
3.3 Adaptation Rules
With the error dynamics defined, the parameter estimate dynamics will now be
derived. The adaptation rule is chosen using a gradient update law, which minimizes

%(Cm€1)2 = %(y — ym)2. The adaptation rule is

r o
G c e M (3.15)

_ ml
0 PH 30

Q)-
II

where

Q)
|Il>
£0) O >

and where PG 6 R2"‘1"2”‘1, 1"” E R’Wxn” are diagonal matrices with positive coeffi—
cients, 70 and m = 751/, respectively. Note that 95 = 90 + 90 and 9” = 9” + 9” are the
estimates of the plant and hysteresis parameters, respectively. Applying the adaptation rule

(3.15) to the error dynamics (3.14) gives the following adaptation law:

 

 

 

 

( - Cpxp - \
6 = -70Cmele W“ , (3-16)
W2
I ( ((3%)th _...,) l

19

where Gm([z(-)]) is the output of G... (s) acting component-wise on the input vector z(-),
0 < 70 << 1 separates the error dynamics and adaptation into two time scales, and O <
Y << 1 further separates the adaptation of the plant and hysteresis parameters into two
time scales. Notice that applying the gradient update law (3.15) to (3.8) gives a factor of
é. By (A3) this can be absorbed into the adaptation gain. All of the signals in (3.16) are
available online. The true parameters 9 are taken as constant, so 9 = 9. As in the case of

error dynamics, the parameter update law will be expressed in terms of error e1 and driving

 

 

 

 

signals.
( ' o ' )
,e Wm Q31 “’11 — W1
0 = —ngmele + + (3.17)
0 O 0
T x
K (Germ) 9......) l

Equations (3. 14) and (3.17) now give a complete expression for the closed-loop dynamics.
From this point on, W” (t) will be considered the approximation of w” (t) when the states
are near the origin, i.e. 3711(1) = wH(t) V t 2 0 if e1 (t) = O, 9(t) = O, V t 2 0. Now make

the following assumption:

— (A9) The signal W” (t) is bounded and continuous.
This assumption is not stringent, provided the input to the Preisach operator is bounded
and continuous. Using the fact that (ef, 5T) = O is an equilibrium point of (3.14) and
(3.17), a local approximation of the dynamics near the origin will be derived, so that local
results can be obtained. Establishing the stability of the approximate system is useful for
establishing the stability of the origin of the complete nonlinear system, at least in some
neighborhood of the true parameter values. The approximate dynamics in the neighbor-

hood of (e{, 9T) = 0 are found by linearizing (3.14) and (3.17) after substituting W” (t)

20

for wH(t). The results are

B ~ w

93 (aerial—Wu

Y ((fiillwggr 91 —WH)

Since the right-hand side now only depends linearly on el and 9, a new signal is defined

Qt.
||

_ I’GCme l Gm

War 3 (% W;)T91—WHl

so that the approximate linear dynamics take the final form

B - w
e. = Ame1+Em eT "’ (3.18)
WHF
.-. Wm
9 = —’)’GCmele . (3.19)
IIWHF

 

In Section 4.3, the validity of this approximation method will be discussed.

2]

4. Averaging Analysis

4.1 An Identification Problem

In [6], Tan and Khalil present an identification problem containing both linear and
bilinearly coupled parameters, as in the closed-loop error dynamics in (3.14)1 . This first
step in the system analysis has promising results, and is repeated here to motivate further

analysis. Consider the error

ez(t) = 2(t) - z(t) (4.1)

and corresponding signals
z(t) = 93wa(t) + Ong(s) [907%,] (t) (4.2)
2(t) = éfwao) + 9364s) [éZwa] (t), (4.3)

where 9,, E R“, wa E R“, 9;, E R"b and Gb (3) is a stable transfer function with dimension
nb. In the original paper there is also a scalar transfer function Operating on the right
hand side of both equations. This transfer function is omitted here because it is shown
therein that the transfer function does not complicate the analysis when slow adaptation is
implemented and the time scales of the parameter adaptation and the transfer function are

sufficiently separated. The parameter adaptation rules are

)0

a = —yaez(wa+Gb(s) [wZ] 9,.) (4.4)
9b = —7besz(s)[90Twa], (4.5)

 

1In order to clearly see the bilinearly coupled parameters in (3.14), one can make a substitution using (3.11).

22

where O < )2, << 7), << 1. To simplify the notation, assume n), = l. The following

 

matrices
m .9. avg(wa(t)w5<r)),
an 9-. avg (Gus) [we] malls) [wz (0]),
A—BT 2 avg (Wa(t)Gb(S) [W5 (t )1) ,
or 2 avg(cb<s)lw.l<r)w5(r)),
(A+9bB)BT 3 .7177 + 9835?,
m 2— avg (mo) + 9.0.8) [w] (0) (war) + sells) [w] (of)

are assumed to exist for arbitrary’ to, where avg(-) = Tlim .1. 120+T(.)dt. For T-periodic
regressor wa(t), all of the matrices are guaranteed to exist. Theorem 5.1 of [6] states that

for the system with error dynamics (4.l)-(4.3) and adaptation rules (4.4) and (4.5), if

93 13—370.. > o (4.6)
(A + 9,3) BT 8.8,,T (A + 9.3) BTT
93231379,

 

 

 

(A + 8.19)2 — > 0 (4.7)

and if Q, (t = O) is sufficiently close to then 9., all parameter estimates will converge to their
true values. The proof is straightforward and uses averaging analyses and linearization of
the slow dynamics. This will give zero output error and complete knowledge of the system

parameters.

4.2 The Local Approximate Closed-Loop System

The identification problem in Section 4.1 is a good motivational example, but the
most interesting aspect of the problem is the behavior of the closed-loop system. In the
closed-loop control case, the regressor vectors cannot be arbitrarily chosen. They will

depend on feedback, and the regressor w” (t) depends nonlinearly on internal signals, as

23

discussed in Section 3.3. Nonlinear analysis of the closed-loop system has not been thor-
oughly completed, but analysis of the approximate system in a sufficiently small neigh-
borhood of zero parameter error is presented. These results are then used to construct a
converse Lyapunov function and show local stability Of the true closed-loop system dy-
namics.

Two-time scale averaging theory is now applied to the approximate system dynamics
given in (3.18) and (3.19). Equation (3.18) is the fast state, while (3.19) is the slow state, as
dictated by the adaptation gain 76. The steady-state value of e1 with the parameter estimate

error 9 frozen is

I
v (1,6) 2 f0 amt-1% [wgn oz...) dré.

The function
.. A Wm
f’ (t,9,v) = Cvam
YWHF
is defined for more compact notation. These assumptions are made:

- (Bl)

1 lo+T W t B
Fwé lim — G... "' cm/ eAmI‘-‘)—"1[w;,w,l,p]ds d:
TENT t° TWHF o 93

exists uniformly in to.

- (32)

3f' (t,§,v) —F
aé av

has zero average value, with convergence function C (T), where C() : R+ —> IL. is

strictly decreasing and

1 to+T ~
ll? f,(Tr9rv)dT-'Fav SC(T)

to

 

 

24

Theorem 4.1 Consider the error dynamics (3.18) and the adaptation rule (3.19). If as-
sumptions (BU-(B2) are satisfied, and E = O is an exponentially stable equilibrium point
of the averaged system,

5=—ma§ up

then (cl = O, 5 = 0) is an exponentially stable equilibrium of the system (3.18)-(3.19) for

75 sufiiciently small.

Proof. From its structure, f’ (t,9,v) is locally Lipschitz in (By). Moreover, since

w... and WHF are piecewise continuous in t and bounded, f’ is continuous in t. Define

Clearly, h(t,0) = O for all t > 0 and ”@3952” is bounded for all t > O. Combining the
above and assumptions (B1)-(B2), the conditions of Theorem 4.4.3 of [8] are satisfied. El
It remains to be shown that F = O is an exponentially stable equilibrium of the aver-

aged system. The averaged system is given by (4.8). Notice that

Cmv (1,6) = 53—6... (p5,, Wm)?

By assumption (B1), these quantities all exist,

 

W 2 avg (Gm(lwml)Gm (14.1)), (49>
HHT 9-— avg (GmflWHFDGmaWiiFD), (W)
W 2 avg(Gm<lwml)Gm(l-wzpl)), (4.11)
W 2 avg (Gm(lWHrl)Gm MD)- (W)

25

Substituting the above quantities into (4.8) the averaged dynamics are

 

_—=_ 70 MM 7 MHT .2.
9 = __ 6 (4.13)
93 yHMT yHHT
Within the averaged system, the parameter estimate dynamics can be further separated into

two time scales. The averaged dynamics, separated by time scales, are

 

 

 

 

3,, = —Z‘1(HHT6H+HMT66) (4.14)
93

i _ _7_G T~ TU

96 _ 93 (MI-I 9H+MM 90), (4.15)

where 35 and F” are the averaged parameter estimate vectors. As a result of averaging,
the system (4.14)-(4.15) is autonomous, and a two-time scale, slow and fast manifold
analysis [13] can be performed. The system will be transformed into the slow time scale,

1' = (76)// 93) t, for convenience:

aai: = — (HHTéH +HMT60) (4-16)
are ———.~ —T:
73? = —(MH 0H+MM 96). (4.17)

The parameter '1/ dictates that the state vector 50 is the fast state and that DH is the slow

state.

Theorem 4.2 (Theorem 3.1 from [9]) Define the slow and fast system matrices as

 

 

|li>

_l —
9, (HH7 —HMT (MMT) MHT) (4.18)

||l>

of (W). (4.19)

_ —1
If (MM T) exists, then as 7/ ——+ 0 the first nH eigenvalues of the system (4.16)-(4.I7) tend

26

to positions in the complex plane defined by the eigenvalues of 9,, namely,
24(95), 1': l,2,...,nH; (4.20)
while the remaining 2n — 1 eigenvalues of the system (4.16)-(4.I 7) tend to infinity, with the

rate of 1/'/, along asymptotes defined by the eigenvalues of D f, namely

1
?A,-(Of), j=1,2,...,2n—l. (4.21)
Furthermore, if the n” eigenvalues of Q, are distinct and the Zn — 1 eigenvalues of 52 f are

distinct, where 2.,- ($23) = ll]- (Qf) is allowed, then the eigenvalues of the original system

(4.16)-(4.I 7) are approximated as

 

2,-=i.-(O,)+0(y), i = 1,2,...,nH; (4.22)
t=lj<flf)y+0(y), i = nH+j,j=l,2,...,2n-1, (4.23)

where 0(1/ ) denotes a term on the order of magnitude of y.

Equations (4.16), (4.17), (4.18) and (4.19) are well known in singular perturbation
analysis of linear systems [9]. Using the results of Theorem 4.2 in Theorem 4.1 shows
that the origin (el = 0, 9 = O) is an exponentially stable equilibrium point of the approxi-
mated system, (3.18)-(3.19), provided that both 9,- and Q, are positive definite (PD) and
y’ is sufficiently small. The next theorem shows that (4.19) being PD is not a stringent

condition.

Theorem 4.3 Consider the matrix

 

_l __
c, = HHT — HMT (MMT) MHT (4.24)

27

where HHT E R’Wxn”, HMT E R"”"2"'1, MHT E len'lx’m andMMT E Rzn—lxzn‘l are

defined as in equations (4.9) through (4.12). Then, for all x E R” 7E 0, xTQSx 2 0.

Proof. Choose x E R"”. Since the vector x is a constant, it can be moved inside the

integral present in the averaging, i.e.

1 m+T 1 m+T
T ° ._ . = I — T ' .
x [711—124 fro ( )dt] x 1113; T to x ( )xdt
Let avg (xTHHTx), avg (xTHMT) and avg (MHTx) be denoted fl, hMT and m, respec-
tively. This gives

__ _ __ -l _
xTc,x = h2 — hMT (MMT) Mh. (4.25)

It is easy to show for z E R" and A = AT 6 RM", zTAz = trace [zzTA] , so we have

 

£52.12 = fi—trace [Kr—h hMT (MMT 1)] .

Notice that zzT is a singular matrix with rank 1, which implies that for a full rank matrix

A, zzTA is also a singular matrix with rank 1. Therefore we have

hZ—tracePl—lhhMT (MMT l)]=fi—am[mhMT (MMT 1)].

 

 

From [14] p. 318,

 

— — — TM—h—hMT
hZ—Amax[MhhMT(MMT1)]=h2— max {y _y .

 

yew-“{0} yTMMT)’

Now if one considers the averaging operator as an inner product and using the Cauchy-

Schwarz Inequality,

_ THE W7 — — TMMT
h2 - max y __.__ y 2 h2 — h2 max y_____y
y€R2"‘l\{0} yTMMTy y€R2"—1\{0} yTMMT)’

= P — 727 = 0. (4.26)

 

28

Taking the first equation (4.25) and last inequality (4.26), it is clear that {9.x 2 0. CI

4.3 Local Lyapunov Stability

In the previous section, it was shown that the approximate local system has an ex—
ponentially stable equilibrium point at the origin. In order to rigorously show that this
analysis implies exponential stability of the nonlinear system, a Lyapunov function will be

constructed. Theorems 4.1 and 4.2 together state that for the approximate local system

 

 

 

 

_ B _ -
A’" 70:1 [ W; "’le ]
é1 e
1 = w". _ , (4.27)
9 —yG f5 34.41—an d’tCm o 9
YWHF
(e17, 9T) = O is exponentially stable. These quantities are defined for future use:
A a(s)
”’"F ‘ (WI mfg“
T
~ A a s _ _
WHF = (Ti MI -th]) 91-(WH —WH)
Equation (4.27) is the local approximation of
. B... 3m T “T ”T
el = Am + E M(GgQ) e1 + a (chm-I- 91 (W11-W1)-— GHWH),
A W11
9 = —yGCmele , (4.28)
W2
I _ rum _ )

 

 

a so-calledfunctional differential equation (FDE). Recall that the signal pr is a function
of the infinite-dimensional history of the state (e17, ST) and Wm: is the model version of

this signal, i.e. how the signal evolves when the true system parameters are known and the

29

state is identically zero. Adding and subtracting convenient terms to the ri ght-hand side of

 

 

 

 

 

(4.28) gives
A... is? [ W5: will? i
é] el
5 ‘70 fl; CmeAmO—VBm Wm dTCm O é
YWHF
.. I
0 £3,” [ 0T WIT” ]
el
+ t 0 ~
”70 lo 6"” 0—103»: dTC’" O 0
YWHF
%(65Qe1+é11.%% [égwu] — aim W516”)
( Cpell
+ 812 + a: 67W” (4.29)
-?’GCmele 1% [ H ]
313
T
a s T ~
_ v((iiw) 91)) .

 

In order to proceed with analysis of the complete nonlinear system, the convolution oper-
ations in (4.29) must be realized in state-space form. Recall that Gm(s) is a single-input
single-output system and that

Wm

I
/ Cme4m(“‘)B,,. dr
0 YWHF

is actually the parallel action of Zn — 1 + my instances of this system. The state-space

realization Of this operation is

(4.30)

UN)
I
:1:-

in

UN)
.1.
W

n

(4.31)

W)
II
«9
if“

30

where

 

 

.. 1 - .
0 Am 0 0 EM 0
AC = r8; = a
i- 'I
Cm 0 0
0 C 0
c; _ "'
0 O Cm

 

 

Assume the steady-state solution of (4.30), denoted C (t), exists and has state-space repre-

sentation

Wm
_ AC; '1' BC
YWHF

QC,

and let 5 = I: — C. Now to illustrate the idea, the first term in (4.29) is derived in the

augmented state-space:

 

 

 

 

 

 

 

 

B _ -
, A... a: I 4. win]
e1 81
.:. : Wm ~
9 —ya féCme4m(’—")Bm drC... 0 9
YWHF J
- 1 _ 9 _ _ - - ,
é1 Am 7.? w; WIT-IF 0 81 0
=> 5 = —75C§§C.., O O 5 -YGCc5Cmel
_ 5 _ o o A; _ _ 5 _ 0 _

31

 

Continuing in this manner on the second and third terms in (4.29), grouping similar result-

ing terms, and augmenting another new state for §—% will lead to the complete expression
x = A(t)x+f(t,x) +g(t,x,) (4.32)

where x is the new augmented state and x, denotes a function of its history. Once favorable
properties of (4.32) are established, the stability analysis can proceed using a variation of
typical Lyapunov techniques. See for instance [15]. Establishing the asymptotic stability
of (4.32) also validates using (4.27) as a local approximation of the nonlinear system. Until
this point it is accepted that substituting W” for w” in (3.14) and (3.17), then linearizing
the standard way, via partial derivatives, results in a legitimate estimate. The first term in
(4.32) is merely an asymptotically stable, linear time-varying (LTV) matrix.

Now the stability of a system of the form (4.32) described by an FDE with desir-
able properties is shown. Define the Banach space <I> consisting Of the set of continuous

functions (1) : R+ -—> R“, equipped with
||¢II = sup||¢(t)llo
:20
Consider the asymptotically stable, LTV system
at = A(t)x, x E Rn" (4.33)

with A(t) continuous in t and bounded. From [13] Theorem 4.12, for a continuous matrix
Q(t) satisfying
0 < c1] < Q(t) < czI, v: > o, Q(t) = QT(t)

32

there exists a continuously differentiable matrix P(t) satisfying

0 < c3] < P(t) < c4], Vt > 0, P(t) = PT(t),

—P(t) = P(t)A(t) +AT(t)P(t) + Q(t).
Equation (4.32) is a perturbed version of (4.33). Let

f : W X 131"" -+ W" ||f(t,X)|| S Klletllzi

8 : R" X G -+ 4’ |I8(tixt)|| S PIIXzIIHXII G C ‘1’-

Now, for the perturbed system (4.32), choose the Lyapunov function candidate V(t,x) =
xTP(t)x. The derivative of V along (4.32) is

V(t,x) = xTPU)(f(tiX)+g(tixt))+(f(tax)+8(tiX:))TP(t)x
+ x1 (P(t)A(t) +AT(t)P(t) +P(t))x
= —xTQ(t)x+2(f(t,x)+8(‘:xt))TP(‘)x
=>V S -¢‘lIIJC||2+2€4KIIJ'C||3*I‘ZCMDIIMIIIIJCIIZ+26405||JC|I2

= —(c1—2C4K||x||-2(:4p||x,||)||x||2. (4.34)
Define the neighborhood <I>x = {x, E <I>: ”x,” < min{c1 / 8C4 K,C1 / 8C4p}}. For x, e (bx

v g —(c1—ZC4KIIXII- 2C410||x:||)IIXII2
< c21x” ||2<--—'V

=>V < ——V, (4.35)
2C4

33

where c1/2c4 > 0. The comparison system

. Cl A
y 264v (y) ( )

is introduced. This system is clearly exponentially stable, as it is merely a linear time-
invariant (LTI) system with lone eigenvalue 2., = —c1/2C4 > 0. Comparing (4.35) and
(4.36) gives

V < h(V). (4.37)

Theorem 6.6.1 from [15] states that when (4.37) is true, that the exponential stability of
(4.36) implies the exponential stability of the original system (4.32).

The above conclusion pertains to (4.32). The last step to verify (4.28) is locally ex-
ponentially stable is to show, when all convolution integrals are eliminated by state vector
augmentation, that (4.28) indeed takes the form (4.32), with the same boundedness and
continuity properties. Continuity properties of the Preisach operator and its right inverse
are discussed in [16]. These results will be useful in showing the continuity of w”; as
in (4.29). The second part of (A5) ensures a sufficient condition on the Preisach density
function therein. Continuity and boundedness properties of many of the signals will follow
from boundedness and continuity of the driving signal, along with the stability of G... (s)

(l
and“?

Remark 4.1 Note that in the context of Section 4.3. |

 

x“ -—> 0 does not imply ||x.|| —-+ O. In
fact, if at any single point in time ”x” 96 0, then ”x. II will always be nonzero, regardless of

what happens to ”x” thereafter.

Remark 4.2 Assumption ( A8) considers d” (t) in (2.11) to be identically zero. Relaxing
this assumption should not complicate the analysis. Since d” (I) will always be bounded,
the output error should be ultimately bounded, with ultimate bound proportional to an

upper bound on d” (t).

34

Remark 4.3 For sums of bounded, periodic signals, assumptions (BU-(B2) will always
be satisfied [13]. Periodic driving signals are common in many applications, e. g. [ I 7], so

often these are not stringent conditions.

35

5. Simulation Results

Simulation results can now be introduced in order to better understand the behavior
of the closed-loop system and to motivate further work. Several different aspects of interest
will be explored. The simulations are performed on the complete closed-loop system, using
software provided by Dr. Xiaobo Tan to simulate the Preisach operator and its approximate
inverse. The algorithm for the inversion of the Preisach Operator in [3] is an iterative
algorithm, so it must be implemented in discrete time. Therefore, all simulations in this
thesis are implemented in discrete-time, and parameter updates are performed according

to the first-order approximation of (3.16):
906+ 1) = 6(k) - roe(k)Gd(k),

where e(k) and Gd (k) are the sampled data versions of Cmel and

l ' a... ')
Wll

Gm . (5.1)

(v((iiillwllYél—wl.) _ )

 

 

 

 

respectively.

In Section 5.1 the persistence Of excitation of driving signals is examined. A vector
v() : R+ ——> R” is said to be persistently exciting (PE) of order n if there exist a1, a2, 5 > 0
such that

to+6
a2] 2 / v(r)vT(r)dr 2 all (5.2)
to

for all to 2 0, where I is the n x n identity matrix [8]. Roughly speaking, a signal that
is PE of order n can drive n adapted parameters to their true values in an MRAC system.

The relationship between the initial conditions of the plant and hysteresis parameters is the

36

topic of Section 5.2. Let qJ(t,xo) be a solution to

x=f(X).

with an asymptotically stable equilibrium point x = 0. Then, the region of attraction is

defined as the set of all points x0 such that (p(t,xo) is defined for all t 2 0 and

lim (p(t,xo) = 0

t—ooo

in [13]. The main point of Section 5.2 is to illustrate how the initial estimate of the plant pa-
rameters affects the possible initial estimates of the hysteresis parameters, and vice versa.
In addition, the region of attraction of zero output error and zero parameter error is in-
vestigated for the case when the hysteresis and plant parameter adaptation time-scales are

sufficiently separated and when they are very near each other.

5.1 Parameter Convergence

Persistence of excitation of driving signals is well understood for open-loop identi-
fication problems. For linear systems with n parameters to be identified, the driving signal
must have n unique nonzero spectral components (provided Zp (s) and Pp (s) are coprime)
[18]. PE conditions for linear MRAC systems depend instead on the internal regressor
vectors, which cannot be explicitly chosen. A rule of thumb is to require the driving signal
to have as many nonzero spectral components as the number of parameters to be estimated
(the sum of the lengths of all regressors) [12]. Ioannou and Sun give sufficient conditions
to guarantee parameter convergence when this type of driving signal is applied. In ad-
dition to assumptions (A1)-(A4), G...(s) must also be strictly positive real for a classical
MRAC with relative degree n* = 1. For an increasingly complex linear system (larger
number of parameters to be identified), the PE condition will require an increasing number

of components in the driving signal.

37

 

Simulation Parameters

 

 

 

6(3) ___ 11213+91650

Linear Dynarmcs s2+l475s+91650

 

0(2) = 22422—18754

Discrete Approxrmatlon 24.95.4533.)

 

 

 

Driving Signal ymf = 22 sin(27t25t) + 22 sin(27t50t)
Sampling Period T = 2 x 10‘3
Adaptation Gains yo = 5 x 10-5, 7” = 1 x lo-2

 

 

 

 

Table 5.1: Parameter convergence simulation settings.

In [7] Tan and Baras introduce PE conditions in the identification of hysteresis pa-
rameters. Instead of depending on the number Of spectral components, the hysteresis input
signal (v(t) in Fig. 2.4) must contain a sufficient number of reversals at specific levels. In
other words, each discretization level should contain a local extremum of the signal v(t).

PE conditions for the proposed closed-loop controller are not yet completely under-
stood. Refer to the simulation parameters in Table 5.1. The amplitude of the driving signal
is chosen in order to drive the hysteresis Operator very near its upper and lower saturation
levels. This is to ensure that each level in the truncated Preisach plane is reached. For each

simulation in this Section 5.1, the parameter errors will be calculated as

2n—l
80(1) = 100 ( )and
i=1

"H A . _ .
em) = 100 2 _—
The first simulation case has discretization level L = 4 and initial parameter estimate 9 =

90:0) — 90:

 

 

 

 

1.059. The matrices (If and Q, are calculated numerically with sampling period T =
2 x 10-8, which gives a million samples for every period Of the driving signal. For L = 4,
both of these matrices have full rank, i.e. rank[§2f] = 3 and rank[£2,.] = 11. Fig. 5.1 shows

38

 

 

 

 

 

 

 

 

 

20 e .
—eG(t)
to ...e 1
§ 15’ H” -
8
95. 10- 4
d)
E
E
3 5~
5“.
4M...“ - ___ 1 4
G0 1 000 2000 3000 4000 5000

iteration k/100

Figure 5.1: Parameter error for case L = 4.

that when the plant and hysteresis parameters are updated using the adaptation gains in
Table 5.1, the signal drives all parameter estimates to their true values. Thus, for this case
the above approach produces a persistently exciting driving signal.

The implication of the previous example is significant. In this particular case, there
are only two frequency components in the driving signal and all 14 parameter estimates
reach their true values. For L = 8, on the other hand, 52, is not positive definite, in partic-
ular, rank[.Qs] = 30 while my = 37. From Fig. 5.2, it is clear that the parameter estimates
do not converge to their true values, but rather some other steady-state. It is encouraging
to note, however, that while the parameter estimates do not converge, the output error,

computed as

eL(t) =108|y(t)-ym(t)|

continues to approach zero as time approaches infinity, as in Fig. 5.3. The downward spike
is when the error y(t) — y... (t) switches signs. This means that under appropriate conditions
it may be possible to guarantee perfect tracking in the absence of a PE driving signal, which

makes the proposed controller much more practical.

39

20 . . . .

 

 

 

 

 

 

 

 

 

 

 

_eG(t)
a) ---e 1
E 15 H( ) .
3
§
0)
E l ______________________
9 ...........................................
CU )-
Q.
G0 10100 20100 30100 4000 5000
iteration k/2000
Figure 5.2: Parameter error for case L = 8.
§
8
‘5
.9-
3
O

 

 

 

 

0 10100 2000 3000 4000 5000
iteration k/2000

Figure 5.3: Output error for case L = 8.

5.2 Region of Attraction

In this section, the effects that the separation of time-scales and initial parameter es-
timates have on convergence are investigated. Since it is impossible to completely classify
the region of attraction of the 2n — 1 + nH-dimensional system, the problem is restricted to
a two degree-of-freedom (DOF) problem instead. Consider two design parameters a and

b. The initial parameter estimates are governed by these parameters as

00(0)=a9(; and 9H(0)=b9... (5.3)

40

When a = b = 1 the initial parameter estimates are the true parameter values, and the
output error is identically zero. Fig. 5.4 shows the region of attraction in the a — b plane.
Two cases are shown: the solid line is for m = 10‘2 and the dashed line is for m = 10".
In the first case, the hysteresis parameters evolve at a rate that is approximately ten times
slower than the plant parameters. This is considered a sufficient separation of time scales
for this particular system, as decreasing 7).) further does not have a significant impact on
the region of attraction. In the second case, the hysteresis parameters evolve at nearly the
same rate as the plant parameters. It is clear from Fig. 5.4 that while there seems to be
little regularity or symmetry to the encircled regions, the separation Of time scales brings
the added benefit of a larger region of attraction. In fact the (dimensionless) area of each

encircled region can be easily computed:

AREAsolid = 3.9453 and AREAdashed = 2.5122.

In this case, the separation of time scales increases the two DOF region of attraction over
57%. This snapshot of the respective regions of attraction are in no way the true maximal
regions of attraction, which can only be represented as a 2n — 1 + nH-dimensional volume.
However, even in this simplified case the benefit is clearly evident.

It is also worth noting that the region in the a — b plane to the immediate right of the
encircled regions still gives output error stability, but the parameters do not converge their
true values. In the rest of the a - b plane, outside of the encircled areas, the parameter esti-
mates and the output error are unstable. As such, one cannot generalize that the complete
nonlinear system will always be stable, without placing sufficiently small limits (e.g. via
parameter projection) on the parameter estimates.

From Fig. 5.4 it is clear that in this simplified investigation of the region of attraction,
that both an adequate separation of time scales and and a more precise initial estimate of the

plant (hysteresis) parameters allows more variation of the initial estimate Of the hysteresis

41

parameter b
fl
.— 'ut

.9
Ut

 

 

 

 

 

 

 

 

(plant) parameters.

 

Figure 5.4: Comparison of regions of attraction.

42

 

6. Future Work

There is still much work to be done regarding this particular control scheme. The
analysis herein is based completely on “linearizing” the error and adaptation dynamics,
then performing averaging. The local exponential stability via a Lyapunov argument is
based on the linear system’s stability being proven by averaging. In the identification
problem in [6] and repeated in Section 4.1, averaging is applied to the complete nonlinear
system. Linearization is only used to show the stability of the averaged nonlinear sys-
tem. This allows the authors to estimate the region of attraction. In general, averaging
the closed-loop nonlinear system in Section 4.2, rather than linearizing prior to averaging,
may provide more useful insight into the system behavior. This will require a thorough
understanding Of the signal wH, so that a new averaging theorem can be derived, or the
work regarding nonlinear averaging by Tee], Moreau and Nesié in [19] or the work re-
garding averaging with hysteresis by Pokrovskii, Rasskazov and Vladimirov in [20] can
be applied.

The PD requirement on {If and S2, only guarantee parameter convergence in a neigh-
borhood of the origin. Finding an analytical expression for the region of attraction will
demonstrate how practical this control scheme will be in systems with poor initial param-
eter estimates. Additionally, one may investigate whether the region of attraction is in any
way dependent on the magnitude of separation of time scales. In this thesis, the aforemen-
tioned PD conditions were satisfied in an ad hoc manner. A signal was chosen, and then
the matrices were computed to show whether or not the signal would lead to convergence.
It will be of interest to find sufficient conditions on the driving signal that will lead to both
matrices being PD. In Section 5.1, simulations showed that when conditions (4.18) and
(4.19) are not satisfied, the output error may still converge to zero. This needs rigorous

analysis to establish under which conditions output error convergence can be guaranteed.

43

References

[1] X. Tan, “Control of smart actuators,” Ph.D. dissertation, University of Maryland,
College Park, 2002, available at http://www.egr.msu.edu/~xbtan/papers.html.

[2] X. Tan and J. S. Baras, “Modeling and control of a magnetostrictive actuator,” in
Proceedings of IEEE Conference on Decision and Control, 2002, pp. 866—872.

 

, “Modeling and cOntrOl of hysteresis in magnetostrictive actuators,” Automat-
ica, vol. 40, no. 9, pp. 1469—1480, 2004.

[3]

[4] G. Tao and P. V. Kokotovic, “Adaptive control of plants with unknown hysteresis,”
IEEE Transactions on Automatic Control, vol. 40, no. 2, pp. 200—212, 1995.

 

[5]

, Adaptive Control of Systems with Actuator and Sensor Nonlinearities. New
York: John Wiley & Sons, Inc., 1996.

[6] X. Tan and H. K. Khalil, “Control of unknown dynamic hysteretic systems using
slow adaptation: Preliminary results,” in Proceedings of the 26th American Control
Conference, New York, NY, 2007, pp. 3294-3299.

[7] X. Tan and J. S. Baras, “Adaptive identification and control of hysteresis in smart
materials,” IEEE Transactions on Automatic Control, vol. 50, no. 6, pp. 827—839,
2005.

[8] S. Sastry and M. Bodson, Adaptive Control. Englewood Cliffs, NJ: Prentice Hall,
1989.

[9] P. Kokotovié, H. K. Khalil, and J. O’Reilly, Singular Perturbation Methods in Con-
trol: Analysis and Design. Philadelphia: Society for Industrial and Applied Mathe-
matics, 1999.

[10] I. D. Mayergoyz, Mathematical Models of Hysteresis and Their Applications. New
York, NY: Elsevier, 2003.

[l 1] J. J. Reynolds, X. Tan, and H. K. Khalil, “Closed-loop analysis of slow adaptation
in the control of unknown dynamic hysteretic systems,” in Proceedings of the 46th
IEEE Conference on Decision and Control, New Orleans, LA, 2007, to appear.

[12] P. A. Ioannou and J. Sun, Robust Adaptive Control. Prentice-Hall, 1995.

[13] H. K. Khalil, Nonlinear Systems, 3rd ed. Upper Saddle River, New Jersey: Prentice
Hall, 2002.

[14] D. S. Bernstein, Matrix Mathematics: Theory, Facts and Formulas with Application
to Linear Systems Theory. Princeton, New Jersey: Princeton University Press, 2005.

[15] A. N. Michel and K. Wang, Qualitative Theory of Dynamical Systems. New York,
NY: Marcel Dekker, Inc., 1995.

[16]

[17]

[18]

[19]

[20]

R. V. Iyer, X. Tan, and P. S. Krishnaprasad, “Approximate inversion of the Preisach
hysteresis operator with application to control of smart actuators,” IEEE Transactions
on Automatic Control, vol. 50, no. 6, pp. 798-810, 2005.

D. Croft, G. Shed, and S. Devasia, “Creep, hysteresis, and vibration compensation for
piezoactuators: Atomic force microscopy application,” Journal of Dynamic Systems,
Measurement, and Control, vol. 123, no. 1, pp. 35—43, 2001.

B. D. 0. Anderson, R. R. Bitmead, C. R. Johnson, Jr., P. V. Kokotovic, R. L. Kosut,
I. M. Y. Mareels, L. Praly, and B. D. Riedle, Stability of Adaptive Systems: Passivity
and Averaging Analysis. Cambridge, MA: The MIT Press, 1986.

A. R. Teel, L. Moreau, and D. Nesic’, “A unified framework for input-tO-state stability
in systems with two time scales,” IEEE Transactions on Automatic Control, vol. 48,
no. 9, pp. 1526—1544, Sept. 2003.

A. Pokrovskii, O. Rasskazov, and A. Vladimirov, “Averaging principle for differential
equations with hysteresis,” University College Cork, Ireland, Tech. Rep., 2004.

45

 

I”lllllllililliii