Load Parameter Identification for Parallel Robot Manipulator Based on Extended Kalman Filter

,


Introduction
Parallel robot manipulator (PRM) entails the advantages of higher precision, faster response, higher rigidity, and stronger carrying capacity over serial robots and, hence, is widely applied in many fields of industry [1,2]. Currently, the standard industrial control technique applied in PRMs is PID control, which neglected the complex dynamic characteristics of robots, so that it can hardly meet the requirement of fast and accurate motion [3,4]. During the past decade, several model-based control methods, such as computed torque control [5], dynamic feedforward control [6], and modal control [7,8], are used in controlling PRMs because of the excellent control accuracy and dynamic performance of these methods. However, these methods are all based on the dynamic model to design, which makes the control performance of these methods strongly depend on the accurate knowledge of the dynamic model.
In practical applications, PRMs always have to hold an uncertain load to work, and the dynamic parameters of load are hard to know a priori and cannot be measured directly, which makes several model-based control schemes unable to be applied [9]. Take modal control as an example, the key of the control strategy is the modal conversion matrix. e calculation of modal conversion needs the system mass matrix, but the uncertainty of the load will make the accurate system mass matrix impossible to obtain [10,11]. Otherwise, the uncertain load, which is a main external disturbance of the system [12], has a significant impact on system dynamic performance; it will cause dynamic coupling among different degrees of freedom (DOFs) [13,14].
us, in order to achieve the high-performance control of the system, the uncertain load disturbance is necessary to be estimated and compensated.
In parameter identification, three kinds of methods have been proposed to estimate the robot dynamic parameters. e first method is using computer-aided design (CAD) techniques to obtain the dynamic parameters of robots. e 3-dimensional (3D) model of robots generally provided by the robot manufacturer and the parameters, such as the inertia tensor and centroid position, can be solved by any CAD software based on these models. However, the parameters obtained from the CAD techniques are not identical to the real robot because of the manufacturing and assembling error. In addition, for the PRMs with uncertain load, it is difficult to build a 3D model for load, which is usually not provided by the manufacturer and is made up of many complex parts. To avoid this problem, the method of the physical experiment is used to perform the determination of the dynamic parameter. is method can better estimate the mass, centroid position, and inertia tensor of the part. However, the necessity of disassembling the robot and measuring by special devices limits the application of this method. Moreover, this method is inappropriate for the measurement of large-sized parts, such as aircraft and submarine, which are common types of loads and are generally too big to be measured and too complex to be disassembled. e last method is the theoretical identification method, which can obtain better identification result, and does not require disassembly of robots and special devices for measurement compared to the above two methods. Until now, many identification algorithms have been presented to estimate the parameters of robots [15], such as the least squares method [16,17], Kalman filtering method [18][19][20], maximum-likelihood method [21,22]. In addition, among these identification algorithms, the most used method is the least squares method [23]. However, most of the research objects of these methods are serial manipulators, while for parallel manipulators, there is little research at present.
Compared to the serial manipulator, the research on parameter identification of PRMs started late, and too little work has been devoted to the load parameters identification. Chen derived the estimation equation in a linear form of identified parameters based on a new structured Boltzmann-Hamel-d'Alembert approach and used the least square method to identify the parameters [24]. Tian proposed an inertial parameter identification method based on sinusoidal vibrations of a six-degree-of-freedom parallel manipulator and used the least square method to identify the parameters [25]. Briot and Gautier used total least squares to identify the parallel robot dynamic parameters [26]. Wu et al. investigated the dynamic parameter identification of a redundantly actuated parallel manipulator and proposed a two-step identification approach based on the least squares method to identify the dynamic parameters of the system [27]. anh used the direct pattern search technique to do the dynamics identification for a redundant 3-(P) RRR manipulator [28].
Normally, the methods to identify the parameters of the PRM adopt the methods of the least squares method and the weighted least squares method. However, the least square method is not suitable for the identification of parallel manipulators. e main reasons are as follows: (1) e LMS needs to establish an inverse dynamic model that is linear with respect to the dynamic parameters. However, due to the complex and coupled dynamics of the parallel manipulators, it is still not easy to rewrite the dynamic equation into a linear form that is suitable for using the parameter identification algorithm [29]. (2) e LSM is sensitive to measurement noise. However, parallel manipulators are generally driven by hydraulic pressure, and the noise of driving force, displacement, speed, and acceleration is relatively large, which seriously limits the identification accuracy and convergence speed of the method [30,31]. (3) For the LSM, the observation matrix in the inverse dynamic model requires the value of displacement, speed, and acceleration in the task space. Since the displacement and speed of each leg of the parallel mechanism are relatively easy to measure, the displacement in the task space can be obtained by the positive kinematics solution, and the speed in the task space can be solved by the speed Jacobian matrix. But, for acceleration, it is not easy to be directly measured.
is paper proposed a method to identify the load of a parallel robot manipulator. Compared with the traditional least square method, the proposed identification approach does not require linearization of the dynamic model and optimization of the excitation trajectory, which are two complex problems to solve. Under the specific signals, the dynamic equation derived by the Newton-Euler method can be simplified and the simplified dynamic equation has two advantages: (1) this equation is weakly nonlinear with respect to the dynamic parameters of load and (2) the parameters to be identified are independent of each other in the equation and there is no product form. According to the simplified model, we designed the excitation trajectory, which has a simple form and can excite the dynamic parameters of the load. Importantly, the EKF algorithm is applied to estimate the load parameters, which is not sensitive to measurement noise and without acceleration measurements in the process of identification. e organization of this paper is as follows. In Section 2, the dynamic model of PRM with uncertain load is established and the influence of load disturbance for the system is also analyzed. Section 3 presents the design of excitation trajectory and the process of load parameters identification based on the EKF algorithm. Section 4 shows the numerical simulation for verifying the proposed method. Finally, the main conclusions of this paper are presented in Section 5.

