Model Predictive Control With Integral Compensation for Motion Control of Robot Manipulator in Joint and Task Spaces

The motion control of robot manipulators is a crucial problem concerning automatically controlled robots. In this work, the model predictive control method with an integral compensation (MPC-I), which compensates for the matched uncertainties due to unmodeled dynamics, is proposed to solve the trajectory tracking problem of robot manipulators in joint and task spaces. First, this paper decouples the joint variables of the robot manipulator using a computed torque control method. The MPC-I method is, thereafter, derived to realize the motion control of the robot manipulators in joint space. To realize the motion control of the robot manipulator in task space, the task space is, thereafter, converted into the joint space, in which the MPC-I method is executed, afterward, to control the robot. Furthermore, an MPC-I variation, in which the inverse kinematics is calculated indirectly, is proposed to achieve the motion control in task space. The novelty of this paper is to propose the MPC-I method and the method of converting task space to joint space with indirect inverse kinematics calculation. The former is suitable for the dynamic control of the robot manipulators in the joint space, and the latter can extend the MPC-I method to dynamic control in task space. To evaluate the performance of the proposed control method, motion control simulations are performed in the task and joint spaces, respectively. Simulation results and comparisons verify the effectiveness of the proposed control approach for the dynamic control of the UR5 robot manipulator.


I. INTRODUCTION
Due to existing demands to improve reliability, accuracy, repeatability, and productivity of industrial processes, the significance of robot technology has been highlighted [1]. Furthermore, robots are being used for different purposes, such as surgery [2], handling hazardous materials [3], and motion assistance [4], and in different environments, such as space [5], underwater [6] etc. In all these cases, robot manipulators are employed; therefore, understanding their motion control is crucial. The conventional control methods, such as proportional integral derivative (PID) [7] and linear quadratic regulator [8], can guarantee the reliability and safety of the system, realizing a stable and effective motion control of The associate editor coordinating the review of this manuscript and approving it for publication was Guilin Yang . the manipulator. Their main advantage is their simplicity, as linear controllers are perfectly suited for numerous applications. However, some disadvantages, such as cumbersome gain tuning and adjustment requirements for implementation on a robot manipulator, are still not addressed. Furthermore, obtaining optimal parameters of classical PID controllers may be difficult, especially for systems comprising of high nonlinearities and couplings [9]. As a result, in the past few decades, advanced control technologies [10]- [16] have attracted more attention. The control techniques like Adaptive Control [10], [13], Sliding Mode Control [11], Robust Control [14]- [16] are very promising for implementation on a robot manipulator. Among them, the most representative [15] dealt with the robust task space control of electrically driven robot manipulators using voltage control strategy, which is simpler, less computational and requires less feedbacks than VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ the conventional robust controllers. However, model predictive control (MPC) is not only robust enough to eliminate the effects of nonlinear behavior, complexity and coupling effects in the dynamics, but also complies to possibly enforced constraints on the state and input variables, which is what the above control methods are not good at. The MPC, also named receding horizon control, is recognized as the most efficient and application potential method, considering the advanced process control methodology [17]. In MPC, a process model is used to predict and optimize the future behavior of the system [18]. Future control inputs in certain control horizon and future plant responses in certain prediction horizon are both predicted using a system model and optimized at regular intervals with respect to a performance index. MPC has generated huge economic profits on thousands of industrial control systems over the world [19]. This kind of control has been first employed in large chemical factories; however, in recent years, it has been widely promoted in other industrial fields. For example, Paluszczyszyn et al. [20] has recently proposed this method to adjust the patient support system position. Meanwhile, Tadokoro et al. [21] uses this technology to adjust the motion control for rigid-body dynamical systems, such as spacecraft or aerial vehicles. Furthermore, other interesting results on MPC can be found in [22]- [26].
In recent years, MPC application in robot manipulator motion control has also been very popular; however, additional difficulties due to the nonlinearities of the manipulator occur. The most common solution of the problem is the linearization method, as illustrated in [27] and [28], in which the feedback of the inverse dynamics of the manipulator and the Taylor series are, respectively, used in the linearization. Furthermore, boundedness and stability convergence of the state tracking have also been proved in [28]. However, the common disadvantage between [27] and [28] is that compensation measures after linearization were not considered. The other solution is to combine MPC with other advanced control methods to realize the motion control of the manipulator. Reference [29] proposes a combination of MPC and H-infinity control, which is also a nonlinear design technique for a robot manipulator. The result of this combination, however, could not combine the respective advantages of these schemes while trying to avoid their shortcomings. Meanwhile, a robust hierarchical multiloop control scheme design to solve motion control problems for robot manipulators was proposed in [30]. The inverse dynamics-based feedback linearized robotic system and the combination of MPC with an integral sliding mode controller are the key elements of this control approach. The internal control loop was used to compensate for the matched uncertainties due to unmodeled dynamics. This method was verified and validated by simulations; however, the predictive control method is only applicable to the motion control in joint space but not to the motion control in task space that is the workspace of the end-effector.
Based on prior research, MPC variation methods that can be applied to both joint and task spaces of robot manipulators have been rarely found. As a result, this paper proposes an model predictive control with integral compensation (MPC-I) for the motion control of the robot manipulator that can be used in the joint and task spaces. Motion control in joint space was first decoupled by a computed torque control method [31]. Afterward, an approximate prediction was performed using Taylor expansion. Due to the truncation error of the Taylor expansion, an integral compensation is introduced to realize the motion control. To realize motion control in task space, the task space is first transformed into the joint space and, thereafter, the MPC-I in joint space is executed. The task space can be transformed into the joint space by two different procedures. In the first one, the inverse kinematics are calculated directly, whereas, in the second one, which is to calculate inverse kinematics indirectly. The latter has a low computational complexity and is not affected by the configuration of the robot manipulator. Both procedures were compared in this work with simulations.
The rest of the paper is organized in the following manner. In Section II, preliminary knowledge, including the robotic dynamics in joint space and exponential coordinate representation of rigid-body motions are discussed. In Section III, the MPC-I methods that realizes the motion control of the robot manipulator in the joint and task spaces are introduced. In Section IV, the simulation results, which were obtained from the co-simulation of the UR5 manipulator on MATLAB and V-rep [32] platforms, are exhibited. Finally, in Section V, conclusions and future work are informed.

