Thermal error compensation on a computer numerical control machine tool considering thermal tilt angles and cutting tool length

Abstract

The present error compensation technology of computer numerical control machine tools ignores radial thermal tilt angle errors of the spindle, while the thermal-induced offset is closely related to the tilt angle and the handle length. To solve this problem, three models of spindle thermal errors are proposed for the thermal yaw, pitch angles and elongation, and error compensation is performed based on the thermal tilt angles and cutting tool length. A five-point method was applied to measure the spindle thermal drifts at different speeds by eddy current sensors, which could effectively analyse the changes in the position-pose of the errors. Fuzzy clustering and correlation analysis were applied to group and optimise the temperature variables and select the variables sensitive to thermal errors in order to depress the multicollinearity of the temperature variables and improve the stability of the model. Finally, the thermal offset compensation was conducted in three directions. The results indicate that back propagation has a better capability for nonlinear fitting, but its generalisation is far less than that of time series. While the structure of multiple linear regression analysis is simple, its prediction accuracy is not satisfied. Time series adequately reflects the dynamic behaviours of the thermal error, and the prediction accuracy can reach 94%, with excellent robustness under different cutting conditions. The thermal error compensation equation that includes thermal tilt angles and cutting tool length is suitable for actual conditions and can accurately describe the space-pose of the thermal deformation and improve the machining accuracy.

Introduction

Precision jig-boring machining has long been applied to manufacture complicated box-type parts. However, the accuracy may decrease and become far lower than the initial design value after the machine is used for a long period of time. This decreased accuracy over time primarily results from inadequate maintenance and accuracy stability and thermal error is the main factor leading to the inadequate accuracy, accounting for 70% of the total errors arising from various error sources. ¹ Thermal error accounts for an increasing proportion of the total error as machining tools become more sophisticated. Moreover, the dynamical characteristics of the spindle also affect the thermal error. Zhang et al. ² proposed a holospectrum-based balancing method to improve machining accuracy. Non-uniform temperature distribution causes thermal errors in computer numerical control (CNC) machine tools, and this distribution becomes nonlinear and non-stationary and varies with time. The mutual coupling of the location and intensity of the heat source, the thermal expansion coefficient and machine structure can create complex thermal characteristics. ³ Donmez et al. ⁴ proposed that changing temperatures produce thermal error, which is a major factor in the reduction in machining accuracy.

To study the thermal characteristics of a spindle, the temperature field and thermal errors must first be accurately measured. Existing measurement methods discussed in the literature mainly focus on the axial thermal elongation, while radial thermal tilt errors have not been sufficiently considered. In this article, five eddy current displacement sensors (five-point measurement) are applied to simultaneously measure axial and radial thermal deformations, and thermal error models are established based on the analysis of the space and pose of the thermal deformation.

There are two methods to establish the thermal error model, the theoretical model based on heat transfer and the experimental model. In recent years, the finite element method (FEM), which is used to analyse temperature fields and the thermal deformation of machine tools, has become a topic of particular interest. Min and Jiang ⁵ established a variety of thermal boundary conditions for a thermal error model based on the Fourier thermodynamic equation and analysed the temperature field distribution of the screw under different heat fluxes. Zhao et al. ⁶ proposed a method to calculate the thermal conductivity coefficient of the spindle surface and simulated and analysed the variation principles of the temperature field and the thermal deformation of the spindle. Creighton et al. ⁷ used the FEM to analyse the temperature distribution characteristics of a high-speed micro-milling spindle and constructed an exponential model of the axial thermal error related to the spindle speed and running time. However, the thermal error of a precision CNC machine tool is a mutual coupling of many complex factors that are affected by numerous variables; therefore, it is extremely difficult to establish a theoretical equation from the perspective of thermoelasticity and heat transfer.

When pure theoretic equations encounter difficulties, many scholars transfer their emphases to experimental research and construct thermal error models based on engineering experience. Neural networks can describe nonlinear mapping relationships very well. Since Rumelhart et al. ⁸ proposed learning methods of multilayer back propagation (BP), many scholars have begun to apply neural networks to thermal error modelling. Yang et al.^9,10 used artificial neural networks (ANNs) to establish a relationship between the temperature and thermal error of a spindle, and the model proved to be useful for making generalisations. Ouafi et al. ¹⁰ constructed an ANN model for the spindle thermal error, with the temperature drawn based on a statistical methodology, and carried out error compensation experiments that effectively improved the machining accuracy. Multiple linear regression analysis (MLRA) with a simple structure is a traditional modelling method for the thermal error. Jenq and Wei ¹¹ and Lee and Yang ¹² established a compensation model for the spindle thermal error using an MLRA method. Hong and Ibaraki ¹³ studied the thermal characteristics of a rotary axis on the five-axis machine, analysed the effect of thermal error on errors in the motion of the rotary axis and calculated how errors in motion can be caused by thermal errors together with geometric errors. Vyroubal ¹⁴ presented a cost-effective method focused on the error compensation of the machine’s thermal deformation in the axis direction of the spindle based on decomposition analysis. Vissiere et al. ¹⁵ measured the thermal drifts of the spindle with a new method in which measurement accuracy can reach even the nanometre scale. Time series (TS) analysis provides a set of approaches to process the dynamic data. The essence of this method is that all types of data are approximately described by mathematical models. Through analysis of these models, the internal structure of the data can be analysed so that we can forecast its trends and apply necessary controls to it. Wang and Yang, ¹⁶ Miao et al. ¹⁷ and Wang et al. ¹⁸ applied the TS analysis method to establish a spindle thermal error model and used it to compensate for thermal errors, thereby obtaining better results. In this article, eddy current displacement sensors are used to measure spindle thermal errors; such sensors have the advantage of contact-free measurement. In recent years, new technologies have been used for the manufacture of such sensors. ¹⁹ PT100 has a high accuracy, and Frohlich et al. ²⁰ quantitatively demonstrated its excellent stability, enabling it to be used for accurate temperature measurements.

With the thermal error model presented, the next work would be compensation for errors during machining. Fu et al. ²¹ and Miao et al. ²² built a model of the spindle axial thermal error using the multivariate linear regression method. Wang and Yang ¹⁶ also proposed a prediction model for the axial thermal deformation and applied the model to compensate a CNC machine. Ouafi et al. ¹⁰ presented a comprehensive modelling approach for real-time thermal error compensation based on multiple temperature measurements, and the thermal-induced errors of the spindle were reduced from 19 µm to less than 1 µm. There are other scholars who investigated the axial thermal error compensation method for the spindle, ²³ and the machining accuracy was effectively improved. Zhang et al. ²⁴ proposed a compensation implementation technique based on machine zero point shift and Ethernet data communication protocol for machine tools and thereby improved the machine accuracy.

Most of the present literature focuses on the axial expansion of the spindle, ignoring the radial tilt angle errors. However, thermal angle errors are key factors in compensating the terminal processing accuracy of the spindle. Thus, the spindle thermal error offset compensation should also consider the tilt angles and the length of the cutting tool. Therefore, in this article, three methods are proposed to modelling the spindle thermal errors on a jig-boring machine, considering the radial yaw and pitch angles and the length of the cutting tool, and the criteria for estimating the prediction goodness of the models are established. Thermal balance experiments were performed using a temperature displacement acquisition system to synchronously measure the temperature distribution and thermal deformation at different spindle speeds, and the effects of different spindle speeds on thermal characteristics were analysed. The fuzzy clustering analysis method was used to optimise the temperature variables, selecting the variables that were thermal error–sensitive. multi-input and multi-output (MIMO) ANN, MLRA and TS models were proposed for the axial elongation and radial thermal tilt angles of the spindle. Subsequently, a new set of sample data was used to validate the model robustness, and the predictive goodness of the three models was compared. Finally, the thermal offset compensation based on the cutting tool length and thermal tilt angles was conducted in three directions on the CNC machine, and the three models’ compensative efficiencies were analysed.

Experimental principles and equipment

Experimental system

The experimental system is shown in Figure 1, which focuses on the spindle of a precision CNC jig-boring machine. The changes in the temperature field and the thermal distortion of the spindle can be analysed by this system. The maximum speed of the spindle is N = 20000 r/min, and the cooling system of the spindle is controlled intelligently by the temperature. The front and rear bearings and motor are cooled.

Figure 1.

Spindle experiment setup: (a) boring machine and measurement system, (b) left view and (c) front view.

The measured equipment and functions are as follows: a synchronous acquisition system (developed by our group based on NI SCXI) is used to determine the temperatures and the thermal deformations. This system uses the precise magnetic temperature sensors of the Pt100 to measure the temperatures of the motor, bearings, pedestal and coolant of the spindle system and environment. High-precision eddy current sensors are utilised to measure the spindle thermal drifts. The system performs real-time synchronous acquisition of temperatures and thermal drifts.

