Abstract

Load is the main external disturbance of a parallel robot manipulator. This disturbance will cause dynamic coupling among different degrees of freedom and make heaps of model-based control methods difficult to apply. In order to compensate this disturbance, it is crucial to obtain an accurate dynamic model of load. However, in practice, the load is always uncertain and its dynamic parameters are arduous to know a priori. To cope with this problem, this paper proposes a novel and simple approach to identify the dynamic parameters of load. Firstly, the dynamic model of the parallel robot manipulator with uncertain load is established and the dynamic coupling caused by load is also analyzed. Then, according to the dynamic model, the excitation signal is designed and a weak nonlinear dynamic model is derived. Furthermore, the identification model is presented and the identification algorithm based on the extended Kalman filter is designed. Lastly, numerical simulation results, obtained using a six-degree-of-freedom Gough–Stewart parallel manipulator, demonstrate the good estimation performance of the proposed method.

1. 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. The 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]. Thus, 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. The first method is using computer-aided design (CAD) techniques to obtain the dynamic parameters of robots. The 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. This 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. The 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–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]. Thanh 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. The main reasons are as follows:(1)The 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)The 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.

This 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.

The 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.

2. Mathematical Modeling

The parallel manipulator studied in this paper is a six-degree-of-freedom Gough–Stewart parallel manipulator, as shown in Figure 1(a). This manipulator is mainly composed of two platforms and six driving legs. The lower platform fixed on the ground is named static platform and the upper platform used to carry loads is named moving platform. The driving leg is the actuator to drive the moving platform to realize translation and rotation. In order to facilitate the analysis of the parallel manipulator, the coordinate system is established and shown in Figure 1(b). The static coordinate system is fixed on the static platform and the origin O_B is in the center of the lower platform. The moving coordinate is fixed on the upper platform and the origin O is in the center of the upper platform. Not only is O_A the centroid point of the upper platform, but also it is the control point C of the system. G_P is the centroid point of load. a_i and b_i are the upper and lower hinge points, respectively.

(a)

(b)

(a)
(b)

Figure 1 

Schematic of the Stewart robot. (a) Simplified structure. (b) Mathematical representation.

2.1. Kinematics Modeling

According to the space geometry theory, the length vector of the six legs can be expressed bywhere is the coordinate matrix of six lower hinges in the static coordinate. is the coordinate matrix of six upper hinges in the moving coordinate. denotes the transformation matrix from the static coordinate to the moving coordinate. is a translation displacement vector. , and are the displacement of platform mass center along the x-axis, y-axis, and z-axis, respectively. is the initial height matrix composed of six initial height vectors, and each vector is . is the initial height when the upper platform motion table in the middle position.

By differentiating equation (1), one obtains

The matrix form of equation (2) iswhere is velocity Jacobian matrix between the velocity of the upper platform and the leg velocity. is the unit direction vector matrix of the six legs’ direction and . is the upper platform velocity, in which is the angular velocity of the moving platform and is the leg velocity. Importantly, matrix 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 , , and . The load eccentric position along three axes of the static coordinate can be expressed aswhere c () and s () are abbreviations for the trigonometric functions cosine and sine; q₄, q₅, and q₆ are the three Euler angles of the system; and , , and are the load eccentric position along three main axes. From equation (4), it is easy to know that the , , and are related to the Euler anglers of the upper platform.

Let the load centroid acceleration in the static coordinate be . According to the rigid body dynamics theory, the relation between the load centroid acceleration and the upper platform centroid acceleration can be expressed by

The relationship between and is shown in Figure 2. Note that the coordinate system in Figure 1 is not the same as the coordinate system in Figure 2. The origin of the 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. coordinate system is the coordinate system used to establish the dynamic equation. In addition, is the position of the load center of mass.

Figure 2 

Kinematic relationship.

2.2. 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 is the mass of load; is the inertia tensor of load in the static coordinate system; and 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 is the force of hydraulic cylinder of ith leg; is the unit direction vector of the ith legs; is inertia force of load; is inertia moment of load; is the mass of moving platform; and 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. The 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 thatwhere is the driving force vector in the task space; is the centrifugal/Coriolis term; is the gravity term; and is the mass of the system.

The 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 aswhere is the system mass matrix in joint space, is the centrifugal/Coriolis term in joint space, is the gravity term in joint space, and

The impact of the centrifugal/Coriolis term is small, which can be ignored, and equation (13) can be rewritten as

2.3. 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. The load will affect the value of the Coriolis force/centripetal force, but this item of parallel manipulator is generally small and generally ignored. The 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. Therefore, we mainly analyze the effect of load on the inertial force of the system.

The inertial forces and torques of the system are as follows:where , , and are the moments of inertia of moving platform in the static coordinate system, and , , and are the products of inertia of moving platform in the static coordinate system; , , and are the moments of inertia of load in the static coordinate system, and , and are the products of inertia of load in the static coordinate system.

Equation (16) is the inertial force equation of the PM with load. The 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

The mass matrix of equation (17) can be expressed by

The off-diagonal elements in are the cause of the dynamic coupling among the six DOFs. Extracting these nondiagonal elements, the coupled mass matrix can be obtained:where is third-order zero matrix and , , and are the third-order submatrix of . First, analyze the submatrix . The nonzero elements in the matrix are the coefficients of angular acceleration in the inertial force equation. Thus, the element value in this matrix can be used to measure the dynamic coupling influence of rotational motion on the translational motion. The 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 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 and 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 describes the dynamic coupling between the three rotational degrees of freedom and the nondiagonal elements in the matrix 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. Therefore, the dynamic identification of the load is very important.

3. 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. The 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.

The inertia tensor of load is a matrix and can be expressed bywhere is the inertial tensor of load in moving coordinate system.

The position of the center of mass is a column vector and its equation iswhere 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.

3.1. 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

The identified parameters in equation (22) are , , , and . 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 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 and . 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 and . Since the last dynamic parameter only exits in the equation of , the Rz direction is selected to apply an excitation signal to excite dynamic parameter . Note that the amplitude of the excitation signal on Rz direction should satisfy the condition of equation (22).

The 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.

3.2. 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 , , .

The state vector of the system is

Establish the system differential equation as follows:

The discrete form of formula (25) is as follows:

The derivative of the state vector is as follows:

From formulas (26) and (27), the state equation of the system is

The observed values of the system are the angular velocities in the directions of the three degrees of freedom of Rx, Ry, and Rz, and the observation equation is as follows:

The observation vector and the observation noise vector are as follows:

Equations (28) and (29) are the state space equation and observation equation of the system, respectively, which are written aswhere is the observation matrix of the system.

Equation (31) is the system model that needs to be identified. The extended Kalman filter algorithm is used to estimate the parameters of the model. The extended Kalman filter algorithm is mainly divided into two parts, the time update part and the measurement update part. The two parts are introduced separately below.

3.2.1. Time Update

The main purpose of the time update equation is to calculate the next prior estimate vector on the basis of the last optimal estimate vector and also to update the prior error covariance matrix. The time update (prediction) equation is

When the system performs state transition, there will be state transition noise, which is Gaussian white noise, and its covariance matrix is in the time update equation.

is the Jacobian matrix of the system, which is obtained by the partial derivative of the state function f on the state vector x. The expression of the Jacobian matrix of the identification model in this paper is as follows:

is the submatrix of , and its expression is as follows:where