Modeling and flatness-based control of a 3d of helicopter laboratory experiment

Copyright © IFAC Nonlinear Control Systems, Stuttgart, Germany, 2004

ELSEVIER

IFAC PUBLICATIONS www.elsevier.com/locare/ifac

MODELING AND FLATNESS-BASED CONTROL OF A 3DOF HELICOPTER LABORATORY EXPERlMENT Thomas Kiefer, Andreas Kugi and Wolfgang KemmetmiiIIer

Chair of System Theory and Automatic Control, Saarland University, PO Box 15 11 50, 66041 Saarbriicken, Germany, E-Mail: [email protected]

Abstract: This paper deals with the modeling and control of a three-degrees-offreedom helicopter laboratory experiment. The helicopter belongs to the class of mechanical systems underactuated by one control. The mathematical model is derived using the Lagrange formalism in combination with the concept of twists and wrenches. It can be proven that the full system is not configuration flat. Nevertheless, by a slight modification of the generalized forces, which can also be interpreted in terms of a constructive change in the experimental set-up, we are able to design a flatnessbased controller. Experimental results demonstrate the effectiveness of the proposed concept. Copyright © 2004 lFAC Keywords: Underactuated Mechanical Systems, Helicopter Laboratory Experiment, Nonlinear Control, Flatness-based Control, Configuration Flatness

1. INTRODUCTION

e.g., (Murray et al., 1994), (Stramigioli, 2001) that allow a global description of rigid body motions which does not suffer from singularities due to the use of local coordinates. Given the kinetic and potential energy of the helicopter we will derive the equations of motion using Lagrange's formalism.

The present contribution is devoted to the modeling and flatness-based control of the laboratoryexperiment helicopter 3DOF of Quanser, see (Quanser, 2002). The helicopter as shown in Fig. 1 consists of the base upon which the arm is mounted. The arm carries the helicopter body containing two propellers driven by dc-motors on the one end and a counterweight on the other end. The three parts are connected via revolute joints, giving the following three degrees-of-freedom: the travel angle ql, the elevation angle q2 and the pitch angle q3. The helicopter can be controlled by means of the thrusts Fb and FI which are generated by the two propellers. In the subsequent modeling and controller design we will neglect the fast dynamics of the dc-motors and use the thrusts as control inputs. For the description of the configuration of the helicopter we will use exponential coordinates expressed by twists and wrenches, see,

The system is higWy nonlinear and underactuated by one control. This is why the design of a controller for the motion control of this laboratory experiment helicopter is much more complicated compared to other helicopters which are mechanically fully actuated. Thus, for instance in (Mullhaupt et al., 1999) the non-flatness is due to a cross-coupling stemming from the propeller acceleration. In absence of this cross-coupling term the helicopter model would be flat. In contrast to this, we have to deal with a mechanical system where the number of degrees-of-freedom is less than the available control inputs. In (Rathinam

207

in a subspace of ~* Q which only depends on the configuration q and thus can be described by a codistribution P c ~* Q. A Riemannian metric G on the manifold Q is a second-rank covariant tensor which assigns, in a differentiable fashion, a positive definite inner product (-,.) in each tangent space Tq Q. In terms of a coordinate basis we denote the components of G by gij and hence we have 1, see, e.g., (Dubrovin et al., 1992), (Frankel, 2001)

and Murray, 1998) the notion of configuration flatness was introduced for Lagrangian systems underactuated by one control. Thereby, a Lagrangian system is said to be configuration flat if it is differentially flat with flat outputs that only depend on configuration variables. Applying this theory to our helicopter, we can show that it is not configuration flat. Furthermore, the system proves to be not exact input-state linearizable. Nevertheless, utilizing the constructive character of the method proposed in (Rathinam and Murray, 1998), we are able to render the resulting system configuration flat by a slight modification of the generalized forces. Moreover, these modifications can even be interpreted in terms of a constructive change in the experimental set-up. Based on these results a flatness-based controller is developed and implemented in a real-time environment. Finally, measurement results show the feasibility of the proposed concept.

