Evolution of the DeNOC-based dynamic modelling for multibody systems

Abstract. Dynamic modelling of a multibody system plays very essential role in its analyses. As a result, several methods for dynamic modelling have evolved over the years that allow one to analyse multibody systems in a very efficient manner. One such method of dynamic modelling is based on the concept of the Decoupled Natural Orthogonal Complement (DeNOC) matrices. The DeNOC-based methodology for dynamics modelling, since its introduction in 1995, has been applied to a variety of multibody systems such as serial, parallel, general closed-loop, flexible, legged, cam-follower, and space robots. The methodology has also proven useful for modelling of proteins and hyper-degree-of-freedom systems like ropes, chains, etc. This paper captures the evolution of the DeNOC-based dynamic modelling applied to di fferent type of systems, and its benefits over other existing methodologies. It is shown that the DeNOC-based modelling provides deeper understanding of the dynamics of a multibody system. The power of the DeNOC-based modelling has been illustrated using several numerical examples.


Introduction
Over the last two decades, applications of multibody dynamics have expanded over the fields of robotics, automobile, aerospace, bio-mechanics, and many others.With continuous development in the above mentioned fields, many complex multibody systems have evolved whose dynamics play a pivotal role in their behaviour.Hence, computer-aided dynamic analysis of multibody systems has been a prime motive to the engineers, as high speed computing facilities are readily available.In order to perform computer-aided dynamic analysis, the actual system is represented with its dynamic model which has the information of its link parameters, joint variables and constraints.The dynamic model is nothing but the equations of motion of the multibody system at hand derived from the physical laws of motions.For a system with fewer links, it is easier to obtain explicit expressions for the equations of motion.However, finding equations of motion for complex systems with many links is not an easy task.Sometimes even with 4 or 5 links, say, a 4-bar mechanism, it is difficult to find an explicit expression for the system's inertia in terms of its link lengths, masses, and joint angles.Hence, development of the equations of motion is an essential step for the dynamic analysis.
There are several fundamental methods for the formulation of equations of motion (Greenwood, 1988).For example, Newton-Euler (NE) formulation, Euler-Lagrange principle, Gibbs-Appel approach, Kane's method, D'Alembert's principle, and similar others.All the above mentioned approaches when applied to multibody systems have their own advantages and disadvantages.For example, NE approach, which is one of the classical methods for dynamic formulation, uses the concept of "free-body diagrams".For coupled systems, constrained forces (which are meant here to include both forces and moments) along with those applied externally are included in the free-body diagrams.Mathematically, the NE equations of motion lead to three translational equations of motion of the Centre-of-Mass (COM), and three equations determining the rotational motion of the rigid body.The NE equations of any two free bodies are related through the constraint forces acting at their interface.The constraint forces arise due to the presence of a kinematic pair, e.g., a revolute or a prismatic, between the two neighbouring bodies.For an open-loop multibody system, these constraints along with other unknowns, i.e., the actuating forces can be easily solved recursively.However, for a closed-loop system, the NE equations generally need to be solved simultaneously in order to obtain the driving and constraint forces together.Hence, the use of the NE equations of motion for closed-loop systems is not as efficient as those for open-loop systems.
Euler-Lagrange (EL) formulation is another classical approach which is widely used for dynamic modelling.The EL formulation uses the concept of generalized coordinates instead of Cartesian coordinates.It is based on the minimization of a functional called "Lagrangian" which is nothing but the difference between kinetic energy and potential energy of the system at hand.For open-loop multibody systems, where typically the number of generalized coordinates equals the degree-of-freedom of a system, the constraint forces do not appear in the equations of motion.For closed-loop multibody systems, however, the forces of constraints appear as Lagrange's multipliers.
Kane's formulation (Kane and Levinson, 1983), which is same as the Lagrange's form of D'Alembert's principle, has also been used by many researchers for the development of equations of motion.It is found to be more beneficial than other formulations when used for systems with nonholonomic constraints.Several other methods of dynamic formulations were also proposed in the literature.For example, Khatib (1987) presented the operational-space formulation, whereas Angeles and Lee (1988) presented the natural orthogonal complement (NOC) based approach.Blajer et al. (1994) have also presented an orthogonal complement based formulation for the constrained multibody systems.Park et al. (1995) presented robot dynamics using a Lie group formulation, while Stokes and Brockett (1996) derived the equations of the motion of a kinematic chain using concepts associated with the special Euclidean group.McPhee (1996) showed how to use linear graph theory in multibody system dynamics.Cameron and Book (1997) described a technique based on Boltzmann-Hamel equations to derive dynamic equations of motion.Comprehensive discussion on dynamic formalisms can be found in the seminal text by Roberson and Schwertassek (1988), Schiehlen (1990Schiehlen ( , 1997)), Shabana (2001), and Wittenburg (2008).Recent trends in dynamic formalisms can also be found in the work by Eberhard and Schiehlen (2006).

