Causal end-effector inversion of a flexible link manipulator

Mechatronics 19 (2009) 1197–1210

Contents lists available at ScienceDirect

Mechatronics journal homepage: www.elsevier.com/locate/mechatronics

Technical note

Causal end-effector inversion of a flexible link manipulator M. Vakil, R. Fotouhi *, P.N. Nikiforuk Department of Mechanical Engineering, University of Saskatchewan, 57 Campus Drive, Saskatoon, Saskatchewan, Canada S7N 5A9

a r t i c l e

i n f o

Article history: Received 4 April 2007 Accepted 22 March 2009

Keywords: Flexible link manipulator Output redefinition Non-minimum systems

a b s t r a c t A new method to obtain the required base torque for the rest-to-rest manoeuvre of a single flexible link manipulator as its end-effector moves along a desired trajectory is introduced. Contrary to available causal and non-causal end-effector inversion techniques, this new approach does not require pre- and/or post-actuation and works even in the presence of the purely imaginary zeros for the transfer function. In this approach, the desired end-effector trajectory is divided into a finite number of segments. In each segment, but the last one, the desired trajectory is redefined so that a bounded continuous torque can be obtained using causal dynamic inversion. For the last segment the desired trajectory is redefined so that not only the causal inversion is achieved but also the final conditions are satisfied; that is, the end-effector reaches its desired position and the manipulator comes to rest at a given time. The redefinition of the desired trajectory at each segment employs summation of the stable exponential functions, which leads to a family of possible answers for the redefined trajectory. To obtain the best member of this family which minimizes the difference between the desired trajectory and redefined trajectory, an efficient search algorithm is proposed and used. Simulation results show the effectiveness of this new method. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction Light flexible link manipulators are potential substitutes for heavy rigid link robot arms provided their performance becomes more reliable and predictable. Since in many industrial applications the time history of the end-effector movement is given beforehand, end-effector trajectory tracking {EETT} of flexible link manipulators {FLM} is of great importance. A possible approach for the EETT of FLM is to use the output regulation method [1]. The feasibility of applying this method to FLM was studied in [2]. However, for the EETT of FLM which have nonminimum phase characteristics [3], applying the method of [1] leads to transient errors at the beginning and end of the manoeuvres [4]. Another alternative approach for the EETT of FLM is to use stabilizing feedback command when the feedforward command is created by the inversion of the system dynamic [4]. Due to the non-minimum phase behaviour of the FLM, causal1 inversion of the dynamic equation for a desired end-effector trajectory can not be achieved and the non-causal inversion technique has to be adopted [4–6]. It is to be noted that the non-casual end-effector inversion method leads to pre- and post-actuations and the method does not work in the presence * Corresponding author. Tel.: +1 306 966 5453; fax: +1 306 966 5427. E-mail addresses: [email protected] (M. Vakil), [email protected] (R. Fotouhi), [email protected] (P.N. Nikiforuk). 1 Causal inversion leads to a causal signal. The causal signal at any time depends on the values of the states up to that time, while the non-causal signal depends on the values of states before and after that time. 0957-4158/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechatronics.2009.03.010

of purely imaginary zeros [7]. Rather than the non-causal input torque a method was proposed in [8] that creates a causal endeffector inversion of a single flexible link manipulator {SFLM} for the point-to-point and rest-to-rest motion. In general, however, it is more desirable to be able to inverse the dynamic equations for the desired end-effector trajectory rather than a rest-to rest and point-to-point manoeuvre. In this paper a causal end-effector trajectory inversion for a SFLM by output redefinition is presented. This method determines the required base torque for the rest-to-rest manoeuvre of a SFLM while its end-effector moves along a desired trajectory. Contrary to the available non-causal end-effector inversion approach [4,5] and the recent approximate causal inversion technique [9], the technique introduced in this paper does not lead to pre- and/or postactuation(s) and works even in the presence of the purely imaginary zeros. For the output redefinition, the summation of the stable exponential functions {SSEF} [10] is employed, as follows:

~d ðtÞ ¼ y

r X

cj emj t ;

mj < 0

ð1Þ

j¼0

~d is the redefinition of the desired end-effector trajectory, yd, where y and the coefficients cj are calculated after selecting r and mj so that a bounded torque with no pre- and post-actuation is obtained; the details of the method are given in Section 4. Comparing a SSEF with the following polynomial yp having order r, which is the same as the number of terms in the SSEF:

1198

yp ðtÞ ¼

M. Vakil et al. / Mechatronics 19 (2009) 1197–1210 r X

pj t j

ð2Þ

j¼0

it is clear that the number of the degree of selection in the exponential functions, that is cj and mj in Eq. (1), is twice the number in the polynomial, that is pj in Eq. (2). Therefore, for a given set of conditions which is equal to the number of the coefficients pj in Eq. (2), there is only one solution for yp, while for the same set of conditions ~d in Eq. (1). Hence, the there is a family of possible solutions for y member of the family of the possible solutions, whose deviation from yd is the smallest, can be obtained. This member is found in Section 6 by using a new search algorithm which is easy to use and enables the near optimum member to be found using a zero-order optimization approach; that is, this new algorithm does not require any derivative calculation. Simulation studies conducted to show the effectiveness of this new inversion method as described in Section 7. In these studies, the required torque for the rest-to-rest manoeuvre, while the end-effector moves along a desired trajectory, was first obtained by using the linear dynamic model of a single flexible link manipulator. The linear dynamic model was derived by using the assumed mode shape method, which is explained in Section 2. The torque obtained based on this linear dynamic model, was then applied to the full-nonlinear finite element model of the single flexible link manipulator. Details of the finite element model used here can be found in [11]. The full-nonlinear finite element model was used to check the accuracy of the dynamic modelling procedure discussed in Section 2 and also to verify the effectives of the new method introduced in this paper. Although, this method was examined for the end-effector trajectory inversion of a SFLM, it can easily be extended to the output inversion of linear single-input single-output non-minimum phase systems with purely imaginary zeros, hyperbolic internal dynamic,2 or non-minimum phase systems without purely imaginary zeros, non-hyperbolic internal dynamic. 2. Dynamic modeling of a SFLM To change the infinite dimensional dynamic equations of a FLM [12] to a finite dimensional one, in this paper a combination of the Lagrange method and the assumed mode shape method {AMM} approximation [13] is employed. In the AMM, the spatial deflection of the flexible link, n in Fig. 1, is described by a finite series composed of spatial pre-defined shape functions /j(z), which are adopted here from [14], multiplied by time varying weight functions kj ðtÞ that is:

nðz; tÞ ¼

n X

/j ðzÞkj ðtÞ

ð3Þ

j¼1

where n is the number of the pre-defined shape functions. By employing the Lagrange equation for a SFLM; that is:

  d @T SM @T SM @U SM þ ¼ Qj  dt @ q_ j @qj @qj

ð4Þ

where TSM is the kinetic energy of the SFLM, USM is the potential (strain) energy of the SFLM, qj is the jth element of the generalized coordinates vector q ¼ ½ h k T ¼ ½ h k1 . . . kn T , and Qj is the generalized force corresponding to qj, the linear dynamic equation can be written as:

€ þ Kq ¼ Hs Mq

ð5Þ

2 The part of the system dynamics which has been rendered ‘‘unobservable” in the input-output linearization is called the internal dynamic. See [15 Chap. 6] for more details.

Fig. 1. Schematic of a SFLM.

Details of the elements of the mass and stiffness matrices, M and K, and the vector H in Eq. (5) are given in the Appendix A. This dynamic model for a SFLM, Eq. (5), with the end-effector displacement as the output can be written also as:

X_ ¼ AX þ Bs

ð6aÞ

y ¼ CX

ð6bÞ

where

X¼

  q ; q_

A¼



