Kinematics and Dynamics of Flexible Robotic Manipulators Using Dual Screws

In this article, we present a new procedure for the derivation of the linearized kinematics and dynamics of a flexible industrial robotic manipulator. We introduce the Lie groups of dual rotation and dual homogeneous transformation matrices, which are the basis for the derivation of dual twists. In addition, dual screws and dual screw transformation matrices are introduced, which are used for the development of a general and systematic linearization procedure for both kinematics and dynamics. This leads to expressions that are linearized in all the states associated with the elastic motion. The dynamic modeling procedure is based on Kane's method, where the partial velocities and partial angular velocities are given as dual screws arranged as columns of dual projection matrices. The elasticity of the links is modeled by the assumed mode method. The presented procedure is implemented numerically for a 4-degrees of freedom manipulator and the simulation results are given.


I. INTRODUCTION
R OBOTIC manipulators are widely used in industrial applications. In some applications manipulators need to be lightweight and need to have large configuration space. Such long-reach lightweight manipulators may have significant elasticity, which has to be accounted for in the derivation of manipulator kinematics and dynamics. Modeling the elasticity of links leads to more complex expressions with additional higher order terms, which have little contribution to the dynamical simulation results. Often it is reasonable to partly linearize the dynamics by removing those higher order terms, which are the products of elastic displacements and velocities. In this article, we present the systematic procedure for derivation of linearized kinematics and dynamics of flexible robotic manipulators.
Dynamic behavior of flexible manipulators has been investigated by many researchers in the recent decades. Flexible links of manipulators posses an infinite number of degrees of freedom (DOFs) and are governed by partial differential equations. The equations are normally solved by discretization, where the Manuscript received February 17, 2020; revised June 22, 2020; accepted July 28, 2020. This work was supported by the Norwegian Research Council, SFI Offshore Mechatronics, under Project 237896. This paper was recommended for publication by Associate Editor C. Rucker and Editor E. Yoshida upon evaluation of the reviewers' comments. (Corresponding author: Andrej Cibicik.) The authors are with the Department of Mechanical and Industrial Engineering, Norwegian University of Science and Technology, NO-7491 Trondheim, Norway (e-mail: andrej.cibicik@ntnu.no; olav.egeland@ntnu.no).
Color versions of one or more of the figures in this article are available online at https://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TRO.2020.3014519 assumed mode method (AMM), finite-element method (FEM), or lumped parameter method [1] are normally used [2]. Modeling flexibility of manipulators by FEM can lead to accurate results, however large number of finite elements and DOFs may be needed. The advantage of FEM is the ability of modeling complex systems with nonuniform link geometry [3]. Alternatively, using AMM for modeling the flexibility may result in a smaller number of DOFs. In [4], the authors compare AMM and FEM formulations for a flexible manipulator and conclude that for manipulators with uniform cross-sectional area of the links AMM leads to a system with fewer equations and is a preferred option. However, the problem of choosing the proper mode shapes might be complex for manipulators with several links and variable cross-sections. In [5], the authors discussed the two main alternatives for mode shapes in terms of boundary conditions, that is pinned-free (or pinned-mass) and clamped-free (or clamped-mass) modes. It was found that clamped boundary condition at the controlled end of a link leads to significant simplifications in the design of the controlled. Similarly in [6], it was concluded that a clamped boundary condition at the controlled end of a link leads to more accurate predictions of closed-loop behavior. In [4], [7], the authors suggested that eigenmodes of a clamped-free Euler-Bernoulli beam are often used for flexible serial manipulators with high hub-to-link mass ratio, however, it was highlighted that using clamped-mass boundary conditions is more correct. AMM can also be applied for modeling flexible robotic manipulators with telescopic joints [8].
The assumption of small deformations is often adopted in flexible manipulator models, which leads to the elimination of high-order terms related to elastic displacements and velocities. The dynamic equation of motion linearized in the terms associated with elastic motion for a multilink planar manipulator was presented in [7], [9], and [10], where the Euler-Lagrange formulation and AMM was used. In most of the reviewed references, the small deformation assumption was adopted, however no explicit discussion on the linearization procedure was presented.
Kane's equation of motion (or Kane's method) [11] is a variation of the Newton-Euler formulation, where a minimal set of ODE's (ordinary differential equations) is obtained by elimination of workless constraint forces. The constraint forces are eliminated by d'Alembert's principle and principle of virtual work. In practice it is done by premultiplication of the equations of motion with the partial velocities, which serve as projection operators mapping the forces on the directions of generalized speeds. Kane's method was formulated for flexible systems This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ [12]- [14], and a particular example of a flexible beams on a moving base was given in [15], where the authors also included the effect of geometric stiffening. Additionally, the authors proposed a simplified method for linearization of dynamic equations of motion in all elastic displacements and velocities.
Lie group theory [16] for modeling of robot kinematics is a well-established field [17]- [19]. This includes the application of the Lie algebra se(3) to describe the differential kinematics and dynamics of multibody systems, where derivations of differential kinematics rely on the fact that the vector representation of se(3) is a twist, which can be transformed by an adjoint map Ad SE(3) . This was used for rigid mechanisms in [20], [21], and [22], and for flexible systems in [23] and [24]. The fact that a twist is a screw given in Plücker coordinates allows for an alternative formulation in terms of screw theory [25], [26], where the geometrical interpretation of the derivations may be more evident.
The review articles [2] and [27] reveal the amount of research done in the field of modeling of flexible manipulators and flexible multibody systems. It is generally accepted that the equations of motion are highly nonlinear due to the coupling between the DOFs of the system. In this article, we propose a method for linearization of equations of motion of flexible manipulators by eliminating higher order terms formed by the products of elastic deflections and velocities. The kinematic framework of the method is based on the Lie groups of dual rotation and dual homogeneous transformation matrices, which are the basis for the derivation of dual twists. The method is an extension of the previous work [25], [28], [29] and is based on the linearization procedure proposed in [12] and [15]. The advantage of the proposed method is that we propose to use screw theory and dual linear algebra for the derivations, which is easier to implement computationally than the original formulation in [12] and [15] where coordinate-free vectors and dyadics were used. This advantage becomes more apparent for large systems with many links. Another advantage of the proposed method is that the expressions are linearized at each step of derivation. The method does not require a readily derived fully nonlinear model before the linearization is performed [30], which is a quite common approach. The method can also be used for the derivation of models which are partly nonlinear in some elastic states or velocities.
In the end of this article, we present a numerical example, where the method is implemented for a four-link manipulator with two flexible and two rigid links.

