New development of the dynamic modeling and the inverse dynamic analysis for flexible robot

Abstract

When a segment of a flexible link of a flexible robot is currently sliding through a prismatic joint, it is usually assumed that the elastic deformation of the segment equals to zero. This is a kind of time-dependent boundary condition when formulating the dynamics model of a flexible robot consisting of prismatic joints. Hence, the dynamic modeling and especially the inverse dynamic analysis of the flexible robots with the prismatic joints are challenging. In this article, we present a new development of the dynamic modeling method for a generic two-link flexible robot that consists of a prismatic joint and a revolute joint. Moreover, a new bisection method-based algorithm is proposed to analyze the inverse dynamic responses of the flexible robots. Since the bisection method is a rapid converging method in mathematics, the proposed algorithm is effectively applicable to solving the inverse dynamic problem of a flexible robot in a robust manner. Last, the numerical simulation results show the effectiveness and the robustness of the proposed method.

Keywords

Flexible robot finite-element method inverse dynamics prismatic joint

Introduction

A very important issue in the dynamic modeling and analysis of the flexible robot arm is the computation of joint torques/forces that are required for driving the arm to move along a prescribed trajectory. This is a kind of inverse dynamic problem of manipulators that have been extensively studied for rigid link robots. Different from the inverse dynamics of a rigid link arm, the inverse dynamics of a flexible arm is an order-of-magnitude complex and intricate because of the elastic deformation of the arm links. In particular, the deformation of the arm links, however, is not negligible as the required positioning accuracy of the arm tip increases. Robot arms tracking curved trajectories at a high speed with a large acceleration in a large space, while the tracking error at the arm tip must be kept within an accurate limitation, need to be considered with the elastic deformation of links. If one tries to make such an arm structure sufficiently stiff such that the total elastic deformation is less than the allowed limitation at the arm tip, it could end up with extremely heavy arm mass and inertia. Since drive motors are, in many cases, limited in peak torques, the heavy arm mass and inertia prevent rapid accelerations. A practical solution will not be to design an extremely stiff arm but to allow the arm structure to deform under high accelerations, and compensate for the dynamic deformations by control of dynamic applied torques/forces. Namely, one can estimate the arm tip deformation and compensate for the tip error by computing and modifying the dynamic applied torques/forces acting on joints so that the deflected arm tip may track the trajectory exactly.

Decades ago, the dynamic modeling and analysis of flexible robots have been well documented.^1

–8 There have also been attempts to generalize the kinematic and dynamic modeling problems for the multilink flexible robots, and most of these investigations focus on the multilink flexible robots consisting of all revolute joints.^4,9

–13 There were also studies considering the use of the prismatic joints for the flexible robots.^10,11,14,15 However, these studies mainly focused on the use of the assumed modes method (AMM) for the dynamic modeling. Al-Bedoor and Khulief¹⁴ investigated the dynamics of a sliding flexible link through a prismatic joint. Pan et al.¹⁵ developed the modeling for a two-link flexible robot with a prismatic joint. Korayem et al.¹¹ addressed the dynamic modeling of a multilink flexible robot, of which each link rotates and reciprocates through a prismatic joint. Khadem and Pirmohammadi¹⁰ formulated the dynamic equations for a multilink flexible robot consisting of the prismatic joints as well. Though there have been several research works working on the dynamic modeling of different kinds of specific flexible robots, little attention has been paid to generalize the finite-element formulation of the dynamic equations for a generalized two links flexible robot that consist of total four specific configurations of kinematic chain. The four configurations include the two revolute joints configuration ( $C_{R R}$ ), the two prismatic joints configuration ( $C_{P P}$ ), the revolute–prismatic joints configuration ( $C_{R P}$ ), and the prismatic–revolute joints configuration ( $C_{P R}$ ).

Notice that the flexible robots are the continuous dynamical systems, which are often characterized by an infinite number of degrees of freedom and are governed by the nonlinear coupled differential equations, and the exact solution of such systems is not practical. To formulate the dynamic equations, the continuous dynamical systems are usually discretized using two main methods, including (i) the AMM and (ii) finite-element method (FEM). Theodore and Ghosal¹⁶ and Dwivedy and Eberhard⁴ showed that the main drawback of AMM is the difficulty in finding the modes for links with the nonregular cross sections and the multilink robots, and FEM normally requires fewer computations; therefore, FEM lends itself ideally suited for the modeling of multilink flexible robots. The advantages of FEM over AMM have also motivated the finite-element formulation of dynamic equations for the generalized two links flexible robot that is presented in this article.

It can be seen that, in recent years, though there have been several studies related to the flexible manipulators, a few attempts focus on the inverse dynamic analysis for the flexible robots, especially the flexible robots with the translational joints. Most of the works only take into account the revolute joints^{17

–28} when considering the dynamic modeling and analysis of a flexible robot. The large deformation of a flexible link was studied by Celentano.¹⁷ The dynamic modeling methods were developed by Li et al.²¹ and Wei et al.²² The vibration of the flexible robots was focused in a number of investigations.^18
–20,23 Also, the control issues were addressed in the recent works.^24

–28

