New identification method for computer numerical control geometric errors

Abstract

To address the problems of the low accuracy of geometric error identification and incomplete identification results of the linear axis detection of computer numerical control (CNC) machine tools, a new 21-item geometric error identification method based on double ball-bar measurement was proposed. The model between the double ball-bar reading and the geometric error term in each plane was obtained according to the three-plane arc trajectory measurement. The mathematical model of geometric error components of CNC machine tools is established, and the error fitting coefficients are solved through the beetle antennae search particle swarm optimization (BAS–PSO) algorithm, in which 21 geometric errors, including roll angle errors, were identified. Experiments were performed to compare the optimization effect of the BAS–PSO and PSO and BAS and genetic particle swarm optimization (GA–PSO) algorithms. Experimental results show that the PSO algorithm is trapped in the local optimum, and the BAS–PSO is superior to the other three algorithms in terms of convergence speed and stability, has higher identification accuracy, has better optimization performance, and is suitable for identifying the geometric error coefficient of CNC machine tools. The accuracy and validity of the identification results are verified by the comparison with the results of the individual geometric errors detected through laser interferometer experiments. The identification accuracy of the double ball-bar is below 2.7 µm. The proposed identification method is inexpensive, has a short processing time, is easy to operate, and possesses a reference value for the identification and compensation of the linear axes of machine tools.

Introduction

The CNC machine-on-machine measurement system is composed of the machine tool body and trigger probe. The probe generates a trigger signal, which prompts the machine tool’s 3D coordinate reading device to obtain a reading and realize 3D coordinate measurement of the measured point. The machining accuracy of the tested workpiece is guaranteed by the hardware installation accuracy and machining parameters of the machine tool. If the machine tool is used to build an on-machine measurement system to measure the machining accuracy of the workpiece, then the positioning and measurement accuracy of the machine tool structure and reading system as the measurement device and reference must be three times higher than the machining accuracy of the machine tool. Thus, only the machine tool hardware is used to ensure that the method cannot satisfy this requirement. The error separation method must be utilized to separate the geometric errors of the machine tool, and a high-precision error model ¹ needs to be established to compensate for the measurement error.

At present, the geometric error identification methods for computer numerical control (CNC) machine tools mainly use dual-frequency laser interferometers, double ball-bar, laser trackers,^2–4 and other equipment. ⁵ Laser interferometers use the principle of interference to measure displacement errors. Through various specific combinations of measuring mirrors, single measurements of positioning errors, straightness, verticality, and pitch and yaw errors of each axis can be derived. Although the laser interferometer has high detection accuracy, its installation and debugging process is complicated, time-consuming, and unable to detect the roll angle error of each axis, and one measurement can only measure the geometric error on one spatial measurement line, and the measurement efficiency is low ⁶ . The laser interferometer can only measure in a limited space, and there is no method to achieve the full spatial coverage and rapid error measurement. These shortcomings can be addressed by adopting the double ball-bar method. The double ball-bar can separate and identify 21 geometric errors, including rolling angle errors, by measuring conventional three-plane arc trajectories. In recent years, domestic and foreign scholars have proposed many methods for identifying spatial errors on the basis of the double ball-bar. Tsutsu et al. ⁶ used a three-axis linkage method to set two different measurement trajectories of the double ball-bar to identify the position-independent geometric errors of the rotation axis with a tilt angle. Chen et al. ⁷ constructed a robust identification model to identify six geometric errors of the tilting rotary table on the basis of the double ball-bar measurement results at three different installation positions. Zhong et al. ⁸ designed a nine-minute circular measurement path for the comprehensive error of each axis of a five-axis machine tool to identify spherical deviations with a double ball-bar. However, the geometric error elements were not identified. Guangsheng et al. ⁹ identified seven geometric errors of the turntable through the four special installation methods of the double ball-bar on the turntable. Gang et al. ¹⁰ utilized the six-step measurement method of the double ball-bar and the least squares method to eliminate the influence of the installation error of the former. They identified 10 geometric errors of the multi-axis CNC machine tool turntable, including verticality and position errors. Compared with the geometric error of the rotating table identified by the double ball-bar, fewer research focus on the linear axis error identification based on the double ball-bar. Wenjie et al. ¹¹ used the ridge regression method to identify the polynomial coefficients of the geometric error elements on the basis of a linear mapping model between the geometric error source and the length error of the double ball-bar. However, this method took the cost of losing some information and reducing the accuracy, and the identification accuracy was not it is high, and the laser interferometer detection results were compared and verified for only one direction of the X-axis. Kai et al. ¹² constructed an identification model for homogeneous coordinate transformation and used least squares method to solve the fitting coefficients to identify 18 position-dependent geometric errors of the linear axis. However, these results were not compared with the measurement results of the high-precision detection instruments for verification, and unable to verify the identification accuracy problem.

The geometric errors of CNC machine tools have complex changing laws and many influencing factors. Therefore, addressing the geometric error identification problem of CNC machine tools on the basis of the double ball-bar can be classified as a multidimensional model optimization problem. Genetic algorithm, particle swarm optimization (PSO), and beetle antennae search algorithm (BAS) are the primary solutions for complex multidimensional models at home and abroad. Wan et al. ¹³ utilized the particle swarm algorithm with strong global search capability to solve the fitting coefficients and the geometric errors identification of CNC machine tools, but the algorithm fell into local extreme values when addressing complex optimization problems up to 36 dimensions, resulting in inaccurate optimal solutions. BAS is a new and efficient optimization algorithm in accordance to the foraging behavior of long beetles; it was proposed by Jiang and Li ¹⁴ in 2017. Compared with other algorithms, BAS is simple and flexible, avoids local optimal solutions, and is suitable for high-dimensional search. However, its search consumes a large amount of time, and its global search capability needs further improvement. Combining the characteristics of the two algorithms, this study proposes a beetle antennae search–particle swarm optimization (BAS-PSO) algorithm, which is applied to solve the identification coefficients of the geometric errors of CNC machine tools on the basis of the double ball-bar to obtain enhanced global search capabilities and algorithm stability, and then compared the results of the three-axis single geometric error detection through laser interferometer for a comparative verification.

