Pseudo-inverse-based adaptive control for uncertain discrete time systems preceded by hysteresis

Automatica 45 (2009) 469–476

Contents lists available at ScienceDirect

Automatica journal homepage: www.elsevier.com/locate/automatica

Brief paper

Pseudo-inverse-based adaptive control for uncertain discrete time systems preceded by hysteresisI Xinkai Chen a,∗ , Takeshi Hisayama b , Chun-Yi Su c a

Department of Electronic and Information Systems, Shibaura Institute of Technology, Minuma-ku, Saitama-city, Saitama 337-8570, Japan

b

Control Systems Group, IHI Scube Co., Ltd., Japan

c

Department of Mechanical and Industrial Engineering, Concordia University, Montreal, Quebec, H3G 1M8, Canada

article

info

Article history: Received 24 July 2007 Received in revised form 28 April 2008 Accepted 4 August 2008 Available online 9 December 2008 Keywords: Hysteresis Prandtl–Ishlinskii model Adaptive control Sliding mode control Disturbance

a b s t r a c t This paper discusses the adaptive control for the uncertain discrete time linear systems preceded by hysteresis nonlinearity described by the Prandtl–Ishlinskii (PI) model. The contribution of the paper is the development of an adaptive algorithm in which a pseudo-inversion is introduced to avoid difficulties of the directly inverse construction for complex hysteresis models, especially for the unknown hysteresis case. In the developed approach, only those parameters in the formulation of the sliding mode controller are adaptively estimated. The stability in the sense that all signals in the loop remain bounded is analyzed. Simulation results show the effectiveness of the proposed algorithm. © 2008 Elsevier Ltd. All rights reserved.

1. Introduction The hysteresis phenomenon can be found in diverse disciplines ranging from smart materials (Banks & Smith, 2000; Moheimani & Goodwin, 2001; Webb, Lagoudas, & Kurdila, 1998), to ferromagnetism and superconductivity (Mayergoyz, 1991), to economics (Cross, Krasnosel’skii, & Pokrovskii, 2001) and geosciences (Guyer, McCall, & Boitnott, 1994). When a plant is preceded by hysteresis, the system usually exhibits undesirable inaccuracies or oscillations and even instability (Tao & Kokotovic, 1995). The development of control techniques to mitigate the effects of hysteresis has been studied for decades and has recently re-attracted significant attention, e.g. Moheimani and Goodwin (2001) and the references therein. Much of the interest is a direct consequence of the importance of hysteresis in numerous new applications. Interest in studying dynamic systems with actuator hysteresis is also motivated by the fact that they are nonlinear systems with nonsmooth nonlinearities for which traditional control methods are

I This paper was not presented at any IFAC meeting. This paper was recommended for publication in revised form by Associate Editor Gang Tao under the direction of Editor Miroslav Krstic. ∗ Corresponding author. Tel.: +81 48 687 5805; fax: +81 48 687 5198. E-mail addresses: [email protected], [email protected] (X. Chen), [email protected] (T. Hisayama), [email protected] (C.-Y. Su).

0005-1098/$ – see front matter © 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.automatica.2008.08.004

insufficient and thus require development of alternate effective approaches (Chen, Su, & Fukuda, 2008; Su, Stepanenko, Svoboda, & Leung, 2000; Su, Wang, Chen, & Rakheja, 2005; Tao & Lewis, 2001; Tan & Baras, 2004; Zhou, Wen, & Zhang, 2004). The development of a general frame for control of an uncertain dynamical system in the presence of unknown hysteresis and disturbance is quite a challenging task. To deal with the control problem of a dynamical system preceded by hysteresis, the thorough characterization of the hysteresis nonlinearities forms the foremost task. Appropriate hysteresis models may then be applied to the formulation of control algorithms. Hysteresis models can be roughly classified into physics-based models and phenomenological models. Physicsbased models are built on first principles of physics (Jiles & Atherton, 1986). Phenomenological models are used to produce behaviors similar to those of the physical systems without necessarily providing physical insight into the problems. The phenomenological models are formulated by the weighted aggregate effect of all possible so-called elementary hysteresis operators. Elementary hysteresis operators are noncomplex hysteretic nonlinearities with a simple mathematical structure. The popular phenomenological models are the Preisach model (Adly, Mayergoyz, & Bergqvist, 1991; Croft, Shed, & Devasia, 2001; Mayergoyz, 1991; Natale, Velardi, & Visone, 2001), Prandtl–Ishlinskii model (Brokate & Sprekels, 1996; Visintin, 1994), and Krasnosel’skii–Pokrovskii model (Krasnosel’skii & Pokrovskii, 1989; Visintin, 1994). The

470

X. Chen et al. / Automatica 45 (2009) 469–476

