Simple method to the dynamic modeling of industrial robot subject to constraint

Abstract

It is a common movement type in industrial production for the end-effector of industrial robot executing tasks according to desired track. High nonlinearity and coupling are shown from the dynamic characteristic of industrial robot because of the constrained relationship brought by desired task. Therefore, it is difficult to establish dynamical equation with traditional Lagrange equation. Aiming at the dynamics modeling of industrial robot subject to constraint, the additional torque and dynamical equation of industrial robot subject to some desired trajectory is acquired based on the famous Udwadia–Kalaba equation in analytical mechanics field. The simple approach overcomes the disadvantage of obtaining dynamical equation from traditional Lagrange equation by Lagrange multiplier. But the numerical results obtained by integrating the constrained dynamic equation yield the errors. A numerical algorithm is proposed to reduce the errors by modifying the derived dynamic equation. The simulation results prove that the modified algorithm can be more precisely utilized in dynamics modeling of industrial robot subject to constraint.

Keywords

Constraint industrial robot dynamics modeling errors reducing Udwadia–Kalaba equation

Introduction

For industrial robot, the accurate kinematic and dynamical models are the basis for completing all kinds of tasks. Compared with mature kinematic model, the dynamical model, especially the industrial robot whose end-effector moves according to desired track, shows high nonlinearity and coupling in its dynamic characteristic for the constrained relationship brought by desired track. Therefore, it is an important basis for realizing high-quality movement control to establish mathematical model which can describe dynamic characteristics of movement limitation system. At present, most of the scholars pay their attention to control algorithm of industrial robot, while few of them concern about the dynamical modeling problem. A constrained robotic system is a mechanical system, where we need to consider kinematic constraint conditions explicitly in dynamic modeling process. MH Korayem et al.¹ used the traditional Lagrange multiplier method to deal with dynamical modeling of mobile manipulator. The disadvantage of the Lagrange multiplier is relying on special approaches to the determination of the multipliers, and it is often very difficult to find the multipliers to obtain the explicit equations for systems that have large numbers of degrees of freedom and a series of non-integrable constraints. P Sanchez-Sanchez and MA Arteaga-Perez² proposed a simple method for obtaining the dynamic model of robot manipulators. But it is only beneficial to establish the dynamic equation for the unconstrained system. Kane³ proposed a general method to obtain the dynamical system’s equation of motion which originated from Gibbs and Appell’s quasi-coordinates. Y-H Jia et al.⁴ derived the equations of motion for the unconstrained planar dual-arm space robot using Kane’s equations. By inserting the constraint equations in terms of accelerations into the equations of motion for the unconstrained system, the reduced-order equations for the constrained system were obtained which contain only independent motion variables. Kane’s equation applies to holonomic constraints and also non-holonomic constraints and does not involve Lagrange multiplier. Computational complexity can be largely reduced for systems with a large number of degrees of freedom compared with Newtonian and Lagrange methods. Certainly, Kane’s method has its own disadvantage. It requires a suitable choice of problem-specific generalized speeds which relies on experience in a large part. SK Saha et al.⁵ reviewed different methods on constrained multi-body dynamics formulations and summarized the advantages and disadvantages, respectively.

It is well known that it is impossible to create a new fundamental principle for the theory of motion and equilibrium of discrete, dynamical systems. However, as early as in 1996, Firdaus Udwadia of University of Southern California comes up with a multi-body dynamical modeling method that is Udwadia–Kalaba equation and obtains many research results. The Udwadia–Kalaba equation has been proposed that is useful to help us understand nature’s law from new points of view. The Udwadia–Kalaba equation describes the dynamics of constrained systems from Gauss’ principle. This method can establish system dynamical equations under holonomic constraint and non-holonomic constraint in a relatively simple way. Using this method, one can get analytical expression of constraint force without Lagrange multiplier.^6–15 Unfortunately, in the hereafter 20 years until now, this method has not been paid too much attention, and its main application field focuses on dynamical modeling and control research of satellite system. It has seldom been used in robot field. Only J Huang et al.¹⁶ used this method to build a parallel robot dynamical model, and H Zhao et al.¹⁷ established the dynamical model of a machine fish.

