# Dynamic Modeling and Development of Symbolic Calculation Software forN-DOF Flexible-Link Manipulators Incorporating Lumped Mass

Authors: Feng Yu; Xinyuan Chen
Journal: Shock and Vibration
 (2019)

---

## Abstract

Considering the complex dynamic modeling of multi-DOF planar flexible manipulators, a general-purpose method for the rigid-flexible coupling dynamic modeling ofN-DOF flexible manipulators is proposed in this paper, and symbolic calculation software is developed. The modeling method is based on the Lagrange equation and assumed mode method (AMM). First, the N-DOF flexible manipulator is divided into two parts, which are assumed to be rigid and flexible. On this basis, the rigid part and the flexible part are coupled, and the calculation process of the model is further simplified. Then, the simplest general symbolic expression of the dynamic model of the N-DOF flexible manipulator is obtained with the induction method. According to the modeling method, “symbolic expression computation software for dynamic equations of N-DOF flexible manipulators” is developed using the symbolic calculation software Mathematica. Finally, the dynamic modeling method and the symbolic calculation software are verified by a trajectory tracking experiment with a PD control applied to a 2-DOF flexible manipulator. Compared with the traditional modeling method, the calculation time can be reduced by more than 90% using the modeling method proposed in this paper, which significantly reduces the complexity of the modeling process.


---

## Body

## 1. Introduction

Robotic manipulators are widely used in environments that are hazardous to humans and in situations where there are substantial amounts of repetitive work to reduce the burden on humans and to assist in accomplishing special tasks. To ensure sufficient positioning accuracy at the end of a manipulator, the links of the manipulator can generally be stiffened by a bulky structure to prevent the links from vibrating and to avoid elastic deflection during operation [1–4]. However, heavy links limit the speed of the manipulator’s movement and consume excessive amounts of energy. In some cases, thickening the structure of the links to ensure the positioning accuracy of the end of the manipulator may not be practical. For example, the weight of a manipulator operating on the space station must be as light as possible because of launch costs [5]. In addition, to ensure that the manipulator can perform in a greater workspace, the boom of a concrete pump truck structure should be more slender [6]. Robotic manipulators are also used in nuclear maintenance [7, 8], for spraying large aircraft [9], etc. Under these circumstances, these robotic manipulators have to be designed as flexible-link manipulators (FLMs) with a slender structure.FLMs have been found in extensive applications because of their advantages in terms of being lightweight and having high load ratios (the ratio of load mass to manipulator mass), high movement speeds, and low energy consumption [10–12]. Moreover, in human-robot interactions, a lightweight FLM is safer to humans and has smaller workload requirements on joint drivers, sensors [13, 14], etc. However, such an FLM experiences vibration and elastic deformation during movement, which will greatly affect the stability of the entire system and the positioning accuracy of the end [15, 16]. Therefore, it is important to study the effects of the dynamic characteristics of FLMs, in terms of the vibration suppression and elastic deformation compensation of the flexible links. For example, the dynamic model of a single-joint flexible manipulator has been developed to suppress the vibration of the flexible manipulator [3, 13].Because multijoint FLMs are highly nonlinear, strongly coupled, multi-degree-of-freedom (DOF) systems, the dynamic modeling is complex. In addition, considering the payload mass at the end of the manipulator and the lumped mass (for example, the mass of the actuators) on each link, multi-DOF rigid-flexible coupling dynamic modeling is more complicated. These challenging topics have attracted a large number of excellent researchers, who have conducted a considerable amount of research on the dynamic modeling of FLMs. In Giorgio and Del Vescovo’s work [17], the lumped parameter method is used to model the dynamics of flexible links, and although this modeling method is simple, the number of elements must be as large as possible to obtain an adequate model accuracy. The modeling and control of a special class of single-link, an arbitrary number of lumped masses, flexible arms is presented in [18], in which the arm consists of massless flexible structures. Subudhi and Morris adopted the Lagrange equation and the assumed mode method (AMM) to establish the dynamic model of a manipulator with flexible links and joints [19]. A dynamic model of 3-PRR parallel manipulator with flexible links based on AMM was proposed in [20]. Wei et al. provided a global mode method for dynamic modeling of a flexible arm flexible joint manipulator with a tip mass [21]. Some other modeling methods, such as the finite element method (FEM), are also commonly used. For example, the dynamic modeling of the flexible links using FEM is presented in [22]. The FEM is used to derive the dynamic equations for finding the maximum allowable dynamic load of flexible-link mobile manipulators [23]. In addition, the Newton–Euler method is also used to dynamically model the flexible arms [24]. Nonlinear dynamics of motion are developed for flexible manipulator links consisting of rotary joints that connect pairs of flexible links [25]. Bascetta et al. developed a closed-form dynamic model of flexible manipulators based on the Newton–Euler equations [26]. However, the dynamic modeling in the above studies only considers flexible arms with a certain number of DOFs or flexible links without payload mass, and the expressions of the dynamic model are still very computationally expensive in practice.Considering the complexity of the dynamic modeling of multi-DOF FLMs with a lumped mass on the links, a general-purpose dynamic modeling method forN-DOF FLMs based on the Lagrange equation is presented in this paper, and the symbolic expression computation software for the dynamic equations of N-DOF flexible manipulators is developed. In a more general situation, where the lumped mass is located not only on the end of the links but also at other positions of the links (for example, the flexible arms composed of functionally graded materials (FGM) in the fields of spacecraft, nuclear industries, robotics, etc., the flexibility of the arms can be increased or changed, or changing the vibration frequency of the flexible arm structure by adding concentrated mass to the intermediate position of flexible arms, see [27]) is also considered in this paper. The modeling method is as follows. First, the flexible arms are divided into imaginary rigid and flexible parts and modeled separately. Then, the rigid and flexible parts are coupled, and the dynamic model is further simplified. Finally, the simplest general expression of the dynamic model of the flexible arm with N-DOFs is obtained by the induction method, and “symbolic expression computation software for dynamic equations of N-DOF flexible manipulators” is developed using the symbolic calculation software Mathematica. The test results indicate that compared with the traditional modeling method, the proposed modeling method can reduce the computation time by more than 90% and the software can quickly calculate and generate the dynamic model of the flexible manipulator, thereby greatly reducing the difficulty of the dynamic modeling of the flexible manipulator. Finally, the position tracking control experiment of a 2-DOF flexible manipulator is taken as an example and the dynamic model is generated automatically by the program. Meanwhile, a physical model is also constructed. In an experiment on the end effector of the manipulator tracking a trajectory of a circle, the numerical simulation results of the dynamic model are compared with the experimental results. The dynamic model and the dynamic modeling method are both verified to be correct.The remainder of this paper is organized as follows. In Section2, the motion of N-DOF FLM systems is described. Then, in Section 3, the dynamic equations of N-DOF FLMs are derived and established. “The symbolic expression computation software for dynamic equations of N-DOF flexible manipulators” is developed in Section 4, in which the dynamic equations of two-DOF FLMs are automatically generated by the program, and the experiments are conducted. Finally, conclusions are drawn in Section 5.


## 2. Motion of Flexible Links

