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

Authors: Shijie Song; Xiaolin Dai; Zhangchao Huang; Dawei Gong
Journal: Complexity
 (2020)

---

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


---

## Body

## 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 Section2, 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 Figure1(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 OB−XBYBZB is fixed on the static platform and the origin OB is in the center of the lower platform. The moving coordinate OA−XAYAZA is fixed on the upper platform and the origin O is in the center of the upper platform. Not only is OA the centroid point of the upper platform, but also it is the control point C of the system. GP is the centroid point of load. ai and bi are the upper and lower hinge points, respectively.Figure 1
Schematic of the Stewart robot. (a) Simplified structure. (b) Mathematical representation.
(a)(b)


### 2.1. Kinematics Modeling

According to the space geometry theory, the length vector of the six legs can be expressed by(1)L=L1,…,L6=R⋅A+p+p0−B,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=xyzT 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. p0 is the initial height matrix composed of six initial height vectors, and each vector is 00h0T. h0 is the initial height when the upper platform motion table in the middle position.By differentiating equation (1), one obtains(2)L˙=dLdt=p˙+R˙⋅RT⋅R⋅A=p˙+ω×RAA.The matrix form of equation (2) is(3)L˙=LnTR⋅A×LnT⋅q˙=Jvq˙,where Jv is velocity Jacobian matrix between the velocity of the upper platform and the leg velocity. Ln is the unit direction vector matrix of the six legs’ direction and Ln=L1/L1,…,L6/L6. q˙=q˙=p˙,ωT is the upper platform velocity, in which ω=ωxωyωzT 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. The load eccentric position along three axes of the static coordinate can be expressed as(4)lx=cq5cq6Δx+sq4sq5cq6cq4sq6Δy+cq4sq5cq6+sq4sq6Δz,ly=cq5sq6Δx+sq4sq5sq6+cq4cq6Δy+cq4sq5sq6−sq4cq6Δz,lz=−sq5Δx+sq4cq5Δy+cq4cq5Δz,where c () and s () are abbreviations for the trigonometric functions cosine and sine; q4, q5, and q6 are the three Euler angles of the system; and lx, ly, and lz are the load eccentric position along three main axes. From equation (4), it is easy to know that the lx, ly, and lz are related to the Euler anglers of the upper platform.Let the load centroid acceleration in the static coordinate bep¨l=x¨ay¨az¨aT. According to the rigid body dynamics theory, the relation between the load centroid acceleration and the upper platform centroid acceleration can be expressed by(5)x¨a=x¨+lzω˙y−lyω˙z,y¨a=y¨+lzω˙x−lxω˙z,z¨a=z¨+lyω˙x−lxω˙y.The relationship betweenp¨l and p¨ is shown in Figure 2. Note that the OC−XCYCZC coordinate system in Figure 1 is not the same as the OB−XBYBZB coordinate system in Figure 2. The origin of the OC−XCYCZC 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. OC−XCYCZC coordinate system is the coordinate system used to establish the dynamic equation. In addition, Gp 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:(6)fL=mlp¨l=mlp¨−mlQω˙,(7)TL=ILBω˙−mlQp¨,where ml is the mass of load; ILB 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:(8)Q=0−lzlylz0−lx−lylx0.Based on the Newton–Euler method, the dynamic model in task space can be obtained:(9)∑i=16ei⋅fi=fL+mpp¨+mp+mlg,(10)∑i=16R⋅ai×ei⋅fi=TL+IPBω˙+ω×ILB+IPBω,where fi is the force of hydraulic cylinder of ith leg; ei is the unit direction vector of the ith legs; fL is inertia force of load; TL is inertia moment of load; mp is the mass of moving platform; and IPB 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 that(11)τ=Mq⋅q¨+Hq˙,q⋅q˙+G,where τ is the driving force vector in the task space; Hq˙,q is the centrifugal/Coriolis term; G is the gravity term; and Mq 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(12)τ=JTf.From equations (11) and (12), the dynamic model of the robot in joint space can be rewritten as(13)f=MLl⋅l¨+HLl,l˙⋅l˙+GL,where MLl is the system mass matrix in joint space, HLl,l˙ is the centrifugal/Coriolis term in joint space, GL is the gravity term in joint space, and(14)MLl=J−1TMqJ−1,HLl,l˙=J−1TMqJ˙−1+Hq˙,qJ−1,GL=J−1TGq.The impact of the centrifugal/Coriolis term is small, which can be ignored, and equation (13) can be rewritten as(15)f=MLl⋅l¨+GL.