In the literature, the inverse dynamic analysis for flexible robots has been studied by a number of authors. Kwon and Book²⁹ used AMM for the modeling of the linear inverse dynamics of a single flexible link robot with one rotational joint. Asada et al.³⁰ employed the coordinates of virtual rigid links to simplify the formulation of the dynamic equations for a flexible robot consisting of two revolute joints. A numerical algorithm was constructed to compute the actuator torques. Bayo et al.³¹ and Trautt and Bayo³² studied the response function of the applied torques imposing on two revolute joints of a flexible robot. The inverse dynamic problem was solved in frequency domain. Boyer and Khalil³³ investigated an AMM formulation incorporated with Newton–Euler approach to derive the dynamic equations. Based on the formulated model, an iterative technique was studied to find out the root of the inverse dynamic equations. Carrera and Serna³⁴ addressed the inverse dynamics of a flexible robot with all rotary joints. Ledesma and Bayo³⁵ studied a nonrecursive algorithm for solving the inverse dynamic problem of a flexible robot with revolute joints also. Zhaocai and Yueqing³⁶ investigated the inverse dynamics of a flexible parallel robot.

It is clearly seen that the inverse dynamics of flexible robots has been investigated; however, most of the investigations use AMM and focus on the robot architecture with all revolute joints. The inverse dynamics of flexible robot with the prismatic joints has been overlooked. Therefore, in this article, we develop the dynamic modeling and analysis of the general flexible manipulator by formulating the dynamic equations with the use of FEM and solving the inverse dynamic problem effectively. A numerical algorithm is proposed to calculate the time-varying values of all joint variables and actuator torques/forces. To construct the algorithm, the effective and robust bisection method is used to find out the feasible solution of the inverse dynamics in time domain. Since no inverse transformation or approximation of the input/output signals is needed for the inversion, the algorithm can be implemented in a simple and effective manner. The numerical simulation results of the inverse dynamics show the effectiveness and the robustness of the proposed method.

In the following, the article is organized with four more sections. The dynamic modeling of the system is detailed in the second section, where the generalized mechanism of the robot is descripted and interpreted. Based on the dynamic model formulated, the inverse dynamic algorithm is proposed in the third section, in which the use of the bisection method is explained. To verify the research results, the numerical simulation is addressed in the fourth section. Conclusions are made in the last section.

Equations of motion

The schematic model of the generalized two-link flexible robot is shown in Figure 1. L ₁ and L ₂ are the lengths of the links, respectively. m ₁ is the mass of link 1, and m ₂ is the mass per length unit of link 2. m_t is the mass of a payload at the end point of link 2. The joints J ₁ and J ₂ are either a prismatic joint $(P)$ or a revolute $(R)$ joint. With all possible combinations of the different joint types, there exist four different configurations of a two-link flexible robot. The configurations as denoted as $C_{P R}$ , $C_{R P}$ , $C_{R R}$ , and $C_{P P}$ , which are shown in Figure 2. For example, in the configuration $C_{P R}$ , the first joint is a P joint, and the second one is a R joint.

Figure 1.

Generalized schematic model of a two-link flexible robot.

Figure 2.

Specific models of a two-link flexible robot. (a) C_RR model, (b) C_RP model, (c) C_PR model, and (d) C_PP model.

In Figure 1, $X O Y$ is a reference frame, $X_{1} O_{1} Y_{1}$ is a link frame attached to the end point of the link 1, and $X_{2} O_{2} Y_{2}$ is another link attached to link 2, where $O_{2} \equiv O_{1}$ .

$q_{i} (t)$ is the generalized coordinate corresponding to the joint J_i . $d_{i} (t)$ and $θ_{i} (t)$ are the translational and rotational joint variables denoted (see Table 1). $ζ_{1} (t)$ and $ζ_{2} (t)$ are the external generalized forces applying on J ₁ and J ₂ correspondingly. For the configurations, $F_{i} (t)$ and $τ_{i} (t)$ are the applied force/torque imposing on the translational and rotational joint. In Table 1, for all the configurations, the joint displacements and applied force/torque are summarized.

Table 1.

Joint displacements and applied force/torque.

	$C_{P R}$	$C_{R P}$	$C_{R R}$	$C_{P P}$
$q_{1}, q_{2}$	$d_{1}, θ_{2}$	$θ_{1}, d_{2}$	$θ_{1}, θ_{2}$	$d_{1}, d_{2}$
$ζ_{1}, ζ_{2}$	$F_{1}, τ_{2}$	$τ_{1}, F_{2}$	$τ_{1}, τ_{2}$	$F_{1}, F_{2}$

Notice that for $C_{P R}$ and $C_{P P}$ , $O X$ and $O_{1} X_{1}$ could not coincide, there could exist an angle $α_{0}$ between the two axes. Also, for $C_{R P}$ and $C_{P P}$ , there could exist an angle $β_{0}$ between $O_{1} X_{1}$ and $O_{2} X_{2}$ .

As usual, to determine the kinematic relationships for the configurations of the robot, the transformation matrices describing the kinematics of the motion of the joints must be derived. For all the configurations, the homogeneous transformation matrices $H_{01}$ and $H_{12}$ are formulated and presented in Appendix. Notice that $H_{01}$ describes the motion of the link 1 relative to the base and $H_{12}$ describes the relative motion of the link 2.