To study the general situation of dynamic modeling for serial multijoint FLMs and to simplify the derivation process, fixed-base and small-deformation flexible manipulators moving in the vertical plane are the focus of this paper. Before we derive and establish the dynamic model of the multi-DOF FLM, the appropriate float frame for each flexible link should be assigned, as shown in Figure1.Figure 1
Deformation of linki represented in the float frame i and inertial reference frame 0.The fixed base of the manipulator is typically treated as the global inertial reference frame and denoted0, and the i-th flexible link is associated with the float frame i. The x-axis of the float frame i passes through the start and end points of the i-th flexible link, and its origin coincides with the end point of the i-th flexible link. Assuming that the length of the i-th flexible link without deformation is Li, the angle between the i-th and i − 1-th links along the x-axis is θi. Thus, according to reference [28], the homogeneous transformation matrix from frame i to the coordinate system i−1 is derived as(1)Aii−1=cosθi−sinθiLicosθisinθicosθiLisinθi001.In equation (1), the superscript i−1 of the homogeneous transformation matrix Aii−1 represents the reference frame i−1 and the subscript i represents the current frame i. Thus, the position vector of an arbitrary point Pi on link i can be expressed in the global inertial frame 0 as(2)ri0=Ai0rii,where Ai0=A10A21⋯Aii−1, rii is the position vector of point Pi on link i in float frame i, and rii=xiyi1T.According to the aforementioned assumption, the flexible link here is considered to be an Euler beam with a constant cross section. For a small deformation of the link, because the deflection of the flexible links in the length direction is very small, we can ignore it and only consider the transverse deformation. Because the exact mathematical description and calculation of the flexible link deformation are complex given that the flexible link deformation is infinite dimensional, a finite-dimensional deformation is generally introduced to approximate the link deformation. The assumed mode method (AMM) is adopted in this paper. In the AMM, the elastic displacementsyi are generally described by a truncated modal series function and time-dependent mode amplitudes (for example, see [29]) as(3)yixi,t=∑k=1sφikxifikt,where φik is a k-th-order mode shape function, fik are generalized coordinates with respect to the φik of link i, and s is the number of modes used to describe the deflection of link i.Generally, this method has sufficient approximation accuracy for small deformations of flexible links when we truncateyi at s=2. Thus, in the following, we only consider the first second-order mode shape function and its generalized coordinates. For the Euler beam with clamped boundary conditions at both ends of the link, we choose the mode shape function φik=sinkπxi/Li,k=1,2. Thus, equation (3) can be rewritten as(4)yi=sinπxiLifi1+sin2πxiLifi2.The velocity vector ofPi in frame 0 can be expressed as(5)r˙i0=ddtri0=ddtAi0rii=A˙10A21⋯Aii−1rii+A10A˙21⋯Aii−1rii+⋯+A10⋯A˙ii−1rii+Ai0r˙ii=∑j=1i∂Ai0∂θjθ˙jrii+Ai0r˙ii.In general, the elements in matrixAi0 are a function of θ1θ2⋯θi and the elements in vector rii are a function of fi1fi2.Note that the partial derivative ofAi0 with respect to θi can be calculated as(6)∂Aii−1∂θi=QiAii−1,where, for a revolute joint, the matrix Qi is defined as(7)Qi=0−10100000.The advantage of usingQi is that we can transform the differential operation of ∂Aii−1/∂θi into a constant matrix Qi multiplying the Aii−1 matrix, which can still be used; then, we can obtain(8)∂Ai0∂θj=A10A21⋯Aj−1j−2QjAjj−1⋯Aij−1,j≤i,0,j>i.To simplify the derivation of dynamic equations of multi-DOF FLMs, we can divide the flexible link system into two parts: an imaginary rigid part and an imaginary flexible part. The expressions include the variablesθi, representing the rigid part, and the variables fi1 or fi2, representing the flexible part. Now, we imagine that the x-axis of frame i from the start point to the end point of flexible link i is a virtual rigid link i, which can be regarded as a rigid link of the flexible link i without deformation; thus, equation (8) can be interpreted as the effect of the motion of joint j on all the points on virtual rigid link i. To simplify the notations, we define(9)Uij=∂Ai0∂θj=Aj−10QjAij−1,j≤i,0,j>i.Thus, equation (5) can be rewritten as(10)r˙i0=∑j=1iUijθ˙jrii+Ai0r˙ii.For multi-DOF manipulators, the dynamic equations are a strongly coupled system; thus, considering the interaction effects between virtual rigid links, we can obtain(11)Uijk=∂Uij∂θk=Aj−10QjAk−1j−1QkAik−1,j≤k≤i,Ak−10QkAj−1k−1QjAij−1,k≤j≤i,0,i<j∪i<k.Equation (11) can be interpreted as the interaction effects of the motions of joint j and joint k on all points on virtual link i.


## 3. Derivation of Dynamic Equations for Flexible Links with Lumped Mass

### 3.1. Calculation of Kinetic Energy

After obtaining the velocity of points on linki, we can calculate the kinetic energy of link i. Assuming that the mass of flexible link i is mi, let Tim be the kinetic energy of link i and let dTim be the kinetic energy of a particle with differential mass dm on link i; then, we can obtain that(12)dTim=12Trr˙i0r˙i0Tdm=12Tr∑j=1iUijθ˙jrii+Ai0r˙ii×∑j=1iUijθ˙jrii+Ai0r˙iiTdm=12∑p=1i∑r=1iTrUipriiriiTdmUirTθ˙pθ˙r+∑j=1iTrUijriir˙iiTdmAi0Tθ˙j+12TrAi0r˙iir˙iiTdmAi0T,where Tr is a trace operator and TrX≜∑i=1nxii. Now, let us compute the kinetic energy of lumped mass TiM on flexible link i. We assume that the lumped mass can locate at any position on flexible link i. Similarly, we use the same method as equations (10) and (12) and obtain(13)TiM=12TrR˙i0R˙i0TMi=12∑p=1i∑r=1iTrUipRiiRiiTMiUirTθ˙pθ˙r+∑j=1iTrUijRiiR˙iiTMiAi0Tθ˙j+12TrAi0R˙iiR˙iiTMiAi0T,where Mi is the lumped mass and Ri0 and R˙i0 are the position vector and velocity vector of the lumped mass Mi in the global inertial reference frame 0, respectively. Rii and R˙ii are the position vector and velocity vector of lumped mass Mi in float frame i, respectively. Thus, we let Rii=XiYi1T, as shown in Figure 2.Figure 2
Position of lumped massMi with different value of μi.Because we ignore the deformation of the flexible link in the longitudinal direction, the length of flexible linki is always Li. To simplify the derivation, we define Xi=μiLi, μi∈−1,0. When μi=0, lumped mass Mi locates on the end of flexible link i if μi=−1/2 represents Mi at the middle of flexible link i. Thus, we can obtain Rii=μiLiYi1T, and according to equation (4), we can obtain(14)Yi=sinπμifi1+sin2πμifi2=∑k=12fikSki.To simplify the notation, we defineSki=sinkπμi, k=1,2. The total kinetic energy of flexible link i and lumped mass Mi can be computed as(15)Ti=∫−Li0dTim+TiM=12∑p=1i∑r=1iTrUipJuiUirTθ˙pθ˙r+∑j=1iTrUijJviAi0Tθ˙j+12TrAi0JwiAi0T,where Jui=Juim+JuiM, Jvi=Jvim+JviM, and Jwi=Jwim+JwiM. Additionally, Juim is the pseudoinertia matrix of flexible link i, which, for a homogeneous Euler beam, can be calculated as(16)Juim=∫−Li0riiriiTdm=∫−Li0xiyi1xiyi1dm=∫−Li0xi2xiyixixiyiyi2yixiyi1dm=Li23Li2fi1−fi22π−Li2Li2fi1−fi22π12fi12+fi22−2fi1π−Li2−2fi1π1mi.As shown, the pseudoinertia matrixJuim is a function of generalized coordinates f. The variation in the generalized coordinates f can reflect the deformation of the flexible link. Similarly, JuiM is calculated as(17)JuiM=RiiRiiTMi=XiYi1XiYi1Mi=Xi2XiYiXiXiYiYi2YiXiYi1Mi=Li2μi2Liμi∑k=12fikSkiLiμiLiμi∑k=12fikSki∑k=12fikSki2∑k=12fikSkiLiμi∑k=12fikSki1Mi.Note thatx˙i=0, and we can easily obtain Jvim and Jwim:(18)Jvim=∫−Li0riir˙iiTdm=0Li2f˙i1−f˙i22π0012∑k=12fikf˙ik00−2f˙i1π0mi,Jwim=∫−Li0r˙iir˙iiTdm=000012f˙i12+f˙i220000mi.J
v
i
M and JwiM are computed in the same manner:(19)JviM=RiiR˙iiTMi=0Liμi∑k=12fikSki00∑k=12fikSki∑k=12f˙ikSki00∑k=12f˙ikSki0Mi,JwiM=R˙iiR˙iiTMi=0000∑k=12f˙ikSki20000Mi.