The modeling process of the Udwadia–Kalaba equation has three steps. First, in terms of the generalized coordinates, the unconstrained dynamical system whose equations of motion can be written using traditional Lagrange method is considered. Then, the constraint equations are formed. In the end, the additional generalized forces of constraint are imposed on the system. That is to say, add the force resulted from the presence of constraints to the unconstrained system’s force.

In this article, the motion path of the end-effector is regarded as a constraint, and the dynamic equation is obtained based on the Udwadia–Kalaba equation. However, the results using numerical integration of the dynamic equations do not satisfy the constraint equations precisely. A modified approach is proposed to reduce the errors and obtain more precise constrained motion.

Udwadia–Kalaba equation

To industrial manipulator without constraint, the general form of its dynamical equation can be shown as

M (q) . \overset{\cdot\cdot}{q} + C (q, \overset{\cdot}{q}) . \overset{\cdot}{q} + G (q) = τ

(1)

in which $q$ is the generalized coordinate n-vector. $M (q)$ is the $n \times n$ mass matrix. $C (q, \overset{\cdot}{q}) . \overset{\cdot}{q}$ is the $n \times 1$ centrifugal force and Coriolis forces vector. $G (q)$ is the gravitational vector of $n \times 1$ which is related to morpheme $q$ of manipulator. $τ$ is the joint torque vector of $n \times 1$ .

The uncontrolled system is now subjected to the $p$ constraints given by

\begin{matrix} φ_{j} (q, \overset{\cdot}{q}, t) = 0, & j = 1, 2, \dots, p \end{matrix}

(2)

These $p$ constraints include both geometric and trajectory requirements in this article. Differentiating equation (2) with respect to time twice yields the following constraint equation

A (q, \overset{\cdot}{q}, t) \overset{\cdot\cdot}{q} = b (q, \overset{\cdot}{q}, t)

(3)

in which $A (q, \overset{\cdot}{q}, t)$ is a $p \times n$ matrix and $b (q, \overset{\cdot}{q}, t)$ is a $p \times 1$ vector. $t$ is time variable. The dynamical equation of the constrained system is obtained as

M (q) . \overset{\cdot\cdot}{q} = Q (q, \overset{\cdot}{q}, t) + Q_{c} (q, \overset{\cdot}{q}, t)

(4)

in which $Q = τ - C (q, \overset{\cdot}{q}) . \overset{\cdot}{q} - G (q)$ . $Q_{c}$ is the control force n-vector that arises to ensure that the control requirements of the form of equation (3) are satisfied. Using the Udwadia–Kalaba equation,⁶ one can give the constraint force $Q_{c}$ explicitly as

\begin{array}{l} Q_{c} (q, \dot{q}, t) = M {(q)}^{\frac{1}{2}} B^{+} (q, t) \\ [b (q, \dot{q}, t) - A (q, \dot{q}, t) M^{- 1} (q) Q (q, \dot{q}, t)] \end{array}

(5)

in which $B (q, t) = A (q, \dot{q}, t) M^{- \frac{1}{2}} (q)$ , and the superscript “+” represents the Moore–Penrose generalized inverse.

That is to say, the dynamical equation of constraint to multi-body mechanical system is shown in the following general form

\begin{array}{l} M (q) . \overset{\cdot\cdot}{q} = Q (q, \dot{q}, t) + M {(q)}^{\frac{1}{2}} B^{+} (q, t) \\ [b (q, \dot{q}, t) - A (q, \dot{q}, t) M^{- 1} (q) Q (q, \dot{q}, t)] \end{array}

(6)