Without loss the generality, it is assumed that the first link is rigid. The second link is flexible, which is divided into n beam elements. The elements are assumed interconnected at nodes. Each element j, $j = 1 \div n$ , has four elastic displacements at the two nodes that are the flexural displacements $(u_{2 j - 1}, u_{2 j + 1})$ and the slope displacements $(u_{2 j}, u_{2 j + 2})$ .

For the generalized model, the dynamic equation of motion relies on the Lagrange’s equations with Lagrange function $L = T - P$ that is given by

\frac{d}{d t} {(\frac{\partial L}{\partial \dot{Q} (t)})}^{T} - {(\frac{\partial L}{\partial Q (t)})}^{T} = F (t)

where T and P are the total kinetic energy and the total potential energy of the system. The vector $F = {[\begin{matrix} ζ_{1} & ζ_{2} & 0 & . & . & 0 & 0 \end{matrix}]}^{T}$ is the external generalized forces acting along with components of the vector of the generalized coordinates $Q = {[\begin{matrix} q_{1} & q_{2} & u_{1} & u_{2} & ... & u_{2 n + 1} & u_{2 n + 2} \end{matrix}]}^{T}$ . The size of F and Q is $(2 n + 4) \times 1$ .

The equation of motion can be expressed as

M (Q) \ddot{Q} + C (Q, \dot{Q}) \dot{Q} + D \dot{Q} + K Q = F (t)

where the structural damping matrix D and the Coriolis and centrifugal matrix C are correspondingly calculated as

C (Q, \dot{Q}) \dot{Q} = \dot{M} (Q) \dot{Q} - \frac{1}{2} (\frac{\partial}{\partial Q} {({\dot{Q}}^{T} M (Q) \dot{Q})}^{T})

D = α M + β K

The notations α and β are the damping ratios of the system (Zhang et al.).³⁷ The size of the global mass matrix M and the global stiffness matrix K is $(2 n + 4) \times (2 n + 4)$ .

Consider an element j of the link 2, a point $r_{2 j} (x)$ represented in $X_{2} O_{2} Y_{2}$ , where $0 \leq x \leq l_{e}$ and $l_{e} = \frac{L_{2}}{n}$ , can be computed for the configurations as follows