### 3.2. Calculation of Potential Energy

Let the gravitational potential energy of flexible linki with lumped mass Mi be Pi¯; thus,(20)Pi¯=−migr¯i0+MigRi0=−migAi0r¯ii+MigAi0Rii,where g=gxgy0=0−g0 is a gravity vector expressed in the inertial reference frame 0 and g is the gravitational constant. Because of the small deformation of the flexible link, the variation in the gravitational potential energy of the flexible link caused by the flexible link’s deformation can be neglected and is small compared with the total gravitational potential energy of the flexible link. Let r¯i0 and r¯ii be the mass center vector of flexible link i without deformation in inertial frame 0 and float frame i, respectively.Meanwhile, the elastic potential energy of flexible linki caused by deformation can be expressed as(21)Ki=EiIi2∫−Li0∂2yi∂xi22dxi=π4EiIifi12+16fi224Li3,where Ei and Ii are Young’s modulus and the second moment of the area of flexible link i, respectively. Now, from equations (12) to (21), we define the Lagrangian function of FLMs with N-DOFs as(22)L=∑i=1nTi−Ki−Pi¯=12∑i=1n∑p=1i∑r=1iTrUipJuiUirTθ˙pθ˙r+∑i=1n∑j=1iTrUijJviAi0Tθ˙j+12∑i=1nTrAi0JwiAi0T+∑i=1nmigAi0r¯ii+MigAi0Rii−12∑i=1nπ4EiIifi12+16fi222Li3.We now apply the Lagrange–Euler formulation to the Lagrangian function of theN-DOF FLM equation (22) to obtain the necessary generalized torque to drive the manipulators:(23)ddt∂L∂q˙−∂L∂q=τ,where q and q˙ are a vector of generalized coordinates and its derivative with respect to time, respectively. We define generalized coordinates q=θfT, θ=θ1⋯θn, and f=f11f12⋯fn1fn2, and τ is a 3n×1 generalized torque vector.Before we apply the Lagrange–Euler formulation to the Lagrangian function of equation (22), we should note that the transformation matrix Ai0 and the matrix Uij are only functions of generalized coordinates θi,i=1,2,…,n; thus, the derivation of Ai0 or Uij with respect to fiλ,i=1,2,…,n,λ=1,2 will be zero. The matrices Jui, Jvi, and Jwi are only functions of the generalized coordinates fiλ,λ=1,2 and generalized velocity f˙iλ,λ=1,2; thus, the derivatives of Jui, Jvi, and Jwi with respect to θi,i=1,2,…,n will be zero.First, we only consider the generalized coordinateθ; thus, take θ1⋯θnT as q and substitute equation (22) into equation (23) to obtain the rigid generalized torque vector Fθ:(24)Fθi=∑j=in∑k=1jTrUjkJujUjiTθ¨k+∑j=in∑λ=12TrUjiJvjλAj0Tf¨jλ+∑j=in∑k=1j∑m=1jTrUjkmJujUjiTθ˙kθ˙m+∑j=in∑k=1j∑λ=12TrUjkJujλUjiTθ˙kf˙jλ−∑j=inmjgUjir¯jj−∑j=inMjgUjiRjj,where i=1,2,…,n. Then, we derive the flexible generalized torque vector Ff by substituting equation (22) into equation (23) when considering f11f12⋯fn1fn2T as q(25)Ffiλ=∑j=1iTrUijJviλAi0Tθ¨j+12mif¨iλ+∑j=1i∑k=1iTrUijkJviλAi0T+UijJviλ−12JuiλUikTθ˙jθ˙k+π4EiIi2−λfiλ+16λ−1fiλ2Li3−MigAi0Riλ.Fori=1,2,…,n,λ=1,2, where Riλ is a partial derivative of Rii with respect to the generalized coordinate fiλ and Riλ=0sinλπμi0T.Equations (24) and (25) are the dynamic equations of FLMs with N-DOFs. Because the entire process of deriving the dynamic equations from equations (23) to (25) is very complex, we divide the equations into two parts: a rigid part with respect to the generalized coordinate θ and a flexible part with respect to the generalized coordinate f. Therefore, the complexity is greatly reduced. Note that in equations (24) and (25), we should calculate the derivative of the matrices Jui, Jvi, and Jwi with respect to time t because(26)dJxjdt=∂Jxj∂fjλf˙jλ,where x=u,v,w and j=1,2,…,nandλ=1,2 because the ∂Jwj/∂fjλ is canceled during the calculation, and thus, we define(27)Jujλ=∂Juj∂fjλ=∂Jujm∂fjλ+∂JujM∂fjλ,Jvjλ=∂Jvj∂fjλ=∂Jvjm∂fjλ+∂JvjM∂fjλ.It is easy to calculateJujλ and Jvjλ according to the expressions of Juj and Jvj for j=1,2,…,n,λ=1,2:(28)Jujλ=2−λ0Ljπ0Ljπfjλ−2π0−2π0mj+2−λ0Ljμj0Ljμj2∑k=12fjkSkj1010S1j+λ−10−Lj2π0−Lj2πfjλ0000mj+λ−10Ljμj0Ljμj2∑k=12fjkSkj1010S2j,Jvjλ=2−λ0Ljπ00fjλ200−2π0mj+2−λ0Ljμj00∑k=12fjkSkj0010S1j+λ−10−Lj2π00fjλ20000mj+λ−10Ljμj00∑k=12fjkSkj0010S2j.From equations (24) and (25), we find that there are not only rigid motion terms with θ¨i but also flexible motion terms with f¨iλ, in addition to the coupling terms between rigid motion θ˙iθ˙j and rigid motion and flexible link deformation θ˙if˙iλ. Notably, there are no flexible deformation coupling terms, which should have been f˙iλf˙jλ. This makes solving the dynamic equations easier, and this is an important reason for establishing such float frames to describe the FLMs. In essence, the two equations also reflect the vibration of the FLMs during movement.


### 3.3. Relationship between the Generalized Forces and Actual Driving Forces

We continue to discuss the relationship between the generalized forces corresponding to the generalized coordinates and the actual forces applied to the manipulator’s joints. As shown in Figure3, τi,i−1 is the actual force applied to joint i, where the first subscript i denotes the link that the force acts on and the second subscript i−1 denotes the link that the force reacts on. According to the principle of virtual work, the virtual work acting on link i is(29)δWi=τi,i−1δθi+δαi+τi+1,iδβi,where αi=∂yi/∂xixi=0 and βi=∂yi/∂xixi=Li; thus, the total virtual work on all joints is(30)δW=∑i=1nτi,i−1δθi+∑λ=12τn,n−1φiλ′0δfnλ+∑i=1n−1∑λ=12τi,i−1φiλ′0+τi,i+1φiλ′Liδfiλ,where φiλ′0=dφiλ/dxixi=0 and φiλ′Li=dφiλ/dxixi=Li, and the generalized forces corresponding to the generalized coordinates θi are(31)Fθi=τi,i−1,i=1,2,…,n.Figure 3
Torques applied at linki.The generalized forces corresponding to the generalized coordinatesfiλ are(32)Ffiλ=τi,i−1φiλ′0+τi+1,iφiλ′Li,Ffnλ=τn,n−1φnλ′0,where i=1,2,…,n−1,λ=1,2.


