Decomposition and Cross-Product-Based Method for Computing the Dynamic Equation of Robots

Abstract

This paper aims to demonstrate a clear relationship between Lagrange equations and Newton-Euler equations regarding computational methods for robot dynamics, from which we derive a systematic method for using either symbolic or on-line numerical computations. Based on the decomposition approach and cross-product operation, a computing method for robot dynamics can be easily developed. The advantages of this computing framework are that: it can be used for both symbolic and on-line numeric computation purposes, and it can also be applied to biped systems, as well as some simple closed-chain robot systems.

1. Introduction

When looking at robotics research, dynamic modelling is basically a form of artwork. The complexity of a dynamic model concerns the degrees of freedom (DOFs) of a machine or a robot and the dynamic model of a low DOFs robot can be easily developed, but dynamic models for a high DOFs robot are difficult to derive by hand. The drawbacks to creating the models by hand are: difficulty of maintenance, complexity of proofreading the dynamic model and inflexibility. The biped robot is the target of a complex dynamic system study [1] that also reveals a demand for a computer-aided systematic derivation tools for robotics.

Using an auto-deriving tool to obtain a dynamic model of a robot is very economical. Tošic [2] proposed the benefits of a symbolic simulation of an engineering system, while Cetinkunt [3] proposed a symbolic modelling of the dynamics of robotic manipulators on the numerical tool REDUCE. Since the coding of REDUCE is based on C and Fortran, modifying the system model by simply changing the symbol is rather difficult. Vibet demonstrated the symbolic tool FORM [4], which automatically generates the dynamic equations of the manipulator. The MACSYMA tool was then introduced to mechanical engineering education [5], providing the steps for determining a system dynamics model, but it does not support online simulations that integrate both system models and control rules.

Numerous contributions in both algorithms and their computational efficiency have been made in the field of robot dynamics [6–9]. The need for an online computing [6] method is motivated by simulation studies of the dynamics of a robotic system with large DOFs. While computation efficiency continues to be crucial for the simulation and control of increasingly complex mechanisms operating at higher speeds, other aspects of the dynamics problem are also vital. Algorithms should be formulated under a simple framework to enable ease of development, modification, extension and understanding.

This paper aims to demonstrate a clear relationship between Lagrange equations and Newton-Euler equations regarding the computational methods for robot dynamics, from which we derive a systematic method using either symbolic or on-line numerical computations. First, this study reviews the Lagrange formulation, as well as its symbolic generation method, which is suitable for studying a low-DOF robotic system. Second, based on the decomposition approach and cross-product operation, a computing framework can be easily developed to calculate the inertia matrix, the Coriolis and centrifugal matrix, and the gravity force vector of robot dynamics equations. The advantages of this computing method are that: it can be used for both symbolic and the on-line numeric computation purposes, and it also can be applied to a closed-chain robot system. We present several robotic system examples to demonstrate the versatility of the proposed method.

2. Decomposition method

Assume that a robot system has a generalized coordinate with n DOFs and let q = [q₁ q₂ ṫ q_n]^T. The equation of robot motion [9] can be written as

$D(q)\stackrel{..}{q}+C(q,\stackrel{.}{q})\stackrel{.}{q}+G(q)=u,$

(1)

where the inertia matrix D is an n×n positive-definite symmetric matrix, C is an n×n Coriolis and centrifugal matrix, and G is an n×1 gravity force vector. The control input u = u₁ u₂ ṫ u_n] ^T is the vector of joint input forces or torques, and u_i ≡ 0, if the i-th joint is not actuated. Significantly, C is not unique, and S = Ḋ − 2C is a skew-symmetric matrix caused by the null energy property.

2.1 Lagrange formulation approach

The total kinetic energy K of a robot consists of the rotational kinetic energy K_ωi and the translational kinetic energy K_vi of each link. The total potential energy P consists of the potential energy P_i of each link; and therefore,

$K={\displaystyle \sum _{i=1}^n{K}_{\omega i}}+{\displaystyle \sum _{i=1}^n{K}_{vi}},$

(2)

and

$P={\displaystyle \sum _{i=1}^n{P}_i}.$

(3)

The kinetic energies K_ωi and K_vi, and potential energy P_i of each link can be calculated as stated below. Let m_i be the mass of link i and I_i be the inertia of link i expressed in the body-attached frame. Let R_i be the rotational matrix, ω_i be the angular velocity and p_ci be the CoM (Centre of Mass) of link i. The angular velocity of link i can thus be expressed as

${\omega }_i=J_{\omega i}\stackrel{.}{q},$