Table 1 presents and describes the positions of the sensors on the machine. Temperatures are denoted as T ₁,…, T ₁₁, and the eddy current displacement sensors are defined as S ₁,…, S ₅.

Table 1.

Temperature sensor and displacement sensor mounting position.

Temperature sensors/displacement sensors	Installation location
T 1	Rear bearing
T	Spindle base
T 3	Rear bearing cooling out
T 4	Front bearing cooling out
T 5	Ambient temperature
T 6	Front bearing oriented X−
T 7	Front bearing oriented Y+
T 8	Motor oriented Y+
T 9	The cooling in
T 10	Motor cooling out
T 11	Motor oriented X−
S 1	Radial near X-axis
S 2	Radial near Y-axis
S 3	Radial distal X-axis
S 4	Radial distal Y-axis
S 5	Z axial direction

Measurement principle

The spindle thermal drifts are measured by a five-point method, ²⁵ with the displacement sensor fixture and measurement diagram shown in Figure 2. The spindle is parallel to the Z-axis, and the axial thermal expansion can be obtained by S ₅. The radial thermal yaw ${\theta }_x$ partial X-direction is measured by S ₁ and S ₃, and the radial thermal pitch ${\theta }_y$ partial Y-direction is measured by S ₂ and S ₄. The acquisition system recorded data each second.

Figure 2.

Spindle five-spot installation diagram: (a) the sensor fixture and (b) five-point measurement sketch map.

After the spindle has run for a long period, thermal elongation expands the axial direction, and thermal angle change causes a tilt in the radial direction due to the uneven gradient distribution of the temperatures, as shown in Figure 3, where the thermal yaw angle ${\theta }_x$ is defined as follows

$$\tan {\theta }_x=\frac{\Delta L}{D}$$

(1)

Figure 3.

Spindle thermal inclination sketch.

where i denotes the number of measurements. The thermal yaw angle is too small in this experiment (i.e. ${\theta }_x\to 0$ ), so

$${\theta }_x~\tan {\theta }_x$$

(2)

The thermal yaw can be obtained by applying equations (1) and (2)

$${\theta }_x=\frac{({L_3}^i-{L_1}^i)-({L_3}^0-{L_1}^0)}{D}$$

(3)

where ${L_3}^0$ and ${L_1}^0$ are the radial displacements between the sensor probe and the spindle, respectively, measured by S ₃ and S ₁ in the initial state, and ${L_3}^i$ and ${L_1}^i$ are the transient displacements during operation. D is the distance between S ₁ and S ₃ and between S ₂ and S ₄; D = 120 mm.

Similarly, the thermal pitch angle ${\theta }_y$ in the Y-direction can be obtained as follows

$${\theta }_y=\frac{({L_4}^i-{L_2}^i)-({L_4}^0-{L_2}^0)}{D}$$

(4)

Experimental results and analysis

The spindle speed affects the distribution of the temperature field and the magnitude of the thermal drifts. This section focuses on a set of experimental data. To simulate the actual changes in the spindle speeds during processing, the velocity was varied in the experiment, as shown in Figure 4.

Figure 4.

Distribution of the step speeds.

Time-domain analysis of thermal errors

The thermal drifts of the spindle are shown in Figure 5. The measurements indicate that the displacement trends increased gradually over time, eventually reaching thermal equilibrium. The axial elongation increased over time, and the direction was negative, which indicates that the elongation of the spindle expanded in the negative direction on the Z-axis. The time to equilibrium was approximately 385 min, with a maximum elongation of 39.6 µm. The thermal error was positive in the X-axis direction, which indicates that the spindle moved away from displacement sensors S ₁ and S ₃ during operation, deviating from the Z-axis and swinging to the negative direction of the X-axis on the XZ plane. The thermal yaw angle with the Z-axis is ${\theta }_x$ , and the maximum value of the heat offset error was 35.8 µm. The thermal error in the Y-direction was negative, which indicates that the spindle was closer to displacement sensors S ₂ and S ₄ during operation, deviating from the Z-axis. In the YZ plane, the spindle pitched to the negative direction of the Y-axis, and its thermal pitch angle with the Z-axis is ${\theta }_y$ , so that the maximum thermal offset was 20.2 µm.

Figure 5.

Thermal drifts of the spindle: (a) thermal errors and (b) the position and pose of the thermal spindle.

Temperature field distribution

The temperature variations in the spindle system are shown in Figure 6. The temperature trends of all measured points increased with time, but the increasing curves of the temperatures exhibited cyclical changes because the temperatures of the spindle were controlled by an intelligent cooling system, which set a threshold to reduce the temperature. If the temperature of the components reached this threshold, the cooling system activated. Therefore, the temperature exhibited fluctuation changes.

Figure 6.

Temperatures of the spindle.

As the front bearing, rear bearing and the motor were cooled separately, temperatures could be inhibited from unbounded increase. However, the heat removed by the cooling fluid was less than the heat generated by the spindle, so the overall trend was still an increase to thermal equilibrium. Approximately 320 min was required for the temperature to reach thermal equilibrium, at which point the rear bearing had the highest temperature of 30.4 °C due to its large capacity, heavy load and severe friction, which generated more heat. The motor temperature was 26.7 °C. The temperatures of the other measured points were approximately 26 °C.

Temperature fuzzy cluster grouping and optimisation

The famous scholar Ruspini ²⁶ was the first to propose the concept of fuzzy partition. Specifically, Ruspini introduced fuzzy set theory into cluster analysis. Other researchers then presented a variety of the fuzzy clustering analysis method based on fuzzy graph theory, including the biggest tree method based on fuzzy graph theory. ²⁷ This article groups the temperature variables of 11 measuring points using fuzzy clustering, applies statistical correlation to optimise the measuring points and calculates the correlation coefficient between each temperature variable and thermal error. Finally, the measurement points with the greatest correlation coefficients in each group are taken as the typical temperature variables.

Temperature fuzzy clustering

As the temperature variables have small values, system cluster analysis was used, and the variable packet flow is shown in Figure 7.

Figure 7.

Fuzzy clustering grouping.

Assuming the temperature variable T = {T ₁, T ₂,…, T_m } is the object to be analysed by fuzzy clustering, each object in T is T_k (k = 1, 2,…, m), whose characteristics can be described by a limited number of values. Therefore, there is a corresponding vector P(T_k ) = (T_{k 1}, T_{k 2},…, T_ks ) matched with the object T_k. T_kj (j = 1,2,…, s) is the jth characteristic value of T_k. P(T_k ) is the set of eigenvectors for T_k . The sample T is divided into c fuzzy subsets ${\tilde{T}}_1,{\tilde{T}}_2,\dots ,{\tilde{T}}_c$ by applying fuzzy clustering analysis, according to the similarity between the feature vectors.

Cluster analysis is also known as the hierarchical clustering method. The feature vectors are gradually clustered according to the distance criteria. The classification moves from greater to lesser distance until it reaches the desired classification. The general steps of system clustering are as follows:

1. Initialise the data. Assuming that the sample set T contains m subsets T ₁ ⁽⁰⁾, T ₂ ⁽⁰⁾,…, T_m ⁽⁰⁾, which form one class, calculate the distance between each pair of subsets and obtain a distance matrix $D^{(b)}$ with $m\times m$ dimensions.

2. Find the smallest element in the distance matrix $D^{(b)}$ (except diagonal elements). If the minimum element is the distance between T_i ^(b) and T_j ^(b), they will be merged into T_ij ^(b+1). Finally, obtain a new classification T ₁ ^(b+1), T ₂ ^(b+1),…, T_{m −1} ^(b+1).

3. Compute distances between the new categories to obtain the distance matrix $D^{(b+1)}$ after merging clusters.

4. Repeat the second step until the classification meets the requirements.

For m = 11, the number of packets is set to C = 4. After calculating, the optimal grouping is obtained using the clustering algorithm of the combination of Euclidean-centroid, as shown in Figure 8. The temperature variables are divided into Groups {T ₁}, {T ₅}, {T ₂, T ₃, T ₄, T ₆, T ₉, T ₁₀}, {T ₇, T ₈, T ₁₁}.

Figure 8.

Clustering dendrogram.

Optimisation of temperature variables

Based on the results of the above groups, correlation coefficients between the axial thermal error E and the temperature T_i can be calculated

$${\rho }_{T_{i}E}=\frac{{\displaystyle \sum _{j=1}^n}(T_{\mathrm{ij}}-{\overline{T}}_i)(E_j-{\overline{E}}_j)}{\sqrt{{\displaystyle \sum _{j=1}^n}{(T_{\mathrm{ij}}-{\overline{T}}_i)}^2}\sqrt{{\displaystyle \sum _{j=1}^n}{(E_j-{\overline{E}}_j)}^2}}$$