\begin{array}{l} \overset{\cdot\cdot}{q} = M^{- 1} (q) Q (q, \dot{q}, t) + M {(q)}^{- \frac{1}{2}} B^{+} (q, t) \\ [b (q, \dot{q}, t) - A (q, \dot{q}, t) M^{- 1} (q) Q (q, \dot{q}, t)] \end{array}

The standard Lagrange multiplier formulation of the problem of constrained motion can be shown as

[\begin{matrix} M (q) & - A {(q, \overset{\cdot}{q}, t)}^{T} \\ A (q, \overset{\cdot}{q}, t) & 0 \end{matrix}] [\begin{matrix} \overset{\cdot\cdot}{q} \\ λ \end{matrix}] = [\begin{matrix} Q (q, \overset{\cdot}{q}, t) \\ b (q, \overset{\cdot}{q}, t) \end{matrix}]

(7)

in which $λ$ is the Lagrange multiplier. Equation (7) can be shown as

\overset{\cdot\cdot}{q} = M^{- 1} Q + M^{- 1} A^{T} λ

(8)

A \overset{\cdot\cdot}{q} = b

(9)

Now, substituting equation (8) into the practical acceleration level constraints (9), one can solve for the Lagrange multipliers by direct inversion

λ = {(A M^{- 1} A^{T})}^{- 1} (b - A M^{- 1} Q)

(10)

Substituting equation (10) into equation (8)

\overset{\cdot\cdot}{q} = M^{- 1} Q + M^{- 1} A^{T} {(A M^{- 1} A^{T})}^{- 1} (b - A M^{- 1} Q)

(11)

in which $(A M^{- 1} A^{T})^{- 1} = (A (M^{\frac{1}{2}} M^{\frac{1}{2}})^{- 1} A^{T})^{- 1} = (A^{T})^{- 1} M^{\frac{1}{2}} (A M^{- \frac{1}{2}})^{- 1}$ . That is to say, equation (11) can be shown in the following form

\overset{\cdot\cdot}{q} = M^{- 1} Q + M^{- \frac{1}{2}} {(A M^{- \frac{1}{2}})}^{- 1} (b - A M^{- 1} Q)

(12)

According to the matrix theory, for the matrices $A$ and $M$ to be nonsingular, the following relationship can be set up

{(A M^{- \frac{1}{2}})}^{- 1} = {(A M^{- \frac{1}{2}})}^{+}

(13)

So equation (12) has the same form as equation (6). There are two points to highlight:

It is impossible to create a new fundamental principle for the theory of dynamic modeling of the constraint systems. Even though the fundamental theory is completely different (The standard Lagrange multiplier formulation is based on D’Alembert’s principle, and the Udwadia–Kalaba equation is based on Gauss’ rule), the dynamic equations have the same form on the premise that the matrices $A$ and $M$ are nonsingular. But the Udwadia–Kalaba equation provides a new perspective of understanding the laws of nature. So, in the process of the subsequent discussion about the computing error, it is not necessary to compare the Lagrange multiplier formulation and the Udwadia–Kalaba equation.

Comparing with the Lagrange multiplier formulation, Udwadia–Kalaba equation eliminates the computing complexity, which comes from the Lagrange multiplier, and enlarges the applied range of the dynamic modeling because of the application of the Moore–Penrose generalized inverse.

Dynamic modeling and simulation of the industrial robot subject to constraint

Dynamic modeling on manipulator without constraint

As shown in Figure 1, it is a schematic diagram of two-link robot manipulator moving on the vertical plane. The system generalized coordinate is $q = [\begin{matrix} q_{1} & q_{2} \end{matrix}]^{T}$ . It is assumed that lengths of each link are $l_{1}, l_{2}$ , qualities of each link are $m_{1}, m_{2}$ , torques of each joint are $τ_{1}, τ_{2}$ , and rotational inertias of each link are $I_{1}, I_{2}$ . $XOY$ is the global coordinate system.

Figure 1.

Sketch map of the two-link planar manipulator subject to constraint.