Geometric error model of a three-axis vertical machining center

When establishing the position-dependent geometric error model of machine tools, the coordinates of the workpiece and tool should be transformed into a fixed coordinate system (i.e. a reference coordinate system) because the positions of both change with the movement of the motion pair. ¹⁵ In general, the machine tool base or the lathe bed is set as the reference coordinate system, and the motion coordinate system is established on each moving part relative to the reference coordinate system (i.e. X-, Y-, and Z-slide rails). The entire machine tool’s motion circuit consists of workpiece kinematic chain and tool kinematic chain according to the cutting tool tips (Figure 1).

The workpiece kinematic chain and tool kinematic chains of the three-axis vertical machine tool are respectively sequenced as follows:

1. Bed R → Y-axis → X-axis → workpiece W

2. Bed R → Z-axis → spindle S → tool T

Figure 1.

Schematic diagram of machine tool kinematic chain.

In this study, the machine tool is regarded as a rigid body. Hence, no deformation can occur during processing. Based on the rigid body assumption, the machine tools can be regarded as a rigid chain link composed of revolute joints and prismatic joints. The transfer relationship between adjacent parts can be expressed as a 4 × 4 homogeneous transformation matrix (HTM), and the entire motion model of the machine tool can be derived finally through transforming the point coordinates from adjacent local coordinate systems to reference coordinate systems.

Two kinds of motion models are established on the basis of whether or not geometric errors are considered, namely, the ideal kinematic model ¹⁶ (IKM) and the real kinematic model (RKM). Both kinematic chains have been established for three-axis vertical CNC machine tools (workpiece kinematic chain and tool kinematic chain). The workpiece kinematic chain is taken as an illustration in the subsequent discussion.

The HTM of the Y-axis in IKM can be represented by

$$\mathrm{Trans}(\mathrm{y})=\left[\begin{array}{cccc}1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & y_m\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right],$$

(1)

Ignoring higher-order terms, the HTM of the Y-axis geometric error can be represented by

$$\mathbf{E}(\mathrm{y})=\left[\begin{array}{cccc}1\hfill & -{\epsilon }_z(y)\hfill & {\epsilon }_y(y)\hfill & {\delta }_x(y)\hfill \\ \multicolumn{1}{c}{{\epsilon }_z(y)}& 1\hfill & -{\epsilon }_x(y)\hfill & {\delta }_y(y)\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y(y)}& {\epsilon }_x(y)\hfill & 1\hfill & {\delta }_z(y)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]$$

(2)

According to equations (1) and (2), the HTM of the Y coordinate system in bed R, including each geometric error, can be obtained as

$$\begin{array}{c}\multicolumn{1}{c}{{}_{\mathrm{Y}}{}^{R}T=\mathrm{Trans}(\mathrm{y})\times \mathbf{E}(\mathrm{y})}\\ \multicolumn{1}{c}{=\left[\begin{array}{cccc}1\hfill & -{\epsilon }_z(y)\hfill & {\epsilon }_y(y)\hfill & {\delta }_x(y)\hfill \\ \multicolumn{1}{c}{{\epsilon }_z(y)}& 1\hfill & -{\epsilon }_x(y)\hfill & {\delta }_y(y)+y_m\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y(y)}& {\epsilon }_x(y)\hfill & 1\hfill & {\delta }_z(y)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]}\end{array}$$

(3)

Similarly, the HTM of X in IKM can be represented by

$$\mathrm{Trans}(\mathrm{x})=\left[\begin{array}{cccc}1\hfill & 0\hfill & 0\hfill & x_m\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right],$$

(4)

Ignoring higher-order terms, the HTM of the X-axis geometric error can be represented by

$$\mathbf{E}(\mathrm{x})=\left[\begin{array}{cccc}1\hfill & -{\epsilon }_z(x)\hfill & {\epsilon }_y(x)\hfill & {\delta }_x(x)\hfill \\ \multicolumn{1}{c}{{\epsilon }_z(x)}& 1\hfill & -{\epsilon }_x(x)\hfill & {\delta }_y(x)\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y(x)}& {\epsilon }_x(x)\hfill & 1\hfill & {\delta }_z(x)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right],$$

(5)

Using the assembly method, with the Y-axis as the mounting reference, a verticality error between the X and Y axes can be obtained.

$$\begin{array}{c}\mathrm{\Delta }x=x_m-x_m\times \cos \left({\gamma }_{x,y}\right)=x_m\times (1-\cos \left({\gamma }_{x,y}\right))\approx 0\hfill \\ \multicolumn{1}{c}{\mathrm{\Delta }y=-x_m\times \sin \left({\gamma }_{x,y}\right)\approx -x_m·{\gamma }_{x,y}}\end{array},$$

(6)

Then, the error transfer matrix can be expressed as