(5)

In the equation above, i = 1, 2,…, m is the ith measurement point of the temperatures, and j = 1, 2,…, n is the number of measurements. T_ij is the temperature of the ith measured point at the jth measurement. E_j is the thermal elongation at the jth measurement, ${\overline{T}}_i$ is the average temperature of the ith measured points, and ${\overline{E}}_j$ is the average value of the thermal elongation. Coefficients are shown in Table 2. The temperature variable with the highest coefficient is selected as the typical variable in each clustering group. T ₁₀ is the temperature of the outlet liquid of the motor coolant, and as the temperature has a more significant influence on the temperature of the motor, this value is designated as a key variable. Finally, T ₁, T ₅, T ₆, T ₇ and T ₁₀ are chosen as the typical temperature variables.

Table 2.

Correlation coefficients between temperature and axial thermal error.

Cluster group	Temperature sensors	Correlation coefficient	Key point
1	T 7	0.9902	T 7
	T 8	0.9739
	T 11	0.9737
2	T 2	0.8593	T 6 T 10
	T 3	0.9054
	T 4	0.9242
	T 6	0.9706
	T 9	0.8948
	T 10	0.9344
3	T 5	0.9546	T 5
4	T 1	0.9651	T 1

Thermal error modelling and prediction

Based on the five typical temperature variables determined by the packet optimisation, three spindle thermal error models are established, and the performance of the models is estimated through criteria provided to evaluate the goodness of the fitting.

Linear regression analysis

Linear regression analysis is a valid statistical method of determining an interdependent quantitative relationship between two or more variables. According to the principle of minimum error, a mathematical model describing the statistical relationship of the variables is determined.

MLRA

Assuming that there is a linear relationship between the dependent variable y and independent variable $x=(x_1,x_2,\dots ,x_m)$ ²⁸

$$y=b_0+b_1{x}_1+b_2{x}_2+\cdots +b_m{x}_m+\epsilon$$

(6)

where $b_0,b_1,\dots ,b_m$ are undetermined parameters and $\epsilon$ is the random error. Equation (6) is called the MLRA model.

For n sets of experimental data $(x_{i1},x_{i2},\dots ,x_{\mathrm{im}};y_i)$ , $i=1,2,\dots ,n$ , equation (6) can be expressed as

