A model-based control law for vibration reduction of serial robots with flexible joints

The presence of flexibilities in rotational joints can limit the kinematic performances of manipulators doing high speed tasks as Pick and Place. The problem addressed in this work concerns the vibration control of serial robots with flexible joints performing Pick and Place tasks in order to improve productivity. Based on a dynamic model of a robot with flexible joints, a model-based control law is proposed with its associated tuning methodology. The robot dynamic model is then the key point of ourmethodology. This dynamic model considers stiffness and damping of each flexible joint. To guarantee its accuracy, a geometrical and dynamic identification procedure is realized. The objective is to show the relevancy of the proposed approach which integrates joint flexibilities in the control law. Theoretical results based on a representative model are used to illustrate the benefit of this model-based control law compare to two other control strategies (Feedforward control and control dedicated to rigid structures). Finally, a sensitivity analysis of this control law is realized to quantify the impact of modelling error and conclude on the criticality of joint damping value on vibration decreasing.


Introduction
Nowadays, the demand of productivity requires robots to present a control behaviour in terms of geometric and dynamic accuracy and task execution time. In this paper, we are focusing on serial robots for Pick and Place tasks.
Pick and Place tasks requires pose accuracy, repeatability, orientation accuracy, reduced stabilization time, position overshoot and static compliance [1]. One of the current challenges for this type of task is to increase productivity while ensuring high accuracy for the final pose.
To achieve this objective, lightweight robots are developed; however resulting flexibilities lead to a high sensitivity to different kind of loaded [2]. Thus, during high acceleration motions, vibratory phenomena appears [3]. Main causes of these vibrations are inertial forces, external loads coupled with the mechanical behaviour of the robot structure and control system. Several researches are carried out to improve the accuracy of trajectory tracking by considering flexibilities in joints and links according to the influential elements [4,5].
Robot vibratory decreasing can be achieved with different viewpoints: task pose in the robot workspace [6,7], tool path interpolation [3,8] or control law [9]. In this article, we focus on the influence of control law and identification of influent parameters in vibratory behaviour of Scara robot during a Pick and Place task.
In practice, this vibratory behaviour is due to the flexibility in motorized joints transmission system or/and the flexibility induced by the lightening of moving arms [10,11]. Moreover, these elements may cause low resonance frequencies [12,13].
This issue can be observed with a first experimental test realised on a Scara S600 robot, set in motion by an industrial control law, using an external measurement sensor (Laser tracker) with an accuracy level of 15 mm and an acquisition time period of 5 ms. The same end-effector path (movement of the two first axes) was realised with different speed levels from 10 to 80 rad/s. Vibrations appear at the end of the movement phase. The amplitude of the vibrations at the end of the movement is related to the robot's movement speed with a maximum of 0.15 mm as shown in Figure 1.
Different researches are carried out depending on the type of flexibility parameters and models and compensated by the control law [9]. In this article, we work on industrial robot with flexible joints. Indeed, the presence of flexibility in the motion transmission systems of industrial robots is classic due to the use of belts, long shafts, cables, harmonic drive reducers or gears [9]. In very small to medium-sized industrial robots, harmonic drives are frequently used; these drives have low backlash [9]. Thus, backlash is neglected, and classical control law only compensates friction. The difficulty in setting up a control law to compensate joint flexibility is based on the fact that the system for measuring the position and speed of each joint is generally placed before the system of transmitting movement [9]. Thus, there is no direct measurement of the elastic deformation of the joints and its compensation requires to develop dedicated models.
Multi-body models were established by adding generalized coordinates, corresponding to the deformation of each spring, to the rigid body model [15]. Moberg [16] explained that the improvement of tracking accuracy of a heavy robot requires considering the non-linear behaviour of the transmission system (hysteresis, backlash, friction and non-linear flexibility), and radial and axial joint flexibilities. In addition, the consideration of non-linear behaviour in the control requires the development of complex strategy [17]. Hamon, added a friction model at the transmission system level [18]. Different friction models are currently available in the literature [19]. Thus in many research projects on light weight robots, the flexibility of the joints is often modelled by a torsion spring taking into account dry friction, damping or viscous friction [9,20,21]. To identify stiffness of the joints, Dumas suspended weights at the end effector and measured the deformation using a laser tracker [4]. Other methods based on the analysis of vibration responses were proposed in Gunnarson [22]. However, there is no discussion on the influence of model flexibility parameters and identification errors on vibration decreasing. Thus, the aim of this paper is to identify the critical flexibility parameters for vibration decreasing of a robotic cell (robot and control law).
From this literature review and according to Scara robot behaviour, it seems relevant to authors to model joint flexibilities by a linear torsion spring associated to a damping, a dry friction and a viscous friction. The dynamic model is then established using Lagrange's formalism [23,24]. This model aims to study the influence of control strategy on the vibratory behaviour of the robot.  This article aim is based on the study of the impact of damping and stiffness parameters on the control law performance regarding vibration decreasing and identification errors.
Several control laws have been proposed to solve the problem of controlling trajectory of flexible joint robots. Control strategies proposed are within the field of nonlinear and multivariable system control. The particularity of such robots is that elastic deformations are not directly measured and can generate vibration frequencies that limit stabilization time [25,26]. Control techniques like those used for rigid robots can be proposed. Thus, Tomei varied gains of a conventional PD controller to reduce final pose accuracy [27]. An extension of the PD controller taking into account the dynamics of the actuators and friction is proposed by Lozano [23]. In the literature, other more or less complex control strategies are developed to improve the performance of Pick and Place operations. The singular disturbance method can be applied to take joint flexibility into account in the control scheme [26]. Flexibility is then integrated into the control diagram as a disturbance. The loop linearization method allows a non-linear model to be transformed into a linear and decoupled system. Other work focused on the implementation of advanced control strategies such as adaptive, optimal, robust predictive control laws based on neural networks [22].
The contribution of our work is to propose a modelbased control strategy in order to quantify the impact of the particular behaviour of the motion transmission chain (stiffness and damping) on the accuracy of the end effector in terms of precision and stabilization time. The robot under study is a Scara robot S600 illustrated on Figure 3.
In order to be able to extend our work to all anthropomorphic robots, we are interested in the first two axes, which are rotational joints and allow the end effector to move in a horizontal plane.
This paper is organized as follows: modelling and identification of the geometric and dynamic parameters including damping, friction and stiffness is described in Section 2. The vibration control approach is developed in Section 3. Comparison with two other control strategies is also achieved in order to show the effectiveness of the proposed control strategy. Finally, a sensitivity analysis is carried out in Section 4. Conclusion of this research work are presented in Section 5.
2 Modelling and identification of geometric and dynamic parameters