Natural Orthogonal Complement (NOC)
It is pointed out here that the Newton-Euler (NE) equations of motion are still found to be popular in the literature of dynamic modelling and analyses.However, it requires solution of the constraint forces which do not play any role in the motion of a system.Hence, extra calculations are required in motion studies.To avoid such extra calculations, there are formulations proposed in the literature where the equations of motion in the Euler-Lagrange (EL) form are obtained from the NE equations.Huston and Passerello (1974) were first to introduce a computer oriented method to reduce the dimension of the unconstrained NE equations by eliminating the constraint forces.Later, Kim and Vanderploeg (1986) derived the equations of motion in terms of relative joint coordinates from Cartesian coordinates through the use of velocity transformation matrix.Velocity transformation matrix relates linear and angular velocities of the links with joint velocities.It is worth noting here that the vector of constraint forces is orthogonal to the columns of the velocity transformation matrix.More precisely, the columns of the velocity transformation matrix span the nullspace of the matrix of velocity constraints.Hence, the said velocity transformation matrix is also referred to as an "orthogonal complement matrix".The phrase "orthogonal complement" was first coined by Hemami and Weimer (1981) for the modelling of nonholonomic systems.Orthogonal complements are not unique.In some approaches, it was obtained numerically, e.g., using singular value decomposition or treating it as an eigen value problem (Wehage and Haug, 1982;Kamman andHuston, 1984, Mani et al., 1985), which are computationally inefficient.
Alternatively, Angeles and Lee (1988) presented a methodology where they derived an orthogonal complement naturally from the velocity constraints.Hence, the name Natural Orthogonal Complement (NOC) was attached to their methodology.The NOC matrix, when combined with the NE equations of motion, leads to the minimal-order constrained dynamic equations of motion by eliminating the constraint forces.This facilitates the representation of the equations of motion in Kane's form that is suitable for recursive computation in inverse dynamics or in the EL form that is suitable for forward dynamics and integration.Later, Angeles andMa (1988), Cyril (1988), Angeles et al. (1989), and Saha and Angeles (1991) showed the effectiveness of the use of the NOC matrix while applied to systems with holonomic and nonholonomic constraints.

The Decoupled NOC (DeNOC)
Subsequently, Saha (1995Saha ( , 1997) ) presented the decoupled form of the NOC for the serial multibody systems.The two resulting block matrices, namely, an upper block triangular and a block diagonal matrices, are referred to as the Decoupled NOC (DeNOC) matrices.In contrast to the NOC, the DeNOC matrices allow one to recursively obtain the analytical expressions of the vectors and matrices appearing in the equations of motion (Saha, 1999a).This in turn helps to analytically decompose the Generalized Inertia Matrix (GIM) arising out of the constrained equations of motion of the system at hand, allowing one to obtain analytical inverse of the GIM (Saha, 1999b) and a recursive algorithm for forward dynamics (Saha, 2003).Later, Saha and Schiehlen (2001) showed the power of the DeNOC matrices in obtaining recursive algorithms for the dynamics analyses of closed-loop parallel systems.Subsequently, Khan et al. (2005) illustrated the effectiveness of the DeNOC-based methodology in modelling parallel manipulators.Inspired by the concept of the DeNOC matrices, Dimitrov (2005) used a similar method for dynamic analysis, trajectory planning, and control of space robots.Garcia de Jalon et al. (2005) have also derived matrices which they have pointed out to be similar to the De-NOC matrices of Saha (1995Saha ( , 1997)).The DeNOC matrices have also found an application in the architecture design of a manipulator through its dynamic model simplifications (Saha et al., 2006).More recently, Chaudhary and Saha (2007) have applied the concept of the DeNOC matrices for the dynamic analyses of general closed-loop systems.They have also introduced the concepts like "determinate" and "indeterminate" subsystems which helped to achieve subsystem-level recursions for the inverse dynamics of a general closed-loop system.Systems with closedloops which are used in automobile steering systems were analyzed by Hanzaki et al. (2009), whereas fuel injection pumps of diesel engines with rolling contacts were analyzed by Sundarranan et al. (2012).Extending the concept of the DeNOC matrices to other type of systems, Mohan and Saha (2007) showed how to derive the DeNOC matrices for a rigid-flexible multibody system.The methodology not only provided efficient dynamic algorithms but also produced numerically stable results.Very recently, Shah et al. (2012a) introduced a concept of "kinematic module" to a tree-type multibody system and derived module-level De-NOC matrices, which provided macroscopic purview of the multibody systems.Moreover, intra-and inter-modular recursive algorithms were derived for the analyses and control of legged robots (Shah, 2011;Shah et al., 2013).It was shown that the concept of Euler-angle-joints (EAJs) (Shah et al., 2012b) coupled with the module-level DeNOC matrices provided very efficient dynamic algorithms for the multibody system consisting of multiple branches and multipledegrees-of-freedom joints.The algorithms have been implemented in a free software called ReDySim (acronym for Recursive Dynamic Simulator), which can be downloaded free from http://www.redysim.co.nr.ReDySim can be easily used by the students and researchers of multibody dynamics.Note here that the DeNOC-based algorithm was also used by the researchers from other domain, e.g., Patriciu et al. (2004) have adopted the concept for the analysis of conformational dependence of mass-metric tensor determinants in serial polymers with constraints.
The main motivation behind this paper is to bring forth the developments of the DeNOC-based dynamic modelling for multibody systems, which have taken place over more than one and half decades.The paper explains the fundamental principles of the DeNOC-based formulation, their benefits and applications.Rest of the paper is organized as follows: Sect. 2 presents the DeNOC-based dynamic modelling for serial-chain systems, which forms the basis for the dynamic modelling of other type of systems, e.g., tree-type systems explained in Sect.3. Application to closed-loop systems is explained in Sect.4, whereas two software, namely, Robo-Analyzer and ReDySim, developed for the use by the students and researchers of multibody dynamics are explained in Sect. 5.The computational aspects are provided in Sect.6.Finally, conclusions are given in Sect.7.

DeNOC-based dynamic modelling for serial-chain systems
The Natural Orthogonal Complement (NOC) matrix proposed by Angeles and Lee (1988) relates the angular and linear velocities of the rigid bodies in a mechanical system to its associated joint-rates.It is used to develop a set of independent equations of motion from the unconstrained or uncoupled Newton-Euler (NE) equations using free-body diagrams.These independent set of equations was referred by the authors as the Euler-Lagrange equations of motion.Unlike the NOC, its decoupled form, i.e., the DeNOC, proposed by Saha (1995Saha ( , 1997)), allows one to write the expressions of each element of the matrices and vectors associated with the dynamic equations of motion in analytical recursive form.