Mathematical Modeling
e parallel manipulator studied in this paper is a six-degree-of-freedom Gough-Stewart parallel manipulator, as shown in Figure 1(a). is manipulator is mainly composed of two platforms and six driving legs. e lower platform fixed on the ground is named static platform and the upper platform used to carry loads is named moving platform. e driving leg is the actuator to drive the moving platform to realize translation and rotation. In order to facilitate the 2 Complexity analysis of the parallel manipulator, the coordinate system is established and shown in Figure 1

Kinematics Modeling.
According to the space geometry theory, the length vector of the six legs can be expressed by where A is the coordinate matrix of six lower hinges in the static coordinate. B is the coordinate matrix of six upper hinges in the moving coordinate. R denotes the transformation matrix from the static coordinate to the moving coordinate. p � x y z T is a translation displacement vector. x, y and z are the displacement of platform mass center along the x-axis, y-axis, and z-axis, respectively. p 0 is the initial height matrix composed of six initial height vectors, and each vector is 0 0 h 0 T . h 0 is the initial height when the upper platform motion table in the middle position.
By differentiating equation (1), one obtains e matrix form of equation (2) is where J v is velocity Jacobian matrix between the velocity of the upper platform and the leg velocity. L n is the unit direction vector matrix of the six legs' direction and T is the angular velocity of the moving platform and _ L is the leg velocity. Importantly, matrix J is the key to derive the dynamic equation of joint space.
In addition, kinematic analysis of the load is also required, such as centroid position and acceleration. Let the load eccentric position along three axes of the moving coordinate system be Δx, Δy, and Δz. e load eccentric position along three axes of the static coordinate can be expressed as l x � c q 5 c q 6 Δx + s q 4 s q 5 c q 6 c q 4 s q 6 Δy + c q 4 s q 5 c q 6 + s q 4 s q 6 Δz, l y � c q 5 s q 6 Δx + s q 4 s q 5 s q 6 + c q 4 c q 6 Δy + c q 4 s q 5 s q 6 − s q 4 c q 6 Δz, where c () and s () are abbreviations for the trigonometric functions cosine and sine; q 4 , q 5 , and q 6 are the three Euler angles of the system; and l x , l y , and l z are the load eccentric position along three main axes. From equation (4), it is easy to know that the l x , l y , and l z are related to the Euler anglers of the upper platform.

