Calibration Optimization of Kinematics and Dynamics for Delta Robot Driven by Integrated Joints in Machining Task

Abstract

For the application of Delta robots with a 3-R(RPaR) configuration in machining tasks, this paper constructed a 54-parameter kinematic error model and a simplified dynamic model incorporating an integrated joint’s position error and friction, respectively. Utilizing Singular Value Decomposition (SVD) of the Linear Model Coefficient Matrix (LMCM) and the coefficient chart, a criterion for identifiability of error components is established. For good identification results, the optimal measurement surface with Fourier series form is obtained using a combination of the Hook–Jeeves Direct Search Algorithm (DSA) and Inner Point Method (IPM). The friction coefficients and other dynamic parameters are obtained through fitting the integrated joint torque-angle pairs measured along specific trajectories using nonlinear least squares regression. The validation of the calibration process is conducted through simulations and experiments. The calibration results provide a foundation for the precise control of integrated joints and the high-precision motion of robots.

1. Introduction

In comparison to serial robots, parallel robots offer significant advantages, including enhanced stiffness and improved repeatability. Delta robots, as a highly successful configuration of parallel robots, have found extensive applications in various fields such as Pick-and-Place [1], 3D printing [2], medical fields [3], and machining [4]. For machining applications, Delta robots demand exceptional motion precision. Precision in robotic systems is typically achieved through two approaches: precise analysis during the design stage and error compensation following assembly. Error calibration and compensation, facilitated using algorithms and software, substantially enhance the precision of deployed robots [5,6].

Kinematic modeling serves as the cornerstone of robot error calibration. Various methods are employed, including the Denavit–Hartenberg (DH) method [7], exponential product [8], and vector loop method [9]. The vector loop method formulates closed-loop motion equations through position vectors at the link ends or joints within kinematic chains and is highly suitable for parallel robots. Vischer et al. [10] utilized the vector loop method to construct two kinematic models, “Model 54” and “Model 24”, for the 3-R(2-SS) configuration. “Model 54” considers geometric errors in all mechanical components except for the spherical joints; if the spatial parallelogram is assumed to remain perfect, “Model 54” can be simplified to “Model 24”. In “Model 24”, the direction of the Delta robot’s moving platform is assumed to be constant. Li et al. [11], also focusing on the 3-R(2-SS) configuration, demonstrated a significant influence of parallel mechanism errors on motion precision, especially in the Z-direction, with an increase in the tilt angle of the moving platform. Shen et al. [12], utilizing the variational method, studied the effect of error components (including the parallelogram mechanism) on the motion precision of Delta-RU robots, significantly enhancing the accuracy of the end effector through compensation. While scholars have recently introduced data-driven approaches for kinematic calibration, the model-based method continues to be predominant. This preference stems from its clear physical significance and the extensive measurement data demanded by the former [13].

Kinematic parameter identification is typically based on redundant equations formed using kinematic models and multiple measurement points. Considering noises in measurement data, the final parameter solution is obtained through nonlinear least squares solutions, with the Levenberg–Marquardt algorithm being the most commonly used algorithm [14]. Other probabilistic methods, such as maximum likelihood estimation and Kalman filtering, are less commonly used. Different parameter error components have varying degrees of influence on the total pose error, typically characterized using sensitivity analysis [15,16]. A local parameter sensitivity can be defined as the ratio of the sensitivity coefficient of a single error component to the combined one of all errors at a measurement point. A global sensitivity can be defined as the integration of local sensitivity over the robot workspace. Sensitivity analysis can identify the dominant error components of the robot. In addition to the calibration for link dimensions and assembly errors, scholars have also investigated the calibration of joint clearance. Kim et al. [17] discussed the relationship between the end effector pose error of the Stewart platform and the tolerance limits of its joints. Tannous et al. [18] conducted a sensitivity analysis on joint clearance using interval analysis methods, taking three manipulators as examples.

The problem of identifying dynamic parameters of robots involves several aspects, including specific identification parameters, dynamic models, excitation trajectories, and solution algorithms. Jan et al. [19] proposed a statistical framework for the dynamic parameter identification of multi-joint serial robots. The excitation trajectories were defined as finite Fourier series on each joint, with the optimal criterion being the uncertainty on the estimated parameters or a lower bound for it. The parameters identified included mass, inertia parameters, and friction parameters. In another study, the coefficients of the Fourier series were optimized to minimize the sensitivity of the identification to measurement disturbances [20]. Zafer et al. [21] used the least squares (LS) method and particle swarm optimization (PSO) technique to estimate distinct inertia parameters of the Staubli RX-60 robot, showing that the PSO method had a smaller identification error. Grotjahn et al. [22] investigated the identification of friction and rigid-body dynamics in parallel kinematic structures, achieving their separation via point-to-point motions. Han et al. [23] considered measurement errors in the process of optimizing excitation trajectories, integrating weighted least squares (WLS) and nonlinear error models, and iteratively reweighted least squares to unify the modeling process of friction and rigid-body dynamics. Dong et al. [24] first identified the friction parameters in the nonlinear Tustin friction model using the LS method during forward and reverse rotations, then estimated the remaining dynamic parameters of the Stäubli TX-90 robot using the symbiotic organisms search (SOS) algorithm. Song et al. [25] implemented the load parameter identification of parallel robots based on the Extended Kalman Filter. Ohno et al. [26] designed target trajectories for the detection of joint clearances based on actuation torque fluctuation. Diaz-Rodriguez et al. [27] developed a simplified dynamic model for a 3-DOF prismatic-revolute-spherical parallel manipulator to improve trajectory tracking precision and reduce the real-time computational burden. Abed Azad et al. [28] obtained the base inertial parameters through singular value decomposition, constructing an optimized path containing several harmonics within the desired bandwidth, achieving offline–online identification of the base inertial parameters.

Measurement processes can be categorized into two types: with or without external measurement devices. Without external measurement device support, the end effector of the robot is virtually constrained to a point, line, or surface [29,30], commonly applied in serial robots. The calibration process of parallel robots often employs external measurement devices, including laser trackers, laser interferometers, coordinate measurement machines (CMMs), vision-based measurement systems, ball-bar systems, etc. [10,31,32,33]. The selection of measurement points directly affects the quality of error parameter identification, leading to an optimization problem, typically aiming for an observation metric as the optimization objective along with constraints. Specific solving methods include local search methods represented by gradient descent and global search algorithms represented by meta-heuristic algorithms. Daney et al. [34] proposed a constrained optimization method for iteratively selecting measurement poses to minimize indicators based on the Jacobi matrix, utilizing the tabu search metaheuristic algorithm to avoid local optima.

This study tackles the relatively unexplored challenge of kinematic error modeling and dynamic calibration for the 3-R(RPaR) configuration Delta robot integrated into a self-built five-axis blade grinding device. It specifically investigates how geometric errors affect the poses of parallelogram mechanism while considering the constraints imposed by the robot’s DoF. Departing from conventional measurement point selection methods, this research reframes the issue as an optimization problem aimed at sha** the measurement workpiece surface, under the assumption that pose errors can be inferred from surface machining errors. Given the complexity of existing parameter sensitivity approaches, the paper seeks to establish criteria for parameter identifiability through the singular value decomposition (SVD) of the LMCM. For machining tasks with smooth and uniform motion, a simplified dynamic model with friction coefficients and other dynamic parameters is established and identified. The objective of this study is to devise a Delta robot calibration methodology tailored for machining environments.

The structure of the subsequent sections is as follows: Section 2 provides an introduction to the blade grinding device and details the kinematic and dynamic models of the Delta robot. Section 3 examines the influence of geometric errors on the parallel mechanism and describes the development of the kinematic error model tailored for the 3-R(RPaR) configuration. In Section 4, the identifiability criteria for error components are proposed, along with the definition of periodic error measurement surfaces and the optimization models under constraints. Section 5 delves into solution algorithms for the optimization problem and validates the kinematic calibration process through simulation. Section 6 validates the dynamic calibration method through experiments. The paper concludes with a summary of the findings and final remarks.

2. Kinematic and Dynamic Modeling

The blade grinding equipment, depicted in Figure 1, includes an optical platform, an aluminum alloy frame, a grinding head unit with a C and B axis, and a Delta robot unit suspended at the top of the frame. The Delta robot’s static platform is secured to the frame’s top, while a six-axis force sensor is mounted on its moving platform and connected to the workpiece. The C axis rotates vertically on the optical platform, and the B axis is perpendicular to its rotation axis. A spindle, mounted perpendicular to the B axis, holds an abrasive flap wheel via an extension rod. In this paper, only the Delta robot unit of the equipment is addressed, with the angular positions of the C and B axis held fixed during the calibration process.