Preliminaries and notation
An open-loop serial-chain system, e.g., a robotic manipulator shown in Fig. 1, has a fixed-base, denoted by #0, and n moving rigid bodies or links, indicated with #1, ..., #n, coupled by n single degree-of-freedom (DOF) kinematic pairs or joints numbered as 1, ..., n.The joints are generally revolute or prismatic.In presence of higher-DOF joints, they are modelled as combinations of single-DOF joints.For example, a spherical joint can be modelled as three intersecting revolute joints, whereas a cylindrical joint is modelled as a combination of revolute and prismatic joints.Few terms are defined below which will be used throughout the paper for the derivation of the dynamic models.The 6-dimensional vectors, twist (t i ) of the i-th rigid link undergoing motion in the 3-dimensional Cartesian space and wrench (w i ), acting on the i-th link are defined by: where ω i is the 3-dimensional vector of angular velocity, and v i is the 3-dimensional vector of linear velocity of the mass center (C i ) of the i-th link, whereas n i and f i are the 3dimensional vectors of the moment and force applied about and at C i , respectively.The 6 × 6 matrices of mass M i , and angular velocity W i , of the i-th body are represented by: Endeffector Base, #0 where ω i × 1 is the 3 × 3 cross-product tensor associated with the angular velocity vector ω i which when operates on any 3-dimensional Cartesian vector x leads to the cross-product vector between ω i and x, i.e., (ω i × 1) x ≡ ω i × x.Also, 1 and O are the 3 × 3 identity and zero matrices, respectively, whereas I i and m i are the 3 × 3 inertia tensor about C i , and the mass of the i-th link, respectively.For the serial-chain mechanical system shown in Fig. l, the method to obtain the dynamic equations of motion using the DeNOC matrices is as follows: -Derive the DeNOC matrices.
-Obtain the unconstrained NE equations of motion from the free-body diagrams of each link, and -Couple the DeNOC matrices with the unconstrained NE equations to obtain a set of constrained independent equations of motion which are same as the system's EL equations of motion.
The above steps are explained next in the following subsections.

Kinematic constraints
The kinematic constraints in terms of the velocities of two neighbouring links, say, #i and #j, coupled by a revolute joint, as shown in Fig. 2, are given by where ω j and v j are the angular velocity and velocity of the mass of link j, i.e., C j , respectively.Similarly, ω i and v i are defined for the neighbouring link i, whereas θi is the jointrate of the i-th joint.The above six scalar equations can be written in a compact form as where B i j is the 6 × 6 matrix and p i is the 6-dimensional vector which are given by Here, c i j is the 3-dimensional position vector from C i to C j given by c i j ≡ −d i − r j , and c i j × 1 is the cross-product tensor associated with vector c i j .It is defined similar to ω i × 1 of Eq. (2).Moreover, e i is the unit vector parallel to the axis of rotation of the i-th revolute joint.Interestingly, matrix B i j and vector p i have the following interpretations: -If links #i and # j are rigidly attached, B i j propagates twist or velocities of # j to #i.Hence, B i j is termed in Saha (1999a) as the twist-propagation matrix, which satisfies -On the other hand the vector p i takes into account the motion of the i-th joint.Hence, vector p i is termed as the joint-rate-propagation vector.The vector p i in Eq. ( 5) is defined for a revolute joint.For a prismatic joint, it is given by Equation ( 4) can be written for i = 1, ..., n, as where 1 is the 6n × 6n identity matrix, and the 6n × 6n matrix B has the following representation: It is now simple matter to invert the 6n × 6n matrix, (1 − B), and hence, Eq. ( 8a) can be rewritten as In Eq. ( 9a), the matrix N is the 6n × n Natural Orthogonal Complement (NOC) matrix, as introduced by Angeles and Lee (1988), whereas N l and N d are the decoupled form of the NOC or the DeNOC matrices proposed first time in Saha (1995).The 6n × 6n matrix N l and the 6n × n matrix N d are given by and Note that in Eq. (9b), N l is a lower block-triangular matrix, whereas N d is a block-diagonal matrix, as indicated through their subscripts "l" and "d", respectively.Moreover, O and 0 are the 6 × 6 matrix of zeros and the 6-dimensional vector of zeros, respectively.The n-dimensional vector θ is defined as which contains the joint-rates of all the joints in the serialchain system shown in Fig. 1.

Unconstrained Newton-Euler (NE) equations
The unconstrained or uncoupled Newton-Euler (NE) equations of motion for the i-th rigid-link (Saha, 1999a) can be written from its free-body diagram, Fig. 3, as where ωi and vi are the angular acceleration and acceleration of the mass center C i , respectively.Moreover, I i is the 3 × 3 inertia tensor of i-th link about its mass center C i , and m i is its mass.Other variables were defined after Eq. ( 1).The above six scalar equations can be put in a compact form as Revolute Joint 5 Revolute Joint 6 Gravity where t i , w i and W i , M i are defined in Eqs. ( 1) and ( 2), respectively.Moreover, ṫi is the time derivative of the twist t i of the i-th link.For the whole system of n rigid links, the 6n scalar equations (for i = 1, ..., n, where n is the number of moving rigid links in the serial chain system) can be written as In Eq. ( 13), ṫ is the time derivative of the generalized twist, t.Moreover, M and W are the 6n × 6n generalized mass matrix and generalized matrix of angular velocities, respectively, i.e., Moreover, w and t are the 6n-dimensional vectors of generalized wrench and twist, respectively.They are defined as