### 3.4. Matrix Form of Dynamic Equations for FLM

To make the dynamic equations easy to understand, we rewrite the equations from equations (24) to (31) in matrix form as(33)MθθMθfMfθMffθ¨f¨+HθθHθfHf02n×2nθ˙f˙+CrKf+Cf=FθFf,where Mθθ∈ℜn×n, Mθf∈ℜn×2n, Mfθ∈ℜ2n×n, and Mff∈ℜ2n×2n are mass matrices and Hθθ∈ℜn×n, Hθf∈ℜn×2n, and Hf∈ℜ2n×n are Coriolis and centrifugal matrices. Cr∈ℜn×1 and Cf∈ℜ2n×1 are gravity vectors. K∈ℜ2n×2n is a diagonal stiffness matrix. The detailed expressions are as follows:(34)Mθθi,k=∑j=maxi,knTrUjkJujUjiT,where i=1,2,…,nandk=1,2,…,n and(35)Mθfi,j,λ=TrUjiJvjλAj0T,where i=j=1,2,…,n,λ=1,2. Because we only take the first two mode shape functions and generalized coordinates, we can regard every two-column vector as one group in the matrix Mθf∈ℜn×2n, and “λ” means the λ-th column vector in the j-th group. Hence, Mθfi,j,λ represents the i-th row, j-th group, and λ-th element in the matrix Mθf:(36)Mfθi,j,λ=TrUijJviλAi0T,where i=j=1,2,…,n,λ=1,2. From the expression, we find that the matrix Mfθ is the transpose of Mθf:(37)Mffi,i=12mj,where i=1,2,…,2n,j=2i+1+−1i+1/4, and mff∈ℜ2n×2n is a diagonal matrix.(38)Hθθi,k=∑m=1nhθθi,k,mθ˙m,where i=1,2,…,n,k=1,2,⋯,n, with(39)hθθi,k,m=∑j=maxi,k,mnTrUjkmJujUjiT,Hθfi,j,λ=∑k=1jTrUjkJujλUjiTθ˙k,where i=j=1,2,…,n,λ=1,2. The computation of the matrix Hθf is the same as that of Mθf:(40)Hfi,j,λ=∑k=1iTrUijkJviλAi0T+UijJviλ−12JuiλUikTθ˙k,where i=j=1,2,…,n,λ=1,2. The calculation of the matrix Hf is the same as that of the matrix Mfθ:(41)Cri=−∑j=inmjgUjir¯jj−∑j=inMjgUjiRjj,where i=1,2,…,n.(42)Ki,i=diagπ4EiIi2Li3,16π4EiIi2Li3,where i=1,2,…,n. K∈ℜ2n×2n is a diagonal matrix, and the element Ki,i is also a 2×2 diagonal matrix.(43)Cfi,λ=−MigAi0Riλ,where i=1,2,…,n,λ=1,2. Here, we also regard every two-row vector as one group, similar to equation (35). Thus, Cfi,λ represents the i-th group, λ-th element of the matrix Cf.We now rewrite equation (33) in a more compact and general matrix form as(44)Mq¨+Hq,q˙q˙+Nq=F,where M∈ℜ3n×3n is the mass matrix, the matrix H∈ℜ3n×3n consists of the centrifugal force and Coriolis force, the column vector N∈ℜ3n×1 denotes the gravitational and elastic force, and F∈ℜ3n×1 is a generalized force.


### 3.5. Symbolic Expression Calculation Procedure

The complete symbolic expression of the rigid-flexible coupling dynamic equation of theN-DOF flexible manipulator with the lumped mass on an arbitrary location of links has been obtained. According to the whole derivation process of this paper and from equations (33) to (44), software that can automatically generate the rigid-flexible coupling dynamic equation of an N-DOF flexible manipulator is developed by using Wolfram Mathematica symbolic calculation software. The operating interface of the software is shown in Figure 4.Figure 4
Operating interface of symbolic expressions calculation software.In the program interface, only the number ofN-DOFs, the arm length, the section inertia moment, the elastic modulus, the mass of the flexible arm, the symbolic parameters of the end mass, and the location of the lumped mass are input. According to the program, the symbolic expression of the dynamic equation of the flexible manipulator with N-DOFs can be obtained conveniently and quickly. Here, we rearrange the symbolic expression computation procedures as follows:(i)
Step 0. Initialization of the algorithmThe number of DOFs and the parameters of the FLM, such asN,mi,Mi,Li,Ei,Ii,θi,μi, are predefinedComputation ofJui,Jvi,Jwi,Juiλ,Juiλ,Q,g are predefinedAii−1 and Ai0 from (1); ri0 from (2)Uij from (9); Uijk from (11)(ii)
Step 1. Computation of the global mass matrixComputation ofMθθ from (34); Mθf from (35)Mfθ from (36); Mff from (37)(iii)
Step 2. Computation of the global matrices of the centrifugal and Coriolis forcesComputation ofHθθ from (38); Hθf from (39)Hf from (40)(iv)
Step 3. Computation of the global vector of the gravity, stiffness matrix, and generalized forceComputation ofCr from (41); Cf from (43)K from (42)Fθ from (31); Ff from (32)(v)
Step 4. Postprocessing of the resultSimplicity of the dynamic model: some terms can be ignored for special situationsOutput of the symbolic expression for the model.The flowchart of the symbolic expression computation procedures of the dynamic equation of the flexible manipulator with N-DOFs is shown in Figure5.Figure 5
Flowchart of symbolic expressions calculation procedures.


### 3.6. Comparison with Traditional Modeling Methods

The traditional dynamic modeling methods of flexible manipulators include the Lagrange method, Newton–Euler method, and Kane method, which have their advantages and disadvantages. In this paper, the traditional modeling method is compared with the proposed method, which is also based on the Lagrange equation. The traditional modeling method is as shown in [19].The traditional modeling method is relatively simple, but it is quite complex in terms of the actual calculation process for solving the dynamic model expression and there is no general formula for flexible manipulators with different DOFs. Therefore, it often takes substantial time to calculate the dynamic model of various types of flexible manipulators in practice.Here, the program developed in Section3.5 is tested. Additionally, according to the modeling method in [19], the calculation programs for flexible manipulators with different DOFs developed using Mathematica are also tested. The time required for the two modeling methods to run the computational program (on a PC with a i3-7100 3.90 GHz CPU with 8 GB of RAM) was tested when N = 1, 2, …, 5. The comparison results are shown in Table 1.Table 1
Required computation time of two modeling methods for different DOF models.
DOF
Proposed method (s)
Traditional method (s)
Time saved (%)
N
=
1
t
=
0.36
t
=
41.63
99.12
N
=
2
t
=
0.92
t
=
106.28
99.17
N
=
3
t
=
4.21
t
=
212.57
98.01
N
=
4
t
=
21.3
t
=
365.28
94.14
N
=
5
t
=
114.6
t
=
1316.45
91.29As shown in Table1, the modeling method proposed in this paper can reduce the modeling time by more than 90% compared with the traditional modeling method. This result is mainly because in the traditional modeling process, the Lagrange energy function of the manipulator system requires a large number of integral operations; then, in the derivation of the generalized coordinates, there are many differential operations. Moreover, there are many repeated intermediate calculations in the differential calculations. The method proposed in this paper eliminates the repetitive intermediate calculation process theoretically and directly derives the simplest expression of the dynamic model equation.Note that to explain the modeling method proposed in this paper, the simplest planar flexible manipulator is selected as the research object. However, from the entire modeling process shown above, the flexibility characteristics of rigid-flexible coupling manipulators are related only to the matricesJui, Jvi, and Jwi and their derivatives. Thus, for the dynamic modeling of complex spatial flexible manipulators, large-deformation flexible arms, or even irregularly shaped flexible links, we only need to modify the matrices Jui, Jvi, and Jwi and their derivatives. The structure of the symbolic expressions remains unchanged; thus, this dynamic modeling approach also has good portability. It can be predicted that to satisfy the increasingly complex working conditions, the structure of the flexible arm will be more complex, and it will be a difficult problem to establish a sufficiently accurate dynamic model. Therefore, the modeling method proposed in this paper will have substantial advantages in terms of modeling convenience and computational efficiency.