v, wE TqQ .

(v, w) = %viwi

(2)

Since G is positive definite, we may also calculate its inverse C- 1 , in coordinates denoted by gik, where gik gkj = 8j. A manifold with a Riemannian metric is called a Riemannian manifold. Given a metric G the covariant derivative 'V provides an intrinsic way to differentiate objects on the manifold. In local coordinates q the covariant derivative 'V v w of a vector field w = wi a~j along a vector field v

=

v j a~) reads as

.

k

.ow k

0

.

0

'Vvw = vJw rjk~ + vJ~7lk' uq'

uqJ uq

(3)

where rjk are the Christoffel symbols of the second kind, i, j, k = 1, ... , n, (Dubrovin et al., 1992). Note that the covariant derivative 'V vi of a function i along the vector field v = v j a~) corresponds

fJ4.

to the well-known Lie derivative 'V vi = v j An operation of covariant differentiation is afso called a connection. If the metric G is non-singular (det (gij) :/= 0), then there exists a unique symmetric connection = r~j' which can be computed from the metric G in the form

r;k

ri _ ~ g il Jk - 2

2. PRELIMINARIES

.

(4)

Let us describe the mechanical system under consideration on an n-dimensional Riemannian manifold Q. Then the metric G is defined in a natural way by the kinetic energy T in the form 2T = (q, q), q E TqQ. Clearly, the components gij of G correspond to the entries in the generalized inertia matrix D, see (5). Further, we will only consider so-called natuml mechanical systems whose Lagrangian has the form kinetic minus potential energy and whose potential energy function U is solely a function of the generalized coordinates q, (Arnold, 1989). Then the equations of motion can be written in matrix notation in the well-known form

Let us consider a mechanical system with n degrees-of-freedom which can be locally described by the generalized coordinates q = (ql, ... , qn) of the configuration manifold Q. The tangent bundle T Q is the collection of all tangent spaces Tq Q at all points q E Q. Similarly, we call the space of all cotangent spaces ~* Q to Q the cotangent bundle and write T* Q. The Lagrangian L (q, q) can be regarded as a real valued function from the tangent bundle T Q to IR with the generalized velocities q E TqQ. The well-known Euler-Lagrange equations read as

Q,

Ogjl _ o9jk } oql

+ oqk

This connection is also known as the Levi-Civita connection, see, e.g., (Dubrovin et al., 1992), (Frankel, 2001) for more details.

Fig. 1. Laboratory experiment helicopter 3DOF.

:t (~~) - ~~ =