### 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:(16)fI,x=max¨+mpx¨+mplzω˙y−mplyω˙z,fI,y=may¨+mpy¨−mplzω˙x+mplxω˙z,fI,z=maz¨+mpz¨+mplyω˙x+mplxω˙y,MI,rx=Ip,xxB+Il,xxBω˙x−mplzy¨+mplyz¨−mplxlyω˙y−mplxlzω˙z−Il,xyBω˙y−Il,xzBω˙z−Ip,xyBω˙y−Ip,xzBω˙z,MI,ry=Ip,yyB+Il,yyBω˙y+mplzx¨−mplxz¨−mplylzω˙z−mplxlyω˙x−Il,xyBω˙x−Il,yzBω˙z−Ip,xyBω˙x−Ip,yzBω˙z,MI,rz=Ip,zzB+Il,zzBω˙z−mplyx¨+mplxy¨−mplxlzω˙x−mplylzω˙y−Il,xzBω˙x−Il,yzBω˙y−Ip,xzBω˙x−Ip,yzBω˙y,where Ip,xxB, Ip,yyB, and Ip,zzB are the moments of inertia of moving platform in the static coordinate system, and Ip,xxB, Ip,xxB, and Ip,xxB are the products of inertia of moving platform in the static coordinate system; Ip,xxB, Ip,yyB, and Ip,zzB are the moments of inertia of load in the static coordinate system, and Il,xxB, Il,xxB and Il,xxB 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(17)FI=Mqq¨.The mass matrix of equation (17) can be expressed by(18)Mq=m000mplz−mply0m0−mplz0mplx00mmply−mplx00−mplzmplyIl,xxB+Ip,xxB−Ip,xyB−Il,xyB−mplxly−Ip,xzB−Il,xzB−mplxlzmplz0−mplx−Ip,xyB−Il,xyB−mplxlyIl,yyB+Ip,yyB−Ip,yzB−Il,yzB−mplylz−mplymplx0−Ip,xzB−Il,xzB−mplxlz−Ip,yzB−Il,yzB−mplylzIl,zzB+Ip,zzB.The off-diagonal elements inMq are the cause of the dynamic coupling among the six DOFs. Extracting these nondiagonal elements, the coupled mass matrix Mcq can be obtained:(19)MCq=03×3L3×3J3×3K3×3=0000mplz−mply000−mplz0mplx000mply−mplx00−mplzmply0−Ip,xyB−Il,xyB−mplxly−Ip,xzB−Il,xzB−mplxlzmplz0−mplx−Ip,xyB−Il,xyB−mplxly0−Ip,yzB−Il,yzB−mplylz−mplymplx0−Ip,xzB−Il,xzB−mplxlz−Ip,yzB−Il,yzB−mplylz0,where 03×3 is third-order zero matrix and L3×3, J3×3, and K3×3 are the third-order submatrix of MCq. First, analyze the submatrix L. The nonzero elements in the matrix L 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 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. Therefore, the dynamic identification of the load is very important.


## 2.1. Kinematics Modeling