$$\mathbf{E}(\mathrm{x})=\left[\begin{array}{cccc}1\hfill & -{\epsilon }_z(x)\hfill & {\epsilon }_y(x)\hfill & {\delta }_x(x)\hfill \\ \multicolumn{1}{c}{{\epsilon }_z(x)}& 1\hfill & -{\epsilon }_x(x)\hfill & {\delta }_y(x)-x_m·{\gamma }_{x,y}\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y(x)}& {\epsilon }_x(x)\hfill & 1\hfill & {\delta }_z(x)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]$$

(7)

According to equations (4) and (7), the HTM from the X-axis coordinate system to the Y-axis coordinate system in RKM is expressed as

$$\begin{array}{c}\multicolumn{1}{c}{{}_X{}^{Y}T=\mathrm{Trans}(\mathrm{x})\times \mathbf{E}(\mathrm{x})}\\ \multicolumn{1}{c}{=\left[\begin{array}{cccc}1\hfill & -{\epsilon }_z(x)\hfill & {\epsilon }_y(x)\hfill & x_m+{\delta }_x(x)\hfill \\ \multicolumn{1}{c}{{\epsilon }_z(x)}& 1\hfill & -{\epsilon }_x(x)\hfill & {\delta }_y(x)-x_m·{\gamma }_{x,y}\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y(x)}& {\epsilon }_x(x)\hfill & 1\hfill & {\delta }_z(x)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right],}\end{array}$$

(8)

Setting the position coordinate vector of the theoretical processing point of the workpiece in the workpiece coordinate system as ${}_W{}^{X}T=[\begin{array}{cccc}{X}_w\hfill & Y_w\hfill & Z_w\hfill & 1\hfill \end{array}{]}^T$ , the HTM of workpiece W with the geometric error relative to machine bed R is expressed as

$$\begin{array}{c}\multicolumn{1}{c}{{}_W{}^{R}T={}_Y{}^{R}T·{}_X{}^{Y}T·{}_W{}^{X}T}\\ \multicolumn{1}{c}{=\left[\begin{array}{cccc}1\hfill & -{\epsilon }_z(y)\hfill & {\epsilon }_y(y)\hfill & {\delta }_x(y)\hfill \\ \multicolumn{1}{c}{{\epsilon }_z(y)}& 1\hfill & -{\epsilon }_x(y)\hfill & {\delta }_y(y)+y_m\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y(y)}& {\epsilon }_x(y)\hfill & 1\hfill & {\delta }_z(y)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]·}\\ \multicolumn{1}{c}{\left[\begin{array}{cccc}1\hfill & -{\epsilon }_z(x)\hfill & {\epsilon }_y(x)\hfill & x_m+{\delta }_x(x)\hfill \\ \multicolumn{1}{c}{{\epsilon }_z(x)}& 1\hfill & -{\epsilon }_x(x)\hfill & {\delta }_y(x)-x_m·{\gamma }_{x,y}\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y(x)}& {\epsilon }_x(x)\hfill & 1\hfill & {\delta }_z(x)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]·\left[\begin{array}{c}{X}_w\hfill \\ \multicolumn{1}{c}{Y_w}\\ \multicolumn{1}{c}{Z_w}\\ \multicolumn{1}{c}{1}\end{array}\right]}\end{array}$$

(9)

Based on the similar principal, the RKM of the Z-axis ${}_Z{}^{R}T$ and the HTM between spindle S and the Z axis and between tool T and spindle S can be established. Similarly, the HTM of tool T relative to machine bed R, which is used to set up the entire kinematic model of the machine tool, can be expressed as

$${}_T{}^{R}T={}_Z{}^{R}T·{}_S{}^{Z}T·{}_T{}^{S}T·T_t,$$

(10)

where $T_t$ is the position coordinate vector of the tool tip in the tool coordinate system T and is set to $T_t=[\begin{array}{cccc}0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}]$ .

In IKM, the tool position should be consistent with the machining point of workpiece. However, considering the non-negligible geometric errors, these positions do not coincide in most cases, thus affecting the machining accuracy and surface quality.

The difference between equations (9) and (10) is the relative displacement error matrix between the tool tip and the workpiece, that is, the geometric error matrix of the machine tool.

$${}_T{}^{W}E={}_T{}^{R}T-{}_W{}^{R}T$$

(11)

Measurement principle and identification method

Double ball-bar’s measurement principle

Figure 2 shows a schematic of the measuring principle of the double ball-bar. In the diagram, the coordinates of the center of the precision ball adsorbed on the worktable mounting base are $O(0,0,0)$ , whereas that adsorbed on the main shaft mounting base rod are $P(x,y,z)$ . When the workbench performs a theoretical circular interpolation movement, the coordinates of the center of the precision ball absorbed on the main axis are $P\prime (x\prime ,y\prime ,z\prime )$ due to errors. Then, Δx, Δy, and Δz are the components of the error between points $P$ and $P\prime$ in the X, Y, and Z directions, respectively.

Figure 2.

Measuring principle diagram of the double ball-bar method.

The true length containing the geometric error of the double ball-bar can be represented by

$$(r+\mathrm{\Delta }r{)}^2=(x+\mathrm{\Delta }x{)}^2+(y+\mathrm{\Delta }y{)}^2+(z+\mathrm{\Delta }z{)}^2,$$

(12)

where Δx, Δy, and Δz are all within a micrometer range, and the quadratic term is ignored. The measuring error of the double ball-bar can be expressed as

$$\mathrm{\Delta }r=\frac{1}{r}(\mathrm{\Delta }x\times x+\mathrm{\Delta }y\times y+\mathrm{\Delta }z\times z),$$

(13)