$$\{\begin{array}{c}{y}_1=b_0+b_1{x}_{11}+b_2{x}_{12}+\cdots +b_m{x}_{1m}+{\epsilon }_1\hfill \\ \multicolumn{1}{c}{y_2=b_0+b_1{x}_{21}+b_2{x}_{22}+\cdots +b_m{x}_{2m}+{\epsilon }_2}\\ \multicolumn{1}{c}{\cdots }\\ \multicolumn{1}{c}{y_n=b_0+b_1{x}_{n1}+b_2{x}_{n2}+\cdots +b_m{x}_{\mathrm{nm}}+{\epsilon }_n}\end{array}$$

(7)

Then, presume

$$Y=\{\begin{array}{c}{y}_1\hfill \\ \multicolumn{1}{c}{y_2}\\ \multicolumn{1}{c}{\cdots }\\ \multicolumn{1}{c}{y_n}\end{array},X=\left[\begin{array}{ccccc}1\hfill & x_{11}\hfill & x_{12}\hfill & \cdots \hfill & x_{1m}\hfill \\ \multicolumn{1}{c}{1}& x_{21}\hfill & x_{22}\hfill & \cdots \hfill & x_{2m}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{1}& x_{n1}\hfill & x_{n2}\hfill & \cdots \hfill & x_{\mathrm{nm}}\hfill \end{array}\right],B=\{\begin{array}{c}{b}_0\hfill \\ \multicolumn{1}{c}{b_1}\\ \multicolumn{1}{c}{\cdots }\\ \multicolumn{1}{c}{b_m}\end{array},\epsilon =\{\begin{array}{c}{\epsilon }_1\hfill \\ \multicolumn{1}{c}{{\epsilon }_2}\\ \multicolumn{1}{c}{\cdots }\\ \multicolumn{1}{c}{{\epsilon }_n}\end{array}$$

(8)

Equation (7) becomes

$$Y=\mathrm{XB}+\epsilon$$

(9)

Equation (9) is a matrix form of the MLRA model.

Least-square estimation of parameters

Assuming that parameters $\hat{b}=({\hat{b}}_0,{\hat{b}}_1,\dots ,{\hat{b}}_m)$ are the estimated values of the MLRA model, and there is

$$\hat{y}={\hat{b}}_0+{\hat{b}}_1{x}_1+{\hat{b}}_2{x}_2+\cdots +{\hat{b}}_m{x}_m$$

(10)

$\hat{b}=({\hat{b}}_0,{\hat{b}}_1,\dots ,{\hat{b}}_m)$ are the least-square estimations of the model parameters, making the sum of residual squares $Q({\hat{b}}_0,{\hat{b}}_1,\dots ,{\hat{b}}_m)={\displaystyle \sum _{i=1}^n}{(y_i-{\hat{y}}_i)}^2$ minimum, and they can obtained through the extremum principle; thus, the partial differential is 0

$$\{\begin{array}{c}\frac{\partial Q}{\partial {\hat{b}}_0}=-2{\displaystyle \sum _{i=1}^n}[y_i-({\hat{b}}_0+{\hat{b}}_1{x}_{i1}+{\hat{b}}_2{x}_{i2}+\cdots +{\hat{b}}_m{x}_{\mathrm{im}})]=0\hfill \\ \multicolumn{1}{c}{\frac{\partial Q}{\partial {\hat{b}}_1}=-2{\displaystyle \sum _{i=1}^n}[y_i-({\hat{b}}_0+{\hat{b}}_1{x}_{i1}+{\hat{b}}_2{x}_{i2}+\cdots +{\hat{b}}_m{x}_{\mathrm{im}})]x_{i1}=0}\\ \multicolumn{1}{c}{\cdots }\\ \multicolumn{1}{c}{\frac{\partial Q}{\partial {\hat{b}}_m}=-2{\displaystyle \sum _{i=1}^n}[y_i-({\hat{b}}_0+{\hat{b}}_1{x}_{i1}+{\hat{b}}_2{x}_{i2}+\cdots +{\hat{b}}_m{x}_{\mathrm{im}})]x_{\mathrm{im}}=0}\end{array}$$

(11)

That is

$$\{\begin{array}{c}n{\hat{b}}_0+b_1{\displaystyle \sum _{i=1}^n}{x}_{i1}+\cdots +b_m{\displaystyle \sum _{i=1}^n}{x}_{\mathrm{im}}={\displaystyle \sum _{i=1}^n}{y}_i\hfill \\ \multicolumn{1}{c}{{\hat{b}}_0{\displaystyle \sum _{i=1}^n}{x}_{i1}+b_1{\displaystyle \sum _{i=1}^n}{x}_{i1}^2+\cdots +b_m{\displaystyle \sum _{i=1}^n}{x}_{i1}{x}_{\mathrm{im}}={\displaystyle \sum _{i=1}^n}{x}_{i1}{y}_i}\\ \multicolumn{1}{c}{\cdots }\\ \multicolumn{1}{c}{{\hat{b}}_0{\displaystyle \sum _{i=1}^n}{x}_{\mathrm{im}}+b_1{\displaystyle \sum _{i=1}^n}{x}_{\mathrm{im}}{x}_{i1}+\cdots +b_m{\displaystyle \sum _{i=1}^n}{x}_{\mathrm{im}}^2={\displaystyle \sum _{i=1}^n}{x}_{\mathrm{im}}{y}_i}\end{array}$$

(12)

Equation (12) can be written in the form of a matrix formula

$$X^T\mathrm{XB}=X^{T}Y$$

(13)

If the matrix $X^{T}X$ is invertible, equation (13) has a unique solution

$$B=\left(\begin{array}{c}{\hat{b}}_0\hfill \\ \multicolumn{1}{c}{{\hat{b}}_1}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\hat{b}}_m}\end{array}\right)={\left(X^{T}X\right)}^{-1}{X}^{T}Y$$

(14)

Vector B is the least-square estimation of the MLRA model.

The maximum coefficient R ² is applied to determine a combination of independent variables

$$R^2=1-\frac{{\displaystyle \sum _{i=1}^n}{(y_i-{\tilde{y}}_i)}^2}{{\displaystyle \sum _{i=1}^n}{(y_i-{\overline{y}}_i)}^2}$$

(15)

The identified temperature variables $x=[\Delta T_1,\Delta T_5,\Delta T_6,\Delta T_7,\Delta T_{10}]$ are independent inputs, $i=1,\dots ,n$ and n = 89; then, the comprehensive thermal error model is as follows, including thermal elongation E, thermal yaw angle error ${\theta }_x$ and thermal pitch angle error ${\theta }_y$

$$\left(\begin{array}{c}E\hfill \\ \multicolumn{1}{c}{{\theta }_x}\\ \multicolumn{1}{c}{{\theta }_y}\end{array}\right)=\left(\begin{array}{ccccc}0\hfill & 0\hfill & 19.533\hfill & -21.205\hfill & -5.138\hfill \\ \multicolumn{1}{c}{-0.289}& 0\hfill & -3.878\hfill & -1.887\hfill & -2.961\hfill \\ \multicolumn{1}{c}{0.555}& -1.159\hfill & 0\hfill & -1.456\hfill & -1.685\hfill \end{array}\right)\left(\begin{array}{c}\Delta T_1\hfill \\ \multicolumn{1}{c}{\Delta T_5}\\ \multicolumn{1}{c}{\Delta T_6}\\ \multicolumn{1}{c}{\Delta T_7}\\ \multicolumn{1}{c}{\Delta T_{10}}\end{array}\right)+\left(\begin{array}{c}2.343\hfill \\ \multicolumn{1}{c}{-0.990}\\ \multicolumn{1}{c}{-1.919}\end{array}\right)$$

(16)

MIMO neural network modelling

BP model

The relationship of the spindle thermal error to temperature is strongly nonlinear. An MIMO model of the thermal error is established based on a feedback neural network (BP)

$$P=\psi (T)$$

(17)

where $P=(E,{\theta }_x,{\theta }_y)$ is an output vector of the model, including axial thermal elongation E, thermal yaw angle error ${\theta }_x$ on XZ plane and thermal pitch angle error ${\theta }_y$ on YZ plane. $T=(\Delta T_1,\Delta T_5,\Delta T_6,\Delta T_7,\Delta T_{10})$ is the input vector, $\psi$ is the neuron weight matrix and $T\to P$ is the transfer mapping function.

Model structure

Setting a desired target, the efficiency of learning and the number of the iteration, the network is trained using the trainlm algorithm. Figure 9 illustrates the structure of the model: the hidden layer contains two layers made up of 5 × 1 neurons. The network constitutes five temperature variable inputs and three thermal drift variable outputs.

Figure 9.

Neural network structure.

TS analysis

The basic idea of the TS analysis is that a mathematical model is established through the analysis of the time sequence of samples based on limited data of the observation system, accurately reflecting the dynamic dependency of the system, which is applied to predict and monitor the future behaviour of the system.

Stationarity determination and standardisation of thermal error sequence

Given that the sequence $\{X_t,t=0,1,2,\dots \}$ is a discrete stochastic process, the stationarity of the TS is determined by the augmented Dickey–Fuller (ADF) test. Assuming that three TS of the thermal errors of the spindle system are $E_t,{\theta }_{x,t},{\theta }_{y,t}$ , and the results demonstrate that the thermal error sequences are non-stationary, then the sequences are normalised by Gaussian standardisation so that the sequences exhibit smoothness

$$\mu =\frac{1}{n}{\displaystyle \sum _{i=1}^n}{x}_i$$

(18)

$${\sigma }^2=\frac{1}{n-1}{\displaystyle \sum _{i=1}^n}{(x_i-\mu )}^2$$

(19)

where ${\mu }_i$ and ${\sigma }_i^2$ represent the mean and variance of axial elongation and radial thermal declination angle errors, respectively.

The original sequences of the spindle thermal errors undergo standardising processing

$$E_t~N({\mu }_1,{\sigma }_1^2)\Rightarrow E{\prime }_t=\frac{E_t-{\mu }_1}{{\sigma }_1}$$

(20)

$${\theta }_{\mathrm{xt}}~N({\mu }_2,{\sigma }_2^2)\Rightarrow \theta {\prime }_{\mathrm{xt}}=\frac{{\theta }_{\mathrm{xt}}-{\mu }_2}{{\sigma }_2}$$

(21)

$${\theta }_{\mathrm{yt}}~N({\mu }_3,{\sigma }_3^2)\Rightarrow \theta {\prime }_{\mathrm{yt}}=\frac{{\theta }_{\mathrm{yt}}-{\mu }_3}{{\sigma }_3}$$

(22)

where $E{\prime }_t,\theta {\prime }_{\mathrm{xt}}\mathrm{and}\theta {\prime }_{\mathrm{yt}}$ are the new TS of thermal errors, and the new sequences approximate standard normal distribution.

Model identification and parameter estimation

The autocorrelation function (ACF) and the partial autocorrelation function (PACF) are applied to calculate the tailing of the series, and the results indicate that the new TS of the standardised thermal errors are the autoregressive and moving average hybrid model autoregressive moving average (ARMA)(p, q), which is as follows ²⁹

$$X_t={\varphi }_1{X}_{t-1}+\cdots +{\varphi }_p{X}_{t-p}+{\epsilon }_t-{\theta }_1{\epsilon }_{t-1}-\cdots -{\theta }_q{\epsilon }_{t-q}$$

(23)

where ${\varphi }_i(i=1,2,\dots ,p)$ is the autoregressive parameter and ${\theta }_j(j=1,2,\dots ,q)$ is the moving average parameter. Defining the latter operator as B

$$B{X}_t=X_{t-1},B^p{X}_t=X_{t-p}$$

(24)

the ARMA model is transformed into

$$\varphi (B)X_t=\theta (B){\epsilon }_t$$

(25)

Box et al. ²⁹ suggested that if the ARMA(p, q) model contains p order autoregressive AR(p) and q order moving average MA (q), its ACF is a pattern mixed exponential and attenuation sine wave after p-q order delay. Correspondingly, the partial correlation function is not an exact exponential form but, rather, is controlled by a mixture of index and decaying sine wave. The covariance is

$$\begin{array}{c}\multicolumn{1}{c}{{\gamma }_{t+k,t}=\mathrm{Cov}(X_{t+k},X_t)}\\ \multicolumn{1}{c}{}& =E[(X_{t+k}-{\mu }_{t+k})(X_t-{\mu }_t)]={\gamma }_k,\hfill \\ \multicolumn{1}{c}{k=(0,\pm 1,\pm 2,\dots ,n-1)}& \hfill \end{array}$$

(26)

The ACF is

$${\rho }_k=\frac{{\gamma }_k}{{\gamma }_0}$$

(27)

According to statistical theory, the covariance function of TS with stationarity and 9 mean is estimated as follows

$${\hat{\gamma }}_k=\frac{1}{N}{\displaystyle \sum _{t=k+1}^N}{X}_t{X}_{t-k}$$

(28)

Thus, the ACF is estimated as follows

$${\hat{\rho }}_k=\frac{{\hat{\gamma }}_k}{{\hat{\gamma }}_0}=\frac{{\displaystyle \sum _{t=k+1}^N}{X}_t{X}_{t-k}}{{\displaystyle \sum _{t=1}^N}{X}_t^2}$$

(29)

Yule–Walker equations can obtain autoregression coefficients; set $k=1,2,\dots ,p$ to obtain linear equations

$$\{\begin{array}{c}{\rho }_1={\varphi }_1+{\varphi }_2{\rho }_1+\cdots +{\varphi }_p{\rho }_{p-1}\hfill \\ \multicolumn{1}{c}{{\rho }_2={\varphi }_1{\rho }_1+{\varphi }_2+\cdots +{\varphi }_p{\rho }_{p-2}}\\ \multicolumn{1}{c}{⋮⋮⋮}\\ \multicolumn{1}{c}{{\rho }_p={\varphi }_1{\rho }_{p-1}+{\varphi }_2{\rho }_{p-2}+\cdots +{\varphi }_p}\end{array}$$

(30)

Replacing the theoretical autocorrelation ${\rho }_k$ with the estimated autocorrelation ${\hat{\rho }}_k$ , the autoregression coefficients of the Yule–Walker estimation can be obtained. Remembering that

$$\varphi =\left(\begin{array}{c}{\varphi }_1\hfill \\ \multicolumn{1}{c}{{\varphi }_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\varphi }_p}\end{array}\right),{\rho }_p=\left(\begin{array}{c}{\rho }_1\hfill \\ \multicolumn{1}{c}{{\rho }_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\rho }_p}\end{array}\right),P_p=\left(\begin{array}{ccccc}1\hfill & {\rho }_1\hfill & {\rho }_2\hfill & \cdots \hfill & {\rho }_{p-1}\hfill \\ \multicolumn{1}{c}{{\rho }_1}& 1\hfill & {\rho }_1\hfill & \cdots \hfill & {\rho }_{p-2}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\rho }_{p-1}}& {\rho }_{p-2}\hfill & {\rho }_{p-3}\hfill & \cdots \hfill & 1\hfill \end{array}\right)$$

(31)

the parameter $\varphi$ in equation (30) can be written as follows

$$\varphi =P_p^{-1}{\rho }_p$$

(32)

Defining $w{\prime }_t=\varphi (B)X_t$ and dealing with the ARMA model as a moving average process

$$w{\prime }_t=\theta (B){\epsilon }_t$$

(33)

the covariance $\gamma {\prime }_j$ of $w{\prime }_t$ can be represented by covariance ${\gamma }_j$ of $X_t$ , and $j=0,1,\dots ,q$

$$\begin{array}{cc}\multicolumn{1}{c}{\gamma {\prime }_j}& ={\displaystyle \sum _{i=0}^p}{\varphi }_i^2{\gamma }_j\hfill \\ \multicolumn{1}{c}{}& +{\displaystyle \sum _{i=1}^p}({\varphi }_0{\varphi }_i+{\varphi }_1{\varphi }_{i+1}+\cdots +{\varphi }_{p-i}{\varphi }_p)({\gamma }_{j+i}+{\gamma }_{j-i})\hfill \end{array}$$

(34)

Using the convention that ${\varphi }_0=-1$ , the covariance function of the MA (q) process is

$$\gamma {\prime }_0=(1+{\theta }_1^2+\cdots +{\theta }_q^2){\sigma }_{\epsilon }^2$$

(35)

$$\gamma {\prime }_k=(-{\theta }_k+{\theta }_1{\theta }_{k+1}+\cdots +{\theta }_{q-k}{\theta }_q){\sigma }_{\epsilon }^2$$

(36)

Estimations of the parameters ${\sigma }_{\epsilon }^2,{\theta }_q,{\theta }_{q-1},\dots ,{\theta }_1$ can be calculated under q order moving average, using iteration

$${\sigma }_{\epsilon }^2=\frac{{\gamma \prime }_0}{1+{\theta }_1^2+\cdots +{\theta }_q^2}$$

(37)

$${\theta }_j=-(\frac{{\gamma \prime }_j}{{\sigma }_{\epsilon }^2}-{\theta }_1{\theta }_{j+1}-{\theta }_2{\theta }_{j+2}-\cdots -{\theta }_{q-j}{\theta }_q)$$

(38)

${\theta }_0=0$ , ${\theta }_q,{\theta }_{q-1},\dots ,{\theta }_1$ can be set to 0 at the beginning of the iterations.

Setting the order range $p,q\in [0,5]$ and calculating the ARMA(p, q) model, the optimal orders p and q are determined by the Akaike information criterion (AIC). ³⁰ The basic idea is to construct a criterion function AIC(p, q) that considers both the original data fitting and the number of unknown parameters in the model. When the parameters and fitting residual variance ${\hat{\sigma }}_{\epsilon }^2$ make the AIC value reach a minimum, orders p and q are determined

$$\mathrm{AIC}(p,q)=\ln {\hat{\sigma }}_{\epsilon }^2+\frac{2(p+q+1)}{n}$$

(39)

where ${\hat{\sigma }}_{\epsilon }^2$ is the variance estimation of fitting residual errors, n is the total number of samples and $p,q=1,2,\dots ,L$ is highest order of the model.

After calculation, the new TS $E{\prime }_t$ of the thermal elongation is ARMA(5, 3) model, the thermal yaw angle error $\theta {\prime }_{x,t}$ is ARMA(2, 5) and the thermal pitch angle error $\theta {\prime }_{y,t}$ is ARMA(2, 1), and the corresponding models are as follows

$$\{\begin{array}{c}{E\prime }_t=1.583{E\prime }_{t-1}-0.471{E\prime }_{t-2}+0.644{E\prime }_{t-3}-1.246{E\prime }_{t-4}+0.488{E\prime }_{t-5}+{\epsilon }_t^1+0.263{\epsilon }_{t-1}^1-0.267{\epsilon }_{t-2}^1-0.959{\epsilon }_{t-3}^1\hfill \\ \multicolumn{1}{c}{{\theta \prime }_{x,t}=0.18{\theta \prime }_{x,t-1}+0.779{\theta \prime }_{x,t-2}+{\epsilon }_t^2+1.139{\epsilon }_{t-1}^2+0.163{\epsilon }_{t-2}^2-0.139{\epsilon }_{t-3}^2+0.231{\epsilon }_{t-4}^2+0.376{\epsilon }_{t-5}^2}\\ \multicolumn{1}{c}{{\theta \prime }_{y,t}=1.996{\theta \prime }_{y,t-1}-0.997{\theta \prime }_{y,t-2}+{{\epsilon }^3}_{y,t}-{{\epsilon }^3}_{y,t-1}}\end{array}$$

(40)

Assuming that vector $\mu$ includes the means of the three TS of the spindle thermal errors, and the vector $\sigma$ includes the standard deviations

$$\mu =\left(\begin{array}{c}{\mu }_1\hfill \\ \multicolumn{1}{c}{{\mu }_2}\\ \multicolumn{1}{c}{{\mu }_3}\end{array}\right)=\left(\begin{array}{c}-24.681\hfill \\ \multicolumn{1}{c}{-9.845}\\ \multicolumn{1}{c}{-10.095}\end{array}\right),\sigma =\left(\begin{array}{c}{\sigma }_1\hfill \\ \multicolumn{1}{c}{{\sigma }_2}\\ \multicolumn{1}{c}{{\sigma }_3}\end{array}\right)=\left(\begin{array}{c}14.11\hfill \\ \multicolumn{1}{c}{6.272}\\ \multicolumn{1}{c}{6.571}\end{array}\right),$$

(41)

the new TS are reversed according to the following transformation, and the final spindle system thermal error model is

$$\left(\begin{array}{c}{E}_t\hfill \\ \multicolumn{1}{c}{{\theta }_{x,t}}\\ \multicolumn{1}{c}{{\theta }_{y,t}}\end{array}\right)=\left(\begin{array}{ccc}{\sigma }_1\hfill & \hfill & \hfill \\ \multicolumn{1}{c}{}& {\sigma }_2\hfill & \hfill \\ \multicolumn{1}{c}{}& \hfill & {\sigma }_3\hfill \end{array}\right)\left(\begin{array}{c}{E\prime }_t\hfill \\ \multicolumn{1}{c}{{\theta \prime }_{x,t}}\\ \multicolumn{1}{c}{{\theta \prime }_{y,t}}\end{array}\right)+\left(\begin{array}{c}{\mu }_1\hfill \\ \multicolumn{1}{c}{{\mu }_2}\\ \multicolumn{1}{c}{{\mu }_3}\end{array}\right)$$

(42)

Model predictions and comparison

The experimental data in Figure 4 are taken as the training sample, and the sample data number is 89. Then, the three models – MIMO neural network (BP), MLRA and TS – are used to predict the thermal drifts of the spindle. The predicted and the measured values are compared in Figures 10–12.

Figure 10.

Axial thermal elongation: (a) prediction and measurement and (b) residual error.

Figure 11.

Radial thermal yaw angle error: (a) prediction and measurement and (b) residual error.

Figure 12.

Radial thermal pitch angle error: (a) prediction and measurement and (b) residual error.

Next, we establish the evaluation criteria for the model fitting. Assuming the absolute value of the residual error is $\left|e_i\right|$ , we set its minimum as ${\left|e_i\right|}_{min}$ , its maximum as ${\left|e_i\right|}_{max}$ and its mean value as $\overline{\left|e_i\right|}$ . The root mean square error is RMSE, the determination coefficient is R ² and the predictive ability is $\eta$

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n}{(y_i-{\tilde{y}}_i)}^2}$$