II. PRELIMINARIES
In this section, the preliminary knowledge used in this article was introduced. First, the robotic dynamics in joint space was described. The exponential coordinate representation of rigid-body motions was, thereafter, described, which lays the foundation for the conversion of the task space to the joint space in the Section III-B2, and facilitates the definition of rotation and translation errors in the Section IV-B.

A. THE ROBOTIC DYNAMICS IN JOINT SPACE
The dynamic equation for n degree of freedom robot manipulators can be written as: where M(θ ) ∈ R n×n is the positive definite mass matrix, C(θ,θ) ∈ R n×n is called the Coriolis matrix, g(θ) ∈ R n×1 is the vector specifying the effects due to gravity, τ ∈ R n×1 is the vector of torques (forces) of each joint actuator, and θ = [θ 1 θ 2 . . . θ n ] T ∈ R n×1 is the vector of joint variables. Furthermore,θ can be expressed as the following to calculate the acceleration of the joint variables: where h(θ,θ) is defined as:

B. EXPONENTIAL COORDINATE REPRESENTATION OF RIGID-BODY MOTIONS
The product-of-exponentials formula used to describe the kinematics of robot manipulators is introduced by [33]. Meanwhile, the Chasles-Mozzi theorem states that every rigid-body displacement can be expressed as a displacement along a fixed screw axis in space [34]. By combining this theorem with the product of the exponential formula, a six-dimensional exponential coordinate (S) of a homogeneous transformation (T h ) can be defined as: where where S a and S p represent the rotation and translation, respectively. The matrix are therefore defined by: where T h belongs to Special Euclidean Group, named SE(3), and [S], which belongs to se(3), the Lie algebra of the Lie group SE (3), is defined as: The specific calculation process of (7) and (8) is described by [33].