where $\mathrm{\Delta }r$ is the change in the length of the double ball-bar. The above equation is the theoretical basis for measuring the geometric error of the CNC machine tool through the double ball-bar method, from which the error separation model of the CNC machine tool can be deduced.

Geometric error identification method

According to theoretical derivation, equation (10) can be written as

$$\begin{array}{c}\multicolumn{1}{c}{{}_T{}^{R}T=\left[\begin{array}{cccc}1\hfill & -{\epsilon }_z(z)\hfill & {\epsilon }_y(z)\hfill & {\delta }_x(z)-z_m·{\gamma }_{x,z}\hfill \\ \multicolumn{1}{c}{{\epsilon }_z(z)}& 1\hfill & -{\epsilon }_x(z)\hfill & {\delta }_y(z)-z_m·{\gamma }_{y,z}\hfill \\ \multicolumn{1}{c}{-{\epsilon }_y(z)}& {\epsilon }_x(z)\hfill & 1\hfill & z_m+{\delta }_z(z)\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]·}\\ \multicolumn{1}{c}{\left[\begin{array}{cccc}1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]·\left[\begin{array}{cccc}1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill & d\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \end{array}\right]·\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{1}\end{array}\right],}\end{array}$$

(14)

where $d$ is the length of the machine tool (i.e. length of the mounting base rod).

According to equations (9), (11), and (14), let ${}_T{}^{W}E=\left[\begin{array}{cccc}\mathrm{\Delta }x\hfill & \mathrm{\Delta }y\hfill & \mathrm{\Delta }z\hfill & 1\hfill \end{array}\right]$ in the solution process and discard the quadratic and high-order terms of the errors. Moreover, the simplified geometric error model of the CNC machine tool moving along the coordinate axis is obtained as

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{\Delta }x={\delta }_x(z)-z·{\gamma }_{x,z}+d·{\epsilon }_y(z)-X_w+Y_w({\epsilon }_z(x)+{\epsilon }_z(y))}\\ \multicolumn{1}{c}{-Z_w({\epsilon }_y(x)+{\epsilon }_y(y))-x-{\delta }_x(x)-{\delta }_x(y),}\end{array}$$

(15)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{\Delta }y={\delta }_y(z)-z·{\gamma }_{y,z}-d·{\epsilon }_x(z)-X_w({\epsilon }_z(y)+{\epsilon }_z(x))-Y_w+}\\ \multicolumn{1}{c}{Z_w({\epsilon }_x(x)+{\epsilon }_x(y))-x·{\epsilon }_z(y)-{\delta }_y(x)}\\ \multicolumn{1}{c}{+x·{\gamma }_{x,y}-y-{\delta }_y(y),}\end{array}$$

(16)

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{\Delta }z=z+{\delta }_z(z)+d+X_w({\epsilon }_y(y)+{\epsilon }_y(x))}\\ \multicolumn{1}{c}{-Y_w({\epsilon }_x(y)+{\epsilon }_x(x))}\\ \multicolumn{1}{c}{-Z_w+x·{\epsilon }_y(y)-{\delta }_z(x)-{\delta }_z(y)}\end{array},$$

(17)

where x, y, and z are the theoretical coordinate values of the measurement points.

Substituting the above equations into equation (13) yields

$$\mathrm{\Delta }r=\frac{1}{r}\left\{\begin{array}{c}x({\delta }_x(z)-z·{\gamma }_{x,z}+d·{\epsilon }_y(z)-X_w+Y_w({\epsilon }_z(x)+{\epsilon }_z(y))-Z_w({\epsilon }_y(x)+{\epsilon }_y(y))\hfill \\ \multicolumn{1}{c}{-x-{\delta }_x(x)-{\delta }_x(y))+y({\delta }_y(z)-z·{\gamma }_{y,z}-d·{\epsilon }_x(z)-X_w({\epsilon }_z(y)+{\epsilon }_z(x))}\\ \multicolumn{1}{c}{-Y_w+Z_w({\epsilon }_x(x)+{\epsilon }_x(y))-x·{\epsilon }_z(y)-{\delta }_y(x)+x·{\gamma }_{x,y}-y-{\delta }_y(y))+z(z+{\delta }_z(z)}\\ \multicolumn{1}{c}{+d+X_w({\epsilon }_y(y)+{\epsilon }_y(x))-Y_w({\epsilon }_x(y)+{\epsilon }_x(x))-Z_w+x·{\epsilon }_y(y)-{\delta }_z(x)-{\delta }_z(y))}\end{array}\right\},$$

(18)

Experiments indicate that the geometric error component is correlated with the movement position of the CNC machine tool, revealing that the geometric error value during the movement varies with the position and is not constant. This phenomenon should be considered using a simple method that involves fitting the geometric error components with a cubic polynomial function.

According to the quadratic analytical model of the geometric error of the machine tool, ¹⁷ the non-negligible angle deviation around the central axis perpendicular to the translation axis always exists, resulting in errors in straightness, and the yaw angle error is calculated based on the partial differential of the straightness error. Considering the Y-axis as an example, for the simplicity of calculation, the higher-order terms of errors are discarded, and the positioning, yaw angle, and X-direction straightness errors can be expressed as

$${\delta }_y(y)=\mathrm{dyy}\frac{y^2}{r^2},$$

(19)

$${\epsilon }_z(y)=2\mathrm{dxy}\frac{y}{r^2}$$

(20)

$${\delta }_x(y)=\mathrm{dxy}\frac{y^2}{r^2},$$

(21)

