THEORETICAL & APPLIED MECHANICS LETTERS 2, 063011 (2012)

Recursive dynamics simulator (ReDySim): A multibody dynamics solver

Suril V. Shah,1, a) Paramanand V. Nandihal,2, b) and Subir K. Saha2, c)
1)
Department of Mechanical Engineering McGill University, QC H3A0C3, Canada
2)
Department of Mechanical Engineering Indian Institute of Technology, New Delhi 110016, India
(Received 10 September 2012; accepted 25 September 2012; published online 10 November 2012)
Abstract Recursive formulations have significantly helped in achieving real-time computations and
model-based control laws. The recursive dynamics simulator (ReDySim) is a MATLAB-based recur-
sive solver for dynamic analysis of multibody systems. ReDySim delves upon the decoupled natural
orthogonal complement approach originally developed for serial-chain manipulators. In comparison
to the commercially available software, dynamic analyses in ReDySim can be performed without
creating solid model. The input parameters are specified in MATLAB environment. ReDySim has
the capability to incorporate any control algorithm with utmost ease. In this work, the capabilities
of ReDySim for solving open-loop and closed-loop systems are shown by examples of robotic gripper,
KUKA KR5 industrial manipulator and four-bar mechanism. ReDySim can be downloaded for free
from http://www.redysim.co.nr and can be used almost instantly. ⃝ c 2012 The Chinese Society of
Theoretical and Applied Mechanics. [doi:10.1063/2.1206311]
Keywords ReDySim, multibody systems, dynamic modeling, recursive dynamics, DeNOC

Multibody dynamics find applications in robotics, have significantly helped in achieving real-time com-
automobile, aerospace and many other streams for anal- putations and model-based control laws. ReDySim6–8
ysis, simulation and control. It has evolved a lot in the is a MATLAB-based solver, which is a general pur-
last two decades and there is a huge scope of research pose platform, essentially consisting of very efficient re-
for scientists and engineers. Computer-aided dynamic cursive order (n) inverse and forward dynamics algo-
analysis of multibody systems has been the prime mo- rithms for simulation and control of a tree-type system.
tive of the engineers since the evolution of high speed ReDySim delves upon the decoupled natural orthogo-
facilities using computers. Dynamic analysis involves nal complement (DeNOC) approach9 originally devel-
force and motion analyses. Force or inverse dynamics oped for serial-chain manipulators. The specialty of
analysis attempts to find the driving and reactive forces this solver exists in its recursive nature and flexibility in
for the given input motion, whereas the motion analysis solving complex problems. The gain in computational
or forward dynamics, obtains system’s configuration un- time is more as the number of links and joints in the
der the input forces. Force analysis helps in design and system increases.
control of multibody systems, whereas motion analysis This paper mainly addresses analysis of fixed-base
allows one to study and test a design virtually without open- and closed-loop systems using ReDySim. The
really building a real prototype. An efficient framework floating-base module of RedySim is also available on
is essential for the dynamic analysis of complex multi- http://www.redysim.co.nr, and can be used for analyz-
body systems. ing legged robots7 and space robots with mobile or free
There are several software and toolboxes avail- base. Details of floating-base module are not given due
able for simulating multibody systems. They are in to space limitation.
the form of toolkit for MATLAB1,2 and LabVIEW.3 Even though, the algorithms in ReDySim are meant
The commercially available tools such as ADAMS4 and for tree-type systems as shown in Fig. 1, it can effec-
RecurDyn5 are the general purpose software used for tively be used to solve closed-loop systems by simply
dynamic analysis. These software are well suited for providing constrained Jacobian matrix resulting out of
some standard industrial applications, however, they loop-closure equations and letting ReDySim know how a
fail to attract research community mainly due to their closed loop system is cut-open. The flow chart showing
inflexibility at user’s end for solving complex problems. inverse and forward dynamics algorithms6 of the tree-
The generic nature of this software many times com- type systems are given in Fig. 2. The physical configu-
promise with accuracy of the results even for smaller ration of the system is mainly defined by the Denavit-
system for longer simulation time. Hartenberg (DH)10 parameters, as proposed by Khalil
In this context, recursive dynamics algorithms play and Kleinfinger.11 The user inputs for both the inverse
an important role. They are attractive due to simplicity and forward dynamics are discussed.
and computational uniformity regardless of ever grow- In order to perform inverse dynamics, the following
ing complex multibody systems. Recursive formulations input parameters are required:
Model parameters:
(1) number of links (nl );
a) Correspondingauthor. Email: [email protected]. (2) type of system (type), i.e., open-loop or closed-
b) Email: [email protected]. loop;
c) Email: [email protected].
(3) degrees-of-freedom (dof) of the system;
063011-2 S. V. Shah, P. V. Nandihal, and S. K. Saha Theor. Appl. Mech. Lett. 2, 063011 (2012)