(43)

$$R^2=1-\frac{{\displaystyle \sum _{i=1}^n}{(y_i-{\tilde{y}}_i)}^2}{{\displaystyle \sum _{i=1}^n}{(y_i-{\overline{y}}_i)}^2}$$

(44)

$$\eta =1-\frac{{\displaystyle \sum _{i=1}^n}|y_i-{\tilde{y}}_i|}{{\displaystyle \sum _{i=1}^n}\left|y_i\right|}$$

(45)

where $y_i$ is the measurement value, ${\tilde{y}}_i$ is the predicted value, ${\overline{y}}_i$ is the average value of the measurement and $i=1,\dots ,n$ , n are the data points. The fitting performance parameters of BP, MLRA and TS are shown in Tables 3–5.

Table 3.

Fitting goodness parameters for thermal elongation E.

Output	${\left|e_i\right|}_{min}(\mu \mathrm{m})$	${\left|e_i\right|}_{max}(\mu \mathrm{m})$	$\overline{\left|e_i\right|}(\mu \mathrm{m})$	RMSE	$R^2$	$\eta (\%)$
BP	0.01	3.42	0.581	0.688	0.998	97.7
MLRA	0.04	3.12	0.863	1.089	0.994	96.5
TS	0.00	1.50	0.335	0.473	0.999	98.7