Let us establish the dynamical equation without constraint using Lagrange equation.

The total kinetic energy of system is

\begin{array}{l} E = \frac{1}{8} m_{1} l_{1}^{2} {\dot{q}}_{1}^{2} + \frac{1}{2} I_{1} {\dot{q}}_{1}^{2} + \frac{1}{2} m_{2} \\ [{(l_{1} {\dot{q}}_{1})}^{2} + {(\frac{1}{2} l_{2} {\dot{q}}_{2})}^{2} + l_{1} l_{2} \cos (q_{1} - q_{2}) {\dot{q}}_{1} {\dot{q}}_{2} + \frac{1}{2} I_{2} {\dot{q}}_{2}^{2}] \end{array}

(14)

The potential energy of system is

U = - \frac{1}{2} m_{1} l_{1} g \sin q_{1} - m_{2} g (l_{1} \sin q_{1} + \frac{1}{2} l_{2} \sin q_{2})

(15)

The dynamical equation of the system can be acquired by Lagrange equation

M (q) . \overset{\cdot\cdot}{q} + C (q, \overset{\cdot}{q}) . \overset{\cdot}{q} + G (q) = τ

(16)

in which

M (q) = [\begin{matrix} \frac{1}{4} m_{1} l_{1}^{2} + m_{2} l_{1}^{2} + I_{1} & \frac{1}{2} m_{2} l_{1} l_{2} \cos (q_{1} - q_{2}) \\ \frac{1}{2} m_{2} l_{1} l_{2} \cos (q_{1} - q_{2}) & \frac{1}{4} m_{2} l_{2}^{2} + I_{2} \end{matrix}]

C (q, \overset{\cdot}{q}) = [\begin{matrix} 0 & \frac{1}{2} m_{2} l_{1} l_{2} \sin (q_{1} - q_{2}) {\overset{\cdot}{q}}_{2} \\ \frac{1}{2} m_{2} l_{1} l_{2} \sin (q_{2} - q_{1}) {\overset{\cdot}{q}}_{1} & 0 \end{matrix}]

G (q) = [\begin{matrix} (\frac{1}{2} m_{1} + m_{2}) g l_{1} \cos q_{1} \\ \frac{1}{2} m_{2} g l_{2} \cos q_{2} \end{matrix}]

τ = {[\begin{matrix} τ_{1} & τ_{2} \end{matrix}]}^{T}

The dynamical equation of the system without constraint is shown as

M (q) \overset{\cdot\cdot}{q} = Q

(17)

in which

Q = τ - C (q, \overset{\cdot}{q}) \overset{\cdot}{q} - G (q)

Because the system has not been constrained by now, $τ = [\begin{matrix} 0 & 0 \end{matrix}]^{T}$ .

Then, we have

Q = - C (q, \overset{\cdot}{q}) \overset{\cdot}{q} - G (q)

Dynamics modeling subject to constraint

The trajectory constraint of the end-effector is summarized as

X_{B} = L

Y_{B} = 0.5 (1 - \cos t)

(18)

in which $[\begin{matrix} X_{B} & Y_{B} \end{matrix}]$ is the coordinate of the end-effector in the global coordinate system. As shown in Figure 1, one also has the following equation

l_{1} \cos q_{1} + l_{2} \cos q_{2} = X_{B}

l_{1} \sin q_{1} + l_{2} \sin q_{2} = Y_{B}

(19)

Differentiating equation (19) with respect to time twice

- l_{1} \sin q_{1} {\overset{\cdot\cdot}{q}}_{1} - l_{2} \sin q_{2} {\overset{\cdot\cdot}{q}}_{2} = l_{1} \cos q_{1} {\overset{\cdot}{q}}_{1}^{2} + l_{2} \cos q_{2} {\overset{\cdot}{q}}_{2}^{2} + {\overset{\cdot\cdot}{X}}_{B}