C ¼ ½ D 01nþ1 ;

0nþ1nþ1 M 1 K

Inþ1nþ1 0nþ1nþ1

 ;

B¼



0nþ1 1

M H

D ¼ ½ L /1 ðLÞ . . . /n ðLÞ 

 ;

ð7Þ

and y is shown in Fig. 1. 3. Output inversion and internal dynamics Assumption 1. It is assumed that the desired end-effector trajectory, yd, is at least three times differentiable and satisfies the following conditions: (1) yd = 0 for t < 0. €d , the desired end-effector acceleration, is bounded. (2) y (3) yd(t) = yd(tf) for t P tf, where tf is the final manoeuvre time and yd(tf) is the desired final position. Consequently €d ðtÞ ¼ 0 for t > tf . y_ d ðtÞ ¼ y Definition 1. The causal end-effector inversion of a SFLM, Eqs. (6a) and (6b) with A, B and C defined in Eq. (7), for a desired trajectory yd satisfying Assumption 1 means that the bounded causal pair (s, X) has to be found so that: 

_ XðtÞ ¼ AXðtÞ þ BsðtÞ for 0  t  t f yd ðtÞ ¼ CXðtÞ (2) s(t) = 0, X(t) = 0 for t < 0 _ ¼ 0, and y(t) = yd(tf) for t > tf where (3) sðtÞ ¼ 0; kðtÞ ¼ kðtÞ   _ ¼ k_1 . . . k_ n : kðtÞ ¼ ½ k1 . . . kn  and kðtÞ (1)

Definition 2. Consider the linear single-input single-output system of the form given in Eqs. (6a) and (6b). This system has a well-defined relative degree rw if and only if CAðrw 1Þ B – 0. In other words, in the rwth time differentiation of Eq. (6b), the input s appears explicitly. r r That is yðrw Þ ¼ CArw X þ CArw 1 Bs; where yðrw Þ ¼ d w y=dt w . It is to be noted that in the absence of a modeling error, initial condition error and external disturbance, applying the offline computed inverse torque s as the input to the system given in Eq. (6a) results in the desired end-effector movement. Even in the presence of the above errors and disturbance, stabilization of the system around the nominal states can be done for instance using the online state-feedback technique [16]. Thus the main challenge of the end-effector regulation for flexible link manipulators with

1199

M. Vakil et al. / Mechatronics 19 (2009) 1197–1210

the inverse dynamic command as the feedforward signal, which are a non-minimum phase system, is to find the nominal state, X, and input torque, s [4,7]. For a SFLM rw = 2 as noted in [5,6], which means CAB – 0 where C, A and B are given in Eq. (7). Thus the output, that is the endeffector displacement y, has to be differentiated twice so that the base torque, s, appears explicitly; thus:

€ ¼ CA2 X þ CABs y

ð8Þ

inversion based on the assumption that there are no purely imaginary eigenvalues for AI, non-causal inversion [4–6] and approximate causal inversion [9] can be adopted which lead to pre- and/or postactuation(s). However, in this paper, to have a bounded solution for XI, yd is redefined so that causal integration of Eq. (16) provides a bounded XI and the corresponding torque does not have post-actuations. Moreover, it is to be noted that the existence of the purely imaginary eigenvalues for AI is not a restriction for the implementation of the new method introduced here.

Assuming the base torque:

s¼

4. Causal inversion by output redefinition

1 €d  CA2 XÞ ðy CAB

ð9Þ

where yd is the desired end-effector trajectory and substituting Eq. (9) in Eq. (8) leads to:

€¼y €d y

ð10Þ

Since the order of the inversed system in Eq. (10) is two while the order of the original system given in Eq. (5) is 2(n + 1) there is an internal dynamic [15] of the order 2(n + 1)  2. To obtain the associated internal dynamic, the transformation:

w ¼ TX

ð11Þ

is considered where:



T

w ¼ y y_ k1 . . . kn k_ 1 . . . k_ n ; 2 3 D1nþ1 01nþ1  Inn 6 D1nþ1 7 T ¼ 4 01nþ1 5; Io ¼ 0nnþ1 02n1 ðIo Þ2n2nþ1

0nnþ1



Inn

The transformation in Eq. (11), changes the dynamic model given in Eq. (6a) into:

w_ ¼ TAT 1 w þ TBs

ð12Þ

Substituting the torque given in Eq. (9) into Eq. (12) results in:

€d w_ ¼ Aw w þ Bw y

ð13Þ

where Aw and Bw are:

Aw ¼ T 1 AT 

TBCA2 T 1 ; CAB

Bw ¼

TB CAB

ð14Þ

and as mentioned earlier the scalar CAB – 0. According to the definition of w and using the torque given in Eq. (9), for a SFLM the matrices Aw and Bw are:

2

0

1

6 Aw ¼ 4 0 02n1

0 02n1

3

2

0

3

7 012n 5; ðAI Þ2n2n

6 Bw ¼ 4

1

7 5

012n

ð15Þ

ðBI Þ2n1

umns of the matrix G are the eigenvectors of the matrix AI, and as a result G1AIG is a diagonal matrix, Eq. (16) can be written as:

"

X_ sI X_ u I

G1 BI

G1 AI G

# ¼