A. Dual Numbers, Dual Vectors, and Dual Matrices
A dual number [31], [32] is given byx = r x + ε d x, where r x and d x are real numbers, and the dual unit ε satisfies ε = 0 and ε 2 = 0. A dual vector of dimension n is given byv where r v and d v are vectors of dimension n. In the same way, a dual n × n matrix is defined byD = r D + ε d D, where r D and d D are n × n matrices. Calculations on dual numbers, vectors and matrices are done by treating the expressions as polynomials in ε, and using ε 2 = 0.
Duality can be applied in linearization. In this article, we will in particular use this by defining dual matrices, where the real part is a zero-order term and the dual part is a first-order term. Then, the product of two dual matrices will have a real part which is of order zero, and a dual part which is of order one.
In this article, we define elastic coordinates and velocities as dual parts. Then, the products of elastic coordinates and velocities vanish from the equations, after each multiplication operation during the model derivation. For example, letp 1 = r 1 + εe 1 andp 2 = r 2 + εe 2 be some position vectors given in the dual from, where r i is the rigid displacement, and e i is the elastic deformation. Then, the product of those two terms is found asp 1p2 = r 1 r 2 + ε(r 1 e 2 + e 1 r 2 ), which is linear in e i as the term ε 2 e 1 e 2 is vanished.

B. Matrix Lie Groups
Consider a matrix Lie group G ⊂ GL(n, R), where GL(n, R) is the linear group of invertible n × n matrices with real elements [16], [33]. Then, if G,Ḡ ∈ G, it follows that GḠ ∈ G, and there is an inverse G −1 ∈ G so that GG −1 = I ∈ G, where I is the identity matrix of the appropriate dimension. Let U be an element of the associated Lie algebra g ⊂ M(n, R), where M(n, R) is the algebra of n × n matrices with real elements. Let U ∈ g be the logarithm of G, which means that G = exp(U). The logarithmic map is the inverse map of the exponential map, and gives U = log(G). The logarithm U can be represented by the vector u ∈ R n . The conversion between U ∈ g and u ∈ R n is done with the notation U = The time derivative of a Lie group element G ∈ G is given byĠ = V l G = GV r , where V l ∈ g is the left velocity [33], which is also called the spatial velocity [18], and V r ∈ g is the right velocity, which is also called the body velocity. The corresponding velocity vectors v l = [V l ] ∨ ∈ R n and v r = [V r ] ∨ ∈ R n are v l = Ψ l (ad(u))u and v r = Ψ r (ad(u))u, where Ψ l is the left Jacobian and Ψ r is the right Jacobian of G [34].

C. Twists and Screw Transformations
The displacement of Frame b relative to Frame a is described by the homogeneous transformation matrix T a b ∈ SE(3) [17]. The time derivative of the homogeneous transformation matrix is given byṪ a (3) is the right velocity. The transformation from the right to the left velocity is given by The left and right velocities can be represented by the twist vectors t a ab/a ∈ R 6 and t b ab/b ∈ R 6 from a to b, which transform according to is a matrix representation of the adjoint map in SE (3). A twist vector is in general given by [34] t i ab/r = where ω i ab is the angular velocity of b relative to a, and v i ab/r is the velocity of b relative to a referenced to r. The superscript i indicated that the vectors are given in the coordinates of Frame i.
The twist vector is a screw [32], [34], which satisfies the screw transformation with the screw transformation matrix V ij rs =R i j U j rs . This transformation involves a change of coordinates, which is due tō and a change of the reference point, which is due to where p j rs is a vector from r to s given in the coordinates of Frame j. It is straightforward to verify that the twists of the composite displacement T a c = T a b T b c can be added according to [25] t i ac/r = t i ab/r + t i bc/r (5) under the condition that the twists are referenced to the same point r, and are given in the coordinates of the same Frame i.

D. Wrenches
A wrench is a screw representation of forces and torques acting on a rigid body [32], [35]. The wrench with a line of action through the origin of Frame j, referenced to j and given in the coordinates of j is given as where f j j is a force acting at the origin of Frame j and n j j/j is a torque acting on Link j, both given in the coordinates of j. A wrench is a screw and it satisfies screw transformations. The wrench (6) can be referenced to the COG (center of gravity) of link i and expressed in the coordinates of i as