to these input parameters, the following parameters are

required for the purpose of integration:
Alink or (1) initial conditions y0 = [q T q̇ T Eact ]T where q,
body q̇ and Eact are initial joint positions, joint rates, and
actuator energy, respectively;
(2) initial time (ti ) and final time (tf ) of simulation,
and step size (dt);
(3) relative tolerance (rtol ) and absolute tolerance
(atol ), and the type of integrator, note that one may use
either adaptive solver ode45 (for non-stiff problem) or
ode15s (for stiff problem) or fixed step solver ode5 by
specifying the index 0, 1 or 2, respectively.
Base 0
These parameters are entered in the function file
named initials.m. The joint torque on each joint, which
Fig. 1. Tree-type system. are input for forced simulation, can be entered in the
function file torque.m. The default value of torque at
each joint is kept zero which lead to the free simulation.
(4) actuated joint in the system (aj ); However, in the case of force simulation, the user can de-
(5) vector containing number of joint variables as- fine proportional (P), proportional and derivative (PD)
sociated with each joint (nj ); or model-based control law for torque input as the cur-
(6) constant Denavit-Hartenberg (DH) parameters rent time (t), number of links (nl ), and vectors of joint
for revolute joints (a, α and b); positions (θ) and joint velocities (θ̇) are passed to the
(7) parent of each link (β); function torque.m. One can also integrate any user de-
(8) vector dk measured from origin Ok to the center- fined control algorithm in this function file. Dynamic
of-mass (COM) Ck of the kth link; simulation of closed-loop system also requires entering
(9) mass of each link (mk ); Jacobain and its time derivative in the function file ja-
(10) vector of gravitational acceleration (g) in the cobian.m. It is worth noting that the vectors of current
inertial frame; joint angels (θ) and joint velocities (θ̇) are passed as the
(11) inertia tensor of each link about COM and rep- inputs to this function and outputs are Jacobian matrix
resented in body-fixed frame (IkC ); and its time derivative.
(12) time period (Tp ) and step size (dt). Finally, simulation is performed by running pro-
The above input parameters are entered in the func- tected function file runfor.p. The output joint motions
tion file named inputs.m. The joint torques are then are stored in data file statevar.dat whereas the time his-
obtained by calling protected function file runinv.p. tory is stored in timevar.dat. The joint motions can be
User may call function file plot tor.m to see the out- plotted by using function file plot motion.m. The func-
put torques. The function file runinv.p also calculates tion file for kine.m can be used to calculate the position
Lagrange multipliers at the cut joints for closed-loop of the link origin and COM, velocity, angular velocity
systems. and the tip position. Moreover, the total energy can
In order to generate inverse dynamics results, cy- be calculated by running the file energy.p and plotted
cloidal joint trajectory is as default. However a user using plot en.m. This can be used for the purpose of
can define any trajectory through the function file tra- validating the simulation results. The system can also
jectory.m. Using the defined trajectory, the calcu- be animated using the file animate.m.
lated joint torque values are stored in data file tor.dat. Dynamic analyses of serial, tree-type and closed-
The numerical results of joint trajectories are stored in loop systems are presented next using ReDySim. The
traj.dat. One can also visualize input joint trajectory by planar 4-dof robotic gripper and a spatial 6-dof indus-
using file animate.m. The file animate.m is not generic trial manipulator KUKA KR512 are selected as open-
and need to be modified for animating different systems. loop systems, whereas a 1-dof planar four-bar mecha-
In order to solve the inverse dynamics of a closed- nism is selected as an example of closed-loop system.
loop system, first the trajectories of independent joints A tree-type robotic gripper, as shown in Fig. 3, can
are entered in the function trajectory.m whereas the hold objects to be manipulated by a robotic manipu-
relationship between independent and dependent joint lator. Numerical results for inverse dynamics, i.e., to
trajectories and Jacobian matrix are entered in the func- find the joint torques for a given set of input motions,
tion file inv kine.m. The inv kine.m file is specific to a were obtained using the inverse dynamics module of the
given system and the user is required to modify it de- ReDySim. The detailed steps are shown in Fig 2(a)
pending on the type of system to be analyzed. where the definitions of all variables are available in
In order to perform forward dynamics, input pa- Ref. 6. The motion for each joint for this purpose was
rameters are provided in the file inputs.m. These input computed using
parameters are nothing but the first eleven entities of θ (T ) − θ (0)
model parameters, as discussed in prevous. In addition θ (t) = θ (0) + ·
T
063011-3 Recursive dynamics simulator Theor. Appl. Mech. Lett. 2, 063011 (2012)