## 3.1. Calculation of Kinetic Energy

After obtaining the velocity of points on linki, we can calculate the kinetic energy of link i. Assuming that the mass of flexible link i is mi, let Tim be the kinetic energy of link i and let dTim be the kinetic energy of a particle with differential mass dm on link i; then, we can obtain that(12)dTim=12Trr˙i0r˙i0Tdm=12Tr∑j=1iUijθ˙jrii+Ai0r˙ii×∑j=1iUijθ˙jrii+Ai0r˙iiTdm=12∑p=1i∑r=1iTrUipriiriiTdmUirTθ˙pθ˙r+∑j=1iTrUijriir˙iiTdmAi0Tθ˙j+12TrAi0r˙iir˙iiTdmAi0T,where Tr is a trace operator and TrX≜∑i=1nxii. Now, let us compute the kinetic energy of lumped mass TiM on flexible link i. We assume that the lumped mass can locate at any position on flexible link i. Similarly, we use the same method as equations (10) and (12) and obtain(13)TiM=12TrR˙i0R˙i0TMi=12∑p=1i∑r=1iTrUipRiiRiiTMiUirTθ˙pθ˙r+∑j=1iTrUijRiiR˙iiTMiAi0Tθ˙j+12TrAi0R˙iiR˙iiTMiAi0T,where Mi is the lumped mass and Ri0 and R˙i0 are the position vector and velocity vector of the lumped mass Mi in the global inertial reference frame 0, respectively. Rii and R˙ii are the position vector and velocity vector of lumped mass Mi in float frame i, respectively. Thus, we let Rii=XiYi1T, as shown in Figure 2.Figure 2
Position of lumped massMi with different value of μi.Because we ignore the deformation of the flexible link in the longitudinal direction, the length of flexible linki is always Li. To simplify the derivation, we define Xi=μiLi, μi∈−1,0. When μi=0, lumped mass Mi locates on the end of flexible link i if μi=−1/2 represents Mi at the middle of flexible link i. Thus, we can obtain Rii=μiLiYi1T, and according to equation (4), we can obtain(14)Yi=sinπμifi1+sin2πμifi2=∑k=12fikSki.To simplify the notation, we defineSki=sinkπμi, k=1,2. The total kinetic energy of flexible link i and lumped mass Mi can be computed as(15)Ti=∫−Li0dTim+TiM=12∑p=1i∑r=1iTrUipJuiUirTθ˙pθ˙r+∑j=1iTrUijJviAi0Tθ˙j+12TrAi0JwiAi0T,where Jui=Juim+JuiM, Jvi=Jvim+JviM, and Jwi=Jwim+JwiM. Additionally, Juim is the pseudoinertia matrix of flexible link i, which, for a homogeneous Euler beam, can be calculated as(16)Juim=∫−Li0riiriiTdm=∫−Li0xiyi1xiyi1dm=∫−Li0xi2xiyixixiyiyi2yixiyi1dm=Li23Li2fi1−fi22π−Li2Li2fi1−fi22π12fi12+fi22−2fi1π−Li2−2fi1π1mi.As shown, the pseudoinertia matrixJuim is a function of generalized coordinates f. The variation in the generalized coordinates f can reflect the deformation of the flexible link. Similarly, JuiM is calculated as(17)JuiM=RiiRiiTMi=XiYi1XiYi1Mi=Xi2XiYiXiXiYiYi2YiXiYi1Mi=Li2μi2Liμi∑k=12fikSkiLiμiLiμi∑k=12fikSki∑k=12fikSki2∑k=12fikSkiLiμi∑k=12fikSki1Mi.Note thatx˙i=0, and we can easily obtain Jvim and Jwim:(18)Jvim=∫−Li0riir˙iiTdm=0Li2f˙i1−f˙i22π0012∑k=12fikf˙ik00−2f˙i1π0mi,Jwim=∫−Li0r˙iir˙iiTdm=000012f˙i12+f˙i220000mi.J
v
i
M and JwiM are computed in the same manner:(19)JviM=RiiR˙iiTMi=0Liμi∑k=12fikSki00∑k=12fikSki∑k=12f˙ikSki00∑k=12f˙ikSki0Mi,JwiM=R˙iiR˙iiTMi=0000∑k=12f˙ikSki20000Mi.


## 3.2. Calculation of Potential Energy

Let the gravitational potential energy of flexible linki with lumped mass Mi be Pi¯; thus,(20)Pi¯=−migr¯i0+MigRi0=−migAi0r¯ii+MigAi0Rii,where g=gxgy0=0−g0 is a gravity vector expressed in the inertial reference frame 0 and g is the gravitational constant. Because of the small deformation of the flexible link, the variation in the gravitational potential energy of the flexible link caused by the flexible link’s deformation can be neglected and is small compared with the total gravitational potential energy of the flexible link. Let r¯i0 and r¯ii be the mass center vector of flexible link i without deformation in inertial frame 0 and float frame i, respectively.Meanwhile, the elastic potential energy of flexible linki caused by deformation can be expressed as(21)Ki=EiIi2∫−Li0∂2yi∂xi22dxi=π4EiIifi12+16fi224Li3,where Ei and Ii are Young’s modulus and the second moment of the area of flexible link i, respectively. Now, from equations (12) to (21), we define the Lagrangian function of FLMs with N-DOFs as(22)L=∑i=1nTi−Ki−Pi¯=12∑i=1n∑p=1i∑r=1iTrUipJuiUirTθ˙pθ˙r+∑i=1n∑j=1iTrUijJviAi0Tθ˙j+12∑i=1nTrAi0JwiAi0T+∑i=1nmigAi0r¯ii+MigAi0Rii−12∑i=1nπ4EiIifi12+16fi222Li3.We now apply the Lagrange–Euler formulation to the Lagrangian function of theN-DOF FLM equation (22) to obtain the necessary generalized torque to drive the manipulators:(23)ddt∂L∂q˙−∂L∂q=τ,where q and q˙ are a vector of generalized coordinates and its derivative with respect to time, respectively. We define generalized coordinates q=θfT, θ=θ1⋯θn, and f=f11f12⋯fn1fn2, and τ is a 3n×1 generalized torque vector.Before we apply the Lagrange–Euler formulation to the Lagrangian function of equation (22), we should note that the transformation matrix Ai0 and the matrix Uij are only functions of generalized coordinates θi,i=1,2,…,n; thus, the derivation of Ai0 or Uij with respect to fiλ,i=1,2,…,n,λ=1,2 will be zero. The matrices Jui, Jvi, and Jwi are only functions of the generalized coordinates fiλ,λ=1,2 and generalized velocity f˙iλ,λ=1,2; thus, the derivatives of Jui, Jvi, and Jwi with respect to θi,i=1,2,…,n will be zero.First, we only consider the generalized coordinateθ; thus, take θ1⋯θnT as q and substitute equation (22) into equation (23) to obtain the rigid generalized torque vector Fθ:(24)Fθi=∑j=in∑k=1jTrUjkJujUjiTθ¨k+∑j=in∑λ=12TrUjiJvjλAj0Tf¨jλ+∑j=in∑k=1j∑m=1jTrUjkmJujUjiTθ˙kθ˙m+∑j=in∑k=1j∑λ=12TrUjkJujλUjiTθ˙kf˙jλ−∑j=inmjgUjir¯jj−∑j=inMjgUjiRjj,where i=1,2,…,n. Then, we derive the flexible generalized torque vector Ff by substituting equation (22) into equation (23) when considering f11f12⋯fn1fn2T as q(25)Ffiλ=∑j=1iTrUijJviλAi0Tθ¨j+12mif¨iλ+∑j=1i∑k=1iTrUijkJviλAi0T+UijJviλ−12JuiλUikTθ˙jθ˙k+π4EiIi2−λfiλ+16λ−1fiλ2Li3−MigAi0Riλ.Fori=1,2,…,n,λ=1,2, where Riλ is a partial derivative of Rii with respect to the generalized coordinate fiλ and Riλ=0sinλπμi0T.Equations (24) and (25) are the dynamic equations of FLMs with N-DOFs. Because the entire process of deriving the dynamic equations from equations (23) to (25) is very complex, we divide the equations into two parts: a rigid part with respect to the generalized coordinate θ and a flexible part with respect to the generalized coordinate f. Therefore, the complexity is greatly reduced. Note that in equations (24) and (25), we should calculate the derivative of the matrices Jui, Jvi, and Jwi with respect to time t because(26)dJxjdt=∂Jxj∂fjλf˙jλ,where x=u,v,w and j=1,2,…,nandλ=1,2 because the ∂Jwj/∂fjλ is canceled during the calculation, and thus, we define(27)Jujλ=∂Juj∂fjλ=∂Jujm∂fjλ+∂JujM∂fjλ,Jvjλ=∂Jvj∂fjλ=∂Jvjm∂fjλ+∂JvjM∂fjλ.It is easy to calculateJujλ and Jvjλ according to the expressions of Juj and Jvj for j=1,2,…,n,λ=1,2:(28)Jujλ=2−λ0Ljπ0Ljπfjλ−2π0−2π0mj+2−λ0Ljμj0Ljμj2∑k=12fjkSkj1010S1j+λ−10−Lj2π0−Lj2πfjλ0000mj+λ−10Ljμj0Ljμj2∑k=12fjkSkj1010S2j,Jvjλ=2−λ0Ljπ00fjλ200−2π0mj+2−λ0Ljμj00∑k=12fjkSkj0010S1j+λ−10−Lj2π00fjλ20000mj+λ−10Ljμj00∑k=12fjkSkj0010S2j.From equations (24) and (25), we find that there are not only rigid motion terms with θ¨i but also flexible motion terms with f¨iλ, in addition to the coupling terms between rigid motion θ˙iθ˙j and rigid motion and flexible link deformation θ˙if˙iλ. Notably, there are no flexible deformation coupling terms, which should have been f˙iλf˙jλ. This makes solving the dynamic equations easier, and this is an important reason for establishing such float frames to describe the FLMs. In essence, the two equations also reflect the vibration of the FLMs during movement.