E. Projection Matrices for Rigid Manipulators
Consider two Frames j − 1 and j that are fixed in two robotic links connected by a 1-DOF revolute joint. Let the direction vector of the joint be the unit vector a j j through the origin of Frame j. Then, the line through the joint axis referenced to the origin of Frame j will be L j j/j = [(a j j ) T 0 T ] T , where the moment is zero since the line passes through the reference point. The twist of the link displacement T j−1 j can be given in terms of L j j/j as t j j−1,j/j = L j j/j u j , where u j =q j is the generalized speed of the rotation in Joint j [25], [26].
The line L j j/j can be referenced to the COG of any subsequent Link i and expressed in the coordinates of i as L i j/m i = V i,j m i ,j L j j/j . Then, the twist of Link i relative to the inertial frame, referenced to the COG of Link i and expressed in the coordinates of i can be given as t i 0i/m i = P i u, where u = [u 1 . . . u n ] T is a vector of generalized speeds [11] and the projection matrix [25], [26], [28] is given by The projection matrix (7) have two main functions in the derivation of the equations of motion. First, it serves as a geometric Jacobian [34], which maps the robot joint space velocities into Cartesian velocities of the point m i and angular velocities of Link i. Second, it serves as an operator for projecting the external and inertial wrenches on the direction of the joint lines. In the following sections of this article, we will introduce the Jacobian form and the projecting form of projection matrices, which are different in the case of linearized dynamics of flexible manipulators.

III. LINEARIZED KINEMATICS OF FLEXIBLE MANIPULATORS IN TERMS OF DUAL SCREWS
In this section, we present a novel derivation procedure for linearized kinematics of flexible manipulators, where the equations include only the first-order effects of elastic deflections and velocities. This includes the introduction of two Lie Groups: the group of dual rotation matrices and the group of dual homogeneous transformations matrices. The group velocity for the group of dual homogeneous transformation matrices is derived as a dual twist, which can be given in terms of dual screws.
The development of this section starts with the introduction of the Lie groups. At this stage, the Lie groups are given as mathematical objects without any particular reference to flexible manipulators. Then, the concepts of dual rotation and dual homogeneous transformation matrices with application to flexible manipulators are first introduced in the dual matrix form. Later it is demonstrate that dual rotation and dual homogeneous transformation matrices can alternatively be represented as members of the proposed Lie groups, that is matrices in the dual form are represented by real matrices which allows using basic linear algebra for the derivations. In addition, Lie algebras associated with the proposed Lie groups are derived in both vector and matrix forms, which are also used in the derivations.
We consider serial link manipulators with elastic links. The elastic links are described as Euler-Bernoulli beams with the assumed mode approach, where elastic deformations are modeled by a set of eigenmodes for a constrained beam with clamped-free or clamped-mass boundary conditions [7]. IEEE TRANSACTIONS ON ROBOTICS

A. Definition of the Matrix Lie Group G 1
Consider the matrix Lie group G 1 , where the group element is given as where A ∈ SO(3) and θ × ∈ R 3×3 is skew-symmetric. The notation X = (A, θ × A) ∈ G 1 will also be used for the group element in (8). Let X 1 = (A 1 , θ × 1 A 1 ) and X 2 = (A 2 , θ × 2 A 2 ) be elements of G 1 . The group operation is matrix multiplication and the product is given by is the identity element of the group, which means that XX I = X I X = X. The derivative of (8) with respect to time is given byẊ where the right velocity From where Ad G 1 (X) = X.
Consider the two matrices be the corresponding vectors. The following property is verified by direct computation As well-established Lie groups used in robotics (SO(3) and SE(3)), the group G 1 has the exponential and logarithmic maps. Although these maps are not directly used in this article, they are still derived for the sake of generality in Appendix A. It is noted that the Lie group G 1 has the same structure as the adjoint representation of SE(3). The difference will however become more apparent later in this article, where the group will be used as a dual rotation matrix for coordinate transformation of dual vectors. Then, the vector θ is a rotation angle between two frames, rather than a distance vector in SE (3). Then, the derivations and expressions in this chapter obtain a new kinematic interpretation and, therefore, are presented. Equivalent derivations for the adjoint representation of SE(3) are given in [34]. It is also noted that the R 4×4 representation of G 1 cannot be interpreted as the Euclidean transformation since the vector θ is an angle.

B. Definition of the Matrix Lie Group G 2
Consider the matrix Lie group G 2 , where the group element is given as where X ∈ G 1 andv ∈ R 6 . The notation Y = (X,v) ∈ G 2 will also be used. Let The time derivative of the group element (14) is given bẏ where the right velocity is which is in the Lie algebra g 2 . A vector representation of the right velocity ist where the notation [[t] ∧ ] ∨ =t will also be used. The adjoint map is given by This is verified from and in combination with (12) and [v] ∧ Xω = −[Xω] ∧v , which follows from (13). Similarly as for G 1 , the group G 2 has the exponential and logarithmic maps, which are not directly used in this article, they are also derived for the sake of generality in Appendix A.