l_{1} \cos q_{1} {\overset{\cdot\cdot}{q}}_{1} + l_{2} \cos q_{2} {\overset{\cdot\cdot}{q}}_{2} = l_{1} \sin q_{1} {\overset{\cdot}{q}}_{1}^{2} + l_{2} \sin q_{2} {\overset{\cdot}{q}}_{2}^{2} + {\overset{\cdot\cdot}{Y}}_{B}

(20)

Differentiating equation (18) with respect to time twice

{\overset{\cdot\cdot}{X}}_{B} = 0

{\overset{\cdot\cdot}{Y}}_{B} = 0.5 \cos t

(21)

According to equations (20) and (21), one can have the following constraint equation

A \overset{\cdot\cdot}{q} = b

(22)

in which

A = [\begin{matrix} - l_{1} \sin q_{1} & - l_{2} \sin q_{2} \\ l_{1} \cos q_{1} & l_{2} \cos q_{2} \end{matrix}]

b = [\begin{matrix} l_{1} \cos q_{1} {\overset{\cdot}{q}}_{1}^{2} + l_{2} \cos q_{2} {\overset{\cdot}{q}}_{2}^{2} \\ l_{1} \sin q_{1} {\overset{\cdot}{q}}_{1}^{2} + l_{2} \sin q_{2} {\overset{\cdot}{q}}_{2}^{2} + 0.5 \cos t \end{matrix}]

Then, the additional torques for system to satisfy the given trajectory is explicitly determined by equation (5) and is shown as

Q_{c} = M^{\frac{1}{2}} {(A M^{- \frac{1}{2}})}^{+} (b - A M^{- 1} Q)

(23)

Finally, the dynamic equation is obtained by equations (4), (5), and (23)

M \overset{\cdot\cdot}{q} = Q + Q_{c}

(24)

Result and analysis to simulation

It is assumed that $m_{1} = m_{2} = 1 kg$ , $l_{1} = l_{2} = 1 m$ , $I_{1} = I_{2} = (1 / 3) kg m^{2}$ , $q_{1} = 60 \circ$ , $q_{2} = - 60 \circ$ , ${\overset{\cdot}{q}}_{1} = {\overset{\cdot}{q}}_{2} = 0 m / s$ , and $L = 1 m$ . The simulation time is 20 s. The initial configuration is shown in Figure 2.

Figure 2.

Sketch map of the initial configuration of the two-link planar manipulator subject to constraint.

The theoretical curve of joint angles and joint angular velocities can be obtained by solving equation (19). It can be regarded as reference curve for the following study. Figure 3 represents the joint angles at each joint of the robotic system subjected to the constraints of the vertical plane while the end-effector stays at the initial position, in which the solid line represents the theoretical value and the dashed line represents the numerical integration value. It is obvious that the dashed line approximates the solid line. The first link moves from the initial configuration to 90° and then comes back to the initial position. The second link moves from the initial configuration to the horizontal position and then comes back to the initial position. So, the link’s movement is periodic which leads to the periodic moving of the end-effector from the lowest position to the highest position and then back to the lowest position. However, obvious errors between the solid line and the dashed line still exist. Figure 4 represents the errors between the solid line and the dashed line. It is observed that the errors gradually increase with time. Figure 5 represents the joint angular velocities. Figure 6 represents the angular velocities errors between the theoretical value and the numerical integration value. It is also observed that the errors gradually increase with time. But the trend is gentle than the joint angular errors. Thus, it is necessary to consider the numerical method for reducing the errors; this is handled in the following section.

Figure 3.

Sketch map of the joint angles.

Figure 4.

Errors of the joint angles.

Figure 5.

Sketch map of the joint angular velocities.

Figure 6.

Errors of the joint angular velocities.

Numerical method to reduce errors

It is well known that numerical integration scheme on the second-order differential equation needs a set of two first-order differential equations. The second-order differential equation of equation (6) is expressed as a set of the first-order differential equations as

\overset{\cdot}{q} = V

(25a)