Geometric and dynamic modelling
This section provides a brief description on the modelling of the robot. The robot segments are perfectly considered rigid and joints are considered geometrically ideal. Geometric parameters of the robot are presented in Figure 4. O 1 and O 2 are centers of joint 1 and joint 2, G 1 and G 2 are centers of mass of solid 1 and solid 2, q 1 and q 2 are articular values of joint 1 and joint 2, L 1 and L 2 are the length of arm 1 and arm 2, and Lc 1 and Lc 2 are the distance between arm center of mass and joint center of arm 1 and arm 2.
The direct geometric model of the Scara robot is written as: This model is linear with respect to dynamic parameters [28]: where G mid ; q; _ q; € q respectively represent vectors of motor torques, positions, speeds and joint accelerations. A(q) and C q; _ q ð Þare two square matrices, the inertia matrix and the Coriolis matrix. F s and F v are vectors containing dry and viscous friction coefficients of both joints.
For this case study, we obtain the following model:

See equation below.
where m 1 and m 2 are the respective masses of links 1 and 2. I 1 and I 2 are the inertia of solids 1 and 2 along the z-axis at their centre of mass G 1 and G 2 .

Dynamic modelling with flexible articulation of the Scara S600 robot
The robot is modelled as a combination of rigid solid set in motion between them by motorized joints whose transmission system is represented by a spring with a stiffness k i and a damping coefficient c i . We note q i the instruction sent to the articulation motor and q ei the movement transmitted to the solid after the flexible motion transmission system (Fig. 5).
The equations of motion in matrix form are written as: where j and K are respectively damping and stiffness matrices such that: The system's natural frequencies v i and mode quation: The accuracy of these models is obtained from a modal identification procedure. In a first time, a geometric identification is realized before the identification of dynamic parameters.