C. Linearized Dual Rotation Matrices
Consider an elastic robotic link given in Fig. 1. Define the orthogonal rotation matrix R j j ∈ SO(3) from the body-fixed Frame j attached at the controlled end of Link j to the body-fixed Frame j attached to the distal end of Link j. The rotation matrix given in terms of the angle-axis parameters (θ j , k j j ) is [17] It is assumed that the elastic deformations are small θ j → 0, and therefore a first-order approximation of the rotation matrix will be used. This leads to sin(θ j ) → θ j , cos(θ j ) → 1, and linearization of (22) givesR j j = I + θ j k j× j . Assume that this linear rotation matrix is given as a product of three consecutive linearized rotations about the x j , y j , and z j axes as where θ j j = [θ x,j θ y,j θ z,j ] T and the skew-symmetric form is θ j× j ∈ so(3). A similar approximation of the rotation matrix was given in [36]. The higher order terms (h.o.t) in (23) can be omitted, and the linearized rotation matrix can be formulated as a dual matrixR j j = I + εθ j× j (24) where the zero-order term I ∈ SO(3) is the real part, the firstorder term θ j× j is the dual part. The inverse of the rotation matrix will have the angle-axis parameters (−θ j , k j j ), and it follows that the first-order approximation of the inverse is It is noted that The time derivative of the rotation matrix and its first-order approximation is given bẏ whereθ j j = ω j j,j . The first-order vector εω j j,j can be expressed in the coordinates of j as which leads to the first-order approximation ω j j,j = ω j j,j . Similarly, the first-order vector εθ j j can be expressed in the coordinates of j as which leads to another first-order approximation θ j × j = θ j× j . This elastic deflection can be combined with rigid motion described with the second dual rotation matrixR j i = R j i from j to i, where i is a frame fixed at the controlled end of Link i, which is attached to Link j through a single-DOF rotational joint. The dual rotation matrix from j to i is then given aŝ where (29) is used and

D. Linearized Dual Rotation Matrices as Elements of G 1
The dual composite rotation matrixR j i can be given as It is noted that The inverse of (31) is given as and the time derivative of (31) is given as The right velocity (10) for the dual rotation matrix (31) is found as The right velocity (34) can be given in a vector form aŝ which corresponds to the dual formω i ji = ω i j ,i + εω i j,j , where the real term is of order zero and the dual term is of order one in θ j j orθ j j . Denote the real term as r ω i ji = ω i j ,i and the dual term as d ω i ji = ω i j,j , then (35) can alternatively be written

E. Coordinate Transformation of Dual Vectors
Consider the dual vector from the origin of Frame i to some point a i fixed in Link i given in the coordinates of i asâ The term d i i,a i describes the displacement due to elastic deformations, which are assumed to be small. The dual vectorâ i i,a i can be given in R 6 aŝ and can be given in a matrix form as The vector (37) can be expressed in the coordinates of j aŝ whereR j i is given by (31). The real term of (39) is of order zero and the dual term is of order one in either θ j j or d i i,a i . It is straightforward to verify that the inverse transformation iŝ where [R j i ] −1 is given by (32). Alternatively, the dual skewsymmetric matrix (38) can be expressed in the coordinates of j as which follows from (12)

F. Linearized Dual Homogeneous Transformation Matrices
Let T = (R, p) ∈ SE(3) be the pose of an undeformed elastic link which has a body-fixed frame. Suppose that the link has an elastic deflection, such that the elastic displacement and rotation of the frame is described by the vector δ = [θ T , d T ] T , so that the deflected pose is If a first-order approximation of the elastic deflection is used, then T δ ∈ SE(3) is approximated byT δ = (I + [δ] ∧ )T, which is formulated in a dual form aŝ where the real part T ∈ SE(3) is a homogeneous transformation matrix, which is the rigid pose, and the dual part is the linearized elastic deflection [δ] ∧ T ∈ T T SE(3). The matrixT will be referred to as a dual homogeneous transformation matrix in the following.
be two dual homogeneous transformation matrices. Then, the composite dual homogeneous transformation matrix will bê Consider the dual homogeneous transformation matrixT j j from j to j given aŝ where (29) is used and the reference for the notations is made to Fig. 1. The dual vectorp j ji = r j ji + εd j ji can alternatively be expressed by a dual coordinate transformation then (45) can be written aŝ The dual homogeneous transformation matrixT j i from j to i isT Then, the dual homogeneous transformation matrixT j i is found from (45) and (48) to bê

G. Linearized Dual Homogeneous Transformation Matrices as Elements of G 2
Consider the dual homogeneous transformation matrixT j i in (49) given as The inverse of (50) is found to be and the time derivative is is used, see (27) and (28). The Lie where the following dual coordinate The Lie algebra (53) is alternatively written as which can be given in a vector form as a dual twist

IV. LINEARIZED DIFFERENTIAL KINEMATICS
In this section, we propose a method where dual screws and dual screw transformations are used to describe the differential kinematics of flexible manipulators. All dual matrices in this section are given in a matrix form, andÂ notation will for simplicity be used instead of [Â] notation as in previous sections.

