# Intelligent Manipulator with Flexible Link and Joint: Modeling and Vibration Control

Authors: Minqiang Shao; Yuming Huang; Vadim V. Silberschmidt
Journal: Shock and Vibration
 (2020)

---

## Abstract

This paper presents a finite-element (FE) model of a manipulator with a flexible link and flexible joint as well as embedded PZT actuators and proposes a corrected rebuilt reduced model (CRRM) to make its dynamic characteristics more consistent with reality and facilitate control design. The CRRM considers the holding torque of the manipulator driving motor and eliminates the response divergence induced by a fault of the mass matrix of the FE model. In order to reduce the dimensions and maintain the precision of the model, an iterated improved reduction system (IIRS) method is adopted. Additionally, a LQR controller is designed based on the output function of the improved model. The simulation results demonstrate that the CRRM is consistent with reality and the active controller has good performance in suppressing vibration of the manipulator with both the flexible link and the flexible joint.


---

## Body

## 1. Introduction

Lightweight and large-scale flexible structures are utilized in many fields. Some of them exhibit excellent dynamic performance when integrated with intelligent components. Especially in space applications, manipulators are always designed as lightweight structures [1]. Such structures have the advantages of lower energy consumption in operation and higher speed in motion. However, they always have small damping and low natural frequencies that can lead to residual vibration with large amplitude [2]. Additionally, the vibration and rigid motion of manipulators are highly coupled [3]. So, it is a challenge to control the vibration and maintain the rigid motion of the system at the same time. Dynamics of manipulators was studied for a long time, with first works traced back to the 1970s. Based on the differences between modeling, control and experimental studies, Dwivedy and Eberhard [4] summarized the original works in the area of dynamics of flexible robotic manipulators performance. The main focus of scholars in the early studies was on dynamic modeling. Currently, dynamics and control of flexible manipulators is still a hot topic; in addition to development of new models, more attention is paid to the control of flexible manipulators, with applications to different robotic platforms [5, 6]. Control of vibration of flexible manipulators is also concerned, with active control methods and advanced intelligent materials used in some research studies [7–10].A high-precision mathematical model is crucial for design of a controller to avoid unwanted vibration. The two main classes of modeling approaches are an assumed mode method (AMM) and finite-element (FE) method. Junkins and Kim [11] described the theory of these two methods systematically. Shin and Choi [12] developed a model of a two-link flexible manipulator using the Lagrangian method associated with the AMM. Abe [13] proposed a method for prediction of a trajectory of two-link rigid-flexible manipulators to suppress residual vibration, also employing a consideration of the Lagrangian approach and AMM. The order of the assumed mode models is low, and it is easy to apply them to the multilink manipulators [14, 15]; the AMM is also suitable for structures with available mode shapes. However, it is difficult to assess the real dynamic characteristics of a structure with complex geometry or constraints using the AMM. On the other hand, the FE method can deal with complex structures and has a wider of applications. Orszulik and Shan [16, 17] deduced the dynamic functions of a single-link manipulator based on the FE method in detail. Dubay et al. [2] presented an approach for active vibration control of a single-link manipulator based on the FE model. Generally, single-link manipulators were studied widely with the FE method [18], and the dynamic characteristics of the models were verified by experiments [19]. Furthermore, manipulators with two or more links were also analyzed employing the FE method [20, 21]. An important purpose of dynamical analysis of flexible manipulators is to suppress the vibration of flexible structures [22, 23]. Many studies showed that the control schemes were effective for vibration control of flexible manipulators. For example, a linear quadratic Gaussian (LQG) method was applied to control a flexible-beam structure subjected to external disturbance [24], an optimal tracking control theory was used for position control of a flexible hub-beam system [25], and an optimal placement of sensors and actuators distributed on flexible space structures [26]. Other control algorithms, such as robust control [21, 27], sliding-surface constraint scheme, velocity feedback control, and fuzzy control [28–30], were also broadly used in controller design for manipulators.Almost all lightweight and large-scale manipulators have flexible arms, and joints of real-life manipulators are always flexible. Most of the research studies in previous papers only refer to the flexible link of manipulators [31, 32]. There are a few studies considering both flexible links and joints within the same approach. In order to develop an effective model adequately reflecting the reality, it is necessary to consider both flexible links and flexible joints in the same model. This paper considers this case and studies the vibration suppression of a manipulator with both flexible links and joints based on the FE method. A pair of embedded-PZT actuators is also considered in the FE model. The IIRS method [33] is introduced to reduce the order of the FE model. The high-precision character of the suggested method helps the reduced model to maintain consistent performance in the effective frequency range comparable to that of the original. Furthermore, in order to fit the dynamic characteristics of the model to reality, excess inertia instead of a holding torque is added to the rotor of the drive motor. A type of the mass matrix of the developed FE model induces its response to diverge. So, a reconstructed method is presented and deals with this problem effectively. The function of the reconstructed model is transformed into an expression with output variables; then, it is easy to design an optimal controller based on the output function. The obtained simulation results show good performance of the suggested approach in vibration suppression of manipulators.


## 2. Dynamic Modeling

A schematic of a flexible manipulator with embedded actuators is shown in Figure1. This manipulator has general dynamical characteristics of most realistic systems with both a flexible link and a flexible joint. The manipulator is primarily composed of an electric motor, an elastic joint, a follow-up unit, PZT actuators, a mounting base for the actuators, a flexible link, and a tip mass. The manipulator is driven by the electric motor, and the follow-up unit is connected to the rotor of the motor by the elastic joints. The link is a flexible beam. A pair of actuators was embedded in the mounting base, which was fixed to the root segment of the beam.Figure 1
Manipulator with flexible link and flexible joint.In some cases, the mounting base can be considered as a part of the flexible beam. At the end of the link, a tip mass was considered. This manipulator demonstrates a general case for conditions with one link. The outline of the constituent parts of the manipulator can be designed arbitrarily. Therefore, the model could describe a vast majority of structures with slender rotating beams. The motion of the manipulator was planar. The typical Euler–Bernoulli assumption was made for the flexible beam.


### 2.1. Lagrange Function of the System