The Delta robot is composed of a stationary platform (base), a moving platform, and three interconnected kinematic chains. This study focuses on the 3-R(RPaR) configuration, illustrated in Figure 2. Each chain connects the base to an active arm via an active revolute joint denoted as A_i (i = 1, 2, 3). The active arm is linked to a parallelogram mechanism, consisting of four revolute joints (D_i, E_i, F_i, G_i) and four links, through another passive revolute joint labeled as B_i. The moving platform connects to the parallel mechanism through a passive revolute joint denoted as C_i. The three chains are evenly distributed around the stationary platform, spaced 120° apart, with the Z₀ axis defined as the central axis. The origin S of the stationary coordinate system {W₀} is at the intersection of the plane containing A₁A₂A₃ with the X₀ axis, and Y₀ = Z₀ × X₀. Coordinate systems {W_i,1} and {W_i,2} are established at A_i, fixed to the stationary platform and the active arm, respectively. Z_i,1 is oriented along the rotation axis of joint A_i, while Y_i,1 is parallel to Z₀ but in the opposite direction. Z_i,2 is aligned with Z_i,1, and X_i,2 is oriented along the vector A_iB_i. The angle between X_i,2 and X_i,1 is denoted as θ_i, representing the rotation angle of the active arm. The unit vector of D_iE_i is denoted as Z_i,3, and the unit vector of C_iG_i is denoted as Z_i,4. Z_i,4 also serves as the Z-axis in {W_i,4}, fixed to point C_i and the moving platform. The direction of Y_i,4 is consistent with Z₀. The lengths SA_i and TC_i represent the radius r₁ of the stationary platform and the radius r₂ of the moving platform, respectively. l_p and l_n denote the lengths of active arm A_iB_i and the vector B_iC_i, respectively, defining the four basic parameters of the Delta robot. Φ₁ = 0°, Φ₂ = 120°, and Φ₃ = 240° indicate the angles between the three active arm rotation planes and Z₀X₀ plane for chains 1, 2, and 3, respectively. Under ideal conditions, the moving platform only undergoes translation, maintaining a fixed orientation. A local coordinate system {W₅} is established with T as the origin, aligned with {W₀}. {W₇}, {W₈} and {W₉} in Figure 2 are, respectively, associated with the C axis, B axis, and the grinding tool, which are not involved in the subsequent sections.

T is the center of C₁, C₂ and C₃ and is denoted as ⁰T (x_t, y_t, z_t) and ^i,1T(^i,1x_t, ^i,1y_t, ^i,1z_t) in {W₀} and {W_i,1}, respectively. In {W_i,1},

$$B_i{C}_i=T{C}_i-\left(S{A}_i+A_i{B}_i\right),$$

(1)

where $S{A}_i={\left[-r_1,0,0\right]}^T$, $A_i{B}_i={[l_p\cos {\theta }_i,l_p\sin {\theta }_i,0]}^T$, $T{C}_i={\left[{}^{i,1}x_t,{}^{i,1}y_t,{}^{i,1}z_t\right]}^T+{\left[r_2,0,0\right]}^T$. Substituted into Equation (1),

$$B_i{C}_i={\left[{}^{i,1}x_t+r_2,{}^{i,1}y_t,{}^{i,1}z_t\right]}^T-{\left[l_p\cos {\theta }_i,l_p\sin {\theta }_i,0\right]}^T.$$

(2)

B_iC_i always satisfies:

$$‖B_i{C}_i‖=l_n \mathrm{or} B_i{C}_i\cdot B_i{C}_i=l_n^2.$$

(3)

The relationship of T coordinates in {W₀} and {W_i,1} can be represented as the following:

$${\left[{}^{i,1}x_t,{}^{i,1}y_t,{}^{i,1}z_t\right]}^T=Rot(X,-{90}^{\circ })Rot(Z,-{\Phi }_i){\left[x_t,y_t,z_t\right]}^T-{\left[r_1,0,0\right]}^T.$$

(4)

where Rot(axis, angle) represents the rotation matrix around axis by angle. The fundamental kinematic equations of the Delta robot are formed by substituting Equation (3) for i = 1,2,3, respectively. Given θ₁, θ₂, and θ₃, if the goal is to determine x_t, y_t, z_t, it constitutes a forward kinematics problem. Conversely, it represents an inverse kinematics problem. By substituting Equation (2) into Equation (3),

$${\left(l_p\cos {\theta }_i-{}^{i,1}x_t-r_2\right)}^2+{\left(l_p\sin {\theta }_i-{}^{i,1}y_t\right)}^2+{\left({}^{i,1}z_t\right)}^2=l_n^2.$$

(5)

Beginning with the inverse kinematics, Equation (5) yields the following:

$$\cos {\psi }_i\cos {\theta }_i+\sin {\psi }_i\sin {\theta }_i=\cos \left({\theta }_i+{\psi }_i\right)=M_i/2l_p{N}_i.$$

(6)

where ${\psi }_i=\mathrm{atan}2\left({}^{i,1}y_t,{}^{i,1}x_t+r_2\right)$, $M_i=l_p^2-l_n^2+{\left({}^{i,1}x_t+r_2\right)}^2+{}^{i,1}y_t^2+{}^{i,1}z_t^2$, $N_i=\sqrt{{\left({}^{i,1}x_t+r_2\right)}^2+{}^{i,1}y_t^2}$. Solving Equation (6),

$${\theta }_i=\mathrm{acos}\left(M_i/2l_p{N}_i\right)+{\psi }_i, \mathrm{or} {\theta }_i=-\mathrm{acos}\left(M_i/2l_p{N}_i\right)+{\psi }_i.$$

(7)

Proceeding to obtain the forward solution, similarly derived from Equation (5),

$$x_t^2+y_t^2+z_t^2+a_i{x}_t+b_i{y}_t+c_i{z}_t+d_i=0.$$

(8)

where $a_i=-2l_p\cos {\theta }_i\cos {\Phi }_i-2\Delta r\cos {\Phi }_i$, $b_i=-2l_p\cos {\theta }_i\sin {\Phi }_i-2\Delta r\sin {\Phi }_i$, $c_i=-2l_p\sin {\theta }_i$, $d_i=l_p^2+\Delta r^2+2\Delta r{l}_p\cos {\theta }_i-l_n^2$, and Δr = r₁ − r₂. The solutions for x_t, y_t and z_t are obtained through the resolution of quadratic equations, the details of which are omitted here for brevity. In both forward and inverse solutions, the selection among multiple solutions is necessary. The viable solution must adhere to inequality (9).

$$\left({}^0\left(A_i{B}_i\right)\times {}^0\left(B_i{C}_i\right)\right)\cdot {}^{0}Z_{i,3}>0$$

(9)

where the left superscript “0” denotes the vector value in {W₀}, and ${}^{0}Z_{i,3}=Rot({\rm Z},f_i)Rot(X,{90}^{\circ }){[0,0,1]}^T$, ${}^{0}A_i=Rot(Z,f_i)Rot(X,{90}^{\circ }){[r_1,0,0]}^T$, ${}^{0}B_i=Rot(Z,f_i)Rot(X,{90}^{\circ }){[l_p\cos {\theta }_i+r_1,l_p\sin {\theta }_i,0]}^T$, ${}^{0}C_i={}^{0}T_i+Rot(Z,f_i)Rot(X,{90}^{\circ }){[r_2,0,0]}^T$.

In accordance with the principle of virtual work, the virtual work carried out using all non-inertial forces should equal the virtual work carried out using inertial forces. This principle is applied to driving integrated joints:

$$\left(\Gamma +{\Gamma }_c+J^T{F}_G+{\Gamma }_{CM}\right)\delta \Theta =\left(I\ddot{\Theta }+J^T{F}_A\right)\delta \Theta ,$$

(10)

where Γ represents the motor torque, Γ_c represents frictional torque, I denotes the moment of inertia of the integrated joints (including active arms), J signifies the Jacobian matrix, F_A stands for the inertial force generated by the moving platform mass (m_e), F_G represents the gravitational force acting on m_e, Γ_CM accounts for the gravitational moment produced by the unbalanced mass of the integrated joints (including active arms), δΘ represents virtual displacement, and $\ddot{\Theta }$ denotes joint acceleration. Here, it is assumed that the mass of the parallelogram linkage segment is negligible, or it can be separately accounted for at either end of integrated joints or the moving platform.

3. Kinematic Error Modeling for the 3-R(RPaR) Configuration

Geometric errors of the Delta robot unit result in pose deviations of the workpiece mounted on the moving platform during grinding operations. Under error conditions, the relationships between coordinate systems are adjusted. Variables labeled with the right superscript “0” represent their nominal values.

3.1. Transfer Relationships between Coordinate Systems or Points

The orientation and origin of {W_i,1} in {W₀} are expressed as Equations (11) and (12), respectively.

$${}_{i,1}{}^{0}R=Rot\left(Z,{\phi }_i\right)Rot\left(X,{\gamma }_i\right)$$

(11)

$${}^{0}A_i={\left[x_{ai},y_{ai},z_{ai}\right]}^T$$

(12)

where the nominal values of φ_i, γ_i, x_ai, y_ai, z_ai are denoted as ${\phi }_i^0$, ${\gamma }_i^0$, $x_{ai}^0$, $y_{ai}^0$, $z_{ai}^0$, respectively, and ${\phi }_i^0={\Phi }_i$, ${\gamma }_i^0={90}^{\circ }$, $x_{ai}^0=r_1\cos {\Phi }_i$, $y_{ai}^0=r_1\sin {\Phi }_i$, $z_{ai}^0=0$.

The origin of {W_i,2} is the same as that of {W_i,1} and the relative orientation is the following:

$${}_{i,2}{}^{i,1}R=Rot\left(Z,{\theta }_i\right),$$

(13)

where the nominal value of ${\theta }_i$ is ${\theta }_i^0$. $\Delta {\theta }_i={\theta }_i-{\theta }_i^0$ is the zero position offset of the active revolute joint i.