A. Dual Twists of Rotary Joints and Dual Screw Transformations
Consider the system given in Fig. 2 where the adjacent elastic Links j − 1 and j are connected with Joint j, and Links j and j + 1 are connected with Joint j + 1. Link j + 1 is referred to as Link i and Joint j + 1 is referred to as Joint i to simplify the notation. The rotary Joint j is a rotation from Frame (j − 1) fixed in the distal end of Link j − 1 to Frame j fixed in the controlled end of Link j. The rotary Joint i is a rotation from Frame j fixed in the distal end of Link j to Frame i fixed in the controlled end of Link i. The dual twist of the rigid motion due to rotation in Joint j is found as in the rigid case given in [25] to bet where u j is the generalized speed of Joint j and the dual line of The real parts of the dual line are r a j j = [0 0 1] T and r m j j/j = 0, which are also the axis and the moment of the real line L j j/j as given in Section II-E, while the dual parts are d a j j = 0 and d m j j/j = 0. The dual line (58) can be transformed to be referenced to the origin of Frame j bŷ is a dual screw reference transformation matrix, which moves the point of reference from j to j , and [p j j ,j ] ∧ is a matrix form (38) of the dual vectorp j j ,j , which is given aŝ Here, r p j j ,j is a vector from j to j when the link is undeformed, and d p j j ,j is a vector from j in the deformed state to j in the undeformed state. This means that the total distance from j to j in the deformed state is p j j ,j = r p j j ,j + d p j j ,j . The dual line (59) can be expressed in the coordinates of Frame j byL is the dual screw coordinate transformation matrix, andR j j is the dual rotation matrix. The dual line (62) can be further transformed aŝ Then, the dual twist of Frame j relative to Frame (j − 1) , reference to the origin of Frame i e and expressed in the coordinates of i can be given ast

B. Dual Twists Due to Elastic Deformation
The dual twist of Frame j relative to Frame j due to elastic deformation, referenced to j and expressed in the coordinates of j is given ast where u bj = [ u bj 1 . . . u bj n ] T is the vector of generalized speeds associated with elastic motion of Link j, and is the dual screw matrix, which is given as a product of the dual screw basiŝ and the matrix of assumed mode shape functions Ψ j evaluated at the origin of Frame j . The dimensions of n = [1 . . . 1] in (69) depend on the number of assumed modes as given by Ψ j . The auxiliary matrix N interchanges the real and dual parts of the vectors in a screw, and it is defined as where it is noted that N 2 = I. For an elastic link, where only flexural and torsional elasticity is taken into account, Ψ j can be given as where (·) notation refers to a partial derivative with respect to the spatial coordinate, where n is a number of assumed modes. The dual linesL j ·/j and dual screwsN j ·/j in (69) have zero dual partŝ It is noted that assigning the coordinate axes directions to the real parts in (72) is a convention we use for the given notation. In fact, velocities associated with the elastic deformations should be defined as dual parts in (68). Therefore, N matrix is introduced in (69) to interchange the real and dual parts. Since all columns of (69) are dual screws, and Ψ j are scalar values of functions evaluated at the defined point, it follows that the columns of the dual screw matrixŜ j j /j are dual screws, which means that the matrix satisfies dual screw transformations. The dual screw matrix (68) can be transformed to be referenced to the origin of Frame i e and expressed in the coordinates of i aŝ then the dual twist of Frame j relative to Frame j, reference to the origin of Frame i e and expressed in the coordinates of i can be given ast where the dual basisÊ i i/i is defined as (69). The matrix of assumed mode functions of Link i Ψ i e is defined as in (71) and the assumed mode functions are evaluated at the origin of Frame i e . It is straightforward to extend this procedure for a case oft i 0,i e /i e , which is a dual twist of Frame i e fixed in Link i relative to the inertial frame, reference to i e and given in the coordinates of i. This will be done in the next section.

C. Derivatives of Dual Lines and Dual Screws
The time derivative of a dual line is given as an extension of the time derivative of a real line [25]. The geometric differentiation procedure given in this section is more computationally efficient for symbolic expressions than using the build-in differentiation commands in Matlab and the chain rule. The comparison of the calculation times is given in Section VI. The time derivative of the dual line of Joint j referenced to the point a i on Link i and expressed in the coordinates of i (see Fig. 3) can be expressed in terms of the dual line aṡ where will be referred to as the dual differentiation operator, andυ i a i ,j will be referred to as the auxiliary dual velocity vector, which is given aŝ (78) It is noted that if Link i is rigid, thenω k i ,j =ω k i,j , which means that the rotation of the frames at the start and the end of the link is the same.
The derivative of the dual screw matrix with respect to time can be expressed in terms of the dual screw matrix aṡ In this case, the dual differentiation operator [d i ji ] is given by (77) with the termv i j,j in (78) being equal to zero.

D. Dual Wrenches
A dual wrench is a dual screw representation of forces and torques, which is an extension of (6). Consider a dual wrencĥ which can alternatively be given aŝ where the matrixL j j/j is a dual screw basis, which is given aŝ and ρ r is a vector of magnitudes of forces and torques given as The dual lines and screws corresponding to the directions of coordinate frame axes in (82) are given by (72). The dual screw basis (82) can be transformed to be referenced to the point a i on Link i and expressed in the coordinates of i aŝ where the dual screw transformation matrix is defined in a similar manner as in (66). The dual wrenchŵ i j/a i referenced to the point a i and expressed in the coordinates of i can then be given asŵ In this article, application of dual screw transformations on dual wrenches is limited to coordinate transformations only.