RMSE: root mean square error; BP: back propagation; MLRA: multiple linear regression analysis; TS: time series.

Table 4.

Fitting goodness parameters for thermal yaw angle ${\theta }_x$ .

Output	${\left|e_i\right|}_{min}$ "	${\left|e_i\right|}_{max}$ "	$\overline{\left|e_i\right|}$ "	RMSE	$R^2$	$\eta (\%)$
BP	0.01	6.14	0.652	0.829	0.982	93.5
MLRA	0.01	3.10	1.090	1.301	0.956	89.1
TS	0.00	3.44	0.441	0.730	0.986	95.6

RMSE: root mean square error; BP: back propagation; MLRA: multiple linear regression analysis; TS: time series.

Table 5.

Fitting goodness parameters for thermal pitch angle ${\theta }_y$ .

Output	${\left|e_i\right|}_{min}$ "	${\left|e_i\right|}_{max}$ "	$\overline{\left|e_i\right|}$ "	RMSE	$R^2$	$\eta (\%)$
BP	0.00	2.93	0.374	0.485	0.995	96.3
MLRA	0.01	3.03	0.820	1.069	0.973	92.0
TS	0.00	1.62	0.242	0.416	0.996	97.6

RMSE: root mean square error; BP: back propagation; MLRA: multiple linear regression analysis; TS: time series.

The determination coefficient R ² can then be calculated with a boundary region of [0, 1], which characterises the extent to which the original data fall within the curve at a certain confidence level, with a larger R ² indicating that the curve fits the original data better, and that the argument for explaining the model is more powerful. The root mean squared error of the fitted and measured values is denoted by RMSE, which indicates how closely the curve fits the original measurement data; a smaller value (closer to 0) indicates that the original data are more concentrated around the fitted curve on both sides. Analysing the parameters of fitting goodness in the three tables shows that the prediction behaviour of the three models is very good, and the maximum residual error, minimum error and average absolute error are very small, which indicates that the accuracy of the three models is high. The RMSE is near 0, and the determination coefficients are also close to 1. The prediction accuracy can reach 90%, and it further exhibits a strong fitting ability. The TS model is better than the BP neural network, which is superior to the MLRA.

Model validation

To verify the universality and generalisation of the models, a new data sample is used to predict the spindle thermal errors. Figure 13 describes the spindle speeds map and experimental results.

Figure 13.

New air cutting condition: (a) random distribution of the speed and (b) thermal drifts.

Then, using the BP, MLRA and TS models to predict thermal errors in the three directions of the spindle, the predicted values are compared to the measurements, as shown in Figures 14–16.

Figure 14.

Axial thermal elongation: (a) prediction and measurement and (b) residual error.

Figure 15.

Radial thermal yaw angle error: (a) prediction and measurement and (b) residual error.

Figure 16.

Radial thermal pitch angle error: (a) prediction and measurement and (b) residual error.

The prediction-goodness parameters of the three models are presented in Tables 6–8. For E, ${\theta }_x$ and ${\theta }_y$ , the predictive abilities of the TS for the three errors are greater than 94%, which are much higher than for BP and MLRA. The BP predictive performance is better than that of MLRA only for the prediction of E, while it is far lower for the prediction of the radial thermal tilt angle errors.

Table 6.

Prediction-goodness parameters for thermal elongation E.

Output	${\left|e_i\right|}_{min}(\mu \mathrm{m})$	${\left|e_i\right|}_{max}(\mu \mathrm{m})$	$\overline{\left|e_i\right|}(\mu \mathrm{m})$	RMSE	$R^2$	$\eta (\%)$
BP	0.03	6.48	2.63	3.200	0.944	89.8
MLRA	0.01	7.95	4.13	4.671	0.880	84.0
TS	0.00	6.13	1.19	1.790	0.982	95.4

RMSE: root mean square error; BP: back propagation; MLRA: multiple linear regression analysis; TS: time series.

Table 7.

Prediction-goodness parameters for thermal yaw angle ${\theta }_x$ .

Output	${\left|e_i\right|}_{min}$ "	${\left|e_i\right|}_{max}$ "	$\overline{\left|e_i\right|}$ "	RMSE	$R^2$	$\eta (\%)$
BP	0.27	6.14	2.58	2.934	0.720	78.6
MLRA	0.02	6.30	1.67	2.171	0.847	86.0
TS	0.00	6.07	0.70	1.131	0.958	94.1

RMSE: root mean square error; BP: back propagation; MLRA: multiple linear regression analysis; TS: time series.

Table 8.

Prediction-goodness parameters for thermal pitch angle ${\theta }_y$ .

Output	${\left|e_i\right|}_{min}$ "	${\left|e_i\right|}_{max}$ "	$\overline{\left|e_i\right|}$ "	RMSE	$R^2$	$\eta (\%)$
BP	0.07	9.82	3.56	4.656	0.496	71.5
MLRA	0.05	3.73	1.38	1.621	0.939	88.9
TS	0.00	2.13	0.48	0.771	0.986	96.1

RMSE: root mean square error; BP: back propagation; MLRA: multiple linear regression analysis; TS: time series.

Comparison of the prediction-goodness parameters of the three models shows that when the spindle speed condition changes, the predictive abilities of all models for new data decrease. However, the forecasting accuracy of the TS model remains above 90%, indicating that TS can adequately reflect the dynamic characteristics of spindle thermal errors, and it provides excellent generalisation under different conditions. In contrast, BP has a powerful fitting capability and can perfectly describe the nonlinear relationship of temperatures with spindle thermal errors, but the model is strongly dependent on the training sample. Under new conditions, its prediction-goodness parameters drop sharply, even becoming far lower than the traditional MLRA model.

Model universality analysis

We now further explain the robustness of the TS model for application to different cutting conditions and its capability to adequately reflect the thermal error dynamic characteristics.

• 1. The temperature rate of the spindle system has a great influence on the thermal error, affecting the model prediction accuracy.

The essential reason for thermal errors is the existence of a temperature gradient in the spindle system. The uneven distribution of the temperature field of the spindle causes different thermal expansions of the components, which causes a relative deviation of the geometric positions among the parts, leading to thermal deformation. The influence factors of the spindle thermal characteristics are very complex, as the rotor and stator of the motor, bearings and other components can generate different amounts of heat per unit time in different operating conditions so that the temperature rate $d\Delta T/\mathrm{dt}$ is different, which results in different rates $d\Delta E/\mathrm{dt}$ of thermal deformation.

In the two groups of experiments with different speeds, even if the temperatures are the same, the corresponding thermal deformations are not entirely equivalent, as shown by a comparison between the working condition for training models (1) shown in Figure 4 and the working condition for the validation of the models (2) shown in Figure 13(b). Figures 17–19 show that E, ${\theta }_x$ and ${\theta }_y$ change with temperature increases of the rear and front bearings, with the prediction curves of the E for condition 2 defined as TS2, BP2 and MLRA2, respectively.

Figure 17.

Axial thermal elongation with temperature: (a) $\Delta T_1=3{}^{\circ }\mathrm{C}$ and (b) $\Delta T_7=2.6{}^{\circ }\mathrm{C}$ .

Figure 18.

Radial thermal yaw angle error with temperature: (a) $\Delta T_1=3{}^{\circ }\mathrm{C}$ and (b) $\Delta T_7=2.6{}^{\circ }\mathrm{C}$ .

Figure 19.

Radial thermal pitch angle error with temperature: (a) $\Delta T_1=3{}^{\circ }\mathrm{C}$ and (b) $\Delta T_7=2.6{}^{\circ }\mathrm{C}$ .

In condition 2, when t = 165 min, the temperature rises of the rear/front bearings are $\Delta T_1=3{}^{\circ }\mathrm{C},\Delta T_7=2.6{}^{\circ }\mathrm{C}$ , respectively, and the measured and the predicted results are presented in Table 9. For example, the thermal elongations under conditions 1 and 2 are $E_1=-29.3\mu \mathrm{m}$ and $E_2=-20.6\mu \mathrm{m}$ , respectively, which are clearly not equal. The models’ predicted values in condition 2 are $E_{\mathrm{BP}}=-25.3\mu \mathrm{m}$ , $E_{\mathrm{MLRA}}=-23.8\mu \mathrm{m}$ and $E_{\mathrm{TS}}=-21.3\mu \mathrm{m}$ , which are not equal to the measured E ₂ but closer to E ₁. The reason for this deviation is that the independent variables of BP and MLRA are temperatures, and the models do not consider the influence of the temperature rates $d\Delta T/\mathrm{dt}$ of the different working conditions on the thermal error models. The model parameters are identified by condition 1, so the prediction accuracies of BP and MLRA for condition 1 are better than for condition 2. The analyses of ${\theta }_x$ and ${\theta }_y$ reach similar conclusions.