Kinetic and potential energies of the entire system was generated by the rotor of the electrical motor, the flexible joint, the follow-up unit, the flexible beam, and the tip mass. The kinetic energies of each part were denotedTr, Tf, Tb, and Tm while the potential energies, Vr, Vf, Vb, and Vm, respectively. The mounting base and piezoelectric actuators were considered as constituent parts of the beam. Therefore, the kinetic energy function is given by(1)T=Tr+Tf+Tb+Tm,where Tr=1/2Jrβ˙2, Tf=0, Tb=∑i=1NTbi, Tm=1/2mmvm2+1/2Jmωm2, Jr and β are the rotational inertia and the rotation angle of the rotor, respectively, Tbi is energy of the ith finite element of the beam, N is the total number of elements of the beam, mm and Jm are the mass and rotational inertia of the tip mass, respectively, and vm and ωm are the translational and angular velocities of the tip mass, respectively.Similarly, the potential energy can be described as(2)V=Vr+Vf+Vb+Vm,where Vr=0, Vf=1/2kfθ−β2, Vb=∑i=1NVbi, Vm=0, kf is the coefficient of elasticity of the flexible joint, θ is the rotation angle of the mounting base, and Vbi is the energy of the ith finite element of the beam.Defining(3)L1=Lr=Tr−Vr,L2=Lf=Tf−Vf,L3=Lb=bb−Vb,L4=Lm=Tm−Vm.The Lagrangian equation of the system can be described as(4)∑j=14ddt∂Lj∂q˙k−∂Lj∂qk=Qk,k=1,…,Na.The abovementioned equation indicates that the model of the system can be established by combining all constituent parts of the system. Nevertheless, the responses of movement and vibration of the system are coupled to each other, and the total degrees of freedom depend on FEs of the beam.


### 2.2. FE Model of Flexible Structure with Rotation

The finite-element method (FEM) was applied for flexible components. The beam with the mounting base was divided into finite segments. The material of the beam with densityρb was assumed to be homogeneous and isotropic. The cross-section area Abi of the ith element was approximately constant. Therefore, the kinetic energy of the element is given by(5)Tbi=∫0li12ρbAbivi2x,tdx,where vix,t=a+xi+xθ˙+w˙ix,t, x∈0,li, wi is the transverse displacement of ith element, li is the length of the ith element, a is the distance between the rotational axis and the first node of the beam, and xi is the distance between the first and the ith nodes of the beam. The potential energy is defined as(6)Vbi=∫0liEbIbi∂2wix,t∂x2dx,where Eb is Young’s modulus for the beam and Ibi is the inertia of the element of the beam.By combining equations (5) and (6) and Lagrangian equation, the coefficient matrices and equation of the ith element can be derived as(7)mθθiMθqiMθqiTMqqiθ¨q¨i+000Kqqiθqi=QθiQqi,where qi is the displacement vector of ith element,(8)mθθi=ρbAbia+xi2+a+xili+li23li,Mθqi=ρbAbia+xili126,li,6,−li+ρbAbili2609,2li,14,−3li,Mqqi=ρbAbili420156−22li5413li−22li4li2−13li−3li254−13li15622li13li−3li222li4li2,Kqqi=EbIbili3126li−126li6li4li2−6li2li2−12−6li12−6li6li2li2−6li4li2.According to the positions of nodes of the elements, the total mass and stiffness matrices were assembled.w1j and w2j were defined as the displacement and rotation angle of the jth node of the beam FE, and Q¯1j and Q¯2j are the external force and the moment acting on the jth node, respectively. Therefore, the dynamic equation of the finite-element model of the beam is given by(9)mθθMθqMθqTMqqθ¨q¨+000Kqqθq=QθQ¯q,where mθθ, θ, and Qθ are scalar quantities, q and Qq are 2(N + 1) dimensional vectors, q=w11,w21,w12,w22,…,w1N+1,w2N+1T, Q¯q=Q¯11,Q¯21,Q¯12,Q¯22,…,Q¯1N+1,Q¯2N+1T, matrix Mθq∈R1×2N+1, Mqq,Kqq∈R2N+1×2N+1, and N is the total number of elements of the FE model of the beam. The matrices Mqqi and Kqqi shown in equation (7) are rewritten in a partitioned form:(10)Sqqi=S11iS12iS12iTS22i,Sqqi=Mqqi,Kqqi,where S11i,S12i,S22i∈R2×2. By combining equations (7) and (10), the elements of the coefficient matrices in equation (9) can be expressed by(11)mθθ=∑i=1NρbAbia+xi−12li,Mθq1=Mθq11,Mθq2=Mθq12,Mθq2i+1=Mθqi3+Mθqi+11,Mθq2i+1=Mθqi4+Mθqi+12,Mθq2N+1=MθqN3,MθqN2N+1=MθqN4,i=1,2,…,N−1,Sqq=S111S121S221+S112⋱⋱S12i−1S22i−1+S11i⋱sym⋱S12N−1S22N−1+S11NS12NS22N.


### 2.3. Model of Piezoelectric Actuators