r_{2 j} (x) = [\begin{matrix} {[\begin{matrix} (j - 1) l_{e} + x & w_{j} (x) \end{matrix}]}^{T}, for C_{P R} and C_{R R} \\ {[\begin{matrix} (j - 1) l_{e} - L_{2} + d_{2} + x & w_{j} (x) \end{matrix}]}^{T}, for C_{R P} and C_{P P} \end{matrix}

Notice that $w_{j} (x)$ is the flexural displacement at the point x on the element j. The value of $w_{j} (x)$ is determined via the vector of the shape functions $N_{j} (x) = {[\begin{matrix} ϕ_{1} (x) & ϕ_{2} (x) & ϕ_{3} (x) & ϕ_{4} (x) \end{matrix}]}^{T}$ presented in the work by Zhang et al.³⁷

w_{j} (x) = N_{j} (x) Q_{j}

where $Q_{j} = {[\begin{matrix} u_{2 j - 1} & u_{2 j} & u_{2 j + 1} & u_{2 j + 2} \end{matrix}]}^{T}$ is the vector of elastic displacements of the element j.

Therefore, the point $r_{2 j}$ represented in $X O Y$ can be expressed as

r_{02 j}^{​} = H_{01} H_{12} r_{2 j}^{​}

Notice that the end point of the link 2, $χ = {[\begin{matrix} x_{E} & y_{E} \end{matrix}]}^{T}$ , can be computed with equation (7), given that $j = n$ and $x = l_{e}$ .

The kinetic energy of the element j of the link 2 is determined as

T_{2 j}^{​} = \frac{1}{2} \int_{0}^{l_{e}} m_{2} {\dot{r}}_{02 j}^{T} {\dot{r}}_{02 j} d x = \frac{1}{2} {\dot{Q}}_{j g}^{T} M_{j}^{​} {\dot{Q}}_{j g}^{​}

where $Q_{j g}^{} = {[\begin{matrix} q_{1} & q_{2} & u_{2 j - 1} & u_{2 j} & u_{2 j + 1} & u_{2 j + 2} \end{matrix}]}^{T}$ . The matrix $M_{j}$ can be computed as

M_{j}^{​} (s, e) = \int_{0}^{l_{e}} m_{2} {[\frac{\partial r_{02 j}}{\partial Q_{j s}}]}^{T} [\frac{\partial r_{02 j}}{\partial Q_{j e}}] d x; s, e = 1 \div 6

Notice that $Q_{j s}$ and $Q_{j e}$ are the s’th and e’th element of $Q_{j g}^{}$ , respectively. It can be shown that $M_{j}^{}$ is of the form

M_{j}^{​} = [\begin{matrix} m_{11} & m_{12} & m_{13} & m_{14} & m_{15} & m_{16} \\ m_{21} & m_{22} & m_{23} & m_{24} & m_{25} & m_{26} \\ m_{31} & m_{32} & ​ & ​ & ​ & ​ \\ m_{41} & m_{42} & ​ & M_{base} & ​ & ​ \\ m_{51} & m_{52} & ​ & ​ & ​ & ​ \\ m_{61} & m_{62} & ​ & ​ & ​ & ​ \end{matrix}]

where

M_{base}^{​} = [\begin{matrix} \frac{13}{35} m_{2} l_{e} & \frac{11}{210} m_{2} l_{e}^{2} & \frac{9}{70} m_{2} l_{e} & - \frac{13}{420} m_{2} l_{e}^{2} \\ \frac{11}{210} m_{2} l_{e}^{2} & \frac{1}{105} m_{2} l_{e}^{3} & \frac{13}{420} m_{2} l_{e}^{2} & - \frac{1}{140} m_{2} l_{e}^{3} \\ \frac{9}{70} m_{2} l_{e} & \frac{13}{420} m_{2} l_{e}^{2} & \frac{13}{35} m_{2} l_{e} & - \frac{11}{210} m_{2} l_{e}^{2} \\ - \frac{13}{420} m_{2} l_{e}^{2} & - \frac{1}{140} m_{2} l_{e}^{3} & - \frac{11}{210} m_{2} l_{e}^{2} & \frac{1}{105} m_{2} l_{e}^{3} \end{matrix}]

For $C_{P R}$ model, the other entries of the elemental matrix $M_{j}^{}$ are calculated as follows

\begin{array}{l} m_{11} = m_{2} l_{e}; \\ m_{12} = - \frac{1}{12} m_{2} l_{e} [(6 u_{2 j - 1} + 6 u_{2 j + 1} + l_{e} u_{2 j} - l_{e} u_{2 j + 2}) sin θ_{2} + 6 l_{e} (1 - 2 j) cos θ_{2}]; \\ m_{13} = m_{15} = 0.5 m_{2} l_{e} cos θ_{2}; m_{14} = m_{16} = \frac{1}{2} m_{2} l_{e}^{2} cos θ_{2}; \\ m_{21} = m_{12}; m_{23} = \frac{1}{20} m_{2} l_{e}^{2} (10 j - 7); m_{24} = \frac{1}{60} m_{2} l_{e}^{3} (5 j - 3); \\ m_{25} = \frac{1}{20} m_{2} l_{e}^{2} (10 j - 3); m_{26} = - \frac{1}{60} m_{2} l_{e}^{3} (5 j - 2); \\ m_{22} = \frac{1}{210} m_{2} l_{e} (\begin{array}{l} 210 l_{e}^{2} j (j - 1) + l_{e}^{2} (2 u_{2 j}^{2} - 3 u_{2 j} u_{2 j + 2} + 2 u_{2 j + 2}^{2}) + \\ 22 l_{e} (u_{2 j - 1} u_{2 j} - u_{2 j + 1} u_{2 j + 2}) + + 13 l_{e} (u_{2 j} u_{2 j + 1} - u_{2 j - 1} u_{2 j + 2}) + \\ 78 (u_{2 j - 1}^{2} + u_{2 j + 1}^{2}) + 70 l_{e}^{2} + 54 u_{2 j - 1} u_{2 j + 1} \end{array}); \\ m_{31} = m_{13}; m_{32} = m_{23}; m_{41} = m_{14}; m_{42} = m_{24}; m_{51} = m_{15}; m_{52} = m_{25}; m_{61} = m_{16}; m_{62} = m_{26} \end{array}

For $C_{R P}$ model, the other entries of the elemental matrix $M_{j}^{}$ are calculated as follows

\begin{array}{l} m_{11} = f (m_{2}, l_{e}, j, L_{1}, L_{2}, d_{2}, Q_{j}); \\ m_{12} = m_{21} = - \frac{1}{12} m_{2} l_{e} (u_{2 j} l_{e} - u_{2 j + 2} l_{e} + 6 u_{2 j - 1} l_{e} + 6 u_{2 j + 1} l_{e}); \\ m_{13} = m_{31} = - \frac{1}{20} m_{2} l_{e} (10 L_{1} - 10 L_{2} + 10 d_{2} - 7 l_{e} + 10 j l_{e}); \\ m_{14} = m_{41} = - \frac{1}{60} m_{2} l_{e} ​^{2} (5 L_{1} - 5 L_{2} + 5 d_{2} - 3 l_{e} + 5 j l_{e}); \\ m_{15} = m_{51} = - \frac{1}{20} m_{2} l_{e} (10 L_{1} - 10 L_{2} + 10 d_{2} - 3 l_{e} + 10 j l_{e}); \\ m_{16} = m_{61} = - \frac{1}{60} m_{2} l_{e} ​^{2} (5 L_{1} - 5 L_{2} + 5 d_{2} - 2 l_{e} + 5 j l_{e}); \\ m_{22} = m_{2} l_{e}; m_{23} = m_{32} = m_{24} = m_{42} = 0; \\ = m_{25} = m_{52} = m_{26} = m_{62} = 0 \end{array}

For $C_{R R}$ model, the other entries of the elemental matrix $M_{j}^{}$ are calculated as follows

\begin{array}{l} m_{11} = f (m_{2}, l_{e}, L_{1}, j, Q_{j}, θ_{2}) \neq m_{12} = f (m_{2}, l_{e}, L_{1}, j, Q_{j}, θ_{2}) \\ m_{13} = m_{31} = \frac{1}{20} m_{2} l_{e} (10 L_{1} cos θ_{2} + 10 j l_{e} - 7 l_{e}); \\ m_{14} = m_{41} = \frac{1}{60} m_{2} l_{e}^{2} (5 L_{1} cos θ_{2} + 5 j l_{e} - 3 l_{e}) \\ m_{15} = m_{51} = \frac{1}{20} m_{2} l_{e} (10 L_{1} cos θ_{2} + 10 j l_{e} - 3 l_{e}); \\ m_{16} = m_{61} = - \frac{1}{60} m_{2} l_{e}^{2} (5 L_{1} cos θ_{2} + 5 j l_{e} - 2 l_{e}) \\ m_{21} = m_{12}; m_{22} = f (m_{2}, l_{e}, j, Q_{j}); m_{23} = m_{32} = \frac{1}{20} m_{2} l_{e}^{2} (10 j - 7); \\ m_{24} = m_{42} = \frac{1}{60} m_{2} l_{e}^{3} (5 j - 3) \\ m_{25} = m_{52} = \frac{1}{20} m_{2} l_{e}^{2} (10 j - 3); m_{26} = m_{62} = - \frac{1}{60} m_{2} l_{e}^{3} (5 j - 2) \end{array}

For $C_{P P}$ model, the other entries of the elemental matrix $M_{j}^{}$ are calculated as follows

\begin{array}{l} m_{11} = m_{2} l_{e}; m_{12} = m_{21} = 0; m_{13} = m_{31} = \frac{1}{2} m_{2} l_{e}; m_{14} = m_{41} = \frac{1}{12} m_{2} l_{e}^{2}; \\ m_{15} = m_{51} = \frac{1}{2} m_{2} l_{e}; m_{16} = m_{61} = - \frac{1}{12} m_{2} l_{e}^{2}; \\ m_{22} = m_{2} l_{e}; m_{23} = m_{32} = m_{24} = m_{42} = m_{25} = m_{52} = m_{26} = m_{62} = 0 \end{array}

The kinetic energy of the link 2 is yielded as

T_{e}^{​} = \sum_{j = 1}^{n} T_{2 j}^{​} = \frac{1}{2} {\dot{Q}}^{T} M_{e}^{​} {\dot{Q}}^{​}

The matrix $M_{e}^{}$ is assembled from all the elemental matrices $M_{j}^{}$ . The pseudocodes of the algorithm for the construction of $M_{e}^{}$ are presented in Appendix.

The total kinetic energy of the system is determined as

T^{​} = T_{1}^{​} + T_{e}^{​} + T_{P}^{​} = \frac{1}{2} {\dot{Q}}^{T} M \dot{Q}

where $T_{1}^{}$ and $T_{P}^{}$ are the kinetic energy of the first link and the payload that can be easily determined via the rigid model.

Let us denote E and I as Young’s modulus and the inertial moment of the link 2. The elastic potential energy of the element j, $P_{j}^{}$ , is computed as

where the elemental stiffness matrix $K_{j}^{}$ is presented as

P_{j}^{​} = \frac{1}{2} \int_{0}^{l_{e}} E I {(\frac{\partial^{2} w_{j} (x, t)}{\partial x_{​}^{2}})}^{2} d x = \frac{1}{2} Q_{j}^{T} K_{j}^{​} Q_{j}

K_{j}^{​} = [\begin{matrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{12 E I}{l_{e}^{3}} & \frac{6 E I}{l_{e}^{2}} & - \frac{12 E I}{l_{e}^{3}} & \frac{6 E I}{l_{e}^{2}} \\ 0 & 0 & \frac{6 E I}{l_{e}^{2}} & \frac{4 E I}{l_{e}} & - \frac{6 E I}{l_{e}^{2}} & \frac{2 E I}{l_{e}} \\ 0 & 0 & - \frac{12 E I}{l_{e}^{3}} & - \frac{6 E I}{l_{e}^{2}} & \frac{12 E I}{l_{e}^{3}} & - \frac{6 E I}{l_{e}^{2}} \\ 0 & 0 & \frac{6 E I}{l_{e}^{2}} & \frac{2 E I}{l_{e}} & - \frac{6 E I}{l_{e}^{2}} & \frac{4 E I}{l_{e}} \end{matrix}]

The total elastic potential energy of the whole system is yielded as

P = \sum_{j = 1}^{n} P_{j} = \frac{1}{2} Q^{T} K Q

The global stiffness matrix K is constructed from all matrices $K_{j}^{}$ .

Substituting M and K into equations (3) and (4) yields C and D accordingly. Finally, substituting all components M, C, D, and K into equation (2) obtains the dynamic equations for the generalized model.

Boundary conditions

For $C_{P R}$ and $C_{R R}$ configurations, the link 2 connects with the link 1 via a rotational joint. It is assumed that there are no elastic displacements at the first node of the element 1 of the link 2 ( $u_{1} = u_{2} = 0$ ). Therefore, the third row, the fourth row, the third column, and the fourth column of the matrices $M, K, D$ , and C are eliminated. The third and the fourth components of $F (t)$ and $Q (t)$ are eliminated as well. For $C_{R P}$ and $C_{P P}$ configurations, when an element k of the second link is sliding inside the translational joint, the displacements of the element are equal to zero, because the prismatic joint hub is rigid.

Inverse dynamics analysis

The inverse analysis of the flexible robots is a complex problem since there exists $2 n$ nodal coordinates in the vector of the generalized coordinates, $Q = {[\begin{matrix} q_{1} & q_{2} & u_{1} & u_{2} & ... & u_{2 n + 1} & u_{2 n + 2} \end{matrix}]}^{T}$ . Given a trajectory of the end-effector, $χ (t) = {[\begin{matrix} x_{E} (t) & y_{E} (t) \end{matrix}]}^{T}$ , represented in the workspace, the joint trajectory $Q (t)$ and the force/torque $F (t)$ could not be determined with the use of the existed techniques that have been using for the inverse kinematics and inverse dynamics of the rigid robots. Therefore, in this section, an effective numerical algorithm is constructed to find a feasible and robust solution, $F (t)$ , to the inverse dynamic problem. The idea of the algorithm construction is that the inverse dynamic solution for a flexible robot could be obtained by modifying the inverse dynamic solution for a corresponding rigid robot that has the same system parameters. The modification is calculated iteratively in such a way that the error between the desired trajectory $χ (t)$ and the computed trajectory $χ^{C} (t)$ minimizes.

The algorithm is proposed, which is based on the bisection method to compute $F (t)$ approximately, given the trajectory $χ (t)$ and the allowable error ε of the trajectory. For a given interval of time $[0, T]$ , a time stepwise is selected as $Δ t = T / m$ , where m is an integer number.

At a time $t_{i} \in [0, T]$ , for i from 0 to m, the computational steps of the algorithm are as follows.

Step 1. Compute the inverse dynamics of the robot without any elastic deformation, based on the inputs $χ (t_{i})$ . In other words, $F^{r} (t_{i})$ is computed with equation (2), in which the effect of the elastic deformation is neglected.

Step 2. Compute the forward dynamics of the flexible model. Based on $F^{r} (t_{i})$ computed in step 1, integrating equation (2) yields $Q (t_{i})$ ; and the end-effector trajectory $χ^{C} (t_{i})$ is computed with equation (7).

Step 3. Compute the trajectory error δ.

The trajectory error δ is approximated as the distance from the point $χ^{C} (t_{i})$ to the velocity vector $\dot{χ} (t_{i})$ (tangent vector) of the curve $χ (t_{i})$ at the time t_i . Therefore, it can be calculated as

δ = \frac{|Δ χ \times \dot{χ} (t_{i})|}{|\dot{χ} (t_{i})|}, where Δ χ = χ (t_{i}) - χ^{C} (t_{i})

Step 4. Estimate the interval $F^{U}$ , in between which the feasible value of $F (t)$ can be found. Corresponding to $Δ χ$ , compute $Δ Θ = {[\begin{matrix} Δ q_{1} & Δ q_{2} \end{matrix}]}_{} ^{T}$ and its time derivatives $Δ \dot{Θ}$ and $Δ \ddot{Θ}$ via the inverse kinematics of the rigid model. Substituting $Δ Θ$ , $Δ \dot{Θ}$ , and $Δ \ddot{Θ}$ into equation (2) (neglecting the elastic deformation) yields $Δ F$ . Estimate $F^{U}$ as $F^{U} = F^{r} (t_{i}) + κ Δ F$ , where the weight κ is a positive integer. To ensure that $F^{U}$ is a large enough interval, the selected value of κ should be greater than 2.

Step 5. Modify $F^{r} (t_{i})$ to refine $F (t_{i})$ .

Loop the following substeps:

– Compute $F_{M} = \frac{1}{2} [F^{U} + F^{r} (t_{i})]$ ,

– Implement step 2 to recompute $χ^{C}$ with the input $F_{M}$ ,

– Set $δ_{f} : = δ$

– Implement step 3 to recompute δ,

– If $δ_{f} < δ$ then $F^{U} : = F_{M}$ else $F^{r} (t_{i}) : = F_{M}$

Until $δ < ε$ .

Step 6. Finally, $F (t_{i}) = F_{M}$ .

All the values $F (t_{1}), F (t_{2}), \dots, F (t_{m})$ are computed while iincreasing from 1 to m, where m is the number of the time steps. The time next to t_i can be calculated as $t_{i + 1} = t_{i} + Δ t$ , where $Δ t = \frac{T}{m}$ . Figure 3 describes the algorithm implementation steps.

Figure 3.

The bisection algorithm diagram.

Numerical simulations

In this section, a numerical simulation of the inverse dynamics for the two configurations of the robot $C_{P R}$ and $C_{R P}$ is presented. The structural parameters of the robots, the desired trajectories, and the other simulation parameters are given in Table 2.

Table 2.

Simulation parameters.

Parameters	$C_{P R}$	$C_{R P}$
Length of link 1 (m)	$L_{1} = 0.2$	$L_{1} = 0.2$
Mass of link 1 (kg)	$m_{1} = 0.157$	$m_{1} = 0.226$
Fixed angle between links	$α_{0} = 0$	$β_{0} = 0$
Length of link 2 (m)	$L_{2} = 0.8$	$L_{2} = 0.8$
Number of elements of the link 2	$n = 5$	$n = 5$
Mass per meter of link 2 (kg/m)	$m_{2} = 0.0785$	$m_{2} = 0.0785$
Young’s modulus (N/m²)	$E = 2 \times 10^{10}$	$E = 2 \times 10^{10}$
Inertial moment of cross section (m⁴)	$I = 8.33 \times 10^{- 13}$	$I = 8.33 \times 10^{- 13}$
Damping ratios	$α = 0.005; β = 0.007$	$α = 0.005; β = 0.007$
Mass of payload (g)	$m_{t} = 20$	$m_{t} = 20$
Desired trajectory and simulation time
Desired trajectory	$χ (t) = [\begin{matrix} 0.5 - 0.3 * sin (t - \frac{π}{2}) \\ 0.3 * sin (t) \end{matrix}]$	$χ (t) = [\begin{matrix} 0.4 + 0.2 * sin (\frac{π}{3} t - \frac{π}{2}) \\ 0.2 * sin (\frac{π}{3} t) \end{matrix}]$
Simulation time (s)	$T_{d} = 4$	$T_{d} = 6$

The results of the simulation for $C_{P R}$ are shown in Figures 4 –6. The simulation results of the corresponding rigid robot of the same configuration are depicted in the figures as well. Figure 4 shows the translational joint displacement $d_{1} (t)$ and the rotation joint displacement $θ_{2} (t)$ for both the rigid robot and the flexible robot.

Figure 4.

Displacement of joints of $C_{P R}$ configuration. (a) The first joint displacement and (b) the second joint displacement.

Figure 5.

The applied force/torque of $C_{P R}$ configuration. (a) The applied force imposing on the first joint and (b) the applied torque imposing on the second joint.

Figure 6.

(a) The flexural displacement at the end-effector and (b) slope displacement at the end-effector of $C_{P R}$ configuration.

Figure 5 shows the resulting curves $F_{1} (t)$ and $τ_{2} (t)$ . $F_{1}^{r} (t)$ and $τ_{2}^{r} (t)$ are calculated with the rigid robot.

The elastic displacements at the end-effector are shown in Figure 6. The maximum value of the flexural displacement is $0.245 (m)$ and the maximum value of the slope displacement is $0.45 (rad)$ .

The simulation results for the configuration $C_{R P}$ are shown in Figures 7 –9.

Figure 7.

The displacement of joints of $C_{R P}$ configuration. (a) The first joint displacement and (b) the second joint displacement.

Figure 8.

The applied torque/ force of $C_{R P}$ configuration. (a) The applied force/torque imposing on the joints and (b) the deviation of the torque/force in comparison with the rigid robot.

Figure 9.

(a) The flexural displacement at the end-effector and (b) slope displacement at the end-effector of $C_{R P}$ configuration.

Figure 7 shows the value of the joint 1 and joint 2 variables with the difference between the two models does not exceed $0.2 (rad)$ . Figure 8 shows the results of the inverse dynamics problem based on the proposed algorithm. These results showed that the deviation value of force/torque is very small on both driving joints. It is easy to realize that the difference between the two rigid models and the flexible model C_RP is very slight due to the flexible link 2 sliding along the translational joints. This is also shown in Figure 9. The elastic displacement value at the last operating point is extremely tiny.

Conclusions

A new development of the dynamic modeling and the inverse dynamic analysis for all four configurations of a generalized two-link flexible robot that consists of a prismatic joint and a revolute joint was presented in this study. A recursive finite-element Lagrangian formulation of the dynamic equations for all the configurations of the robot has been addressed. Finally, the numerical simulation results demonstrated the effectiveness and the robustness of the proposed method. In comparison with the previous works, the proposed investigation has two main advantages. The first one is that the dynamic modeling method newly proposed in this article is addressed for a generic two-link flexible robot, which is constructed with different joint types and different joint orders. Hence, the proposed method is generalized and applicable to all four types of two-link flexible robots. The second advantage is that the inverse dynamics of the generic two-link flexible robot is effectively analyzed using a simplified and robust algorithm.

Note that the accuracy and the speed of the bisection algorithm mainly depend on the chosen allowable error ε between the computed trajectory and the desired trajectory. When choosing a small ε, the accuracy of the inverse dynamic responses is increased, however, the computation speed is decreased.

Footnotes

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) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported by Vingroup Innovation Foundation (VINIF) in project code VINIF.2019.DA08.

ORCID iD

Chu A My

Appendix

The homogeneous transformation matrices

References

Alandoli

Sulaiman

Rashid

, et al. A review study on flexible link manipulators. J Telecomm Elect Comput Eng 2016; 8(2): 93–97.

Ata

Johar

. Dynamic simulation of task constrained of a rigid-flexible manipulator. Int J Adv Robot Syst 2008; 1(2): 61–66.

Benosman

Le Vey

. Control of flexible manipulators: a survey. Robotica 2004; 22(5): 533–45.

Dwivedy

Eberhard

. Dynamic analysis of flexible manipulators, a literature review. Mech Mach Theory 2006; 41(7): 749–777.

Etxebarria

Sanz

Lizarraga

. Control of a lightweight flexible robotic arm using sliding modes. Int J Adv Robot Syst 2005; 2(2): 11.

Halalchi

Laroche

Bara

. Flexible-link robot control using a linear parameter varying systems methodology. Int J Adv Robot Syst 2014; 11(3): 46.

Rahimi

Nazemizadeh

. Dynamic analysis and intelligent control techniques for flexible manipulators: a review. Adv Robot 2014; 28(2): 63–76.

Wang

Yin

Yang

, et al. Motion planning of rigid chain for rigid–flexible coupled robot. Int J Adv Robot Syst 2018; 15(3): 1–10. DOI: 1729881418772815.

Chen

. Dynamic modeling of multi-link flexible robotic manipulators. Comput Struct 2001; 79(2): 183–195.

10.

Khadem

Pirmohammadi

. Analytical development of dynamic equations of motion for a three-dimensional flexible link manipulator with revolute and prismatic joints. IEEE Trans Syst Man Cybernet Part B (Cybernetics) 2003; 33(2): 237–249.

11.

Korayem

Shafei

Dehkordi

. Systematic modeling of a chain of N-flexible link manipulators connected by revolute–prismatic joints using recursive gibbs-appell formulation. Arch Appl Mech 2014; 84(2): 187–206.

12.

Usoro

Nadira

Mahil

. A finite element/Lagrange approach to modeling lightweight flexible manipulators. J Dyn Syst Meas Control 1986; 108(3): 198–205.

13.

Zhang

. Recursive Lagrangian dynamic modeling and simulation of multi-link spatial flexible manipulator arms. Appl Math Mech 2009; 30(10): 1283–1294.

14.

Al-Bedoor

Khulief

. General planar dynamics of a sliding flexible link. J Sound Vib 1997; 206(5): 641–661.

15.

Pan

Scott

Ulsoy

. Dynamic modeling and simulation of flexible robots with prismatic joints. J Mech Des 1990; 112(3): 307–314.

16.

Theodore

Ghosal

. Comparison of the assumed modes and finite element models for flexible multilink manipulators. Int J Robot Res 1995; 14(2): 91–111.

17.

Celentano

. Modeling of flexible robots with varying cross section and large link deformations. J Dyn Syst Meas Control 2016; 138: 1–12.

18.

Karagulle

Malgaca

Dirilmis

, et al. Vibration control of a two-link flexible manipulator. J Vib Control 2015; 23(12): 2023–2034.

19.

Wang

, et al. Two-time scale virtual sensor design for vibration observation of a translational flexible-link manipulator based on singular perturbation and differential games. Sensor 2016; 16: 1–14.

20.

Bian

Gao

, et al. Theoretical and experimental study on vibration control of flexible manipulator based on internal resonance. J Vib Control 2017; 24(15): 3321–3337.

21.

Guo

, Preferable workspace for fatigue life improvement of flexible-joint robots under percussive riveting. J Dyn Syst Meas Control 2017; 139: 1–7.

22.

Wei

Cao

Liu

, et al. Global mode method for dynamic modeling of a flexible-link flexible-joint manipulator with tip mass. Appl Math Model 2017; 48: 787–805.

23.

Zhang

. The impulse spectrum method for vibration suppression of a flexible multi-link robot arm. J Vib Control 2017; 1–17.

24.

Zhang

Liu

. Adaptive boundary control for a flexible manipulator with state constraints using a barrier Lyapunov function. J Dyn Syst Meas Control 2018; 140: 1–7.

25.

Kim

Boo

. Use of a simple mechanical analogy to analytically tune the PD controller of a flexible manipulator system. Shock Vib 2018; 2018: 15.

26.

Lochan

Dey

Roy

, et al. Tracking control with vibration suppression of a two-link flexible manipulator using singular perturbation with composite control design. Adv Intel Syst Comput 2018; 634: 365–377.

27.

Ravandi

Khanmirza

Daneshjou

. Hybrid force/position control of robotic arms manipulating in an uncertain environment based on the adaptive fuzzy sliding mode control. Appl Soft Comput 2018; 70: 864–874.

28.

Mishra

Singh

. Hybrid vibration control of a two-link flexible manipulator. Appl Sci 2019; 715: 1–25.

29.

Kwon

Book

. A time-domain inverse dynamic tracking control of a single-link flexible manipulator. J Dyn Syst Meas Control 1994; 116(2): 193–200.

30.

Asada

Tokumaru

. Inverse dynamics of flexible robot arms: modeling and computation for trajectory control. J Dyn Syst Meas Control 1990; 112(2): 177–185.

31.

Bayo

Papadopoulos

Stubbe

, et al. Inverse dynamics and kinematics of multi-link elastic robots: an iterative frequency domain approach. Int J Robot Res 1989; 8(6): 49–62.

32.

Trautt

Bayo

. Inverse dynamics of flexible manipulators with Coulomb friction or backlash and non-zero initial conditions. Dyn Control 1999; 9(2): 173–195.

33.

Boyer

Khalil

. An efficient calculation of flexible manipulator inverse dynamics. Int J Robot Res 1998; 17(3): 282–293.

34.

Carrera

Serna

. Inverse dynamics of flexible robots. Math Comput Simul 1996; 41(5–6): 485–508.

35.

Ledesma

Bayo

. A non-recursive Lagrangian solution of the non-causal inverse dynamics of flexible multibody systems: the planar case. Int J Numer Methods Eng 1993; 36(16): 2725–2741.

36.

Zhaocai

Yueqing

. Dynamic modeling and inverse dynamic analysis of flexible parallel robots. Int J Adv Robot Syst 2008; 5(1): 13.

37.

Zhang

Liu

. Vibration control for a nonlinear three dimensional flexible manipulator trajectory tracking. Int J Control 2016; 89(8): 1641–1663.