where y is the theoretical coordinate value of the measurement point, $r$ is the measuring radius, and $\mathrm{dyy},\mathrm{dxy}$ are the coefficients of the polynomial; these coefficients are the only requirements to calculate the geometric error of each component.

Apart from the three verticality errors, ¹⁸ position-dependent geometric component errors can be obtained to express the correlation equation into equation (18). An equation with 12 unknown parameters relevant to the length deviation of the double ball-bar, can also be obtained.

$$\mathrm{\Delta }r=\frac{1}{r}\left\{\begin{array}{c}x(\mathrm{dxz}{(\frac{\mathrm{z}}{r})}^2-z·{\gamma }_{x,z}+2d·\mathrm{dxz}(\frac{z}{r^2})-X_w+Y_w(2\mathrm{dyx}(\frac{x}{r^2})+2\mathrm{dxy}(\frac{y}{r^2}))\hfill \\ \multicolumn{1}{c}{-Z_w(2\mathrm{dzx}(\frac{x}{r^2})+\mathrm{eyy}(\frac{y}{r}))-x-\mathrm{dxx}(\frac{x}{r})-\mathrm{dxy}{(\frac{y}{r})}^2)}\\ \multicolumn{1}{c}{+y(\mathrm{dyz}{(\frac{z}{r})}^2-z·{\gamma }_{y,z}-2d·\mathrm{dyz}(\frac{z}{r^2})-X_w(2\mathrm{dxy}(\frac{y}{r^2})+2\mathrm{dyx}(\frac{x}{r^2}))}\\ \multicolumn{1}{c}{-Y_w+Z_w(\mathrm{exx}(\frac{x}{r})+2\mathrm{dzy}(\frac{y}{r^2}))-2x·\mathrm{dxy}(\frac{y}{r^2})-\mathrm{dyx}(\frac{x}{r}{)}^2+x·{\gamma }_{x,y}}\\ \multicolumn{1}{c}{-y-\mathrm{dyy}(\frac{y}{r}))+z(z+\mathrm{dzz}(\frac{z}{r})+d+X_w(\mathrm{eyy}(\frac{y}{r})+2\mathrm{dzx}(\frac{x}{r^2}))}\\ \multicolumn{1}{c}{-Y_w(2\mathrm{dzy}(\frac{y}{r^2})+\mathrm{exx}(\frac{x}{r}))-Z_w+x·\mathrm{eyy}(\frac{y}{r})-\mathrm{dzx}{(\frac{x}{r})}^2-\mathrm{dzy}{(\frac{y}{r})}^2)}\end{array}\right\},$$

(22)

where $\mathrm{dxx},\mathrm{dyy},\mathrm{dzz}$ are the fitting coefficients of positioning error when moving along the x, y, and z coordinate axis, $\mathrm{dyx},\mathrm{dzx},\mathrm{dxy},\mathrm{dzy},\mathrm{dxz},\mathrm{dyz}$ are the fitting coefficients of straightness error in other directions when x, y, and z coordinate axis is moving, $\mathrm{exx},\mathrm{eyy},\mathrm{ezz}$ are the fitting coefficients of rolling error when moving along the x, y, and z coordinate axis. ${\gamma }_{x,z},{\gamma }_{y,z},{\gamma }_{x,y}$ are the three verticality errors measured by the double ball-bar; $x,y,z,\mathrm{\Delta }r$ are the theoretical coordinate and radial error values of each measurement point; and $r$ is the measurement radius.

From the measurements, data with thousands of sets of radial error values $\mathrm{\Delta }r$ and corresponding angles $\theta$ can be obtained. By substituting these data into equation (22), a system of overdetermined equations with unknown parameters can be generated. This system is solved to calculate various geometric error values at different positions. The overdetermined equations do not have exact solutions, only approximate ones in most cases, and the least square solution is a commonly used method. This study used the PSO algorithm based on beetle antennae search algorithm (BAS) to solve the approximate solution of the overdetermined equations.

BAS–PSO algorithm

BAS–PSO

PSO is a random algorithm. When solving some complicated optimization problems, falling into local extreme values is easy, resulting in inaccurate optimal solutions. The BAS–PSO is a PSO algorithm based on BAS, which can accelerate the parallel calculation of global search and local search, and quickly obtain the global optimal solution by quickly approximating the fitness value of the objective function. The basic idea of this algorithm is to transplant the antenna search operation of BAS into the speed and position update process of PSO to enhance the global search capability. Using this algorithm to solve the approximate solution of the double ball-bar’s overdetermined equations can resolve the shortcomings of the PSO. The flow chart of the BAS–PSO is displayed in Figure 3.

Figure 3.

Flow diagram of the BAS–PSO algorithm.

The antenna search operation follows the following rules:

$$\begin{array}{c}\multicolumn{1}{c}{v{b}_i=-{\delta }^t·\overrightarrow{b}·\mathrm{sign}(f(x_{\mathrm{rt}})-f(x_{\mathrm{lt}}))}\\ \multicolumn{1}{c}{v_i^{j+1}=v_i^j+c_1·\mathrm{rand}·({\mathrm{Pb}}_i^j-x_i^j)+c_2·\mathrm{rand}·({\mathrm{Pg}}_i^j-x_i^j)}\\ \multicolumn{1}{c}{+c_3·\mathrm{rand}·v{b}_i}\\ \multicolumn{1}{c}{q_i^{j+1}=q_i^j+v_i^{j+1}}\end{array}$$

(23)

where j represent the number of iterations, $v_i^{j+1},q_i^{j+1}$ represent the ith speed and position vector after the jth iteration, and $v{b}_i$ represents the update rate generated by BAS-PSO, $c_3$ is the learning factor, sign() is a sign function.