Two piezoelectric actuators shown in Figure1 were embedded in the mounting base, which was considered as a part of the beam. The actuators, denoted ACT1 and ACT2, were placed near the flank and symmetrical as to the central plane of the mounting base. Forces generated by the actuators are considered as external forces. Deformations of the flexible structure induce interactions between the beam and the actuators. Strain along the axis of the actuator was generated by structural deformation of the beam and active deformation induced by external voltage. The active deformation Δlpv of the actuator without constraint is proportional to voltage, and it can be expressed as(12)Δlpv=kvu,where u is the voltage and kv is the proportionality coefficient.A pretightening load should be applied to the actuators, since they cannot take a tension force. Assuming the cross-section of the actuator to be constant, the deformation induced by the pretightening load can be expressed as(13)Δlpt=−FtlpEpAp,where Ft is the pretightening force and lp, Ap, and Ep are the geometrical and material parameters of the PZT actuator.Another important factor affecting deformation of the actuators was induced by a nodal rotational angle of the FE beam. According to the position of the actuators, deformation of ACT1 induced by nodal rotation was expressed by(14)Δlpr=−w2p2−w2p1h,where p1 and p2 are nodal numbers of the finite element of the beam, on which the actuators act, and h is the distance between the acting point of ACT1 to the central plane of the beam.Therefore, the total deformation of ACT1 from its free state without constraint to the state with pretightening force and interaction of the beam is given by(15)Δlp1=Δlpt+Δlpr−Δlpv=−FtlpEpAp−w2p2−w2p1h−kvu.ACT2 was placed on the other side of the beam, symmetric to ACT1. So, the direction of deformation induced by rotation is opposite to that of ACT1. The driven voltage also was set opposite to that of ACT1. Therefore, the total deformation of ACT2 is expressed by(16)Δlp2=−FtlpEpAp+w2p2−w2p1h+kvu.Deformation of the actuators generated forces acting on the beam; according to Hooke’s law, the forces can be expressed as(17)Fpi=EpεiAp=EpApΔlpilp,i=1,2,where εi is the strain of the ith actuator. Substituting equations (15) and (16) into equation (17), the moments acting on the nodes p1 and p2 can be obtained. By defining an anticlockwise direction as positive, one obtains(18)Mp2=−Mp1=Fp2h−Fp1h=c¯pw2p2−w2p1+k¯vu,where c¯p=2EpAph2/lp and k¯v=2EpAphkv/lp.Equation (18) was substituted into dynamic equation (9), and the items including w2p1 and w2p2 were moved to the left-hand part of the equation. The updated elements of the stiffness matrix of the beam are expressed by(19)K¯qqij=Kqqij+c¯p,i=j=2p1or2p2,Kqqij−c¯p,i=2p1,j=2p2ori=2p2,j=2p1,Kqqij,else,and the force vector of the beam can be derived as(20)Q¯qi=−k¯vu,i=2p1,k¯vu,i=2p2,0,else,where i=1,2,…,N+1.


### 2.4. System Assemble

By combining equations (1)–(3) and equation (4), the dynamic equation of the rotor and the flexible joint can be derived as(21)Jr000β¨θ¨+kJ−kJ−kJkJβθ=Mr0,where Mr is the torque acting on the rotor and kJ is the elastic rotation coefficient of the flexible joint.Assuming the tip mass set on the end point, the numbered (N + 1)th node of the beam, the velocity of the tip mass is(22)vm=a+lθ˙+w˙1N+1,where l is the length of the beam and w1N+1 is the deflection of the (N + 1)th node of the beam.The angular velocity of the tip mass is expressed as(23)ωm=w˙2N+1Substituting equations (24)-(25) into equation (1) and combining equations (3) and (4) obtains mass and stiffness matrices:(24)Mm=mma+l2mma+l0mma+lmm000Jm,Km=03×3,where Mm and Km correspond to vector θ¨,w¨1N+1,w¨1N+1T and θ,w1N+1,w2N+1T, respectively.Combining equations (9), (19)–(21), and (24), assembling the mass and stiffness matrices, the entire model of the system is derived as(25)Jr000m¯θθM¯θq0M¯θqTM¯qqβ¨θ¨q¨+kJ−kJ0−kJkJ000K¯qqβθq=Mr0Q¯q,where K¯qq and Q¯q are clarified in equations (19) and (20), m¯θθ=mθθ+mma+l2,(26)M¯θqj=Mθqj+mma+l,j=2N+1,Mθqj,else,M¯qqij=Mqqij+mm,i=j=2N+1,Mqqij+Jm,i=j=2N+2,Mqqij,else.


## 2.1. Lagrange Function of the System

Kinetic and potential energies of the entire system was generated by the rotor of the electrical motor, the flexible joint, the follow-up unit, the flexible beam, and the tip mass. The kinetic energies of each part were denotedTr, Tf, Tb, and Tm while the potential energies, Vr, Vf, Vb, and Vm, respectively. The mounting base and piezoelectric actuators were considered as constituent parts of the beam. Therefore, the kinetic energy function is given by(1)T=Tr+Tf+Tb+Tm,where Tr=1/2Jrβ˙2, Tf=0, Tb=∑i=1NTbi, Tm=1/2mmvm2+1/2Jmωm2, Jr and β are the rotational inertia and the rotation angle of the rotor, respectively, Tbi is energy of the ith finite element of the beam, N is the total number of elements of the beam, mm and Jm are the mass and rotational inertia of the tip mass, respectively, and vm and ωm are the translational and angular velocities of the tip mass, respectively.Similarly, the potential energy can be described as(2)V=Vr+Vf+Vb+Vm,where Vr=0, Vf=1/2kfθ−β2, Vb=∑i=1NVbi, Vm=0, kf is the coefficient of elasticity of the flexible joint, θ is the rotation angle of the mounting base, and Vbi is the energy of the ith finite element of the beam.Defining(3)L1=Lr=Tr−Vr,L2=Lf=Tf−Vf,L3=Lb=bb−Vb,L4=Lm=Tm−Vm.The Lagrangian equation of the system can be described as(4)∑j=14ddt∂Lj∂q˙k−∂Lj∂qk=Qk,k=1,…,Na.The abovementioned equation indicates that the model of the system can be established by combining all constituent parts of the system. Nevertheless, the responses of movement and vibration of the system are coupled to each other, and the total degrees of freedom depend on FEs of the beam.


## 2.2. FE Model of Flexible Structure with Rotation