(4)

where $J_{\omega i}=\frac{\partial {\omega }_i}{\partial \stackrel{.}{q}}$ is the angular Jacobian matrix of link i, and the rotational kinetic energy of link i can be obtained through

$K_{\omega i}=\frac{1}{2}{\omega }_i^T\left(R_i{I}_i{R}_i^T\right){\omega }_i=\frac{1}{2}{\stackrel{.}{q}}^T\left(J_{\omega i}^T{R}_i{I}_i{R}_i^T{J}_{\omega i}\right)\stackrel{.}{q}.$

(5)

The potential energy of link i can be obtained by

$P_i=m_{i}g{p}_{ci}^T{e}_3,$

(6)

where g is the acceleration of gravity and e₃ = [0 0 1]^T. The velocity of the CoM of link i can be computed by

${\stackrel{.}{p}}_{ci}=J_{vi}\stackrel{.}{q},$

(7)

where $J_{vi}=\frac{\partial p_{ci}}{\partial q}$ is the linear Jacobian matrix of link i, and the translational energy of link i is via

$K_{vi}=\frac{1}{2}{m}_i{\stackrel{.}{p}}_{ci}^T{\stackrel{.}{p}}_{ci}=\frac{1}{2}{\stackrel{.}{q}}^T\left(m_i{J}_{vi}^T{J}_{vi}\right)\stackrel{.}{q}$

(8)

From the Lagrange formulation, we have

$\begin{array}{c}{D}_{\omega i}=\frac{\partial }{\partial \stackrel{.}{q}}\frac{\partial K_{\omega i}}{\partial \stackrel{.}{q}}\\ C_{\omega i}=\frac{\partial }{\partial q}\frac{\partial K_{\omega i}}{\partial \stackrel{.}{q}}-\frac{1}{2}\frac{\partial }{\partial q}\left(D_{\omega i}\stackrel{.}{q}\right)\end{array}$

and

$\begin{array}{c}{D}_{vi}=\frac{\partial }{\partial \stackrel{.}{q}}\frac{\partial K_{vi}}{\partial \stackrel{.}{q}}\\ C_{vi}=\frac{\partial }{\partial q}\frac{\partial K_{vi}}{\partial \stackrel{.}{q}}-\frac{1}{2}\frac{\partial }{\partial q}\left(D_{vi}\stackrel{.}{q}\right)\\ G_i=\frac{\partial P_i}{\partial q}.\end{array}$

The D, C and G matrices of the robot's equations of motion are then

$D={\displaystyle \sum _{i=1}^n{D}_{\omega i}}+{\displaystyle \sum _{i=1}^n{D}_{vi}}$

(9)

$C={\displaystyle \sum _{i=1}^n{C}_{\omega i}}+{\displaystyle \sum _{i=1}^n{C}_{vi}}$

(10)

and

$G={\displaystyle \sum _{i=1}^n{G}_i}.$

(11)

2.2 Newton-Euler formulation approach

The decomposition of the robot's equations of motion can also be conducted using the Newton-Euler approach. The angular acceleration of link i is expressed as

${\stackrel{.}{\omega }}_i=J_{\omega i}\stackrel{..}{q}+{\stackrel{.}{J}}_{\omega i}\stackrel{.}{q},$

(12)

where ${\stackrel{.}{J}}_{\omega i} =\frac{\partial }{\partial q}\left(J_{\omega i} \stackrel{.}{q}\right)$ is the differentiation of the angular Jacobian matrix J_ωi. The Euler equation of link i is

$R_i{I}_i{R}_i^T{\stackrel{.}{\omega }}_i+{\omega }_i\times \left(R_i{I}_i{R}_i^T\right){\omega }_i={\tau }_i,$

(13)

and hence,

$R_i{I}_i{R}_i^T\left(J_{\omega i}\stackrel{..}{q}+{\stackrel{.}{J}}_{\omega i}\stackrel{.}{q}\right)+{\omega }_i\times \left(R_i{I}_i{R}_i^T\right){\omega }_i={\tau }_i.$

Since $J_{\omega i}^T {\tau }_i=u$ , we now have the equation

$\begin{array}{l}{J}_{\omega i}^T{R}_i{I}_i{R}_i^T{J}_{\omega i}\stackrel{..}{q}+J_{\omega i}^T{R}_i{I}_i{R}_i^T{\stackrel{.}{J}}_{\omega i}\stackrel{.}{q}+\\ +J_{\omega i}^T{\omega }_i\times \left(R_i{I}_i{R}_i^T\right)J_{\omega i}\stackrel{.}{q}=J_{\omega i}^T{\tau }_i=u\prime \end{array}$