Table 9.

Comparison between the measured and the predicted results at $\Delta T_1=3^{\circ }\mathrm{C},\Delta T_7=2.6^{\circ }\mathrm{C}$ .

	Measurement 1	Measurement 2	BP 2	MLRA 2	TS 2
E (µm)	−29.3	−20.6	−25.3	−23.8	−21.3
${\theta }_x$ "	−6.9	−11.7	−6.5	−7.4	−11.5
${\theta }_y$ "	−8.6	−10.4	−8.3	−7.5	−10.7

BP: back propagation; MLRA: multiple linear regression analysis; TS: time series.

The BP neural network has a strong nonlinear mapping ability, but its physical explanation is weak, its predictive ability is heavily dependent on the typical training samples, and the approximation and generalisation ability of the BP model are closely related to the typicality of the learning samples. When the sample changes, the predictive ability decreases. Fortunately, the TS model has no temperature variables, thus avoiding this difficulty.

• 2. The intrinsic properties of thermal errors on machines fitted by the TS model.

The thermal errors of the spindle are low-frequency signals that change slowly, show an increasing trend in processing due to the temperature rise and grow approximately linearly. Meanwhile, influenced by cooling fluid, heat conduction inside the shaft and air convection, there is a fluctuation with a certain periodicity in the series of the thermal error, which is obviously nonlinear. Moreover, the thermal error sequence is a random process that is influenced by outside factors, such as environmental temperature, radiation and vibration of the machine, so that the series of the thermal error inevitably contains random noise. In summary, the thermal error series usually includes three parts: a trend term, a periodic component and random elements, and Figures 10(a), 11(a) and 12(a) and Figures 14(a), 15(a) and 16(a) show these features clearly.

The $\mathrm{AR}(p)$ part of the proposed ARMA model can present the trend and periodic characteristics of the thermal error, while the $\mathrm{MA}(q)$ part can effectively characterise the random item of the thermal deformation. Therefore, the ARMA model $\mathrm{ARMA}(p,q)$ provides an appropriate description of the behaviours of the spindle thermal error. Because BP and MLRA cannot take so many features of the spindle thermal error into account, the prediction accuracy of the TS model is better.

• 3. The nonlinear relationship between the temperature rise and the thermal error is an important consideration in model selection.

The machine thermal error is nonlinear not only in the time domain but also in the relationship with temperature rise, $\Delta T$ , which is shown in Figures 17–19. Although BP has a strong nonlinear processing ability, the generalisation ability is weak and it cannot adapt to new samples. The MLRA model is apparently not very good at dealing with nonlinear problems in engineering practice, while the TS model has a strong nonlinear prediction ability.

• 4. The hysteresis effect of the thermal error has a great impact on the selection of the model.

Because an effect follows its cause, the change in thermal deformation of the spindle lags behind changes in the temperature field. In the case of thermal elongation, the experimental data are $(E_i:\Delta {T^i}_1,\dots ,\Delta {T_{11}}^i)$ at the ith measurement, where the thermal stress caused by temperature could lead to thermal deformation, and this process requires a period. Hence, the acquired temperature and thermal deformation are not equivalent at this moment. The thermal deformation E_i is the accumulation of the effect of the temperature rises $\Delta {T_1}^{i-1},\dots ,\Delta {T_{11}}^{i-1},\dots ,\Delta {T_1}^{i-k},\dots ,\Delta {T_{11}}^{i-k}$ within the first k moments, and the present temperature rises $\Delta {T^i}_1,\dots ,\Delta {T_{11}}^i$ will affect the next moment E_{i +1}.

This conclusion was proven in the experiment. While the spindle speeds decreased, as shown in Figures 4 and 13(b), the dissipated heat of the spindle was greater than the generated energy, so the temperatures declined slightly, while the thermal deformations continued to increase slowly, as shown in Figures 17–19. Both BP and MLRA only consider the present temperatures and thermal errors, ignoring the delayed effect between them, so the two models do not fully express the essential hysteresis characteristics of the thermal errors. Contrasting finely with BP and MLRA, the TS model reflects the dynamic dependency of the thermal errors series by the following equation

$$X_t={\varphi }_1{X}_{t-1}+\cdots +{\varphi }_p{X}_{t-p}+{\epsilon }_t-{\theta }_1{\epsilon }_{t-1}-\cdots -{\theta }_q{\epsilon }_{t-q}.$$

• 5. TS model can better characterise the cumulative effect of the thermal errors.

The acquired thermal error at each moment in the experiment is the accumulation of the thermal errors at the previous stage; namely, there is an intrinsic dependency between adjacent observations. $\mathrm{AR}(p)={\varphi }_1{X}_{t-1}+\cdots +{\varphi }_p{X}_{t-p}$ in model $\mathrm{ARMA}(p,q)$ describes the essential dynamic dependency of $E_t\mathrm{and}{E}_{t-i},{\theta }_{x,t}\mathrm{and}{\theta }_{x,t-i},{\theta }_{y,t}\mathrm{and}{\theta }_{y,t-i}$ , while the models of BP and MLRA have no similar ability.

• 6.  $\mathrm{ARMA}(p,q)$ can better describe non-stationary characteristics of the thermal error series.

TS in practical engineering are commonly non-stationary, which can be realised as a stochastic process. The thermal drifts and temperatures of the spindle in this article are no exception, and they are proven to be non-stationary series by the ADF test algorithm (ADF).

In this article, when the original sequences $E_t,{\theta }_{x,t},{\theta }_{y,t}$ of the spindle thermal drifts are translated into new TS $E{\prime }_t,\theta {\prime }_{\mathrm{xt}},\theta {\prime }_{\mathrm{yt}}$ through Gaussian standardisation, the ADF test algorithm is performed to identify the $E{\prime }_t,\theta {\prime }_{\mathrm{xt}},\theta {\prime }_{\mathrm{yt}}$ stationarity. The calculated results indicate that the characteristic roots of the new series are within the unit circle, demonstrating that the standardised series are stationary and do not need to be smoothed further. Therefore, $\mathrm{ARMA}(p,q)$ can better adapt to the non-stationary characteristic of the thermal error series.

In summary, the most important feature of MLRA is that the model has a simple structure with ordinary prediction accuracy, although the adaptive ability of the model is not strong and not suitable for nonlinear problems. The BP neural network has a strong nonlinear mapping ability, but it depends heavily on the training samples, and the approximation and generalisation abilities of the network model are closely related to the typicality of the learning samples. If the difference between new experimental samples and the sample for estimating the parameters of models is great, the predictive accuracy values of BP and MLRA are poor.

The TS analysis provides a set of approaches to process dynamic data. The basic idea of the TS analysis is that a mathematical model, which accurately reflects the system dynamic dependency, is established through the analysis of the time sequence samples based on a limited sample of the observation system. The method’s primary meaning is that all types of data are approximately described by mathematical models. Through the model analysis, the internal structure of the thermal error series can be mined so that we can forecast its trends and apply necessary controls to it. Hence, the TS model can adequately reflect the dynamic characteristics of the spindle thermal error, and it has excellent generalisation ability under different conditions.

Thermal error compensation

Cutting tool compensation offset

Figure 20 shows the thermal drift of the spindle in plane XOZ, with the spindle parallel to the Z-axis. During the experiments, the deviation mainly arose from spindle deformation due to the slight temperature rise of the detection probe. Therefore, it can be approximately considered that there were no changes in the shape of the probe and the tool. Supposing that the detecting probe was AP, where P was the end and the probe had length L, the displacement sensors S ₁ and S ₃ in the X-direction were placed at points B and H, respectively, with a distance between them of D, but the measurement point S ₃ was not at the end of the probe. AT represents the handle. After the spindle deformed, the location of the probe translated to $A\prime P\prime$ , so that the cutting tool position was $A\prime T\prime$ . The axial thermal elongation PK was E, and the cutting tool tip drifted from T to $T\prime$ . Eventually, the thermal drifts of the tool were $\Delta O_x$ and $\Delta O_z$ in plane XOZ.

Figure 20.

Thermal error offset map of the cutting tool.