The finite-element method (FEM) was applied for flexible components. The beam with the mounting base was divided into finite segments. The material of the beam with densityρb was assumed to be homogeneous and isotropic. The cross-section area Abi of the ith element was approximately constant. Therefore, the kinetic energy of the element is given by(5)Tbi=∫0li12ρbAbivi2x,tdx,where vix,t=a+xi+xθ˙+w˙ix,t, x∈0,li, wi is the transverse displacement of ith element, li is the length of the ith element, a is the distance between the rotational axis and the first node of the beam, and xi is the distance between the first and the ith nodes of the beam. The potential energy is defined as(6)Vbi=∫0liEbIbi∂2wix,t∂x2dx,where Eb is Young’s modulus for the beam and Ibi is the inertia of the element of the beam.By combining equations (5) and (6) and Lagrangian equation, the coefficient matrices and equation of the ith element can be derived as(7)mθθiMθqiMθqiTMqqiθ¨q¨i+000Kqqiθqi=QθiQqi,where qi is the displacement vector of ith element,(8)mθθi=ρbAbia+xi2+a+xili+li23li,Mθqi=ρbAbia+xili126,li,6,−li+ρbAbili2609,2li,14,−3li,Mqqi=ρbAbili420156−22li5413li−22li4li2−13li−3li254−13li15622li13li−3li222li4li2,Kqqi=EbIbili3126li−126li6li4li2−6li2li2−12−6li12−6li6li2li2−6li4li2.According to the positions of nodes of the elements, the total mass and stiffness matrices were assembled.w1j and w2j were defined as the displacement and rotation angle of the jth node of the beam FE, and Q¯1j and Q¯2j are the external force and the moment acting on the jth node, respectively. Therefore, the dynamic equation of the finite-element model of the beam is given by(9)mθθMθqMθqTMqqθ¨q¨+000Kqqθq=QθQ¯q,where mθθ, θ, and Qθ are scalar quantities, q and Qq are 2(N + 1) dimensional vectors, q=w11,w21,w12,w22,…,w1N+1,w2N+1T, Q¯q=Q¯11,Q¯21,Q¯12,Q¯22,…,Q¯1N+1,Q¯2N+1T, matrix Mθq∈R1×2N+1, Mqq,Kqq∈R2N+1×2N+1, and N is the total number of elements of the FE model of the beam. The matrices Mqqi and Kqqi shown in equation (7) are rewritten in a partitioned form:(10)Sqqi=S11iS12iS12iTS22i,Sqqi=Mqqi,Kqqi,where S11i,S12i,S22i∈R2×2. By combining equations (7) and (10), the elements of the coefficient matrices in equation (9) can be expressed by(11)mθθ=∑i=1NρbAbia+xi−12li,Mθq1=Mθq11,Mθq2=Mθq12,Mθq2i+1=Mθqi3+Mθqi+11,Mθq2i+1=Mθqi4+Mθqi+12,Mθq2N+1=MθqN3,MθqN2N+1=MθqN4,i=1,2,…,N−1,Sqq=S111S121S221+S112⋱⋱S12i−1S22i−1+S11i⋱sym⋱S12N−1S22N−1+S11NS12NS22N.


## 2.3. Model of Piezoelectric Actuators

Two piezoelectric actuators shown in Figure1 were embedded in the mounting base, which was considered as a part of the beam. The actuators, denoted ACT1 and ACT2, were placed near the flank and symmetrical as to the central plane of the mounting base. Forces generated by the actuators are considered as external forces. Deformations of the flexible structure induce interactions between the beam and the actuators. Strain along the axis of the actuator was generated by structural deformation of the beam and active deformation induced by external voltage. The active deformation Δlpv of the actuator without constraint is proportional to voltage, and it can be expressed as(12)Δlpv=kvu,where u is the voltage and kv is the proportionality coefficient.A pretightening load should be applied to the actuators, since they cannot take a tension force. Assuming the cross-section of the actuator to be constant, the deformation induced by the pretightening load can be expressed as(13)Δlpt=−FtlpEpAp,where Ft is the pretightening force and lp, Ap, and Ep are the geometrical and material parameters of the PZT actuator.Another important factor affecting deformation of the actuators was induced by a nodal rotational angle of the FE beam. According to the position of the actuators, deformation of ACT1 induced by nodal rotation was expressed by(14)Δlpr=−w2p2−w2p1h,where p1 and p2 are nodal numbers of the finite element of the beam, on which the actuators act, and h is the distance between the acting point of ACT1 to the central plane of the beam.Therefore, the total deformation of ACT1 from its free state without constraint to the state with pretightening force and interaction of the beam is given by(15)Δlp1=Δlpt+Δlpr−Δlpv=−FtlpEpAp−w2p2−w2p1h−kvu.ACT2 was placed on the other side of the beam, symmetric to ACT1. So, the direction of deformation induced by rotation is opposite to that of ACT1. The driven voltage also was set opposite to that of ACT1. Therefore, the total deformation of ACT2 is expressed by(16)Δlp2=−FtlpEpAp+w2p2−w2p1h+kvu.Deformation of the actuators generated forces acting on the beam; according to Hooke’s law, the forces can be expressed as(17)Fpi=EpεiAp=EpApΔlpilp,i=1,2,where εi is the strain of the ith actuator. Substituting equations (15) and (16) into equation (17), the moments acting on the nodes p1 and p2 can be obtained. By defining an anticlockwise direction as positive, one obtains(18)Mp2=−Mp1=Fp2h−Fp1h=c¯pw2p2−w2p1+k¯vu,where c¯p=2EpAph2/lp and k¯v=2EpAphkv/lp.Equation (18) was substituted into dynamic equation (9), and the items including w2p1 and w2p2 were moved to the left-hand part of the equation. The updated elements of the stiffness matrix of the beam are expressed by(19)K¯qqij=Kqqij+c¯p,i=j=2p1or2p2,Kqqij−c¯p,i=2p1,j=2p2ori=2p2,j=2p1,Kqqij,else,and the force vector of the beam can be derived as(20)Q¯qi=−k¯vu,i=2p1,k¯vu,i=2p2,0,else,where i=1,2,…,N+1.


## 2.4. System Assemble

By combining equations (1)–(3) and equation (4), the dynamic equation of the rotor and the flexible joint can be derived as(21)Jr000β¨θ¨+kJ−kJ−kJkJβθ=Mr0,where Mr is the torque acting on the rotor and kJ is the elastic rotation coefficient of the flexible joint.Assuming the tip mass set on the end point, the numbered (N + 1)th node of the beam, the velocity of the tip mass is(22)vm=a+lθ˙+w˙1N+1,where l is the length of the beam and w1N+1 is the deflection of the (N + 1)th node of the beam.The angular velocity of the tip mass is expressed as(23)ωm=w˙2N+1Substituting equations (24)-(25) into equation (1) and combining equations (3) and (4) obtains mass and stiffness matrices:(24)Mm=mma+l2mma+l0mma+lmm000Jm,Km=03×3,where Mm and Km correspond to vector θ¨,w¨1N+1,w¨1N+1T and θ,w1N+1,w2N+1T, respectively.Combining equations (9), (19)–(21), and (24), assembling the mass and stiffness matrices, the entire model of the system is derived as(25)Jr000m¯θθM¯θq0M¯θqTM¯qqβ¨θ¨q¨+kJ−kJ0−kJkJ000K¯qqβθq=Mr0Q¯q,where K¯qq and Q¯q are clarified in equations (19) and (20), m¯θθ=mθθ+mma+l2,(26)M¯θqj=Mθqj+mma+l,j=2N+1,Mθqj,else,M¯qqij=Mqqij+mm,i=j=2N+1,Mqqij+Jm,i=j=2N+2,Mqqij,else.