The coordinate of B_i in {W_i,2} is as follows:

$${}^{i,2}B_i={\left[l_{pi},0,0\right]}^T,$$

(14)

where the nominal value of the active arm length $l_{pi}$ is $l_p$.

The relationship between ${}^{i,2}Z_{i,3}$ and ${}^{i,2}Z_{i,2}$ is as follows:

$${}^{i,2}Z_{i,3}=Rot\left(X,\Delta {\alpha }_i\right)Rot\left(Y,\Delta {\beta }_i\right){}^{i,2}Z_{i,2},$$

(15)

where Δα_i and Δβ_i are approximate infinitesimal errors, Z_i,2 and Z_i,3 are approximately parallel under error conditions.

The orientation and origin of {W_i,4} in {W₅} are the following:

$${}_{i,4}{}^{5}R=Rot\left(Z,{\epsilon }_i\right)Rot\left(X,{\tau }_i\right),$$

(16)

$${}^{5}C_i={\left[x_{ci},y_{ci},z_{ci}\right]}^T,$$

(17)

where the nominal values of ${\epsilon }_i$, ${\tau }_i$, x_ci, y_ci, z_ci are denoted as ${\epsilon }_i^0$, ${\tau }_i^0$, $x_{ci}^0$, $y_{ci}^0$, $z_{ci}^0$, respectively, and ${\epsilon }_i^0={\Phi }_i$, ${\tau }_i^0={90}^{\circ }$, $x_{ci}^0=r_2\cos {\Phi }_i$, $y_{ci}^0=r_2\sin {\Phi }_i$, $z_{ci}^0=0$.

The orientation and origin of {W₅} in {W₀} are the following:

$${}_5{}^{0}R=Rot\left(X,{\alpha }_5\right)Rot\left(Y,{\beta }_5\right)Rot\left(Z,{\gamma }_5\right),$$

(18)

$${}^{0}T={\left[x_t,y_t,z_t\right]}^T.$$

(19)

3.2. The Influence of Parallelogram Mechanism Errors

As shown in Figure 3, in the part of the parallelogram mechanism, the nominal points of D_i, E_i, F_i and G_i are denoted as $D_i^0$, $E_i^0$, $F_i^0$ and $G_i^0$, respectively. The following relations are satisfied: $‖D_i{E}_i‖=l_{bi}$, $‖F_i{G}_i‖=l_{ci}$, $\Delta l_{bi}=l_{bi}-l_b^0$, $2\Delta l_{ci}=l_{ci}-l_{bi}$, $‖D_i^0{F}_i^0‖=‖E_i^0{G}_i^0‖=l_n^0$, $‖D_i{F}_i‖=l_{ni}+\Delta l_i$, $‖E_i{G}_i‖=l_{ni}-\Delta l_i$, $\Delta l_{ni}=l_{ni}-l_n^0$. In the RPaR configuration, even if there are errors Δl_bi and Δl_ci, it can be assumed that B_i and C_i are located at the midpoint of D_iE_i and F_iG_i, respectively. Establishing a coordinate system {W_i,3} with B_i as the origin, where its Z-axis is Z_i,3, the X axis is perpendicular to the plane D_iE_iG_iF_i and points outward from the paper, denoted as X_i,3, and Y_i,3= Z_i,3 × X_i,3. The unit vector along B_iC_i is denoted as ${Y^{\prime }}_{i,3}$, ${Z^{\prime }}_{i,3}=X_{i,3}\times {Y^{\prime }}_{i,3}$. The angle between ${Y^{\prime }}_{i,3}$ and $Y_{i,3}$ is ξ_i.

$${\xi }_i=\mathrm{asin}‖Y_{i,3}\times {Y^{\prime }}_{i,3}‖$$

(20)

Due to the presence of Δl_ci and Δl_i, an angle Δξ_i exists between E_iG_i and ${Y^{\prime }}_{i,3}$, and an angle Δσ_i exists between Z_i,4 and Z_i,3. In Figure 3, the new point obtained by rotating G_i around C_i by Δσ_i is denoted as H_i. C_iH_i has the same direction as Z_i,3. ${G^{\prime }}_i$ lies on the line C_iH_i, and $E_i{G^{\prime }}_i$ is parallel to ${Y^{\prime }}_{i,3}$. The intersection point of $E_i{G^{\prime }}_i$ and F_iG_i is ${H^{\prime }}_i$. In the area enclosed by ${H^{\prime }}_i{G}_i{H}_i{G^{\prime }}_i$, assuming errors Δl_ci, Δσ_i, Δl_i, Δξ_i are small, and $\sin \Delta {\sigma }_i\approx \Delta {\sigma }_i$, $\sin \Delta {\xi }_i\approx \Delta {\xi }_i$, then

$$\Delta l_{ci}{Z}_{i,3}-\Delta {\sigma }_i{l}_{ci}{Y}_{i,3}=-\Delta l_i{Y^{\prime }}_{i,3}+\Delta {\xi }_i\left(l_{ni}-\Delta l_i\right){Z^{\prime }}_{i,3}.$$

(21)

Equation (21) depicts the extension of $C_i{G^{\prime }}_i$ by the distance Δl_ci to reach point H_i, followed by the rotation of H_i around C_i by Δσ_i to reach point G_i on the left side. Conversely, on the right side, it illustrates the shortening of l_ni by Δl_i to form segment $E_i{H^{\prime }}_i$, followed by the rotation of ${H^{\prime }}_i$ around E_i by Δξ_i to reach point G_i. Considering the relationship between ${Y^{\prime }}_{i,3}$ and $Y_{i,3}$, as well as ${Z^{\prime }}_{i,3}$ and $Z_{i,3}$, and $\Delta {\sigma }_i\Delta l_{ci}\approx 0$, $\Delta {\xi }_i\Delta l_i\approx 0$, it is derived that,

$$\Delta l_{ci}{Z}_{i,3}-\Delta {\sigma }_i{l}_{bi}{Y}_{i,3}=-\Delta l_{i}Rot\left(X,{\xi }_i\right)Y_{i,3}+\Delta {\xi }_i{l}_{ni}Rot\left(X,{\xi }_i\right)Z_{i,3},$$

(22)

where Rot(X, ξ_i) represents the rotation matrix around X axis by ξ_i. Equations (23) and (24) can be derived along Y_i,3 and Z_i,3.

$$\Delta {\sigma }_i{l}_{bi}=\Delta l_i\cos {\xi }_i+\Delta {\xi }_i{l}_{ni}\sin {\xi }_i,$$

(23)

$$\Delta l_{ci}=-\Delta l_i\sin {\xi }_i+\Delta {\xi }_i{l}_{ni}\cos {\xi }_i.$$

(24)

Therefore, it can be solved that,

$$\Delta {\sigma }_i=\frac{\Delta l_i+\Delta l_{ci}\sin {\xi }_i}{l_{bi}\cos {\xi }_i}.$$

(25)

Equation (25) reveals that the deviations Δl_i of links D_iF_i and E_iG_i, and the Δl_ci of links D_iE_i and F_iG_i cause relative pose changes in link F_iG_i. The value of ξ_i also affects Δσ_i when Δl_i and Δl_ci remain constant. Notably, l_bi and Δσ_i exhibit an inverse relationship, where larger l_bi values mitigate the impact of installation and manufacturing errors on link poses.

Although this analysis focuses on point G_i, a similar approach can be applied to point F_i. Assuming small Δl_ci, Δσ_i, Δl_i, Δξ_i, Equation (25) still holds. Based on these considerations, the Delta robot necessitates the calibration of 54 parameters (

Table 1. Parameters for calibration in Delta robot.

Index	Parameter	Meaning
1	$x_{ai}$	the X component of Ai coordinate in {W0}
2	$y_{ai}$	the Y component of Ai coordinate in {W0}
3	$z_{ai}$	the Z component of Ai coordinate in {W0}
4	${\phi }_i$	the Z rotation component transforming from {W0} to {Wi,1}
5	${\gamma }_i$	the X rotation component transforming from {W0} to {Wi,1}
6	${\theta }_i$	the rotation angle of the active revolute joint i
7	$l_{pi}$	the length of the active arm link
8	$∆{\alpha }_i$	the X rotation component transforming from i,2Zi,2 to i,2Zi,3
9	$∆{\beta }_i$	the Y rotation component transforming from i,2Zi,2 to i,2Zi,3
10	$l_{ni}$	the length of the vector BiCi within parallelogram mechanism
11	$∆l_i$	the lengthening or shortening amount of the links DiFi and EiGi relative to lni
12	$l_{bi}$	the length of link DiEi
13	$∆l_{ci}$	the half-deviation between the lengths of FiGi and DiEi
14	$x_{ci}$	the X component of Ci coordinate in {W5}
15	$y_{ci}$	the Y component of Ci coordinate in {W5}
16	$z_{ci}$	the Z component of Ci coordinate in {W5}
17	${\epsilon }_i$	the Z rotation component transforming from {W5} to {Wi,4}
18	${\tau }_i$	the X rotation component transforming from {W5} to {Wi,4}

), all of which can be assumed to have small errors.

3.3. The Kinematic Equations under Error Conditions

Under error conditions, the kinematic chain of the Delta robot satisfies the following equations:

$$‖B_i{C}_i‖=l_{ni},$$