The experimental results showed that the spindle expanded in the negative Z-direction after being heated and elongated from point T to M, so that the thermal offset of the cutting tool tip in the Z-direction was $\Delta O_z$ , and the direction was negative

$$\Delta O_z=(L_t-L+\frac{E+L_d}{\cos {\theta }_x})\cos {\theta }_x-(L_t-L+L_d)$$

(46)

where $L_t$ is the cutting tool length, $L$ is the length of detection probe, $L_d$ is the distance between the measured point H of S ₃ sensor and the end point P of the probe, E is the thermal elongation and ${\theta }_x$ is the thermal yaw angle of the spindle in the XOZ plane.

Because the thermal yaw angle was small (i.e. ${\theta }_x\to 0$ )

$$\cos {\theta }_x\to 1$$

(47)

$$\Delta O_z=E$$

(48)

This indicates that the thermal offset of the cutting tool in the Z-direction was not related to the length of the tool.

The thermal offset of the cutting tool in the X-direction was $\Delta O_x$ , and the direction was negative

$$\Delta O_x=(L_t-L+\frac{E+L_d}{\cos {\theta }_x})\sin {\theta }_x+\Delta L_3$$

(49)

The thermal yaw angle was small, so

$$\begin{array}{cc}\multicolumn{1}{c}{}& \cos {\theta }_x\to 1\hfill \\ \multicolumn{1}{c}{}& \sin {\theta }_x\to {\theta }_x\hfill \end{array}$$

(50)

Applying equation (60) to equation (59), the thermal offset of the tool in the X-direction can be obtained

$$\Delta O_x=(L_t-L+E+L_d){\theta }_x+\Delta L_3$$

(51)

where $\Delta L_3$ is a real-time thermal offset HN measured by sensor S ₃. This indicates that the cutting tool thermal offset in the X-direction was closely related to the length of the tool.

Similarly, the thermal offset of the cutting tool in the Y-direction $\Delta O_y$ can be obtained

$$\Delta O_y=(L_t-L+\frac{E+L_d}{\cos {\theta }_y})\sin {\theta }_y+\Delta L_4$$

For $\begin{array}{cc}\multicolumn{1}{c}{}& \cos {\theta }_y\to 1\hfill \\ \multicolumn{1}{c}{}& \sin {\theta }_y\to {\theta }_y\hfill \end{array}$

$$\Delta O_y=(L_t-L+E+L_d){\theta }_y+\Delta L_4$$

(52)

where $\Delta L_4$ is a real-time thermal offset HN measured by sensor S ₄. This indicates that the cutting tool thermal offset in the Y-direction was closely related to the length of the tool. ${\theta }_y$ is the thermal pitch angle of the spindle in the YOZ plane.

Because the thermal offsets of the tool in the X-, Y- and Z-directions were all negative, we can set the vector $\overrightarrow{T\prime T}$ of the tool thermal offset

$$\overrightarrow{T\prime T}\{-\Delta O_x,-\Delta O_y,-\Delta O_z\}$$

(53)

Assuming that the original coordinate value of point W on the workpiece is $\mathrm{OW}(P_x,P_y,P_z)$ , the compensation offset is $-\overrightarrow{T\prime T}$ , and the new machining coordinate value is $\mathrm{OW}\prime (P{\prime }_x,P{\prime }_y,P{\prime }_z)$

$$\mathrm{OW}\prime =\mathrm{OW}-\overrightarrow{T\prime T}$$

(54)

That is

$$\{\begin{array}{c}{P\prime }_x=P_x+(L_t-L+E+L_d){\theta }_x+\Delta L_3\hfill \\ \multicolumn{1}{c}{{P\prime }_y=P_y+(L_t-L+E+L_d){\theta }_y+\Delta L_4}\\ \multicolumn{1}{c}{{P\prime }_z=P_z+E}\end{array}$$

(55)

The spindle axial thermal elongation E, radial thermal yaw angle ${\theta }_x$ in the X-direction and thermal pitch angle ${\theta }_y$ in the Y-direction can be predicted by the multiple linear regression model (MLRA), BP neural network model and TS model (TS). Equations (19), (20) and (49) can be substituted into equation (65), through which we obtain the spindle thermal error ultimate compensation equations. This is the basis of error compensation.

Real-time compensation control

A sketched map of the spindle thermal error compensation is provided in Figure 21, and the CNC system is the Siemens 840D. The compensation system is composed of three parts: the processing module, HMI interface and PLC compensation module. The processing module acquires the temperatures and displacements from the sensors installed on the machining tool and can establish the thermal error model. While the real-time compensation is carried out, this module can calculate the thermal offsets and send the ultimate compensation values to the HMI based on the 840D secondary developments. Subsequently, the offsets are accepted by the PLC, and the offset value is inserted into the processing instructions to implement compensation for the thermal errors.

Figure 21.

Thermal error compensation control.

Compensation result

The spindle thermal error compensation is implemented in three directions, and Figures 22–24 show the comparison of the measured value of the thermal error before and after compensation. After compensation in the X-direction, the displacement sensor S ₃ is measured again and compared with the measurements before compensation. Similarly, displacement sensor S ₄ in the Y-direction is compared before and after compensation, while the compared object in the Z-direction is S ₅. The compensation parameters of the three models are shown in Table 10. Maximum and mean are the maximum and average absolute value of the thermal error after thermal compensation, respectively, and $\eta$ is the improved accuracy percentage.

Figure 22.

Compensation in the X-direction.

Figure 23.

Compensation in the Y-direction.

Figure 24.

Compensation in the Z-direction.

Table 10.

Compensation result.

Axis	Model	Maximum, µm	Mean, µm	$\eta (\%)$
Z-direction elongation	Measurement	39.6	24.6
	BP	8.0	7.4	70
	MLRA	24.0	9.1	63
	TS	7.0	2.8	89
X-direction S3	Measurement	26.3	14.5
	BP	15.1	5.2	64
	MLRA	9.3	3.6	75
	TS	9.7	3.3	77
Y-direction S4	Measurement	12.1	8.3
	BP	3.4	1.2	85
	MLRA	8.2	2.0	74
	TS	3.7	1.1	88

BP: back propagation; MLRA: multiple linear regression analysis; TS: time series.

Table 10 indicates that the TS model performed better than BP or MLRA. After compensation for the thermal errors, the error was reduced significantly in the axial direction. For TS, the maximum error decreased from 39.6 to 7 µm in the axis, and the average error decreased from 24.6 to 2.8 µm. The axial accuracy improved by 89%, which effectively demonstrates the proposed measurement and modelling method. BP and MLRA improved the accuracy by 70% and 63%, respectively. The absolute average value of the thermal error S ₃ decreased from 14.5 to 3.3 µm in the X-direction, and the accuracy improved by 77% based on TS. In contrast, BP could improve the accuracy by 64% in the X-direction, while MLRA improved it by 75%. The accuracy was enhanced by 88%, 85% and 74% in the Y-direction using TS, BP and MLRA, respectively.

Conclusion

In this article, three types of thermal-induced error models with prediction-goodness evaluation criteria were constructed, with the aim of solving the problem that traditional spindle thermal error compensation ignores the radial tilt angle errors and the cutting tool length and does not improve the accuracy of the spindle in three dimensions. The error offsets in three directions were compensated with perfect accuracy improvement. A system of technical solutions was proposed for the measurement, modelling and compensation of the thermal error.

1. A synchronous acquisition system was developed, including a thermal deformation collection unit based on five eddy current displacement sensors and a multi-point temperature measurement unit with platinum resistance, which could simultaneously obtain real-time status information of the position-pose of the thermal error and the temperatures of the critical points. This information provides an in-depth analysis of the thermal behaviour of the motorised spindle. In addition, the five-point measurement method could effectively analyse changes in the axial expansion and radial thermal tilt angle errors.

2. A comprehensive thermal offset equation was presented, considering the radial tilt angles, elongation and cutting tool length, which can accurately describe the real position and pose of the spindle thermal deformation and efficiently improve the machine accuracy through compensation.

3. A compensation module was developed, and the compensated strategy was provided based on the feedback integral method and the secondary development of Siemens 840D, realising real-time on-line compensation.

4. A novel method was proposed to group temperature variables, using fuzzy cluster theory to packet temperature variables, and typical temperature variables were selected through statistical correlation analysis. This method can improve the efficiency of the model prediction with minimum independent variables and reduce modelling costs. Comparison of the three thermal error models demonstrated that the BP neural network has a desired fitting ability, but its generalisation is far less than that of the TS model. While the traditional MLRA model has a simple structure, it is a compromise choice between forecast accuracy and modelling cost. In contrast, the TS model can fully exploit the inherent dynamic characteristics of the thermal errors, which is more suitable to be used on the spindle thermal error compensation.