## 3. Model Reduction and Dynamical Response

Based on the parameters of the manipulator system presented in Table1, a full FE model, in which flexible components were described into 100 elements, was developed. The follow-up unit and the mounting base of the actuators were divided into 5 elements based on their cross-sections. The first five elements were special; they are marked e1 to e5 in Figure 2. Parameters of the special elements and PZT actuators are shown in Tables 1 and 2. The manipulator was driven by the electrical motor, with its rotor rotating together with the follow-up unit. The rotor was locked by a holding torque [34] after energizing the motor. The holding torque was always present, no matter the rotor was in motion or at rest. Therefore, the rotor remained in a fixed rotation constraint state at any instant. Another constraint was the deflection and rotation of node 1, w11=0 and w12=0. According to the constraints mentioned, the columns and rows corresponding to variable β, w11 and w12 of coefficient matrices K, and M should be eliminated.Table 1
Parameters of the manipulator and special elements.
Following up unit
Mounting base
Flexible beam
ρ
f kg/m3
E
f (GPa)
L
f (mm)
ρ
m (kg/m3)
E
m (GPa)
L
m (mm)
ρ
b (kg/m3)
E
b (GPa)
L
b (mm)
b
b (mm)
h
b (mm)
7.9 × 103
200
50
7.9 × 103
200
82
7.9 × 103
200
500
20
1
e1
e2, e5
e3
e4
Rotor
Flex joint
Tip mass
A
e1 (mm2)
I
e1 (mm4)
A
e2,5 (mm2)
I
e2,5 (mm4)
A
e3 (mm2)
I
e3 (mm4)
A
e4 (mm2)
I
e4 (mm4)
J
r (Kgm2)
k
J (Nm/degree)
m
t (kg)
J
t (Kgm2)
525
5.3 × 104
2624
1.47 × 106
1563
7.1 × 105
384
4608
3 × 10−2
200
0.5
1.25 × 10−4Figure 2
Elements of flexible components.Table 2
Parameters of the piezoelectric actuator and acting position.
E
p (GPa)
A
p (mm2)
L
p (mm)
h
p (mm)
28.6
490
70
26In order to improve the computational efficiency, it is necessary to reduce the dimension of the system. The IIRS method [31] has good performance and is easy to work with. Combined with the full FE model, assembling positions of the actuators and the constraints the master coordinates were selected:(27)xm=θ,w32,w52,w31,w51,w61,w111,w211,…,w1011T.The detailed numerical results of frequencies and shapes of the modes of the reduced model are compared to those of the full model in Table3 and Figure 3. It is obvious that the natural frequencies of the reduced model converge with the full model within few iterations. The first two mode shapes show that the rotation angle θ is coupled with the flexible link. In one case, the mode shape of the flexible link is above the balance position of the link, and in the other case, it is below it. In both cases, the rotation angle of the flexible joint is positive (counterclockwise is considered positive). Comparing the higher modes of the link, the joint is more flexible, so the high-order mode shapes are similar to those of the cantilever beam. The calculated mode characteristics reproduce the real response and they agree with reality of the flexible-joint and flexible-link structures.Table 3
Mode frequencies of the reduced model compared to the full model.
Iteration no.
Mode frequencies (Hz)
1
2
3
4
5
6
7
8
9
10
1
1.588
1.987
13.468
43.931
91.882
157.376
240.498
341.385
460.403
600.949
5
—
—
—
—
—
—
—
341.369
460.121
596.996
40
—
—
—
—
—
—
—
—
460.119
596.919
282
—
—
—
—
—
—
—
—
—
596.913
Full model
1.588
1.987
13.467
43.931
91.882
157.376
240.498
341.369
460.119
596.913
Reduced error (%)
0
0
0.007
0
0
0
0
0
0
0Figure 3
First six mode shapes of the manipulator system.The mode results were calculated for the state of rotors fixing. In fact, the rotation angleβ changed constantly when the manipulator rotated. So, the reduced model should be rebuilt with equation (21). This results is(28)Jr00MIRSβ¨mx¨m+kJK1mK1mTKIRSβxm=MrFIRS,where K1m=−kJ,0,…,0, MIRS and KIRS are the reduced model matrices with rotation constraint, and FIRS is the force acting on the master nodes.By substituting parameters of Table1 into equation (28), the mode characteristics of the rebuilt reduced model (RRM) can be obtained; these results are presented in Table 4. The lowest frequency of the RRM in Table 4 is close to zero; it shows the rigid motion of the manipulator. Other frequencies of the RRM agree with those of the RM except for the first two modes. The RRM releases the rotational constraint; so, the first two modes are not correct since the RRM does not.Table 4
Natural frequencies of the rebuilt reduced model (RRM) compare the reduced model (RM).
Natural/Mode frequencies (Hz)
Rigid motion
1
2
3
4
5
6
7
8
9
RRM
0.0008
0.697
12.768
13.648
43.931
91.882
157.376
240.498
341.369
460.119
RM
—
1.588
1.987
13.468
43.931
91.882
157.376
240.498
341.369
460.119
Error (%)
—
56.1
542
0
0
0
0
0
0
0Consider the holding torque of the motor. How can rotation of the manipulator and the constraint induced by the holding torque be considered at the same time? It is seemingly a paradox. Another approach is needed to deal with it. The effect of the holding torque is to keep the rotor in the current state. As known, it is harder to move an object while a greater mass; it is also for rotation. A local correction, increasing the rotational inertia of the rotor, was implemented to ensure that the model had the same dynamic characteristics as the RM with the rotational constraint. The simulation results of the corrected RRM (CRRM), with inertia of the rotor changed toJrnew=100kg⋅m2, as shown in Table 5.Table 5
Natural frequencies of the corrected RRM (CRRM) compare reduced model (RM).
Natural/Mode frequencies (Hz)
Rigid motion
1
2
3
4
5
6
7
8
9
CRRM
2.7 × 10−5
1.595
1.980
13.648
43.931
91.882
157.376
240.498
341.369
460.119
RM
—
1.588
1.987
13.468
43.931
91.882
157.376
240.498
341.369
460.119
Error (%)
—
0.41
0.35
0
0
0
0
0
0
0The holding moment considered in CRRM is equivalent to the rotational inertia moment of the motor’s rotor, which simulates effectively the dynamical state of the flexible-joint and flexible-link system under the external moment. The proposed CRRM is different from other models of flexible links and joints mentioned in [3, 35], and to the author’s knowledge, similar models were not published. The frequencies of the first two modes of the CRRM agree well with those of the RM, while other modes still fit the RM. The proposed model effectively solves the problem of considering rotation and constraint at the same time.