Preisach model and Krasnosel’skii–Pokrovskii (KP) model are parameterized by a pair of threshold variables (Mayergoyz, 1991), whereas the Prandtl–Ishlinskii (PI) model is a superposition of elementary stop operators which are parameterized by a single threshold variable (Visintin, 1994). With the developments in various hysteresis models, it is natural to seek means to fuse these hysteresis models with the available control techniques to mitigate the effects of hysteresis, especially when the hysteresis is unknown, which is a typical case in many practical applications. The most common approach in coping with hysteresis in the literature is to construct an inverse operator, which is pioneered by Tao and Kokotovic (1995), and the reader may refer to, for instance, Huhnen and Janocha (1999), Krejci and Huhnen (2001) and Tan and Baras (2004) and the references therein for the further development. Essentially, the inversion problem depends on the phenomenological modeling methods. Some hysteretic nonlinearities are very complicated with multi-values and nonsmooth features, such as piezo-electric actuators and magnetostrictive actuators, where phenomenological hysteresis models are required to describe the hysteretic nonlinearity. In this case, the analytic inversions of the phenomenological hysteresis models are generally very complicated, causing the difficulty of stability analysis for the hysteretic system. Focusing on the phenomenological Prandtl–Ishlinskii (PI) model, this paper proposes a new approach to deal with the inversion problem. Instead of directly constructing the inversion from the phenomenological hysteresis model, a so-called pseudoinversion is introduced. The inversion is obtained by searching the optimal value of the inversion based on a performance index, thus avoiding the mathematical difficulties in the inversion construction, especially for the unknown hysteresis case. Such a pseudo-inversion technique is developed for uncertain linear discrete time dynamical systems preceded by the hysteresis with the phenomenological PI representation. The considered uncertain linear discrete time dynamical system contains unknown parameters and disturbance, where the disturbance is composed of model uncertainties, nonlinearities, external disturbances, etc. (Chen, 2006). In the developed approach, only those parameters in the formulation of the controller are adaptively estimated online using adaptive algorithms with dead zone. The adaptive sliding mode controller is synthesized by using the estimated parameters. The stability in the sense that all signals in the loop remain bounded is analyzed, and the simulation results show the effectiveness of the proposed algorithm. The remainder of this paper is organized as follows. Section 2 describes the problem and the PI-type hysteresis model. In Section 3, first, the parameters in the formulation of the controller are adaptively estimated. Then, the adaptive sliding mode controller is formulated and the stability of the closed system is analyzed. In Section 4, the simulation results are presented. Section 5 concludes this paper.

called ‘‘disturbance’’ in this paper. H [v] is the hysteresis operator which will be described later; q−1 is the delay operator. d is the time delay; A(q−1 ) and B(q−1 ) are polynomials defined by A(q−1 ) = 1 + a1 q−1 + · · · + an q−n , B(q

−1

) = b0 + b1 q

−m

+ · · · + bm q

(3)

,

b0 6= 0, m < n,

(4)

and A(q−1 ) and B(q−1 ) are assumed to be relatively prime. w(k) is the hysteresis output described in (2). The control purpose is to drive the output y(k) to track a uniformly bounded signal yd (k) for the uncertain system (1) preceded by hysteresis (2). For system (1), we make the following assumptions. Assumption 1. The time delay d and the plant order n are known. Assumption 2. The plant is in minimum phase. Assumption 3. The parameters {ai , bi } in the polynomials A(q−1 ) and B(q−1 ) are unknown. Furthermore, the sign of b0 6= 0 is known. Without loss of generality, it is assumed b0 > 0.

2.2. Hysteresis model In this paper, the Prandtl–Ishlinskii (PI) model is adopted. The basic element of the PI model is the so-called stop operator ω(k) = Er [v](k) with threshold r. For an arbitrary piece-wise monotone function v(k), define er : R¯ → R¯ (where R¯ denotes the space of real numbers) as er (v) = min(r , max(−r , v)).

(5)

Suppose the function v(k) is monotone on ki ≤ k ≤ ki+1 (for i = 0, 1, . . .) with k0 = 0. For any initial value w−1 ∈ R¯ and r ≥ 0, the stop operator Er [·; w−1 ](k) is defined as Er [v; w−1 ](0) = er (v(0) − w−1 ),

(6)

Er [v; w−1 ](k) = er (v(k) − v(ki ) + Er [v; w−1 ](ki )),

(7)

for ki < k ≤ ki+1 (Brokate & Sprekels, 1996). The stop operator is mainly characterized by the threshold parameter r ≥ 0 which determines the height of the hysteresis region in the (v, w) plane. For simplicity, denote Er [v; w−1 ](k) by Er [v](k) in the following of this paper. It should be noted that the stop operator Er [v](k) is rate-independent. The PI hysteresis model is defined by

w(k) =

2. Problem statement

−1

∞

Z

p(r )Er [v](k)dr ,

(8)

0

2.1. System description Consider a system composed of an uncertain linear plant preceded by hysteresis described by A(q−1 )y(k) = q−d B(q−1 )w(k) + σ (k),

(1)

w(k) = H [v](k),

(2)

where y(k) and v(k) are the output and input, respectively; σ (k) is an unknown signal composed of model uncertainties, nonlinearities, external disturbances, etc. For simplicity, σ (k) is

where p(r ) is the density R ∞function which is usually unknown, satisfying p(r ) ≥ 0 with 0 rp(r )dr < ∞ (Visintin, 1994; Webb et al., 1998). Since the density function p(r ) vanishes for large values of r, it is reasonable to assume that there exists a constant R such that p(r ) = 0 for r > R (Brokate & Sprekels, 1996; Visintin, 1994). Thus, model (8) gives

w(k) =

R

Z

p(r )Er [v](k)dr . 0

(9)

X. Chen et al. / Automatica 45 (2009) 469–476

3. Adaptive control design

with

θ = [f0 , f1 , . . . , fn−1 ]T .

3.1. Some preliminaries

s(k + d) = C (q−1 )(y(k + d) − yd (k + d)),

(10)

where C (q−1 ) is a Schur polynomial defined by (11)