which is the same as

$D_{\omega i}\stackrel{..}{q}+C_{\omega i}\stackrel{.}{q}=u.$

Therefore, D_ωi and C_ωi can be obtained respectively using:

$D_{\omega i}=J_{\omega i}^T\left(R_i{I}_i{R}_i^T\right)J_{\omega i}$

(14)

and

$C_{\omega i}=J_{\omega i}^T\left(R_i{I}_i{R}_i^T\right){\stackrel{.}{J}}_{\omega i}+J_{\omega i}^T\left[{\omega }_i\times \right]\left(R_i{I}_i{R}_i^T\right)J_{\omega i},$

(15)

where $\left[\begin{array}{c}x\\ y\\ z\end{array}\right]\times \equiv \left[\begin{array}{ccc}0& -z& y\\ z& 0& -x\\ -y& x& 0\end{array}\right]$ is the skew-symmetric matrix representation of the cross-product operation.

The acceleration of CoM of link i is

${\stackrel{..}{p}}_{ci}=J_{vi}\stackrel{..}{q}+{\stackrel{.}{J}}_{vi}\stackrel{.}{q},$

(16)

where ${\stackrel{.}{J}}_{vi}=\frac{\partial }{\partial q}\left(J_{vi}\stackrel{.}{q}\right)$ is the differentiation of the linear Jacobian matrix J_vi. The Newton equation of the CoM of link i is

$m_i{\stackrel{..}{p}}_{ci}+m_{i}g{e}_3=f_i$

(17)

and hence,

$m_i\left(J_{vi}\stackrel{..}{q}+{\stackrel{.}{J}}_{vi}\stackrel{.}{q}\right)+m_{i}g{e}_3=f_i.$

Since $J_{vi}^T f_i=u$ , we formulate the equation

$m_i{J}_{vi}^T{J}_{vi}\stackrel{..}{q}+m_i{J}_{vi}^T{\stackrel{.}{J}}_{vi}\stackrel{.}{q}+m_{i}g{J}_{vi}^T{e}_3=J_{vi}^T{f}_i=u,$

which is also the same as

$D_{vi}\stackrel{..}{q}+C_{vi}\stackrel{.}{q}+G_i=u.$

Therefore, the D_vi, C_vi and G_i matrices can be obtained by

$D_{vi}=m_i{J}_{vi}^T{J}_{vi}$

(18)

$C_{vi}=m_i{J}_{vi}^T{\stackrel{.}{J}}_{vi}$

(19)

and

$G_i=\frac{\partial P_i}{\partial q}.$

(20)

2.3 Extended coordinate system

The Cartesian coordinates of the base point, [x₀ y₀ z₀] ^T , are frequently added to form an extended set of generalized coordinates for the robot systems. Let q_e = [x₀ y₀ z₀ q₁ ṫ q_n] ^T be the extended coordinate system and next define the total mass $M={\displaystyle \sum _{i=1}^n{m}_i}$ , Jacobian matrix $J_{CoM}=\frac{1}{M}{\displaystyle \sum _{i=1}^n{m}_i}{J}_{vi}$ , and its derivative ${\stackrel{.}{J}}_{CoM}=\frac{1}{M}{\displaystyle \sum _{i=1}^n{m}_i}{\stackrel{.}{J}}_{vi}$ . The robot's equations of motion in this extended coordinate system can thus be expressed as

$D_e(q_e){\stackrel{..}{q}}_e+C_e(q_e,{\stackrel{.}{q}}_e){\stackrel{.}{q}}_e+G_e(q_e)=B_{e}u,$

(21)

where

$\begin{array}{l}{D}_e=\left[\begin{array}{cc}M{I}_{3\times 3}& M{J}_{CoM}^{}\\ M{J}_{CoM}^T& D\end{array}\right]\\ C_e=\left[\begin{array}{cc}{0}_{3\times 3}& M{\stackrel{.}{J}}_{CoM}\\ 0_{n\times 3}& C\end{array}\right]\\ G_e=\left[\begin{array}{c}M{e}_3\\ G\end{array}\right]\end{array}$

and

$B_e=\left[\begin{array}{c}{0}_{3\times n}\\ I_{n\times n}\end{array}\right].$

3. Procedures for computing the dynamic equations of a serial-link robot