## 4. Response of CRRM

### 4.1. Decoupled Method

Large-scale motivation and small-scale vibration, which are also called a slow system and a fast system, are simultaneously present in the response of the CRRM. Two types of mechanical behavior are coupled, so precision of the solution cannot be guaranteed. Traditional numerical algorithms are not suitable for this situation. It is well known that the modal method can decouple a vibration system into independent-mode functions. Therefore, the method can also deal with the rigid motion as an independent mode. By considering each independent-mode function separately, the slow and fast systems can be divided effectively.The dynamic function of the manipulator system expressed in equation (28) was rewritten in a general form:(29)Mx¨+Cx˙+Kx=f,where M,K∈Rn×n are the mass and stiffness matrices, C=αcM+βcK is the damping matrix, x=xt∈Rn×1 is the system’s displacement vector, and f=ft∈Rn×1 is the excitation vector of the system.Φ
=
φ
1
,
φ
2
,
…
,
φ
n is assumed to be the eigenvector matrix of equation (29). The displacement vector can be described in the independent-mode coordinates q=q1,q2,…,qnT as(30)xt=∑i=1Nφixqit,where φi is also the ith mode shape and N is the order of the system. Therefore, one obtains(31)miq¨i+ciq˙i+kiqi=fqi,i=1,…,N,where mi=φiTMφi, ki=φiTKφi, ci=αcmi+βcki, and fqi=φiTf.Equation (31) is a group of independent mode functions of the manipulator system.


### 4.2. Model Reconstitution

In general, it is easy to get the independent functions in equation (31) and substitute the solutions into equation (30); then, the solution of the system is obtained. But there is a problem that complicates attainment of the correct solution. Some eigenvalues in equation (29) are negative, causing the divergence of the manipulator system. This problem never appears in pure vibration systems, but in case, when the large-scale variables of the manipulator are considered in the function, the matrices cannot be well-defined. According to the decoupled function of equation (31), the mode mass mi is negative corresponding to respective eigenvalues. A simple method for dealing with the problem and maintaining the systems convergence is to reconstitute the model. The main process of the method was to make positive the negative mass of the mode. The new mass matrix can be expressed as Mqrec=diagmi, i=1,…,N.Therefore, the reconstituted mass matrix of the original system is(32)Mrec=Φ−TMqrecΦ−1.The reconstituted function was derived as(33)Mrecx¨+Crecx˙+Krecx=f,where Krec=K and Crec=αcMrec+βcKrec.The mode vectors of equation (33) are the same as those of the original system shown in equation (29) and whether the mode frequencies conform to reality can be validated experiment or experience.


### 4.3. Simulation Results

The manipulator was loaded with the external momentMr, generated by the driving motor. The damping coefficients were αc=0.1 and βc=10−4. Different types of the applied moment are listed in Table 6, and the response results are shown in Figures 4–7.Table 6
Different external excitations.
Forms of the torque
Mathematical expression
No. of response figures
Impulse excitation
M
r
=
100
,
t
=
0.5
s
0
,
else
Figure4
Step excitation
M
r
=
0
,
t
<
0.5
s
10
,
else
Figure5
Harmonic excitation
M
r
=
sin
2
π
⋅
1.6
⋅
t
Figure6
Hybrid excitation
M
r
=
sin
2
π
⋅
1.6
⋅
t
+
sin
2
π
⋅
2
⋅
t
+
w
t
,
where
w
t
is random excitation
Figure7Figure 4
Response with impulse excitation.
(a)
(b)Figure 5
Response with step excitation.
(a)
(b)Figure 6
Response with harmonic excitation.
(a)
(b)Figure 7
Response with hybrid excitation.
(a)
(b)The bandwidth of the impulse excitation was 2 ms; the responses in Figure4 include rigid rotation and flexible vibration. The damping coefficient of the manipulator system was small, so the rotating angle of the response could not reach a constant value, and it was also hard to attenuate the amplitude of vibration. Figure 5 exhibits the response of the manipulator for the step excitation.Vibration frequencies were similar to those of the impulse response, including the first two modes of the system. The step excitation caused the rigid rotation of the system to accelerate. At the same time, the vibration was changed to a DC (direct current) bias. The harmonic and hybrid excitations included the first two eigenfrequencies of the system. As shown in Figures6 and 7, the responses of the system with the flexible joint and flexible link diverge, when the frequencies of excitations are close to the first two natural frequencies of the system. This phenomenon is in accordance with characteristics of a small damping system. The system also demonstrated its ability to isolate high-frequency vibration induced by the motor; the response of the system to the hybrid excitation that included white noise illustrates this.The simulated responses of the system indicated that the vibration of the flexible beam could be induced easily by most types of excitations, and the response was hard to attenuate. The vibration regime almost included the first two eigenmodes of the system. According to the character of these responses of the system, it is necessary to design a vibration controller for it.


## 4.1. Decoupled Method