Complexity
Let the load centroid acceleration in the static coordinate be € p l � € x a € y a € z a T . According to the rigid body dynamics theory, the relation between the load centroid acceleration and the upper platform centroid acceleration can be expressed by e relationship between € p l and € p is shown in Figure 2. Note that the O C − X C Y C Z C coordinate system in Figure 1 is not the same as the O B − X B Y B Z B coordinate system in Figure 2. e origin of the O C − X C Y C Z C coordinate system is always coincident with the origin of the moving coordinate system, and the axes are always in the same direction as the axes in the static coordinate system. O C − X C Y C Z C coordinate system is the coordinate system used to establish the dynamic equation. In addition, G p is the position of the load center of mass.

Dynamic
Modeling. First, we analyze the load dynamics. Based on Newton's second law and angular momentum theory, inertia force and moment of load can be described: where m l is the mass of load; B I L is the inertia tensor of load in the static coordinate system; and Q is the antisymmetric matrix composed of each element of the eccentric vector of the center of mass loaded in the static coordinate system and is as follows: Based on the Newton-Euler method, the dynamic model in task space can be obtained: where f i is the force of hydraulic cylinder of ith leg; e i is the unit direction vector of the ith legs; f L is inertia force of load; T L is inertia moment of load; m p is the mass of moving platform; and B I P is the inertia tensor of moving platform in the static coordinate system. In addition, the third term on the right of equation (9) is the gravity of the system. e second term on the right of equation (9) represents the inertia force of the moving platform and the second term on the right of equation (10) represents the inertia force of the moving platform. Moreover, the third term on the right of equation (10) is Coriolis/centrifugal forces.
From equations (9) and (10), it can be obtained that where τ is the driving force vector in the task space; H( _ q, q) is the centrifugal/Coriolis term; G is the gravity term; and M(q) is the mass of the system. e relationship between the driving force vector in the task space and the driving force vector in the joint space can be expressed by From equations (11) and (12), the dynamic model of the robot in joint space can be rewritten as where M L (l) is the system mass matrix in joint space, H L (l, _ l) is the centrifugal/Coriolis term in joint space, G L is the gravity term in joint space, and e impact of the centrifugal/Coriolis term is small, which can be ignored, and equation (13) can be rewritten as

Dynamic Effects of Load.
Before establishing the identification model and designing the excitation trajectory, it is necessary to analyze the impact of the load on the dynamics of the manipulator. Analyzing the dynamic equation (11), the equation mainly includes three items, Coriolis/centripetal force, gravity force, and inertial force. e load will affect the value of the Coriolis force/centripetal force, but this item of parallel manipulator is generally small and generally ignored. e gravity term is only related to the mass of the system. For parallel manipulators, the mass of the load is usually a constant value and can be easily measured, so it is not necessary to identify the load mass. erefore, we mainly analyze the effect of load on the inertial force of the system. e inertial forces and torques of the system are as follows: where B I p,xx , B I p,yy , and B I p,zz are the moments of inertia of moving platform in the static coordinate system, and B I p,xx , B I p,xx , and B I p,xx are the products of inertia of moving platform in the static coordinate system; B I p,xx , B I p,yy , and B I p,zz are the moments of inertia of load in the static coordinate system, and B I l,xx , B I l,xx and B I l,xx are the products of inertia of load in the static coordinate system. Equation (16) is the inertial force equation of the PM with load. e equation shows that the load dynamic parameters, such as mass, inertia tensor, and position of load centroid have an effect on the inertial force or moment of the system. Analyzing equation (16), we can see that the inertial force or moment on a certain degree of freedom is not only determined by the acceleration on that degree of freedom, but also determined by the acceleration on the other degrees of freedom. So the PM with load has strong dynamic coupling characteristics among the six degrees of freedom. In addition, it is easy to know that the dynamic coupling is mainly caused by load.
Equation (16) can be rewritten as e mass matrix of equation (17) e off-diagonal elements in M(q) are the cause of the dynamic coupling among the six DOFs. Extracting these nondiagonal elements, the coupled mass matrix M c (q) can be obtained: where 0 3×3 is third-order zero matrix and L 3×3 , J 3×3 , and K 3×3 are the third-order submatrix of M C (q). First, analyze the submatrix L. e nonzero elements in the matrix L are the coefficients of angular acceleration in the inertial force Complexity 5 equation. us, the element value in this matrix can be used to measure the dynamic coupling influence of rotational motion on the translational motion. e larger the element value, the stronger the coupling effect. For example, if the matrix is a zero matrix, the coefficients representing the angular acceleration in the inertial force term are all zero, indicating that the rotational motion does not affect the translational motion. In the same way, the matrix L is composed of the coefficients of each translational acceleration in the moment of inertia, which represents the influence of translational motion on rotational motion. It can be seen from equation (19) that the elements in J and L are only related to the load, and the value of each nonzero element is determined by the load mass and the position of the centroid, so it is important to know their values a priori. Matrix K describes the dynamic coupling between the three rotational degrees of freedom and the nondiagonal elements in the matrix K are mainly determined by the load dynamic parameters.
In summary, the load has a greater impact on the dynamics of the system, especially the inertial force of the load. erefore, the dynamic identification of the load is very important.