(26)

$$‖Z_{i,3}\times Z_{i,4}‖=\Delta {\sigma }_i.$$

(27)

Equations (26) and (27) can be further expressed as follows:

$$‖\left({}^{0}A_i+{}_{i,1}{}^{0}R{}_{i,2}{}^{i,1}R{}^{i,2}B_i\right)-\left({}^{0}T_i+{}_5{}^{0}R{}^{5}C_i\right)‖=l_{ni},$$

(28)

$$‖\left({}_{i,1}{}^{0}R{}_{i,2}{}^{i,1}R\cdot Rot\left(X,\Delta {\alpha }_i\right)Rot\left(Y,\Delta {\beta }_i\right){}^{i,2}Z_{i,2}\right)\times \left({}_5{}^{0}R{}^{i,4}R_5\right){}^{i,4}Z_{i,4}‖=\Delta {\sigma }_i.$$

(29)

Expand Equations (28) and (29), and the results can be formally represented as follows:

$$f_i\left(x_{ai},y_{ai},z_{ai},{\phi }_i,{\gamma }_i,{\theta }_i,l_p,l_n,x_{ci},y_{ci},z_{ci},x_t,y_t,z_t,{\alpha }_5,{\beta }_5,{\gamma }_5\right)=0,$$

(30)

$$g_i\left({\phi }_i,{\gamma }_i,\Delta {\alpha }_i,\Delta {\beta }_i,{\epsilon }_i,{\tau }_i,\Delta {\sigma }_i,{\alpha }_5,{\beta }_5,{\gamma }_5\right)=0.$$

(31)

When the Delta robot with the 3-R(RPaR) configuration satisfies nominal geometric parameters, according to the modified Chebyshev–Grübler–Kutzbach formula [35], the mechanism’s DoF is as follows:

$$M_{3-R\left(RPaR\right)}=6\left(n-g-1\right)+{\displaystyle \sum _{q=1}^g{f}_q}+v-\kappa =6\left(17-21-1\right)+21+12-0=3,$$

(32)

where n is the number of bodies, g is the number of kinematic pairs, f_q is the DoF of the qth joint, and κ represents the internal mobilities or idle DoF. ν = ν_PM + ν_CJ = 3 + 9 = 12; here, ν_PM is the number of overconstraints assembling the kinematic chains to a parallel manipulator and ν_CJ is the number of overconstraints within the complex joints.

In Figure 3, deviations in parameter errors may lead to non-parallelism between F_iG_i and D_iE_i. Connecting the three kinematic chains to the moving and stationary platforms alters the overconstraint situation, resulting in a v_PM value of 0. Thus, v = ν_PM + ν_CJ = 0 + 9 = 9, and M_3-R(RPaR) = 0, indicating zero DoF for the Delta robot. Forcefully driving all three active arms induces deformation in the kinematic chains, creating surplus DoF. This paper defines structural deformation in line with the characteristics of the 3-R(RPaR) structure and referencing motion under error in 3-R(2-SS) as the rotation of link F_iG_i (Figure 3) about the Y_i,3 axis passing C_i. To clarify this process, in the planar state depicted in Figure 3, F_iG_i is redefined as $F_i^{″}{G}_i^{″}$ in Figure 4. Upon rotating $F_i^{″}{G}_i^{″}$ about the Y_i,3 axis by Δη_i, the final F_iG_i is obtained.

The pose relationship described by Equation (21) is redefined as follows:

$$Rot({}^{0}Y_{i,3},\Delta {\eta }_i)Rot({}^{0}X_{i,3},\Delta {\sigma }_i){}_{i,1}{}^{0}R{}_{i,2}{}^{i,1}R\cdot Rot\left(X,\Delta {\alpha }_i\right)Rot\left(Y,\Delta {\beta }_i\right){}^{i,2}Z_{i,2}=\left({}_5{}^{0}R{}_{i,4}{}^{5}R\right){}^{i,4}Z_{i,4},$$

(33)

where Rot(⁰Y_i,3, Δη_i) represents the rotation matrix around ⁰Y_i,3 by Δη_i and Rot(⁰X_i,3, Δσ_i) represents the rotation matrix around ⁰X_i,3 by Δσ_i. Both sides of Equation (33) represent unit vectors. Given the small Z-component of these unit vectors in {W₀} (attributable solely to errors), this paper formulates equations based on their X and Y components. With three chains (i = 1, 2, 3), Equation (33) yields six independent equations, formally expressed as follows:

$$q_{xi}\left(x_{ai},y_{ai},z_{ai},{\phi }_i,{\gamma }_i,{\theta }_i,l_{pi},\Delta {\alpha }_i,\Delta {\beta }_i,l_{ni},\Delta {\sigma }_i,x_{ci},y_{ci},z_{ci},{\epsilon }_i,{\tau }_i,\Delta {\eta }_i,x_t,y_t,z_t,{\alpha }_5,{\beta }_5,{\gamma }_5\right)=0,$$

(34)

$$q_{yi}\left(x_{ai},y_{ai},z_{ai},{\phi }_i,{\gamma }_i,{\theta }_i,l_{pi},\Delta {\alpha }_i,\Delta {\beta }_i,l_{ni},\Delta {\sigma }_i,x_{ci},y_{ci},z_{ci},{\epsilon }_i,{\tau }_i,\Delta {\eta }_i,x_t,y_t,z_t,{\alpha }_5,{\beta }_5,{\gamma }_5\right)=0.$$

(35)

Equations (30), (34) and (35) collectively depict the kinematic equations of the Delta robot in the 3-R(RPaR) configuration under error conditions. These equations entail multiple nonlinear terms involving trigonometric functions and higher-order factors, rendering the solution process cumbersome. To streamline these equations, parameters are initially expressed as the sum of nominal values and errors, as presented below:

$$x_{ai}=x_{ai}^0+\Delta x_{ai}, y_{ai}=y_{ai}^0+\Delta y_{ai}, z_{ai}=z_{ai}^0+\Delta z_{ai},$$

(36)

$${\phi }_i={\phi }_i^0+\Delta {\phi }_i, {\gamma }_i={\gamma }_i^0+\Delta {\gamma }_i, {\theta }_i={\theta }_i^0+\Delta {\theta }_i,$$

(37)

$$l_{pi}=l_p+\Delta l_{pi}, l_{ni}=l_n+\Delta l_{ni},$$

(38)

$$x_{ci}=x_{ci}^0+\Delta x_{ci}, y_{ci}=y_{ci}^0+\Delta y_{ci}, z_{ci}=z_{ci}^0+\Delta z_{ci},$$

(39)

$${\epsilon }_i={\epsilon }_i^0+\Delta {\epsilon }_i, {\tau }_i={\tau }_i^0+\Delta {\tau }_i.$$

(40)

The right superscript “0” denotes the nominal value of a parameter, while the prefix Δ represents the parameter’s error. For example, $x_{ai}^0$ and Δx_ai represent the nominal and error values of x_ai, respectively. All parameter error values are considered small quantities. Taking φ_i as an example, the simplification is as follows:

$$\cos \left(\Delta {\phi }_i\right)\approx 1, \sin \left(\Delta {\phi }_i\right)\approx \Delta {\phi }_i,$$

(41)

$$\cos \left({\phi }_i^0+\Delta {\phi }_i\right)\approx \cos {\phi }_i^0-\Delta {\phi }_i\sin {\phi }_i^0,$$

(42)

$$\sin \left({\phi }_i^0+\Delta {\phi }_i\right)\approx \sin {\phi }_i^0+\Delta {\phi }_i\cos {\phi }_i^0.$$

(43)

The other angle parameters γ_i, θ_i, ε_i, τ_i, follow a similar pattern of φ_i. Additionally, all second-order and higher small quantities are omitted.

$$\Delta x_{ai}\cdot \Delta x_{ai}\approx 0, \Delta x_{ai}\cdot \Delta x_{ai}\cdot \Delta x_{ai}\approx 0,\dots$$

(44)

$$\Delta x_{ai}\cdot \Delta y_{ai}\approx 0, \dots , \Delta x_{ai}\cdot \Delta {\phi }_i\approx 0,\dots$$

(45)

Equations (30), (34) and (35) can be expressed in matrix form after simplification as follows:

$$F_{p}v+h_p=0,$$

(46)

$$G_{ox}v+h_{ox}=0,$$

(47)

$$G_{oy}v+h_{oy}=0,$$

(48)