Large-scale motivation and small-scale vibration, which are also called a slow system and a fast system, are simultaneously present in the response of the CRRM. Two types of mechanical behavior are coupled, so precision of the solution cannot be guaranteed. Traditional numerical algorithms are not suitable for this situation. It is well known that the modal method can decouple a vibration system into independent-mode functions. Therefore, the method can also deal with the rigid motion as an independent mode. By considering each independent-mode function separately, the slow and fast systems can be divided effectively.The dynamic function of the manipulator system expressed in equation (28) was rewritten in a general form:(29)Mx¨+Cx˙+Kx=f,where M,K∈Rn×n are the mass and stiffness matrices, C=αcM+βcK is the damping matrix, x=xt∈Rn×1 is the system’s displacement vector, and f=ft∈Rn×1 is the excitation vector of the system.Φ
=
φ
1
,
φ
2
,
…
,
φ
n is assumed to be the eigenvector matrix of equation (29). The displacement vector can be described in the independent-mode coordinates q=q1,q2,…,qnT as(30)xt=∑i=1Nφixqit,where φi is also the ith mode shape and N is the order of the system. Therefore, one obtains(31)miq¨i+ciq˙i+kiqi=fqi,i=1,…,N,where mi=φiTMφi, ki=φiTKφi, ci=αcmi+βcki, and fqi=φiTf.Equation (31) is a group of independent mode functions of the manipulator system.


## 4.2. Model Reconstitution

In general, it is easy to get the independent functions in equation (31) and substitute the solutions into equation (30); then, the solution of the system is obtained. But there is a problem that complicates attainment of the correct solution. Some eigenvalues in equation (29) are negative, causing the divergence of the manipulator system. This problem never appears in pure vibration systems, but in case, when the large-scale variables of the manipulator are considered in the function, the matrices cannot be well-defined. According to the decoupled function of equation (31), the mode mass mi is negative corresponding to respective eigenvalues. A simple method for dealing with the problem and maintaining the systems convergence is to reconstitute the model. The main process of the method was to make positive the negative mass of the mode. The new mass matrix can be expressed as Mqrec=diagmi, i=1,…,N.Therefore, the reconstituted mass matrix of the original system is(32)Mrec=Φ−TMqrecΦ−1.The reconstituted function was derived as(33)Mrecx¨+Crecx˙+Krecx=f,where Krec=K and Crec=αcMrec+βcKrec.The mode vectors of equation (33) are the same as those of the original system shown in equation (29) and whether the mode frequencies conform to reality can be validated experiment or experience.


## 4.3. Simulation Results

The manipulator was loaded with the external momentMr, generated by the driving motor. The damping coefficients were αc=0.1 and βc=10−4. Different types of the applied moment are listed in Table 6, and the response results are shown in Figures 4–7.Table 6
Different external excitations.
Forms of the torque
Mathematical expression
No. of response figures
Impulse excitation
M
r
=
100
,
t
=
0.5
s
0
,
else
Figure4
Step excitation
M
r
=
0
,
t
<
0.5
s
10
,
else
Figure5
Harmonic excitation
M
r
=
sin
2
π
⋅
1.6
⋅
t
Figure6
Hybrid excitation
M
r
=
sin
2
π
⋅
1.6
⋅
t
+
sin
2
π
⋅
2
⋅
t
+
w
t
,
where
w
t
is random excitation
Figure7Figure 4
Response with impulse excitation.
(a)
(b)Figure 5
Response with step excitation.
(a)
(b)Figure 6
Response with harmonic excitation.
(a)
(b)Figure 7
Response with hybrid excitation.
(a)
(b)The bandwidth of the impulse excitation was 2 ms; the responses in Figure4 include rigid rotation and flexible vibration. The damping coefficient of the manipulator system was small, so the rotating angle of the response could not reach a constant value, and it was also hard to attenuate the amplitude of vibration. Figure 5 exhibits the response of the manipulator for the step excitation.Vibration frequencies were similar to those of the impulse response, including the first two modes of the system. The step excitation caused the rigid rotation of the system to accelerate. At the same time, the vibration was changed to a DC (direct current) bias. The harmonic and hybrid excitations included the first two eigenfrequencies of the system. As shown in Figures6 and 7, the responses of the system with the flexible joint and flexible link diverge, when the frequencies of excitations are close to the first two natural frequencies of the system. This phenomenon is in accordance with characteristics of a small damping system. The system also demonstrated its ability to isolate high-frequency vibration induced by the motor; the response of the system to the hybrid excitation that included white noise illustrates this.The simulated responses of the system indicated that the vibration of the flexible beam could be induced easily by most types of excitations, and the response was hard to attenuate. The vibration regime almost included the first two eigenmodes of the system. According to the character of these responses of the system, it is necessary to design a vibration controller for it.


## 5. Design and Simulation of the Controller

### 5.1. Controller Design

According to the simulated response, vibration was inevitable when the manipulator was driven by the motor. The flexible-beam vibration may reduce the tracking precision of the system; so, it is necessary to design a controller to suppress vibration. The function of the system under control was derived from equation (33):(34)Mrecx¨+Crecx˙+Krecx=f+Buu,where u=ut∈R3×1, Bu∈RNr×3 is the distribution matrix of the control force and Nr is the dimension of the reduced system. Defining x1=x and x2=x˙1 and substituting into equation (34), one obtains the following form:(35)X˙=AX+Bu+w,where X=x1x2, w=0Mrec−1f, A=0I−Mrec−1Krec−Mrec−1Crec, B=0Mrec−1Bu, 0 is the Nr-order zero matrix, and I is the Nr-order unit matrix.Vibration control focuses on suppressing the deformation induced by fast variables; so, there is no need to consider the slow variables—the rotation angleβ and θ—in the control function. Therefore, the observation equation for vibration can be presented:(36)y=HX,where H=0Nr−2×2,INr−2,0Nr−2×2,INr−2, with the subscript presenting the matrix order.The solution of equation (36) is the eigenform without rigid rotation. It can also be described in the form of partial states. H+ is the generalized inverse matrix of H; therefore, the eigenvector can be described as X=H+y; substituting it into equation (35), one obtains(37)y˙=A¯y+B¯u+Hw,where A¯=HAH+ and B¯=HB.In order to facilitate numerical simulation or experimental application, it is necessary to obtain the discrete form of equation (37). Assuming T as the discrete sample time, one obtains(38)yk+1=A¯dyk+B¯duk+Hdwk,where A¯d=expA¯⋅T, B¯d=∫0TexpA¯τdτB¯, and Hd=∫0TexpA¯τdτH.The form of equation (38) is a standard eigenfunction; it facilitates designing a state feedback controller with the LQR optimal method. One obtains(39)uk=−Lyk,where L is the matrix of the feedback control coefficients, which can be calculated by Matlab function [L, S, e] = dlqr(A¯d, B¯d, Q, R) and Q and R are weighted matrices of y and u, respectively.