Identification Process
In order to eliminate the influence of the load on the dynamic characteristics of the parallel manipulator and ensure the application of some advanced model-based control strategies, it is necessary to identify the parameters of the dynamic model of the load. Since the mass of the load is generally easier to determine, this paper only studies the parameter identification of the position of the centroid and the inertia tensor of the load. e inertia tensor and the position of the center of mass used in equation (11) are referred to as the dynamic parameters in the static coordinate system. e inertia tensor of load is a 3 * 3 matrix and can be expressed by where A I L is the inertial tensor of load in moving coordinate system. e position of the center of mass is a column vector and its equation is From equations (20) and (21), it is easy to know that the dynamic parameters are not constant and they change with the excitation trajectory of three rotational DOFs. However, the parameter identification of variables is very difficult and the corresponding real-time identification algorithm needs to be designed. To cope with this problem, the inertia tensor and the position of the center of mass relative to the moving coordinate system are selected as a base parameters set because these two parameters are constant and not affected by the trajectory of excitation, which can reduce the difficulty of identification.

Excitation Trajectory.
Although there is no variable in the basic parameters set, the equation of inertial force/ moment with regard to the identified parameters is not simple. In equation (16), the identification parameters are coupled with each other, which is impossible to separate. In order to reduce the coupling degree of the identified parameters in the dynamic equation, it is necessary to simplify the dynamic equation with regard to the identified parameters. When the given signal is a translational degree of freedom signal with a small amplitude, the coupling amplitude of each rotational degree of freedom is relatively small, and the system can be regarded as a small translational movement near the zero position, and the rotation matrix R can be treated as the identity matrix. In addition, when the given signal is a rotation signal with a small amplitude (less than 1 degree), the rotation matrix can also be regarded as the identity matrix. When the rotation matrix R is the identity matrix, formula (16) can be simplified to formula e identified parameters in equation (22) are A I l,xx , A I l,yy , A I l,zz , and Δz. In order to obtain sufficient information on these parameters, the exciting trajectory should be well designed.
Equation (3) requires that the amplitude of the given excitation signal on rotational direction should be small enough. Because when the rotation degree of freedom excitation signal is applied, the simplified model will 6 Complexity have a large modeling error relative to the original model. Due to the dynamic coupling characteristics between Dx and Ry, the motion on Dx and Ry both can excite dynamic parameters A I l,yy and Δz. For these two parameters, the dx direction is selected to apply an excitation signal because the translation signal satisfies equation (6) better than the rotation signal. Similarly, the dx direction is selected to apply an excitation signal to excite dynamic parameters A I l,xx and Δz. Since the last dynamic parameter A I l,zz only exits in the equation of M I,rz , the Rz direction is selected to apply an excitation signal to excite dynamic parameter A I l,zz . Note that the amplitude of the excitation signal on Rz direction should satisfy the condition of equation (22). e excitation trajectory of the system is as follows: It can be seen from the above equation that the excitation signals of the system are two translational signals, one rotational signal, and the rotational signal with a small amplitude, so it satisfies the system dynamics equation (22). Note that the excitation signal in the rotation direction cannot be too small; otherwise, the dynamic characteristics of the system cannot be well excited.