Geometric identification
The geometric parameters of the model are arm lengths L 1 and L 2 . Measurement of the end effector position is carried out using a Leica AT901 Laser Tracker, the target being positioned at the end of the robot at axis 4. According to the Scara geometrical behaviour, we choose to perform the identification using points measured statically during two types of displacements: -Displacement of q 1 with fixed q 2 (Fig. 6a). The end effector follows a circular trajectory of center O 1 . 20 positions are measured over an angle near 45°. A least squares plane P of all these measuring points is then determined as the displacement plane of the end effector. The measured points are then projected into this plane before calculating the least squares circle of these points. We thus determine the position of its center point O 1 in the measuring frame and its radius R.
-Displacement of q 2 with fixed q 1 (Fig. 6b). Trajectory followed by the end-effector is a circle of center O 2 and radius L 2 . 14 positions are measured over an angle near 40°. By projecting these points into the plane P and calculating the associated least square circle, we determine the length L 2 and coordinates of O 2 in the measurement frame.

Dynamic parameters identification
Dynamic parameters are different kind of quantity (stiffness, damping, inertia, friction and position of mass centre). These parameters have to be identified with different type of data.

Identification of joint stiffness and damping
The identification of the joint stiffness and damping can be done from the determination of system's natural frequencies by modal analysis (Fig. 7). The modal analysis is performed using LMS-TAB software.
Three configurations have been chosen to measure natural frequencies of the robot: extended arm (q 2 = 0°), arm bent at 45°(q 2 = 45°) and arm bent at 90°(q 2 = 90°). Once the accelerometers are in place, the measurement protocol is done in two steps: Step 1: we apply an excitation force at the end effector using an impact hammer.
Step 2: After choosing the bandwidth for signal acquisition, the frequencies are then selected manually by using PolyMax identification procedure. The software proposes unstable and stable modes representative of numerical or structural modes for the results. Then, only the so-called "stable" modes are selected.
Observation of the modal shape of the model shows that the first 2 modes are associated to the deformation of joints 1 and 2.
For the identification of stiffness and damping parameters of joints 1 and 2, the first 2 modes are considered (Tab. 1).
The others dynamic parameters are then identified from motor torque measurement.

Dynamic identification
According to equations (2), (10) parameters must be identified for the rigid dynamic model: The identification procedure is done in four steps using the same trajectory performed at different speeds: Step 1: The robot goes from position A to position B in the joint space.  Step 2: The torque is not directly measurable. It is deduced from the measurement of the motor current by the following relationship for axis j: where v j is the value of the motor current, R j the reduction ratio of the joint, g j the static gain of the amplifier, K T the torque coefficient of the motor.
Step 3: We calculate the torque G mdi from equation (2). To do so, we measure the angular position with the encoder of each motorized axis. The estimation of acceleration and speed is then computed by numerical derivation and Butterworth filtering of this last measurement [29,30].
Step 4: The identification procedure is achieved by minimizing the following cost function: First of all, the dry friction parameters were identified using a low speed movement (axis by axis). The trajectory is composed of two movements: firstly, movement of q 1 with fixed q 2 and secondly, movement of q 2 with fixed q 1 (Fig. 8).
For the rest of the parameters, trajectory is the same as before but with higher speed. Measured torques during movement of axis 1 and then axis 2 are isolated and used to perform the identification. After optimization (Fig. 9), the mean error is 0.648 N m and the standard deviation is 1.676 N m.   To identify the stiffness of the joints, we couple the cost function expressed by equation (6) with another cost function which compares the natural frequencies measured with the modal analysis and estimated by the model (Eq. (4)). Values of these modes varies according to q 2 as the inertia matrix is a function of q 2 (Eq. (2)). Identified parameters are found in Table 2. The parameters already identified in the previous section are validated with a random trajectory made by the end effector. The mean error is 1.4 N.m and standard deviation is 3 N m.
To finalize the identification process, damping parameters values have to be determined.

Damping identification
The modal analysis allowed us to determine the two modes corresponding to deformations of joint 1 and joint 2. Damping ratios are evaluated for these two modes. We only consider the first extended arm configuration for identifying joint damping because all accelerometers and hammer chock are oriented tangent to the displacement of joint 1 and 2. Thus, in this configuration, joints are directly loaded by hammer chock and their vibration amplitudes are directly measured by accelerometers. This makes it possible to link the measured damping to the torsional damping at each joint. Finally, damping coefficients of joint 1 and joint 2 are deduced from measured damping ratios (z) for these two modes according to robot configuration using the following relation: where & is the measured damping ratio matrix and C d , K d , A d respectively the diagonalized matrixes of joint damping, stiffness and inertia. Damping parameters are finally determined as follows, c 1 =57.484 N.m.s. rad −1 and c 2 = 12.983 N.m.s. rad −1 .
After the identification process, the identified model is used to define a model-based control law.