Let m_i be the mass of link i, I_i be the inertia of link i expressed in the body-attached frame and then q = [q₁ q₂ ṫ q_n] ^T be the joint coordinate vector, where n is the number of joints. For simplicity, it is assumed that the relative motion of link i with respect to link (i-1) is either the rotational (revolute joint) or the translational (prismatic joint) motion with one of the basic axes in the local body frame from the set {± x-axis, ± y-axis, and ± z-axis} and is denoted by symbol ρ_i, where ρ_i ∈ {+Rx, −Rx, +Ry, −Ry, +Rz, −Rz, +Tx, −Tx, +Ty, −Ty, −Tz, −Tz}.

Define three basic rotational matrices: R_x(q_i), R_y(q_i), and R_z(q_i), as shown below

$\begin{array}{l}{R}_x(q_i)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos q_i& -\sin q_i\\ 0& \sin q_i& \cos q_i\end{array}\right]\\ R_y(q_i)=\left[\begin{array}{ccc}\cos q_i& 0& \sin q_i\\ 0& 1& 0\\ -\sin q_i& 0& \cos q_i\end{array}\right]\end{array}$

and

$R_z(q_i)=\left[\begin{array}{ccc}\cos q_i& -\sin q_i& 0\\ \sin q_i& \cos q_i& 0\\ 0& 0& 1\end{array}\right].$

The procedures for computing the dynamic equation of a serial-link robot are thus as follows:

Step 1: Computing the rotational matrix of link i, R_i, i = 1, 2, ṫ, n

Let R₀ = I_3×3, and