III. MOTION CONTROL WITH TORQUE OR FORCE INPUTS
Due to the nonlinearities of the robot manipulator dynamics, motion control with torque or force inputs is more difficult than motion control with velocity inputs. However, the latter is generally limited to applications with low or predictable force/torque requirements, the former is suitable for more applications. In this section, therefore, a control method that generates joint torques or forces to try to track the desired trajectory in joint and task spaces is discussed.

A. MPC-I METHOD OF ROBOT MANIPULATORS IN JOINT SPACE
In this section, the MPC-I method for the motion control of the robot manipulators in joint space is introduced, and the prediction equation required by this method is easily expressed in this space. Trajectories are also naturally described by the joint variables, and there are no issues of singularities or redundancy [33]. First, the joint variables are decoupled by a computed torque method. Afterward, an equation that predicts the behavior of the joint variables is obtained. The objective function is, thereafter, calculated, and the optimization problem is solved. Moreover, the stability analysis of the control method is performed. Finally, integral compensation is introduced to eliminate the steady state error.
The following assumptions were made before calculating the appropriate solution: 1) The MPC controller used in this paper is based on the solution of the so-called finite horizon optimal control problem. Therefore, the prediction and control horizon were defined within p and m steps (where p ≥ m), respectively, and the time interval of each step is t. 2) For discrete systems, the acceleration of joint variables between steps k and k + 1 was considered constant. 3) Beyond the control horizon, the control quantity was preserved.
where k + i|k represents the prediction of the k + i step at k step.

1) PREDICTION OF JOINT VARIABLES
According to these assumptions, the prediction formula of joint variables could be written as where where I and 0 denote the (n × n) identity matrix and (n × n) null matrix, respectively, and the matrix subscript represents the number of elements rather than the dimension of the matrix. VOLUME 8, 2020 According to (2),θ is coupled with τ . These variables were decoupled with a computed torque control method. As a result, the torque or force at each step was defined as: By combining (16) and (2), the following expression can be obtained:θ As a result, (10) can be written as: where The optimization problem can be described as where the objective function can be defined as where nonnegative weighting matrix and the expected target joint vector were evaluated by equations (22), (23), and (24), respectively.
As a result, the solution of (20) was obtained: where H and F(k + 1|k) were defined by (26) and (27), respectively. According to the MPC method, the first element of W (k) * should be applied to the system: where E θ (k + 1|k) is defined by (29) and where As a result, according to (16), the control torque/force was where w(k) * and h(k) are defined as (33) and (3), respectively.

3) STABILITY ANALYSIS OF THE CONTROL METHOD IN JOINT SPACE
This section discusses the stability of the methods proposed above. The authors adopted a similar concept as reported by [28] to conduct this analysis.
Using (36) and (29) in (2), the acceleration of the joint variables could be calculated as: where d (k + 1) can be written as , (38) where I a , I b , and A are defined as (13), (14), and (15), respectively, and O( t 3 ) is truncation error. Furthermore, the joint variable error θ e is defined as Therefore, the joint variable acceleration error (θ e ) can be calculated as ). (40) By combining (38) with (40), the joint variable error dynamic could be obtained: where e 1 = θ e and e 2 =θ e . Moreover, the error state (e) can be defined as: As a result, the tracking problem of θ e andθ e is reduced to the regulation problem of e that was treated as a linear system with non vanishing perturbation [35], which are described by the following equations: where the perturbed term is . (44) The origin may not be an equilibrium point of system (41); therefore, the stability of the system at the origin was not analyzed. Furthermore, the solution of the perturbed system may not approach the origin when t → ∞. However, the boundedness of the system (41) could be analyzed. Lyapunov function candidate is defined as where both QK mpc I a and Q are positive definite matrixes. By deriving the Lyapunov function concerning time, (45) was transformed into: After combining (41) with (46) the following expression was obtained: There exist an lower bound δ 1 (θ d ), δ 2 ( t), and positive reals t p1 , t p2 for all e 2 : The time derivative of V is bounded as follows: As a result, if λ min (QK mpc I b ) + (t p1 + t p2 )λ min (Q) > 0, the solutions of the system expressed by (41) are uniformly bounded. The dynamic of the tracking error (41) is specified by the matrix Q and parameters t,θ d . Because this bound could be made small by reducing the penalty of the control signal w → 0 in (21). With w = 0, that will lead to As a result, t p1 = 0 and the truncation error O( t 3 ) can be eliminated by introducing integral compensation, which will be introduced in the next section. Considering this equation, the origin is, therefore, the equilibrium point of the system (41) and the time derivative of the Lyapunov function could be modified to: As QK mpc I b is positive definite matrix, the equilibrium point e o = [0 0] T is globally asymptotically stable, according to La Salle's invariant theorem [35].