V. LINEARIZED DYNAMICS
In this section, we use the proposed kinematic formulation to derive the linearized equations of motion for a flexible robotic manipulator. This is based on Kane's method [11] where dual projection matrices are used to represent the partial velocities. The model is based on linear elastic kinematics, where the linear Euler-Bernoulli beam model and engineering theory of torsion is used. The final equations of motion are derived such that they include the terms up to the order.
1) O( q e q r q e ), O( q e q r q r ), which stem from the mass matrix. 2) O( q e q n r q 2 r ), O( q n r q r q e ), which stem from the Coriolis and centrifugal force matrix, and n is any power of a rigid coordinate.
3) O( q e q n r ), which stem from the vector of gravity forces. 4) O( q e ), which stem from the stiffness matrix. In this definition of the order, we have used symbolic notation of q r for rigid DOF and q e for elastic DOF. It is noted the order above is defined for the final terms obtained by multiplication of the system matrices with the vectors of generalized accelerations, speeds, or coordinates.
The discussion on the necessity of inclusion of the higher order terms in a model based on linear elastic kinematics will be addressed in Section VI.
The model is similar to the standard models used in majority of dynamic codes and software [37]. Such models can give sufficiently accurate simulation results for flexible manipulators with slower actuation and smaller end-effector load. In the case when a manipulator is exposed to high velocities, accelerations or end-effector loads, the so-called effect of dynamic stiffening needs to be introduced to the equations of motion [36], [37]. In this article, we consider slower manipulators where dynamic stiffening effects can be neglected.
Often the dynamical models presented in the literature are linearized to a certain order. A good discussion on type of nonlinearities was given in [36]. However, what is common for linearized models is that: first, the nonlinear model is derived and, then, linearized to a desired order [30]. Some methods suggest that linearization can be done in several steps during the derivation, e.g., [15]. In this article, we propose a procedure, where both kinematic and dynamic expressions are linear to a given above order at each step of derivation.

A. Dual Projection Matrices
Consider an open-chain robotic manipulator consisting of n q links, where some of the links are assumed to be rigid, and some are assumed to be elastic. The configuration space of the manipulator is defined by the vector of generalized coordinates where q i are rigid rotations in the joints, q bk = [ q bk 1 . . . q bk nk ] T are the elastic coordinates of elastic Link k (which is the first elastic link in the chain), and n k is the number of assumed modes for Link k. The corresponding vector of generalized speeds is given as where u i =q i . The dual line of Joint j referenced to the point a i on Link i and expressed in the coordinates of i iŝ The dual screw matrix of Frame j referenced to the point a i on Link i and expressed in the coordinates of i iŝ The dual twist of Link i relative to the inertial frame, referenced to the point a i on Link i and expressed in the coordinates of i can be given as a sum The dual twist (90) can be alternatively given in terms of the Jacobian form of the projection where a i is the COG point m i for a rigid link and the origin of Frame i e for an elastic link. The projection matrixP i for Link i is defined aŝ where n q is a number of links. In addition to the Jacobian form of the dual projection matrix (92), we define the projecting form of the dual projection matrix, which serves as a projection operator of wrenches on the directions of generalized speeds. The projecting form of the dual projection matrix for Link i is defined as The dual screw basis matrix in this case isȆ i i/i = NÊ i i/i , wherê E i i/i is given by (69) and N is given by (70). The derivative of the twist (91) is given aṡ where the time derivative of the projection matrix is found by the time differentiation of the dual lines and screws contained in it, while the time derivation procedure for dual lines and screws is discussed in Section IV-C.

1) Rigid Links:
The equation of motion of a rigid Link i can be derived using the Newton-Euler approach where v i 0i/m i is the linear velocity of the COG of Link i relative to the inertial frame, ω i 0i is the angular velocity of Link i relative to the inertial frame, and M i i is the inertia tensor of Link i. The external gravity force f and which satisfies the identityÎ[ω] ∧ = [ω] ∧Î . The equation of motion (98) contains linearized form of equations (97) and, in addition, the equations formed by only real parts of the terms. The advantage of this formulation will become evident when the linearized equation of motion of the whole system will be derived.
The equation of motion (98) can alternatively be formulated using the dual screw notations aŝ The dual screwŵ 2) Elastic Links: The equation of motion of an elastic Link i is derived using the Newton-Euler approach (105) The equation of motion (105) can be formulated using the dual screw notations as

3) Linearized Equation of Motion for the Whole System:
In further derivations of equations of motion (101) and (106), it should be taken into account that product of two dual terms d (·) d (·) is equal to zero due to ε 2 = 0. This property again leads to cancellation of all higher order terms associated with the elastic displacements and velocities.
The linearized vector of generalized forces is formulated as the following sum: (109) where the wrenches of constraint are cancelled from the equation of motion due to D'Alembert's principle. The linearized equation of motion for the entire system is formulated as Substitution of (91) and (96) into (110) leads to the following matrix formulation of the equation of motion: where the mass matrix M, the centrifugal and Coriolis force matrix C, and the stiffness matrix K are defined as where Θ j,n = ∂(θ j,n )/∂q. The proposed modeling procedure is summed up in Algorithm 1. The goal of this section is to give a reader more details on how the proposed method can be implemented. The discussion on the best practices how the assumed modes and mass distribution have to be selected is outside the scope of this work.