(12)

It is obvious that limk→∞ s(k) = 0 implies limk→∞ (y(k) − yd (k)) = 0. Now, consider the polynomial equation C (q−1 ) = A(q−1 )Ω (q−1 ) + q−d F (q−1 ),

(13)

) and F (q ) are in the following form

Ω (q ) = h0 + h1 q

R

Z

g1 p(r )Er [v](k − 1)dr 0

R

gm+d−1 p(r )Er [v](k − m − d + 1)dr

− ··· − 0

s(k + d) = 0.

−1

−1

g0 p(r )Er [v](k)dr = −φ (k)θ − T

Z

Let the sliding surface be defined by

−1

R

Z 0

C (q−1 ) = 1 + c1 q−1 + · · · + cn q−n .

where Ω (q

(23)

In the case that the parameters are all known and the disturbance is absent, if v(k) is chosen such that

Define the variable

−1

471

−d+1

+ · · · + hd−1 q

+ C (q−1 )yd (k + d) + δ · s(k),

(24)

where δ is the weighting factor in the range of 0 < δ < 1, then it can be easily concluded that limk→∞ s(k) = 0 and thus the output tracking can be asymptotically achieved. In the following development, it is assumed that the v(k) satisfying (24) exists for all k > 0. For the uncertainty ω(k), let

ω(k) = ω(1) (k) + ω(2) (k),

,

F (q−1 ) = f0 + f1 q−1 + · · · + fn−1 q−n+1 .

(25)

(14)

where ω(1) (k) and ω(2) (k) are respectively defined by

(15)

ω(1) (k) = Ω (q−1 )σ (1) (k),

Thus, the parameters in Ω (q−1 ) and F (q−1 ) can be determined uniquely and h0 = 1 (Astrom & Wittenmark, 1989). It is assumed that the disturbance σ (k) is composed of model uncertainties, nonlinearities, external disturbances, etc.

ω(2) (k) = Ω (q−1 )σ (2) (k).

(26)

From (17) and (26), it yields

(2) ω (k) ≤ ρ,

(1) ω (k) ≤ α · β T (k − d),

(27)

where T (k − d) is defined by

Assumption 4. The disturbance σ (k) is of the following form

σ (k) = σ (1) (k) + σ (2) (k), (1) σ (k)

(16)

+

! 21

R

XZ i=0

(Er [v](k − d − i))2 dr

,

(17a) (17b)

where α is a very small unknown nonnegative constant; ρ¯ is an unknown nonnegative constant; the regressor φ(k) and the norm kφ(k)k2 are respectively defined as

φ(k) = [y(k), y(k − 1), . . . , y(k − n + 1)]T

(18)

and

 1 kφ(k)k2 = φ T (k)φ(k) 2 .

(19)

If m is unknown, m should be replaced by n in the proposed formulation. Now, multiplying Eq. (13) with y(k) and employing Eq. (1) yield C (q−1 )y(k + d) = G(q−1 )w(k) + F (q−1 )y(k) + ω(k + d)

(20)

where G(q−1 ) and ω(k) are respectively defined by

C (q

)y(k + d) = φ (k)θ +

Z 0

(29)

i=0

3.2. Parameter estimation Since the parameters {ai , bi } of the polynomials A(q−1 ) and B(q−1 ) are unknown, the parameters in G(q−1 ) and F (q−1 ) are not available. In the following, we will try to estimate these parameters together with the hysteresis density function p(r ) by using adaptive algorithms with dead zone. Let

h iT θˆ (k) = fˆ0 (k), fˆ1 (k), . . . , fˆn−1 (k)

(30)

denote the estimate of the unknown parameter θ at the kth step and let pˆ i (r , k) be the estimate of gi p(r ) at the kth step for a fixed r. Define the estimation error as

−

R

X Z

pˆ i (r , k − 1)Er [v](k − d − i)dr .

(31)

0

e(k) = φ T (k − d){θ − θˆ (k − 1)} + ω(k)

R

g0 p(r )Er [v](k)dr + · · ·

m+d−1

+

R

gm+d−1 p(r )Er [v](k − m − d + 1)dr + ω(k + d),

(28)

Then, by substituting (22) into (31), e(k) can also be expressed as

0

+

i =0

i =0

Substituting (2) into (20) yields T

.

0

β and ρ are unknown positive constants defined by ! d−1 d−1 X X |hi | , |hi | . β= ρ = ρ¯

m+d−1

(21)

ω(k) = Ω (q−1 )σ (k). Z

(Er [v](k − d − τ − i)) dr 2

e(k) = C (q−1 )y(k) − φ T (k − d)θˆ (k − 1)

G(q−1 ) = Ω (q−1 )B(q−1 )

= g0 + g1 q−1 + · · · + gm+d−1 q−m−d+1 ,

! 12

R

XZ i =0

0

(2) σ (k) ≤ ρ, ¯

−1

kφ(k − d − τ )k22

max 0≤τ ≤d−1

m−1

m−1

≤ α kφ(k − d)k22 +

T (k − d) =

i =0

(22)

R

X Z

gi p(r ) − pˆ i (r , k − 1) Er [v](k − d − i)dr .




0

(32)

472

X. Chen et al. / Automatica 45 (2009) 469–476

Since β and ρ are unknown, their corresponding estimates ˆ k − 1) and ρ( β( ˆ k − 1) are used to construct the dead-zone function. Now, define

ˆ k − 1)T (k − d) + ρ( η( ˆ k − 1) = γ · β( ˆ k − 1),

(33)

where γ > 0 is a design parameter which is very small (Chen, 2006), the initial values of w(k) (for k = 3 − 2d − m, . . . , 0) and y(k) (for k = 3 − 2d − n, . . . , 0) can be assigned to any small values. The dead-zone function is defined as

λ(k) =

 

1−

0

if |e(k)| > η( ˆ k − 1)

R

(35)

0

i=0

The estimates θˆ (k) and pˆ i (r , k) are updated by the following adaptation laws with constraint

λ(k)e(k)φ(k − d) , 1 + D(k − d) λ(k)e(k)Er [v](k − d − i) pˆ 0i (r , k) = pˆ i (r , k − 1) + χ(k) , 1 + D(k − d)  0 0 pˆ 0 (r , k) if pˆ 0 (r , k) ≥ 0 pˆ 0 (r , k) = ˆ k) = θˆ (k − 1) + χ (k) θ(

0

otherwise

pˆ i (r , k) = pˆ 0i (r , k),

+ C (q−1 )yd (k + d) + δ · s(k) − e(k),

(45)

Remark 2. W (k) has a relation with the values of v(k − 1), . . . , v(k − m − 2d + 1). For an assigned admissible error ε , we try to find the inverse V ∗ (k) such that R

(37)

Let [vmin , vmax ] be the practical input range to the hysteresis operator, which is a strict subset of [−R, R]. Let

(38a) (38b)

nonnegative constants. Remark 1. The parameter adaptation gain χ(k), which is introduced to adjust the adaptation speed, should be chosen in the range κ ≤ χ(k) ≤ 2 − κ with 0 < κ < 1. χ(k) should be increased in the range if a fast adaptation of the parameters is needed, and vice versa. Some further restrictions for χ(k) and the initial values θˆ (0), ˆ 0) and ρ( pˆ i (r , 0), β( ˆ 0) will be discussed in Remark 5. Lemma 1. For the adaptation algorithm in (36)–(40), the following properties hold.


2

pˆ m+d−1 (r , k)Er [v](k − m − d + 1)dr

0

(36)

λ(k) |e(k)| ˆ k) = β( ˆ k − 1) + χ (k) β( γ T (k − d), (39) 1 + D(k − d) λ(k) · |e(k)| ρ( ˆ k) = ρ( ˆ k − 1) + χ (k) . (40) 1 + D(k − d) ˆ 0) and ρ( The initial values β( ˆ 0) can be chosen as very small

R

R

Z

ˆ k) and ρ( The estimates β( ˆ k) are respectively updated by the following adaptation laws

dr are bounded for

λ2 (k)e2 (k)

(P3) limk→∞ 1+D(k−d) = 0

(P4)

Z

pˆ 0 (r , k)Er [V ](k)dr − W (k) ≤ ε. ∗

0

(46)

R

Z

pˆ 0 (r , k)Er [vmin ](k)dr = W sat (k),

(47)

pˆ 0 (r , k)Er [vmax ](k)dr = W sat (k).

(48)

0

i = 1, 2, . . . , m + d − 1.

R ˆ k), β( ˆ k), ρ( (P1) θ( ˆ k) and 0 pˆ i (r , k) − gi p(r ) all k > 0. P∞ λ2 (k)e2 (k) (P2) k=1 1+D(k−d) < ∞.

pˆ 1 (r , k)Er [v](k − 1)dr

where δ is the weighting factor in the range of 0 < δ < 1, e(k) defined in (31) is added to compensate the disturbance.

T

(Er [v](k − d − i))2 dr .

R

Z 0

(34)

D(k − d) = (γ · T (k − d)) + φ (k − d)φ(k − d) 2

X Z

W (k) = −φ T (k)θˆ (k) −

otherwise.

For simplicity, define

m+d−1

In this section, the control input is determined so that the sliding mode exists along the sliding surface s(k) = 0. By referring to the right-hand side of (24), define

− ··· −

η( ˆ k − 1) |e(k)|

+

3.3. The adaptive sliding mode control design

2 P∞

ˆ ˆ k=υ θ (k) − θ (k − υ) < ∞, 22 P∞  ˆ ˆ β( k ) − β( k − υ) < ∞, k=υ
2 P∞ ρ( ˆ k) − ρ( ˆ k − υ) < ∞,
2 Pk∞=υ R R ˆ p ( r , k) − pˆ i (r , k − υ) dr < ∞ i k=υ 0 for any positive finite integer υ .

Proof. The lemma can be similarly proved by referring to Chen (2006) and Goodwin and Sin (1984). 

R

Z 0

Since pˆ 0 (r , k) ≥ 0, for vmin ≤ v ≤ vmax , it yields W sat (k) ≤

R

Z

pˆ 0 (r , k)Er [v](k)dr ≤ W sat (k).

(49)

0

Discretize [0, R] uniformly into L sub-intervals and define ∆r = R/L, where L is an even number. V ∗ (k) should be derived as follows. If W (k) < W sat (k), V ∗ (k) = vmin ; If W (k) > W sat (k), V ∗ (k) = vmax ; Otherwise, V ∗ (k) is derived as follows.

Step 1: V (0) (k) := V ∗ (k − 1), l := 0. Step 2: x(l) (k) :=

RR pˆ (r , k)Er [V (l) ](k)dr. 0 0 (l) If x (k) − W (k) ≤ ε , go to Step 6; Else if x(l) (k) < W (k) − ε , let V (l+1) (k) := V (l) (k) + ∆r and l := l + 1, then go to Step 3. Else (i.e. x(l) (k) > W (k) + ε ), let V (l+1) (k) := V (l) (k) − ∆r and l := l + 1, then go to Step 4.

Step 3: x(l) (k) :=

RR pˆ (r , k)Er [V (l) ](k)dr. 0 0 (l) If x (k) − W (k) ≤ ε , go to Step 6; Else if x(l) (k) < W (k) − ε , let V (l+1) (k) := V (l) (k) + ∆r and l := l + 1, then go to Step 3. Else (i.e. x(l) (k) > W (k) + ε ), let V (l) (k) := V (l−1) (k) and (l)

V (k) := V (l) (k), then go to Step 5.

Step 4: x(l) (k) :=

RR pˆ (r , k)Er [V (l) ](k)dr. 0 0 (l) If x (k) − W (k) ≤ ε , go to Step 6; Else if x(l) (k) > W (k) + ε , let V (l+1) (k) := V (l) (k) − ∆r and l := l + 1, then go to Step 4. Else (i.e. x(l) (k) < W (k) − ε ), let V (l) (k) := V (l) (k) and (l)

V (k) := V (l−1) (k), then go to Step 5.

X. Chen et al. / Automatica 45 (2009) 469–476

Step 5: x(l) (k) := (l)

x¯ (k) :=

RR 0

pˆ 0 (r , k)Er [V (l) ](k)dr,

R

Z

where $ (k) is defined by

  $ (k) = φ T (k) θˆ (k + d − 1) − θˆ (k) Z R
pˆ 0 (r , k + d − 1) − pˆ 0 (r , k) Er [v](k)dr +

(l)

pˆ 0 (r , k)Er [V ](k)dr , 0

 W (k) − x(l) (k)  (l) (l) . V (l+1) (k) := V (l) + (l) V − V x¯ (k) − x(l) (k) (l)

RR

0

(l−1)

(l)

V

(k) and

V (k) := V (l) (k), then go to Step 5. Step 6: V ∗ (k) := V (l) (k) and stop.

Remark 4. The algorithm is started from the value V (k − 1) obtained in the previous sampling time. This value is adjusted by increasing or decreasing with ∆r until it falls in a certain range. Then, by applying the algorithm in Step 5 (which is an interpolation method), the value of V ∗ (k) is derived. Lemma 2. V ∗ (k) can be found by finite steps of operations, i.e. l is finite when the operations stop.

In this paper, the adaptive sliding mode control input is considered as (50)

3.4. Stability analysis In the following, it is supposed that (51)

Let

ρ(k) =

R

Z

pˆ 0 (r , k)Er [v](k)dr − W (k).

(52)

0

Then,

|ρ(k)| ≤ ε.

(53)

From (45) and (52), it has

φ (k)θˆ (k) +

pˆ 0 (r , k)Er [v](k)dr 0

R

Z

pˆ 1 (r , k)Er [v](k − 1)dr

+ 0

Z

R

pˆ m+d−1 (r , k)Er [v](k − m − d + 1)dr

+ ··· + 0

= C (q−1 )yd (k + d) + δ s(k) + ρ(k) − e(k).

(54)

φ T (k)θ +

R 0 R

Z



"

#



0

" # + 1

C (q−1 )yd (k + d) − ω(k + d) .




(58)

Theorem 1. Consider system (1) satisfying Assumptions 1–4 controlled by the adaptive controller (50). If there exists an integer K > 0 such that W sat (k) ≤ W (k) ≤ W sat (k) for all k > K , then there is a constant γ ∗ > 0 such that, for 0 < γ < γ ∗ , the existence of a constant α ∗ > 0 is guaranteed such that, for the class of disturbances described in (16)–(17) satisfying 0 < α < α ∗ , (i) limk→∞ λ(k) |e(k)| = 0, (ii) all the signals in the closed-loop remain bounded. Proof. See the Appendix.



ˆ 0) and ρ( Remark 5. The initial conditions θˆ (0), pˆ i (r , 0), β( ˆ 0) in (36)–(40) should be chosen such that pˆ i (r , 0) ≥ 0, 0 < RR r pˆ i (r , 0)dr < ∞, and W sat (0) ≤ W (0) ≤ W sat (0). In almost 0 all cases, the parameter adaptation gain χ (k) (for k ≥ 1) in (36)–(40) can be chosen in the range κ ≤ χ (k) ≤ 2 − κ with 0 < κ < 1 such that W sat (k) ≤ W (k) ≤ W sat (k).

Remark 7. In the implementation of the proposed algorithm, by discretizing [0, R] uniformly into L sub-intervals, the integrals can be calculated by the well-known ‘‘Simpson Method’’. Thus, a good calculation of the integrals can be obtained even though L is not very large. As a matter of fact, instead of (1), it is equivalent to considering the control problem of the system

σ¯ (2) (k) = σ (2) (k) + q−d B(q−1 )(w(k) − w( ¯ k)).

gm+d−1 p(r )Er [v](k − m − d + 1)dr

(2)

0

= C (q−1 )yd (k + d) + δ s(k) + e(k + d) + ρ(k) − e(k) − ω(k + d) + $ (k),

0

y(k) −δ  w(k) s(k) 1 − δ q−d 0

0 1 0 σ (k) +  1  (e(k + d) − e(k) + ρ(k) + $ (k)) 0 q−d

" # =

−q−d B(q−1 ) G(q−1 )

(59)

where w( ¯ k) is the calculated value of w(k) defined in (9) using the Simpson Method, σ¯ (2) (k) is defined as

g0 p(r )Er [v](k)dr

+ ··· +

A(q−1 ) F (q−1 ) 0

A(q−1 )y(k) = q−d B(q−1 )w( ¯ k) + σ (1) (k) + σ¯ (2) (k),

By equations (32) and (54), we have

Z

(57)

Remark 6. The parameter δ (0 < δ < 1) in (45) is introduced in order to modify the system response and to control the size of the control signal (Chen, 2006).

R

Z

T

s(k + d) = δ s(k) + e(k + d) + ρ(k) − e(k) + $ (k).

0



W sat (k) ≤ W (k) ≤ W sat (k).

(56)

By substituting (55) into (22) and using (10), the variable s(k + d) becomes

 ∗

v(k) = V ∗ (k).

× Er [v](k − m − d + 1)dr .

The closed-loop system equation can be obtained by combining (1), (55) and (57)

Remark 3. V ∗ (−1) can be chosen as V ∗ (−1) = vmin .

Proof. The proof is omitted.




0

Else if x(l) (k) < W (k) − ε , let V (l) (k) := V (l) (k) and (l−1)

pˆ m+d−1 (r , k + d − 1) − pˆ m+d−1 (r , k)

+ ··· +

(l)

V (k) := V (k), then go to Step 5; Else (i.e. x(l) (k) > W (k) + ε ), let V (l) (k) :=

R

Z

Let l := l + 1 and x (k) := 0 pˆ 0 (r , k)Er [V ](k)dr. If x(l) (k) − W (k) ≤ ε , go to Step 6; (l)

473

(55)

(60)

The boundedness of σ¯ (k) still holds. By replacing the integrals in the formulation with their corresponding calculated values, the result in Theorem 1 is still valid.

474

X. Chen et al. / Automatica 45 (2009) 469–476

4. Simulation results Consider the system described by y(k) + a1 y(k − 1) + a2 y(k − 2)

= b0 w(k − 1) + b1 w(k − 2) + σ (k), Z R p(r )Er [v](k)dr w(k) = 0

with a1 = −1.3, a2 = 0.42, b0 = 1, b1 = −0.65, p(r ) = e−0.067(r −1) , the initial value of w−1 is 0, and 2

h

σ (k) = 0.04 (cos(0.02k))  Z × y(k − 1) + y(k − 2) + 2

ˆ k) = fˆ0 (k), fˆ1 (k) Fig. 1. The estimate θ( R

p(r )Er [v](k − 1)dr

iT

.



0

+ 0.1(sin(0.02k)) qR   R 2 E [v]( k − 1 )) dr kφ(k − 1)k2 + ( r 0  qR × 1 + R 2 E [v]( k − 1 )) dr 1 + kφ(k − 1)k2 + ( r 0 T

with φ(k − 1) = yT (k − 1), yT (k − 2) .



The polynomial C (q−1 ) is chosen as C (q−1 ) = 1 + q−1 + 0.25q−2 . Since d = 1, it is easy to see Ω (q−1 ) = 1. The unknown polynomials F (q−1 ) and G(q−1 ) are respectively F (q−1 ) = 2.3 − 0.17q−1 and G(q−1 ) = 1 − 0.65q−1 .

Fig. 2. The estimates of pˆ 0 (r , k) and pˆ 1 (r , k) with r = 0.1.

φ(0) is set to φ(0) = [0.15, 0]T . θˆ (0) is set to θˆ (0) = [0, 0]T . The initial conditions of pˆ 0 (r , 0) and pˆ 1 (r , 0) are all set to 0.1. The initial value of w−1 is simply assumed as 0. The parameter γ in (33) is chosen as γ = 0.1. The adaptation gain χ (k) is set to χ (k) = 0.8. The weighting factor δ in W (k) is chosen as δ = 0.1. In the simulation, R is chosen as R = 60, L is chosen as L = 6000. The admissible error ε is set to ε = 0.005. The desired trajectory kπ of the output is yd (k) = sin( 100 ). Fig. 1 shows the estimate θˆ (k)

=

h

fˆ0 (k), fˆ1 (k)

iT

of the

parameter θ = [f0 , f1 ] = [2.3, −0.17] . Fig. 2 shows the estimates of pˆ 0 (r , k) and pˆ 1 (r , k) with r = 0.1. It can be seen that the estimates converge very fast. The convergence of pˆ 0 (r , k) and pˆ 1 (r , k) with different r is also confirmed. Fig. 3 shows the control input described in (50), Fig. 4 shows the output tracking error. It can be seen that a very good result is obtained in a relatively short time. In the simulation, Step 3–Step 5 in the algorithm of deriving V ∗ (k) only need to be run for at most one time for a definite k, and the relation W sat (k) ≤ W (k) ≤ W sat (k) holds for all k > 0. T

T

Remark 8. If R and L are chosen as R = 20 and L = 2000, the steady state response of the controlled system does not differ much from that in the case R = 60 and L = 6000; however, the calculation speed becomes much faster. Finally, if the amplitude of the disturbance σ (k) becomes smaller, then, by using the same design parameters, the estimated parameters fˆi (k) and pˆ i (r , k) (for i = 0, 1) converge much faster and the output tracking error can be controlled to be much smaller. Particularly, if the disturbance σ (k) is absent, then the estimated parameters converge even faster and the output tracking error can be controlled to be very small which is almost zero by using the same design parameters.

Fig. 3. The control input.

Fig. 4. The output tracking error.

5. Conclusions This paper discusses the design of an adaptive sliding mode controller for a general class of uncertain discrete time linear systems preceded by hysteresis and disturbance, where the hysteresis is represented by the Prandtl–Ishlinskii model, and the disturbance is composed of model uncertainties, nonlinearities,

X. Chen et al. / Automatica 45 (2009) 469–476

external disturbances, etc. The proposed method fuses the hysteresis output control method with the adaptive control theory for linear discrete time dynamical systems with uncertainties. Only the parameters directly needed in the formulation of the sliding mode controller are updated by adaptive algorithms with dead zone. The robust sliding mode tracking controller is synthesized by using the estimated parameters. The stability of the controlled system is guaranteed in the sense that all signals remain bounded. The proposed method is well suited for controlling the systems actuated by smart-material-based actuators such as piezo-electric actuators and magnetostrictive actuators, where complex hysteresis models are required to describe the hysteretic nonlinearity.

475

By applying (33) and property (P1) of Lemma 1 to relation (67), it can be seen that there exist positive constants C9 and C10 such that

kφ(k)k2 ≤ C9 + C8 max {λ(τ + d) |e(τ + d)|} 1≤τ ≤k

+ C10 γ max kφ(τ )k2 .

(68)

1≤τ ≤k

Then, it can be seen that there exist positive constants C11 , C12 and γ ∗ = 1/C10 such that, if 0 < γ < γ ∗ , then the existence of a constant α ∗ > 0 is guaranteed such that

kφ(k)k2 ≤ C11 + C12 max {λ(τ + d) |e(τ + d)|}

(69)

1≤τ ≤k

The work for C.Y. Su was partially supported by the 111 Project of China (B08015).

for the class of disturbances described in (16)–(17) satisfying 0 < ∗ α< R Rα . By the definition of D(k − d) in (35) and the boundedness of 0 pˆ j (r , i)Er [v](k)dr, it can be proved that there exist positive constants C13 , C14 and γ ∗ such that, if 0 < γ < γ ∗ , then the existence of a constant α ∗ > 0 is guaranteed such that

Appendix. Proof of Theorem 1

p

Acknowledgement

D(k − d) ≤ C13 + C14 max {λ(τ ) |e(τ )|}

(70)

1≤τ ≤k

Since A(q−1 )  det F (q−1 ) 0

−q−d B(q−1 ) G(q−1 )



0



0 0 1 − δq

−d



= (1 − δ q )B(q )C (q ) −d

−1

−1

(61)

which is a Schur polynomial, from (58), there exist positive constants C1 , C2 and C3 such that

kφ(k)k2 ≤ C1 + C2





max |e(τ + d)| + max |$ (τ )|

1≤τ ≤k

1≤τ ≤k

+ α · C3 max kφ(τ )k2 ,

(62)

1≤τ ≤k

RR

where the boundedness of 0 pˆ j (r , i)Er [v](k)dr for arbitrary i, j and k is employed. So, there exists a positive constant α ∗ = 1/C3 such that C1 max kφ(τ )k2 ≤ 1≤τ ≤k 1 − α C3

+



C2 1 − α C3



max |e(τ + d)| + max |$ (τ )|

1≤τ ≤k

(63)

1≤τ ≤k

for the class of disturbances described in (16)–(17) satisfying 0 < α < α ∗ . Therefore, by substituting (63) into (62), it can be concluded that there exist positive constants C4 and C5 such that

kφ(k)k2 ≤ C4 + C5





max |e(τ + d)| + max |$ (τ )|

1≤τ ≤k

(64)

1≤τ ≤k

for the class of disturbances (16)–(17) satisfying 0 ≤ α < α ∗ . By the definition of $ (k) in (56) and Lemma 1, it can be easily shown that there exists a positive constant C6 such that C5 max |$ (τ )| ≤ C6 + 1≤τ ≤k

1

max kφ(τ )k2 . 2 1≤τ ≤k

(65)

By substituting (65) into (64), it can be proved that there exist positive constants C7 and C8 such that

kφ(k)k2 ≤ C7 + C8 max |e(τ + d)| .

(66)

1≤τ ≤k

From (66), it also gives

kφ(k)k2 ≤ C7 + C8 max

|e(τ + d)| − η(τ ˆ + d − 1) 1≤τ ≤k + η(τ ˆ + d − 1) ≤ C7 + C8 max {λ(τ + d) |e(τ + d)|} 




1≤τ ≤k

+ C8 max η(τ ˆ + d − 1). 1≤τ ≤k

(67)

for the class of disturbances described in (16)–(17) satisfying 0 < α < α∗ . From (70) and property (P3) of Lemma 1, result (i) can be proved by employing the famous ‘‘Key Technical Lemma’’ in Goodwin and Sin (1984). Furthermore, from (69), it can be concluded that all the signals in the loop are uniformly bounded. Therefore, the theorem is proved. References Adly, A. A., Mayergoyz, I. D., & Bergqvist, A. (1991). Preisach modeling of magnetostrictive hysteresis. Journal of Applied Physics, 69(8), 5777–5779. Astrom, K. J., & Wittenmark, B. (1989). Adaptive control. Reading, MA: AddisonWesley. Banks, H. T., & Smith, R. C. (2000). Hysteresis modeling in smart material systems. International Journal of Applied Mechanics and Engineering, 5, 31–45. Brokate, M., & Sprekels, J. (1996). Hysteresis and phase transitions. New York: Springer-Verlag. Chen, X. (2006). Adaptive sliding mode control for discrete-time multi-input multioutput systems. Automatica, 42(3), 427–435. Chen, X, Su, C. Y., & Fukuda, T. (2008). Adaptive control for the systems preceded by hysteresis. IEEE Transactions on Automatic Control, 53(4), 1019–1025. Croft, D., Shed, G., & Devasia, S. (2001). Creep, hysteresis, and vibration compensation for piezoactuators: Atomic force microscopy application. ASME Journal of Dynamic Systems, Measurement, and Control, 123, 35–43. Cross, R., Krasnosel’skii, M. A., & Pokrovskii, A. V. (2001). A time-dependent Preisach model. Physica B, 306, 206–210. Goodwin, G. C., & Sin, K. S. (1984). Adaptive filtering, prediction and control. Englewood Cliffs, NJ: Prentice-Hall, Inc. Guyer, R. A., McCall, K. R., & Boitnott, G. N. (1994). Hysteresis, discrete memory and nonlinear propagation in rock. Physical Review Letters, 74, 3491–3494. Huhnen, K., & Janocha, H. (1999). Adaptive inverse control of piezoelectric actuators with hysteresis operators. In Proceedings of the European control conference (paper F0291). Jiles, D. C., & Atherton, D. L. (1986). Theory of ferromagnetic hysteresis. Journal of Magnetism and Magnetic Materials, 61, 48–60. Krasnosel’skii, M. A., & Pokrovskii, A. V. (1989). Systems with hysteresis. New York: Springer-Verlag. Krejci, P., & Huhnen, K. (2001). Inverse control of systems with hysteresis and creep. Proceedings of the Institution of Electrical Engineers Control Theory and Applications, 148, 185–192. Moheimani, S. O. R., & Goodwin, G. C. (2001). Guest editorial introduction to the special issue on dynamics and control of smart structures. IEEE Transactions on Control Systems Technology, 9, 3–4. Mayergoyz, I. D. (1991). Mathematical models of hysteresis. New York: SpringerVerlag. Natale, C., Velardi, F., & Visone, C. (2001). Identification and compensation of preisach hysteresis models for magnetostrictive actuators. Physica B, 306, 161–165. Su, C. Y, Stepanenko, Y., Svoboda, J., & Leung, T. P. (2000). Robust adaptive control of a class of nonlinear systems with unknown backlash-like hysteresis. IEEE Transactions on Automatic Control, 45(12), 2427–2432. Su, C. Y., Wang, Q., Chen, X., & Rakheja, S. (2005). Adaptive variable structure control of a class of nonlinear systems with unknown Prandtl–Ishlinskii hysteresis. IEEE Transactions on Automatic Control, 50(12), 2069–2074. Tao, G., & Kokotovic, P. V. (1995). Adaptive control of plants with unknown hysteresis. IEEE Transactions on Automatic Control, 40, 200–212.

476

X. Chen et al. / Automatica 45 (2009) 469–476

Tao, G., & Lewis, F. L. (Eds.) (2001). Adaptive control of nonsmooth dynamic systems. New York: Springer-Verlag. Tan, X., & Baras, J. S. (2004). Modeling and control of hysteresis in magnetostrictive actuators. Automatica, 40(9), 1469–1480. Visintin, A. (1994). Differential models of hysteresis. New York: Springer-Verlag. Webb, G. V., Lagoudas, D. C., & Kurdila, A. J. (1998). Hysteresis modeling of SMA actuators for control applications. Journal of Intelligent Material Systems and Structures, 9(6), 432–448. Zhou, J., Wen, C., & Zhang, Y. (2004). Adaptive backstepping control of a class of uncertain nonlinear systems with unknown backlash-like hysteresis. IEEE Transactions on Automatic Control, 49(10), 1751–1757.

Takeshi Hisayama received his B.E. and M.E. degrees in Electronics from Shibaura Institute of Technology, Japan, in 2006 and 2008, respectively. He is now serving in the Control Systems Group, IHI Scube Co., Ltd. His research interests include robust nonlinear control, adaptive control and motion control.

Xinkai Chen received his B.S. and M.S. degrees in mathematics from Hebei University, China, in 1986 and 1989, respectively, M.S. degree in electrical engineering from Mie University, Japan, in 1995, and Ph.D. degree in engineering from Nagoya University, Japan, in 1999. He is currently an Associate Professor in the Department of Electronic and Information Systems, Shibaura Institute of Technology, Japan. His research interests include adaptive control, robust nonlinear control, motion control and machine vision. Now, he is serving as an Associate Editor of IEEE Transactions on Automatic Control.

Chun-Yi Su received his Ph.D. from South China University of Technology in 1990. He is the Professor and Concordia Research Chair at Concordia University (Canada). His current research interest is in control of systems involving hysteresis nonlinearities. Dr. Su has served on the editorial boards of several journals, including IEEE Transactions on Automatic Control and IEEE Transactions on Control Systems Technology. He has also served as an organizing committee member for many conferences.