EKF Identification Process.
Select the state variable of the system: angular velocity in Rx, Ry, and Rz directions, centroid eccentricity z, and load moment of inertia A I l,xx , A I l,yy , A I l,zz . e state vector of the system is e a priori estimated value and the actual value of the system are not necessarily equal, so the a priori estimated value needs to be corrected by the measured value.
e main purpose of the measurement update equation is to find the optimal estimated value of the current iteration of the system. Update the Kalman gain, the optimal estimate of the system, and the posterior error covariance. e measurement update equation is In the equation, K is the Kalman gain and P is the posterior error covariance matrix. e system will produce measurement noise during measurement. is noise is white noise. R in the measurement update equation is the covariance matrix of the measurement noise. Now, we summarize the EKF identification algorithm for the parallel robot manipulator as shown in Algorithm 1.

Numerical Simulations
In this section, the simulation analysis is carried out in MATLAB/Simulink. For the simulation model, the mechanical model is built using Multibody and the sampling time is set to 1 ms. Moreover, the main parameters of the Stewart robot with uncertain load are shown in Table 1.
Corresponding simulation strategies are designed for load dynamics parameter identification, as shown in Figure 3.
First, a specific excitation signal is used to excite the dynamic characteristics of the load. en, the displacement, velocity, and output force of each leg were collected. According to the kinematic analysis of the parallel manipulator, the displacement, velocity, and force in the DOF space can be obtained by the kinematic positive solution or the velocity Jacobian matrix transformation of these measured values. Finally, the displacement, velocity, and force in the DOF space are input into the estimator and the estimator can calculate the dynamic parameters of the load based on the previously designed identification algorithm.
According to Figure 4, the corresponding simulation system was built in MATLAB/Simulink, as shown in Figure 4. e first is a signal generator block used to generate a target trajectory. e second is the kinematics block and its main function is to convert the task space signal to the joint space signal and calculate the speed Jacobian matrix. e third is the PID controller block. e fourth is the hydraulic system block.
e fifth is the mechanical model of the Stewart parallel mechanism. e sixth is the EKF estimator.
In the simulation, the initial estimated value of the state variables is all 0 and the initial error covariance matrix is diag 1 1 1 1 10 500 500 . Moreover, the value of the dynamic parameters to be identified is shown in Table 2.
e identification results of load dynamics parameters are shown in Figure 4.
It can be seen from Figure 5 that the proposed method in this paper can estimate the dynamic parameters of loads. It can be seen from Figure 5(a) that this identification method has a high estimation accuracy for centroid eccentricity along the z-axis direction, and the identification of Δz curve in the figure almost overlaps with the real value. Figures 5(b)-5(d) are identification curves of load dynamics parameters A I l,xx , A I l,yy , and A I l,zz , respectively. Although there is a certain error between the values of the three inertial parameters and the real values, the error is not large. Table 3 shows the true value X r , identification value X id , and relative error e r of the load dynamic parameters obtained by the proposed method. e relative error e r is defined as Algorithm extended Kalman filter algorithm (1) Initialize the estimated value of the state variable and the error covariance matrix (2) Repeat (3) Calculate prior state estimates, x − k � f(x k−1 ) (4) Calculate the Jacobian matrix, J k−1 (5) Calculate the prior error covariance,

Conclusion
e proposed method has been successfully applied to a sixdegree-of-freedom Gough-Stewart parallel manipulator for load dynamic parameters estimation. e identification algorithm is simple and easy to implement. Compared with the traditional least square method, the proposed identification approach does not require linearization of the dynamic model and optimization of the excitation trajectory. Moreover, this method is not sensitive to measurement noise and without acceleration measurements in the process of identification.
Note that the method proposed in this paper is an offline identification method. How to study a load real-time identification algorithm not limited to trajectory needs further study.

Data Availability
All data included in this study are available upon request to the corresponding author.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.