Joint torques (-
τi), inertia parameters
-
(Mi), twist and motion propagations
- -
(Ai,j and Ni), and initial conditions
- .
qi and -
qi

-- . -.. i/.s
Joint motions (q i, qi,qi), inertia .
-/ - - --

recursion
-

Forward
ti Ai,βtβi +Niθi
parameters (Mi), twist and . . . .-.
- -- -- -
motion propagations ti/Ai,βtβi +Ai,βtβi+Niθi
- and N - --. - - --
(q -
~ i*/M
i,β i) w iti +ΩiMiEiti

i/.s i/.s
-/ - - -- .
-
^ i/M -
^ iN-i - - - - - ^-
, I^i/NiT Ψ ^, Ψ ^
recursion

ti Ai,βtβi +Niθi Ψ ΨiIi -1

i/Ψ
Forward

-. .
-- -- . .-
- . -- .. i

ti/Ai,βtβi +Ai,βtβi+Niθi+Niθi -
^
ϕ i/ϕ - -
ϕi -Ni ηi T-~

recursion
Backward
--. - - --
-
w~ i/M iti +ΩiMiEiti ~ i/-
-
ϕ ^ -1-
I i
^i
ϕ
-
^ -
^ -
^- T - -- ^i +-~i
M i,i/Mi -Ψ i Ψi , η i/ΨΨi ϕ η
-
^ -^ -T -
^ -
Mβ(i)/Mβ(i)+Ai,β(i)M i,iAi,β(i)
- - -
~β(i)/η
η η~β(i)+ATi,β(i)-
ηi
i/s.
recursion
Backward

- -- ~i
τi/Niw
i/.s

recursion
- - - -
~ β(i)+ATi,β(i)w
~ β(i)/w ~i
Forward
w - - - .. -
~i/Ai,β(i)(N
µ ~
β(i)qβ(i)+µβ(i))
.. - -
- ϕi-Ψi T-
~
qi/ϕ ~i
µ

..
Joint torques (-
τi) -
Independent accelerations q

(a) Inverse dynamics (b) Forward dynamics

Fig. 2. Recursive dynamics algorithms.6

X

θ O
O .
ϑ  .
ϑ
ϑ
θ X ϑ
X
X  O
O
ϑ ϑ . ϑ
θ
ϑ
Y  O O
Y .

.
O,O  θ X O,O
ϑ

Fig. 3. Robotic gripper.

[ ( )]
T 2π ple joint torques for the robotic gripper are then plot-
t− sin t , (1)
2π T ted in Fig. 4. In order to validate the results, a CAD
where θ(0) and θ(T ) denote the initial and final joint an- model of the gripper was developed in ADAMS4 soft-
gles, respectively, as given in Table 1. The joint trajec- ware and used for the computation of the joint torques.
tory in Eq. (1) is so chosen that the initial and final joint The results were superimposed in Fig. 4, which show
rates and accelerations for all the joints are zero. These close match with the values obtained using the proposed
ensure smooth joint motions. The lengths and masses of inverse dynamics algorithm of the ReDySim. ReDySim
all links were taken as l1 = 0.1 m, l2 = l3 = l4 = 0.05 m, took only 0.025 s on Intel [email protected] GHz computing
m1 = 0.4 kg, and m2 = m3 = m4 = 0.2 kg. The sam- system. The ADAMS software, however, took 1.95 s,
063011-4 S. V. Shah, P. V. Nandihal, and S. K. Saha Theor. Appl. Mech. Lett. 2, 063011 (2012)

1.0 0.3
ReDySim
0.8 ADAMS 0.2

0.1

τ2/Nm
τ1/Nm 0.6

0.4 0

0.2 -0.1

0 -0.2
0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0
Time/s Time/s

Fig. 4. Joint torques for robotic gripper.

-60 20
ReDySim
ADAMS
10
-80

θ2 /(Ο)
θ1 /(Ο)

0
-100
-10

-120 -20
0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0
Time/s Time/s