zfflfflfflfflfflfflffl}|fflfflfflfflfflfflffl{   " s # AI 0 XI 0

AþI

X uI

þ

zfflffl}|fflffl{ " # BsI BuI

€d y

ð17Þ

where A I is a diagonal matrix with diagonal elements that are the eigenvalues of AI with negative real part, and Aþ I is a diagonal matrix with diagonal elements that are the eigenvalues of AI with positive real part. Also, if AI has purely imaginary eigenvalues, the corresponding Jordan block will be included in Aþ I . Due to the existence of Aþ I the forward integration of Eq. (17) for €d , generally results in an una bounded desired acceleration, y bounded answer for X uI and thus XI will be unbounded. However, in this paper, for a given set of initial conditions on X uI , the value ~d so that the causal solution of: of yd is redefined by y

€~d X_ uI ¼ AþI X uI þ BuI y

ð18Þ

~ ~€d €d is the redefined acceleration.4 After finding y is bounded where y so that X uI from Eq. (18) is bounded, the bounded causal X sI will be found solving:

~€ X_ sI ¼ AI X sI þ BsI y d Having bounded h ~d bounded XI, w ¼ y

ð19Þ X uI

and

~_ d y

X sI

X TI

~ €d , the for the redefined acceleration y

iT

and X (from Eq. (11)) can be calcu-

~ €d and the corresponding X, the required causal lated. Thus having y ~d is obtained from Eq. (9). To clarify the torque for the inversion of y basic concept of calculating bounded X uI by output redefinition the following example is provided. 4.1. Example for the stable inversion by output redefinition

and the internal dynamic is:

€d X_ I ¼ AI X I þ BI y

In this section the basic concept of the causal inversion by output redefinition is clarified and a simple example is provided. h iT where the colUsing the transformation X I ¼ G ðX sI ÞT ðX uI ÞT

ð16Þ



T where X I ¼ k1 . . . kn k_ 1 . . . k_ n . If the causal XI from  T Eq. (16) is bounded for a bounded €d ; X ¼ T 1 yd y_ d X TI from Eq. (11) is also bounded. Thus, the y base torque s from Eq. (9) for the causal inversion of yd is bounded. However, since the SFLM is a non-minimum phase system, some of the eigenvalues of AI in Eq. (16) have positive real parts ([15] Ch. 6.1.3, [6]) and thus the internal dynamic, Eq. (16), is unstable. This means that the feedforward integration of Eq. (16) even for a €d generally leads to an unbounded XI. For the stable bounded3 y

3 The bounded signal does not grow to infinity as time goes to infinity while the unbounded signal will. Also, a signal which has an ever-oscillating behavior is classified here as an unbounded signal.

€d Find the bounded xuI by the output redefinition for x_ uI  axuI ¼ by where xuI is scalar. Without the loss of generality, it is assumed that the following differential equation has to be solved for xuI :

€d x_ uI  axuI ¼ by

ð20Þ

€d is the accelwhere a > 0 and b are arbitrary constant scalars and y eration of the desired trajectory. Also the initial condition on xuI is:

xuI ðti Þ ¼ ðxuI Þi

ð21Þ

4 The velocity and acceleration which correspond to the redefined trajectory are called redefined velocity and acceleration, respectively.

1200

M. Vakil et al. / Mechatronics 19 (2009) 1197–1210

where ti is the initial time and ðxuI Þi is the initial value of xuI . The xuI from Eq. (20) is composed of two parts, the complementary part, ðxuI Þc , and the particular part, ðxuI Þp , that is:

xuI ¼ ðxuI Þc þ ðxuI Þp

ð22Þ ðxuI Þp ,

€d and can be found by the The particular part, depends on y convolution integral. The complementary part, ðxuI Þc , is:

ðxuI Þc ¼ de

at

ð23Þ

where d is constant and is found by imposing the initial condition on xuI . Since a > 0 if d – 0, ðxuI Þc will be unbounded and thus generally an unbounded xuI will exist (See Eq. (22)). To have a bounded ~ €d is replaced by it redefinition, y €d , so that: xuI from Eq. (20), y (1) The constant d in Eq. (23), which is obtained after imposing the initial condition on xuI , given in Eq. (21), becomes zero. ~ €d is bounded. (2) The particular solution ðxuI Þp corresponds to y

Zk

ðC T Þk zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl3{ ck 2 2 m1 tf m t m t m t zfflffl}|fflffl{ bm1 e bm22 e 2 f bm23 e 3 f bm24 e 4 f 2 3 zfflfflfflfflfflffl}|fflfflfflfflfflffl{ 2 u 3 ðx c 1 m2 a m3 a m4 a I Þf 7 6 m1 a 6 bm2 em1 ti bm2 em2 ti bm2 em3 ti bm2 em4 ti 7 6 c 7 6 ðxu Þ 7 2 3 4 7 6 27 6 I i 7 6 1 m2 a m3 a m4 a 7 6 6 m1 a 7¼6 7 7 4 c3 5 4 yd ð0Þ 5 6 mt m t m t m t 1 2 3 4 i i i i 5 4 e e e e yd ðtf Þ c4 em1 tf em2 tf em3 tf em4 tf 1 ck ¼ Z k ðC T Þk )

ð30Þ

It is to be noted that the addition of the condition given in Eq. (29) necessitates that r in Eq. (24) increases from 3 to 4. In addition to the conditions presented in Eq. (26), to have the same velocity and acceleration for the redefined trajectory at initial time ti and final time tf as the desired ones, the following conditions have also to be considered:

~_ ðti Þ ¼ y_ ðti Þ; y d d ~€ ðti Þ ¼ y € ðti Þ; y d

d

~_ ðt Þ ¼ y_ ðt Þ y d f d f ~€ ðt Þ ¼ y € ðt Þ y d

f

d

f

€d are redefined by: For this purpose, yd and y

~d ¼ y

r X

cj emj t ;

~€ ¼ y d

j¼1

r X

m2j cj emj t ;

mj < 0

ð24Þ

j¼1

~ €d in Eq. (20) is replaced with y €d , giving: Then y

~€ x_ uI  axuI ¼ by d

ð25Þ

~d , to have the same In addition, to force the redefined trajectory, y values as the desired trajectory, yd, at initial time ti and final time tf the following conditions must be imposed:

~d ðti Þ ¼ yd ðt i Þ; y

~d ðt f Þ ¼ yd ðt f Þ y

ð26Þ

The xuI from Eq. (25) is given in Eq. (22) where ðxuI Þc is given in Eq. ~ €d , is: (23). Also ðxuI Þp for the redefined acceleration, y

ðxuI Þp ¼

2 r X bmj cj m t e j mj  a j¼1

ð27Þ

and by selecting mj < 0 and mj – a, ðxuI Þp is bounded. Therefore, to have a bounded solution for xuI ; ðxuI Þc must be zero, which means that d = 0. Thus, assuming that the mj in Eq. (24) are known,5 the unknown cj are chosen to make d = 0 in Eq. (23) and to meet the conditions given in Eqs. (21) and (26). To have a unique solution for cj, the sum of the number of the imposed conditions, Eqs. (21) and (26), has to be the same as the number of the unknowns, cj. Since three conditions are given in Eqs. (21) and (26), r has to be 3 and to force d to be zero, cj(j = 1. . .3) is calculated as: Zk

ck ðC T Þk zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ 2 2mt 3 zfflffl}|fflffl{ 2 3 zfflfflfflfflfflffl}|fflfflfflfflfflffl{ 2 u 3 bm1 e 1 i bm22 em2 ti bm23 em3 ti ðx c 1 I Þi m2 a m3  7 6 m1 a 6 7 6 7 6 m1 ti m2 t i m3 t i 7 4 c 2 5 ¼ 4 yd ð0Þ 5 e e 4 e 5 yd ðt f Þ c3 em1 tf em2 tf em3 tf

4.2. Satisfaction of the final conditions imposed on X sI ~ €d will be found from Eq. (19). X sI which corresponds to y Since all the eigenvalues of A I have a negative real part, the forward integration of Eq. (19) from any initial condition leads to bounded X sI . Thus satisfaction of any initial condition on X sI is always possible. However if it is required that X sI ðtf Þ ¼ fðX sI Þgf , which means that the final value of X sI at tf is fðX sI Þgf , the redefinition of the desired trajectory has to be such that this condition is met. This can be achieved similar to satisfying the final condition on X uI , which is explained in the example provided in Section 4.1 whose required final condition is given in Eq. (29). 5. Piece-wise trajectory inversion by output redefinition For the piece-wise trajectory inversion the following steps have to be taken. (1) Divide the desired end-effector trajectory, yd, into several consecutive equal segments such that:

)

ck ¼ Z 1 k ðC T Þk ð28Þ

After finding cj(j = 1. . .3) for the given mj(j = 1. . .3) from Eq. (28), ~ ~d and y €d are known as given in Eq. (24). Also xuI ¼ ðxuI Þp , as given y in Eq. (27), which is bounded. If in the above example the following condition is required:

xuI ðt f Þ ¼ ðxuI Þf

ð29Þ

where ðxuI Þf is given, the coefficients of the redefined trajectory, cj in Eq. (24), have to be selected so that the settings given in Eqs. (21), (26), and (29) are satisfied and the condition d = 0 in Eq. (23) is met. Thus cj(j = 1. . .4) must be calculated from:

5

~d approach yd one can add the conditions that Moreover to make y ~ €d , up to order h be the derivatives of the redefined acceleration, y equal to the corresponding derivative of the desired acceleration, €d , at the initial time ti and final time tf, if these derivatives of the y desired acceleration exist. This will increase the value of r in Eq. (24).

A method for the selection of mj < 0 is given in Section 6.

ðyd Þk ¼ yd ;

t ik  t  t fk ;

k ¼ 1...v

ð31Þ

where v is the number of the segments and tik and tfk; are the initial and final times of the kth segments. Note that ti1 = ti and tfv = tf are the initial and final manoeuvre times, respectively. (2) In the kth segment, redefine the desired trajectory, (yd)k, by the summation of the stable exponential functions {SSEF}, that is:

~d Þk ¼ ðy

r X

cjk emjk t ;

mjk < 0

ð32Þ

j¼0

where the selection of mjk and r are explained in detail in Section 6. (3) Find the contribution of each stable exponential function, cjk in Eq. (32), for the kth segment so that the following conditions are met:

1201

M. Vakil et al. / Mechatronics 19 (2009) 1197–1210

3.1.

Assuming yd e Cw (w P 3 according to assumption 1 in Section 3) the displacement and its derivatives up to order h where 3 6 h 6 w for the redefined trajectory at the beginning and end of the kth segment must be equal to the desired trajectory; that is: ðiÞ ~ðiÞ y d ðt ik Þ ¼ yd ðt ik Þ

ð33Þ

ðiÞ ~ðiÞ y d ðt fk Þ ¼ yd ðt fk Þ

3.2.

3.3.

where k = 1. . .v, i = 0. . .h and y(h) = dhy/dth. The complementary part of the solution of Eq. (18) must be zero. This condition guarantees that a bounded solution for the internal dynamic exists. The continuity of X uI at the beginning of the kth segment must be satisfied. This condition guarantees that the states X and torque s are continuous. From ~_d and y ~ ~d ; y €d are Eq. (33), it can be seen that y h iT ~_ and continuous. Since X ¼ T 1 y ~d y X TI d  s T  T u T ; to have a continuous X, XI has X I ¼ G ðX I Þ ðX I Þ

to be continuous. This means X sI and X uI must be continuous. Since X sI ; which is calculated by the forward integration of Eq. (19) is continuous, if X uI is also continuous, then X and consequently the corresponding torque s are also continuous. Thus, to have a continuous solution for X uI the final value of X uI at the (k  1)th segment is considered as the initial condition for X uI at the kth segment. 3.4. For the last segment, not only the three conditions of 3.1, 3.2 and 3.3 must be satisfied but also the condition X uI ðt f Þ ¼ X sI ðt f Þ ¼ 0 must be met where tf is the final manoeuvre time. This condition ensures that at the end of the manoeuvre, the deflection, n, and its _ are zero.6 velocity, n, ~ €d and the corresponding bounded X uI ; find (4) After calculating y s X I by feedforward integration from Eq. (19). ~ €d , the bounded Having X uI and X sI for the redefined acceleration, y h iT T _ in Eq. (11), and X can be calculated. Then, the XI, w ¼ y ~d y ~d X I ~d is available from Eq. (9). required causal torque for the inversion of y Theorem. Given the linear dynamic model of a SFLM, Eq. (5) or Eq. (6a) and a desired trajectory yd which satisfies Assumption 1 given in Section 3, the causal continuous torque for the end-effector inversion by the output redefinition with no post-actuation can be found provided conditions (3.1) to (3.4) are satisfied. Proof. See Appendix B. h

6. Selection of the variables mjk, r and

v

In this section the required number of exponential functions, r, that should be used for the output redefinition in the kth segment is explained. Moreover, the selection of the near optimum set of mjk(j = 1. . .r) in the kth segment is discussed. This set minimizes the difference between the redefined trajectory and desired trajectory, as defined in Eq. (37). Finally, the maximum number of segments, v, that can be used for the redefinition of the desired trajectory is detailed.

6

See Eq. (3) for n ¼

Pn

j¼1 kj /j ,

also note X I ¼ ½k

kT and X I ¼ G½ðX sI ÞT

ðX uI ÞT T .

6.1. Required number of exponential functions, r After the selection of mjk for the kth segment, which will be explained in Section 6.2, cjk is found by solving a set of linear algebraic equations (see Eq. (36)). To find a unique cjk from the linear algebraic equations, the number of unknowns, r in P ~d Þk ¼ rj¼0 cjk emjk t , must be equal to the number of equations ðy which resulted from conditions 3.1 to 3.4 for the last segment and conditions 3.1 to 3.3 for the rest of the segments. These conditions are discussed in Section 5. Therefore, the required number of exponential functions for the all segments except the last segment is:

r ¼ 2ðh þ 1Þ þ nu

ð34Þ

Here 2(h + 1) equations come from the condition 3.1 where h is highest derivative of the redefined displacement in Eq. (33), and the nu equations are due to the continuity of X uI at the start of each segment, in which nu is the size of the vector X uI in Eq. (18). For the last segment the required number of the exponential functions is:

r ¼ 2ðh þ 1Þ þ 2nu þ ns

ð35Þ X sI

where ns is the size of the vector in Eq. (19). The addition of nu + ns to Eq. (34) to calculate the required number of exponential functions for the last segment, Eq. (35), is due to the addition of the condition 3.4 in Section 5. 6.2. Selection of the best family member in the kth segment To satisfy the conditions for the piece-wise stable inversion, which are discussed in Section 5, the following equation between the cjk and mjk ðj ¼ 1 . . . rÞ of the kth segment has to be solved (see Eq. (30) for example):

Kðt ik ; tfk ; m1k ; . . . ; mrk ÞC k ¼ Ck

ð36Þ

T where K(tik, tfk, mj1, . . . , mjr) is a square matrix, C k ¼ ½ c1k . . . crk  1 and Ck is a known vector. From Eq. (36), Ck = K Ck and conseP ~d Þk ¼ rj¼0 cjk emjk t is available in quently the redefined trajectory ðy ~d Þk is available; terms of mjk (j = 1. . .r) and tik and tfk. Given that ðy the overall difference between the redefined and desired trajectories for the kth segment is defined as:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Z tfk 2 ~ errork ðt ik ; tfk ; m1k ; . . . ; mrk Þ ¼ ððyd Þk  ðyd Þk Þ dt

ð37Þ

t ik

Since tik and tfk are known, errork in the above equation only depends on mjk (j = 1. . .r). To minimize the error defined in Eq. (37), a method to find the near optimum set of mjk (j = 1. . .r) is introduced in the following. Maximum and minimum values for the mjk in the kth segment, are first selected as (mmax)k and (mmin)k, respectively. Then, to maximize the difference between the consecutive mjk, for the numerical stability of the inverse calculation of matrix K defined in Eq. (36), it is assumed that mjk(j = 1. . .r) are equally spaced between (mmax)k and (mmin)k. Therefore:

mjk ¼ ðmmax Þk þ

ðmmin Þk  ðmmax Þk ðj  1Þ; ðr  1Þ

j ¼ 1...r

ð38Þ

where

ðmmin Þk < ðmmax Þk  0

ð39Þ

By combing Eqs. (38) and (37), errork will be a function of (mmax)k and (mmin)k. Therefore, it is possible to take partial derivatives from errork with respect to (mmax)k and (mmin)k and equal them to zero to find the optimum (mmax)k and (mmin)k. However, deriving the analytical expression for errork, taking partial derivatives from it, and solving the nonlinear algebraic equations, which are obtained by

1202

M. Vakil et al. / Mechatronics 19 (2009) 1197–1210

equalling these derivatives to zero, are not easy. To overcome these difficulties a new search algorithm, which is a zero-order optimization algorithm, is introduced and used in this paper. Based on this algorithm for the kth segment, the 3D graph of errork versus (mmax)k and (mmin)k is plotted. Then from this graph, (mmax)k and (mmin)k which correspond to the minimum point is searched and obtained. That is the plane (mmax)k and (mmin)k is first meshed as shown in Fig. 2. Then at each node the value of errork from Eq. (37) is calculated. Finally (mmax)k and (mmin)k which give the minimum value for errork are selected as the best pair. After obtaining this best pair of (mmax)k and (mmin)k, their corresponding mjk (j = 1. . .r) are found from Eq. (38). As shown in Fig. 2, there is no need to mesh the quadrants one, two and four and also the area above the line b in quadrant three, since these domains are infeasible because of Eq. (39). Moreover, by approaching the line b the difference between the (mmax)k and (mmin)k becomes smaller and so do the differences between mjk (j = 1. . .r). As a result, the matrix K approaches singularity, consequently it is impossible to find C k ¼ ½ c1k . . . crk T from Eq. (36), and thus redefinition is not achievable. This area, where K is singular, is called the singular domain in Fig. 2. To find the feasible domain in quadrant three, the condition rank(K) = r is imposed where r is defined in Section 6.1. To check this condition in the simulation studies presented in Section 7, the MATLAB command rank (K,e) = r is used where e was selected to be in the range 1013 to 1011. The command rank (K,e) returned the number of singular values of matrix K that are bigger than e. 6.3. Number of segments, v By increasing the number of segments, v in Eq. (31), the difference between the redefined and desired trajectories will decrease and the efficiency of the method will improve, at the expense of higher computational cost. This is shown in the simulation section, Section 7. However, as v increases the feasible domains for the segments shrink, and for the last segments this reduction in the feasible domain is the greatest. Following paragraph explain this statement.