According to the space geometry theory, the length vector of the six legs can be expressed by(1)L=L1,…,L6=R⋅A+p+p0−B,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=xyzT 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. p0 is the initial height matrix composed of six initial height vectors, and each vector is 00h0T. h0 is the initial height when the upper platform motion table in the middle position.By differentiating equation (1), one obtains(2)L˙=dLdt=p˙+R˙⋅RT⋅R⋅A=p˙+ω×RAA.The matrix form of equation (2) is(3)L˙=LnTR⋅A×LnT⋅q˙=Jvq˙,where Jv is velocity Jacobian matrix between the velocity of the upper platform and the leg velocity. Ln is the unit direction vector matrix of the six legs’ direction and Ln=L1/L1,…,L6/L6. q˙=q˙=p˙,ωT is the upper platform velocity, in which ω=ωxωyωzT 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. The load eccentric position along three axes of the static coordinate can be expressed as(4)lx=cq5cq6Δx+sq4sq5cq6cq4sq6Δy+cq4sq5cq6+sq4sq6Δz,ly=cq5sq6Δx+sq4sq5sq6+cq4cq6Δy+cq4sq5sq6−sq4cq6Δz,lz=−sq5Δx+sq4cq5Δy+cq4cq5Δz,where c () and s () are abbreviations for the trigonometric functions cosine and sine; q4, q5, and q6 are the three Euler angles of the system; and lx, ly, and lz are the load eccentric position along three main axes. From equation (4), it is easy to know that the lx, ly, and lz are related to the Euler anglers of the upper platform.Let the load centroid acceleration in the static coordinate bep¨l=x¨ay¨az¨aT. According to the rigid body dynamics theory, the relation between the load centroid acceleration and the upper platform centroid acceleration can be expressed by(5)x¨a=x¨+lzω˙y−lyω˙z,y¨a=y¨+lzω˙x−lxω˙z,z¨a=z¨+lyω˙x−lxω˙y.The relationship betweenp¨l and p¨ is shown in Figure 2. Note that the OC−XCYCZC coordinate system in Figure 1 is not the same as the OB−XBYBZB coordinate system in Figure 2. The origin of the OC−XCYCZC 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. OC−XCYCZC coordinate system is the coordinate system used to establish the dynamic equation. In addition, Gp 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:(6)fL=mlp¨l=mlp¨−mlQω˙,(7)TL=ILBω˙−mlQp¨,where ml is the mass of load; ILB 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:(8)Q=0−lzlylz0−lx−lylx0.Based on the Newton–Euler method, the dynamic model in task space can be obtained:(9)∑i=16ei⋅fi=fL+mpp¨+mp+mlg,(10)∑i=16R⋅ai×ei⋅fi=TL+IPBω˙+ω×ILB+IPBω,where fi is the force of hydraulic cylinder of ith leg; ei is the unit direction vector of the ith legs; fL is inertia force of load; TL is inertia moment of load; mp is the mass of moving platform; and IPB 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 that(11)τ=Mq⋅q¨+Hq˙,q⋅q˙+G,where τ is the driving force vector in the task space; Hq˙,q is the centrifugal/Coriolis term; G is the gravity term; and Mq 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(12)τ=JTf.From equations (11) and (12), the dynamic model of the robot in joint space can be rewritten as(13)f=MLl⋅l¨+HLl,l˙⋅l˙+GL,where MLl is the system mass matrix in joint space, HLl,l˙ is the centrifugal/Coriolis term in joint space, GL is the gravity term in joint space, and(14)MLl=J−1TMqJ−1,HLl,l˙=J−1TMqJ˙−1+Hq˙,qJ−1,GL=J−1TGq.The impact of the centrifugal/Coriolis term is small, which can be ignored, and equation (13) can be rewritten as(15)f=MLl⋅l¨+GL.