This study refers to the related literature ¹³ and compares the search results of algorithms under multiple sets of inertia weight values. The optimization effect is better when $w_{min}=0.5,w_{max}=0.9$ . The inertia weight of the algorithm is adjusted from 0.5 to 0.9 during the dynamic change process, aimed at obtaining the optimal capability for global and local searching, respectively at the beginning and the later stage of the algorithm.

$$w=w_{min}+\frac{w_{max}-w_{min}}{m},w_{min}=0.5,w_{max}=0.9,$$

(24)

Algorithm test for identifying geometric errors

To test the performance of this algorithm, this study used the PSO, BAS, GA-PSO, ¹³ and BAS–PSO algorithms to calculate the double ball-bar readings of the XY plane machine tool clockwise and applied the fitting coefficients of the overdetermined equations determined by the objective function of equation (22). By calculating and comparing the sum of squares of the least squares errors of the four algorithms, the smaller the result, the closer the fitness value of the objective function is, and the result is the optimal solution of the overdetermined equation system. All algorithms are run on MATLABR2016B. The experimental platform uses AMD Athlon (tm) IIX4 641 quad-core processor, CPU at 2.79 GHz, and a PC with 12 GB of running memory.

PSO algorithm parameter settings: particle size PopSize = 36, maximum iterations maxgen = 30, maximum and minimum speed V_max = 0.5 and V_min =−0.5, acceleration factor c1 = c2 = 1.49445.

BAS algorithm parameter settings: maxgen = 30, step factor: eta = 0.95, initial step: step = 22.

GA–PSO algorithm parameter settings: PopSize=36, maxgen=30, c₁=c₂=1.49445, crossover probability CrossPercent=0.7, and the inertia weight selection method shown in equation (24).

BAS–PSO algorithm parameter settings: PopSize=36, maxgen=30, c₁=c₂=c₃=1.49445, step factor: eta=0.95, initial step: step=22, and the inertia weight selection method is the same as above.

PSO, BAS, GA-PSO, and BAS–PSO are repeatedly run for 10 times each. The single optimization results of the four algorithms are presented in Figure 4, and the test results are summarized in Table 1. Figure 4 shows the variation of the fitness function value with the changes of the number of iterations. In accordance to the measurement data of the XY plane, the fitness function value of the BAS–PSO and the GA-PSO algorithms both reache the minimum value in the ninth generation. Regardless of the number of calculations, the minimum value remains almost unchanged. The experimental results show that the BAS–PSO algorithm possesses fast convergence speed, high identification accuracy, satisfactory optimization performance, and suitability for the geometric error coefficients identification of CNC machine tools.

Figure 4.

Comparison of the optimization results of the four algorithms.

Table 1.

Test results of the four algorithms.

Algorithmname	Single optimizationoperation time/s	Number of iterations	Find the optimalvalue (mm)(two significant digits)	Whether to fallinto a localoptimum	Average of10 runs	Standard deviationof 10 runs
PSO	174.975	13	0.25	Yes	0.2357	0.0033
BAS	132.604	18	0.21	Yes	0.2096	0.0068
GA-PSO	127.064	9	0.20	No	0.1998	0.0028
BAS–PSO	107.082	9	0.20	No	0.1942	0.0014

From the comparison of the optimization results in Table 1, the following inferences are obtained. (1) The PSO algorithm and the BAS algorithm are trapped in the local optimal solution after multiple iterations, and the optimal values are 0.25 and 0.21 mm. Neither the GA-PSO nor the BAS-PSO algorithm is trapped in the local optimal solution, but the former can only obtain the suboptimal approximate solution of the problem, not the optimal one, and the BAS-PSO and the GA-PSO algorithm optimization results are similar. (2) The convergence speed of the BAS-PSO algorithm is faster than that of the BAS and PSO algorithms. The BAS-PSO and the GA-PSO algorithm both converge in the ninth generation, whereas the BAS and PSO algorithms converge in the 18th and 13th generations, respectively. (3) The single optimization operation time of the BAS–PSO algorithm is shorter than that of the BAS, PSO and the GA-PSO algorithms. Compared with the PSO, BAS and the GA-PSO algorithms, the BAS–PSO algorithm improves the speed of a single optimization operation by 63.4%, 23.8%, and 18.7%, respectively. This is also the biggest difference between the BAS-PSO and GA-PSO algorithms. (4) Based on the average value and standard deviation obtained by each algorithm after 10 independent runs, the accuracy and stability of the BAS–PSO algorithm are better than those of the PSO, BAS, and the GA-PSO algorithms.

Experiment and result comparison

The identification method was applied to the Shenyang Machine Tool VMC850E three-axis vertical machining center. Renishaw’s QW20 double ball-bar was used to measure and identify the geometric error of CNC machine tools. Several parameters are shown in Table 2.

Table 2.

Experiment parameters.

Item		Numerical value	Unit
Doubleball-bar	Nominal length	100	mm
	Resolution	0.1	μm
	Measuring range	−100 to +100	mm
	Measuring speed	480	mm/min
	Sampling rate	19.231	Hz
Testingmachinetool	X axis stroke	850	mm
	Y axis stroke	500	mm
	Z axis stroke	540	mm
	Cutting feed rate	1–20,000	mm/min
	Spindle speed	50–8000	r/min
	Driving type	Ball screw andserver motorfor driving
Experimentaltemperature		23.2 ± 2	°C

Three-plane test and error identification