4) INTEGRAL COMPENSATION
In the previous sections, the MPC method was introduced and its stability was proved. It is obvious that the truncation error is caused by the matched uncertainties due to unmodeled dynamics. To eliminate the truncation error, this section demonstrates the necessity of introducing integral compensation in the previously proposed MPC method.
To introduce the integral compensation, (29) is modified as follow, where the integral gain K I is a positive definite matrix: and e (i) can be calculated as Therefore, theθ e can be calculated as Then we calculated the derivative of (56) concerning time, the joint variable error dynamic could be obtained: where e 3 =θ e . Moreover, the error state with integral compensation (e I ) can be defined as: Then (43) can be modified aṡ where the perturbed term is The proof of the stability of (59) is similar to the method in Section III-A3, which will not be repeated here. Below we will mainly analyze the perturbation term (60). Whenθ d is a constant, we can get ... θ d = 0. Whenθ d is a variable, we can set w = 0, that will lead to K mpc A − I = 0.
t is very small, O( t 4 ) in (60) is much smaller than O( t 3 ) in (44), so O( t 4 ) can even be ignored. Therefore, the integral VOLUME 8, 2020 compensation makes the perturbed term (60) much smaller than (44), which effectively eliminates the truncation error.
Therefore, the final equation that represents the MPC-I method can be obtained as: where h(k) is defined as (3) and where K mpc and E * θ (k + 1|k) are defined as (34) and (53), respectively.

5) CONTROL METHOD OF MOTION IN JOINT SPACE
In this section, the detailed process of the MPC-I method to realize the motion control of the robot manipulators in joint space is described.
The MPC-I controller block diagram is illustrated in Fig. 1. It should be noted that K is to transform θ e into e . Furthermore, the detailed control process is shown in Algorithm 1. In line 1, the integral compensation S θ e is the integral of joint variables errors. In line 2, N is the maximum number of control steps. In lines 5 and 6, the joint variables error θ e (k) and e (k) are calculated as (39) and (55), respectively. From line 8 to 10, E * θ (k + 1|k), K mpc , and w * I (k) are calculated as (53), (34), and (62), respectively. In line 11, the dynamics parameters h(k) is calculated as (3), but in order to reduce the inaccuracy of the estimation, the θ andθ used to calculate M(k) and h(k) are replaced by the follows In line 12, the control torque/force τ * I (k) is calculated as (61).

B. MPC-I METHOD OF ROBOT MANIPULATORS IN TASK SPACE
As the robot interacts with the external environment and multiple objects, expressing the motion of the robot as a trajectory of the end-effector in task space is more convenient.
To perform this task, the trajectory in task space must be converted to joint space.
Provided the corresponding trajectory in joint space is feasible, this conversion may be done by direct or indirect inverse kinematics calculation. Afterward, the control method in Section III-A5 be used to realize the motion of the robot.

Algorithm 1
The MPC-I Method to Control the Robot Manipulator in Joint Space 1: Initialization: Determine the time step t, the prediction horizon p, the control horizon m, the nonnegative weighting matrix θ , w , K I , and set the integral compensation S θ e = 0 2: for each k ∈ [0, N ] do 3: Read the current the robot manipulator states θ(k) andθ(k) 4: Determine the expected target joint variables vector d (k + 1) 5: Calculate the joint variables error θ e (k) 6: Determine e (k) 7: Calculate the error integral S θ e = S θ e + e (k) t 8: Determine E * θ (k + 1|k) 9: Calculate K mpc 10: Determine w * I (k) 11: Calculate the dynamics parameters M(k) and h(k) 12: Determine the control torque/force τ * I (k) and act on the robot manipulator 13: end for