## 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:(16)fI,x=max¨+mpx¨+mplzω˙y−mplyω˙z,fI,y=may¨+mpy¨−mplzω˙x+mplxω˙z,fI,z=maz¨+mpz¨+mplyω˙x+mplxω˙y,MI,rx=Ip,xxB+Il,xxBω˙x−mplzy¨+mplyz¨−mplxlyω˙y−mplxlzω˙z−Il,xyBω˙y−Il,xzBω˙z−Ip,xyBω˙y−Ip,xzBω˙z,MI,ry=Ip,yyB+Il,yyBω˙y+mplzx¨−mplxz¨−mplylzω˙z−mplxlyω˙x−Il,xyBω˙x−Il,yzBω˙z−Ip,xyBω˙x−Ip,yzBω˙z,MI,rz=Ip,zzB+Il,zzBω˙z−mplyx¨+mplxy¨−mplxlzω˙x−mplylzω˙y−Il,xzBω˙x−Il,yzBω˙y−Ip,xzBω˙x−Ip,yzBω˙y,where Ip,xxB, Ip,yyB, and Ip,zzB are the moments of inertia of moving platform in the static coordinate system, and Ip,xxB, Ip,xxB, and Ip,xxB are the products of inertia of moving platform in the static coordinate system; Ip,xxB, Ip,yyB, and Ip,zzB are the moments of inertia of load in the static coordinate system, and Il,xxB, Il,xxB and Il,xxB 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(17)FI=Mqq¨.The mass matrix of equation (17) can be expressed by(18)Mq=m000mplz−mply0m0−mplz0mplx00mmply−mplx00−mplzmplyIl,xxB+Ip,xxB−Ip,xyB−Il,xyB−mplxly−Ip,xzB−Il,xzB−mplxlzmplz0−mplx−Ip,xyB−Il,xyB−mplxlyIl,yyB+Ip,yyB−Ip,yzB−Il,yzB−mplylz−mplymplx0−Ip,xzB−Il,xzB−mplxlz−Ip,yzB−Il,yzB−mplylzIl,zzB+Ip,zzB.The off-diagonal elements inMq are the cause of the dynamic coupling among the six DOFs. Extracting these nondiagonal elements, the coupled mass matrix Mcq can be obtained:(19)MCq=03×3L3×3J3×3K3×3=0000mplz−mply000−mplz0mplx000mply−mplx00−mplzmply0−Ip,xyB−Il,xyB−mplxly−Ip,xzB−Il,xzB−mplxlzmplz0−mplx−Ip,xyB−Il,xyB−mplxly0−Ip,yzB−Il,yzB−mplylz−mplymplx0−Ip,xzB−Il,xzB−mplxlz−Ip,yzB−Il,yzB−mplylz0,where 03×3 is third-order zero matrix and L3×3, J3×3, and K3×3 are the third-order submatrix of MCq. First, analyze the submatrix L. The nonzero elements in the matrix L 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 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. 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 a3∗3 matrix and can be expressed by(20)ILB=R⋅ILA⋅RT,where ILA is the inertial tensor of load in moving coordinate system.The position of the center of mass is a column vector and its equation is(21)bB=R⋅bA,where bA=ΔxΔyΔzT is bB=lxlylzT.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(22)fI,xt=max¨t+mpx¨t+mpΔzω˙yt,fI,yt=may¨t+mpy¨t−mpΔztω˙x,fI,zt=maz¨t+mpz¨t,MI,rxt=Ip,xxA+Il,xxAt+mpΔzt2ω˙xt−mpΔztyt,MI,ryt=Ip,yyA+Il,yyAt+mpΔzt2ω˙yt+mpΔztx¨t,MI,rzt=Ip,zzA+Il,zzAtω˙zt.The identified parameters in equation (22) are Il,xxA, Il,yyA, Il,zzA, 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 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 Il,yyA 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 Il,xxA and Δz. Since the last dynamic parameter Il,zzA only exits in the equation of MI,rz, the Rz direction is selected to apply an excitation signal to excite dynamic parameter Il,zzA. 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:(23)dx=0.2⋅sinπt,dy=0.2⋅sinπt,rz=0.2⋅sinπt.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 inRx, Ry, and Rz directions, centroid eccentricity z, and load moment of inertia Il,xxA, Il,yyA, Il,zzA.The state vector of the system is(24)x=vrxtvrytvrztztIl,xxAtIl,yyAtIl,zzAtT.Establish the system differential equation as follows:(25)v˙rxt=mpztfI,yt+mMI,rxtmIp,xxA+Il,xxAt+mpzt2−mp2zt2,v˙ryt=mpztfI,xt−mMI,rytmp2zt2−mIp,yyA+Il,yyAt+mpzt2,v˙rzt=MI,rztIp,zzA+Il,zzAt,z˙t=0,I˙l,xxAt=0I˙l,yyAt=0,I˙l,zzAk=0.The discrete form of formula (25) is as follows:(26)v˙rxk=mpzk−1fI,yk−1+mMI,rxk−1mIp,yyA+Il,yyAtk−1+mpzk−12−mp2zk−12,v˙ryk=mpzk−1fI,xk−1−mMI,ryk−1mp2zk−12−mIp,xxA+Il,xxAtk−1+mpzk−12,v˙rzk=MI,rzk−1Ip,zzA+Il,zzAtk−1,z˙k=0,I˙pc,xxk=0,Ipc,yyk=0,I˙pc,zzk=0.The derivative of the state vector is as follows:(27)x˙k=v˙rxkv˙rykv˙rzkz˙kI˙l,xxAkI˙l,yyAkI˙l,zzAkT.From formulas (26) and (27), the state equation of the system is(28)xk=x˙k−1Δt+xk−1+Wk−1.The observed values of the system are the angular velocities in the directions of the three degrees of freedom ofRx, Ry, and Rz, and the observation equation is as follows:(29)zrxk=vrxk+Vrxk,zryk=vryk+Vryk,zrzk=vrzk+Vrzk.The observation vector and the observation noise vector are as follows:(30)zk=zrxkzrxkzrxkT,Vk=VrxkVrxkVrxkT.Equations (28) and (29) are the state space equation and observation equation of the system, respectively, which are written as(31)xk=fxk−1+Wk−1,zk=Hxk+Vk−1,where H=E3×304×3 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(32)x^k−=fx^k−1,Pk−=Jk−1Pk−1Jk−1T+Qk−1.When the system performs state transition, there will be state transition noise, which is Gaussian white noise, and its covariance matrix isQ in the time update equation.J 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:(33)J=∂f1∂x1∂f1∂x2…∂f1∂x7∂f2∂x1∂f2∂x2…∂f2∂x7⋮⋮⋱⋮∂f7∂x1∂f7∂x2…∂f7∂x7=E3×3N04×3E4×4.N is the submatrix of J, and its expression is as follows:(34)N=j11j1200j210j220000j34,where(35)j11=ΔtfI,xk−1mpmp2zk−12−mIp,xxA+Il,xxAtk−1+mp2z2k−1,−Δt2mp2zk−1−2mmpzk−1mpzk−1fI,xk−1−mMI,ryk−1mp2z2k−1−mIp,xxA+Il,xxAtk−1+mpz2k−12,j12=mmpzk−1fI,xk−1−mMI,ryk−1mp2z2k−1−mIp,xxA+Il,xxAtk−1+mpz2k−12Δt,j21=fI,yk−1mpΔtmIp,yyA+Il,yyAt+mpzk−12−mp2zk−12−2mmpzk−1−2mp2zk−1mpzk−1fI,yk−1+mMI,rxk−1ΔtmIp,yyA+Il,yyAtk−1+mpzk−12−mp2zk−122,j22=mmpzk−1fI,yk−1+mMI,rxk−1ΔtmIp,yyA+Il,yyAk−1+mpzk−12−mp2zk−122,j34=MI,rzk−1Ip,zzA+Il,zzAk−12Δt.