### 5.2. Simulation of the Controlled System

The control forces and moments were generated by the piezoelectric actuators and the electric motor. The control force vector can be defined asu=u1,u2,u3T. According to the installation position of the actuators and the motor, the distribution matrix of control force in equation (34) is confirmed as Bu=B10, where B1=100000010001. By combining equations (35), (37), and (38), the coefficient matrices of the output controlled function can be confirmed.Valuating the weighted matrices asQ=103⋅I2Nr−2 and R=I3, the controller coefficient L∗ can be calculated with the LQR optimal equation solver. Then, substituting equation (39) with L∗ into the discrete form of equation (34), the responses of the controlled system can be obtained (Figures 8–11).Figure 8
Response of the controlled system with impulse excitation.
(a)
(b)
(c)Figure 9
Response of the controlled system with step excitation.
(a)
(b)
(c)Figure 10
Response of the controlled system with harmonic excitation.
(a)
(b)
(c)Figure 11
Response of the controlled system with hybrid excitation.
(a)
(b)
(c)The simulation results indicate good performance of the controlled system. According to Figures8 and 9, vibration decayed much faster with control when the system was exposed to impulse or step excitations. The rigid rotation also agreed with that in the response of the system without control. Figures 10 and 11 show the response of the system to the harmonic and hybrid excitations, respectively. It is obvious that the controller can improve the stability of the system and keep the response converging to a low-vibration state.


## 5.1. Controller Design

According to the simulated response, vibration was inevitable when the manipulator was driven by the motor. The flexible-beam vibration may reduce the tracking precision of the system; so, it is necessary to design a controller to suppress vibration. The function of the system under control was derived from equation (33):(34)Mrecx¨+Crecx˙+Krecx=f+Buu,where u=ut∈R3×1, Bu∈RNr×3 is the distribution matrix of the control force and Nr is the dimension of the reduced system. Defining x1=x and x2=x˙1 and substituting into equation (34), one obtains the following form:(35)X˙=AX+Bu+w,where X=x1x2, w=0Mrec−1f, A=0I−Mrec−1Krec−Mrec−1Crec, B=0Mrec−1Bu, 0 is the Nr-order zero matrix, and I is the Nr-order unit matrix.Vibration control focuses on suppressing the deformation induced by fast variables; so, there is no need to consider the slow variables—the rotation angleβ and θ—in the control function. Therefore, the observation equation for vibration can be presented:(36)y=HX,where H=0Nr−2×2,INr−2,0Nr−2×2,INr−2, with the subscript presenting the matrix order.The solution of equation (36) is the eigenform without rigid rotation. It can also be described in the form of partial states. H+ is the generalized inverse matrix of H; therefore, the eigenvector can be described as X=H+y; substituting it into equation (35), one obtains(37)y˙=A¯y+B¯u+Hw,where A¯=HAH+ and B¯=HB.In order to facilitate numerical simulation or experimental application, it is necessary to obtain the discrete form of equation (37). Assuming T as the discrete sample time, one obtains(38)yk+1=A¯dyk+B¯duk+Hdwk,where A¯d=expA¯⋅T, B¯d=∫0TexpA¯τdτB¯, and Hd=∫0TexpA¯τdτH.The form of equation (38) is a standard eigenfunction; it facilitates designing a state feedback controller with the LQR optimal method. One obtains(39)uk=−Lyk,where L is the matrix of the feedback control coefficients, which can be calculated by Matlab function [L, S, e] = dlqr(A¯d, B¯d, Q, R) and Q and R are weighted matrices of y and u, respectively.


## 5.2. Simulation of the Controlled System

The control forces and moments were generated by the piezoelectric actuators and the electric motor. The control force vector can be defined asu=u1,u2,u3T. According to the installation position of the actuators and the motor, the distribution matrix of control force in equation (34) is confirmed as Bu=B10, where B1=100000010001. By combining equations (35), (37), and (38), the coefficient matrices of the output controlled function can be confirmed.Valuating the weighted matrices asQ=103⋅I2Nr−2 and R=I3, the controller coefficient L∗ can be calculated with the LQR optimal equation solver. Then, substituting equation (39) with L∗ into the discrete form of equation (34), the responses of the controlled system can be obtained (Figures 8–11).Figure 8
Response of the controlled system with impulse excitation.
(a)
(b)
(c)Figure 9
Response of the controlled system with step excitation.
(a)
(b)
(c)Figure 10
Response of the controlled system with harmonic excitation.
(a)
(b)
(c)Figure 11
Response of the controlled system with hybrid excitation.
(a)
(b)
(c)The simulation results indicate good performance of the controlled system. According to Figures8 and 9, vibration decayed much faster with control when the system was exposed to impulse or step excitations. The rigid rotation also agreed with that in the response of the system without control. Figures 10 and 11 show the response of the system to the harmonic and hybrid excitations, respectively. It is obvious that the controller can improve the stability of the system and keep the response converging to a low-vibration state.


## 6. Conclusion

This paper solves the vibration problem for a manipulator with a flexible link and a flexible joint. The main stages of research can be summed up as follows: (1) establishing a FE model of the manipulator with a flexible link and joint; (2) reducing the dimension of the FE model with a reduction method; (3) suggesting the CRRM to update the reduced FE model and fitting the dynamic characteristics of the system to reality; (4) presenting a model reconstruction method to correct the CRRM mass matrix and ensuring the nondivergent response of the system; (5) establishing the state function of the output system, designing a LQR optimal controller for the reconstructed model, and performing simulations of the controlled system.All simulation results indicated good performance in terms of dynamic characteristics and controlled effect. The response of the reconstructed model without control showed that the joint of the manipulator could prevent high-frequency vibration being transmitted to the flexible link. The vibration was almost induced by the first two eigenmodes of the reconstructed model. The vibration phenomenon was in accordance with reality. The controller design only considered the vibration signal, ignoring the rigid motion of the manipulator. Therefore, the response of the controlled system showed good vibration suppression and maintained excellent performance in terms of the rigid motion. The method suggested in this paper can be used for practical applications. It is worth noting that the FE model should be adjusted according to the specific model, mainly by mode adjustment. The mode adjustment for manipulators and their experimental analysis will be studied next.



---

*Source: 4671358-2020-01-31.xml*