SciELO - Scientific Electronic Library Online

vol.22 issue3Development of a three-dimensional nonlinear viscoelastic constitutive model of solid propellantA study of the facial aging - a multidisciplinary approach author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand



Related links


Journal of the Brazilian Society of Mechanical Sciences

Print version ISSN 0100-7386

J. Braz. Soc. Mech. Sci. vol.22 n.3 Campinas  2000 

Nonlinear Dynamics and Control of Multibody Systems


Fernando A. Rochinha
Mechanical Engineering Department. EE-COPPE-UFRJ. Cx. Postal 68503. 21945-970 Rio de Janeiro. Brazil
Rubens Sampaio
Mechanical Engineering Department. PUC-Rio. Rua Marquês de São Vicente 225. 22453-900 Rio de Janeiro. Brazil



The dynamics of flexible systems, such as robot manipulators , mechanical chains or multibody systems in general, is becoming increasingly important in engineering. This article deals with some nonlinearities that arise in the study of dynamics and control of multibody systems in connection to large rotations. Specifically, a numerical scheme that adresses the conservation of fundamental constants is presented in order to analyse the control-structure interaction problems.




Recently there has been a great interest in the study of nonlinear dynamics of structures formed by the connection of ''simple bodies'' and its applications to a wide variety of engineering problems: robotics, spacecraft dynamics and attitude control,vehicle dynamics and large structures vibration. Indeed, there is a growing literature (Garcia de Jalón and Bayo 1993, Hughes 1986 and references therein) in such kind of systems usually called multibody systems. In particular, the approach to the rotational degrees of freedom plays a crucial role in the modeling and numerical approximation of rigid bodies, rods and shells, which are the most common elements forming a multibody. Thus, in order to better establish and assess the performance of some specific approach, the dynamics of a single rigid body is often used as a prototype problem. From a numerical standpoint it represents a significant example due to its high degree of nonlinearity.

The dynamic response associated with the rotational degrees of freedom leads to a problem evolution in the rotation group, which can be parametrized in several different ways (Hughes 1986). In the present work the rotations are described by means of directors (Simo and Wong 1991), what will play a crucial role in the developments presented here.

Mechanical control systems provide an important and challenging research area that falls between the study of classical mechanics and control theory. A very illustrative example is the control of multibody systems in the context of spacecraft attitude. Often one part of a spacecraft (a communication antenna) must point toward the Earth, while another part (a solar panel) must face the sun. To achieve such mission objectives, it is evident that an attitude stabilization and control system must be an important part of the design.

The present article should be understood as the first attempt of the authors to conjugate their own experiences in the nonlinear analysis of structures with control. Thus, in the present work, a simple Proportional Derivative (PD) feedback law is proposed with a respective mathematical analysis based on the use of Lyapunov's formalism. The novelty of the proposed law relies on the treatment of the rotational degrees of freedom, which represents, as was mentioned earlier, a source of difficulties. In order to better understand the role played by the interaction of control and nonlinear mechanics only the rigid body modeling is adressed. Although, it is believed that all experience obtained here can be directly transported for nonlinear rods.

The numerical simulation seems to play a very important role in the design of nonlinear control systems (Bullo and Murray 1997, Fjellstad and Fossen 1994 and Slotine 1991). Thus, the proposed modeling is tested by means of two numerical simulations, in which a conserving algorithm proposed by the authors (Rochinha and Sampaio 2000) for the time integration of the nonlinear dynamics of rigid bodies is applied.


Mechanical Modelling

In the following, a brief summary of the dynamics of a rigid body is presented. Further details can be found in the standard literature (see for example Goldstein 1980).

Let Î be the reference placement of a solid body, with particles labeled by X. So, the motion of the body is described by the mapping :


where t denotes the time, T represents the total time of observation and is the ordinary 3D-Euclidian space. Throughout this paper, bold letters are used to designate vectors.

As in a rigid body the distance between any two points must remain the same during the motion, the position of each particle can be put in the following form:

a08frm02.gif (2657 bytes), (2)

where r defines the position of the center of mass and di are vector fields called directors (Simo and Wong 1991) forming at each instant t an orthonormal basis attached to the body. Within the rigid body literature, the basis {di} is known as the body frame.

In (2), r takes into account the translation of the body and its rotation is described by the comparison among the directors and a fixed (inertial) basis { ei}. Thus, the configuration of the body is represented by the pair:


where dij repesents de delta of Kronecker. Consequently, one refers to x K as the abstract configuration manifold of the rigid body.

Remark 1: Indeed, K and the special group of rotations S0(3) = {Q Î / QQT = QTQ = 1, det(Q) = +1}can be identified by the isomorphism:

di = Q Di (i = 1,3), (4)

where Q is a rotation belonging to S0(3), 1 is the identity in and Di are the directors in the reference placement (t = 0), which, without loss of generality, can be taken coincident with the principal directions of inertia. Hence, the rotation can be recast as Q = di Ä ei, where Ä represents the tensorial product.

The velocities and accelerations fields in the present theory are:

· Translational Motion              · Rotational Motion

a08frm05.gif (1920 bytes)

where a08img01.gif (580 bytes) stands for the time derivative of ( ). The first two fields are, respectively, the translational velocity and the translational acceleration of the center of mass. The others describe the velocity and acceleration of the directors and are related to the angular velocity w and the angular acceleration a of the body by means of :



where denotes the ordinary vector product.

The linear and angular momenta which are given, respectively, by the following expressions:

where r is the mass density,play a crucial role i n the understanding of the dynamics of a rigid body. Using the relation (2) in the above expressions, the linear and angular momenta, in the present context, are done respectively by:

L = Mv, (7)

where M is the total mass given by and: = and:

a08frm09.gif (10109 bytes)

where the scalars Iij = Iji are the components of the moment of inertia tensor, e.g.:

a08frm10.gif (2642 bytes)(9)

defined by:

a08frm11.gif (5087 bytes)

Throughout this paper, the directors will be taken, without loss of generality, coincidents with the principal directions of inertia, which implies that Iij introduced above are 0 for i ¹ j.

Remark 2: By using (5) and the relations di a08img07.gif (993 bytes) = w - (w.di) di and di a08img09.gif (932 bytes)j = - (w.di ) dj , the above expression for the angular momentum can be cast in the classical form of the rigid body dynamics.

J = L + p,

where p = Iw is the spatial angular momentum relative to the center of mass.

Let f and mc be, respectively, the applied force and the applied torque with respect to the center of mass. In the present context, the classical equations of balance of angular and linear momentum take the form:

a08frm12.gif (2144 bytes)

In particular, the applied torque will be written in the following form, which is convenient for the numerical technique used in the present work.

where fi should be understood as the force which generates the torque mc.

The concepts introduced above lead to the following nonlinear system of ordinary differential equations governing the motion of a rigid body:

Problem: Find F =   and    such that:

Ma = f. (10)

Ia + w Iw = mc. (11)

a08img07.gif (993 bytes) = w di (i = 1,3), (12)

where dK is the tangent space of K (Le Tallec et all. 1992).

The above problem is recast in a variational formulation inspired in Le Tallec et all.(1992) and Rochinha (1990) and numerically solved by means of a techinque presented in Rochinha and Sampaio (1996-2000).


The Feedback Control Law

In the present section two basic categories of control task systems are treated, namely: the regulation (or stabilization) and tracking. In regulation problems a control system, called a regulator, is to be designed so that the state of the closed-loop system will be stabilized around an equilibrium point. In tracking control problems, the design objective is to construct a controller so that the system output tracks a given time-varying trajectory. In both cases, nonlinear control laws exploring the description of large rotations by means of directors and the Lyapunov's formalism (see Kelkar and Joshi, 1996, Slotine and Li, 1991 and Junkins and Kim,1993) are presented. In order to isolate difficulties, those laws are designed only for pure rotational motions.


The Regulation Problem

A proportional derivative feedback control law (PD) for the regulation problem in the context of the presented modeling is introduced. For this purpose, the tracking error is defined as:

a08frm16.gif (1776 bytes) (13)

The stability of the closed-loop system is explored with the aid of the following Lyapunov function:

a08frm17.gif (2442 bytes), (14)

where Kp > 0 is the proportional gain associated to the feedback law. It is straightforward to verify that the scalar function defined in (14) is positive-definite "t Î [0, ). Taking the time derivative of V and using the definitions of angular acceleration and the inertia tensor yields.

a08frm18.gif (3474 bytes)

Using the balance of angular momentum (11), the term multiplied by the angular velocity in the right side of the equation above is substituted by the applied momentum, resulting in:

a08frm19.gif (3588 bytes).

After rearranging the terms, the time derivative of V is rephrased as:

a08frm20.gif (3661 bytes).

Now, adopting the following proportional-derivative feedback law in the equation above:

a08frm21.gif (1684 bytes),


a08frm22.gif (2395 bytes),

where Kv is the derivative gain. The expression above implies that for each trajectory of the rigid body, the function V will not increase. Thus, as V(t) ³ 0 and a08img06.gif (928 bytes)(t) ³ 0, "t Î [0, ¥ ), given any initial state represented by a08frm23.gif (2019 bytes) in which V (0) ¹ 0 the final state will be an equilibrium point such that V(tf) = 0, where tf is associated with the final state. Hence, the main question to be answered is that if there is an unique equilibrium point of the system to be reached by the rigid body. In fact, it will be shown next that the state a08frm24.gif (5197 bytes) which is the target of the regulation point is the only attractor for the closed loop system. By the way there are few points (precisely seven points) which satisfy a08img06.gif (928 bytes)(t) = 0. Those points are considered isolated equilibrium points as if they are chosen as initial states the rigid body will stay at them, but any pertubation will lead to an scape from them. This fact, which implies that the regulation point is attained unless those isolated points are initial states, is proved for one point and for the other five analogous arguments can be used. So, choosing the start point represented by P1 : d1 = e2, d2 = -e1, d3 = -e3 and di = 0 (w = 0) yelds a08img06.gif (928 bytes)(t). Now, observing the form of V, it is straightforward to conclude that in the neighbourhood U of the rest configuration P1 one has VP1 > VP, "P Î U, which, together with the fact that a08img06.gif (928 bytes) is a not an incresing function, implies that the rigid body will be repelled as it was claimed before.


The Tracking Problem

In order to design a feedback control law for a tracking problem the following scalar error tracking s is introduced.

a08frm25.gif (2539 bytes).


a08frm26.gif (3473 bytes),

where Ei and VEi are, respectively, the desired rotational trajectory and its velocity.

For the tracking problem, the control law will have two main components: the feedback and the feedforward. Roughly speaking, the feedforward component is responsible for providing the necessary signal to make the required motion and the feedback component mff  is used to cancel the effects of disturbances. Indeed, the feedback stabilizes the tracking error dynamics.

The feedforward control input is determined by the following set of equations.

a08frm27.gif (2094 bytes) (15)

a08frm28.gif (1732 bytes) (16)

where a08img02.gif (588 bytes) and a08img03.gif (541 bytes) are, respectively, the acceleration and velocity reference input signals. Please note that they are not necessarily the angular acceleration and velocity of the desired trajectory. Indeed, a08img02.gif (588 bytes) is defined through the expression:

a08frm29.gif (2733 bytes).

Considering the equations (15), (16), (11) and (12) the following error system of equations is obtained:

a08frm30.gif (3054 bytes). (17)

a08frm31.gif (2113 bytes) (18)

The feedback component of the control is designed inspired in the regulation problem, so the total control input is given by

a08frm32.gif (2252 bytes), (19)

where fi = Kp SiKv VSi

The next step is to demonstrate that the closed-loop system resulting from the use of the above control law is stable and tracks the desired trajectory. For this purpose a Lyapunov like function is introduced below:

a08frm34.gif (4128 bytes)

which can be easily shown to be positive, "t Its time derivative is given by:

a08frm35.gif (6063 bytes). (20)

Now, substituting (17), (18) and (19) in (20), one has:

a08frm36.gif (4552 bytes). (21)

A crucial detail in the present development is to recognize that a08img04.gif (783 bytes), in which W denotes the skew-symmetric matrix associated to the angular velocity, is a skew-symmetric matrix. The demonstration of this this statement is presented in Appendix A.