#### 3.2.2. Measurement Update

The 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. The 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. The measurement update equation is(36)Kk=Pk−HTHPk−HT+Rk−1,x^k=x^k−+Kkzk−Hx^k−,Pk=I−KkHPk−.In the equation,K is the Kalman gain and P is the posterior error covariance matrix. The system will produce measurement noise during measurement. This 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.Algorithm 1: Identification algorithm.
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−=fx^k−1(4)
Calculate the Jacobian matrix,Jk−1(5)
Calculate the prior error covariance,Pk−=Jk−1Pk−1Jk−1T+Qk−1(6)
Calculate Kalman gain,Kk=PK−HTHPK−HT+Rk−1(7)
Update the posterior state estimate,x^k=x^K−+Kkzk−Hx^K−(8)
Update posterior error covariance,Pk=I−KkHPK−(9)
Output the best estimate of this iteration(10)
Until Simulation stopped(11)
End


## 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(22)fI,xt=max¨t+mpx¨t+mpΔzω˙yt,fI,yt=may¨t+mpy¨t−mpΔztω˙x,fI,zt=maz¨t+mpz¨t,MI,rxt=Ip,xxA+Il,xxAt+mpΔzt2ω˙xt−mpΔztyt,MI,ryt=Ip,yyA+Il,yyAt+mpΔzt2ω˙yt+mpΔztx¨t,MI,rzt=Ip,zzA+Il,zzAtω˙zt.The identified parameters in equation (22) are Il,xxA, Il,yyA, Il,zzA, 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 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 Il,yyA 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 Il,xxA and Δz. Since the last dynamic parameter Il,zzA only exits in the equation of MI,rz, the Rz direction is selected to apply an excitation signal to excite dynamic parameter Il,zzA. 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:(23)dx=0.2⋅sinπt,dy=0.2⋅sinπt,rz=0.2⋅sinπt.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 inRx, Ry, and Rz directions, centroid eccentricity z, and load moment of inertia Il,xxA, Il,yyA, Il,zzA.The state vector of the system is(24)x=vrxtvrytvrztztIl,xxAtIl,yyAtIl,zzAtT.Establish the system differential equation as follows:(25)v˙rxt=mpztfI,yt+mMI,rxtmIp,xxA+Il,xxAt+mpzt2−mp2zt2,v˙ryt=mpztfI,xt−mMI,rytmp2zt2−mIp,yyA+Il,yyAt+mpzt2,v˙rzt=MI,rztIp,zzA+Il,zzAt,z˙t=0,I˙l,xxAt=0I˙l,yyAt=0,I˙l,zzAk=0.The discrete form of formula (25) is as follows:(26)v˙rxk=mpzk−1fI,yk−1+mMI,rxk−1mIp,yyA+Il,yyAtk−1+mpzk−12−mp2zk−12,v˙ryk=mpzk−1fI,xk−1−mMI,ryk−1mp2zk−12−mIp,xxA+Il,xxAtk−1+mpzk−12,v˙rzk=MI,rzk−1Ip,zzA+Il,zzAtk−1,z˙k=0,I˙pc,xxk=0,Ipc,yyk=0,I˙pc,zzk=0.The derivative of the state vector is as follows:(27)x˙k=v˙rxkv˙rykv˙rzkz˙kI˙l,xxAkI˙l,yyAkI˙l,zzAkT.From formulas (26) and (27), the state equation of the system is(28)xk=x˙k−1Δt+xk−1+Wk−1.The observed values of the system are the angular velocities in the directions of the three degrees of freedom ofRx, Ry, and Rz, and the observation equation is as follows:(29)zrxk=vrxk+Vrxk,zryk=vryk+Vryk,zrzk=vrzk+Vrzk.The observation vector and the observation noise vector are as follows:(30)zk=zrxkzrxkzrxkT,Vk=VrxkVrxkVrxkT.Equations (28) and (29) are the state space equation and observation equation of the system, respectively, which are written as(31)xk=fxk−1+Wk−1,zk=Hxk+Vk−1,where H=E3×304×3 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(32)x^k−=fx^k−1,Pk−=Jk−1Pk−1Jk−1T+Qk−1.When the system performs state transition, there will be state transition noise, which is Gaussian white noise, and its covariance matrix isQ in the time update equation.J 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:(33)J=∂f1∂x1∂f1∂x2…∂f1∂x7∂f2∂x1∂f2∂x2…∂f2∂x7⋮⋮⋱⋮∂f7∂x1∂f7∂x2…∂f7∂x7=E3×3N04×3E4×4.N is the submatrix of J, and its expression is as follows:(34)N=j11j1200j210j220000j34,where(35)j11=ΔtfI,xk−1mpmp2zk−12−mIp,xxA+Il,xxAtk−1+mp2z2k−1,−Δt2mp2zk−1−2mmpzk−1mpzk−1fI,xk−1−mMI,ryk−1mp2z2k−1−mIp,xxA+Il,xxAtk−1+mpz2k−12,j12=mmpzk−1fI,xk−1−mMI,ryk−1mp2z2k−1−mIp,xxA+Il,xxAtk−1+mpz2k−12Δt,j21=fI,yk−1mpΔtmIp,yyA+Il,yyAt+mpzk−12−mp2zk−12−2mmpzk−1−2mp2zk−1mpzk−1fI,yk−1+mMI,rxk−1ΔtmIp,yyA+Il,yyAtk−1+mpzk−12−mp2zk−122,j22=mmpzk−1fI,yk−1+mMI,rxk−1ΔtmIp,yyA+Il,yyAk−1+mpzk−12−mp2zk−122,j34=MI,rzk−1Ip,zzA+Il,zzAk−12Δt.