A. Modeling
Consider a robotic arm system given in Fig. 4. The open-chain system consists of 4 links. Link 1 is modeled as a rigid body, Link 2 and Link 3 are modeled as elastic bodies, and Links 4 is modeled as a rigid body. The joints between the links are assumed to be rigid. Each Link i has a body-fixed Frame i attached in the controlled end of the link. Elastic links additionally have a body-fixed Frame i attached at the distal end of the link and a body-fixed Frame i e attached at the point located at a distance x from i along x i axis. The elasticity is modeled by the Euler-Bernoulli beam model and the AMM, where elastic axial elongation is neglected. Two modes are considered for each type of vibrations, where the modes are derived using the clamped-mass boundary conditions [7] for flexural modes and clamped-free for torsional modes.
The vector of generalized coordinates (86) is given as where q i is a rigid-body rotation in Joint i, q bx i for i = 1, 2 are the generalized coordinates for torsional oscillations, q bx i for Consider that the dual rotation matrix can be alternatively represented asR = (R 1 , R 2 ), then the coordinate transformations between the frames are defined by the following dual rotation matrices. The dual rotation matrix from Frame 0 to Frame 1 isR 0 1 = (R z (q 1 ), 0), where q 1 is the rotation angle about the z 1 axis. The dual rotation matrix from Frame 1 to Frame 2 isR 1 2 = (R x (−π/2)R z (−π/2)R z (q 2 ), 0), where q 2 is the rotation angle about the z 2 axis. The dual rotation matrix from Frame 2 to Frame 2 isR 2 2 = (I, θ 2 × 2 ), where θ 2 2 = [θ 1,2 θ 2,2 θ 3,2 ] T is a vector of rotation angles due to elastic deformations about x 2 , y 2 and z 2 axes. The dual rotation matrix from Frame 2 to Frame 3 isR 2 3 = (R z (q 3 ), 0), where q 3 is the rotation angle about the z 3 axis. The dual rotation matrix from Frame 3 to Frame 3 isR The dual screw matrix for elastic motion of Frame j can be expressed in the coordinates of Frame i and referenced to the COG of Link i for rigid links, or to the origin of Frame i e for elastic links by the dual screw transformationŝ while the projecting form is defined as The Jacobian form of the dual projection matrix for the rigid Link 1 is given asP The Jacobian form of the dual projection matrix for the elastic Link 2 is given aŝ The Jacobian form of the dual projection matrix for the elastic Link 3 is given aŝ where g = 9.81 m/s 2 is the acceleration of gravity. The dual wrenches of the input torques of Link 1 expressed in the coordinates of 1 and referenced to the COG of Link 1 arê formulated by (109). The rest of the modeling procedure is given in Section V.

B. Simulation
In this section, the simulation of the proposed manipulator is carried out. To highlight the advantage of the linearized model, we have implemented the same model without any linearization and compared the simulation results and the simulation times. In addition, to demonstrate the effect of linearization for the symbolic expressions, we have presented the explicit nonlinear expression for the velocities v 0,2 e and v 0,3 e and their linearized versions in Appendix B. These specific velocities were selected as an example.
The geometrical and mechanical constants of the manipulator which were used in the simulation are given in Table I. The numerical simulation was performed in the following order, see Fig. 5. The manipulator was initiated at (q 1 , q 2 , q 3 , q 4 ) = (0, π/6, 2π/3, π/6) rad, while q i = 0 for i = 5 . . . 16 andq i = 0 for i = 1 . . . 16. The implemented PD (proportional derivative)  controller was commanded the desired position of q i for i = 1 . . . 4 equal to the initial value. The control of flexible manipulators is not the aim of this work, therefore the details of the controller are omitted. Then, at time t = 0.5 s, the desired position of q i for i = 2, 3, 4 was set to (q 2 , q 3 , q 4 ) = (π/3, π/3, π/3) rad and the manipulator started the manoeuvre. In addition, at time t = 3 s, the desired position of q 1 was set to q 1 = π/4 rad. The simulation was stopped at time t = 6 s. The values of generalized coordinates and generalized speeds throughout the simulation are given in Fig. 6, where the units are converted to deg and deg/s for convenience.
The manipulator was initiated in the underformed position, and after time t ≥ 0 some osculations were initiated in local y i direction of Link i due to the force of gravity, see Fig. 7. After time t ≥ 0.5 s, the manipulator started moving in Joints 2, 3, and 4, which can be seen in Fig. 6. Since the slew motion was absent, the manoeuvre was performed in xy planes of Links 2, 3, and 4. Such motion additionally excited flexural oscillations in  the local y i direction of Links 2 and 3, see Fig. 7, while torsional osculations, as well as flexural oscillations in z i direction of Link i were equal to zero, see Figs. 8 and 9. After time t ≥ 3 s, the slew motion started, which excited both torsional oscillations of Link 2 and 3, as well as flexural oscillations in the local z i direction of the same links, see Figs. 8 and 9. The control input torques throughout the simulation are given in Fig. 10.
The results presented in Figs. 6-10 have shown that the difference in simulation results between the nonlinear and linearized models is small, while the linearized model is more computationally efficient. For the comparison of computational efficiency the CPU time of simulations was measured. Simulation of the nonlinear model took 3901 s of the CPU time in Matlab/Simulink, while the same implementation of the linearized model took 314 s to simulate. In addition, for the case of the linearized model optimization of the code for symbolic functions using matlabFunction allowed for the more computationally efficient implementation which took 14 s of the CPU time to simulate. An attempt of equivalent code optimization of the symbolic functions for the nonlinear model was taking impractically long time and was not possible to finalize with the available personal computer.