Constrained equations using the DeNOC matrices
The kinematic constraints in velocities, i.e., Eq. (9a), then can be incorporated into the unconstrained NE equations of motion, Eq. ( 13).This is done by pre-multiplying N T with the 6n unconstrained NE equations of motions of Eq. ( 13), i.e., where w is substituted as, w ≡ w E + w C , in which w E and w C are the 6n-dimensional vectors of external and constraint wrenches, respectively.Since the constraint wrenches do www.mech-sci.net/4/1/2013/Mech.Sci., 4, 1-20, 2013 not do any work, N T w C vanishes (Angeles and Lee, 1988).Hence, N T w C = 0. Substituting the expression of t from Eq. ( 9a) and its time derivative, ṫ = N T θ + Ṅθ into Eq.( 16), one can get the n independent scalar dynamic equations of motion, namely, where, I ≡ N T MN: the n × n generalized inertia matrix (GIM); C ≡ N T (MN + WMN): the n × n matrix of convective inertia terms (MCI); and τ ≡ N T w E : the n-dimensional vector of generalized forces of driving, and those resulting from gravity, dissipation, and other external forces like footground interaction of a walking robot, etc., if any.

Analytical expression of the GIM
The analytical expression of the generalized inertia matrix (GIM) appearing in Eq. ( 17) plays an important role in simplifying, mainly, the forward dynamics agorithm (Saha, 1999a(Saha, , 2003)).In this section, the GIM I is derived using the expressions of the DeNOC matrices (Saha, 1995(Saha, , 1997(Saha, , 1999a(Saha, , b, 2003)).Substituting the expressions of the DeNOC matrices given by Eq. ( 9b) into the expression of the GIM appearing after Eq. ( 17), one gets The 6n × 6n symmetric matrix M can be written as where the 6 × 6 matrix, Mi , for i = 1, • • • , n, can be obtained recursively, i.e., in which Mi+1 ≡ O, because there is no (n+1)st link in the serial-chain.Hence, Mn ≡ M n .The matrix, Mi , is interpreted as the mass matrix of the Composite Body, i, that consists of rigidly connected links #i, ..., #n, as indicated in Fig. 1.Finally, the n × n GIM I can be expressed as , where i i j ≡ p T i Mi B i j p j (21) for i = 1, ..., n; j = 1, ..., i.The term i i j is a scalar and "sym" denotes symmetric elements of the GIM I.

Recursive inverse dynamics algorithm
The inverse dynamics of a serial-chain system is defined as the process of determining the joint forces/torques when the joint motions of the system are known.The inverse dynamics algorithm calculates the joint torque, τ i , for i = 1, ..., n, in two recursive steps, namely, forward and backward recursions.They are given below.

Step 1: forward recursion
First, the 6-dimensional twist and twist-rate vectors of each link, i.e., t i and ṫi , respectively, are calculated, for i =1, ..., n, using the following relations: In the above equations, t 0 = 0 and ṫ0 = 0, as link #0 is fixed without any motion.

Step 2: backward recursion
The 6-dimensional vector, wi , and the scalar, τ i , for i = n, ..., 1, are calculated using the following relations: wi = B T i+1,i wi+1 , and where for i = n, wn+1 = 0, as there is no (n + 1)st link in the system.Hence, wn = w n .The effect of gravity can also be taken into account by providing negative acceleration due to gravity, g, to the twist-rate of the first link as an additional term (Kane and Levinson, 1983), i.e., ṫ1 = p 1 θ1 + ṗ1 θ1 + ρ, where ρ ≡ 0 T , − g T (26) Note that Eqs. ( 22)-( 26) were reported in Saha (1999a) with different notations, which actually have the same interpretations as given above, i.e., twist (t i ), twist-rate ( ṫi ), wrench of composite body ( wi ), etc.Based on the above mentioned recursive inverse dynamics algorithm, a computer program was developed in C++ which was called RIDIM (Recursive Inverse Dynamic for Industrial Manipulators) (Saha, 1999a).
Recently, a similar algorithm has been rewritten in Visual C# and implemented in the "IDyn" module of the newly developed software called RoboAnalyzer (Rajeevlochana and Saha, 2011;Rajeevlochana et al., 2012) which also has 3dimensional visualisation of the system under study.It is explained in Sect.6.1, and available free from http://www.roboanalyzer.comfor the benefits of students and researchers of multibody dynamics community.

Recursive forward dynamics algorithm
Forward dynamics of a serial-chain system is defined as the process of determining the joint accelerations when the jointactuator torques/forces of the system are known.In order to compute the joint accelerations θ recursively, the GIM, I of Eq. ( 17), is decomposed as I ≡ UDU T (Saha, 1995(Saha, , 1997(Saha, , 1999b) ) based on the Reverse Gaussian Elimination (RGE) method, where U and D are upper triangular and diagonal matrices, respectively.The UDU T decomposition results in an efficient order n, i.e., O(n), computational algorithm in contrast to O(n 3 ) computations required by the Cholesky decomposition of the GIM (Strang, 1998).
For the development of recursive O(n) forward dynamics algorithm, the constrained dynamics equations of motion, Eq. ( 17), are rewritten as where ϕ ≡ τ−C θ.Then, three recursive steps are used to calculate the joint accelerations, which are given below.

Step 1
Solution for τ, where τ ≡ DU T θ ≡ U −1 ϕ.It is found as follows: where η i,i+1 is the 6-dimensional vector obtained recursively as in which η n,n+1 = 0, and the 6-dimensional vector ψ i+1 is evaluated using the following relations: In Eq. ( 30), the 6 × 6 matrix, Mi is obtained recursively as The 6 × 6 symmetric matrix Mi is the mass matrix of Articulated Body, i, defined as the links #i, ..., #n, coupled by the joints i+1, ..., n.This is in contrast to the definition of the Composite Body, i, given after Eq. ( 20), where the links are rigidly connected, i.e., the joints are locked.Note that the mass matrix of the i-th Articulate Body Mi is nothing but the Articulated-Body-Inertia (ABI) of Featherstone (1987).

Step 3
Solution for θ , where, θ ≡ U −T τ.It is found as follows: For i = 2, ..., n, where , and for i = 1, µ 10 ≡ 0. Based on the above mentioned forward dynamics algorithm, another C++ program RFDSIM (Recursive Forward Dynamic and Simulation of Industrial Manipulators) was written which was reported in Saha (1999a).A similar algorithm was rewritten in Visual C# and implemented in the "FDyn" module of RoboAnalyzer software (Rajeevlochana et al., 2012; http://www.roboanalyzer.com) with which one can see animation of the systems under study.The numerical integrator used in RoboAnalyzer for the simulation purposes is based on the Runge-Kutta 4th order method (Bathe and Wilson, 1976).

Numerical example: a 6-DOF Stanford arm
The dynamic analyses of the 6-link 6-DOF serial-chain system with both revolute and prismatic joints, namely, the Stanford arm as shown in Fig. 4, were carried out using RoboAnalyzer.The Denavit and Hartenberg (DH) paramters, which were proposed by Denavit and Hartenberg (1955), and the mass and inertia propoerties are taken from Saha (1999a) as per the notations explained there and in Saha (2008).The numerical values are not reproduced here since the focus of this paper is to review the DeNOC-based formulations and their applicability.However, the joint torques (Joints 1-2, 4-6) and force (Joint 3) obtained from the "IDyn" module of RoboAnalyzer software for the following joint input motions are plotted in Fig. 5: where θ i (0) = 0, for i = 1-2, 4-6 and b 3 (0) = 0 are the variable DH parameters (Saha, 2008) or the joint variables at time T = 0, whereas the total time of motion is, T = 10 s.Gravity was acting in the negative Z 1 -direction.The variable, τ i , for i = 1-2, 4-6, and f 3 in Fig. 5 are the joint torques and force, respectively.The results were verified with those reported in Saha (2008).
The forward dynamics and simulation of the Stanford arm was also performed using "FDyn" module of RoboAnalyzer.The Stanford manipulator was assumed to fall freely under gravity without any external torques and force at the actuating joints.The initial positions were taken same as in the inverse dynamics analysis given after Eq. ( 35).The results are plotted in Fig. 6, where the variations of the joint motions with respect to time are shown.The results were also verified with those reported in Saha (2008).

Tree-type systems
A tree-type system has a set of links connected by kinematic pairs, typically, a revolute or a prismatic joint, as shown in Fig. 7. Other type of joints, say, a universal or spherical, and a cylindrical, can be modelled as a combination of two or three intersecting revolute joints, and a pair of revoluteprismatic joints, respectively, as mentioned in the beginning of Sect.2.1.Based on the modelling of serial-chain systems, Shah et al. (2011Shah et al. ( , 2013) ) extended the methodology to model a tree-type system.For this, the tree-type system was assumed to be a combination of several serial-chain systems called "kinematic modules".Consequently, multi-modular recursive algorithms for the tree-type systems were presented against "full-body-level" recursive dynamics algorithms of Featherstone (1987) and Rodriguez (1992).Each "module" of the tree-type architecture was defined as a set of serially connected links emerges from the last link of its parent module.For example, as indicated in Fig. 8, the parent module of M i is module M β .
For the analyses purposes, the tree-type system was first kinematically modularized before its kinematic constraints were derived.The modules are denoted with M 0 , M 1 , M 2 , etc., where a child module bears a number higher than its parent module.Moreover, the links inside any module, say, M i , are denoted as #1 i , ..., #k i , ..., #η i , where the superscript i signifies the module number.Considering the treetype system, there are s number of modules in the system, and there are η i number of links in the i-th module.The total number of links in the whole system is then obtained by n = and inter-modular (between the modules) recursions, as presented in Fig. 9.

Intra-modular kinematic constraints
Intra-modular kinematic constraints are effectively the velocity constraints between the links of a serial-chain system derived in Eqs.(3-10).Here, however, a little modification is proposed in the definition of each link's linear velocity v i .In contrast to the definition of the velocity of the mass center of the i-th link, C i , as v i of Eq. (1), it is defined in this section as the velocity of point O i where the i-th joint couples the j-th link with the i-th one, as indicated in Fig. 2. Such definition of v i in twist expression of Eq. ( 1) was necessitated mainly to take care of the branching issue of the serial-modules in the tree-type system, as shown in Fig. 8.The velocity of the i-th link defined here with respect to O i (sometimes referred to as the origin of the i-th link).It is actually the velocity of the previous link at its connection point, namely, the last link of the parent module where the first link of the child module is coupled.Hence, where branching occurs no additional computations are required for the calculation of the velocity of the first link belonging to the child module.This was not the case with the definition of the velocity of the i-th link with respect to its mass center C i in which additional computations would be required to calculate the velocity of the mass center C i from the origin O i .Moreover, as the main objective of dynamic analyses is to calculate either joint torques or joint motions, selection of O i as a reference point, instead of the C i , can lead to efficient recursive inverse and forward dynamics algorithms, as shown by Shah et al. (2011Shah et al. ( , 2013)).In fact, for the serial-chain systems considered in Sect.2, the same definition with respect to O i could have been adopted.This was done with the "IDyn" and "FDyn" modules of the RoboAnalyzer software.In Sect.2, however, it was shown how the simplest form of the NE equations of motion given by Eq. ( 1) can be used with the definition of the DeNOC matrices, as demonstrated in the original work of Saha (1995Saha ( , 1997Saha ( , 1999aSaha ( , b, 2003)).Now, with the new definitions of v i with respect to O i , Eq. ( 4) is rewritten as where t i and t j are the 6-dimensional twist vectors defined in Eq. ( 1) but with respect to (w.r.t.) the new definition of v i , i.e., w.r.t.point O i .Accordingly, the 6 × 6 matrix A i j is the new twist-propagation matrix.A different notation is used here to distinguish it from B i j which was defined after Eq. ( 4) w.r.t. the definition of the velocity of C i .The 6 × 6 matrix A i j , and the 6-dimensional joint-rate-propagation vector, p i , are given by , and where the 3-dimensional vector a i j is shown in Fig. 2. Notice the change in the expression of p i in Eq. (36b) in comparison to the same in Eq. ( 5) where v i was defined w.r.t.C i .For serially connected rigid links in the i-th serial-chain module, one can write the expression for the generalized twist, t i , similar to Eq. (9a), as In Eq. ( 37), the 6η i -dimensional generalized twist vector t i and the η i -dimensional generalized joint-rates vector θi are defined as follows: where a bar ("−") over an entity in Eqs. ( 37) and ( 38) signifies that the quantity is related to a module and the superscript, i, outside the brackets identifies the module.As a consequence, the generic notation t k (or t k i ) in Eq. ( 38) is the 6-dimensional twist vector for the k th link in the i-th module.The 6η i ×6η i and 6η ,  (Shah, 2011;Shah et al., 2013).

Inter-modular kinematic constraints
Having obtained the intra-modular kinematic constraints in the velocity-level, it is now possible to derive the intermodular kinematic (velocity) constraints, i.e., between two neighbouring serial-chain modules.In a way, each module has been treated similar to a link in a serial-chain module presented in Sect. 2 or Sect.3.1.For this, module M β is con-sidered as the parent of module M i , as shown in Fig. 8.This is similar to link j of Fig. 2 which is the parent of link i.The 6η i -dimensional generalized twist t i is then obtained from the 6η β -dimensional generalized twist t β as where A i,β is the 6η i × 6η β module-twist-propagation matrix which propagates the generalized twist of the parent module (β) to the child module (i) and N i is the 6η i ×η i module-jointrate propagation matrix, which are given by Mech.Sci., 4, 1-20, 2013 www.mech-sci.net/4/1/2013/and The vectors t i and θi are defined in Eq. ( 38).Next, the 6n-dimensional generalized twist vector t, and the ndimensional generalized joint-rate vector θ, for the whole tree-type system which comprises of s modules and n links are defined as where t 0 and θ0 correspond to the base module M 0 which may not be fixed.For example, in the case of a spacecraft carrying a manipulator, the spacecraft floats with motion of 6-degrees-of-freedom (DOF).For the analysis purposes, its motion need to be specified for further motion analyses of other modules, e.g., the manipulator of the above system.Upon substitution of the expressions of N i from Eq. ( 41b), for i = 1,..., s, in Eq. ( 40), and manipulating the expressions like Eqs. ( 8)-( 9), one obtains the expression of the 6ndimensioanl generalized twist t for the whole tree-type system as in which, N l and N d are the 6(n + 1) × 6(n + 1) and 6(n + 1) × (n + n 0 ) matrices, respectively, as the tree-type system was assumed to have module M 0 with one-link with n 0 DOF.Matrices N l and N d for the tree-type system are given by and In Eq. ( 44), 1 i is the 6η i × 6η i identity matrix, whereas γ i stands for the array of all modules including module M i and outward to it, as shown within dashed line of Fig. 10.The  matrices N l and N d are the desired Decoupled Natural Orthogonal Compliment (DeNOC) matrices for the whole treetype system at hand.Note here that the matrices, N l and N d of Eq. (9b), and N l i and N d i of Eq. ( 39), are the special cases of the DeNOC matrices derived in Eqs. ( 44) and ( 45), where each module has only one link without any branching.

Newton-Euler (NE) equations for tree-type systems
In contrast to the expressions for the Newton-Euler (NE) equations of the i-th link given by Eq. ( 11) or ( 12), a deviation in their expressions will be observed.This is due to the modified definition of the velocity of the i-th link, i.e., v i , with respect to point O i .This was mentioned in Sect.3.1.The NE equations of motion of the k-th link (as the letter "i" will be used to denote module) of the i-th module with respect to point O k can be expressed as (Shah et al., 2011(Shah et al., , 2013) Combining Eqs. ( 46a)-(46b), one can obtain an expression equivalent to Eq. ( 12) as where the 6 × 6 matrices M k , Ω k , and E k are defined as and Note in Eq. (47b), that I k is the 3 × 3 mass moment of inertia tensor of the k-th link about O k .Combining Eq. ( 47a) for all η i links of the i-th module and for all s modules, one can write a compact expression equivalent to Eq. ( 13) as (Shah et al., 2011(Shah et al., , 2013)  where matrices M, Ω, and E are the 6(n + 1) × 6(n + 1) blockdiagonal matrices defined similar to Eq. ( 14).For details, readers are referred to the Ph.D. thesis of Shah (2011) or the book by Shah et al. (2013).

Constrained equations for tree-type systems using the DeNOC matrices
The constrained equations of motion for the tree-type systems are derived in this subsection in a similar manner to that of the serial-chain system of Sect.2, i.e., pre-multiply N T d N T l of Eq. ( 43) to the unconstrained NE equations given by Eq. ( 48) to obtain a set of constrained independent equations of motion by eliminating the constraint wrenches.These constrained equations are also referred to as the Euler-Lagrange equations of motion of the tree-type system at hand.They are given by where I is generalized inertia matrix (GIM), C is the matrix of convective inertia terms (MCI), and τ is the vector of generalized driving forces, and due to gravity, dissipation, external forces, etc., which have expressions similar to those after Eq. ( 17).Note that the expression of Eq. ( 49) is same as Eq. ( 17) but the sizes of the corresponding matrices and vectors are different because they represent two different architectures of the multibody systems.Based on Eq. ( 49), recursive inverse and forward dynamics algorithms for tree-type systems were developed by Shah (2011) and implemented in a software called ReDySim (Recursieve Dynamics Simulator) (Shah et al., 2012c).ReDySim was written in MATLAB environment and available free from http://www.redysim.co.nr.(Shah, 2011;Shah et al., 2013).

Numerical example: a spatial biped
In order to illustrate the recursive dynamics algorithms presented in this section, ReDySim was used to analyze a spatial biped shown in Fig. 11.The model parameters were taken from Shah (2011) which will appear in the book by Shah et al. (2013) also.They are not reproduced here due to the reasons cited in Sect.2.8.However, the designed input motions of the trunk's centre-of-mass (COM) and ankle for stable walking (Shah, 2011) are shown in Figs. 12 and 13, respectively.Based on the inputs of Figs. 12 and 13, the inverse dynamics results were obtained which are shown in Fig. 14.Forced simulation was performed next, as reported in Shah (2011), where the motion of the biped was studied under the application of joint torques calculated above, i.e.,

Mech
(b) Swing leg (Module 2) Figure 13.Joint trajectories of the biped obtained from the trajectories of trunk and ankle (Shah, 2011;Shah et al., 2013) Figure 13.Joint trajectories of the biped obtained from the trajectories of trunk and ankle (Shah, 2011;Shah et al., 2013).
those shown in Fig. 14.The joint motions were calculated using the forward dynamics module of ReDySim.The plots for the simulated joint angles are shown in Fig. 15, along with the desired one.It can be seen that the simulated joint angles match with the desired joint angles up to 0.1 s, i.e., until 0.1 s movement of the biped.After this, the system behaves unexpectedly as evident from the divergent plots of the simulated angles in Fig. 15a.The deviation in the simulated angles is mainly attributed to what is known as zero eigen-value effect (Saha and Schiehlen, 2001).The physical system may also not behave as expected due to disturbances caused by unmodelled parameters like friction, backlash, etc., and nonexact geometrical and inertia parameters.Hence, a control scheme must be considered, as this forms a part and parcel of achieving proper walking.These aspects were explained in detail in Shah (2011) and Shah et al. (2013), and not elaborated further due to space limitation of the paper.Note several advantages of the concept of the kinematic modules in the dynamics modelling of tree-type systems consisting of serially connected links (Shah, 2011;Shah et al., 2013), which are as follows: -Extension of the body-to-body velocity transformation relationship to module-to-module velocity transformation relationship.
-Compact representation of the system's kinematic and dynamic models.
-Uniform development of the inverse and forward dynamics algorithms with inter-and intra-modular recursions.
-Module-level analytical expressions of the matrices and vectors appearing in the equations of motion.
-Ease of investigation of any inconsistency in the results of modules without the need to investigate the whole system.
-Possibility of hybrid recursive-parallel algorithms, where each module can be analyzed using recursive relations in parallel.

Closed-loop systems
The DeNOC-based dynamic modelling of serial-chain and tree-type open-loop systems presented in Sects. 2 and 3, respectively, can be extended to closed-loop systems provided one cuts the closed-loops of a system at suitable locations to make it open.Note that, one needs to use suitable constraint forces at the cut-joints to represent the actual presence of the joints.Such constraint forces are known in the literature as Lagrange multipliers (Chaudhary andSaha, 2007, 2009, andothers).The multipliers need to be evaluated from the loop-closure constraints before they can be used as external forces to the resulting open-loop systems.In this section, a planar 4-bar mechanism shown in Fig. 16 is considered to illustrate the concept.However, the methodology is applicable to any general multi closed-loop systems, as shown in Chaudhary and Saha (2007), Shah (2011), andShah et al. (2013).
Note that to model an open-loop system resulting from a closed-loop system, one needs to re-write Eqs. ( 13) or ( 48) as where w λ is the 6n-dimensional vector of generalized wrench due to Lagrange multipliers acting at the cut joints.For the 4-bar mechanism shown in Fig. 16a, the two cut-open serialchain subsystems are shown in Fig. 16b.The resulting openloop tree-type subsystems have one and two links, respectively, connected by one and two one-DOF revolute joints.
Other terms have same meaning as in Sects. 2 and 3.In Eq. ( 50), N T w C = 0 for the reason given after Eq. ( 16), but N T w λ 0. These terms are now the new unknowns to the inverse and forward dynamics problems that need to be evaluated with the help of loop-closure constraints.
For the closed-loop 1-2-3-4 of the 4-bar mechanism shown in Fig. 16a, one can write where 2-dimensional vectors of the planar system, a i for i = 0, 1, 2, 3, represent the relative position vectors of the joints in the 4-bar mechanism, Fig. 16a.Differentiating Eq. ( 51) with respect to time, one obtains J θ = 0, where θ ≡ θ1 θ2 θ3  Figure 15.Simulated joint angles for the biped (Shah, 2011;Shah et al., 2013) Figure 15.Simulated joint angles for the biped (Shah, 2011;Shah et al., 2013).and the 2×3 Jacobian matrix for the 4-bar mechanism at hand can be given by where s 12 ≡ sin(θ 1 + θ 2 ) , s 123 ≡ sin(θ 1 + θ 2 + θ 3 ), and similarly c 12 and c 123 etc. Equation (52b) was also derived in Chaudhary et al. (2007Chaudhary et al. ( , 2009) ) as , where J I ≡ A I en N I and J II ≡ A II en N II (53) In Eq. ( 53), N I and N II are the 6 × 1 and 12 × 2 NOC matrices for the two open-chain subsystems I and II, respectively, shown in Fig. 16b, whereas A I en and A II en are the 2 × 6 and 2 × 12 twist-propagation matrices for the last link from the point of contact to its previous link to the point where the joint is cut, i.e., Joint 4 of Fig. 16b.Such Jacobians using the DeNOC matrices were also derived in Saha (2008).The resulting constrained dynamic equations of motion for the two subsystems are then written as Depending on the type of dynamics problem, i.e., inverse or forward, Eqs.(54a)-(54b) can be solved using recursive "subsystem" or "system" approach.For inverse dynamics, "subsystem recursion" provides a better efficiency, as pointed out by Chaudhary and Saha (2009).   of the crank, output link and coupler were taken as 1.5 kgs, 3 kgs and 5 kgs, respectively.The input joint angle and the joint torque at joint 1 are plotted in Fig. 17.The forward dynamics of the 4-bar mechanism was carried out using the same initial configuration, as specified in the inverse dynamics.The simulation was done for the free-fall of the mechanism under gravity without any external torque applied at joint 1.The joint angles and rates are plotted in Fig. 18.The results of inverse dynamics and forward dynamics were validated with MATLAB's SimMechanics model, as reported in Shah et al. (2012c).

Computational efficiency
In this section, computational efficiencies of the DeNOCbased algorithms are investigated.Figure 18 shows comparisons of computational complexities required by several inverse dynamics algorithms, when a system has 1-DOF (Fig. 19a), 2-DOF (Fig. 19b), 3-DOF (Fig. 19c) and equal numbers of 1-2-and 3-DOF joints (Fig. 19d).It may be seen that the recursive inverse dynamics algorithm given in Fig. 9 of Sect. 3 performs as fast as the fastest algorithm available in the literature when the system has only 1-DOF joints, as evident from Fig. 19a.However, when multiple-DOF joints are introduced in the system, the algorithms of Sect. 3 (Shah, 2011), which have been implemented in ReDySim, outperforms the other algorithms available in the literature.This is clear from Fig. 19b-d.
From Fig. 20, it is also clear that the forward dynamics algorithm of Fig. 9 explained in Sect. 3 (Shah, 2011) performs better than any other algorithm available in the literature.More the number of multiple-DOF joints, more the improvement in the computational efficiency, as shown in Fig. 20b-d.This is mainly due the implicit inversion of the GIM based on the Reverse Gaussian Elimination (Saha, 1995(Saha, , 1997) ) of the GIM, and simplification of the expressions associated with multiple-DOF joints (Shah et al., 2012b).In the tree-type robotic systems, such as biped, quadruped, etc., where the DOF of the system is more than 30, and the system consists of many multiple-DOF joints, the DeNOC-based algorithms significantly improve the computational efficiency.

Software for students and researchers
In order to build up the interest in the areas of multibody dynamics and to provide efficient tools to perform the dynamic analyses, the following multibody simulation tools  (Shah, 2011;Shah et al., 2013).
were developed for the students and researchers using the DeNOC-based formulations presented in Sects.2-4: -RoboAnalyzer: for serial-chain open-loop systems -ReDySim (Recursive Dynamic Simulator): for general tree-type systems These are explained briefly in the following subsections.

RoboAnalyzer
RoboAnalyzer (Rajeevlochana et al., 2012) is a 3dimensional model-based software to solve kinematics and dynamics problems of an serial-chain open-loop system.It was developed using Visual C# and OpenGL that take the description of a serial-chain system using the DH parameters (Denavit and Hartenberg, 1955;Saha, 2008), and the mass and inertia properties of each link.In RoboAnalyzer, one can also see the animation of the analyzed systems.For the benefit of the users, CAD models of the standard systems like KUKA, PUMA robots, and others were made available for analysis.The software is freely downloadable from http://www.roboanalyzer.com.

ReDySim (Recursive Dynamics Simulator)
Recursive Dynamic Simulator (ReDySim) (Shah et al., 2012c) is a multibody dynamics simulation tool which was developed in MATLAB environment.It was developed based on the concept of the DeNOC matrices, and kinematic modules of a tree-type system, as explained in Sect.3. ReDySim can be used to perform inverse dynamics and simulation of multibody systems.It has two modules, namely, fixed-base and floating-base modules.The latter was not presented in this paper due to limited page restriction of a paper.However, the interested readers can refer to Shah (2011) or Shah et al. (2013) for the dynamics analyses of biped, quadruped and six-legged walking robots using the floating-base concept.Note here that the architectural information of the treetype systems were provided using the modified-DH parameters (MDH) parameters, as proposed by Khalil and Kleinfinger (1986), instead of the DH parameters defined in Saha (2008).This was done mainly to improve the computational efficiency of the tree-type systems.
ReDySim was also used to solve flexible systems like ropes, etc. which were modelled as hyper-degrees-offreedom rigid-link systems (Shah, 2011).In fact, the simulation of long chains with the aid of ReDySim showed considerable improvement over commercial software like Recur-Dyn in terms of the computational time and correctness of www.mech-sci.net/4/1/2013/Mech.Sci., 4, 1-20, 2013  (Shah, 2011;Shah et al., 2013).
the results.These results were separately communicated to a journal for publication (Agarwal et al., 2012).ReDySim can be downloaded free from http://www.redysim.co.nr,where the user's manual and some demos are also made available for the benefits of the users.

Conclusions
Recursive algorithms are popular due to their efficiency, computational uniformity, and numerical stability.This paper gave insight to the evolution of the Decoupled Natural Orthogonal Complement (DeNOC) matrices for more than one and half decades.It was shown that the use of DeNOC matrices in the dynamic analyses of multibody systems led to the development of several recursive dynamics algorithms for serial-chain, tree-type, and closed-loop systems.Several numerical examples, e.g., the Stanford arm (serial-chain), spatial biped (tree-type), and 4-bar mechanism (closed-loop), were given to explain the concepts presented in this paper.
Advantages and computational efficiencies of the DeNOCbased methodologies suggest that it should be used when the DOF of a system is large and the application is real-time.Finally, two types of software meant for serial-chain open-loop systems, and general tree-type systems with abilities to solve closed-loop systems, i.e., RoboAnalyzer and ReDySim, re-spectively, were explained.It is expected that the algorithms, and more importantly, the software will benefit immensely the students and researchers of multibody dynamics community.

Figure
Figure 1.A robot manipulator

1Figure
Figure 3. Free-body diagram of the ith link 2 3

Figure 3 .
Figure 3. Free-body diagram of the i-th link.

Figure 4 .
Figure 3. Free-body diagram of the ith link 2 3

Figure 8 .
Figure 7.A tree-type system

Figure 8 .
Figure 7.A tree-type system
Figure 10.Definition of  i