### 3.2.2. Measurement Update

The 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. The 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. The measurement update equation is(36)Kk=Pk−HTHPk−HT+Rk−1,x^k=x^k−+Kkzk−Hx^k−,Pk=I−KkHPk−.In the equation,K is the Kalman gain and P is the posterior error covariance matrix. The system will produce measurement noise during measurement. This 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.Algorithm 1: Identification algorithm.
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−=fx^k−1(4)
Calculate the Jacobian matrix,Jk−1(5)
Calculate the prior error covariance,Pk−=Jk−1Pk−1Jk−1T+Qk−1(6)
Calculate Kalman gain,Kk=PK−HTHPK−HT+Rk−1(7)
Update the posterior state estimate,x^k=x^K−+Kkzk−Hx^K−(8)
Update posterior error covariance,Pk=I−KkHPK−(9)
Output the best estimate of this iteration(10)
Until Simulation stopped(11)
End


## 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(32)x^k−=fx^k−1,Pk−=Jk−1Pk−1Jk−1T+Qk−1.When the system performs state transition, there will be state transition noise, which is Gaussian white noise, and its covariance matrix isQ in the time update equation.J 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:(33)J=∂f1∂x1∂f1∂x2…∂f1∂x7∂f2∂x1∂f2∂x2…∂f2∂x7⋮⋮⋱⋮∂f7∂x1∂f7∂x2…∂f7∂x7=E3×3N04×3E4×4.N is the submatrix of J, and its expression is as follows:(34)N=j11j1200j210j220000j34,where(35)j11=ΔtfI,xk−1mpmp2zk−12−mIp,xxA+Il,xxAtk−1+mp2z2k−1,−Δt2mp2zk−1−2mmpzk−1mpzk−1fI,xk−1−mMI,ryk−1mp2z2k−1−mIp,xxA+Il,xxAtk−1+mpz2k−12,j12=mmpzk−1fI,xk−1−mMI,ryk−1mp2z2k−1−mIp,xxA+Il,xxAtk−1+mpz2k−12Δt,j21=fI,yk−1mpΔtmIp,yyA+Il,yyAt+mpzk−12−mp2zk−12−2mmpzk−1−2mp2zk−1mpzk−1fI,yk−1+mMI,rxk−1ΔtmIp,yyA+Il,yyAtk−1+mpzk−12−mp2zk−122,j22=mmpzk−1fI,yk−1+mMI,rxk−1ΔtmIp,yyA+Il,yyAk−1+mpzk−12−mp2zk−122,j34=MI,rzk−1Ip,zzA+Il,zzAk−12Δt.