Fig. 5. Simulated Joint angles for robotic gripper (Joints 1 and 2).

Fig. 6. Joint torques for KUKA KR5.

Fig. 7. Simulated joint angles of the KUKA KR5 (Joints 1, 2 and 3).
063011-5 Recursive dynamics simulator Theor. Appl. Mech. Lett. 2, 063011 (2012)

150 Total Potential Actuator Kinetic

3
100 θ λ2

50
Energy/J
λ1
0 ϑ
ϑ
-50
θ
-100 2 θ 1
-150
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ϑ
Time/s
Fig. 9. A four-bar mechanism.
Fig. 8. Energy plot of KUKA KR5.
0.5
Table 1. Initial and final joint angles for robotic gripper. 0.4
Joints 1 2 3 4 0.3
θ(0) 0◦ 0◦ 0◦ 90◦

Driving torque/Nm
θ(T ) 60◦ 80◦ 80◦ 120◦ 0.2

0.1

which is longer and is expected as it is general purpose 0

software. -0.1
The forward dynamics simulation of the robotic
gripper was then carried out using forward dynamics -0.2
module of the ReDySim. Numerical results for the ac-
-0.3
celeration were obtained for the free-fall of the grip-
per, i.e., it was left to move under gravity without any -0.4
external torques applied at the joints. The accelera- 0 0.5 1.0 1.5
Time/s
tions were then numerically integrated twice using the
ordinary differential equation (ODE) solver ode45 (de-
fault in ReDySim) of MATLAB. The ode45 solver is Fig. 10. Inverse dynamics results for four-bar mechanism.
based on the explicit Runge-Kutta formula given in Dor-
mand and Prince.13 The initial joint angles and rates
were taken as θ1 = −60◦ , θ2 = θ3 = θ4 = 0◦ , and Next, forward dynamics of the KUKA KR5 robot
θ̇1 = θ̇2 = θ̇3 = θ̇4 = 0 rad/s respectively. Figure 5 was performed for the free-fall under gravity. For this,
shows the comparison of the simulated joint angles and both initial joint angles and velocities are assumed to
the same obtained in ADAMS (by using RKF45 solver) be zero. Figure 7 shows the plot of joint angles for the
over the time duration of 1 s with the step size of 0.001 s. time period of 1 s. The law of conservation of energy
The results obtained using ReDySim is in close match was used to validate the simulated results. Since there
with those obtained using ADAMS. The ReDySim took is no dissipation in the system, the total energy is plot-
0.29 s in contrast to 39 s required by ADAMS. ted in Fig. 8, which remained unchanged throughout
Dynamic analysis of a spatial 6-dof industrial ma- the simulation time. Hence, the simulation results are
nipulator KUKA KR5 was performed next. To study validated.
the dynamic behavior of the KUKA KR5, its kinematic Dynamic analysis of a closed-loop planar four-bar
architecture12 was used with user-defined mass and in- mechanism was carried out next. First, inverse dynam-
ertia properties. The DH parameters considered for ics was carried out using ReDySim’s inverse dynamics
KUKA KR5 and its assumed mass and inertia prop- module. The link length of crank, output link, cou-
erties are given in the Table 2. pler and fixed base were taken as 0.038 m, 0.115 2 m,
First, inverse dynamics of the KUKA KR5 robot 0.115 2 m and 0.089 5 m respectively. Also, the mass of
was performed for initial and final joint angles of all the crank, output link and coupler were taken as 1.5 kg,
the six joints as 0◦ and 60◦ , respectively. Once again 3 kg and 5 kg respectively. Inertia tensor of each link
cycloidal trajectory of Eq. (1) was used as the input about its COM was calculated by assuming each link as
joint motions. The joint torques at first three joints are a slender rod. In order to solve the four-bar mechanism,
shown in Fig. 6. Note that the torque required at joint 1 it was cut at an appropriate joint, as shown in Fig. 9, to
is zero at the beginning and at the end, as joint 1 is not form a tree-type system.14 The opened joint was sub-
affected by the gravitational acceleration. Moreover, stituted with suitable constraint forces (λ) known as
the torque requirement at joint 2 is maximum, which Lagrange multipliers. For this purpose, the Jocobian
is very obvious as this joint is required to overcome the matrix was provided in inv kine.m function file. The
effect of the gravity on links 2, 3, 4, 5 and 6. constant angular velocity of 4.712 4 rad/s was provided
063011-6 S. V. Shah, P. V. Nandihal, and S. K. Saha Theor. Appl. Mech. Lett. 2, 063011 (2012)