Model based flexible control (MBC) law
In this chapter we develop a computed torque control law taking into consideration joint flexibilities. Validity and benefit of this control law for flexible joint robot is evaluated through advanced simulations. All the values of identified parameters in previous chapter is integrated in the simulator in order to have a representation as close as possible to the real behaviour of robot. Before reaching this stage, we define a trajectory generator in the articular space, and develop a methodology for tune gains of the used corrector in order to compare similarly different kind of control laws.

Control law scheme
The equation of motion governing the variable qe is obtained from the relationship (3): where G m and G e are respectively the motor torque and the elastic torque, and Dq e = q À q e .
On the other hand, the equation of motion is written as: where q e = q À Dq e .
These two relationships (Eqs. (7) and (8)) are the basis of the new control law based on the dynamic model of a manipulator with flexible joints. Flexibility of joints is presented in the control diagram as a feedback loop (orange block). In absence of gravity effect, PD control strategies are asymptotically stable and do not require addition of an integrator [31]. The control scheme corresponding to this computed torque control law is defined in Figure 10.
To validate and to illustrate the benefit and limit of the proposed control law, firstly, we develop an advanced simulation with Adams' and Simulink. Secondly, we make a comparative analysis with two control strategies applied to same simulator. Simulations results are presented in a comparative table.

Adams simulink co-simulation unit
The numeric model of a robot with 2 degrees of freedom was developed in Adams software. The idea is to use its mechanism solver to simulate robot motion generated by discretised setpoint computed by the control law. The dimensions and dynamic characteristics correspond to the one identified on the Scara S600 robot of SIGMA Clermont's technological platform. The control laws are implemented in Simulink. The numeric model defined in Adams is set in motion via a co-simulation. The flexibility of the connections is represented under Adams by a linear torsion spring modelled by stiffness and damping. Two segments of negligible mass compared to that of the two robot axes (m 1 = m 2 = 0.15 kg) are added to perform the simulation with flexible connections. The movement of these two segments is controlled. The segments are in rotational connection with the robot elements and the movement is transmitted by the springs. The robot model is defined geometrically and dynamically using the parameters identified above.