So, using in (21 the fact that a08img04.gif (783 bytes) is a skew-symmetric operator yields:

a08frm37.gif (6928 bytes),


a08frm38.gif (3134 bytes)

and finally

a08frm39.gif (3289 bytes). (22)

The above result is analogous to that one appearing in the regulation problem, which means that for each trajectory of the rigid body the function V will not increase. Thus, as V(t) ³ 0 and a08img06.gif (928 bytes)(t) £ 0 , "t Î [ 0 ,¥ ), given any initial state represented by {di(0), di(0)} in which V(0) ¹ 0 the rigid body will track a trajectory in which a08img06.gif (928 bytes)(t) = 0. So, the central issue is to understand if the tracked trajectory is the desired one, e.g. V(t) = 0, "t . In order to understand this situation, the first step consists in recognizing that from a08img06.gif (928 bytes)(t) = 0 one has: a08img07.gif (993 bytes) = VEi, which in turn implies w = a08img08.gif (1036 bytes) and a = a08img02.gif (588 bytes). Thus, from (17).

a08frm40.gif (2509 bytes)

which has some solutions, as it was considered in the analysis of the regulation problem. By the way, any solution different of the desired trajectory (e.g. di ¹ Ei) implies that w = 0 Þ VEi = (0, 0, 0), which often is not coherent with the velocity of the desired trajectory. So, only in the case that the desired trajectory tends to a steady position, those points not coincident with the goal could be reached by the closed-loop system. Using similar arguments as those used in the case of the regulation problem, one shows that those points are not attractors and thus cannot be reached by a continuous path. Finally, this means that the presented control law is globally asymptotically convergent to the desired trajectory.


Numerical Simulations

In the present section two numerical simulations are presented. Both simulations deals with the heavy top problem (Simo and Wong, 1991), which is a purely rotational situation. The first simulation consists in a regulation problem and the second in a tracking problem. In both cases a fully actuated system is supposed, which means that an independent actuation is associated for each degree of freedom. The PD laws introduced in the previous item are used with an additional term that accounts for the gravity compensation.


Regulation Simulation

In the present simulation, the set point is defined by rest position and the inertial frame, e.g. di = ei and di = (0, 0, 0). The initial conditions are given by:

di (0) = exp ([0.3 0 0]) eI,

w (0) = 30 d3 (0).

where exp is the exponential operator (Simo and Wong, 1991), which maps onto SO(3). Besides, the following parameters were adopted: Mg = 20, l = 1, I11 = I22 = 5 and I33 = 1, where Mg is the norm of the gravitational force and the parameter I denotes the distance between the bottom of the top and its center of gravity.

In figures 1 and 2, the norm of the tracking error defined in (13) and the angular velocity associated with the rigid body are plotted against time. Those figures show the asymptoticaly convergence to the regulation set point.




In order to get some insight in the role played by the gains Kp and Kv , the norm of the applied moment a, the norm of the tracking error and the energy of the closed-loop system are depicted for two different values of the gains in figures 3, 4 and 5.





The behavior of the proposed PD law is tested in the neighborhood of one equilibrium points not coincident with the regulation set point. The initial conditions were choosen as:

di(0) = exp([0, 0, 0.9p]) ; w(0) = 0,

which are near the point P1 introduced in the previous item. As it can be observed in figure 6, the systems tends toward the desired goal and not to the equilibrium point given by P1.



Tracking Simulation

The same heavy top is treated here in a tracking context. The desired trajectory is obtained by the superposition of two rotations, namely:

Qd = Q1Q2,

where Q2 consists of a spinning around (d3) and Q1 is a rotation about a fixed axis (e1), e.g.:

Q2 = exp([0, 0, a t]) ; Q1 = exp ([b t, 0, 0]),

with a and b two constant positive scalars. From the prescribed trajectory given above, the angular velocity and acceleration are computed as:

wd = a d3 + be1.,

ad = a a08img05.gif (613 bytes) = a b (e1 d3).

In order to test the proposed control law, two situations were analysed. The first one has its initial condition coincident to those prescribed by the desired trajectory, which implies that only the feedforward component was active resulting in a very small value of the tracking error. In the second situation, the initial conditions of the top are taken slightly different (di = ei; w(0) = (b, 0, 1.01a) in order to assess the feedback capabilities of the PD law. As it is observed, for the second situation, in figures 7 and 8 the tracking is achieved after a while. For the first situation, where no pertubations were introduced, the body had followed exactly the prescribed trajectory.




As it is seen in the results of the numerical experiment depicted in the figure 8, the stability provided by the control law is achieved even for the long term dynamics as it was theoretically predicted. A lot of other experiments involving different kind of pertubations, like, for instance, in the mass, geometry or gravity, were performed and the same asymptotic behavior was found.


Concluding Remarks

In the present work, a proportional-derivative feedback law for controling the movement of rigid bodies in the context of multibody dynamics is presented and explored. The main characterisct of the propose law is the way it deals with the description of the large rotations by means of directors. It is very well known that large rotations, at least in the computational domain, can give rise to problems like singularities or strong nonlinear problems, which implies in serious drawbacks for the control implementation.

A mathematical proof of stability provided by the control law is presented based on the Lyapunov's formalism and this proof is numerically assessed. No systematic ways for determining the optimum control parameters were introduced. However, the intensive use of numerical simulation has showed up as a practical tool so to obtain those parameters.


Appendix A:

In the present appendix the proposition a08img04.gif (783 bytes), in which W denotes the skew-symmetric matrix associated to the angular velocity, is a skew-symmetric matrix

Dem.: Taking the time derivative of (9) yields

a08frm41.gif (3091 bytes)

which can be rewritten with the aid of (5) as:

a08frm42.gif (3724 bytes) (23)

where wk = w. dk.

Also from (5) one has:

a08frm43.gif (3646 bytes)

leading to:

a08frm44.gif (4143 bytes). (24)

Thus, from (23) and (24):

a08frm45.gif (5709 bytes). (25)

Remembering in the above equation that ekji = -eijk and ekji = eijk one cocludes that:

a08frm46.gif (4795 bytes)

proving the initial claim.



Bullo, F. and Murray, R.M., 1997, ''Tracking for Fully Actuated Mechanical Systems: A Geometric Framework'', CIT-CDS Technical report 97-001 - California Institute of Technology (Pasadena).        [ Links ]

Campos, A.D. and Rochinha, F.A, 1996, ''Dynamics of Rods: Numerical Formulation'', Proceedings of Advances in computational Techniques for Structural Engineering, Budapest, Hungary, pp.239-246.        [ Links ]

Fjellstad, O. and Fossen, T.I., 1994, ''Singularity-Free Tracking of Unmanned Underwater Vehicles in 6 DOF '', Proceedings of the 33rd IEEE Conference on Decision and Control, Florida.        [ Links ]

Garcia de Jalón, J. and Bayo, E., 1993, ''Kinematic and Dynamic Simulation of Multibody Systems'', Springer Verlag, New York.        [ Links ]

Goldstein, H., 1980, Classical Mechanics, Second Edition, Addison-Wesley, Reading, MA.        [ Links ]

Hughes, P.C. 1986, ''Spacecraft Attitude Dynamics'', Wiley, New York.        [ Links ]

Junkins, J.L. and Kim, Y., 1993, ''Introduction to Dynamics and Control of Flexible Structures'', AIIA - Education Series, Washington.        [ Links ]

Kelkar, A. and Joshi, S., 1996, ''Control of Nonlinear Multibody Flexible Space Structure '', Lecture Notes in Control and Information Science, 221, Springer, London.        [ Links ]

Le Tallec, P., Mani,S.A. and Rochinha,F.A., 1992, ''Finite element computation of hyperelastic rods in large displacements'', RAIRO Math. Model. Num. Anal. 26, pp. 325-346.        [ Links ]

Rochinha, F.A.,1990, ''Modeling and numerical simulation of rods'' (in portuguese), D.Sc. Thesis, PUC-Rio, Rio de Janeiro, Brazil.        [ Links ]

Rochinha, F.A. and Sampaio, R., 1996, ''Numerical Simulation of Multibody Systems: The Rigid Body Problem '', Proceedings of the Second ECCOMAS Conference on Numerical Methods in Enginnering, Paris, França, pp. 605-610.        [ Links ]

Rochinha, F.A. and Sampaio, R., 2000, ''Non-Linear Dynamics of Rigid Body: Energy and Momentum Conserving Algorithm, ''Computer Modeling in Engineering & Science. 1, no 2, pp. 7-18.        [ Links ]

Simo,J.C. and Wong,K.K., 1991, ''Unconditionally stable algorithms for rigid body dynamics that exactly preserve energy and momentum'', Int. J. Numer. Methods Eng. 31, pp. 19--52.        [ Links ]

Slotine, J.J.E. and Li,W., 1991, ''Applied Nonlinear Control'', Prentice Hall , Englewood Cliffs, New Jersey.        [ Links ]



Manuscript received: September 1999. Technical Editor: Átila P.Silva Freire.