where $F_p=\left[\begin{array}{c}{f}_{1,1},\cdots ,f_{1,j},\cdots ,f_{1,18}\\ 0\\ 0\end{array}\begin{array}{c}0\\ f_{2,1},\cdots ,f_{2,j},\cdots ,f_{2,18}\\ 0\end{array}\begin{array}{c}0\\ 0\\ f_{3,1},\cdots ,f_{3,j},\cdots ,f_{3,18}\end{array}\right]$, $h_p=\left[\begin{array}{c}{h}_{p1}\\ h_{p2}\\ h_{p3}\end{array}\right]$, $f_{i,j}=f_{i,j}\left(x_{ai}^0,y_{ai}^0,z_{ai}^0,{\phi }_i^0,{\gamma }_i^0,{\theta }_i^0,l_p,l_n,{\epsilon }_i^0,{\tau }_i^0,x_{ci}^0,y_{ci}^0,z_{ci}^0,x_t,y_t,z_t,{\alpha }_5,{\beta }_5,{\gamma }_5\right)$, $h_{pi}=h_{pi}{\left(x_{ai}^0,y_{ai}^0,z_{ai}^0,{\phi }_i^0,{\gamma }_i^0,{\theta }_i^0,l_p,l_n,{\epsilon }_i^0,{\tau }_i^0,x_{ci}^0,y_{ci}^0,z_{ci}^0,x_t,y_t,z_t,{\alpha }_5,{\beta }_5,{\gamma }_5\right)}_{18\times 1}$, $v={\left[v_{c1},v_{c2},v_{c3}\right]}^T$, $v_{ci}={\left[\Delta x_{ai},\Delta y_{ai},\Delta z_{ai},\Delta {\phi }_i,\Delta {\gamma }_i,\Delta {\theta }_i,\Delta l_{pi},\Delta {\alpha }_i,\Delta {\beta }_i,\Delta l_{ni},\Delta l_i,\Delta l_{bi},\Delta l_{ci},{\epsilon }_i,{\tau }_i,\Delta x_{ci},\Delta y_{ci},\Delta z_{ci}\right]}^T$, $G_{ox}=\left[\begin{array}{c}{g}_{x1,1},\cdots ,g_{x1,j},\cdots ,g_{x1,18}\\ 0\\ 0\end{array}\begin{array}{c}0\\ g_{x2,1},\cdots ,g_{x2,j},\cdots ,g_{x2,18}\\ 0\end{array}\begin{array}{c}0\\ 0\\ g_{x3,1},\cdots ,g_{x3,j},\cdots ,g_{x3,18}\end{array}\right]$, $h_{ox}=\left[\begin{array}{c}{h}_{ox1}\\ h_{ox2}\\ h_{ox3}\end{array}\right]$, $g_{xi,j}=g_{xi,j}\left(x_{ai}^0,y_{ai}^0,z_{ai}^0,{\phi }_i^0,{\gamma }_i^0,{\theta }_i^0,l_p,l_n,l_b,x_{ci}^0,y_{ci}^0,z_{ci}^0,{\epsilon }_i^0,{\tau }_i^0,\Delta {\eta }_i,x_t,y_t,z_t,{\alpha }_5,{\beta }_5,{\gamma }_5\right)$, $h_{oxj}=h_{oxj}{\left(x_{ai}^0,y_{ai}^0,z_{ai}^0,{\phi }_i^0,{\gamma }_i^0,{\theta }_i^0,l_p,l_n,l_b,x_{ci}^0,y_{ci}^0,z_{ci}^0,{\epsilon }_i^0,{\tau }_i^0,\Delta {\eta }_i,x_t,y_t,z_t,{\alpha }_5,{\beta }_5,{\gamma }_5\right)}_{18\times 1}$, $G_{oy}=\left[\begin{array}{c}{g}_{y1,1},\cdots ,g_{y1,j},\cdots ,g_{y1,18}\\ 0\\ 0\end{array}\begin{array}{c}0\\ g_{y2,1},\cdots ,g_{y2,j},\cdots ,g_{y2,18}\\ 0\end{array}\begin{array}{c}0\\ 0\\ g_{y3,1},\cdots ,g_{y3,j},\cdots ,g_{y3,18}\end{array}\right]$, $h_{oy}=\left[\begin{array}{c}{h}_{oy1}\\ h_{oy2}\\ h_{oy3}\end{array}\right]$, $g_{yi,j}=g_{yi,j}\left(x_{ai}^0,y_{ai}^0,z_{ai}^0,{\phi }_i^0,{\gamma }_i^0,{\theta }_i^0,l_p,l_n,l_b,x_{ci}^0,y_{ci}^0,z_{ci}^0,{\epsilon }_i^0,{\tau }_i^0,\Delta {\eta }_i,x_t,y_t,z_t,{\alpha }_5,{\beta }_5,{\gamma }_5\right)$, $h_{oyj}=h_{oyj}{\left(x_{ai}^0,y_{ai}^0,z_{ai}^0,{\phi }_i^0,{\gamma }_i^0,{\theta }_i^0,l_p,l_n,l_b,x_{ci}^0,y_{ci}^0,z_{ci}^0,{\epsilon }_i^0,{\tau }_i^0,\Delta {\eta }_i,x_t,y_t,z_t,{\alpha }_5,{\beta }_5,{\gamma }_5\right)}_{18\times 1}$.

