Calibration of UR10 Robot Controller through Simple Auto-Tuning Approach

: This paper presents a calibration approach of a manipulator robot controller using an auto-tuning technique. Since the industry requires machines to run with increasing speed and precision, an optimal controller is too demanding. Even though the robots make use of an internal controller, usually, this controller does not fulﬁll the user speciﬁcation with respect to their applications. Therefore, in order to overcome the user requirements, an auto-tuning method based on a single sine test is employed to obtain the optimal parameters of the proportional–integral–derivative PID controller. This approach has been tested, validated and implemented on a UR10 robot. The experimental results revealed that the performances of the robot increased when the designed controller, using the auto-tuning technique, was employed.


Introduction
Manipulator robots can be described as automated electromechanical systems with flexible functionalities (i.e., programming according to the environmental conditions in which they operate).The main role of the manipulator robot is to help the human operators in fulfilling repetitive tasks with increased risk in an industrial environment.From the perspective of robot control, focus is put on the structure and functionality of the robot command unit.Based on the dynamic and geometrical model and the tasks that need to be performed, appropriate commands are established.In order to provide these commands to the actuators, the hardware and the software, it is required to use the feedback signals obtained from the sensory system unit [1].Due to the complexity of the manipulator robot, the control architecture has a hierarchical structure.The upper level makes the decision with respect to the actions that need to be taken (e.g., simple test based on "if-then" actions), while the lower level carries out the control of the joints.The typical structure of the control system consists of a computer on the upper level and a system with one or more microcontrollers to drive the actuators of the joints on the lower level.
Mechatronic systems, in special robots, are very popular for control applications due to their interdisciplinary nature [2][3][4].For linear mechatronic systems, the proportional-integral-derivative (PID) controller has been widely used given its simple structure and robustness [5].Usually, the design of conventional integer-order PID controllers is based on the model of the system.
The manipulator robots are also used in medicine for precision surgery [6][7][8].Surgery has become an extremely precise practice, and professionals often have to make movements in the order of hundreds of microns.The inherent limits of human dexterity are a constraint for new advances in this field.In this context, high-precision robotics are being developed and implemented for medical use [9,10].From the observation of these devices, it is possible to acquire an idea of their versatility.Not only in medicine, but in the high-tech industry, as well, robotic automation is highly competitive and requires very high precision and better performance.Therefore, the performance specifications of control have also become extremely demanding.
In these real applications of the UR10 robot, the characteristics of the system may be subject to changes (e.g., the mass attached to the end effector, the speed of the joints (TCP respectively), the required precision, etc.).The controller, however, has to keep its performance constantly optimal.Usually, the robots already have an internal control loop, making the robot stable and performing quite well.Nevertheless, this controller cannot fulfill all the user specifications and does not perform optimally in all situations.Therefore, a new controller is employed in order to overcome the user requirements [11][12][13][14].Hitherto, PID control is widely used in control engineering and industry.The most challenging step in employing PID controllers, for non-experts in control engineering, is the process of parameter tuning.Nowadays, the self-tuning PID digital controller provides much convenience in engineering [15,16].Optimal control of a plant (in this particular application, the robot arm) is highly dependent on the plant behavior.
Throughout the years, several methods have been introduced [17][18][19] to improve the tuning principles proposed by Ziegler and Nichols.Astrom and Hagglund proposed an AMIGO (approximate M-constrained integral gain optimization) tuning rule for PI controllers [20], based on the ideas presented in [21], where the tuning of the controllers is done so as to maximize the integral gain subject to constraints on the maximum sensitivity.The same approach was later extended to PID controllers [22] and finally perfected by adding additional constraints [23] to make the tuning method suitable for a wide class of processes.Skogestad proposed a novel tuning rule based on the ideas behind IMC (internal model control) that have achieved widespread industrial acceptance [24].These novel tuning rules, referred to as SIMC (Skogestad IMC), work well for both integrating and pure time delay processes, for both setpoint changes and load disturbances and are suitable for first-order or second-order time delay models.
In this paper, we propose simple and efficient calibration of the robot controller based on an auto-tuning procedure.The auto-tuning method attempts to design PID controllers that ensure a minimum phase margin PM = 45 • and a minimum gain margin GM = 2.The method is based on defining a forbidden region in the Nyquist plane, using the performance specifications, and then searches for the optimal PID controller that ensures that the difference between the slope of the region border and the loop frequency response slope is minimum.This requirement leads to PID controllers that are robust to gain, as well as to phase and delay variations.Although the method has already been compared to various other popular PID tuning methods [25], its main features have not been fully revealed.To design the PID controller, a single sine test is applied to the robot to determine its frequency response and its frequency response slope [26], needed in the auto-tuning procedure.The original elements of this paper include, apart from a detailed description of the auto-tuning method, also the particularities of using it in a closed loop setting.Moreover, for the first time, the auto-tuning method is validated experimentally.
This paper is structured as follows: The next section gives a description of the robotic system used in this study.Section 3 discusses the automatic calibration of the robot using an auto-tuning approach.Section 4 presents the implementation and the corresponding results of the tuning procedure.The conclusion is formed in Section 5.