A large stroke cannot be measured easily because the length of the cue bar is divided into 50, 100, and 150 mm and restricted by installation. The comprehensive geometric errors on the XY, YZ, and ZX planes of the machine tools were measured with the double ball-bar, and the length of the rod was set to 100 mm. On the XY plane, the movement stopped in the z direction. Similarly, the movement stopped in the X and Y directions on the YZ and ZX planes, respectively. A 360° (overrange 45°) arc measurement was performed on the XY plane, and a 220° (overrange 2°) arc measurement was performed on the ZX and YZ planes, as shown in Figure 5.

Figure 5.

Trace pattern of the double ball-bar test.

The QW20 double ball-bar was installed on the workbench, and the experiment was performed under a no-load condition of the machine tool. Figure 6 displays the field test chart of the double ball-bar. The spindle of the machine tool performed clockwise and counterclockwise circular path interpolation motions with a radius of 100 mm in each plane at a speed of 480 mm/min. The QW20 double ball-bar was rotated by 10° to select data once, and three measurements were performed in each plane experiment. The average value, which is taken as the final value, was detected by the double ball-bar to reduce the random errors. Take the clockwise operation of the XY plane machine tool as an example. Table 3 lists the double ball-bar measurement data at different positions. This method was used to solve the unknown coefficient of the geometric error model.

Figure 6.

Double ball-bar experiment on CNC machine tools.

Table 3.

Radial error values of the XY plane machine tools at different positions when running clockwise.

Erroroccurrenceposition (°)	Radial error $\mathrm{\Delta }\mathrm{r}$ (μm)	Erroroccurrenceposition (°)	Radial error $\mathrm{\Delta }\mathrm{r}$ (μm)	Erroroccurrenceposition (°)	Radial error $\mathrm{\Delta }\mathrm{r}$ (μm)	Erroroccurrence position (°)	Radial error $\mathrm{\Delta }\mathrm{r}$ (μm)
0	−32.0	90	−19.2	180	37.0	270	6.4
10	−27.2	100	−15.5	190	41.3	280	−0.6
20	−28.9	110	−7.9	200	41.8	290	−7.1
30	−29.0	120	−3.2	210	37.8	300	−12.4
40	−29.0	130	3.8	220	36.7	310	−18.0
50	−26.7	140	13.8	230	34.9	320	−19.2
60	−26.3	150	20.2	240	30.0	330	−23.0
70	−25.8	160	27.3	250	26.3	340	−27.0
80	−22.8	170	31.7	260	19.2	350	−30.3

According to the BAS-PSO algorithm flow chart shown in Figure 3, the corresponding parameters of each geometric error model are solved through programming (Table 4).

Table 4.

Parameters of the geometric error model.

Positioning error parameter	Straightness error parameter	Rolling error parameter	Verticality error parameter
$\mathrm{dxx}=-0.0049$	$\mathrm{dyx}=0.2341$	exx = 0.3824	${\gamma }_{x,y}=-74.5\mu m/m$
	$\mathrm{dzx}=0.2762$
$\mathrm{dyy}=0.0016$	$\mathrm{dxy}=0.3174$	eyy = 0.3271	${\gamma }_{y,z}=-131.7\mu m/m$
	$\mathrm{dzy}=0.2043$
$\mathrm{dzz}=0.0063$	$\mathrm{dxz}=0.7426$	ezz = 0.4136	${\gamma }_{x,z}=-65.7\mu m/m$
	$\mathrm{dyz}=0.8935$

Each parameter and related coordinate values were substituted into equations (19) to (21) to obtain the identification values of the machine tool geometric errors. Figures 7 to 9 present the identification results of the six position-dependent geometric errors of the three axes.

Figure 7.

Identification results of the X-axis position-dependent geometric errors.

Figure 8.

Identification results of the Y-axis position-dependent geometric errors.

Figure 9.

Identification results of the Z-axis position-dependent geometric errors.

Therein, three verticality error parameters were directly obtained by measurement of the double ball-bar.

According to the geometric error identification results of the machine tool shown in Figures 7 to 9, among all the geometric errors of CNC machine tools, the positioning error of each axis exerts the greatest impact on the overall geometric error and far exceeds the impacts of angle and straightness errors. The variation ranges of the positioning error of the X, Y, and Z axes are −20–0 μm, 0–6.4 μm, and 0–25 μm, respectively. The variation range of the Y axis positioning error is smaller than that of the two other moving axes. The absolute value of each geometric error increases with the increase in displacement because of the accumulation of errors and reaches the maximum at the displacement limit position.

Comparative analysis of results

To verify the accuracy of the identification method, the individual errors of the laser interferometer were determined within the corresponding measurement range in the three directions of the double ball-bar experiment. The measurement area was set to 200 mm × 200 mm × 200 mm, and data were collected once every 5 mm. Each measurement point was unidirectionally measured three times. The unit of the angular error was arc seconds, and the unit of the linear error was μm. Thus, for a comparison using the same result graph, the unit of the angular error was converted according to the technical information of the laser interferometer, as shown in equation (25). We measured the pitch and yaw Abbe arm in the X, Y, and Z directions, as shown in Table 5. After conversion, the straightness error caused by the corresponding angle error was obtained.

$$1\mathrm{arc}\sec =5\mathrm{\mu }\mathrm{m}/\mathrm{m}$$

(25)

Table 5.

Experimental measurement of Abbe arm.

Item	Numerical value (mm)
Pitch Abbe arm in Z direction	343
Yaw Abbe arm in Z direction	390
Pitch Abbe arm in Y direction	950
Yaw Abbe arm in Y direction	290
Pitch Abbe arm in X direction	620
Yaw Abbe arm in X direction	140

