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

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.


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][2][3][4][5][6][7][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][10][11][12][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 14 investigated the dynamics of a sliding flexible link through a prismatic joint. Pan et al. 15 developed the modeling for a two-link flexible robot with a prismatic joint. Korayem et al. 11 addressed the dynamic modeling of a multilink flexible robot, of which each link rotates and reciprocates through a prismatic joint. Khadem and Pirmohammadi 10 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 RR ), the two prismatic joints configuration (C PP ), the revolute-prismatic joints configuration (C RP ), and the prismatic-revolute joints configuration (C PR ).
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 16 and Dwivedy and Eberhard 4 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][18][19][20][21][22][23][24][25][26][27][28] when considering the dynamic modeling and analysis of a flexible robot. The large deformation of a flexible link was studied by Celentano. 17 The dynamic modeling methods were developed by Li et al. 21 and Wei et al. 22 The vibration of the flexible robots was focused in a number of investigations. [18][19][20]23 Also, the control issues were addressed in the recent works. [24][25][26][27][28] In the literature, the inverse dynamic analysis for flexible robots has been studied by a number of authors. Kwon and Book 29 used AMM for the modeling of the linear inverse dynamics of a single flexible link robot with one rotational joint. Asada et al. 30 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. 31 and Trautt and Bayo 32 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 33 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 34 addressed the inverse dynamics of a flexible robot with all rotary joints. Ledesma and Bayo 35 studied a nonrecursive algorithm for solving the inverse dynamic problem of a flexible robot with revolute joints also. Zhaocai and Yueqing 36 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 1 and L 2 are the lengths of the links, respectively. m 1 is the mass of link 1, and m 2 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 1 and J 2 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 PR , C RP , C RR , and C PP , which are shown in Figure 2. For example, in the configuration C PR , the first joint is a P joint, and the second one is a R joint.
In Figure 1, X OY is a reference frame, X 1 O 1 Y 1 is a link frame attached to the end point of the link 1, and ð Þ and q i t ð Þ are the translational and rotational joint variables denoted (see Table 1). z 1 t ð Þ and z 2 t ð Þ are the external generalized forces applying on J 1 and J 2 correspondingly. For the configurations, F i t ð Þ and t 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.  Notice that for C PR and C PP , OX and O 1 X 1 could not coincide, there could exist an angle a 0 between the two axes. Also, for C RP and C PP , there could exist an angle b 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 ¼ 1On, has four elastic displacements at the two nodes that are the flexural displacements u 2jÀ1 ; u 2jþ1 À Á and the slope displacements u 2j ; u 2jþ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 where T and P are the total kinetic energy and the total potential energy of the system. The vector F ¼ z 1 z 2 0 : : 0 0 ½ T is the external generalized forces acting along with components of the vector of the generalized coordinates Q ¼ q 1 q 2 u 1 u 2 :::

The equation of motion can be expressed as
where the structural damping matrix D and the Coriolis and centrifugal matrix C are correspondingly calculated as The notations a and b are the damping ratios of the system (Zhang et al.). 37 The size of the global mass matrix M and the global stiffness matrix K is 2n þ 4 ð ÞÂ 2n þ 4 ð Þ. Consider an element j of the link 2, a point r 2j x ð Þ represented in X 2 O 2 Y 2 , where 0 x l e and l e ¼ L 2 = n , can be computed for the configurations as follows 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 where Q j ¼ u 2jÀ1 u 2j u 2jþ1 u 2jþ2 ½ T is the vector of elastic displacements of the element j. Therefore, the point r 2j represented in X OY can be expressed as Notice that the end point of the link 2, χ ¼ x E y E ½ 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 where Notice that Q js and Q je are the s'th and e'th element of Q jg , respectively. It can be shown that M j is of the form  Table 1. Joint displacements and applied force/torque.
For C PR model, the other entries of the elemental matrix M j are calculated as follows m 11 ¼ m 2 l e ; m 12 ¼ À 1 12 m 2 l e 6u 2jÀ1 þ 6u 2jþ1 þ l e u 2j À l e u 2jþ2 À Á sinq 2 þ 6l e 1 À 2j ð Þcosq 2 Â Ã ; For C RP model, the other entries of the elemental matrix M j are calculated as follows m 2 l e u 2j l e À u 2jþ2 l e þ 6u 2jÀ1 l e þ 6u 2jþ1 l e À Á ; For C RR model, the other entries of the elemental matrix M j are calculated as follows

My and Bien
For C PP model, the other entries of the elemental matrix M j are calculated as follows The kinetic energy of the link 2 is yielded as 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 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 The total elastic potential energy of the whole system is yielded as 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 PR and C RR 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 RP and C PP 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 2n nodal coordinates in the vector of the generalized coordinates, Q ¼ q 1 q 2 u 1 u 2 ::: u 2nþ1 u 2nþ2 ½ T . Given a trajectory of the end-effector, χ t ð Þ ¼ x E t ð Þ y E t ð Þ ½ 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 e of the trajectory. For a given interval of time 0; T ½ , a time stepwise is selected as At a time t i 2 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 d.
The trajectory error d is approximated as the distance from the point χ C t i ð Þ to the velocity vector _ χ t i ð Þ (tangent vector) of the curve χ t i ð Þ at the time t i . Therefore, it can be calculated as Step 4. Estimate the interval F U , in between which the feasible value of F t ð Þ can be found. Corresponding to Dχ, compute DΘ ¼ Dq 1 Dq 2 ½ T and its time derivatives D _ Θ and D € Θ via the inverse kinematics of the rigid model. Substituting DΘ, D _ Θ, and D € Θ into equation (2) (neglecting the elastic deformation) yields DF. Estimate F U as where the weight k is a positive integer. To ensure that F U is a large enough interval, the selected value of k should be greater than 2.
Step 5. Modify F r t i ð Þ to refine F t i ð Þ. Loop the following substeps: Step 6. Finally, F t i ð Þ ¼ F M . All the values F t 1 ð Þ; F t 2 ð Þ; . . . ; 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 Figure 3 describes the algorithm implementation steps.

Numerical simulations
In this section, a numerical simulation of the inverse dynamics for the two configurations of the robot C PR and C RP is presented. The structural parameters of the robots, the desired trajectories, and the other simulation parameters are given in Table 2.
The results of the simulation for C PR 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 q 2 t ð Þ for both the rigid robot and the flexible robot. Figure 5 shows the resulting curves F 1 t ð Þ and t 2 t ð Þ. F r 1 t ð Þ and t r 2 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 RP are shown in Figures 7-9. 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     Simulation time (s) T d ¼ 4 T d ¼ 6   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 e between the computed trajectory and the desired trajectory. When choosing a small e, the accuracy of the inverse dynamic responses is increased, however, the computation speed is decreased.

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.