System Description
A robot manipulator, depicted in Figure 1, can be defined as a kinematic chain consisting of multiple rigid bodies, which are interconnected by joints.Every joint of the kinematic chain consists of one degree of freedom, which can be translational or rotational.In the work of Spong et al. [1], it is mentioned that a simple representation of a robot manipulator joint can be modeled as a spring with linear constant stiffness.Based on the Lagrangian formulation, the dynamics of a global manipulator robot with n revolute joints are given by [27]: Q = M(q) q + C(q, q) q + F( q) + G(q) + J(q) T η (1) where q, q, q ∈ n represent the joint variables (the generalized joint coordinates, the velocity and the acceleration), M(q) ∈ n×n is the inertia matrix, C(q, q) is the n × 1 vector of the Coriolis and centrifugal forces, F( q) is the friction torque, G(q) is the gravity force and Q is a vector representing the actuator forces associated with the joint coordinates q.The matrix J(q) T is the transpose of the Jacobian matrix of the robot, while η is the joint force vector applied at the robot end-effector.The matrices M, C, F and G are composed of complex functions, which describe the kinematic link dynamics through a set of parameters (θ j , d j , a j , α j ), j = 1, ..., n, also known as Denavit-Hartenberg parameters.
The considered UR10 robot studied in this paper has the following Denavit-Hartenberg (D-H) parameters; see Table 1.The highest variation in the robot dynamics is introduced by the gravity matrix and the inertia matrix.The gravity term is present even when the robot is stationary or moving slowly.The torque exerted on a joint due to gravity acting on the robot is strongly dependent on the position and orientation of the robot.
Due to the construction of the robot (three large joints and three small joints), we expect to have a larger variation introduced by the shoulder joint and elbow joint.For example, the gravity torque acting on the elbow joint is higher when the robot is in a position such that the elbow joint has to support the shoulder-lift joint and the wrist joint.To illustrate this aspect, we have investigated the variation of the gravity load for Joint 2 (shoulder-lift joint) and Joint 3 (elbow joint) with respect to the robot configuration.The obtained results are illustrated in Figure 2. As observed, Joint 2 has a larger variation (±100 Nm) compared to Joint 3 (±40 Nm).This analysis is important for robot design, as well as for control design in order to ensure a feasible control input for the motors.
The gravity matrix is not the only term that varies with respect to the position of the robot.The other matrix that is highly dependent on the robot pose is the inertia matrix.By inspecting the elements of this matrix, we observe that it is symmetric with diagonal terms M jj describing the inertia related to joint j.The other elements M ji of the inertia matrix represent the coupling between joint j and joint i.In Figure 3, we can see the variation of the Joint 1 inertia with respect to the pose of the robot (here, Joint 2 and Joint 3 are varying between −180 and +180).In most of the robotics applications, one needs to control the pose of the end-effector, which can be done through the direct kinematics formulation.In this research study, a six-degree of freedom manipulator robot is considered; thus, the direct kinematic equation is given as: where i−1 A i (q i ), i = 1, 6 represent the homogeneous transformation matrices.The direct kinematics can be calculated via the Denavit-Hartenberg (D-H) formulation [28].
Starting from the assumption that the pose, the velocity and the acceleration are known and based on Equation (1), it is possible to calculate the required joint torques.
Given the pose configuration of the end-effector, the joint variables corresponding to this configuration can be calculated using the inverse kinematics formulation.The inverse kinematics is a very challenging problem since it is almost impossible to obtain a unique solution.In the literature, there are two approaches used to calculate the inverse kinematics: the analytical method and the numerical method [28].Here, a numerical method was considered [27].The existence of solutions is guaranteed only if the given end-effector position and orientation belong to the manipulator's dexterous workspace.
Nowadays, a large number of the industrial robots are driven by brushless servo motors.A schematic overview of a classical robot joint is illustrated in Figure 4. Several studies in the area of the control of manipulator robots consider as control variables the Cartesian coordinates.However, in order to obtain the desired Cartesian trajectory for the end-effector attached to the robot, each joint axis must follow a specific trajectory.The common approach to control the robot joints is based on the decentralized techniques, which consider an independent controller for each joint.The industrial robot considered for our experiments is electrically actuated.If the robot joint drivetrain is driven by the current: then the torque generated by the motor is proportional to the current: where K a is the transconductance of the amplifier, K m is the motor torque constant and u is the applied control voltage.Based on this assumption, the dynamics of a motor attached to a joint j can be defined as: with τ c Coulomb's law of friction, B the viscous friction and J m the total inertia seen by the motor for joint j, computed as: Since in real life, disturbance torques such as gravity and friction can act on the joints, the common approach where independent control systems are considered to control the robot joints is no longer efficient.The control performance decreases, leading to high overshoot and large steady-state error.A classical approach to overcome these shortcomings is to use a nested control structure composed of two loops: one outer loop and one inner loop, as depicted in Figure 5.The outer loop is used to control the position and to provide the velocity of the joints in order to minimize the position error, while the inner loop is used to minimize the error between the actual velocity of the joint and the velocity demanded by the outer loop.For simplification purposes, if the Coulomb friction from Equation ( 5) is ignored, the equivalent Laplace transform of ( 5) is given by: where Ω(s) and U(s) are the Laplace transform of their corresponding signal from time domain ω and u.Thus, the transfer function of a motor drivetrain attached to a robot joint can be written as: Typically, a PI controller is used to drive the robot joint from the actual velocity to its demanded velocity.However, in this experiment, it is no need for good closed loop dynamics; therefore, a P controller has been selected (note that the P controller was used only for tests in order determine the process frequency response and process frequency response slope) for simplicity in the computations in the next section, and thus, we have:

A Robust Method to Estimate the Process Frequency Response
The frequency response of a process at a specific frequency ω can be written as: where A y stands for the amplitude of the output signal y(t), A u is the amplitude of the input signal u(t) = A ω sin(ωt) and ϕ y is the phase of the output signal y(t).To determine experimentally the frequency response in (10), a robust procedure is needed to evaluate correctly, even in the presence of severe disturbances, the amplitude A y and phase ϕ y of the output signal y(t).Such a procedure is based on the transfer function analyzer-discrete Fourier transform (TFA-DFT) [29,30].
For the process with the frequency response given in (10), the steady-state response can be written as: where n(t) is the stochastic disturbance with zero mean value and b is a non-zero bias term (b = 0 in the case of a nonzero average disturbance).
Consider the multiplication of the measured output y(t) by sine and cosine functions of the frequency of the system input ω and then integrated over the measurement period T m = k 2π ω , with k an integer value [29,30].Then, the following results are obtained:  (13) An analysis of the right-hand side terms in ( 12) and ( 13) leads to a simplification, provided that there is an increase in the averaging time.In this case, the contribution of the last terms can be neglected compared to the first term, which is growing with T m .Furthermore, it is assumed that for a long integration time, the noise n(t) will be filtered out (i.e., zero average) [30].In this case, (12) reduces to: The result in ( 14) is valid even in the case of low signal-to-noise ratios provided that n(t) is a stochastic disturbance with zero mean value and a sufficient number of test signal periods is selected [29].
After some mathematical rearrangement of the terms in ( 14), the following result is obtained: Assuming that T m = k 2π ω , as mentioned above, namely the measurement (also integration) time is an integer number of periods, the last two terms in (15) will be zero.Then, a simplified result for (12) is finally obtained as: A similar analysis can be performed for the signal in (13), leading to: Using ( 12) and ( 13) leads to: The integral T m 0 y(t)e −jωt dt in ( 18) can be calculated via the DFT: where T s is the sampling period, chosen such that T m = NT s .Using now ( 16) and ( 17), the next result holds: which leads to: A y e jϕ y = y c − jy s −j0.5T m Using ( 18) and ( 19) in ( 21) leads to: A y e jϕ y = T s ∑ N−1 0 y(kT s )e −jωkT s −j0.5T m (22) Taking into account that T m = NT s and rearranging (22) leads to the following solution for determining robustly the amplitude A y and phase ϕ y of the output signal y(t) [26]:

A Robust Method to Estimate the Process Frequency Response Slope
For the purpose of designing the PID controller using the KC auto-tuning method [26], the frequency response slope needs to be determined, as well.For the process P(s), the frequency response slope is given as: A simple strategy to estimate (24) has been developed previously, and it is based on filtering the output signal y(t) as indicated in Figure 6 [31].The strategy is simple, but prone to errors in the case of stochastic disturbances and noise [26].An elegant solution to eliminate this problem has been developed and consists of a modification of the basic scheme in Figure 6 as indicated in Figure 7a, where y TR (t) and y SS (t) are the transient and steady-state parts of the output signal y(t), with the transient component going to zero (or to a constant value for an integrating system) [26].The steady-state component of the output signal can be written simply as: where S = A y cos(ϕ y ) and C = A y sin(ϕ y ).Applying the Laplace transform to (25) leads to: Applying the derivative to (26) leads to the following: The following result is obtained using the inverse Laplace transform on (27): Figure 7a can now be replaced by Figure 7b, using the result in (28).The Laplace transform of the signal x(t), computed based on Figure 7b, is obtained as: Applying the inverse Laplace transform to (29) and considering that ω s 2 +ω 2 = L −1 sin(ωt) leads to: However, the component t * y TR (t) does not influence the steady-state oscillation in y(t) and the result in (30).In this case, Figure 7b reduces to the simplified version in Figure 8.Using the result in Figure 8 leads to a means of computing the amplitude A y and phase ϕ y of the signal y(t) at the specific frequency ω as: A y e jϕ y = A x e jϕ x + A y sin(ϕ y ) Since the test frequency ω is user specified and the amplitude A y and phase ϕ y of the output signal y(t) are determined as mentioned previously, in order to determine the amplitude A y and phase ϕ y , the corresponding amplitude A x and phase ϕ x have to be computed, as well.These can be obtained from the frequency response of the system 2s ω Y TR (s), as indicated in Figure 8.Using the definition of the Laplace transform, the following relation is valid: Then, the amplitude A x and phase ϕ x are computed according to: where the discrete Fourier transform is used to evaluate the integral on the right-hand side of (33) as:

The Closed Loop Implementation
The robust procedures in the previous sections assume that an input signal can be applied to the process.In case the loop cannot be opened to enter a test-sine in the process P(s), a closed loop approach must be considered.For this purpose, a test-sine w(t) = A w sin(ωt) must be used, where w(t) is the reference signal of the closed loop system.Figure 9 illustrates this approach.The output y(t) is measured, and the error signal e(t) is further computed as e(t) = w(t) − y(t).Applying the procedure described in Sections 3. Based on (35), it follows that the process frequency response P(jω) can be easily determined as: where C(jω) = 15 is a P-type controller in the initial robot closed loop system.Computing the derivative of (35) leads to the following result: The result in (37) leads to:

An Auto-Tuning PID Controller Based on the Process Frequency Response and Frequency Response Slope
A simple auto-tuning method is used in this paper to design a PID controller using solely the frequency response and its corresponding frequency response slope of the robot at a specific frequency ω.This auto-tuning method is called the KC auto-tuner and produces robust PID controllers [25].The basic idea of the KC auto-tuning method is based on defining a circular region, including the −1 point in the Nyquist plane.This circular region is defined as a "forbidden region", which is a circle computed according to two design constraints: the gain margin (GM = 2) and the phase margin (PM = 45 • ); see Figure 11.The main design requirement is that the loop frequency response L(s) = P(s)C(s) needs to avoid entering the forbidden region, and most importantly, the robust PID controller parameters are determined in such a way as to ensure that the slope-difference between the circle border and the loop frequency response is minimum: where is the slope of the loop frequency response in the Nyquist plane and dIm dRe α is the slope of the "forbidden region" border, as a function of the angle α, as indicated in Figure 11.
, the slope of the loop frequency response has to be evaluated firstly as: The process frequency response P(jω) and the process frequency response slope dP(jω) dω ω=ω are determined as mentioned in Sections 3.1 and 3.2.Furthermore, for a user-specified frequency ω, the frequency response of the PID controller can be determined according to: whereas dC(jω) dω ω=ω is easily calculated in a numerical way for a known controller C(s) [25].
In this way, the slope of the loop frequency response in (41) can be easily broken into a real L and an imaginary L part as: Then, the slope of the loop frequency response in the Nyquist plane can be computed as the ratio . Regarding Figure 11, it is obvious that any point of the forbidden region circle is defined by the radius R and the angle α.Computing the real and imaginary parts of this point in the Nyquist plane, the following equations are obtained: The equation of the circle in that same point gives: with the derivative equal to: Then slope of the "forbidden region" is determined using the result in (46): where the result in (44) has also been used.With the results in ( 43) and ( 47), the minimization problem in ( 40) is solved with a single for-loop where α varies from 0-α max in 1 • steps.Referring to Figure 11, it is obvious that α max = 90 • .

Tuning of the PID Controller
To tune the controller, the test frequency is selected as ω = 22 rad/s, while the gain and phase margins are GM = 2 and PM = 45 • .Firstly, in the KC auto-tuning method, the process frequency response and frequency response slope at the test frequency are required.For this purpose, a reference signal as indicated next is used and applied to Joints 2 and 3 of the robot, similarly to the procedure described in Section 3.3: with A w = 5.The output signal is measured periodically, and the error signal is computed as: The reference and error signals are indicated in Figure 12. and it is given in Figure 14.Next, the discrete Fourier transform is used to evaluate the following integral, as in (34): Then, using the result in Figure 8 for the derivative of the error signal e(t) leads to a means of computing the amplitude A e and phase ϕ e of the signal e(t) at the specific frequency ω as: A e e jϕ e = A x e jϕ x + A e sin(ϕ e ) ω = −0.109+ 0.0042j with A e = 7.4 and ϕ e = −0.33 rad/s, as indicated in (50).With the result in (54), the following is obtained by modifying (24): The frequency response of S(s) at the test frequency ω is computed similarly to the equation in (10) as:

S(jω) =
A e A w e jϕ e = 1.4 − 0.48j (56) The process frequency response P(jω) can be easily determined using (36) and the result in (56), considering that the controller used is P-type with C(jω) = 15: The process frequency response slope is computed according to (39): The resulting PID controller, computed according to the KC auto-tuning procedure described in Section 3.4, is the following:

Experimental Results
This section reports the experimental results that were achieved in order to validate the auto-tuning method of a PID controller, presented in the previous section.The performance of the obtained controllers has been evaluated in the real-time application using the setup depicted in Figure 15.In order to communicate with the UR10 robot, the Robotics System Toolbox of MATLAB R and the ROS platform were used.The UR10 controller uses Ethernet and the TCP/IP protocol to send and receive commands.The control and calculations of the robot are executed in MATLAB R .Although it is possible to send commands directly from MATLAB R to the controller, it was chosen to use the ROS platform for a more versatile solution.Therefore, the ROS-supported ur_moder_driver runs on a computer with a Linux operating system, which is directly connected to the robot controller through TCP/IP.The ur_modern_driver allows an easy communication with many different operating systems and publishes robot data with the use of ROS messages.As MATLAB R has a built-in ROS support, this ensures a very clear and effortless communication; for more details, please see [32].An overview of the complete communication model is given in Figure 16.The running nodes and topics within the ROS platform related to this case study are illustrated in Figure 17.The procedure used to control the movements of the robot is presented in the flowchart in Figure 18.Notice that the auto-tuning method implies the use of ROS − controller and velocity_based_controller, as well.Besides controlling the movements of the robot, it is very important to have access to the status of the robot's joint states (joint position, joint velocity, joint effort, etc.).Thus, it can be mentioned that implementing a new controller on an industrial robot is strongly dependent on the opening level of the robot controller and the type of the driver used.In order to implement the auto-tuning method presented in this paper it is required to have read/write access to the joint states of the robot.
The speed control was developed and applied to Joint 2 and Joint 3 of the manipulator robot (these two joints have a larger variation compared with other joints).The controller was tested for multiple scenarios and different configurations of the manipulated robot.The first test was conducted using the following configuration [−80, −45, −50, −80, 260, −90], and the obtained results are presented in Figures 19a and 20a.The output of the manipulator robot for the second configuration ([35, −102, −75, 90, −120, 110]) considered in this paper is presented in Figure 19b and 20b.The third test case presented in this paper is illustrated in Figures 21 and 22.As observed, the controller performs equally for different joint configurations, while keeping the interaction between joints at a minimum.It can be seen that the PID controller designed using the proposed auto-tuning procedure converges to the imposed reference for all considered cases and presents no overshoot in comparison with the available controller.Given the problem addressed in this paper, the results provided indicate that the auto-tuning method can be successfully applied in the field of robotic control.

Conclusions
In this paper, a PID based on a single sine test auto-tuning method was designed and implemented for a robot system.Due to the advantage of the auto-tuning approach, which ensures that the difference between the slope of the region border in the Nyquist plane and the loop frequency response slope are minimum, the obtained PID controller is robust to gain, phase and delay variations.This robustness is highly demanded for a manipulator robot taking into account its dynamics (variation with respect to the position and orientation of the robot).The outcome of this paper shows that the proposed technique is a suitable method for tuning PID controllers.The experimental results indicate an increased performance of the manipulator robot when the proposed auto-tuning technique is used.Even more, this method could bring a benefit for both real-life and industrial applications given that it is time saving and is easy to use by non-experts in control engineering.

Figure 1 .
Figure 1.Schematic overview of the UR10 robot: CAD model.

Figure 4 .
Figure 4. Schematic overview of a classical robot joint drivetrain.

Figure 5 .
Figure 5.The structure of the nested control architecture.

Figure 6 .
Figure 6.The basic scheme for determining the frequency response slope.

Figure 7 .
Figure 7.The robust scheme for determining the frequency response slope: (a) intermediary solution; (b) final solution.

Figure 8 .
Figure 8. Simplified robust scheme for determining the frequency response slope.

Figure 9 .
Figure 9. Closed loop scheme for determining the frequency response of a process and the frequency response slope.

Figure 10 .
Figure 10.Simplified closed loop scheme for determining the frequency response of S(s)and the corresponding frequency response slope.

Figure 11 .
Figure 11.The loop frequency response (blue line) and the 'forbidden region' (red circle).GM, gain margin; PM, phase margin.

Figure 12 .
Figure 12.Reference and error signals measured experimentally.

Figure 13 .
Figure 13.Reference and steady state error signals used to determine the robot frequency response and its corresponding slope.

∞ 0 e 1 ∑Figure 14 .
Figure 14.Transient component of the error signals used to determine the robot frequency response and its corresponding slope.

Figure 15 .
Figure 15.The UR10 robot available in our laboratory.

Figure 16 .
Figure 16.Communication model of the data transmission between the UR10 robot and the MATLAB R software.

Figure 17 .
Figure 17.Visualization of the ROS computation graph.

Figure 18 .
Figure 18.The flowchart for controlling the movements of the robot.

Figure 19 .
Figure 19.Step response of the UR10 robot for Joint 2(PID, the designed controller; de f , the default controller available on the robot); (a) first robot configuration and (b) second robot configuration.

Figure 20 .
Figure 20.Step response of the UR10 robot for Joint 3; (a) first robot configuration and (b) second robot configuration.

Figure 21 .
Figure 21.Response of the UR10 robot for Joint 2.

Figure 22 .
Figure 22.Response of the UR10 robot for Joint 3.