Δη_i satisfies small quantity assumption. Using $\left\{\begin{array}{c}{\displaystyle \sum _{j=1}^{18}{g}_{xi,j}\cdot v_j+h_{oxi}}=0\\ {\displaystyle \sum _{j=1}^{18}{g}_{yi,j}\cdot v_j+h_{oyi}}=0\end{array}\right.$ to eliminate Δη_i, Equations (47) and (48) can be rearranged as:

$$G_{o}v+h_o=0,$$

(49)

where $G_o=\left[\begin{array}{c}{g}_{1,1},\cdots ,g_{1,j},\cdots ,g_{1,18}\\ 0\\ 0\end{array}\begin{array}{c}0\\ g_{2,1},\cdots ,g_{2,j},\cdots ,g_{2,18}\\ 0\end{array}\begin{array}{c}0\\ 0\\ g_{3,1},\cdots ,g_{3,j},\cdots ,g_{3,18}\end{array}\right]$, $h_o=\left[\begin{array}{c}{h}_{o1}\\ h_{o2}\\ h_{o3}\end{array}\right]$, $g_{i,j}=g_{i,j}\left(x_{ai}^0,y_{ai}^0,z_{ai}^0,{\phi }_i^0,{\gamma }_i^0,{\theta }_i^0,l_p,l_n,l_b,x_{ci}^0,y_{ci}^0,z_{ci}^0,{\epsilon }_i^0,{\tau }_i^0,x_t,y_t,z_t,{\alpha }_5,{\beta }_5,{\gamma }_5\right)$, $h_{oj}=h_{oj}{\left(x_{ai}^0,y_{ai}^0,z_{ai}^0,{\phi }_i^0,{\gamma }_i^0,{\theta }_i^0,l_p,l_n,l_b,x_{ci}^0,y_{ci}^0,z_{ci}^0,{\epsilon }_i^0,{\tau }_i^0,x_t,y_t,z_t,{\alpha }_5,{\beta }_5,{\gamma }_5\right)}_{18\times 1}$.

Equations (46) and (49) constitute the linear kinematic equations of the Delta robot with the 3-R(RPaR) configuration under error conditions.

$$\left[\begin{array}{c}{F}_p\\ G_o\end{array}\right]v+\left[\begin{array}{c}{h}_p\\ h_o\end{array}\right]=0.$$

(50)

In theory, utilizing Equations (28) and (33) to derive the partial derivatives of x_t, y_t, z_t, as well as α₅, β₅ and γ₅ with respect to each error component, yields a similar error model to Equation (50) containing the Jacobian matrix. However, since x_t, y_t, z_t and α₅, β₅, γ₅ are implicit in Equations (28) and (33), the solving process becomes relatively complex. This paper adopts the assumption of small errors and follows the principles of small quantity operations for a more intuitive process. Given that the error values of x_t, y_t, z_t and α₅, β₅, γ₅ consist of various error components and are relatively large, for enhanced accuracy, the relevant elements in matrices F_p and G_o remain expressed as functions of x_t, y_t, z_t and α₅, β₅, γ₅, rather than their corresponding linear error terms.

4. The Identification of Parameter Errors

4.1. Parameter Error Identifiability

Parameter identifiability comprises two primary aspects: the linear correlation between parameters, potentially impeding accurate individual parameter identification, and parameter sensitivity, indicating the relative impact of parameter errors on overall error. Traditional methods often struggle to effectively address both concerns, resulting in complex analyses. This study employs a simplified approach, the workspace random points + SVD method, to address these challenges. By selecting n random points within the Delta robot’s workspace, their actual poses are computed using Equation (50).

In Equation (25), $l_{bi}=l_{bi}^0+\Delta l_{bi}$; if Δl_bi is regarded as a first-order small quantity, it will be absorbed in calculations. Consequently, this study overlooks the impact of Δl_bi and does not account for it. In Equation (50), the component Δl_bi in the parameter error vector v is consistently treated as zero. To streamline the analysis, Δl_b1, Δl_b2 and Δl_b3 are omitted from v. Utilizing Equation (50), an equation set is formulated at each point, yielding the following:

$$Mv+b=0,$$

(51)

where M is the Linear Model Coefficient Matrix (LMCM). Since Δl_bi is excluded from v, M has a size of 6n × 51, and the bias b has a size of 6n × 1. Performing SVD on matrix M, Equation (51) can be rewritten as follows:

$$U\Sigma W^{T}v+b=0,$$

(52)

where U and W represent orthogonal matrices of size 6n × 6n and 51 × 51, respectively. u_i and w_i (for i = 1, …, 51) denote the row vectors of U and W, respectively. Σ is a 6n × 51 matrix, with the first 51 rows forming a diagonal matrix consisting of the singular values of M arranged in descending order. The remaining rows of Σ are all zeros. To examine the influence of singular values on the parameter solution, Equation (52) can be expressed in terms of u_i and w_i:

$${\lambda }_1{u}_1{w}_1^{T}v+{\lambda }_2{u}_2{w}_2^{T}v+\cdots +{\lambda }_{51}{u}_{51}{w}_{51}^{T}v+b=0.$$

(53)

The structure of Equation (53) indicates that the smaller ${\lambda }_i{u}_i{w}_i^{T}v$ is, the lower the impact on the parameter solution. The magnitude of ${\lambda }_i{u}_i{w}_i^{T}v$ is represented by a vector norm, and

$$‖{\lambda }_i{u}_i{w}_i^{T}v‖\le {\lambda }_i‖u_i‖‖w_i^T‖‖v‖={\lambda }_i‖v‖.$$

(54)

In Equation (54), it is evident that the magnitude of the term ${\lambda }_i{u}_i{w}_i^{T}v$ depends on the size of λ_i when v is given or held constant. When λ_i is sufficiently small, ${\lambda }_i{u}_i{w}_i^{T}v$ can be disregarded in Equation (53). Higher-order small terms are neglected during the derivation of Equation (50) due to the assumption of small errors. Therefore, the accuracy of equations built upon Equation (50) does not surpass the order of $‖v‖\cdot ‖v‖$. By comparison with ${\lambda }_i‖v‖$ in Equation (54), the following criterion can be deduced:

$${\lambda }_i<‖v‖.$$

(55)

When λ_i satisfies Equation (55), the corresponding parameter becomes unidentifiable. The relationship between λ_i and the parameters will be illustrated through examples in the simulation section.

The linearity correlation of parameters can be analyzed from the coefficients of the error parameters. For example, in Equation (46), the coefficient of Δx_ai is as follows:

$$\begin{array}{l}{f}_{i,1}=2x_{ai}-2x_t-2z_{ci}\sin \left({\beta }_5\right)-2x_{ci}\cos \left({\beta }_5\right)\cos \left({\gamma }_5\right)\\ \begin{array}{cc}& \end{array}+2y_{ci}\cos \left({\beta }_5\right)\sin \left({\gamma }_5\right)+2\cos \left({\phi }_i\right)\cos \left({\theta }_i\right)l_p-2\cos \left({\gamma }_i\right)l_p\sin \left({\phi }_i\right)\sin \left({\theta }_i\right)\end{array}.$$

(56)

The coefficient f_i,15 can be derived as well, but due to its lengthy equation, it is not listed here. Without assuming α₅, β₅, γ₅ as small quantities, there exists no fixed linear relationship between f_i,1 and f_i,15. However, by assuming α₅, β₅, γ₅ as small quantities and sin() = 0, cos() = 1, then

$$f_{i,1}=-f_{i,15}=2x_{ai}-2x_t-2x_{ci}+2l_p\cos \left({\phi }_i\right)\cos \left({\theta }_i\right)-2l_p\cos \left({\gamma }_i\right)\sin \left({\phi }_i\right)\sin \left({\theta }_i\right).$$

(57)

Similarly, the coefficients corresponding to Δy_ai and Δy_ci, as well as Δz_ai and Δz_ci, exhibit the same relationship. Given these linear relationships, Δx_ai, Δy_ai and Δz_ai can directly replace Δx_ai − Δx_ci, Δy_ai − Δy_ci, and Δz_ai − Δz_ci. This simplifies the dimensionality of v. Analyzing or determining linear relationships between parameters directly from lengthy coefficient expressions is often challenging. In the simulation section, these linear relationships will be illustrated graphically with examples.

4.2. Measurement Surface Optimization

This paper assumes that the pose error of the moving platform can be determined through measuring the surface form error of the machined workpiece, termed as the measurement surface. Assuming a fixed form of point acquisition on this surface, its shape directly influences the error parameter identification process. This section aims to identify the optimal measurement surface to enhance the identification process.

First, the measurement surface is parameterized. The coordinate system of the measurement surface is defined as {W₆}, with {W₆} being fixed relative to {W₅}. The origin of {W₆} is offset along the Z axis of {W₅} by a distance of H_q, with the axes of {W₆} parallel and opposite to those of {W₅}, respectively. H_q is the height of the measurement surface. Within {W₆}, the measurement surface Q(s, t) is defined. This surface can be extended periodically in both the s and t directions within a finite region, creating an infinitely free-from surface. Alternatively, Q(s, t) can be viewed as a single-period surface on an infinitely periodic surface.

$$Q\left(s,t\right)=Q_r\left(s\right)Q_h\left(t\right),$$

(58)

where Q_r(s) and Q_h(t) are both single-period functions, constructed using multiplication to form the Q(s, t) surface. This construction method guarantees that any iso-parameter curve, such as Q(s₀, t) = Q_r(s₀)Q_h(t) or Q(s, t₀) = Q_r(s)Q_h(t₀), remains periodic.

According to the Fourier series expansion theorem, Q_r(s) and Q_h(t) can be further represented as follows:

$$Q_r\left(s\right)=\frac{a_{u0}}{2}+{\displaystyle \sum _{m=1}^3\left(a_{um}\cos \left(3ms\right)+b_{um}\sin \left(3ms\right)\right)},$$

(59)

$$Q_h\left(t\right)=\frac{a_{v0}}{2}+{\displaystyle \sum _{n=1}^3\left(a_{vn}\cos \left(nt\right)+b_{vn}\sin \left(nt\right)\right)},$$

(60)

The Fourier series is an infinite sum of trigonometric functions, with high-frequency terms typically having a minimal impact on function values. Here, both m and n are set to 3, a sufficient choice for engineering purposes. Additionally, to capture the symmetry of the three kinematic chains in the XY plane of {W₅}, coefficients relative to s are represented in a 3m form. Consequently, the workpiece’s cross-sectional lines in the XY plane of {W₅} consist of three identical periodic curves, each spanning 120°.

Subsequently, the optimization problem formulation is established, with Equation (50) serving as the core of parameter identification. At multiple measurement points, an over-constrained equation set can be derived from Equation (50), as detailed in Equation (54) in Section 4.1. The quality of solving Equation (51) is most aptly reflected in the condition number of M. Hence, the optimization problem in this paper is defined as follows:

$$\left\{\begin{array}{cc}\begin{array}{c}\\ s.t.\\ \\ \end{array}& \begin{array}{c}\min f\left(v_{coff}\right)\\ R_q<R_{w\max }\\ R_q>R_{w\min }\\ H_q=H_{wm}\end{array}\end{array}\right.,$$

(61)

where $v_{coff}=[a_{s0},a_{s1},a_{s2},a_{s3},b_{s1},b_{s2},b_{s3},a_{t0},a_{t1},a_{t2},a_{t3},b_{t1},b_{t2},b_{t3}]$, $f\left(v_{coff}\right)=cond\left(M\right)$, and cond(M) denotes the condition number of M in Equation (51). R_q denotes the maximum radius of the measurement surface in the X₅Y₅ plane. R_q must be smaller than the maximum radius R_wmax of the cylindrical workspace of the Delta robot. Simultaneously, R_q should not be excessively small, as this could compromise the strength of the workpiece, leading to fractures, or result in deep grooves that are challenging to measure. The lower limit of R_q is set to R_wmin. The height H_q of the measurement surface is equal to the maximum height H_wmax of the robot cylindrical workspace by replacing t in Equation (60) with $2\pi t/H_{w\max }-\pi$.

The process of measurement surface optimization and parameter identification is shown in Figure 5. Specific implementation details are described in Section 5. The “triangle” in the upper left corner indicates the start, and the “circle” in the lower right corner indicates the end.

5. Simulation

This section validates the proposed kinematic calibration method for the Delta robot with 3-R(RPaR) configuration through computational cases. Initially, nominal geometric parameters are specified: r₁ = 0.127 m, r₂ = 0.046 m, l_p = 0.22751 m, l_n = 0.55361 m, and $l_{bi}^0$ = 0.1 m. Eliminating Δl_bi, the other 51-dimensional (or positional) error terms are set to 1 × 10⁻⁴ m, and angular (or orientation) error terms are set to 0.01° (1.745 × 10⁻⁴ rad). For instance, Δx_ai = 1 × 10⁻⁴ m and Δθ_i = 0.01°.

1. Verification of parameter identifiability. A total of 50 points are randomly selected in the Delta robot’s workspace, as detailed in Section 4.1 with n = 50. Figure 6 illustrates the singular values of the LMCM M from Equation (51), plotted in descending order on a log10 scale. The singular values of M form three distinct tiers, labeled I, II, and III, with magnitudes of approximately 10⁰, 10⁻⁷, and 10⁻¹⁵, respectively. Given that position and angle errors are on the order of 10⁻⁴, the magnitude of $‖v‖$ is also 10⁻⁴. According to inequality (55), parameters can be ignored when λ_i < 10⁻⁴. Therefore, parameters corresponding to tiers other than tier I will be ignored. The lowest point in tier I corresponds to the singular value λ₂₄.

Rewriting Equation (52) yields the following:

$${\displaystyle \sum _{i=1}^{24}{\lambda }_i{w}_i^{T}v}=-U^{T}b,$$

(62)

where v denotes the 51-component error vector excluding Δl_b1, Δl_b2, and Δl_b3. Figure 7 displays a scatter plot of the components of the column vectors w_i, with each color corresponding to the same i value. The jth component of w_i also serves as the coefficient for the jth error component of v. The horizontal axis in Figure 7 represents each of the 17 error components in each chain, ordered as chains 1, 2, and 3. Notably, the component values in w_i have been preprocessed based on the error scale; specifically, if the absolute value of a component in w_i is less than 10⁻⁴, it is set to zero.

From Figure 7, it is evident that the coefficients corresponding to components 8, 9, 11, 12, 13, and 14 in the error vectors within each chain consistently evaluate to zero. Consequently, these parameters are deemed unidentifiable. Furthermore, the coefficients of components 1, 2, and 3 within each chain exhibit sign symmetry in relation to the coefficients of components 15, 16, and 17, respectively. This indicates that either the error components 1, 2, and 3, or the error components 15, 16, and 17, cannot be individually discerned. What can be discerned are v_ci(1)–v_ci(15), v_ci(2)–v_ci(16), and v_ci(3)–v_ci(17), where v_ci(j) represents the jth component of the v_ci vector. In subsequent optimization, v_ci(1), v_ci(2) and v_ci(3) will be utilized to denote v_ci(1)–v_ci(15), v_ci(2)–v_ci(16), and v_ci(3)–v_ci(17), respectively. The remaining identifiable parameters for each chain are as follows:

$$v_{ci}={\left[\Delta x_{ai},\Delta y_{ai},\Delta z_{ai},\Delta {\phi }_i,\Delta {\gamma }_i,\Delta {\theta }_i,\Delta l_{pi},\Delta l_{ni}\right]}^T.$$

(63)

2. Determining the optimal measurement surface. The cylindrical workspace of the Delta robot is characterized by a maximum radius R_wmax = 0.265 m, a given minimum radius R_wmin = 0.100 m, and a height H_wm = 0.2 m, with its zero point in {W₀} along the Z axis positioned at Z_home = 0.3795 m. When evaluating $f\left(v_{coff}\right)$ in Problem (61), the surface coefficients in Equations (59) and (60) are considered given, and a point acquisition method on the surface needs to be specified. Here, $s\in \left[0,2\pi \right]$, $t\in \left[0,H_{wm}\right]$, and s and t are uniformly divided into 36 and 12 parts, respectively, generating a grid of points. Since the s direction forms a closed curve, it is iterated from i = 0 to 35 and j = 0 to 12, yielding a total of 36 × 13 non-repeating grid points in space. While more points theoretically enhance the accuracy of the grid, an excessive number would escalate measurement and computation times. In this study, the choice of points is balanced; even for trigonometric components with m, and n = 3 in Equations (59) and (60), a minimum of four sampling points per cycle is ensured.

This paper employs a hybrid solving algorithm that combines the Hook–Jeeves Direct Search Algorithm (DSA) with the Inner Point Method (IPM) for the optimization problem outlined in Problem (61). Initially, the constrained optimization problem is converted into an unconstrained one through the application of the interior point method. Here, the constraints refer to “<” and “>” constraints, while equality constraints can be ensured through the definition of the measurement surface. The augmented objective function needed for the interior point method is formulated as follows:

$$f_a\left(v_{coff},\tau \right)=f\left(v_{coff}\right)+\tau f_b\left(v_{coff}\right),$$

(64)

where τ > 0 is the penalty factor.

The essence of the interior point method is to transform a constrained optimization problem into a sequence of unconstrained optimization problems $\min f_a\left(v_{coff},{\tau }_k\right)$, where ${\tau }_{k+1}=\rho {\tau }_k$, as shown in “IPM” part of Figure 8. In the simulation process, ρ = 0.5, τ₀ = 1. The barrier function $f_b\left(v_{coff}\right)$ is defined as follows:

$$f_b\left(v_{coff}\right)=\frac{1}{R_{w\max }-R_q}+\frac{1}{R_q-R_{w\min }},$$

(65)

where R_q is achieved by separately solving the maximum value of $Q_r\left(s\right)$ and $Q_h\left(t\right)$, i.e., $R_q=\max \left(Q_r\left(s\right)\right)\cdot \max \left(Q_h\left(t\right)\right)$. As k increases, ${\tau }_k{f}_b\left(v_{coff}\right)$ will gradually decrease. The termination criterion for the iterative process is ${\tau }_k{f}_b\left(v_{coff}\right)<{\epsilon }_{IM}$. In simulation process, ${\epsilon }_{IM}=0.01$.

Each $f_a\left(v_{coff},{\tau }_k\right)$ undergoes optimization using the Hook–Jeeves DSA to determine the minimum value. While DSAs may not converge as rapidly as derivative-based methods like gradient descent or Newton’s method, the latter are susceptible to local optima. The DSA, being derivative-free, offers global optimality and requires fewer parameters than heuristic algorithms, making it well-suited for engineering applications. The Hook–Jeeves algorithm comprises two key components: Descent Direction Search (DDS) and Step Acceleration (SA), depicted in the “DSA” section of Figure 8. In DDS, the search progresses along the positive and negative directions of each unit component e_j of the v_coff vector sequentially, with a step size of d. The initial value of v_coff is denoted as $v_{coff}^{0,0}$. Upon identifying a direction of function reduction, the SA phase is initiated. During SA, the step length is accelerated using the acceleration factor $\alpha \left(v_{coff}^i-v_{coff}^{i,0}\right)$ to update v_coff and commence a new cycle of DDS. In the simulation process, α= 4. If a direction of function reduction is not detected, d = d/2, and DDS continues. The termination criterion for the Hook–Jeeves DSA is d < ε_DS. In simulations, the initial value of d is set to 0.005, with ε_DS = 10⁻⁶.

Figure 9 illustrates the iterative convergence process of the augmented objective function. The horizontal axis represents the number of iterations of the Hook–Jeeves DSA, which is the ordinal of the sequence minf_a(v_coff, τ_k). The convergence process is very rapid, with the augmented objective function entering a bottom plateau phase after approximately five iterations. The value of the augmented objective function decreases from 1239.6853 to 236.8953, a decrease of more than five times. Figure 10a and 10b, respectively, depict the iterative convergence process of the surface coefficients a_s0, a_s1, a_s3, b_s1, b_s2, b_s3 and a_t0, a_t1, a_t2, a_t3, b_t1, b_t2, b_t3 with the same horizontal axis as Figure 9. From Figure 10, it can be observed that the amplitudes of a_s0 and b_s0 are significantly larger compared to other trigonometric coefficients, and their initial values change noticeably during the initial iterations.

Figure 11 depicts the optimized surface profile. In the 3D representation shown in Figure 11a, the surface demonstrates pronounced periodicity along both the s and t directions. Moreover, the surface’s position and height in the Z direction align with the specified values. As illustrated in Figure 11b, the optimized surface profile adheres to the “<” and “>” constraints outlined in Problem (60), with both constraint boundaries being attained.

Table 2. The optimized result for measurement surface.

Title 1	Initial Values	Optimized Values	Conditions
as0/at0	0.2/2	0.53845/1.6278	D = 0.005, εDS = 10−6, ρ = 0.5, τ0 = 1, εIM = 0.01
as1/at1	0.01/0.2	−0.031114/0.025262
as2/at2	0/0	−0.0066394/0.054048
as3/at3	0/0	0.0078314/−0.025262
bs1/bt1	0/0	0.026959/0.078273
bs2/bt2	0/0	0.031114/−0.028414
bs3/bt3	0/0	0.020146/−0.015076
fa	1239.6853	236.8953

presents the conditions used in the optimization process, the values of the surface coefficients before and after optimization, and the values of the augmented objective function before and after optimization.

3. Identify the parameter errors. Firstly, simulation error data are generated for the measurement surface. Assuming a uniform error distribution, Δl_bi, Δx_ci, Δy_ci and Δz_ci are set to zero, and the remaining 42 positional and angular errors are set to 1 × 10⁻⁴ m and 0.01°, respectively. These errors are then added to the nominal parameter values accordingly. Utilizing the actual parameter values, the pose vectors for the points Q(s_i, t_j) on the measurement surface are computed directly by Equations (30), (34), and (35), without introducing linearization errors. This method ensures a more precise assessment of the comprehensive error magnitude. Additionally, besides acquiring the pose components of the moving platform control point T, the non-planar deformation angles for the three parallelogram mechanisms are also calculated. Notably, in the linearized Equation (50), $\Delta {\eta }_1^{i,j}$, $\Delta {\eta }_2^{i,j}$ and $\Delta {\eta }_3^{i,j}$ have been eliminated as intermediate variables.

Subsequently, parameter identification is performed based on the simulated error data. By employing the linearized Equation (51) and the identifiable error vector in Equation (63), a set of linear over-constrained equations is established. The identifiable error components in Equation (63) are then resolved using the least squares method, and the outcomes are presented in

Table 3. The calibration results of Delta robot’s geometric parameter errors.

	GV (mm or °)	CV i = 1 (mm or °)	Error i = 1 (%)	CV i = 2 (mm or °)	Error i = 2 (%)	CV i = 3 (mm or °)	Error i = 3 (%)
Δxai	0.1	0.099873	−0.13	0. 099746	−0.25	0.100294	0.29
Δyai	0.1	0.100206	0.20	0.099781	−0.22	0.099983	−0.02
Δzai	0.1	0.100100	0.10	0.100066	0.07	0.100170	0.17
Δφi	0.01	0.009941	−0.59	0.009910	−0.90	0.009924	−0.76
Δγi	0.01	0.009999	−0.01	0.010009	0.09	0.010005	0.05
Δθi	0.01	0.010012	0.12	0.009994	−0.06	0.010022	0.22
Δlpi	0.1	0.100067	0.07	0.100059	0.06	0.100026	0.03
Δlni	0.1	0.099801	−0.20	0.099937	−0.06	0.099735	−0.27

. In Table 3, “GV” denotes the error values provided during the construction of the simulation data, “CV” represents the error components computed based on Equations (51) and (63), and “Error” indicates the relative deviation between the two. Notably, Table 3 illustrates that the maximum absolute relative deviation is 0.9%, thereby affirming the efficacy of the proposed method for parameter error identification. This also corroborates the parameter identifiability judgment. Furthermore, it suggests that the error components not listed in Table 3 have a negligible impact on the overall error.

In Ref. [12], Table 3 presented the simulation results of kinematic parameter identification for a Delta robot. It showed that the accuracy of the identification results varied with different spatial distributions of measurement points, with Set 1 achieving the best accuracy. The maximum deviation between the identified and given errors in Ref. [12] is 1.6%, which is larger than 0.9%. In our research, the measurement surface or the distribution of measurement points is optimized, allowing the algorithm to determine a unique set of measurement points.

To verify the parameter identification performance of the algorithm under error conditions, random errors of ±5% ΔT_i,jare added to the pose error vector ΔT_i,j, and the error components are recalculated.

Table 4. Parameter error calibration results with random measurement noises.

	GV (mm or °)	CV i = 1 (mm or °)	Error i = 1 (%)	CV i = 2 (mm or °)	Error i = 2 (%)	CV i = 3 (mm or °)	Error i = 3 (%)
Δxai	0.1	0.102633	2.63	0.095674	−4.33	0.096670	−3.33
Δyai	0.1	0.099947	−0.05	0.101908	1.91	0.097117	−2.88
Δzai	0.1	0.093943	−6.06	0.093356	−6.64	0.094659	−5.34
Δφi	0.01	0.010051	0.51	0.009332	−6.68	0.010415	4.15
Δγi	0.01	0.010105	1.05	0.009738	−2.62	0.010226	2.26
Δθi	0.01	0.008964	−10.36	0.009779	−2.21	0.009993	−0.07
Δlpi	0.1	0.100067	1.27	0.100059	0.059	0.100026	0.026
Δlni	0.1	0.102633	2.63	0.095674	−4.33	0.096670	−3.33

presents a typical result of parameter identification under random measurement errors. The maximum identification error in the identification results is −10.36%. Although this is considerably larger than the results in Table 3, it is still sufficient for engineering applications.

6. Experiments

The dynamic calibration method is validated by experiments. The nominal geometric parameters are specified: r₁ = 0.065 m, r₂ = 0.035 m, l_p = 0.1042 m, l_n = 0.1650 m. The motor used has a rated current of 3.47 A, a torque constant of 25.9 mNm/A, a no-load speed of 8810 rpm at 24 V voltage, and a reduction ratio of 12. In contrast to the top-mounted configuration illustrated in Figure 1, in the experiment, the Delta robot is laterally installed, and a preloaded spring in the upper kinematic chain compensates for the torque of the integrated joint. The spring torque is defined as follows:

$${\tau }_s=k_s{\theta }_1+{\tau }_{s0},$$

(66)

where k_s and τ_s0 are the proportional coefficient and initial torque of the spring, respectively. During the grinding task, the robot usually follows a smooth and slow curve with very small acceleration and deceleration. Therefore, the parameters in the simplified dynamic model in this paper mainly include friction and mass, and the inertia parameters closely related to acceleration are not reflected, i.e., $I\ddot{\Theta }$, J^TF_A, and Γ_CM are negligible. Consequently, these three factors will be disregarded in the identification process. Additionally, accounting for the frictional forces of integrated joints, the motor torque Γ will be adjusted to incorporate the effect of the frictional torque Γ_c. The Coulomb model is employed as the friction model,

$${\tau }_{c1}\left({\dot{\theta }}_i\right)={\tau }_{ci0}sign({\dot{\theta }}_i),$$

(67)

where Γ_c = [τ_c1, τ_c2, τ_c3]^T. Therefore, the dynamic equation can be re-defined as follows:

$$\Gamma =-J^T{F}_G-{\Gamma }_c-{[{\tau }_s,0,0]}^T$$

(68)

The excitation trajectory adopts a uniform synchronous curve in joint space, i.e., θ₁(t) = θ₂(t) = θ₃(t) = vt + θ_min, where v is the moving speed, with a motion range of θ_min = −32° and θ_max = 72°. A suitable velocity, such as 5°/s, is selected. Under the servo driver’s CSP mode, the system moves from the starting point to the ending point at a constant speed. The joint angle is obtained through the encoder at the rear of the motor, and the torque is acquired by reading the motor’s current value and combining it with the torque constant and reduction ratio. Both the encoder values and current values are read in real-time at a 1 ms cycle from the corresponding variables inside the driver using the TwinCAT software (v3.1.4024.11). The same procedure is repeated from the ending point back to the starting point. The joint angle is obtained through the encoder at the rear of the motor, and the torque is acquired by reading the motor’s current value and combining it with the torque constant and reduction ratio. Both the encoder values and current values are read in real-time at a 1 ms cycle from the corresponding variables inside the driver using the TwinCAT software.

This simple curve form can fully meet the identification needs of the dynamic model in this paper. Additionally, this curve can be implemented at the early stage of robot controller development, without requiring interpolation in Cartesian space or worrying whether the T point after synthesizing the three joint motions always remains within the robot workspace.

Figure 12 display the raw torque-angle data and the data after applying moving average filtering for positive and negative movements. Based on the filtered data, nonlinear least squares fitting was conducted to obtain the calibration results of the dynamic parameters, as shown in

Table 5. Dynamic parameter calibration results.

Index	Parameter	Value	Unit	Meaning
1	me	0.9735	kg	the mass of moving platform
2	τc10	78.9	mNm	the friction of joint 1
3	τc20	35.3	mNm	the friction of joint 2
4	τc30	45.1	mNm	the friction of joint 3
5	ks	−152.9	mNm/rad	the proportional coefficient of the spring
6	τs0	−313.0	mNm	the intial torque of the spring

. Figure 12 also depicts the model prediction curves based on the identified dynamic parameters. The torque RMSEs between the filtered data and model-predicted data for the three joints are 14.5529 mNm, 26.2232 mNm, and 31.7745 mNm, which are approximately 1.42%, 2.55%, and 3.09% of the joints’ rated torques, respectively. Apart from slightly larger deviations at the beginning and end, the prediction curves fit the actual curves very well in most regions.

In Ref. [28], Table 6 provided the simulation results of RMSE between the given torques and the predicted torques of the identified model for a Delta robot. Ref. [28] did not provide the rated torques of robot joints. To compare with the results in our research, it is assumed that the maximum output torque in Ref. [28] is 0.7 times the rated torque, making it consistent with the situation our research. Additionally, the units in Ref. [28] are converted to mNm. First, looking at the offline parameter identification results in Table 6 of Ref. [28], the RMSE was 632.7 mNm, which is 4.4% of the rated torque. The identification process in our paper is also offline, and the maximum RMSE between the model-predicted torque and the actual torque is 3.09% of the rated torque, which is better than the offline identification results in Ref. [28]. Additionally, Ref. [28] provides online identification results, with the corresponding torque RMSE being 31.5 mNm, or 0.22% of the rated torque. It is worth noting that the results in Table 6 of Ref. [28] are simulation results, while our research provides experimental results. The results in our research include data noise and other unknown factors.

7. Conclusions

This paper focused on the calibration of kinematic errors and dynamic parameters for Delta robots driven by integrated joints with 3-R(RPaR) configuration in machining tasks, including position errors and frictions of integrated joints.

The influence of parallelogram mechanism dimension errors was analyzed based on the vector loop method for the 3-R(RPaR) configuration. The modeling of the in-plane deviation angle of the connecting platform link was conducted. This paper also defined the spatial deformation angle of the parallelogram mechanism constrained by the robot’s DoF, ultimately forming a 54-parameter kinematic error model and its corresponding linearized model.

By employing the SVD of LMCM, this paper established criteria for the identifiability of error component combinations based on the inequality between the singular value and the error vector norm. A measurement surface with Fourier series form containing 14 coefficients was conducted. The condition number of LMCM derived from the measurement surface was defined as the optimization objective. A hybrid algorithm that combines the Hook–Jeeves DSA with IPM was successfully applied to solve the optimal measurement surface under constraints.

The friction coefficients of integrated joints and other dynamic parameters were obtained by fitting the joint torque-angle pairs measured along specific trajectories using nonlinear least squares regression.

The simulation and experiment results validated the proposed kinematic and dynamic calibration methods, providing a foundation for precise control of integrated joints and high-precision motion of robots.

Based on the actual conditions of the grinding process, this study assumed that the motion process is always smooth and uniform, and based on this, the simplified dynamic model was established. If the motion behavior of the robot during the grinding process differs significantly from the assumption, the simplified model proposed in this paper will no longer be applicable, and the calibration accuracy will also be affected.

Our next step is to implement force–position hybrid control during the grinding process, examine the quality of the actual machined surface, and further optimize the dynamic performance through online optimization.