{091k oqj

D (q) q + C (q, q) q + 9 (q)

(1)

where Q denotes the generalized forces acting on the system. For all subsequent considerations we will assume that the generalized forces Q lie

= Q,

(5)

1 In order to kcep the formulas short and readable, we will subsequently use Einstein's convcntion for sums by summing over double repeated indices.

208

with the generalized inertia matrix D (q), the Coriolis matrix C (q, q), the gravity vector 9 (q) and the generalized forces Q. Note that the Christoffel symbols rjk from (4) are related to the Coriolis matrix by Cij = gilr~kqk and the gravity vector 9 (q) can be derived from the potential energy U (q) by the relation gi = a~' U (q).

referred to as the adjoint transformation Adg associated with g, given by Adg

T

°

= Adg-t = [R -RTd] RT

For open-chain rigid body systems with revolute joints as it the case for the helicopter of Fig. 1, we have a simple rule for determining the twists, sce (Murray et al., 1994). Let us define with = q2 = q3 = the reference configuration of the helicopter with the inertial frame (0o, Xo, Yo, zo). Then, the twist vector belonging to the i th revolute joint, i = 1,2,3, has the form

ql

°

In this section, we will use the concept of twists and wrenches which allows to derive the equations of motion in a systematic and simple way. For more details on this topic the reader is referred to (Murray et al., 1994) and (Stramigioli, 2001). It is known that in general the motion of a rigid body can be described by a combination of a rotation and a translation. The configuration (position and orientation) of a coordinate frame K 1 attached to the body relative to an inertial frame Ko is described by the position vector dli E 1R3 from the origin of Ko to the origin of K 1 and the rotation matrix RA E SO (3) = {R E 1R3x31RRT = I, det (R) = 1 }. Thus, the configuration of the rigid body is given by the pair gJ = (dli, RA) E SE (3) = {(d, R) E 1R3 , RE SO (3)}, with SE (3) the special Euclidean group. The homogenous representation of gJ is written in form of a (4 x 4)matrix 2

(10)

with Wi E 1R3 as the unit vector directed into the rotation axis of the i th coordinate frame (Oi,Xi,Yi,Zi), i = 1,2,3, and ai E 1R3 is any point on the rotation axis, both of them measured in the inertial frame. For example for the revolute joint 3 we have wj = [-l,O,OJ and af = [0,0,ar2 - ahJ. By this procedure we may calculate the twist vectors for the helicopter of Fig. 1 in the form

~r = [0,0,0,0,0,1]' ~r = [0, a23 -

~r = [ah,O,O,O,-l,O]

ah, 0, -1,0, 0J

.

(11) By means of these twist vectors we can directly derive the equations of motion according to (5). For this we define new coordinate frames Ei , i = 1,2,3, each of them attached to the center of mass of the i th subsystem. The configuration of the frame ~ relative to the inertial frame in the case when the helicopter is in its reference configuration will be denoted by g~' (0). The socalled product of exponentials formula g~' (q) =

Id

.

-I

Ad g

(9)

3. MATHEMATICAL MODEL OF THE HELICOPTER

gJ = [ ~ ~~]

= [ ROdR] R'

(6)

Now, a twist is an element of the set se (3) = {(v,w) E 1R3 ,w E so (3) }, where so(3) denotes the set of skew symmetric (3 x 3)-matrices, i.e. so (3) = {S E 1R3X31ST = -S}. For a differential geometric interpretation in terms of a Lie algebra, see (Stramigioli, 2001). By using the so-called wedge operator 1\ applied to a vector W E 1R3

Iv

ql) ...

exp (~I exp (~iqi) g~' (0) describes the actual configuration of the frame E i for q i= 0. Furthermore, let M i be the generalized inertia matrix of the i th subsystem with respect to the frame ~' Then, the generalized inertia matrix D (q) of (5) is given by 3

D (q) =

we can identify so (3) with 1R3 . Similarly, a twist ~ E se (3) can be written in homogenous coordinates in the form of a (4 x 4)-matrix .

~=

v]

°°

[£.21

L (Jt{ Mdt ,

(12)

i=1

with the so-called body Jacobian

Jt(q) =

(8)

[d, ... ,d,o,o,o]

(13)

where

~tJ

or in the form of a twist vector ~T = [v T , w T ] E 1R6 . The (6 x 6)-matrix which transforms twist vectors from one coordinate frame to another is

= Ad-(l (.) (.);. )~j exp f,jqj ... exp f"q' go' (0)

(14)

for j ~ i = 1,2,3, see, e.g., (Murray et al., 1994). The Coriolis matrix C (q, q) can be directly calculated by means of (4) and the gravity vector 9 ( q) follows from the potential energy. Since this is rather straightforward we will omit a detailed derivation here. Remains to show how the thrusts

2 Note that here and subsequently we do not consequently distinguish between the element and its matrix or vector representation.

209

.

of the two propellers, F b and F / , act as generalized forces Q on the system (5). In general, a generalized force F is composed of a linear component (pure force) f E ]R3 and an angular component (pure torque) T E ]R3 and thus can be represented as a vector F = [IT, TT] T in ]R6, also referred to as a wrench. Wrenches naturally combine with twists to define instantaneous work. Therefore, they also constitute a pair of so-called power variables. The relation Fj

= Ad~Fi

to a diagonal matrix with constant entries d j j , = 1,2,3, i.e. D (q) = diag (d u ,d22 , d33) 3.

j

4.1 Non-flatness Let us consider a natural mechanical system with n degrccs-of-freedom underactuated by one control with the configuration manifold Q and a metric (.,.) defined by the kinetic energy, see Section 2. Remember that we have supposed that the generalized forces Q can be described by a codistribution P C ~*Q. Since the system under consideration is underactuated by one control we suppose dim (P) = n -1 at a generic point q E Q. For the subsequent considerations we will further define a distribution V in the form

(15)

J

transforms a wrench F i applied to the origin of the coordinate frame K i into an equivalent wrench applied at the origin of the coordinate frame K j . The thrusts F b and F I define the following wrenches PI = [0,0, F b , 0, 0, 0] and Fl = [0,0, F/, 0,0,0] with respect to the coordinate frames (04,X4,Y4,Z4) and (Os,xs,Ys,zs), respectively. Clearly, the corresponding twists ~4 and ~s are identically zero. Thus, by utilizing (15) the generalized forces QT = [Ql,Q2,Q3] in (5) can be directly calculated in the form

V = span {(, V...L(, i = 1, ... Bo'

,n},

(17)

where ( is any vector field such that pJ. = span {(}, with pJ. the annihilator of P. Now we will restate Proposition 2 of (Rathinam and Murray, 1998) which provides a constructive method for the determination of configuration flat outputs y 1 , ... ,yn - 1 :

Theorem 1. Let U be an open neighborhood of a point q E Q and suppose y : U C Q -+ ]Rn-l is a submersion. If yl, . .. , yn-l are configuration flat outputs, then

(16)

(kerTY, V) where Or denotes the rotational component of J the wrench (torque) in the j-axis, j E {x,y,z}. Note that the whole mathematical model of the helicopter is too large to be presented here. The simulation model implemented in MATLABjSIMULINK contains in addition a stick-slip friction model for all three revolute joints.

= 0,

(18)

where in coordinates ker TY is the kernel of the Jacobian of the map y. Conversely, given a set of outputs yl, ... , yn-l and suppose that the function z completes the outputs to a coordinate system. If (kerTY, 'D) = 0 at a point q E Q and the expressions

/;(((,~)) /;((V*~,())

(, a~J)

4. FLATNESS-BASED CONTROLLER DESIGN

(Va7~'()

(19) for j, k, l = 1, ... , n - 1 are not all identically equal in a neighborhood of q, with U as the potential energy function, then yl, ... , yn-l are configuration flat outputs around q. Thereby, ( and TJ are fixed nonvanishing vector fields such that pJ. = span {(} and ker TY = span {TJ}.

Roughly speaking, a system is differentially flat, if there exists a set of outputs equal in number to the number of inputs such that all states and inputs can be determined from these outputs without integration. A rigorous definition of flatness can be found in e.g., (Fliess et al., 1995), (Rudolph, 2003). For the purpose of a controller design, we try to simplify the mathematical model in such a way that it still contains the essential nonlinearities but it becomes manageable for symbolic manipulations. Neglecting terms in the equations of motion often leads to a mathematical model which is no longer physically meaningful. Therefore, we will rather try to simplify in terms of the energies, where the physical structure of the system is preserved. In this way we find out that the kinetic energy can be shortened such that the generalized inertia matrix D (q) simplifies

For a rigorous proof of this Theorem the reader is referred to (Rathinam and Murray, 1998). Next, we will apply Theorem 1 to the mathematical model of the helicopter. The thrusts Fb and FI as the external control inputs lie in the codistribution P = span bldql - cos (q3) a~3dq2 + a~4dq3, 'Y2dql - cos (q3) a~3dq2 -

aXs dq3}

, (20)

Numerous simulation results justify these simplifications.

3

210

where 11

= -

sin (q3) cos (q2) a~3+

sin (q2) (a23 sin (q3) - a~4) - sin (q3) cos (q2) a~3 + sin (q2) (a23 sin (q3)

(21)

+ a~5)

The annihilator pi- is spanned by the vector field (cos (q2) ai3 - sin (q2) ah) sin(q3) 8 cos (q3) a X 8 q2

8

( = 8q1 -

8

+ sin (l)

8

q

23

3 .

(22) Now it is possible to calculate the distribution V from (17). Note that the simplifications made so far drastically simplify the evaluation of the covariant derivatives according to (3). Since the generalized inertia matrix was reduced to a constant diagonal matrix diag (d ll , d 22 , d33 ) the Christoffel symbols fjk from (4) are all equal to zero 4 . Thus, the covariant derivative of ( along the vector field 8~i takes the form 3

'V .L ( = 8.'

L .

J=1

8

a:;T

(,

'V b

8q

thrusts, Fb and Ft, would always lie in a plane parallel to the zo-direction. The simplification being proposed leads to a codistribution P equal to (20) except for 11 and 12 which change to 11

-8 . -8. . qt qJ

C 'V ~ i)

(}

= 12 = -

sin (q3) cos (q2) a~3 .

(25)

In this case the annihilator pi- is given by

8(J 8

(23)

pi- = span

It can be easily verified that the dimension of V

span {C 'V

Fig. 2. Constructive change to render the helicopter configuration flat.

=

{8~1

-

tan (q3) cos (q2)

8~2}

(26)

and dim (V) = 2 at a generic point q. Furthermore, we can show that the vector field TJ = 8~3 satisfies the condition (TJ, V) = O. Since the onedimensional distribution span {TJ} is always integrable locally, we find two independent functions, yl = Xl (q1,q2) and y2 = X2 (ql,q2), which serve as possible flat outputs. However, choosing the functions yl = q1, y2 = q2 and adding the function z = q3 to give a complete set of local coordinates (y1, y2, z) for the configuration manifold Q, it can be shown that the regularity conditions (19) are trivially satisfied. As a result we may conclude that y1 = ql and y2 = q2 are flat outputs for the modified helicopter system.

is 3 at a generic

point q. As a result the system is not configuration flat because we cannot find a nontrivial vector field TJ such that (TJ, V) = O.

4.2 A Flatness-based Approximation From Theorem 1 we see that the rank condition dim (V) s:: n - 1 at a generic point q provides a necessary condition for the existence of configuration flat outputs. Therefore, we are looking for conditions on the vector field ( (see (22)) of the helicopter model in order to reduce the dimension of the distribution V at least by one. A straightforward calculation shows that this is the case when the components of the vector field ( satisfy the PDE-condition

4.3 Controller Design Considering all the simplifications discussed so far the mathematical model of the helicopter reads as

(24) In addition we want the components (J, j = 1,2,3, to be independent of q1 in accordance with (22). Hence, (3 must either be a function of q3 or (2 a function of q2 only. From a geometric point of view it becomes immediately clear that the first condition is approximately fulfilled in the operating area. Moreover, it is even possible to fulfill this condition exactly by a constructive change in the experimental set-up of the helicopter as shown in Fig. 2. Note that with this modification the

iP = - da3

33

y

cos (q2) sin (q3)

+ a34u2 , d33

(27c)

with the new control inputs U1 = Fb + Ft, U2 = F b - Ft and the coefficients al, a2, a3 depending on the masses and on geometric parameters. The derivation of the flatness-based controller for the flat outputs yl = q1 and y2 = q2 is straightforward. Let qid denote the desired trajectories

In this case the connection is said to be Euclidean, see (Dubrovin et al., 1992).

4

211

for qi, i = 1, 2, and suppose they are sufficiently smooth. In a first step, we choose a linear error dynamics for the tracking error of the elevation angle e2 = q2 - q2d in the form

iF =

e

2 ij2d - k 22 2 - k 21 e - k 20 e2I , e 2I

=

f

400r-,--.,.---,-----,--;:=====il 350 300

e 2dt ,

250

(28) with suitable constant coefficients k 2j > 0, j = 0,1,2. Equating (27b) with (28), we obtain the first part of the flatness-based controller U1

= U1

100 50

(q2d,q2d,ij2d,q2,q2,q3,e 21 ) .

(29)

o

In the next step (27a) is differentiated two times w.r.t. t and by eliminating iF and ij3, we get (q

1) (4)

=

f-2 ('2 ··· ) q ,q2 ,q·3 ,q3 ,Ul,U1,U1,U2

,

-50 0L--:'::-0--2O"e---"'3O,----'4O,---50'---60'----l70 1 tin 5

(30) Fig. 3. Tracking behavior of the travel axis.

Jt: (.)

where the following abbreviation (.)(n) = is used. Again by choosing a linear error system for the tracking error of the travel angle e 1 = q1 _ q1d (q

1)(4) -_

5

-measured

(ld)(4) q - k 14 (1)(3) e - k 13e··1 - k 12e·1 1 -k l l e - klOel/, ell = e 1dt,

- - -desired

~

h o

M

f

..

·5

(31) with suitable constant coefficients k 1j > 0, j = 0, ... ,4, we can directly derive the control law for U2 from (27) - (31) in the form U2 =

I

.. ,

[

oS

"b- -10

...

-15

Id ·ld ··ld (ld)(3) (ld)(4) 2d ·2d ··2d U(2 q , q , q , q ,q ,q,q,q, 1·1 2 ·2 3 ·3 II 21) 2d)(3) (2d)(4) , q ,q ,q ,q ,q ,q ,q,e ,e . (q

-20

(32) Note that the explicit expressions of the flatnes&based controller are too large to be presented here.

l,..J o

10

..

U 20

L..i

30

40

..

,

\-I

50

1.0;,

60

70

tin 5

Fig. 4. Tracking behavior of the elevation axis. Fliess, M., J. Levine, Ph. Martin and P. Rouchon (1995). Flatness and defect of non-linear systems: introductory theory and examples. Internat. J. Control 61, 1327 1361. Frankel, T. (2001). The Geometry of Physics. Cambridge University Press. Cambridge. Mullhaupt, P., B. Srinivasan, J. Levine and D. Bonvin (1999). Cascade control of the toycopter. CD-Proc. of the European Control Conference ECC'99, 31.08. -03.09.1999. Murray, R.M., Z. Li and S.S. Sastry (1994). A mathematical introduction to robotic manipulation. CRC Press. Boca Raton. Quanser (2002). 3D Helicopter System, product information. www.quanscr.com. Ontario. Rathirlam, M. and R. Murray (1998). Configuration flatness of lagrangian systems underactuated by one control. SIAM J. Control and Optimization 36(1), 164-179. Rudolph, J. (2003). Beitmge zur fiachheitsbasierten Folgeregelung linearer und nichtlinearer Systeme endlicher und unendlicher Dimension. Shaker Veriag. Aachen. Stramigioli, S. (2001). Modeling and IPC Control of Interactive Mechanical Systems: A Coordinate Free Approach. Springer. London.

5. MEASUREMENT RESULTS The control concept is implemented in a real-time environment of oSPACE. The sampling time is chosen sufficiently small (1 ms) in order to ensure a quasi time-continuous operation of the controller. The reference trajectories are generated by using Gevrey functions, sce, e.g., (Rudolph, 2003). The measurement results depicted in Fig. 3 show the reference tracking behavior of the travel angle q1 while the elevation angle q2 is kept at a constant value. Fig. 4 presents the measured response to a reference command of the elevation angle q2 holding the travel angle q1 at a constant value. In conclusion we may say that the flatneS&- based controller presented in this paper proves to be very efficient in terms of tracking arbitrary flight trajectories within the possible operating range of the laboratory helicopter, see also the videos on our web page http://www.lsr.unisaarland.de/heli.html. 6. REFERENCES Arnold, V.I. (1989). Mathematical Methods of Classical Mechanics. Springer. New York. Dubrovin, B.A., A.T. Fomenko and S.P. Novikov (1992). Modern Geometry - Methods and Applications. Springer. New York.

212