Fig. 2. Schematic of the plane (mmax)k and (mmin)k, feasible, infeasible and singular domains Feasible domain: Marked by solid lines, Infeasible domains: Quadrants one, two, four and the area above line b in quadrant three, Singular domains: Area below line b in quadrant three minus feasible domain.

If the desired trajectory is divided into v segments, then the length of each segment is tf/v, where tf is the final manoeuvre time and the initial manoeuvre time is assumed to be zero. Therefore, an increase in v leads to a decrease in the length of each segment. That is, the initial time tik and final time tfk of the kth segment approach each other. If tik and tfk are too close to each other, then the possibility of having matrix K, defined in Eq. (36), to be singular will increase. Consequently the feasible domain will decrease. For example, if tf approaches ti in Eq. (30) (assume tf = tfk and ti = tik), the possibility of having the last two rows of Zk (or even its first two rows) being identical increases, and Zk thus becomes singular. However by increasing v, the matrix K which is associated with the last segment, has the highest possibility of being singular compared to K for the rest of the segments. This is due to a smaller absolute value of the derivative of the exponential function emt for the longer time t. That is, since for the last segment the values of tiv and tfv are the greatest, compared to the initial and final times of the other segments, the difference between emtiv and emtf v of the last segment becomes the smallest. The small value of this difference for the last segment is shown schematically in Fig. 3 where DFS ¼ emti1  emtf 1 of the first segment is compared with DLS ¼ emtiv  emtf v of the last segment. Consequently by increasing v, not only does tiv approach tfv which makes emtiv approach emtf v , but also the difference between emtiv and emtf v of the last segment becomes much smaller compared to the rest of the segments. Thus, as v increases the greatest reduction in the feasible domain is for the last segment. Therefore, although the larger number of segments increases the accuracy of the output redefinition, there is a practical limit for the maximum number of segments to be used. This maximum number of segments corresponds to the smallest possible difference between the initial and final times of the last segment, which for that there is a feasible domain, and can be found as follows: (1) Start and set v = 1. (2) Try to find the feasible domain of the last segment for the number of segments v. When v = 1 there is only one segment which will be referred to as the final segment. (3) If there is a feasible domain for the last segment, then increase the number of segment by one, that is set v = v + 1 and go to step 2. (4) If there is no feasible domain for the last segment, go to Step 5. (5) Stop and the maximum number of segments is v  1.

Fig. 3. Graph of an exponential function, FS: First segment, LS: Last segment DFS ¼ emti1  emtf 1 and DLS ¼ emtiv  emtf v .

M. Vakil et al. / Mechatronics 19 (2009) 1197–1210 Table 1 Physical properties of the SFLM. L (m)

EI (N m2)

Ih (kg m2)

q (kg/m)

Itip (kg m2)

mtip (kg)

1.000

16.00

0.1000

0.2000

0.0500

1.000

7. Simulation results

SFLM and the end-effector manoeuvre from the finite element model was compared with the desired trajectory. For the inversion, the flexibility of the link was modeled using one natural mode shape in which the first nonzero natural frequency of the manipulator was f1 = 3.200 Hz. Moreover the zeros of the transfer function, the eigenvalues of matrix AI in Eq. (16), were ±29.6i. For this manipulator, the transfer function between the end-effector displacement, y, and input torque, s, was y

In this section the results of the simulation studies for a SFLM of the type shown in Fig. 1, and with the physical properties given in Table 1 are presented. Here L is the length of the link, EI is the rigidity of the link, Ih is the mass moment of inertia of the hub, q is the mass per unit length of the link, and Itip and mtip are the mass moment of inertia and the mass of the payload at the end-effector, respectively. The method introduced in this paper, output redefinition, was used for the calculation of the inverse torque by employing the linear dynamic model given in Eq. (5). This calculated inverse torque was then applied to the full-nonlinear finite element model of the

1203

s

2

þ332:5 ¼ 0:378s . It is to be noted that since the zeros were purely s2 ðs2 þ404:5Þ

imaginary, the available techniques for trajectory redefinition [4–6,9] were not applicable [7]. The desired end-effector displacement was selected as:

yd ðtÞ ¼

8 9:259t 5 þ 6:944t4 for 0:0  t < 0:4 > > > > 5 4 3 2 > > < 2:314t  6:944t þ 6:250t  1:250t þ :0973

for 0:4  t < 0:8 > > 5 > > þ 48:61t4  100:0t 3 þ 100:0t2 þ 48:00t þ 9:112 9:259t > > : for 0:8  t  1:2 ð40Þ

Fig. 4. Desired end-effector displacement, velocity and acceleration.

Fig. 5. Example 1, redefining the desired trajectory with one segment, (a) singular, feasible, and infeasible domains, (b) error1 versus (mmax)1 and (mmin)1.

1204

M. Vakil et al. / Mechatronics 19 (2009) 1197–1210

The desired end-effector displacement, velocity and acceleration are shown in Fig. 4. From Eq. (33), to have a continuously differentiable redefined ~ €d , h = 3 was selected. Moreover since there were only acceleration y two purely imaginary zeros for the system, thus nu = 2 and ns = 0. Therefore, the number of the exponential functions for all segments was 10 from Eq. (34); expect the last segment which was 12 from Eq. (35). To find the best member of the family of possible solutions by using the method introduced in Section 6.2, the plane mmax and mmin was meshed with the increment of 0.04. That is the mesh size, as shown in Fig. 2, was 0.04. 7.1. Example 1, output redefinition with one segment;

v=1

For the output redefinition, first the whole desired trajectory was redefined in one segment which meant that v = 1.0 with ti1 = 0.00 s and tf1 = 1.2 s. The required number of the exponential functions for output redefinition was r = 12 from Eq. (35). Although Eq. (35) has to be used for the last segment, when v = 1.0 there is only one segment which has to be treated as the last segments (See requirement 3.4 in Section 5). The best set of (mmax)1 and (mmin)1 was obtained according to the method explained in Section 6.2. The singular, feasible and infeasible domains for this case are

shown in Fig. 5a. The variation of the error1 versus (mmax)1 and (mmin)1 in the feasible domain is also presented in Fig. 5b. From Fig. 5b, the best values for (mmax)1 and (mmin)1 were (mmax)1 = 0.0000 and (mmin)1 = 3.520 and the associated error1 was 0.0342. As explained in Section 6.2, these best values correspond to the point in plane (mmax)1 and (mmin)1, where error1 from Eq. (37) is minimum. Using (mmax)1 = 0.0000 and (mmin)1 = 3.520, the redefined end-effector trajectory, difference between the redefined and desired trajectories as defined in Eq. (41), and the required torque are shown in Fig. 6.

~d ðtÞ  yd ðtÞ Difference ¼ y 7.2. Example 2, output redefinition with two segments,

ð41Þ

v=2

In the second example, the desired output was divided into two equal segments for the redefinition process; that is v = 2.0 with ti1 = 0.00 s, tf1 = 0.60 s, ti2 = 0.60 s and tf2 = 1.2 s. The numbers of the required exponential functions for the first and second segments were 10 and 12 respectively as already discussed. For this case, the feasible and singular domains, as well as the variation of errork versus (mmax)k and (mmin)k in their corresponding feasible

