THE PUBLISHING HOUSE   PROCEEDINGS OF THE ROMANIAN ACADEMY, Series A,  OF THE ROMANIAN ACADEMY   Volume 13, Number 3/2012, pp. 235 – 242  INTERNAL JOINT FORCES IN DYNAMICS OF A 3- PRP   PLANAR PARALLEL ROBOT  Stefan STAICU   1 , Damien CHABLAT   2 1   Department of Mechanics, University “Politehnica”   of Bucharest, Romania  2   Institut de Recherche en Communications et Cybernétique de Nantes, France staicunstefan@yahoo.com   Damien.Chablat@irccyn.ec-nantes.fr  Recursive matrix relations for the complete dynamics of a 3 -PRP   planar parallel robot are established in this paper. Three identical planar legs connecting to the moving platform are located in the same vertical plane. Knowing the motion of the platform, we develop first the inverse kinematical problem and determine the positions, velocities and accelerations of the robot. Further, the inverse dynamic problem is solved using an approach based on the principle of virtual work. Finally, some graphs of simulation for the input powers of three actuators and the internal joint forces are obtained.  Key words : Kinematics, Dynamics, Planar parallel robot, Virtual work  1.   INTRODUCTION  Parallel manipulators are closed-loop mechanisms that consist of separate serial chains connecting the fixed base to the moving platform. Compared with serial manipulators, the followings are the potential advantages of parallel architectures: higher kinematical precision, lighter weight and better stiffness, greater load bearing, stabile capacity and suitable position of arrangement of actuators. Equipped with revolute or prismatic actuators, parallel   manipulators have a robust construction and can move bodies of large dimensions with high velocities and accelerations. That is reason why the devices, which produce translation or spherical motion to a platform, technologically are based on the concept of parallel manipulators [1]. Over the past two decades, parallel manipulators have received more attention from researches and industries. Important companies such as Giddings & Lewis, Ingersoll, Hexel and others have developed them as high precision machine tools. Considerable efforts have been devoted to the kinematics and dynamic analysis of fully parallel manipulators. Among these, the class of manipulators known as Stewart-Gough platform focused great attention (Stewart [2]; Merlet [3]). They are used in flight simulators and more recently for Parallel Kinematics Machines. The prototype of Delta parallel robot (Clavel [4]; Tsai and Stamper [5]) developed by Clavel at the Federal Polytechnic Institute of Lausanne and by Tsai and Stamper at the University of Maryland as well as the Star parallel manipulator (Hervé and Sparacino [6]) are equipped with three motors, which train on the mobile platform in a three-degrees-of-freedom general translation motion. Angeles [7], Wang and Gosselin [8] analysed the kinematics, dynamics and singularity loci of Agile Wrist spherical robot with three actuators. Planar parallel robots are useful for manipulating an object on a plane. A mechanism is said to be a  planar robot   if all the moving links in the mechanism perform the planar motions. For a planar mechanism, the loci of all points in all links can be drawn conveniently on a plane. In a planar linkage, the axes of all revolute joints must be normal to the plane of motion, while the direction of translation of a prismatic joint must be parallel to the plane of motion. Bonev, Zlatanov and Gosselin [9] describe several types of singular configurations by studying the direct kinematics model of a 3- RPR   planar parallel robot with actuated base joints. Aradyfio and Qiao [10] examined the inverse kinematics solution for the three different 3-DOF planar parallel robots, while Pennock and Kassner [11] present a kinematical study of a planar parallel robot where a moving platform is connected to a fixed base by three links, each leg consisting of two binary links and three parallel revolute joints. Sefrioui and Gosselin [12] give an interesting numerical solution in the inverse and direct kinematics of this kind of planar robot.
Internal joint forces in dynamics of a 3-PRP planar parallel robot   2  Using dual-number quaternion algebra, Mohammadi-Daniali et al. [13], [14] present a study of velocity relationships and singular conditions for general planar parallel robots. 2.   KINEMATICS ANALYSIS  A recursive method is introduced in the present paper, to reduce significantly the number of equations and computation operations by using a set of matrices for kinematics and complete dynamics of the 3- PRP  planar parallel robot. Having a closed-loop structure, this robot is a symmetrical mechanism composed of three planar kinematical chains with identical topology, all connecting the fixed base to the platform. The points 0 0 0   , ,   C B A   represent the summits of a triangular base and other three points define the geometry of the moving platform. Each leg consists of two links, with one revolute and two prismatic joints. The parallel mechanism with seven links ) 7 ..., , 2 , 1 , (    k T k   consists of three revolute joint and six prismatic joints (Fig. 1).  Fig. 1 The 3- PRP   planar parallel robot  In the actuation scheme   PRP   each prismatic joint is an actively controlled prismatic cylinder. Thus, all prismatic actuators can be installed on the fixed base. For the purpose of analysis, we attach a Cartesian frame ) (   0 0 0 0   T z y Ox   to the fixed base with its origin located at triangle centre O   , the 0  z   axis perpendicular to the base and the 0  x   axis pointing along the direction 0 0   B C   . Another mobile reference frame G G G   z y Gx   is attached to the moving platform. The origin of this coordinate central system is located just at the centre G   of the moving triangle. In what follows we consider that the moving platform is initially located at a   central configuration , where the platform is not rotated with respect to the fixed base and the mass centre G   is at the origin O   of the fixed frame. It is noted that the relative rotation of k  T   body with 1 ,    k k     angle must be always pointing about the direction of k  z   axis. One of three active legs (for example leg A   ) consists of a prismatic joint, which is as well as a piston   1  of mass 1  m   linked at the A A A   z y x A   1 1 1 1   frame, having a rectilinear motion of displacement A  10     . Second element of the leg is a rigid body   2   linked at the A A A   z y x A   2 2 2 2   frame, having a relative rotation about A  z   2   axis with the angle A  21     . It has the mass 2  m   and tensor of inertia 2  ˆ J   with respect to A  T 2   frame. Finally, a prismatic joint is introduced at a planar moving platform as an equilateral triangle with the edge 3 0  l l      , mass 3  m   and inertia tensor 3  ˆ J   with respect to 3  A   , which translate relatively along A  z   3   axis with the displacement A  32     (Fig. 2).  C 0  B 0  A 3  A 0  C 3  B 3  G
3   Stefan Staicu, Damien Chablat  At the central configuration, we also consider that all legs are symmetrically extended and that the angles of orientation of three edges of fixed platform are given by 3 , , 3                C B A  .   (1) In the study of the kinematics of robot manipulators, we are interested in deriving a matrix equation relating the location of an arbitrary body to the joint variables. We call the matrix   1 ,    k k  a   , for example, the orthogonal transformation 3 3      matrix of relative rotation with the angle A k k   1 ,        of link A k  T   around A k  z   axis. Starting from the origin O   and pursuing three independent legs 3 2 1 0   A A A OA   , 3 2 1 0   B B B OB   , 3 2 1 0   C C C OC   , we obtain the following transformation matrices [15] ) , , ( ) , , ( , ,   2 1 32 1 21 21 1 10   C B A i c b a q q q q q   T i                 ,   (2) where              1 0 0 0 cos sin 0 sin cos  i i i i i           ,                0 0 1 0 1 0 1 0 0  1     ,                2 0 0 0 1 3 0 3 1 2 1  2                         1 0 0 0 cos sin 0 sin cos ) (   1 , 1 , 1 , 1 , 1 , 1 ,  i k k i k k i k k i k k i k k k k   rot q          ,          k s s k s k k   q q  1 , 1 0 ) 3 , 2 , 1 (    k   .   (3)  Fig. 2 Kinematical scheme of first leg A   of the mechanism  In the inverse geometric problem, it can be considered that the position of the mechanism is completely given through the coordinates G G   y x   0 0   ,   of the mass centre G   of the moving platform and the orientation angle    of the mobile central frame G G G   z y Gx   . The orthogonal known rotation matrix of the platform from 0 0 0   z y Ox   to G G G   z y Gx   reference system is ) (   rot R      . We suppose that the position vector of G   centre T G G G   y x r   ] 0 [   0 0 0       and the orientation angle    , which are expressed by following analytical functions  A 0  y 0  A 2  z 3 A  G  x 2 A  O  A 3  α A  A 1  x 0  y 1 A  y 3 A  z 1 A  y 2 A  v 10 A   v 32 A  x G  y G  2 1  ω 21 A
Internal joint forces in dynamics of a 3-PRP planar parallel robot   4 t y y x x  G G G G  3 cos 1 * * 0 0 * 0 0            (4) can describe the general absolute motion of the moving platform in its   vertical plane . From the conditions concerning the orientation of the platform ) , , ( , 30 30   c b a q R q q   T        (5) with ) , , ( , ,   2 1 30 10 21 32 30   C B A i q q q q q   i              ,   (6) we obtain the first following relations between angles            C B A  21 21 21  .   (7) Six independent variables A A  32 10   ,      , B B  32 10   ,      , C C  32 10   ,      will be determined by several vector-loop equations as follows G k Gi T i k , k T k i   r r q r q r   0 2 1 3 30 1 0 10                 ) , , ( ) , , (   C B A i c b a q      ,   (8) where 3 10 10 0 00 10   ) 3 / (   u q l r r   T i i i            , T Gi T i i i   l r u q r r   ] 3 / 1 1 0 [ 5 . 0 , , 0   0 3 3 32 32 32 21                T A   l r   ] 0 1 0 [ 0 00       T C T B   l r l r   ] 0 1 3 [ 5 . 0 , ] 0 1 3 [ 5 . 0   0 00 0 00                                                          0 0 0 0 0 1 0 1 0 ~ , 1 0 0 , 0 1 0 , 0 0 1  3 3 2 1   u u u u       .   (9) Actually, these vector equations mean that ) sin( 2 ) cos( 3 2 ) 3 cos( cos ) 3 (   0 0 00 0 32 10 0  i i i G i i i i   l l x x l                        ) cos( 2 ) sin( 3 2 ) 3 sin( sin ) 3 (   0 0 00 0 32 10 0  i i i G i i i i   l l y y l                        ) , , (   C B A i     .   (10) The rotations of the compounding elements of each leg (for example the leg   A ) are characterized by recursive relations of following skew-symmetric matrices A k k A k k A k k T k k A k k k A k   u a a   1 , 1 , 3 1 , 1 , 0 , 1 1 , 0   , ~ ~ ~                       , ) 3 , 2 , 1 (    k   (11) which are   associated   to the absolute angular velocities 0 10        A    , 3 21   u  A              , 0 32        A     , 3 20   u  A              , 3 30   u  A              .   (12) Following relations give the velocities A k v   0     of the joints k  A 3 10 10   u v   A A          , 0 21        A  v   , 3 32 32   u v   A A           , 3 1 , 1 , 0 , 1 1 , 0 , 1 1 , 0  ~   u v r a v a v   A k k A k k A k k k A k k k A k                     .   (13) Equations of geometrical constraints (7) and (8) when differentiated with respect to time lead to the following   matrix conditions of connectivity   [16]  3 10 10   u a u v   T T j A     G T j T T j A   r u u a u v   0 3 30 32          } { ~   3 32 32 3 20  GA T A T T j   r a r u a u            ) 2 , 1 (    j  .   (14) From these equations, we obtain the relative velocities A A A   v v   32 21 10   , ,     as functions of angular velocity of the platform and velocity of mass centre G   . The Jacobian matrix given by these conditions of constraints is a fundamental element for the analysis of the robot workspace and the configurations of singularities where the manipulator becomes uncontrollable [17]. Concerning the first leg A   , the characteristic   virtual velocities   are expressed as functions of the pose of the mechanism by the general kinematical equations (14), where we add the contributions of successive virtual translations during two orthogonal fictitious displacements of the revolute joint 2  A   , as follows: ) 2 , 1 ( } { ~   0 3 30 32 3 32 32 3 20 21 3 10 21 2 10 21 3 10 10            j v u u a u v r a r u a u u a v u a u v u a u v   Gv T j T T j Av GA T A T T j Av T Azv T T j Ayv T T j Av                  . (15) Now, let us assume that the robot has successively some virtual motions determined by following sets of velocities:
5   Stefan Staicu, Damien Chablat 1 10     Av a  v  , 0 10     Bv a  v   , 0 10     Cv a  v   , 0 21     iyv a  v   , 0 21     izv a  v 0 10     iv a  v  , 1 21     Ayv a  v   , 0 21     Byv a  v   , 0 21     Cyv a  v   , 0 21     izv a  v 0 10     iv a  v  , 0 21     iyv a  v   , 1 21     Azv a  v   , 0 21     Bzv a  v   , 0 21     Czv a  v ) , , (   C B A i      .   (16) These virtual velocities are required into the computation of virtual power and virtual work of all forces applied to the component elements of the robot. As for the relative accelerations A A A  32 21 10   , ,       of the robot, new conditions of connectivity are obtained by the derivative of above equations (14):  3 10 10   u a u   T T j A       A  32       G T j T T j   r u u a u   0 3 30            } { ~ ~   3 32 32 3 3 20 2   GA T A T T j   r a r u u a u             } { ~   3 32 32 3 20  GA T A T T j   r a r u a u          3 32 3 20 32  ~ 2   u a u a u v   T T T j A         ) 2 , 1 (    j  .   (17) The following recursive relations give the angular accelerations A k   0       and the accelerations A k   0      of joints k  A 3 10 10   u  A A              , 0 21        A     , 3 32 32   u  A A               , 0 10        A     , 3 21   u  A               , 0 32        A   3 1 , 0 , 1 1 , 1 , 3 1 , 0 , 1 1 , 0  ~   u a a u a   T k k A k k k A k k A k k A k k k A k                                          T k k A k A k A k k k A k A k A k   a a   1 , 0 , 1 0 , 1 0 , 1 1 , 0 0 0  ~ ~ ~ ~ ~ ~        3 1 , 0 , 1 1 , 1 , 3 1 , 3 3 1 , 1 ,  ~ ~ 2 ~ ~ ~   u a a u u u   T k k A k k k A k k A k k A k k A k k                   ) 3 , 2 , 1 ( , ~ 2 ) ~ ~ ~ (   3 1 , 3 1 , 0 , 1 1 , 1 , 1 , 0 , 1 0 , 1 0 , 1 1 , 0 , 1 1 , 0                         k u u a a v r a a   A k k T k k A k k k A k k A k k A k A k A k k k A k k k A k                 .   (18)  3. DYNAMICS EQUATIONS  The dynamics analysis of parallel robots is complicated because the existence of a spatial kinematical structure, which possesses a large number of passive degrees of freedom, dominance of the inertial forces, frictional and gravitational components and by the problem linked to real-time control in the inverse dynamics. In the context of the real-time control, neglecting the frictions forces and considering the gravitational effects, the relevant objective of the complete dynamics is first to determine the input torques or forces, which must be exerted by the actuators in order to produce a given trajectory of the end-effector, but also to calculate all internal joint forces or torques. A lot of works have focused on the dynamics of Stewart platform. Dasgupta and Mruthyunjaya [18] used the Newton-Euler approach to develop closed-form dynamic equations of Stewart platform, considering all dynamic and gravity effects as well as viscous friction at joints. Tsai [1] presented an algorithm to solve the inverse dynamics for a Stewart platform-type using Newton-Euler equations. This commonly known approach requires computation of all constraint forces and moments between the links. Three independent mechanical systems acting along the planar directions A  z A   1 1   , B  z B   1 1   and C  z C   1 1   with the forces , 1 10 10   u f f   A A       1 10 10   u f f   B B         , 1 10 10   u f f   C C         can control the general motion of the moving platform. Knowing the position and kinematics state of each link as well as the external forces acting on the 3- PRP   planar  parallel robot, in the present paper we apply the principle of virtual work for the inverse dynamic problem in order to establish some definitive recursive matrix relations for the calculus of input torques of the actuators and internal forces in the joints. The parallel robot can artificially be transformed in a set of three open chains ) , , (   C B A i C   i      subject to the constraints. This is possible by cutting each joint for moving platform, and takes its effect into account by introducing the corresponding constraint conditions. The first and more complicated open tree system includes the first acting link and comprises also the moving platform. The wrench of two vectors   k  F     and   k  M     evaluates the influence of the action of the weight g m   k     and of other external and internal forces applied to the same element k  T   of the mechanism 2 0 *   81 . 9   u q m F   k k k       , 2 0 *   ~ 81 . 9   u q r m M   k C k k k        , ) , , (   c b a q      .   (19) Now, we compute the force of inertia in k  F     and the resulting moment of inertia forces in k  M     of an arbitrary rigid body k  T   of mass k  m   with respect to the centre of its first joint:
Internal joint forces in dynamics of a 3-PRP planar parallel robot   6 }. ˆ ~ ˆ ~ { }, ) ~ ~ ( {   0 0 0 0 0 2 0 0   k k k k k k C k k in k C k k k k k in k   J J r m M r m F                           (20) Pursuing the first leg A   , for example, two significant recursive relations generate the vectors A k T k k A k k A k T k k A k A k A k T k k A k A k   F a r M a M M F a F F   1 , 1 , 1 1 , 1 0 1 , 1 0  ~ ,                          (21) where one denoted A k inA k A k   F F F             0  , A k inA k A k   M M M             0   .   (22) Considering some independent virtual motions of the robot, all virtual displacements and velocities should be compatible with the virtual motions imposed by all kinematical constraints and joints at a given instant in time. The fundamental principle of the virtual work states that a mechanism is under dynamic equilibrium if and only if the virtual work developed by all external, internal and inertia forces vanish during any general virtual displacement, which is compatible with the constraints imposed on the mechanism. Assuming that frictional forces at the joints are negligible, the virtual work produced by all remaining forces of constraint at the joints is zero. Total virtual work contributed by the inertia forces and moments of inertia forces in k  F     , in k  M     , by the wrench of known external forces *  k  F     , *  k  M     and by the first active force A  f 10     or some internal joint forces, for example, can be written in a compact form, based on the relative virtual velocities. Applying   the fundamental equations of the parallel robots dynamics   [19], the following compact matrix relations results } {   2 21 2 21 3 32 2 21 3 1 3 10  C Cv a B Bv a A Av a A Av a T A T A   M M F v M u F u f                     (23) for the   input force   of first prismatic actuator, } {   2 21 2 21 3 32 2 21 3 2 21 2 21  C Cv B Bv A Av a A Av a T A T T A y   M M F v M u F a u f                     (24) for the first   joint force   and } {   2 21 2 21 3 32 2 21 3 2 21 3 21  C Cv B Bv A Av a A Av a T A T T A z   M M F v M u F a u f                     (25) for the second   joint force   acting in the joint 2  A   . The relations (21) - (25) represent the   complete inverse dynamics model   of the 3- PRP   planar parallel robot.  Fig. 3 Powers A  p 10   , B  p 10   , C  p 10   of three actuators   Fig. 4 Internal joint forces A y  f   21   , B y  f   21   , C y  f   21  As application let us consider a 3- PRP   planar robot, which has the following geometrical   and architectural characteristics: 3 , 3 . 0 , 12 , 025 . 0 , 0   0 0 0 0 0 * 0 * 0   l l m OC OB OA l m y x   G G                 s t kg m kg m kg m   3 , 3 , 75 . 0 , 1   3 2 1         . Assuming that there are no external forces and moments acting on the moving platform, a dynamic simulation is based on the computation of three powers i  p 10   required by eac h actuator during the platform’s
7   Stefan Staicu, Damien Chablat  evolution and the internal joint forces i z i y   f f   21 21   , ) , , (   C B A i      . Using the MATLAB software, a computer program was developed to solve the inverse dynamics of the planar parallel robot. To illustrate the algorithm, it is assumed that for a period of three second the platform starts at rest from a central configuration and rotates or moves along rectilinear directions.  If the platform’s centre G   moves along a   vertical trajectory   without rotation of platform, the input powers of three actuators and some joint forces are calculated by the program and plotted versus time as follows: Fig. 3 - Fig. 5. For the second example we consider the   rotation motion   of the moving platform about 0  z  horizontal axis with variable angular acceleration while all the other positional parameters are held equal to zero (Fig. 6 - Fig. 8).  Fig. 5 Internal joint forces A z  f   21   , B z  f   21   , C z  f   21   Fig. 6 Powers A  p 10   , B  p 10   , C  p 10   of three actuators  Fig. 7 Internal joint forces A y  f   21   , B y  f   21   , C y  f   21   Fig. 8 Internal joint forces A z  f   21   , B z  f   21   , C z  f   21  The simulation through the MATLAB program certify that one of the major advantages of the current matrix recursive formulation is a reduced number of additions or multiplications and consequently a smaller processing time of numerical computation.  4. CONCLUSIONS  Within the inverse kinematics analysis some exact relations that give in real-time the position, velocity and acceleration of each element of the parallel robot have been established in the present paper. The dynamics model takes into consideration the mass, the tensor of inertia and the action of weight and inertia force introduced by all compounding elements of the parallel mechanism.
Internal joint forces in dynamics of a 3-PRP planar parallel robot   8  Based on the principle of virtual work, this approach can eliminate all forces of internal joints and establishes a direct determination of the time-history evolution for the torques required by the actuators and the internal forces or torques in joints. Choosing appropriate serial kinematical circuits connecting many moving platforms, the present method can easily be applied in forward and inverse mechanics of various types of parallel mechanisms, complex manipulators of higher degrees of freedom and particularly   hybrid structures , when the number of components of the mechanisms is increased.  REFERENCES  1. TSAI, L.-W.,   Robot analysis: The Mechanics of Serial and Parallel Manipulators , Wiley, New York, 1999. 2. STEWART, D.,   A platform with six degrees of freedom , Proc. Inst. Mech. Eng., Part 1,   180 , 15, pp. 371-386, 1965. 3. MERLET, J.-P.,   Parallel Robots , Kluwer Academic, Dordrecht, 2000. 4. CLAVEL, R.,   Delta: a fast robot with parallel geometry , Proceedings of 18th International Symposium on Industrial Robots, Lausanne, pp. 91-100, 1988. 5. TSAI, L.-W., STAMPER, R.,   A parallel manipulator with only translational degrees of freedom , ASME Design Engineering Technical Conferences, Irvine, CA, 1996. 6. HERVÉ, J.-M., SPARACINO, F.,   Star: A new concept in robotics , Proceedings of the Third International Workshop on Advances in Robot Kinematics, Ferrara, pp. 176-183, 1992. 7. ANGELES, J.,   Fundamentals of Robotic Mechanical Systems: Theory, Methods and Algorithms , Springer, New York, 2002. 8. WANG, J., GOSSELIN, C.,   A new approach for the dynamic analysis of parallel manipulators , Multibody System Dynamics, Springer,   2 , 3, pp. 317-334, 1998. 9. BONEV, I., ZLATANOV, D.,   GOSSELIN, C.,   Singularity analysis of 3-DOF planar parallel mechanisms via screw theory , Journal of Mechanical Design,   25 , 3, pp. 573-581, 2003. 10. ARADYFIO, D.D., QIAO, D.,   Kinematic simulation of novel robotic mechanisms having closed chains , ASME Mechanisms Conference, Paper 85-DET-81, 1985. 11. PENNOCK, G.R, KASSNER, D.J.,   Kinematic analysis of a planar eight-bar linkage: application to a platform-type robot , ASME Mechanisms Conference, Paper DE - 25, pp. 37-43, 1990. 12. SEFRIOUI, J., GOSSELIN, C.,   On the quadratic nature of the singularity curves of planar three-degrees-of-freedom parallel manipulators , Mechanism and Machine Theory, Elsevier,   30 , 4, pp. 533-551, 1995. 13. MOHAMMADI-DANIALI, H., ZSOMBOR-MURRAY, P., ANGELES, J.,   The kinematics of spatial double-triangular parallel manipulators , ASME Journal of Mechanical Design, Vol.   117 , pp. 658-661, 1995. 14.   MOHAMMADI-DANIALI,   H.,   ZSOMBOR-MURRAY,   P.,   ANGELES,   J.,   Singularity   analysis   of   planar   parallel manipulators , Mechanism and Machine Theory, Elsevier,   30 , 5, pp. 665-678, 1995. 15. STAICU, S.,   Power requirement comparison in the 3-RPR planar parallel robot dynamics , Mechanism and Machine Theory, Elsevier,   44 , 5, pp. 1045-1057, 2009. 16. STAICU, S.,   Dynamics of the 6-6 Stewart parallel manipulator , Robotics and Computer-Integrated Manufacturing, Elsevier,  27 , 1, pp. 212-220, 2011. 17. PASHKEVICH, A., CHABLAT, D., WENGER, P.,   Kinematics and workspace analysis of a three-axis parallel manipulator: the Orthoglide , Robotica, Cambridge University Press,   24 , 1, 2006. 18. DASGUPTA, B., MRUTHYUNJAYA, T.S.,   A Newton-Euler formulation for the inverse dynamics of the Stewart platform manipulator , Mechanism and Machine Theory, Elsevier,   33 , pp. 1135-1152, 1998. 19. STAICU, S., LIU, X-J., LI, J.,   Explicit dynamics equations of the constrained robotic systems , Nonlinear Dynamics, Springer,  58 , 1-2, pp. 217-235, 2009.  Received January 11, 2012