## 3.3. Relationship between the Generalized Forces and Actual Driving Forces

We continue to discuss the relationship between the generalized forces corresponding to the generalized coordinates and the actual forces applied to the manipulator’s joints. As shown in Figure3, τi,i−1 is the actual force applied to joint i, where the first subscript i denotes the link that the force acts on and the second subscript i−1 denotes the link that the force reacts on. According to the principle of virtual work, the virtual work acting on link i is(29)δWi=τi,i−1δθi+δαi+τi+1,iδβi,where αi=∂yi/∂xixi=0 and βi=∂yi/∂xixi=Li; thus, the total virtual work on all joints is(30)δW=∑i=1nτi,i−1δθi+∑λ=12τn,n−1φiλ′0δfnλ+∑i=1n−1∑λ=12τi,i−1φiλ′0+τi,i+1φiλ′Liδfiλ,where φiλ′0=dφiλ/dxixi=0 and φiλ′Li=dφiλ/dxixi=Li, and the generalized forces corresponding to the generalized coordinates θi are(31)Fθi=τi,i−1,i=1,2,…,n.Figure 3
Torques applied at linki.The generalized forces corresponding to the generalized coordinatesfiλ are(32)Ffiλ=τi,i−1φiλ′0+τi+1,iφiλ′Li,Ffnλ=τn,n−1φnλ′0,where i=1,2,…,n−1,λ=1,2.


## 3.4. Matrix Form of Dynamic Equations for FLM

To make the dynamic equations easy to understand, we rewrite the equations from equations (24) to (31) in matrix form as(33)MθθMθfMfθMffθ¨f¨+HθθHθfHf02n×2nθ˙f˙+CrKf+Cf=FθFf,where Mθθ∈ℜn×n, Mθf∈ℜn×2n, Mfθ∈ℜ2n×n, and Mff∈ℜ2n×2n are mass matrices and Hθθ∈ℜn×n, Hθf∈ℜn×2n, and Hf∈ℜ2n×n are Coriolis and centrifugal matrices. Cr∈ℜn×1 and Cf∈ℜ2n×1 are gravity vectors. K∈ℜ2n×2n is a diagonal stiffness matrix. The detailed expressions are as follows:(34)Mθθi,k=∑j=maxi,knTrUjkJujUjiT,where i=1,2,…,nandk=1,2,…,n and(35)Mθfi,j,λ=TrUjiJvjλAj0T,where i=j=1,2,…,n,λ=1,2. Because we only take the first two mode shape functions and generalized coordinates, we can regard every two-column vector as one group in the matrix Mθf∈ℜn×2n, and “λ” means the λ-th column vector in the j-th group. Hence, Mθfi,j,λ represents the i-th row, j-th group, and λ-th element in the matrix Mθf:(36)Mfθi,j,λ=TrUijJviλAi0T,where i=j=1,2,…,n,λ=1,2. From the expression, we find that the matrix Mfθ is the transpose of Mθf:(37)Mffi,i=12mj,where i=1,2,…,2n,j=2i+1+−1i+1/4, and mff∈ℜ2n×2n is a diagonal matrix.(38)Hθθi,k=∑m=1nhθθi,k,mθ˙m,where i=1,2,…,n,k=1,2,⋯,n, with(39)hθθi,k,m=∑j=maxi,k,mnTrUjkmJujUjiT,Hθfi,j,λ=∑k=1jTrUjkJujλUjiTθ˙k,where i=j=1,2,…,n,λ=1,2. The computation of the matrix Hθf is the same as that of Mθf:(40)Hfi,j,λ=∑k=1iTrUijkJviλAi0T+UijJviλ−12JuiλUikTθ˙k,where i=j=1,2,…,n,λ=1,2. The calculation of the matrix Hf is the same as that of the matrix Mfθ:(41)Cri=−∑j=inmjgUjir¯jj−∑j=inMjgUjiRjj,where i=1,2,…,n.(42)Ki,i=diagπ4EiIi2Li3,16π4EiIi2Li3,where i=1,2,…,n. K∈ℜ2n×2n is a diagonal matrix, and the element Ki,i is also a 2×2 diagonal matrix.(43)Cfi,λ=−MigAi0Riλ,where i=1,2,…,n,λ=1,2. Here, we also regard every two-row vector as one group, similar to equation (35). Thus, Cfi,λ represents the i-th group, λ-th element of the matrix Cf.We now rewrite equation (33) in a more compact and general matrix form as(44)Mq¨+Hq,q˙q˙+Nq=F,where M∈ℜ3n×3n is the mass matrix, the matrix H∈ℜ3n×3n consists of the centrifugal force and Coriolis force, the column vector N∈ℜ3n×1 denotes the gravitational and elastic force, and F∈ℜ3n×1 is a generalized force.


## 3.5. Symbolic Expression Calculation Procedure