Trajectory generation and gain correction
To be consistent with a classical robotized path, the trajectory followed in this work is described with a fivedegree polynomial function, which ensures a position, velocity and acceleration movement continuities without excessive dynamic load (Fig. 11). The maximum velocity and acceleration are those of Scara robot (max velocity:  Gain tuning is realized to avoid oscillations, and in accordance with the response time at 5% (t r 5%) of position loop. Thus, an appropriate adjustment [20] is: where K p and K d represents respectively proportional and derivative gains of the PD corrector. In the following paragraph, we develop a robot position control using the co-simulation tool Adams and Simulink. Initially, the robot is positioned in a configuration where arms are perpendicular. Initial and desired operational and Cartesian positions of the robot end effector are defined in Table 3.

Modal based control (MBC) law
This control law is our new computed torque control law for flexible joint manipulator. The gain setting used is the one defined in equation (9). For this law, we choose a response time at 5%(t r 5%) of 0.0143 seconds for an adequate setting. Thus w d ¼ 5 tr 5% ¼ 350, the gains obtained for the two rotation axes are K d1 = K d2 = 700 and Kp 1 = Kp 2 = 122500. Evolution curves of the positioning error in both joint and Cartesian spaces are shown in Figure 12a.
Desired position is reached with a maximum trajectory tracking error of 1 millimetre along x-direction and half a millimetre in y-direction. The mean errors and standard deviation are in the order of 10 À4 meters and 10 À6 radians in Cartesian and articular spaces, respectively. At the end of the task, robot's behaviour is stable, positioning error is 10 À9 meters and an accuracy of 1 mm is achieved after 0.913 s.

Rigid control
In order to show the benefit to consider flexibilities in the control loop, a control law dedicated to rigid robots (Fig. 13) is applied to the same Adams robot model with flexible joints.
The gain tuning used to minimize the error is defined in the equation (9). For a t r 5% of 0.0143 seconds (identical to that of the MBC control law), joint flexibilities are highly stressed, vibrations appear, and desired positions are never reached because the robot never converges to a fixed position. The system stabilizes for a t r 5% of 0.3846 seconds, and the proportional and derivative gains are 169 and 26, respectively.

Feedforward control
In this section, a feedforward control combined with a PD control is applied to the same simulator (Fig. 14). An adequate adjustment [19] of the K p and K d gains to control the transient behaviour of the error is: where the element a 11 designate maximum value of matrix A max and has the value 1.9117. Similarly, c 11 has the value 1.14. For a proper gain adjustment, t r 5% is fixed at 0.103 seconds. Thus, the gains obtained for both rotational axes are: Kp = 182 and K d = 4441.
Results of the simulations are presented in Figure 12 and in the comparative balance sheet table (Tab. 4).
As a quick interpretation, MBC control law (considering joint flexibilities as a feedback loop in the control diagram) provides an additional performance advantage over all other control laws in terms of end effector placement accuracy and stability time with a benefit of 10% and 50% according to feedforward and rigid control respectively (Pick and Place operation requirements).
However, model-based control law can be sensitive to identification errors. A sensitivity analysis ensures to determine critical parameters for such control law strategy.

Sensitivity analysis
This sensitivity analysis ensures to determine critical flexible joint parameters with regard to vibration decreasing. The fact that we use a model-based control law strategy allows us to only focus on the influence of mechanical parameters used to describe joint behaviour.
The torques applied to flexible connections, according to estimated dynamic parametersx, are obtained from the following model: whereÂ ¼ Ax ð Þ etĈ¼ Cx ð Þ and ∧ is the symbol of estimation notion. Thus, control law U is defined by the following relationship: where e ¼ q d À q e ; _ e ¼ _ q d À _ q e ; € e ¼ € q d À € q e Let consider that equation (3) is equal to equation (11):  By substituting U with its expression in equation (12) and eliminating the terms of dry and viscous friction on both sides of this equality, equation (13) Let:Ã ¼ A ÀÂ andC ¼ C ÀĈ equation (14) can finally be written as follows: Normally, the error and its derivative converge towards 0 in a second order determined by the servo gains. However, in equation (15), the second member is not null but depends on identify values of inertial parameters and joint flexibilities (stiffness and damping). Hence, this equation illustrates the need to take into consideration the joint flexibility in the control. These three parameters are subject to estimation inaccuracies. A sensitivity study is conduced to quantify the impact of these errors on the vibration behaviour. This sensitivity analysis is realised by isolating each source of error in order to determine its influence separately by numerical simulations.

Errors in the estimation of inertia
In order to study the influence of inertial values on the performance of our new model-based control (MBC) law, we realise several simulations of the MBC control law (with same settings of PD gains) with 2% error steps on inertial parameters. For errors below 12%, simulation results confirm that the control is able to adapt avoid vibrations at the end of the movement. Starting from 12%, vibrations related to instability phenomena appear.

Errors in joint stiffness
The relevance of MBC is evaluated for a case where values of joint stiffness present defect with regard to implemented one in Adams model. In this purpose, we apply the same test as before but with a 1% error step (since the impact of these errors is then greater). Simulations reveals that vibrations appears from errors of 3% in joint stiffness. These vibrations can be controlled by minimizing the performance of the control law due to an increase in response time at 5% (t r 5%) and consequently a decrease in the corrector gains.

Errors in joint damping
The impact of these errors is significantly very important. Vibrations appear from a 3% error on damping. These vibrations can be controlled by minimizing the performance of the control law due to an increase in response time at 5% (t r 5%) and consequently a decrease in the corrector gains. Table 5 resume all the sensitivity test and the decreasing of performances due to an increasing of t r 5%. Thus, we can show that damping and stiffness is critical parameters for vibration control. However, there is work on stiffness joint robot identification and just few on damping. Thus, this work demonstrates a huge need to develop accurate methodology to identify joint damping parameter.

Conclusion
In this work we have developed a new model-based control law dedicated to manipulators with flexible joints to control vibrations behaviour during Pick and Place operations. Simulations are applied on a Scara robot with 2 Degree of freedom, so that this law is applicable to more complex structures. We have modelled and identified the geometric and dynamic parameters of a Sacra robot including friction, stiffness and damping. MBC law dedicated to robots with joints flexibility is then developed using Simulink after defining the trajectory and PD correction system. Advanced simulations results in terms of response time and maximum and average trajectory errors clearly shows the relevance of the proposed approach integrating joint flexibilities as a feedback loop in control diagram. In parallel, its performance is validated due to comparisons with other control laws. Finally, we concluded this work by a sensitivity analysis to modelling errors. We demonstrate, by calculations, that joint flexibilities and damping are influential on the modelling error which indicates the importance of taking them into consideration for the control accuracy of robot vibrations.