Table 2. DH parameters and inertia properties of KUKA KR5.

◦
i ai /m αi /( ) bi /m θi /(◦ ) mi /kg Ii,xx /(kg · m−2 ) Ii,yy /(kg · m−2 ) Ii,zz /(kg · m−2 )
1 0 0 0.400 θ1 16.038 0.176 9 0.266 5 0.262 2
2 0.180 90 0 θ2 7.988 0.271 6 0.276 8 0.021 8
3 0.600 0 0 θ3 12.937 0.388 9 0.376 5 0.104 1
4 0.120 90 0.620 θ4 2.051 0.004 7 0.010 1 0.012 1
5 0 –90 0 θ5 0.811 0.000 7 0.001 7 0.001 8
6 0 90 0 θ6 0.008 0.000 003 0.000 001 0.000 001

250 (ReDySim) has been developed in MATLAB for

θ analyses of multibody systems. Capability of the
200
ReDySim is shown for serial, tree-type and closed-loop
150
systems whose results are validated with the other
commercial software. The ability of ReDySim in
100 including desired trajectory and control law, provide
θ flexibility to researcher in incorporating their cus-
50 tomized algorithms. The control aspect is not reported
θi/(Ο)

here due to space limitation and will be communicated

0 in future. ReDySim does not require building model in
the software environment before simulation; users can
-50
simply provide the input parameters in the MATLAB
-100
environment. The ReDySim showed considerable im-
provement over commercial software such as ADAMS
-150 and algorithms available in literature in terms of the
θ computational time. The ReDySim can be downloaded
-200 free from http://www.redysim.co.nr. The users are
0 0.2 0.4 0.6 0.8 encouraged to send their comments and suggestions to
Time/s
[email protected] or the authors.
Fig. 11. Simulated joint motions for the four-bar mecha- 1. P. I. Corke, Robot. Automat. Mag. IEEE 3, 24 (1996).
nism. 2. M. Toz, and S. Kucuk, Comput. Appl. Eng. Edu. 18, 319
(2010).
3. T. J. Mateo Sanguino, and J. M. A. Márquez, Simulation
Tool for Teaching and Learning 3D Kinematics Workspaces of
as the input to the joint 1 of the four-bar mechanism. Serial Robotic Arms with up to 5-dof. Computer Applications
The constant angular velocity value was provided in in Engineering Education, 2010.
4. Automated Dynamic Analysis of Mechanical System
the function file trajectory.m. The inverse dynamics (ADAMS), Version 2005.0.0, MSC Software (2004).
was carried out for the time period of 1.33 s. The joint 5. RecurDyn. FunctionBay Inc. 2009.
torque at joint 1 was calculated by executing function 6. S. V. Shah, Modular framework for dynamic modeling and
file runinv.p. The joint torque at joint 1 is plotted in analysis of tree-type robotic systems [PhD Thesis]. (Indian
Institute of Technology Delhi, Delhi, 2011).
Fig. 10.
7. S. V. Shah, S. K. Saha, and J. K. Dutt, ASME J. Nonlinear
Next, simulation of the four-bar mechanism was Computat. Dyn. 7, (2012).
carried out using ReDySim’s forward dynamics mod- 8. S. V.Shah, S. K. Saha, and J. K. Dutt, Mech. Mach. Theory
ule. The simulation was done for free-fall under gravity 49, 234 (2012).
for time period of 1.33 s without any external torque 9. S. K. Saha, ASME J. Appl. Mech. 66, 986 (1999).
10. J. Denavit, and R. S. Hartenberg, ASME J. Appl. Mech. 215
applied at the joint 1. The Jacobian matrix was pro- (1955).
vided in jocobian.m function file. By executing function 11. W. Khalil, and J. Kleinfinger, IEEE Int. Conf. on Robotics
file runfor.p, the joint angles were calculated, which are and Automation, 1174 1986.
shown in Fig. 11. 12. KUKA KR5, Technical specification, KUKA Robotics web-
site, http://www.kuka.com, last accessed on 20 May 2012.
The results of the ReDySim were validated with
13. J. R. Dormand, and P. J. Prince, J. Comput. Appl. Math. 6,
MATLAB’s SimMechanics as shown by circular markers 19 (1980).
in Figs. 10 and 11. 14. H. Chaudhary, and S. K. Saha, ASME J. Mech. Design 129,
An efficient Recursive Dynamics Simulator 1234 (2007).