The complete symbolic expression of the rigid-flexible coupling dynamic equation of theN-DOF flexible manipulator with the lumped mass on an arbitrary location of links has been obtained. According to the whole derivation process of this paper and from equations (33) to (44), software that can automatically generate the rigid-flexible coupling dynamic equation of an N-DOF flexible manipulator is developed by using Wolfram Mathematica symbolic calculation software. The operating interface of the software is shown in Figure 4.Figure 4
Operating interface of symbolic expressions calculation software.In the program interface, only the number ofN-DOFs, the arm length, the section inertia moment, the elastic modulus, the mass of the flexible arm, the symbolic parameters of the end mass, and the location of the lumped mass are input. According to the program, the symbolic expression of the dynamic equation of the flexible manipulator with N-DOFs can be obtained conveniently and quickly. Here, we rearrange the symbolic expression computation procedures as follows:(i)
Step 0. Initialization of the algorithmThe number of DOFs and the parameters of the FLM, such asN,mi,Mi,Li,Ei,Ii,θi,μi, are predefinedComputation ofJui,Jvi,Jwi,Juiλ,Juiλ,Q,g are predefinedAii−1 and Ai0 from (1); ri0 from (2)Uij from (9); Uijk from (11)(ii)
Step 1. Computation of the global mass matrixComputation ofMθθ from (34); Mθf from (35)Mfθ from (36); Mff from (37)(iii)
Step 2. Computation of the global matrices of the centrifugal and Coriolis forcesComputation ofHθθ from (38); Hθf from (39)Hf from (40)(iv)
Step 3. Computation of the global vector of the gravity, stiffness matrix, and generalized forceComputation ofCr from (41); Cf from (43)K from (42)Fθ from (31); Ff from (32)(v)
Step 4. Postprocessing of the resultSimplicity of the dynamic model: some terms can be ignored for special situationsOutput of the symbolic expression for the model.The flowchart of the symbolic expression computation procedures of the dynamic equation of the flexible manipulator with N-DOFs is shown in Figure5.Figure 5
Flowchart of symbolic expressions calculation procedures.


## 3.6. Comparison with Traditional Modeling Methods

The traditional dynamic modeling methods of flexible manipulators include the Lagrange method, Newton–Euler method, and Kane method, which have their advantages and disadvantages. In this paper, the traditional modeling method is compared with the proposed method, which is also based on the Lagrange equation. The traditional modeling method is as shown in [19].The traditional modeling method is relatively simple, but it is quite complex in terms of the actual calculation process for solving the dynamic model expression and there is no general formula for flexible manipulators with different DOFs. Therefore, it often takes substantial time to calculate the dynamic model of various types of flexible manipulators in practice.Here, the program developed in Section3.5 is tested. Additionally, according to the modeling method in [19], the calculation programs for flexible manipulators with different DOFs developed using Mathematica are also tested. The time required for the two modeling methods to run the computational program (on a PC with a i3-7100 3.90 GHz CPU with 8 GB of RAM) was tested when N = 1, 2, …, 5. The comparison results are shown in Table 1.Table 1
Required computation time of two modeling methods for different DOF models.
DOF
Proposed method (s)
Traditional method (s)
Time saved (%)
N
=
1
t
=
0.36
t
=
41.63
99.12
N
=
2
t
=
0.92
t
=
106.28
99.17
N
=
3
t
=
4.21
t
=
212.57
98.01
N
=
4
t
=
21.3
t
=
365.28
94.14
N
=
5
t
=
114.6
t
=
1316.45
91.29As shown in Table1, the modeling method proposed in this paper can reduce the modeling time by more than 90% compared with the traditional modeling method. This result is mainly because in the traditional modeling process, the Lagrange energy function of the manipulator system requires a large number of integral operations; then, in the derivation of the generalized coordinates, there are many differential operations. Moreover, there are many repeated intermediate calculations in the differential calculations. The method proposed in this paper eliminates the repetitive intermediate calculation process theoretically and directly derives the simplest expression of the dynamic model equation.Note that to explain the modeling method proposed in this paper, the simplest planar flexible manipulator is selected as the research object. However, from the entire modeling process shown above, the flexibility characteristics of rigid-flexible coupling manipulators are related only to the matricesJui, Jvi, and Jwi and their derivatives. Thus, for the dynamic modeling of complex spatial flexible manipulators, large-deformation flexible arms, or even irregularly shaped flexible links, we only need to modify the matrices Jui, Jvi, and Jwi and their derivatives. The structure of the symbolic expressions remains unchanged; thus, this dynamic modeling approach also has good portability. It can be predicted that to satisfy the increasingly complex working conditions, the structure of the flexible arm will be more complex, and it will be a difficult problem to establish a sufficiently accurate dynamic model. Therefore, the modeling method proposed in this paper will have substantial advantages in terms of modeling convenience and computational efficiency.


## 4. Experimental Study

To verify the correctness of the modeling method proposed in this paper, the dynamic model of a 2-DOF flexible manipulator is automatically generated using the software. Simultaneously, a flexible manipulator experimental platform is set up, as shown in Figure6. Figure 7 is a schematic diagram of the assembly structure of the flexible arm device.Figure 6
Physical model of 2-DOF flexible manipulators.Figure 7
Schematic diagram of the assembly structure of 2-DOF flexible manipulators.The expressions of the dynamic equations for the 2-DOF flexible manipulator are shown in the appendix. The symbols in the formula are the same as those in the derivation process.In this experiment, the two arms of the flexible manipulator are controlled to follow the given trajectories. The torqueτref output from the controller and the angles q of the joints of the flexible manipulator measured by the encoders are recorded. Simultaneously, the angles are transmitted to the dynamics simulation model, and the estimations of the control torques τest are calculated for comparison with the actual situation. The whole experimental scheme and recorded data are shown in Figure 8. Figure 9 is a schematic diagram of the experimental signal acquisition and control.Figure 8
Experimental scheme and data recording.Figure 9
Schematic diagram of the experimental signal acquisition and control.According to the scheme shown in Figure8, the Simulink control model is established and compiled into a “.so file,” the input and output (IO) interface between the control model and the real manipulator is configured using VeriStand, and finally, all the related control model files are downloaded to the real-time cRIO-9035 controller. The output torque calculated by the controller is converted into a PWM signal by dividing the motor torque constant and the gear ratio of the reducer. A 40 MHz FPGA chip is mounted onto the back of the control slot to acquire the high-speed pulse signal produced by the encoders through the digital IO interface. The motors, gears, and drivers are manufactured by Times Brilliant in China. The type of brushless DC motor 1 with the power of 150 W and brushless DC motor 2 with the power of 100 W is 57BL95S15-230TF9B and 57BL75S10-230TF9, respectively. The type of driver 1 and driver 2 is ZM-6610M. Reduction ratios of PX series planetary gears matched with motors are 64 and 32, respectively. The other specific parameters of the components in the experimental device are listed in Table 2.Table 2
Parameters of components in the experimental devices.
Parameters
Values
Length of armsL1,L2m
0.475
,
0.435
Mass of armsm1,m2kg
0.374
,
0.29
Young’s modulus of armsE1,E2GPa
2.06
×
10
11
Second moment of area of armsI1m4
16
×
10
−
11
Second moment of area of armsI2m4
6.75
×
10
−
11
Lumped mass of jointsM1,M2Kg
2.94
,
0.343
Torque values of motorsτ1,τ2Nm
0.5
,
0.32
Gear reduction ratior1,r2
64,32
Coulomb friction of actuatorskc1,kc2Nm
0.0345
,
0.0226
Viscous friction of actuatorskv1,kv2Nms/rad
0.0037
,
0.0023


### 4.1. Case 1: Free Vibration of the End