Fig. 6. Example 1, redefining the desired trajectory with one segment, (a) desired end-effector trajectory (dashed) and redefined end-effector trajectory (solid), (b) difference between redefined and desired trajectories; (c) required torque.

M. Vakil et al. / Mechatronics 19 (2009) 1197–1210

domains when k = 1, 2 are shown in Fig. 7. The best member of the family of possible solutions for the first segment, that is ti1 = 0.00 s and tf1 = 0.60 s, from Fig. 7b corresponded to (mmax)1 = 2.280 and (mmin)1 = 20.00 which resulted in the (error)1 = 0.3893  104. These values for the second segment, that is the segment from ti2 = 0.60 s to tf2 = 1.2 s, from Fig. 7d were (mmax)2 = 0.000 and (mmin)2 = 9.000 and the corresponding (error)2 was 2.021  103. By employing these best values for mmax and mmin for the first and the second segments, the redefined end-effector trajectory, the difference between the redefined and desired trajectories as given in Eq. (41), and the required torque are shown in Fig. 8. 7.3. Example 3, output redefinition with three segments, v = 3 Finally, the desired trajectory was divided into three equal segments for the redefinition process; that is v = 3.0, with ti1 = 0.00, tf1 = 0.40, ti2 = 0.40, tf2 = 0.80, ti3 = 0.80 and tf3 = 1.2. For the first two segments the required number of the exponential functions was r = 10 from Eq. (34) while for the last segment it was r = 12 from Eq. (35). The search algorithm in Section 6.2 was used to find the best family member in each segment. The singular and feasible domains for this case are shown in Fig. 9. From Fig. 9b, the best family member in the first segment with ti1 = 0.00 s and

1205

tf1 = 0.40 s, was obtained as (mmax)1 = 3.440 and (mmin)1 = 20.00 which lead to the (error)1 = 1.013  105. From Fig. 9d for the second segment with ti2 = 0.40 s and tf2 = 0.80 s, the best member corresponded to (mmax)2 = 3.680 and (mmin)2 = 14.52 which resulted in the (error)2 = 2.066  106. Finally from Fig. 9f, for the last segment with ti3 = 0.80 s and tf3 = 1.2 s, these values were (mmax)3 = 0.00, (mmin)3 = 11.56 and the (error)3 = 4.490  105. Using these best values of mmax and mmin obtained above, the redefined end-effector trajectory, the difference between the redefined and desired trajectories and the required torque, when v = 3 are shown in Fig. 10. The redefined trajectory and desired trajectory can not be distinguished in Fig. 10a since they are very close to each other. However, there is a very small difference between them as shown in Fig. 10b. Remark 1. By comparing Figs. 6b, 8b and 10b of example 1, 2 and 3, respectively, it was observed that increasing the number of segments, v in Eq. (31), reduced the magnitude of difference between the desired and redefined trajectories. The maximum absolute value of the difference between these two trajectories were 6.320  102 m, 4.900  103 m and 1.355  104 m for v = 1, 2 and 3, respectively. However, by increasing the number of segments, the feasible domains of the higher segment became

Fig. 7. Example 2, redefining the desired trajectory with two segments, (a) singular, feasible, and infeasible domains for the first segment (ti1 = 0.00(s) and tf1 = 0.60(s)), (b) (error)1 versus (mmax)1 and (mmin)1 for the first segment, (c) singular, feasible, and infeasible domains for the second segment (ti2 = 0.60(s) and tf2 = 1.2(s)), (d) (error)2 versus (mmax)2 and (mmin)2 for the second segment.

1206

M. Vakil et al. / Mechatronics 19 (2009) 1197–1210

smaller, as explained in Section 6.3. For v = 2, example 2, the feasible domain of the second segment, Fig. 7c, is much smaller compared to the feasible domain of the first segment, Fig. 7a. This can also be seen for v = 3, example 3, where the feasible domain of the third segment, Fig. 9e, is much smaller than the feasible domain of the second segment, Fig. 9c, while the feasible domain of the second segment is also smaller compared to that of first segment in Fig. 9a. Moreover, from Figs. 6b, 8b and 10b, the difference between the desired and redefined trajectories and its first three time derivatives at the beginning and end of each segment were zero. For example in Fig. 8b, these derivatives were zero at times 0.00, 0.60 and 1.2 s. This is a confirmation of the accuracy of simulation which indicated that the imposed conditions 3.1 in Section 5 were satisfied, as expected. Remark 2. It is to be noted that by using the method discussed in Section 6.3, it was observed that the maximum number of segment was three. Using this maximum number of segments was resulted in very small difference between the redefined and desired trajectories as these two trajectories can not be distinguished in Fig. 10a. Hence, although it was not possible to increase the number of segments beyond three, since the deviation of the redefined trajectory from the desired trajectory for v = 3 was very small, there was not even a need to have more than three segments to increase the accuracy of the redefinition.

7.4. Example 4, applying the torques in examples 1, 2, and 3 to the fullnonlinear finite element model In the section, the torques obtained from the linear dynamic by modeling the flexibility with one mode shape, were applied to the full-nonlinear finite element model (FEM) of the SFLM. That is, the torques of example 1, 2 and 3, which are given in Figs. 6c, 8c and 10c and corresponded to v = 1, 2, and 3 in Eq. (31) respectively, were applied to the full-nonlinear FEM.For the FEM, the link was modeled by 20 BEAM3 elements and the hub and end-effector mass were modeled by MASS21 from the ANSYS library [17]. BEAM3 is a uniaxial element with three degrees-of freedom (DOF) at each node; two translations and one rotation. It has the ability to model tension, compression and bending. MASS21 is a point element with three DOF; two translations and one rotation. Details of similar finite element analysis can be found in [11,14]. The differences between the desired end-effector displacement and those obtained from FEM are shown in Fig. 11. These differ~d ðtÞ was replaced ences were obtained from Eq. (41), in which y by the end-effector displacement obtained from the FEM. The maximum difference between the desired trajectory and the FEM trajectory when the torque of example 1, v = 1, was used, was 6.480  102 m. These maximum differences were 5.831  103 m and 1.212  103 m for the torques of example 2 and 3, v = 2 and v = 3 respectively. That is by increasing the number

Fig. 8. Example 2, redefining the desired trajectory with two segments, (a) desired end-effector trajectory (dashed) and redefined end-effector trajectory (solid), (b) difference between redefined and desired trajectories, (c) required torque.

M. Vakil et al. / Mechatronics 19 (2009) 1197–1210

of segments the maximum difference was reduced, as already has been observed. However, compared to the maximum differences obtained from Figs. 6b, 8b and 10b, these values were higher. The maximum differences from the simulation when the linear dy-

1207

namic model was used (the difference between the desired and redefined trajectories) and those of the FEM (the difference between the desired and FEM end-effector trajectories) are presented in Table 2. The main reason for the difference between the results

Fig. 9. Example 3, redefining the desired trajectory with three segments, (a) singular, feasible, and infeasible domains for the first segment (ti1 = 0.00(s) and tf1 = 0.40(s)), (b) (error)1 versus (mmax)1 and (mmin)1 for the first segment, (c) singular, feasible, and infeasible domains for the second segment (ti2 = 0.40(s) and tf2 = 0.8(s)), (d) (error)2 versus (mmax)2 and (mmin)2 for the second segment, (e) Singular, feasible, and infeasible domains for the third segment (ti3 = 0.80(s) and tf3 = 1.2(s)), (f) (error)3 versus (mmax)3 and (mmin)3 for the third segment.

1208

M. Vakil et al. / Mechatronics 19 (2009) 1197–1210