1) THE CONTROL METHOD OF DIRECTLY CALCULATING INVERSE KINEMATICS
The end-effector trajectory is specified by T d , where T d ∈ SE (3). The inverse kinematics can transform the end-effector trajectory into a trajectory based on joint variables by the following expression, where f −1 represents the inverse kinematics calculation. Then proceed with control as in Section III-A5. The MPC-I method in which directly calculates the inverse kinematics (MPC-I_IK), not only needs significant computing power to calculate the inverse kinematics, but also is affected by the configuration of the robot manipulator.

2) THE CONTROL METHOD INDIRECTLY CALCULATING INVERSE KINEMATICS
To overcome the shortcomings of the method described in the previous section, the MPC-I method in which the inverse kinematics are indirectly calculated (MPC-I_NIK) is introduced in this section.
By assuming that the end-effector frame can be expressed with a coordinate vector (x) governed by the kinematics x d = g(θ d ), the forward kinematics can be expressed as the Taylor expansion where J(θ d (k)) is the coordinate Jacobian evaluated at θ d (k). After performing the first-order truncation, (66) can be approximated as: (67) where the J(θ d (k)) † denotes the pseudo-inverse of the Jacobian.
To modify (67) to work with a desired end-effector configuration represented as T d ∈ SE(3) instead of as a coordinate vector x d , the coordinate Jacobian J was replaced with the end-effector body Jacobian J b ∈ R 6×n and x d (k + i) − g(θ d (k + i − 1)) was replaced with a body twist V b , which would generate the motion from T d (k + i − 1) to the desired configuration T d (k + i). Furthermore, V b (k + i − 1) is determined by the matrix logarithm introduced in Section II-B: As a result, (67) can be modified to: The method proposed in this section is fully described by Algorithm 2. Lines 1 to 3 of Algorithm 2 are similar to Algorithm 1. From line 4 to 12, the methods to convert the desired pose of the task space T d into joint variables θ d without calculating kinematics are used. In lines 8, 9, and 12, the body twist V b , the expected joint variables of the next step θ d (k + i), and the expected target joint vector are calculated with (68), (69), and (24), respectively. In line 13, the expected joint variables at the current moment θ d (k) are approximately calculated with (70) and (71), as a preliminary step for the calculation of the integral error. As a result, the estimated value of the expected joint variable for step k + 1 at step k is the expected joint variable for the next step. From line 14 to 20, the MPC-I method of a robot manipulator in joint space is executed.

IV. SIMULATION AND DISCUSSION
To illustrate the effectiveness of the MPC-I method proposed in this article, the motion control of the robot manipulators in joint and task spaces is individually simulated. Considering a UR5 robot manipulator as an example, the co-simulation is conducted on the MATLAB 2015b and V-rep 4.0 platforms. The product-of-exponentials parameters of the UR5 with the current configurations relative to the base frame are shown in Table 1, where [ω T i ν T i ] T are the screw axes expressed in the fixed space frame. T h0 is the end-effector configuration when the UR5 is at its home position. All simulations are performed on a desktop computer (Core i5 3.40GHz, 16GB RAM, MATLAB 2019b and V-rep 4.0 software platforms).