## 3.2.2. Measurement Update

The 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. The 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. The measurement update equation is(36)Kk=Pk−HTHPk−HT+Rk−1,x^k=x^k−+Kkzk−Hx^k−,Pk=I−KkHPk−.In the equation,K is the Kalman gain and P is the posterior error covariance matrix. The system will produce measurement noise during measurement. This 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.Algorithm 1: Identification algorithm.
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−=fx^k−1(4)
Calculate the Jacobian matrix,Jk−1(5)
Calculate the prior error covariance,Pk−=Jk−1Pk−1Jk−1T+Qk−1(6)
Calculate Kalman gain,Kk=PK−HTHPK−HT+Rk−1(7)
Update the posterior state estimate,x^k=x^K−+Kkzk−Hx^K−(8)
Update posterior error covariance,Pk=I−KkHPK−(9)
Output the best estimate of this iteration(10)
Until Simulation stopped(11)
End


## 4. 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 Table1.Table 1
Main parameters of the Stewart robot with load.
ParametersValueRadius of upper/lower platform (m)0.4/0.6Initial length of the linear hydraulic cylinder (m)0.8741Min/max stroke of the hydraulic cylinder (m)−0.3/0.3Mass of load (kg)500Mass of upper platform (kg)100Moment of inertia of the upper platform (kg·m2)diag (25, 25, 44)Corresponding simulation strategies are designed for load dynamics parameter identification, as shown in Figure3.Figure 3
Load dynamic parameter identification system.First, a specific excitation signal is used to excite the dynamic characteristics of the load. Then, 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 Figure4, the corresponding simulation system was built in MATLAB/Simulink, as shown in Figure 4. The first is a signal generator block used to generate a target trajectory. The 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. The third is the PID controller block. The fourth is the hydraulic system block. The fifth is the mechanical model of the Stewart parallel mechanism. The sixth is the EKF estimator.Figure 4
Simulink simulation model of load dynamics parameter identification system.In the simulation, the initial estimated value of the state variables is all 0 and the initial error covariance matrix isdiag111110500500. Moreover, the value of the dynamic parameters to be identified is shown in Table 2.Table 2
Load dynamics parameters to be identified.
SymbolParametersValueUnitΔzCentroid position0.3mdiagIl,xxAIl,yyAIl,zzAMoment of inertia of load (kg·m2)diag (25, 25, 25)kg·m2The identification results of load dynamics parameters are shown in Figure4.It can be seen from Figure5 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 Il,xxA, Il,yyA, and Il,zzA, respectively. Although there is a certain error between the values of the three inertial parameters and the real values, the error is not large.Figure 5
Identification results of load dynamics parameters. (a) The identification value ofΔz. (b) The identification value of Il,xxA. (c) The identification value of Il,yyA. (d) The identification value of Il,zzA.
(a)(b)(c)(d)Table3 shows the true value Xr, identification value Xid, and relative error er of the load dynamic parameters obtained by the proposed method. The relative error er is defined as(37)er=Xr−XidXr×100%.Table 3
Real and identification inertia parameters of the payload.
ParemeterReal valueIdentification valueRelative error (%)Δz0.30.29780.73Il,xxA2525.943.76Il,yyA2524.751.00Il,zzA2528.2012.80


## 5. Conclusion

The proposed method has been successfully applied to a six-degree-of-freedom Gough–Stewart parallel manipulator for load dynamic parameters estimation. The 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.



---

*Source: 8816374-2020-11-27.xml*