Fig. 10. Example 3, redefining the desired trajectory with three segments, (a) desired end-effector trajectory (dashed) and redefined end-effector trajectory (solid), (b) difference between redefined and desired trajectories, (c) required torque.

Fig. 11. Example 4, the difference between FEM end-effector trajectories and desired trajectory when the torque of example 1, 2, and 3 was applied to the FEM (v is number of segments used for redefinition).

of the linear dynamic model with one mode shape, and those of the FEM was due to the contribution of higher modes of vibration. Moreover in the finite element a nonlinear model, where the shear and axial deformations as well as the rotary inertia were considered was used while these deformations and also the rotary inertia were neglected in the linear dynamic model. That is the FEM results are more accurate. At the final manoeuvre time, tf = 1.2 s, there was a slight difference between the FEM end-effector trajectory and the desired endeffector trajectory. The absolute value of this difference for v = 1 was 5.025  104 m. These differences for v = 2 and v = 3 were 1.320  104 m and 7.950  104 m, respectively. Moreover, the velocity of the end-effector from the FEM results at tf = 1.2 s was slightly different than zero. The end-effector velocity at tf = 1.2 s for v = 1, v = 2 and v = 3 was 0.2127  102 m/s, 0.5361  102 m/ s and 0.1812  102 m/s, respectively. Thus, there was a small residual vibration. These differences between the desired and simulated (FEM) end-effector displacements and also the residual vibration were due to the fact that the linear dynamic model, which was used for the calculation of inverse torque, was not exactly the same as the FEM. There are two methods to have the FEM end-effector displacement approach the desired one and also suppress the residual vibration as follows. The first approach is to use an on-line state-feedback controller [16], which was done in

M. Vakil et al. / Mechatronics 19 (2009) 1197–1210

1209

Table 2 Example 4, comparison of the maximum differences between the desired trajectory and the FEM trajectory with the maximum differences between the desired trajectory and linear dynamic model trajectory.

v = 1 (one segment) v = 2 (two segments) v = 3 (three segments)

Finite element model (m)

Linear dynamic model (m)

6.480  102 5.831  103 1.212  103

6.320102 4.900103 1.355104

[6,8,18]. The second approach is to introduce an off-line technique which enhances the robustness of the new inversion method introduced here against the changes in the system dynamic model, similar to the robustness of the input shaping technique discussed in [19]. Addressing the latter approach is a challenging yet rewarding research subject, which its investigation is out of the scope of current paper. More simulation examples for the cases where there are positive zeros for the transfer function of the SFLM, positive eigenvalues for matrix AI in Eq. (16), or there is a combination of positive and purely imaginary zeros, positive and purely imaginary eigenvalues for matrix AI, have also been preformed to further check the validity of the method. The results of these studies were not reported here for briefness and can be found in Chapter 3 of [20]. In the following section, the application of the new inversion method, introduced in this paper, in the presence of the joint friction is discussed. This is important as the joint friction is almost inevitable in practice. This discussion is included here where the simulation study for the nonlinear model is presented, since friction makes the dynamics nonlinear. 7.4.1. Employing new inversion method in presence of joint friction In almost all experimental setups of SFLM, friction exists in the hub (joint). This friction, which is a combination of viscous friction and dry friction, makes the precise positioning cumbersome. Moreover, since dry friction adds nonlinearity to the system [21] and because the inversion method introduced in this paper is for linear systems, it is impossible to completely include the friction in the dynamics and use the new inversion method. To alleviate this problem, a remedy is proposed as follows. It is assumed that the friction torque sfriction can be modeled by the LuGre model [22] which is: 2

_ þ rh_ sfriction ¼ ðsfc þ ðsfst  sfc Þeðh=_ lÞ ÞsgnðhÞ

ð42Þ

In this model, which has been widely used to model friction in robot manipulators [20,23–27], the terms rh_ and ðsfc þ ðsfst  sfc Þ 2 _ _ are the viscous and the dry friction torques, respeceðh=lÞ ÞsgnðhÞ tively. Moreover sfc, sfst, l and r are the coulomb friction, static friction, Stribeck velocity constant, and the viscous damping coefficient, respectively. To identify the friction, the parameters sfc, sfst, l and r should be obtained experimentally, for example by using the approaches in [22,23,26,27]. _ does not add any nonlinearity, it Since the viscous friction, rh, can be included in the dynamic model and the inversion method discussed in this paper can be applied without further modification. An example where piece-wise inversion was applied to a SFLM with a joint viscous friction is available in Chapter 3 of [20]. To compensate for the dry friction, the friction torque 2 _ _ can be added to the torque obtained ðsfc þ ðsfst  sfc Þeðh=lÞ ÞsgnðhÞ using piece-wise inversion; note that sfc, sfst and l are available after identification and h_ is measured on-line. However, the complete compensation of the dry friction by feed-forwarding 2 _ _ is impossible. This is due to the fact ðsfc þ ðsfst  sfc Þeðh=lÞ ÞsgnðhÞ that the friction parameters change because of the change in load,

Fig. 12. Schematic of the SFLM closed-loop system, XInversion: the state vector obtained during inversion, sInversion: the required torque obtained during the the dry friction compensation torque which is inversion,sdry friction : 2 _ _ ðsfc þ ðsfst  sfc Þeðh=lÞ ÞsgnðhÞ.

temperature, humidity or operating condition. To compensate for this difference in the friction parameters and also to add robustness to the system, two approaches can be employed, which are also mentioned in Section 7.4. The first approach is to add a state-feedback controller sfeedback. This state-feedback controller is added to the dry friction compensation torque and the inversion torque, in which the schematic of the closed-loop system is shown in Fig. 12. The second approach is to enhance the robustness of the inversion method to the uncertainties in the dynamic model, similar to that discussed in [19] for the input shaping. Addressing these approaches, which are interesting research subjects, is, however, out of the scope of the current work. It is worth noting that utilizing the first approach (a combination of state-feedback controller with dry friction compensation torque) to enhance the control of rigid-link flexible-joint manipulators in the presence of joint friction has successfully been achieved at the Robotic Laboratory of the University of Saskatchewan [23].

8. Conclusions A novel causal end-effector trajectory inversion of a single flexible link manipulator {SFLM} through output redefinition was introduced which did not require post-actuation. The redefinition of the desired trajectory employed the summation of stable exponential functions which led to a family of possible solutions. Thus, by minimizing the error between the desired trajectory and the redefined trajectory, the best member of the family can be found. To find this best family member, a new and efficient optimization algorithm was introduced and used. The proposed algorithm is a zero-order optimization which finds the best member without any derivative calculation. Several examples were presented to show the effectiveness of this new optimization approach. Contrary to the available causal and non-causal methods for the inversion which do not tolerate the existence of the purely imaginary zeros for the transfer function between the end-effector displacement and applied torque, it was shown that the proposed method worked even in the presence of purely imagery zeros. The simulation results for a SFLM with purely imaginary zeros were presented. Although the method described here was only for the inversion of a SFLM, its extension to linear single-input single-output, nonminimum phase systems with hyperbolic or non-hyperbolic internal dynamics is straight forward. Another candidate for the output redefinition is the combination of the stable exponential functions with polynomial or sinusoidal functions, which is the subject of future research.

1210

M. Vakil et al. / Mechatronics 19 (2009) 1197–1210