A. COMPARISON OF MOTION CONTROL BETWEEN MPC AND MPC-I METHODS IN JOINT SPACE
To prove the benefits of the integral compensation, the trajectory tracking of the UR5 robot manipulator in joint space was Algorithm 2 The MPC-I Method to Control the Robot Manipulator in Task Space 1: Initialization: Determine the time step t, the prediction horizon p, the control horizon m, the nonnegative weighting matrix θ , w , K I , and set the integral compensation S θ e = 0 and θ d0 = θ(0) 2: for each k ∈ [0, N ] do 3: Read the current robot manipulator states θ(k) andθ(k) 4: Calculate the pose transformation matrix T h (k) using forward kinematics and set T d (k) = T h (k) 5: Determine the end-effector body Jacobian J b (k|k) 6: for each i ∈ [1, p] do 7: Read expected pose transformation matrix T d (k +i) 8: Calculate the body twist Determine the expected joint variables θ d (k + i) 10: Calculate the end-effector body Jacobian J b (k +i|k) 11: end for 12: Determine the expected target joint vector d (k + 1) 13: Calculate the expected joint variables at the current moment then set 14: Determine e (k) 15: Calculate the error integral S θ e = S θ e + e (k) t 16: Determine E * θ (k + 1|k) 17: Calculate K mpc 18: Determine w * I (k) 19: Calculate the dynamics parameters M(k) and h(k) 20: Determine the control torque/force τ * I (k) and act on the robot manipulator 21: end for  implemented. The UR5 robot manipulator will take 8s to go from the initial angle θ init to the end angle θ end . The specific values of the initial and the end angle of UR5 are shown in Table 2. The expected motion trajectory is, thereafter, VOLUME 8, 2020  obtained by the fifth-order polynomial interpolation of θ init and θ end . Furthermore, the time interval was defined as t = 0.025s, and the prediction and control horizons were defined as p = m = 5. The other parameters in the MPC and MPC-I motion control methods were θ = diag(150I 6 , 120I 6 , 90I 6 , 60I 6 , 30I 6 ), (72) w = diag(0.5I 6 , 0.4I 6 , 0.3I 6 , 0.2I 6 , 0.1I 6 ), (73) K I = 0.003 diag(I 6 , I 6 , I 6 , I 6 , I 6 where I 6 6 denotes the (6 × 6) identity matrix. The actual and expected joint angles for motion control of the UR5, calculated by the MPC and MPC-I methods, versus time are shown in Fig. 2. These six figures respectively correspond to the six joint angles of UR5. In addition, the joint input torque for motion control of the UR5, calculated by the MPC and MPC-I methods, versus time are shown in Fig. 3. The differences between MPC and MPC-I methods are shown by driven motions. It can be seen from the Fig. 2 and 3 that the six joints controlled by the MPC and MPC-I methods all follow very accurately and the fluctuation of joint input torque of the MPC-I method is more serious than the MPC method, which can be relieved by adjusting K I .
However, to further understand the behavior of these steady-state errors, the joint errors of the UR5 controlled by the MPC and MPC-I methods versus time are shown in Fig. 4. The steady-state errors, especially those from the second, third, and fourth joints, were significant when the six joints were controlled by this MPC method, which are 0.0172 • , −0.1183 • , and −0.0323 • respectively. Whereas,  for the MPC-I method, these values converged to 0.0002 • , 0.0009 • , and 0.0011 • , considering the same simulation time. The errors of the MPC-I method were acceptable and were expected to further decrease over time, as analyzed in Section III-A4. To further illustrate the advantages of the MPC-I method, the objective function values (Equation (21)), calculated by the MPC and MPC-I methods, versus time is shown in Fig. 5. As a result, this simulation proves the advantage of the MPC-I method, in which an integral compensation was introduced.

B. COMPARISON OF BETWEEN THE MPC-I_NIK AND MPC-I_IK METHODS IN TASK SPACE
In this section, the reliability of the MPC-I_NIK is verified by comparing its results with those of an MPC-I_IK. The inverse kinematics was calculated directly by the method proposed by [36]. In this experiment, a trajectory that will allow the UR5 to draw the letter ''F'' in 20s is expected to be obtained. The UR5 model of the V-rep simulation platform is shown in Fig. 6(a) and the desired configuration trajectory of the end-effector of the UR5 is shown in Fig. 6(b).
The initial configuration (T init ) and end configuration (T end ) of end-effector are defined as: To illustrate the configuration trajectory and error analysis of the end-effector, the matrix logarithm was used to convert the configuration to six-dimensional exponential coordinates (S) according to Section II-B. When the expected (S d ) and actual (S c ) exponential coordinates were previously determined, as indicated in (77) and (78), respectively, the rotation (E a ) and translation (E p ) errors in task space are defined by (79) and (80), respectively.
The six components of S d and S c versus time are shown in Fig. 7. The actual tracking trajectory was separately obtained by the MPC-I_NIK and MPC-I_IK methods. From  Fig. 7, the two methods are satisfactory for the accuracy of the trajectory tracking. For further analysis, the rotation and translation errors of the trajectory tracking expressed by the exponential coordinate versus time are shown in Fig. 8. Both methods were able to obtain results whose rotation and   translation errors, in exponential coordinates, were within an acceptable range (order of magnitude of 10 −6 ). To evaluate the accuracy of the space conversion method proposed in this paper, the joint angle estimated by (69) was used to calculate the configuration of the end-effector through forward kinematics. The configuration of the end-effector was, thereafter,  converted into exponential coordinates and compared with the exponential coordinates of the desired pose. The rotation and translation errors of the spatial transformation are illustrated in Fig. 9. The rotation error was within 1 × 10 −5 and the translation error was within 3 × 10 −5 . Finally, to illustrate the computational complexity of the MPC-I_NIK method proposed in this paper, the calculation time of the inverse kinematics of methods MPC-I_NIK and MPC-I_IK is shown in Fig. 10. From Fig. 10(a), the advantages of the MPC-I_NIK method to the MPC-I_IK method can be easily verified. From Fig. 10(b), the medians of the time for calculating inverse kinematics for methods MPC-I_NIK and MPC-I_IK are 0.0014s and 0.0132s, respectively. It can be concluded that the former has a great advantage over the latter in terms of computational complexity.

V. CONCLUSION
This article proposes an MPC-I method for motion control of a robot manipulator considering the joint space and analyzes the advantages of introducing an integral control term. To adapt the method to enable the motion control in task space, a space conversion method for indirect calculation of inverse kinematics was proposed. Two simulations were conducted to illustrate the effectiveness of the proposed method. In the first one, the effect of the MPC and MPC-I methods was compared, considering the motion control in joint space. The results indicated that the steady-state errors of the six joints controlled by the MPC method were 0.0008 • , 0.0172 • , −0.1183 • , −0.0323 • , −0.0042 • , and 0.0022 • , respectively; however, with the proposed MPC-I method, the steady-state errors converged to lower values (−0.0003 • , 0.0002 • , 0.0009 • , 0.0011 • , 0.0002 • , and −0.0005 • , respectively), at the same simulation time. In the second simulation, the MPC-I_NIK and MPC-I_IK methods were compared considering the motion control in task space. Both methods could achieve results with rotation and translation errors in exponential coordinates within the order of magnitude of 10 −6 . Meanwhile, the rotation and translation errors of the space conversion method proposed in this paper were 1 × 10 −5 and 3 × 10 −5 , respectively, which shows that the method is effective. And the median time of the MPC-I_NIK method to calculate the inverse kinematics is 0.0014s, which is a great improvement over the calculation time of MPC-I_IK method.
In future works, actual constraints, such as maximum speed and maximum acceleration, will be added to the MPC-I method, resulting in a more practical method. Furthermore, the motion control method in joint and task spaces proposed in this paper will be applied to practical applications, such as handling of hazardous materials and motion assistance. Moreover, the MPC-I method will be extended to hybrid motion-force control. QINGSHENG LUO received the B.S. degree in precision instrument design, the M.S. degree in mechanical optimization design, and the Ph.D. degree in mechanical engineering from the Beijing Institute of Technology, Beijing, China, in 1982China, in , 1987, and 2001, respectively. His research interests include specialized robot technology, mechatronical product innovative design, mechatronical system control technology, and mechatronical device integration technology.
LIJUN QIAO received the B.S. degree in measurement and control technology and instrument from the Taiyuan University of Science and Technology, Taiyuan, China, in 2014, and the M.S. degree in instrument science and technology from the North University of China, Taiyuan, in 2017. He is currently pursuing the Ph.D. degree in mechanical engineering with the Beijing Institute of Technology. His major research interests include dual-arm robot, robot dynamics control, and robot stability. VOLUME 8, 2020