In the experiment, the initial working positions of the flexible manipulator are selected asθ1=45∘ and θ2=−45∘. At the beginning, the external support makes the manipulator as free as possible from gravity and the second flexible manipulator keeps the horizontal position. Then, the flexible manipulator system is allowed to vibrate freely under the force of gravity, and the acceleration of the end of the flexible manipulator under free vibration is observed and recorded. In addition, the accelerometer type is WT901C485, the acquisition time and sampling time are 10 s and 20 ms, respectively, and the simulation interval of Simulink is 0.001 s. The measured acceleration graph is compared with the simulation graph, as shown in Figure 10(a). It is apparent that there are very good correlations between the measured values and simulated values and the amplitude of simulation is smaller, which can be due to model approximation and initial attitude error between the simulation model and the physical model. The spectrum of the end acceleration curve of the simulation model is obtained by fast Fourier transform, as shown in Figure 10(b). It can be seen that the first-order vibration frequency is ω1 = 12.246 rad/s and the second-order vibration frequency is ω2 = 40.82 rad/s. From the experiment, it is obvious that the simulation curves are basically consistent with the experimental curves, which verifies the correctness of the proposed modeling method to some extent. The similar results of vibration characteristics of flexible manipulators can be referred to [30], which can be used to compare with these experimental results.Figure 10
(a) Acceleration of flexible link’s tip. (b) The spectrum of the flexible link’s tip acceleration.
(a)
(b)


### 4.2. Case 2: End Effector Tracking a Trajectory of Circle

In this experiment, the initial working positions of the flexible manipulator are selected asθ1=45∘ and θ2=−45∘. A circle with a radius of 0.1 m in the XOY plane is selected as the end track of the flexible arm, as shown in Figure 11. The equations of the circle are x=x0+0.1cos0.2πt, y=y0+0.1sin0.2πt, whereas the center coordinate of the circle is x0,y0, x0=0.67m, y0=0.34m. The end of the flexible arm is substituted into the inverse equation of rigid kinematics, and the desired joint angles are obtained by real-time operation. In addition, the simple PD control is adopted in the cRIO-9035 controller to make the joints track the desired angles, in which kp and kd were debugged through repeated simulation experiments and finally determined to be kp=125 and kd=25, and the calculation interval is 0.001 s. Figures 12(a) and 12(b) are the angle curves of joint 1 and joint 2. It can be seen from the figures that when the flexible arm starts to move, the experimental curve and the simulation curve have some errors, and then the errors gradually decreases and the curves tracking effect is better. The input torque curves of the joint 1 and the joint 2 are as shown in Figures 13(a) and 13(b), and the vertical acceleration curves of the flexible arm 1 and the flexible arm 2 are as shown in Figures 14(a) and 14(b). It can be seen that the torque simulation curves and the experimental curves follow the same trend, and there are small tracking errors, due to the inevitable vibration of the flexible mechanical arm during the tracking process, the acceleration will cause frequent fluctuations, but the simulation curve and the experimental curve can also achieve a better tracking effect on the whole. The main cause of some vibration of the manipulator during movement is that, to simplify the control model, the distance between the output shafts of two motors is regarded as the arm length of the flexible arms, the mass of the motor reducers is treated as the lumped mass, and the unavoidable clearance of the motor reducers makes the simulation model and the actual model have certain errors, but these errors are within acceptable ranges.Figure 11
Reference trajectory in the workspace.Figure 12
(a) Angle of joint 1. (b) Angle of joint 2.
(a)
(b)Figure 13
(a) Torque of joint 1. (b) Torque of joint 2.
(a)
(b)Figure 14
(a) Acceleration of tip of link 1. (b) Acceleration of tip of link 2.
(a)
(b)


## 4.1. Case 1: Free Vibration of the End

In the experiment, the initial working positions of the flexible manipulator are selected asθ1=45∘ and θ2=−45∘. At the beginning, the external support makes the manipulator as free as possible from gravity and the second flexible manipulator keeps the horizontal position. Then, the flexible manipulator system is allowed to vibrate freely under the force of gravity, and the acceleration of the end of the flexible manipulator under free vibration is observed and recorded. In addition, the accelerometer type is WT901C485, the acquisition time and sampling time are 10 s and 20 ms, respectively, and the simulation interval of Simulink is 0.001 s. The measured acceleration graph is compared with the simulation graph, as shown in Figure 10(a). It is apparent that there are very good correlations between the measured values and simulated values and the amplitude of simulation is smaller, which can be due to model approximation and initial attitude error between the simulation model and the physical model. The spectrum of the end acceleration curve of the simulation model is obtained by fast Fourier transform, as shown in Figure 10(b). It can be seen that the first-order vibration frequency is ω1 = 12.246 rad/s and the second-order vibration frequency is ω2 = 40.82 rad/s. From the experiment, it is obvious that the simulation curves are basically consistent with the experimental curves, which verifies the correctness of the proposed modeling method to some extent. The similar results of vibration characteristics of flexible manipulators can be referred to [30], which can be used to compare with these experimental results.Figure 10
(a) Acceleration of flexible link’s tip. (b) The spectrum of the flexible link’s tip acceleration.
(a)
(b)


## 4.2. Case 2: End Effector Tracking a Trajectory of Circle

In this experiment, the initial working positions of the flexible manipulator are selected asθ1=45∘ and θ2=−45∘. A circle with a radius of 0.1 m in the XOY plane is selected as the end track of the flexible arm, as shown in Figure 11. The equations of the circle are x=x0+0.1cos0.2πt, y=y0+0.1sin0.2πt, whereas the center coordinate of the circle is x0,y0, x0=0.67m, y0=0.34m. The end of the flexible arm is substituted into the inverse equation of rigid kinematics, and the desired joint angles are obtained by real-time operation. In addition, the simple PD control is adopted in the cRIO-9035 controller to make the joints track the desired angles, in which kp and kd were debugged through repeated simulation experiments and finally determined to be kp=125 and kd=25, and the calculation interval is 0.001 s. Figures 12(a) and 12(b) are the angle curves of joint 1 and joint 2. It can be seen from the figures that when the flexible arm starts to move, the experimental curve and the simulation curve have some errors, and then the errors gradually decreases and the curves tracking effect is better. The input torque curves of the joint 1 and the joint 2 are as shown in Figures 13(a) and 13(b), and the vertical acceleration curves of the flexible arm 1 and the flexible arm 2 are as shown in Figures 14(a) and 14(b). It can be seen that the torque simulation curves and the experimental curves follow the same trend, and there are small tracking errors, due to the inevitable vibration of the flexible mechanical arm during the tracking process, the acceleration will cause frequent fluctuations, but the simulation curve and the experimental curve can also achieve a better tracking effect on the whole. The main cause of some vibration of the manipulator during movement is that, to simplify the control model, the distance between the output shafts of two motors is regarded as the arm length of the flexible arms, the mass of the motor reducers is treated as the lumped mass, and the unavoidable clearance of the motor reducers makes the simulation model and the actual model have certain errors, but these errors are within acceptable ranges.Figure 11
Reference trajectory in the workspace.Figure 12
(a) Angle of joint 1. (b) Angle of joint 2.
(a)
(b)Figure 13
(a) Torque of joint 1. (b) Torque of joint 2.
(a)
(b)Figure 14
(a) Acceleration of tip of link 1. (b) Acceleration of tip of link 2.
(a)
(b)


## 5. Conclusions

A general modeling method for the dynamic equations of anN-DOF flexible manipulator based on the Lagrange equation and the small deformation assumption of an Euler beam is proposed in this paper. According to this modeling method, symbolic calculation software for the dynamic equations of N-DOF flexible manipulators is developed, greatly reducing the computational complexity. Compared with the traditional modeling method, the modeling time can be reduced by more than 90%. Finally, the correctness of the modeling method is verified by experiments.The modeling method proposed in this paper takes the planar flexible manipulator as the research object, but according to the modeling process, it can be applied to the more complex rigid-flexible coupling dynamic modeling of flexible manipulators only through modifying the corresponding matrices of the flexible part. This is also future work that needs to be studied. Simultaneously, the modeling method proposed in this paper is also helpful for the dynamic modeling ofN-DOF spatial flexible manipulators, large deformation flexible manipulators, etc.



---

*Source: 5627271-2019-07-01.xml*