\overset{\cdot}{V} = M^{- 1} Q + M^{- \frac{1}{2}} B^{+} (b - A M^{- 1} Q)

(25b)

However, it is found that the numerical results show unneglectable errors between the theoretical curve and the numerical integration curve as shown in the previous example. This is because the dynamic equation utilizes equation (3) and neglects the velocity-based constraint that differentiates the constraint equations once with respect to time. It can be written as

A (q, \overset{\cdot}{q}, t) \overset{\cdot}{q} = b_{q} (q, \overset{\cdot}{q}, t)

(26)

Then, the dynamic equation can be modified as

\overset{\cdot}{q} = V + M^{- \frac{1}{2}} B^{+} (b_{q} - AV)

(27a)

\overset{\cdot}{V} = M^{- 1} Q + M^{- \frac{1}{2}} B^{+} (b - A M^{- 1} Q)

(27b)

It is important to note that the new term $M^{- \frac{1}{2}} B^{+} (b_{q} - AV)$ is not generated by the constraints but rather is introduced to compensate for numerical errors along the integration. Thus, it has no clear physical meaning and can only be considered as a kinematic correction term.

Differentiating equation (19) with respect to time once

- l_{1} \sin q_{1} {\overset{\cdot}{q}}_{1} - l_{2} \sin q_{2} {\overset{\cdot}{q}}_{2} = 0

l_{1} \cos q_{1} {\overset{\cdot}{q}}_{1} + l_{2} \sin q_{2} {\overset{\cdot}{q}}_{2} = 0.5 \sin t

(28)

According to equation (26), one can have the following constraint equation

A \overset{\cdot}{q} = b_{q}

(29)

in which

A = [\begin{matrix} - l_{1} \sin q_{1} & - l_{2} \sin q_{2} \\ l_{1} \cos q_{1} & l_{2} \cos q_{2} \end{matrix}]

b_{q} = [\begin{matrix} 0 \\ 0.5 \sin t \end{matrix}]

Substituting the dynamic equation of the unconstrained system and the constraint equations into equation (27), the dynamic responses and constraint forces can be obtained. Figure 7 represents the joint angles at each joint of the robotic system subjected to the constraints of the vertical plane depending on the modified numerical integration method, in which the solid line shows the theoretical curve and the dashed line shows the modified curve. Figure 8 represents the errors between the theoretical curve and the modified curve. It is observed that the errors’ growth trend is improved obviously and decreased by one order of magnitude. Figure 9 represents the joint angular velocities, in which the solid line shows the theoretical curve and the dashed line shows the modified curve. Figure 10 represents the errors between the theoretical curve and the modified curve. It is observed that the errors also decreased by one order of magnitude. Figure 11 compares the errors of the constraints depending on the application of the proposed numerical scheme, where Figure 11(a) and (c) represents the errors in the satisfaction of equation (19) and Figure 11(b) and (d) represents the errors in the first time derivative of equation (19). It is observed that the errors obtained by the numerical integration of equation (25) gradually increase with time to veer away from the constrained path. The proposed numerical scheme describes more precise dynamic responses; still, the errors do not perfectly disappear. Figure 12 represents the constraint torques to be required at each joint for its movement.

Figure 7.

Sketch map of the joint angles depending on the modified numerical integration method.

Figure 8.

Errors of the joint angles between the theoretical curve and modified curve.

Figure 9.

Sketch map of the joint angular velocities depending on the modified numerical integration method.

Figure 10.

Errors of the joint angular velocities between the theoretical curve and modified curve.

Figure 11.

Errors of constraints: (a) constraint error of X direction, (b) the first time derivative of constraint error of X direction, (c) constraint error of Y direction, and (d) the first time derivative of constraint error of Y direction.

Figure 12.

Sketch map of the joint torque.

Conclusion and future works

Aiming at the dynamics modeling of industrial robot subject to constraint, this article establishes the dynamic equation based on the Udwadia–Kalaba equation and proposes a modified algorithm to reduce the error caused by numerical integration. The research process reveals some conclusions as follows:

The geometric constraint and the trajectory of the end-effector are considered as the external constraint of the system, and the constrained relationship of the system is integrated into the dynamic equation dexterously based on the idea of the Udwadia–Kalaba equation. The simple method overcomes the disadvantage of obtaining dynamical equation from traditional Lagrange equation by Lagrange multiplier effectively.

The velocity-based constraint is considered in dynamics modeling process to reduce the error. The numerical simulation results prove that the dynamical equation established by this method decreases the joint angle and joint angular velocity errors by one order of magnitude and conforms to the matter of fact precisely.

Future work mainly consists of two tasks:

However, the simulation is executed on the condition which the initial condition satisfies the constrained equation, then the next research direction is the dynamics modeling of the industrial robot subject to constraint which the initial condition does not satisfy the constrained equation and expanding to the three-dimensional industrial robot.

The error is improved to a certain degree, but not perfect. So, the next research is to reduce the error by proper control strategy.

Footnotes

Academic Editor: Chuanzeng Zhang

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Korayem

Azimirad

Nikoobin

. Maximum load-carrying capacity of autonomous mobile manipulator in an environment with obstacle considering tip over stability. Int J Adv Manuf Tech 2010; 46: 811–830.

Sanchez-Sanchez

Arteaga-Perez

MA.

Simplied methodology for obtaining the dynamic model of robot manipulators. Int J Adv Robot Syst 2012; 9: 700–709.

Kane

TR.

Dynamics of non-holonomic systems. J Appl Mech 1961; 28: 574–578.

Jia

Y-H

S-J.

Dynamics and adaptive control of a dual-arm space robot with closed-loop constraints and uncertain inertial parameters. Acta Mech Sinica 2014; 1: 112–124.

Saha

Shah

Nandihal

PV.

Evolution of the DeNOC-based dynamic modelling for multibody systems. Mech Sci 2013; 4: 1–20.

Udwadia

Kalaba

REA

. A new perspective on constrained motion. P Roy Soc Lond A Mat 1992; 439: 407–410.

Udwadia

Schutte

AD.

Equations of motion for general constrained systems in Lagrangian mechanics. Acta Mech 2010; 213: 111–129.

Udwadia

Kalaba

RE.

Explicit equations of motion for mechanical systems with nonideal constraints. J Appl Mech 2001; 68: 462–467.

Udwadia

Kalaba

RE.

On constrained motion. Appl Math Comput 2005; 164: 313–320.

10.

Udwadia

Wanichanon

Control of uncertain nonlinear multi-body mechanical systems. J Appl Mech 2013; 81: 1–11.

11.

Pennestri

Valentini

Falco

DD.

An application of the Udwadia–Kalaba dynamic formulation to flexible multibody systems. J Frankl Inst 2010; 347: 173–194.

12.

Schutte

Udwadia

FE.

New approach to the modeling of complex multi-body dynamical systems. J Appl Mech 2010; 78: 856–875.

13.

Cho

New approach to satellite formation-keeping: exact solution to the full nonlinear problem. J Aerospace Eng 2009; 22: 445–455.

14.

Zhen

Zhao

Huang

. On Kepler’s law: application of the Udwadia–Kalaba theory. Sci Sin-Phys Mech Astron 2014; 44: 24–31.

15.

Zhen

Huang

Zhao

. Why can a free-falling cat always manage to land safely on its feet?Nonlinear Dynam 2015; 79: 1–14.

16.

Huang

Chen

Zhong

ZH.

Udwadia–Kalaba approach for parallel manipulator dynamics. J Dyn Syst: T ASME 2013; 135: 1012–1030.

17.

Zhao

Zhen

Chen

YH.

Dynamic modeling and simulation of multi-body systems using the Udwadia–Kalaba theory. Chin J Mech Eng 2013; 26: 839–850.