Figure 10 shows the on-site detection chart of the laser interferometer with the three verticality errors directly measured by the double ball-bar. The correctness of the three-axis roll angle error identification ( ${\epsilon }_x(x),{\epsilon }_{\mathrm{y}}(y),{\epsilon }_{\mathrm{z}}(z)$ ) of the linear axis cannot be fully proven by the conventional plane arc test method. An error comparison was not conducted. The comparison charts of various geometric errors at different positions on the three axes are illustrated in Figures 11 to 13. The three curves in the figure represent the error values calculated with the identification results of this method, and the three discrete points represent the average laser interferometer measurement values.

Figure 10.

Laser interferometer experiment on CNC machine tools.

Figure 11.

Comparison chart of the geometric errors at different positions on the X-axis.

Figure 12.

Comparison chart of the geometric errors at different positions on the Y-axis.

Figure 13.

Comparison chart of the geometric errors at different positions on the Z-axis.

The possible reasons for the stepped curve along the axis length and a certain interval in the measurement error of the positioning error in Figure 11 are the long travel of the machine tools’X-axis. Some sections of the screw are not parallel to the guide rail, leading to the occurrence of the above curve error. The positioning error in Figure 12 exhibits an upward wavy shape. The ball screw can be twisted, and the guide rail is extremely tight during a long processing period. If the guide rail or the screw nut pair is tight, or the ball screw is insufficiently strong, the ball screw may twist, and the above curve error may occur. The laser interferometer measures the geometric error of each linear axis ball screw, and the measurement accuracy is high. In Figure 13, the Z position error measured by the laser interferometer is not continuous due to the defects in the machine tool; the error has a large jump and fall. However, the double ball-bar identification method cannot predict the defects of the machine tool, which is a smooth curve. Therefore, the results obtained by the laser interferometer in Figures 11 to 13 fluctuate more severely than those of the double ball-bar.

The change trend of the geometric error value identified by the double ball-bar and the corresponding geometric error value detected by the laser interferometer are the same. Except for the positioning errors of machine tool X and Z axis, the identification error variation range of straightness error of machine tool X, Z direction, positioning error, and straightness error of Y axis is 0–2.7 μm, and the identification error variation range of the three-axis pitch and yaw error is 0–2.4 μm (as shown in Figures 11–13). The above-mentioned identification error is relatively small, which meets the accuracy requirements of machine tool error identification. Due to the problem of the X and Z axis ball screws of the machine tool being measured, the positioning error of the machine tool measured by the laser interferometer has discontinuities and mutation points. However, the identification model established in this study cannot identify the discontinuities and mutation points, resulting in relatively large identification errors, but the identified error is consistent with the actual measurement error change trend. It can identify the continuously changing error and is used for accurately identifying the error of the machine tool with normal function. Therefore, the geometric error identification model of the machine tool based on the double ball-bar established in this study has usability and universality.

To verify that the BAS-PSO algorithm is more accurate than the others mentioned in this paper, we used the four algorithms to process the same set of data and compared their results with the results obtained by the laser interferometer. During actual error identification of machine tools, GA-PSO and BAS-PSO algorithms are used to solve the geometric error fitting coefficients. The calculated fitting coefficient results do not differ considerably, and the result of the comparison with the graphical method is not notable. Hence, the geometric error comparison chart of the GA-PSO algorithm is no longer shown here. Figure 14 shows the comparison of the geometric errors of PSO and BAS. Figure 14(k) shows the error identification using the PSO algorithm, where k = 1, 3, 5. Figure 14(v) presents the error identification using BAS, where v = 2, 4, 6. As indicated in Table 6, the processing results of the four algorithms were compared with the results of the laser interferometer.

Figure 14.

Comparison of the corresponding geometric errors of two algorithms.

Table 6.

Comparison of the experimental results of four algorithms.

Algorithm	Maximumdeviation (μm)	Mean deviation(μm)
PSO	9.4	4.9
BAS	5.3	2.8
GA-PSO	3.0	1.7
BAS-PSO	2.7	1.5

The error range between the identification results and the laser interferometer results is obtained by comparing the maximum deviations of the four algorithms, and the stability of the algorithm identification is measured using the average deviation. The comparison of the experimental results in Table 1 showed that the BAS-PSO algorithm is more accurate than the three other algorithms. BAS-PSO has high identification accuracy and is suitable for identifying the geometric error coefficients of CNC machine tools.

Conclusions

1. In this study, based on homogeneous transformation matrix, the position-dependent geometric error model for a three-axis vertical machining center is established, and the cubic polynomial function is used to fit the geometric error components. In addition, the identification equations between the error coefficient to be identified and the radial error value measured by the double ball-bar are established.

2. An error coefficient identification algorithm based on BAS-PSO is proposed, which integrates the antenna search operation of the beetle antennae search algorithm into the speed and position update process of the PSO algorithm to enhance the global search capability. The BAS-PSO algorithm is applied to solve the geometric error fitting coefficients of the CNC machine tools. The comparison of the results of PSO, BAS and the GA-PSO shows that the proposed BAS-PSO algorithm is superior to the other three algorithms in terms of optimization performance.

3. Twenty-one geometric errors, including the rolling angle errors, are identified using the proposed identification method, which are then compared with the single geometric error results of the three axes detected through laser interferometer. The results indicate that the proposed method has high accuracy, a maximum identification error of 2.7 µm at a low cost, consumes less time, and is easy to operate.