VII. CONCLUSION
In this article, we presented a systematic and general method for deriving kinematics and dynamics of flexible manipulators, where the final equations are linearized in all elastic displacements and velocities. To form the mathematical background of the method we introduced dual rotation and dual homogeneous transformation matrices, which were presented as matrix Lie groups. The velocities and angular velocities of links were given in terms of dual twists, which are the vector representation of Lie algebra associated with the group of dual homogeneous transformation matrices. Dual twists were conveniently transformed using dual screw transformations, which led to a systematic procedure. The dynamical formulation was based on Kane's equation of motion, where partial velocities and partial angular velocities were given as dual screws, which were arranged as columns of dual projection matrices. The proposed method can also be used for open-chain multibody systems with both rigid and elastic bodies, including the cases when the rotational inertia of cross-sections in elastic links has to be accounted for. This was demonstrated in the implementation example, where a 4-DOF manipulator with two rigid and two elastic links was modeled. The method provided significant simplification in model derivation process compared to the original formulation, which allows for modeling of more complex systems.
The application of the method is limited to open-chain manipulators, where linear elastic kinematics is used and geometric stiffening effects are not taken into account. In the presented formulation, elasticity can only be modeled by the AMM, however modifications of the method for implementing FEM can be seen as a part of the future work.
APPENDIX A EXPONENTIAL AND LOGARITHMIC MAPS FOR G 1 AND G 2 LIE GROUPS Recall an element of the Lie group G 1 given by (8) as The logarithm [Ω] ∧ = log(X) ∈ g 1 is written in the form where U and W are skew-symmetric matrices, while u = [U] ∨ and w = [W] ∨ are the corresponding vector representations. The expressions for U and W will be derived in the following based on the derivation of the time derivative of the exponential map in [16]. It can be verified by direct computation that the kth power of (128) is given as Let the element (2,1) of the matrix (129) be denoted by Z k where L U W = UW and R U W = WU are the left and right translation maps. From ad(U)W = UW − WU it follows that: it follows that (130) can be written as where Ψ l (ad(u)) is the left Jacobian in SO(3) [33] and [38]. The exponential of (128) is then given as .
It can be concluded that the logarithm of (8) is given by where ad(φ) = φ × and φ × = log(A), and the inverse of the left Jacobian in SO(3) is given in [33]. The alternative derivation of the closed-form solutions for exponential and logarithmic maps (136) and (137) for the case of adjoint representation of SE(3) is given in [34].
Recall an element of the Lie group G 2 given by (14) as The logarithm in G 2 is given by It is seen that the kth power of (139) is given as The exponential representation is then found from the series expansion to be is the left Jacobian of G 1 . The closed-form solution for (142) is found in the same way as the closed-form solution for the left Jacobian in SE (3) [39] where

APPENDIX B LINEARIZATION OF VELOCITY
In this Appendix, the velocity of the origin of Frame 3 e is derived using coordinate-free vectors and then linearized as given in Kane's method [15]. This velocity is selected as an example to demonstrate the reduction of symbolic expressions due to linearization. It is noted, that the reduction of symbolic expressions for the final equations is even more significant. The used notations are given in Section VI.
The linear velocity of the COG of Link 1 relative to the inertial frame is v 0,m 1 = 0.
The velocity of the origin of Frame 2 relative to the inertial frame is v 0,2 = 0. The distance from the origin of Frame 2 to the origin of Frame 2 e is p 2,2 e = x x 2 + r 2 y 2 + r 3 z 2 , where r 2 (x, t) and r 3 (x, t) are the elastic deformations along the y 2 and z 2 axes. Then the velocity of the origin of Frame 2 e relative to the inertial frame is v 0,2 e = v 0,2 + 2 d p 2,2 e dt + w 02 × p 2,2 e =ṙ 2 y 2 +ṙ 3 z 2 − u 1 r 3 x 1 + u 1 (xs 2 + r 2 c 2 ) z 2 where w 02 = u 1 z 1 + u 2 z 2 and c i = cos(q i ), s i = sin(q i ). The velocity of the origin of Frame 2 relative to the inertial frame is found similarly to (146) as v 0,2 =ṙ l2 y 2 +ṙ l3 z 2 − u 1 r l3 x 1 + u 1 (L 2 s 2 + r l2 c 2 ) z 2 where p 2,2 = L 2 x 2 + r l2 y 2 + r l3 z 2 and r li = r i (L 2 , t) are the elastic deformations along the y 2 and z 2 axes at the distal end of Link 2. The distance from the origin of Frame 2 to the origin of Frame 3 is p 2 ,3 = 0, then the velocity of the origin of Frame 3 relative to the inertial frame is v 0,3 = v 0,2 . The linearization of (147) gives ṽ 0,3 = v 0,3 , i.e., in this case, the linear and nonlinear velocities are equal. The rotation from Frame 2 to Frame 2 e due to elastic deformation is described by three angles θ i , for i = 1, 2, 3 about the x 2 , y 2 , and z 2 axes, respectively.