(

Acknowledgements The authors would like to express their gratitude toward the referees, and Professor W. Szyszkowski, for their valuable comments. Appendix A. The elements of the matrices in the linear dynamic equation for a SFLM The linear dynamic equation of the single flexible link manipulator shown in Fig. 1 is:

X_ ¼ AX þ Bs ~d ¼ CX y

where matrices A, B and C are defined in Eq. (7). Furthermore, since X uI and X sI are found by the forward integration of Eqs. (18) and (19) ~d ¼ 0 for t < 0, the pair (s, X) are causal and s(t) = 0, and also since y ~_ ðtÞ ¼ 0 for t > t , ~ €d ðtÞ ¼ y X(t) = 0 for t < 0. Finally, since y d f _ ~d ðtÞ ¼ y ðt f Þ then sðtÞ ¼ 0; kðtÞ ¼ kðtÞ ¼ 0; X u ðt f Þ ¼ X s ðt f Þ ¼ 0 and y I

I

d

_ ¼ and y(t) = yd(tf) for t > tf where kðtÞ ¼ ½ k1 . . . kn  and kðtÞ   _k1 . . . k_ n . Therefore, based on the definition 1 given in Section 3, the causal end-effector inversion for a SFLM is achieved.

€ þ Kq ¼ Hs Mq or

References M

K

zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ 2 3 2 3 zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{ 2 32 3 €h M11 . . . M1nþ1 h K 11 . . . K 1nþ1 6 76 7 6 76 7 ... ... 54...5 þ 4 ... ... ... 54...5 4 ... €kn M nþ11 . . . Mnþ1nþ1 K nþ11 . . . K nþ1nþ1 kn H

zfflffl 2 ffl}|fflfflffl3{ 1 6 7 ¼ 4...5s 0

ðA1Þ

where

M 11 ¼ ðqL3 =3Þ þ Ih þ Itip þ mtip L2 M 1i ¼ M i1 ¼ q/i1;x þ mtip L/i1 ðLÞ þ Itip /0i1 ðLÞ;

i ¼ 2::n þ 1

M ij ¼ M ji ¼ q/i1;i1 þ mtip /i1 ðLÞ/j1 ðLÞ þ Itip /0i1 ðLÞ/0j1 ðLÞ; ( K¼ /i;j ¼

i; j ¼ 2::n þ 1 RL K ij ¼ EI 0 /00i ðxÞ/00j ðxÞdx;

0; Z L

i ¼ 1 or j ¼ 1 /j /i dx;

/i;x ¼

0

Z 0

i ¼ 2 . . . n þ 1; j ¼ 2 . . . n þ 1

L

x/i dx; /0i ðxÞ ¼

@/i ðxÞ ; @x

/00i ðxÞ ¼

ðA2Þ

@ 2 /i ðxÞ @x2

and q is the density per unit length, L is the length of the link, EI is rigidity of the link, Mtip and Itip are the mass and mass moment of inertia of the end-effector respectively and Ih is the mass moment of inertia of the hub. Appendix B. Proof of the theorem given in Section 5 Since the redefinition of the desired trajectory by SSEF is such that the complementary part of X uI from the following equation, Eq. (18) in Section 4:

~€ X_ uI ¼ AþI X uI þ BuI y d is zero, the causal

ðA3Þ X uI

is bounded. Also,

€~d X_ sI ¼ AI X sI þ BsI y

X sI

from Eq. (19), which is:

ðA4Þ

~ €d by SSEF. Thus: €d which is the redefinition of y is bounded for the y

h X ¼ T 1 y ~d

~_ y d

X TI

iT

ðA5Þ

and the corresponding torque:

s¼

1 ~ €d  CA2 XÞ ðy CAB

ðA6Þ

are bounded. Moreover, based on the inversion procedure adopted ~d and in the absence of the modeling error, initial condition erfor y ror and external disturbance, the pair (s, X) satisfy:

[1] Isidori A, Byrnes C. Output regulation of nonlinear system. IEEE Trans Automat Contr 1990;35(2):131–40. [2] De Luca A, Lanari L, Ulivi G. Tip trajectory tracking in flexible arms: comparison of approaches based on regulation theory. Lect Notes Contr Inform Sci 1991;162:190–206. [3] Yuan K, Liu LY. On the stability of zero dynamics of a single-link flexible manipulator for a class of parametrized output. J Robot Syst 2003;20(10): 581–6. [4] Devasia S, Chen D, Paden B. Nonlinear inversion-based output tracking. IEEE Trans Automat Contr 1996;41(7):930–42. [5] Bayo E. A finite-element approach to control the end-point motion of a single link flexible robot. J Robot Syst 1989;4(1):63–75. [6] Kwon DS, Book WJ. Time-domain inverse dynamic tracking control of a singlelink flexible manipulator. J Dyn Syst Measur Contr 1994;116(2):193–200. [7] Devasia S. Approximated stable inversion for nonlinear system with nonhyperbolic internal dynamics. IEEE Trans Automat Contr 1999;44(7): 1419–25. [8] Bensoman M, LeVey G. Stable inversion of SISO nonminimum phase systems through output planning: an experimental application to the one-link flexible manipulator. IEEE Trans Contr Syst Technol 2003;11(4):588–97. [9] Wang X, Chen D. Output tracking control of a one-link flexible manipulator via causal inversion. IEEE Trans Contr Syst Technol 2006;14(1):141–8. [10] Vakil M, Fotouhi R, Nikiforuk PN. Causal Inversion of a single-link flexible link manipulator via output planning. In: Proceedings of the IEEE international conference on mechatronics and automation, Ontario, Canada; 2005. p. 376– 81. [11] Fotouhi R. Dynamic analysis of a very flexible beam. J Sound Vib 2007;305(3):521–33. [12] Zhang X, Xu W, Nair SS, Cellabonia V. PDE modeling and control of a flexible two-link manipulator. IEEE Trans Contr Syst Technol 2005;13(2):301–12. [13] Book WJ. Recursive Lagrangian dynamic of flexible manipulator arms. Int J Robot Res 1984;3(3):87–101. [14] Vakil M, Fotouhi R, Nikiforuk PN, Salmasi H. A constrained Lagrange formulation of multilink planar flexible manipulator. ASME J Vib Acoust 2008;130(3):1–16. [15] Slotine JJ, Li W. Applied nonlinear control. NJ: Prentice Hall; 1991. [16] Khalil HK. Nonlinear Systems. NJ: Prentice Hall; 2002. [17] ANSYS, Release 11, Release 11 Documentation for ANSYS, element reference. . [18] Paden B, Chen D, Ledesma R, Bayo E. Exponentially stable tracking control for multi-joint flexible-link manipulators. ASME J Dyn Syst Measur Contr 1993;115:53–9. [19] Singhose WE. Command generation for flexible systems. Ph.D. thesis. Cambridge, MA, USA: Massachusetts Institute of Technology; 1997. [20] Vakil M. Dynamics and control of flexible manipulators. Ph.D. thesis. Saskatoon, SK, Canada: University of Saskatchewan; 2008. [21] Lawrence J, Singhose W, Hekman K. Friction-compensating command shaping for vibration reduction. ASME J Vib Acoust 2005;127:307–14. [22] de Wit CC, Olsson H, Astrom KJ, Lischinsky P. A new model for control of systems with friction. IEEE Trans Automat Contr 1995;40(3):419–25. [23] Salmasi H, Fotouhi R, Nikiforuk PN. Tip trajectory tracking of flexible-joint manipulators driven by harmonic drives. Int J Robot Automat 2009;24(2). [24] Bona B, Indri M, Smaldone N. Rapid prototyping of model-based control with friction compensation for a direct-drive robot. IEEE/ASME Trans Mech 2006;11(5):576–84. [25] Liu M, Quach NH. Estimation and compensation of gravity and friction forces for robot arms: theory and experiments. J Intell Robot Syst 2001;31:339–54. [26] Kelly R, Llamas J, Campa R. A measurement procedure for viscous and Coulomb friction. IEEE Trans Instrum Measur 2000;49(4):857–61. [27] Kermani MR, Patel RV, Moallem M. Friction identification and compensation in robotic manipulators. IEEE Trans Instrum Measur 2007;56(6):2346–53.