$R_i=\{\begin{array}{cc}{R}_{i-1}{R}_x(\pm q_i)& {\rho }_i=\pm Rx\\ R_{i-1}{R}_y(\pm q_i)& {\rho }_i=\pm Ry\\ R_{i-1}{R}_z(\pm q_i)& {\rho }_i=\pm Rz\\ R_{i-1}& \text{prismatic\hspace{0.17em}joint}\end{array},i=1,2,\cdots ,n.$

(22)

Let a_i and t_i be the vectors of the corresponding rotational and translational axes of link i with respect to link (i-1), and then

$a_i=\{\begin{array}{cc}\pm R_{i-1}(:,1)& {\rho }_i=\pm Rx\\ \pm R_{i-1}(:,2)& {\rho }_i=\pm Ry\\ \pm R_{i-1}(:,3)& {\rho }_i=\pm Rz\\ 0_{3\times 1}& \text{prismatic\hspace{0.17em}joint}\end{array}$

(23)

and

$t_i=\{\begin{array}{cc}\pm R_{i-1}(:,1)& {\rho }_i=\pm Tx\\ \pm R_{i-1}(:,2)& {\rho }_i=\pm Ty\\ \pm R_{i-1}(:,3)& {\rho }_i=\pm Tz\\ 0_{3\times 1}& \text{revolute\hspace{0.17em}joint}\end{array},i=1,2,\cdots ,n,$

(24)

where R(:, k), k = 1,2,3 represents the k-th column of the matrix R.

Step 2: Computing the angular Jacobian matrix of link i, ω_i = J _ωiq̇, i = 1,2, ṫ, n

Let ω₀ = 0_3×1, and

${\omega }_i={\omega }_{i-1}+a_i{\stackrel{.}{q}}_i,i=1,2,\cdots ,n.$

(25)

The angular Jacobian matrix of link i, J_ωi, is

$\begin{array}{ll}{J}_{\omega i}(:,j)=\{\begin{array}{ll}{a}_j\hfill & j\le i\hfill \\ 0_{3x1}\hfill & j>i\hfill \end{array}\hfill & \hfill \\ \hfill & ,\hspace{0.17em}i,\hspace{0.17em}j=1,2,...,n.\hfill \end{array}$

(26)

where J_ωi (:, j) represents the j-th column of matrix J_ωi and the total rotational kinetic energy K _ω is obtained from

$\begin{array}{l}{K}_{\omega }={\displaystyle \sum _{i=1}^n{K}_{\omega i}}={\displaystyle \sum _{i=1}^n\frac{1}{2}{\omega }_i^T{R}_i{I}_i{R}_i^T{\omega }_i}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\frac{1}{2}{\stackrel{.}{q}}^T\left({\displaystyle \sum _{i=1}^n{J}_{\omega i}^T{R}_i{I}_i{R}_i^T{J}_{\omega i}}\right)\stackrel{.}{q}.\end{array}$

(27)

Matrix D_ωi is formed by

$D_{\omega i}=J_{\omega i}^T{R}_i{I}_i{R}_i^T{J}_{\omega i},i=1,2,\cdots ,n.$

Step 3: Differentiation of the angular Jacobian matrices The angular acceleration of link i is obtained as

${\stackrel{.}{\omega }}_i=J_{\omega i}\stackrel{..}{q}+{\stackrel{.}{J}}_{\omega i}\stackrel{.}{q}i=1,2,\cdots ,n,$

where ${\stackrel{.}{J}}_{\omega i}=\frac{\partial }{\partial q}\left(J_{\omega i}\stackrel{.}{q}\right)={\displaystyle \sum _{j=1}^n\frac{\partial J_{\omega i}}{\partial q_j}{\stackrel{.}{q}}_j}$ is found to be

${\stackrel{.}{J}}_{\omega i}(:,j)=\{\begin{array}{cc}{\displaystyle \sum _{k=1}^{j-1}{a}_k\times a_j{\stackrel{.}{q}}_k}& j\le i\\ 0_{3\times 1}& j>i\end{array},i,j=1,2,\cdots ,n.$

(28)

Matrix C_ωi is set up by

$C_{\omega i}=J_{\omega i}^T{R}_i{I}_i{R}_i^T{\stackrel{.}{J}}_{\omega i}+J_{\omega i}^T\left[{\omega }_i\times \right]R_i{I}_i{R}_i^T{J}_{\omega i},$

Step 4: Computing the linear Jacobian matrix of the CoM of link i, ṗ_ci = J_viq̇, i = 1,2, ṫ, n

First, compute the joint position and CoM position of link i, p_i and p_ci, i = 1, 2, ṫ, n. The linear Jacobian matrix of the CoM of link i, $J_{vi}=\frac{\partial p_{ci}}{\partial q}$ , is then found as

$J_{vi}(:,j)=\{\begin{array}{cc}\{\begin{array}{cc}{a}_j\times (p_{ci}-p_j)& \text{revolute\hspace{0.17em}joint}\\ t_j& \text{prismatic\hspace{0.17em}joint}\end{array}& j\le i\\ 0_{3\times 1}& j>i\end{array}\begin{array}{l},i,j=1,2,\cdots ,n.\end{array}$

(29)

The total translational kinetic energy K_v and the total potential energy P are obtained as

$\begin{array}{l}{K}_v={\displaystyle \sum _{i=1}^n{K}_{vi}}=\frac{1}{2}{\displaystyle \sum _{i=1}^n{m}_i{\stackrel{.}{p}}_{ci}^T}{\stackrel{.}{p}}_{ci}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\frac{1}{2}{\stackrel{.}{q}}^T\left({\displaystyle \sum _{i=1}^n{m}_i{J}_{vi}^T}{J}_{vi}\right)\stackrel{.}{q}\end{array}$

(30)

and

$P={\displaystyle \sum _{i=1}^n{P}_i={\displaystyle \sum _{i=1}^n{m}_{i}g{p}_{ci}^T{e}_3}}.$

(31)

The matrices D_vi and G_i are set up as

$D_{vi}=m_i{J}_{vi}^T{J}_{vi},$

and

$G_i=m_{i}g{J}_{vi}^T{e}_3,i=1,2,\cdots ,n.$

Step 5: Differentiation of the linear Jacobian matrices The linear acceleration of the CoM of link i is obtained as

${\stackrel{..}{p}}_{ci}=J_{vi}\stackrel{..}{q}+\frac{\partial }{\partial q}\left(J_{vi}\stackrel{.}{q}\right)\stackrel{.}{q}=J_{vi}\stackrel{..}{q}+{\stackrel{.}{J}}_{vi}\stackrel{.}{q},$

where ${\stackrel{.}{J}}_{vi}=\frac{\partial }{\partial q}\left(J_{vi}\stackrel{.}{q}\right)={\displaystyle \sum _{j=1}^n\frac{\partial J_{vi}}{\partial q_j}{\stackrel{.}{q}}_j}$ is found to be

${\stackrel{.}{J}}_{vi}(:,j)=\{\begin{array}{cc}\{\begin{array}{cc}{\displaystyle \sum _{k=1}^{j-1}{a}_k\times J_{vi}(:,j){\stackrel{.}{q}}_k}+{\displaystyle \sum _{k=j}^i{a}_j\times J_{vi}(:,k){\stackrel{.}{q}}_k}& \text{revolute\hspace{0.17em}joint}\\ {\displaystyle \sum _{k=1}^{j-1}{a}_k\times J_{vi}(:,j){\stackrel{.}{q}}_k}& \text{prismatic\hspace{0.17em}joint}\end{array}& j\le i\\ 0_{3\times 1}& j>i\end{array},i,j=1,2,\cdots ,n$

(32)

The matrix C_vi is equal to

$C_{vi}=m_i{J}_{vi}^T{\stackrel{.}{J}}_{vi},i=1,2,\cdots ,n.$

Step 6: Putting it all together

The D, C, and G matrices of the robot's equations of motion are

$D={\displaystyle \sum _{i=1}^n{D}_{\omega i}}+{\displaystyle \sum _{i=1}^n{D}_{vi}}={\displaystyle \sum _{i=1}^n{J}_{\omega i}^T{R}_i{I}_i{R}_i^T}{J}_{\omega i}+{\displaystyle \sum _{i=1}^n{m}_i{J}_{vi}^T}{J}_{vi}$

(33)

$\begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}C={\displaystyle \sum _{i=1}^n{C}_{\omega i}}+{\displaystyle \sum _{i=1}^n{C}_{vi}},\\ ={\displaystyle \sum _{i=1}^n\left(J_{\omega i}^T{R}_i{I}_i{R}_i^T{\stackrel{.}{J}}_{\omega i}+J_{\omega i}^T\left[{\omega }_i\times \right]R_i{I}_i{R}_i^T{J}_{\omega i}\right)}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\displaystyle \sum _{i=1}^n{m}_i{J}_{vi}^T}{\stackrel{.}{J}}_{vi}\end{array}$

(34)

and

$G={\displaystyle \sum _{i=1}^n{G}_i}={\displaystyle \sum _{i=1}^n{m}_i}g{J}_{vi}^T{e}_3.$

(35)

4. Illustrated examples

4.1 Example 1: A 3 DOFs spherical robot

Figure 1 shows an example of a simple 3-DOF spherical robot system with both revolute and prismatic joints, and generalized coordinate q = [q₁ q₂ r] ^T , and is presented to verify the proposed computing method. The system parameters of the spherical robot are shown in Table 1 and the control input vector u = [u₁ u₂ u₃]^T. For a link with two joints, a virtual link with zero mass and/or zero inertia can be inserted for easy programming. The proposed method and the Lagrange method are shown to obtain the same equations of motion as follows

Table 1.

System parameters of the 3-DOF spherical robot system.

I	Joint ρi	Link mass mi	Link length ri	Link CoM di	Link inertia Ii
1	+Rx	0	[0 0 0] T	[0 0 0] T	zero inertia
2	+Ry	0	[0 0 0] T	[0 0 0] T	diag{Ix, Iy, Iz}
3	+Tz	M	[0 0 r] T	[0 0 r] T	zero inertia

Figure 1.

Example of a 3-DOF spherical robot system.

$D(q)\stackrel{..}{q}+C(q,\stackrel{.}{q})\stackrel{.}{q}+G(q)=u,$

where

$\begin{array}{l}D=\left[\begin{array}{ccc}{I}_z+c_2^2(I_x-I_z+M{r}^2)& 0& 0\\ 0& I_y+M{r}^2& 0\\ 0& 0& M\end{array}\right]\\ C=[\begin{array}{c}{c}_2\left(Mr{c}_2\stackrel{.}{r}-(M{r}^2+I_x-I_z)s_2{\stackrel{.}{q}}_2\right)\\ (M{r}^2+I_x-I_z)c_2{s}_2{\stackrel{.}{q}}_1\\ -Mr{c}_2^2{\stackrel{.}{q}}_1\end{array}\\ \begin{array}{ccc}\hspace{0.17em}& -(M{r}^2+I_x-I_z)c_2{s}_2{\stackrel{.}{q}}_1& Mr{c}_2^2{\stackrel{.}{q}}_1\\ \hspace{0.17em}& Mr\stackrel{.}{r}& Mr{\stackrel{.}{q}}_2\\ \hspace{0.17em}& -Mr{\stackrel{.}{q}}_2& 0\end{array}]\end{array}$

and

$G=\left[\begin{array}{c}-Mgr{s}_1{c}_2\\ -Mgr{c}_1{s}_2\\ Mg{c}_1{c}_2\end{array}\right],$

where c_i = cos(q_i), s_i = sin(q_i), and i = 1,2.

4.2 Example 2: A 3-DOF robot system

Figure 2 presents an example of a 3-DOF robotic pendulum system with the generalized coordinate q = [q₁ q₂ q₃]^T. Table 2 lists the system parameters of this example and the control input vector u = [u₁ u₂ u₃]^T. The proposed method and the traditional Lagrange method are shown to obtain the same equations of motion,

Table 2.

System parameters of a 3-DOF robotic pendulum system.

I	Joint ρi	Link mass mi	Link length ri	Link CoM di	Link inertia Ii
1	+Rx	0	[0 0 0] T	[0 0 0] T	zero inertia
2	+Ry	M 1	[0 0 L1] T	[0 0 L1/2] T	diag{I1x,I1y,I1z}
3	+Rx	M 2	[0 0 L2] T	[0 0 L2/2] T	diag{I2x, I2y,I2z}

Figure 2.

Example of a 3-DOF robotic pendulum system.

$D(q)\stackrel{..}{q}+C(q,\stackrel{.}{q})\stackrel{.}{q}+G(q)=u.$

4.3 Example 3: A 10-DOF biped robot

This example considers a 10-DOF fully actuated 3D biped with feet, as shown in Figure 3. The 3D biped robot consists of 7 links: a torso and two legs and two feet. Each leg has a knee and an equal length for the shin and thigh. The 10 actuators are installed as follows: two at the hip, one at the knee and two at the ankle. Table 3 lists the system parameters of this example.

Table 3.

System parameters of a biped robot.

I	Joint ρi	Link mass mi	Link length ri	Link CoM di	Link inertia Ii
0	base	M 0	[0 0 0] T	[0 0 0] T	zero inertia
1	+Rx	0	[–a 0 L0] T	[0 0 0] T	zero inertia
2	+Ry	M 1	[0 0 0] T	[0 0 L1/2] T	zero inertia
3	+Ry	M 2	[0 0 L1] T	[0 0 L2/2] T	zero inertia
4	+Ry	0	[0 0 L2] T	[0 0 0] T	zero inertia
5	+Rx	M 3	[0 0 0] T	[0 −W L3/2] T	zero inertia
6	−Rx	0	[0 −W 0] T	[0 0 0] T	zero inertia
7	+Ry	M 2	[0 0 0] T	[0 0 −L2/2] T	zero inertia
8	+Ry	M 1	[0 0 −L2] T	[0 0 −L1/2] T	zero inertia
9	+Ry	0	[0 0 −L1] T	[0 0 0] T	zero inertia
10	−Rx	M 0	[0 0 0] T	[0 0 −L0/2] T	zero inertia

Figure 3.

Biped coordinate system when the left foot is the stance foot.

Let n = 10, the generalized coordinate q = [q₁ q₂ ṫ q₁₀] ^T , and the control input vector u = [u₁ ṫ u₁₀]^T. In the single-support phase, the biped's dynamic equation can be obtained as

$D(q)\stackrel{..}{q}+C(q,\stackrel{.}{q})\stackrel{.}{q}+G(q)=u,$

where

$\begin{array}{l}D={\displaystyle \sum _{i=1}^n{m}_i{J}_{vi}^T}{J}_{vi}\\ C={\displaystyle \sum _{i=1}^n{m}_i{J}_{vi}^T}{\stackrel{.}{J}}_{vi}\end{array}$

and

$G={\displaystyle \sum _{i=1}^n{m}_i}g{J}_{vi}^T{e}_3.$

Aside from the equations of motion, significant amounts of important data for the analysis of the biped's walking ability are also computed with minimal additional effort. For example, the biped ground reaction force F = [F_x F_y F_z] ^T and the zero moment point (ZMP) during the single-support phase [x_zmp y_zmp] ^T can be set up as:

$\left[\begin{array}{c}{F}_x\\ F_y\\ F_z\end{array}\right]=\left[\begin{array}{c}{\displaystyle \sum _{i=0}^n{m}_i{\stackrel{..}{x}}_i}\\ {\displaystyle \sum _{i=0}^n{m}_i{\stackrel{..}{y}}_i}\\ {\displaystyle \sum _{i=0}^n{m}_i(g+{\stackrel{..}{z}}_i)}\end{array}\right]$

and

$\left[\begin{array}{c}{x}_{zmp}\\ y_{zmp}\end{array}\right]=\frac{1}{{\displaystyle \sum _{i=0}^n{m}_i(g+{\stackrel{..}{z}}_i)}}\left[\begin{array}{c}{\displaystyle \sum _{i=0}^n{m}_i{x}_i(g+{\stackrel{..}{z}}_i)-{\displaystyle \sum _{i=0}^n{m}_i{z}_i{\stackrel{..}{x}}_i}-{\displaystyle \sum _{i=0}^n{({\tau }_i)}_y}}\\ {\displaystyle \sum _{i=0}^n{m}_i{y}_i(g+{\stackrel{..}{z}}_i)-{\displaystyle \sum _{i=0}^n{m}_i{z}_i{\stackrel{..}{y}}_i}+{\displaystyle \sum _{i=0}^n{({\tau }_i)}_x}}\end{array}\right]$

where g is the gravity acceleration, p_ci = [x_i, y_i, z_i] ^T is the centre of mass of the i-th link, and τ_i = [(τ_i) _x , (τ_i) _y , (τ_i) _z ] ^T is the moment of the i-th link. For authentication, the 10-DOF biped robot online simulation work has been implemented on MATLAB. After symbolic calculation, we transfer the dynamic equations into C-Mex files to improve the efficiency of the biped robot walking simulation.

4.4 Example 4: A simple 2-DOF closed-chain robot

We can also apply the proposed method to calculate the equations of motion for simple closed-chain robotic systems. Figure 4 shows an example of a 2-DOF four-bar linkage robot system, where Link 1 and Link 4 are parallel links and their absolute angle θ₁ is controlled by the input torque u₁. Similarly, Link 2 and Link 3 are parallel links and their absolute angle θ₂ is controlled by the input torque u₂. The four-bar robot system can be treated as 2 serial-chain robot systems. The first serial-chain robot consists of Link 1 and Link 2 with generalized coordinates

Figure 4.

An example of a 2-DOF four-bar linkage robot system.

$q=\left[\begin{array}{c}{q}_1\\ q_2\end{array}\right]=\left[\begin{array}{c}{\theta }_1\\ {\theta }_2-{\theta }_1\end{array}\right]=A_1\left[\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right],$

where

$A_1=\left[\begin{array}{cc}1& 0\\ -1& 1\end{array}\right].$

Similarly, the second serial-chain robot consists of Link 4 and Link 3 with generalized coordinates

$q=\left[\begin{array}{c}{q}_1\\ q_2\end{array}\right]=\left[\begin{array}{c}{\theta }_2\\ {\theta }_1-{\theta }_2\end{array}\right]=A_2\left[\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right],$

where

$A_2=\left[\begin{array}{cc}0& 1\\ 1& -1\end{array}\right].$

For the first robot system, the system matrices D₁, C₁ and G₁ can be calculated using the above method. For the second robot system, the matrices D₂, C₂ and G₂ can be similarly obtained. The equations of motion for the four-bar linkage can thus be determined as:

$\begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left(A_1^T{D}_1{A}_1+A_2^T{D}_2{A}_2\right)\stackrel{..}{\theta }+\\ +\left(A_1^T{C}_1{A}_1+A_2^T{C}_2{A}_2\right)\stackrel{.}{\theta }+\left(A_1^T{G}_1+A_2^T{G}_2\right)=u,\end{array}$

and hence

$\begin{array}{l}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}D=A_1^T{D}_1{A}_1+A_2^T{D}_2{A}_2=\\ \hspace{0.17em}\\ \left[\begin{array}{cc}\frac{1}{2}(m_1+2m_2)L_1^2& \frac{1}{2}(m_1+m_2)L_1{L}_2\cos ({\theta }_1-{\theta }_2)\\ \frac{1}{2}(m_1+m_2)L_1{L}_2\cos ({\theta }_1-{\theta }_2)& \frac{1}{2}(m_2+2m_1)L_2^2\end{array}\right]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}C=A_1^T{C}_1{A}_1+A_2^T{C}_2{A}_2=\\ \hspace{0.17em}\\ \left[\begin{array}{cc}0& \frac{1}{2}(m_1+m_2)L_1{L}_2\sin ({\theta }_1-{\theta }_2){\stackrel{.}{\theta }}_2\\ \frac{-1}{2}(m_1+m_2)L_1{L}_2\sin ({\theta }_1-{\theta }_2){\stackrel{.}{\theta }}_1& 0\end{array}\right]\end{array}$

and

$G=A_1^T{G}_1+A_2^T{G}_2=\left[\begin{array}{c}g(m_1+m_2)L_1\cos {\theta }_1\\ g(m_1+m_2)L_2\cos {\theta }_2\end{array}\right].$

5. Conclusion

This paper has developed a decomposition and cross-product-based method to compute robot dynamics, offering several examples that demonstrate the versatility of the proposed method. The focus of this paper is first to establish a clear relationship between Lagrange equations and Newton-Euler equations, and second to develop an extremely intuitive and systematic method to compute the inertia matrix, the Coriolis and centrifugal matrix, and the gravity force vector of robot dynamics equations.