Thermal error compensation of a 5-axis machine tool using indigenous temperature sensors and CNC integrated Python code validated with a machined test piece

Achieving high workpiece accuracy is the long-term goal of machine tool designers. There are many causes for workpiece inaccuracy, with thermal errors being the most common. Indirect compensation (using prediction models for thermal errors) is a promising strategy to reduce thermal errors without increasing machine tool costs. The modelling approach uses transfer functions to deal with this issue; it is an established dynamic method with a physical basis, and its modelling and calculation speed are suitable for real-time applications. This research presents compensation for the main internal and external heat sources affecting the 5-axis machine tool structure including spindle rotation, three linear axes movements, rotary C axis and time-varying environmental temperature influence, save for the cutting process. A mathematical model using transfer functions is implemented directly into the control system of a milling centre to compensate for thermal errors in real time using Python programming language. The inputs of the compensation algorithm are indigenous temperature sensors used primarily for diagnostic purposes in the machine. Therefore, no additional temperature sensors are necessary. This achieved a significant reduction in thermal errors in three machine directions X , Y and Z during verification testing lasting over 60 h. Moreover, a thermal test piece was machined to verify the industrial applicability of the introduced approach. The results of the transfer function model compared with the machine tool ’ s multiple linear regression compensation model are discussed.


Introduction
The heat generated by moving axes and machining processes causes thermal gradients inside the machine tool (MT) structure, resulting in the thermal elongation and bending of MT elements; this substantially deteriorates the MT accuracy [1]. Generally, a solution to this thermal issue falls into one of three categories [2]: � design of the MT structure to reduce its sensitivity to heat flow (e.g. Ref. [3]), � temperature control of the MT and its environments (e.g., using cooling [4] or external heat loads [5]) or � compensation of thermal errors (e.g. Ref. [6]). Furthermore, two alternatives for the compensation of thermal errors at the tool centre point (TCP) can be classified [7]: � direct compensation and � indirect (software) compensation.
In direct compensation, the thermal displacements at the TCP are intermittently measured and superposed to the desired position value of the particular axis, as shown in Brecher et al. [7]. The disadvantage of the direct approach is the required interruption of the machining process to measure the displacement and the additional costs of the measuring probe (e.g., touch probes). The principle of indirect (software) compensation is the readjustment of the positioning of the axes by the MT's control system, which here is based on corrections computed by a mathematical model. Thermally induced displacements at the TCP cannot be sufficiently eliminated through proper MT structure design, temperature control or direct compensation without high additional costs and (or) modelling effort (see Meo et al. [8]). On the contrary, indirect (software) compensation of thermally induced displacements at the TCP is one of the most widely employed techniques to reduce MT thermal errors because of its cost-effectiveness and ease of implementation.
Numerous mathematical approaches can be applied in software compensation. The latest research on approaches to thermal error modelling can be found in Mayr [1]. A review focused on spindle thermal error compensation, which is ordinarily a dominant heat source influencing MT accuracy, can be found in Li [9]. Basically, the thermal error model for indirect compensation can be based on either the numerical simulation results or empirical modelling, where the model is built on the experimental measurements of the MT [9].
In the case of numerical simulations, partial differential equation (PDE) [10], finite difference method (FDM), finite element method (FEM) or a combination of them (FDEM) can be applied to investigate the thermal deformation at the TCP [11]. For example, an FE model-based controller design for the compensation of thermally induced displacements has been reported in Oetinger et al. [12]. However, the creation of a complex and nonstationary predictive numerical model of MT thermal errors is complicated because of many factors (e.g., heat transfer through the joint, establishing the boundary conditions for ambient influence, etc.) and can be a time-consuming process. Therefore, most indirect (software) compensation strategies are based on measured auxiliary variables (empirical models).
A regression analysis is one of the most widely used methods of thermal error modelling [13]. In Srinivas et al. [14], the multiple linear regression (MLR) model for software compensation of thermal displacements was applied. Generally, these MLR models are established in the form of an empirically calibrated polynomial expression. MLR models are often used by MT control system producers for their simplicity despite their weak reaction to a variety of MT conditions and the activity of major heat sources (working cycles).
A general approach to the thermal error modelling of MTs based on transfer functions (TFs) is employed in the current research. TFs contain the nature of the heat transfer principles [15] and use a minimum of additional gauges when compared with models obtained by MLR, artificial neural networks (ANN) [16] and so forth. Brecher et al. [17] showed that even data from a numerical control (NC) system (such as spindle speed or spindle load) are possible inputs into TF models. Internal NC signals are also employed in thermal error models for the compensation of rotary axes thermal errors on 5-axis machines [18]. However, the developed compensation grey-box model is restricted to a description of rotary and swivelling axes and how they impact thermally induced errors; here, other effects are neglected (e.g., the influence of the environmental temperature change, thermal displacements as a result of heat induced by rotation of the main spindle, linear axes movement, etc.).
Despite the fact that different heat sources are permanently combined in real machining conditions, which cause complex thermal errors at the TCP, much of the research has focused on an approximation of only one active heat source in the overall MT thermal errors. This is a typical drawback of many existing compensation models of thermally induced errors; for example, only spindle thermal error compensation is considered in Li [9], time-varying environmental temperature influence in Zhang [19], the rotary axes movement in Blaser [20] and Bitar-Nehme [21] and so forth.
The current paper presents compensation for the main internal and external heat sources affecting the 5-axis MT structure. Indigenous temperature sensors placed close to the main heat sources, such as spindle bearings, feed drives, the rotary table motor and ambient temperatures (used primarily for machine diagnostic) are used to calculate the resulting thermally induced displacements at the TCP in the three machine directions, X, Y and Z.
MT thermal errors are occasionally termed nonlinear in the available literature (see Zhu [22]), and nonlinear models are used for the descriptions. The method used in the present research is based on the linearisation of the MT thermal issues, with the help of an equation system consisting of linear parametric models, which is contrary to models with nonlinear parts. The basic TF model concept has already been introduced for different MT structures (3-axis or 4-axis MT), e.g. Ref. [23,24]. TF model is extended to include a table rotation (C axis) of 5-axis MT. Thus, the compensation TF model consider spindle rotation, three linear axis movements, rotary C axis and time-varying environmental temperature influence, save for the cutting process. Moreover, the model predicts thermally induced displacements at the TCP in the three machine directions, X, Y and Z. The model is implemented directly into the MT control system via Python programming. Furthermore, the new model was tested using the following two approaches: � a long-term verification experiment carried out under different conditions from the calibration tests and � machining a thermal test piece, comparing the approach with an indigenous MLR model composed by the MT manufacturer.
The tested MT and experimental setup for thermo-mechanical measurements are explained in the next section. Section 3 deals with the mathematical mechanism, thermal error model structure and approximation quality expression. Furthermore, section 3 involves the calibration of TFs and identification procedures, including a described thermomechanical system analysis and model implementation into the MT control system. The verification testing and results are elucidated in section 4. The verification testing consists of two parts, one dedicated to measurement according to international standards and the second to thermal test piece manufacturing. In section 5, a critical discussion is conducted, illustrating the advantages of a compensation model with an open structure and modularity. A brief summary and future scope of the work regarding this approach are presented in section 6.

Experimental setup
The current research is focused on minimising of thermally induced errors at the TCP in machine directions X, Y and Z because of the rotation of the main spindle, linear and rotational axis movements and ambient temperature changes. All experiments were performed on a gantry-type 5-axis milling centre with a rotary table (diameter 630 mm, no turning operations). The MT was equipped with a number of internal temperature probes primarily used for safety and diagnostic purposes and that were installed in the MT control system by the MT manufacturer (the temperatures are recorded in 0.1 � C and 0.5 � C resolutions). Six of the internal probes (placed closed to the spindle bearings T sp , X, Y and Z linear axis motors T X , T Y and T Z , rotary table motor T C along with reference temperature T base ) were utilised in the thermal deformation modelling process. The tested machine is schematically shown on the left side of Fig. 1 alongside the approximate positions of the internal temperature probes used in the compensation model. Eddy current sensors firmly clutched in the measuring fixture were employed for noncontact sensing of the displacements between the TCP (in the directions X, Y and Z shown on the right side of Fig. 1) represented by a test mandrel (length 125 mm, diameter 40 mm) and the rotary table (regular position of the workpiece), per ISO 230-3 international standard [25]. The displacements were sensed in micrometre resolution.
All of the results and conclusions are associated with the following experiment conditions: the model was valid for no-load or finishing conditions, the model was calibrated in one configuration of the MT axes (at the centre of the rotary table), and compensation was realised in linear directions only (an investigation of angular deformations is beyond the scope of the current paper).

Compensation model for thermal errors
The model concept lies in the partial linearisation of the issue. This means isolating particular thermal elements, solving them separately and building an approximation model with their subsequent superposition [17,26,27]. A set of TFs each using only one input signal (temperature measured close to a relevant heat source) was deemed a suitable tool to fulfil the aforementioned intent.
Mathematical and modelling approach. A discrete TF was used to describe the link between the excitation and its response (equation (1) and (2)).
yðtÞ ¼ ε:uðtÞ þ eðtÞ; The vector u(t) in equations (1) and (2) is the TF input in the time domain, y(t) is the output vector in the time domain, ε represents the TF in the time domain, e(t) is the disturbance value (further neglected), a n is the calibration coefficient of the TF input, b m is calibration coefficient of the TF output, n is the order of the TF numerator, m is the order of the TF denominator, and z is a complex number. The difference form of a discrete TF (a generally suitable form for modern MT control systems using their programming languages) in the time domain is introduced in equation (3).
where k represents the examined time instant and k-n (k-m) is the nmultiple (m-multiple) delay in sampling frequency of the measured input vector (simulated output vector). The linear parametric models of ARX (autoregressive with external input) or OE (output error) identifying structures were used for an estimation of the TF coefficients. The stability of each TF and relationship between the thermo-mechanical system input and output described by the TF is examined through a linear time invariant (LTI) step response [28]. Excitations in TFs are temperatures measured close to heat sources, and the responses represent the caused deformations in the examined directions.
A general mathematical term of the thermal error approximation model using the superposition principle is possible to express by the following equation: where index i represents the number of axis directions, thermal errors are solved (X, Y and Z in the presented research), j stands for an indication of the number of described heat sources (spindle rotation, three linear axis movements, one rotational axis movement and ambient temperature impact in the case of the presented research), and l is the total number of described heat sources by the compensation model. The deformational elements δi,j in the overall thermal error δi are solved as follows: The term δ i in equation (4) represents the desired approximation value of the MT thermal error in the relevant direction. The term ΔT j in equation (5) expresses the relevant temperature in relative coordinates measured close to a heat source, and ϵ i,j represents the TF in a time domain describing the relation between the relevant input and output (temperature T j and deformational element δ i,j ) of the thermomechanical system.
Based on equations (4) and (5), the approximation of the MT thermal errors can by described as follows: Taking the prediction for the Z direction as an example, the mathematical expression of the specific model is shown in equation (7).
An approximation quality of the simulated behaviour is expressed in equation (8). This value represents the percentage of the output variation reproduced by the model [28].
The δ mea value in equation (8) is the measured output (thermal deformation), δ sim is the simulated/predicted model output, and δ sim expresses the arithmetic mean over time of the measured output. The vector norm used in equation (8) is generally expressed as follows: δ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where δ is a general vector of the length r.

Calibration and identification of MT thermo-mechanical system
Generally, model calibration should describe a transient characteristic between two thermodynamic equilibria: one with the MT ambient temperature and one with all heat sources meant to be approximated.
The described compensation model based on TFs concentrates on the thermal deformations caused by spindle rotation, the rotation movement of table (C axis), linear movements in the X, Y and Z axes and ambient temperature changes.
Compensation for the thermal errors is realised in the X, Y and Z machine directions. Seven calibration measurements were carried out to obtain the required information about the contribution of each heat source to the overall MT thermal error. The calibration measurement setups are summarised in Table 1. The table shows the maximal values (feed rates and rotation speeds) of the described heat sources taken from the MT datasheet (first row in Table 1). The feed rates and rotation speeds for the calibration tests were suggested by the MT manufacturer from the safety of the MT operation point of view (second row in Table 1). Higher values of axis movements and rotations are generally used for short-term fast feeds (except for the spindle).
The ETVE (environmental temperature variation error) test per international standard ISO 230-3 [25] was used for the calibration measurement of ambient temperature influence.
A spindle rotation could be accompanied by the differences between its heating and cooling phases, which is mainly caused by the presence of forced convection in the heating phase [26]. The forced convection is produced by the rotating parts of the spindle. Therefore, considering both the heating and cooling phase is beneficial when describing the impact of the spindle rotation. A total of 20 h (10 for heating, 10 for cooling; see fourth row in Table 1) are sufficient for an understanding of the potential differences between both phases.
Calibration of the heat source impact caused by the spindle speed allows for continuous positioning of the test mandrel in the measuring fixture, whereas a measurement of the axis movements do not. Calibration measurement of movement in the axes (X, Y and Z) and table (C) rotation is divided into two cyclically repeated sections: measurement of test mandrel thermal displacements (10 s in the measuring fixture) and a 5 min repeated load application (the test mandrel is moved from the measuring fixture to prevent a collision). In the second section, the MT conducts a cyclic linear (X, Y or Z axis) or rotational (C axis) motion in trajectories (third row in Table 1) that are under a specific constant speed. Both, the heating and cooling phases are measured. Contrary to the spindle activity, the linear dependency of the inputs and outputs during X, Y, Z and C axis movements are more expected. The calibration time ranging from 12 to 20 h is examined to find an adequate calibration length of X, Y, Z and C axis movement description in thermal error compensation model.
All measured inputs and outputs during the calibration experiments of the MT thermo-mechanical system are summarised in the left and middle parts of Fig. 2. The simulated outputs of the identified TFs (black dash lines) are presented in the middle parts of Fig. 2 alongside the measured deformational behaviours. The sets of identified TFs used in the approximation model for the Z direction were selected according to the best fit value through the estimation process and examined through LTI step response [28]. The LTI expresses a clear dependency between an excitation (1 � C sudden change in temperature input) and a response (predicted deformation) of the thermo-mechanical system. Moreover, the stability of each TF can be proven. LTI diagrams for the transfers used in equation (7) are shown on the right side of Fig. 2. The measured outputs and predicted results in the other directions (X and Y) are not presented in the current paper. Note that the entire data processing and modelling parts of the research lasted up to 60 h (calibration measurements; as shown in Table 1).
Heat sources caused by the ambient temperature, spindle speed and table rotation (C axis) are more sensitive to temperature changes than other heats sources, as observed from the LTI diagrams (right parts of Fig. 2 a, b and f). The deformational responses of heat sources caused by ambient temperature, spindle speed, Z axis movement and table rotation have a reverse impact on the X and Y axis movements. Thus, a certain autocompensation effect during the combined activity of the sources is expected.
The calibration coefficients of the six identified TFs for the approximation of MT thermal errors in the Z direction are shown in Table 2. The approximation quality expressed by the fit values according to (8) of the identified functions reached during the identification process are presented in the last column of Table 2.
The final compensation algorithm is implemented in the MT controller (Heidenhain TNC 640) using the Python OEM process option (using Python 2.7). The Python option for Heidenhain does not allow the use of any external Python libraries, since they cannot be compiled on the PLC unit. The code is thus implemented using standard libraries only. The compensation algorithm is activated via a callback function. This callback function responds to the change of state of a PLC variable oscillating at a frequency of 1 Hz. Upon activation, the values of PLC variables corresponding to the relevant sensors are loaded into the model and the predicted compensation values are evaluated and loaded into the PLC variables responsible for the actual positional adjustments in each axis. The model is represented by a class instance in memory, which stores current transfer function history information and the time of the last successful evaluation and checks input and output validity. The algorithm evaluation process does not exceed 40 ms, leaving enough time for possible future generalisations of the model and/or more frequent evaluation.

Verification
The setup of the two experiments concerning the verification of the implemented model (equation (7)) into the MT control system was chosen for several reasons:  (7)) � examination of the model's validity during the combined activity of the calibrated heat sources, � the model's sensitivity to changes in heat sources and their magnitude compared with calibration measurements and � a comparison of an alternative method [29] to the thermal deformation measurement according to ISO standards [25].

Combined activity of calibrated heat sources
The progress and setup of the first verification testing of the implemented compensation model is presented in Fig. 3. In the figure, the upper chart shows the input signals' behaviour over time, including spindle speed, during the test. The adjacent table in the same figure summarises the magnitude and activity of the individual heat sources caused by movements (linear and rotational) in the machine axes. The whole experiment can be divided into three parts, each followed by a cooling phase. A particular spectrum is divided into sections combining axis movements (or conjugated axis movements) and spindle rotation. Measurements in the sections are conducted analogously to the calibration measurements of the axis movement impacts described in the previous section. The feed rates in the X, Y and Z and table rotation speed (C axis) are set to lesser or equal values as in the calibration tests (see Table 1). The spectrum of spindle speed varies through the complete rpm range.
The deformations measured after compensation δ Z (comp) (residual error of the TF compensation model) during the first verification testing in the Z direction are depicted in Fig. 4, along with information regarding the uncompensated state δ Z (no-comp). Deformations in the uncompensated state δ Z (no-comp) are calculated as the sum of the measured deformations after compensation δ Z (comp) and deformations Table 2 Calibration coefficients of the identified TFs for the Z direction thermal error approximation.   predicted by the TF model δ Z (sim). The developed model, per (7), can be divided in three parts (submodels). The first submodel describes only the ambient temperature influence δ Z,amb (sim), the second submodel characterises only spindle impact δ Z,sp (sim), and the last submodel represents the influence of axis movements (X,Y, Z and C) δ Z,mov (sim). The approximation model analysis of the contribution of deformational elements in the Z direction thermal error compensation during the first verification testing is depicted also in Fig. 4. The improvement of the MT thermo-mechanical state (expressed by the fit value according to (8)) during the first verification testing was 73% in the Z direction compared with the uncompensated state. Spindle activity had a dominant impact on MT thermo-mechanical behaviour in the Z direction over the other described thermal sources. The submodel describing ambient temperature impact followed the long-term gradient of overall thermal error.
A similar structure to the thermal error model of the Z direction (equation (7)) was applied for the X and Y machine directions with the same temperature inputs. The TF calibration coefficients approximating the X and Y directions were identified using the process described in section 3. Note that the X axis movement has a negligible impact on both the X and Y direction thermal errors. Hence, the deformational elements δ X,X , and δ Y,X were not considered during modelling. Models of both directions were implemented into the MT control system using Python programming, as mentioned above.
The measured resultant deformations and reconstructed MT uncompensated state along with the approximation model analysis of the contribution of deformational elements in the X and Y directions during the first verification testing are shown in Fig. 5 and Fig. 6, respectively.
The MT accuracy improvements were 62% in the X direction and 56% in the Y direction when compared with the uncompensated state. The impact of spindle activity and the influence of axis movements on MT thermo-mechanical behaviour in both the X and Y directions, contrary to the Z direction, were approximately similar.
Thermal errors at the TCP in the X and Y directions of the tested gantry-type 5-axis milling centre were lower compared with thermally induced displacements at the TCP in the Z direction. However, thermal errors in the X and Y directions were still significant and should not be neglected during the development of a thermal error compensation model.

Thermal test piece
Thermal test piece machining [29] is an alternative to the thermal displacement measurement given by international standards [25] (i.e., using displacement sensors). The thermal test piece material was an aluminium alloy. The surface was covered with an eloxal coating for better visibility of cutting tool imprints on the surface caused by thermally induced displacements at the TCP. The thermal test piece fixed on the table in the MT workspace is depicted in Fig. 7. The verification of the compensation model at a different position in the machine work envelope using the thermal test piece is somewhat limited. However [30], demonstrated that the influence of spindle rotation on volumetric errors changes is very small (the maximum difference in the magnitude of volumetric errors was 10 μm). Small differences in volumetric errors are probably due to a limited spread of heat in the headstock structure. Given that the spindle rotation is the heat sources with the relatively most significant impact on the MT structure's thermal deformations in Z direction, the validity of the compensation model in the Z direction can be expected to be sufficient in the entire MT work space. This assumption might be addressed in future studies.
The procedure for creation of the thermal test piece is described in Fig. 8. The process was divided into two cyclically repeated phases:   machining tool imprints in a column and conducting a test with simultaneous spindle rotation and movements along the Z axis (see the specification in Fig. 8). One column of tool imprints was machined with a programmed Z tool positioned, as shown on the left side of Fig. 8 (with 10 μm step up to 130 μm). Thereafter, the test with simultaneous spindle rotation and movements in the Z axis was carried out (2 h). The procedure was repeated until the entire matrix of imprints was machined (10 columns; 18 h), as depicted on the right side of Fig. 8.
The complete matrix would have contained only three rows in the ideal case (without any thermally induced displacements at the TCP). This ideal case corresponds to the blue imprints on the right side of Fig. 8. The red imprints reflect a thermal distortion in the Z direction caused by loading periods over time. The vertical scale of the thermal test piece provides a visual imprint of the TCP thermal deformations under the specific test with simultaneous spindle rotation and movements in the Z axis (with 10 μm resolution, as mentioned above).
The thermal test piece constitutes a means of comparison between two approaches of thermal error minimisation. The test piece consists of three parts: without compensation, with compensation by the developed thermal error model, per (7) and with indigenous compensation based on the MLR method provided by the MT manufacturer. The MLR compensation model was built based on the measured temperatures T sp and T base which were obtained during only one calibration test with constant spindle rotation (8000 rpm, 8 h, without a cooling phase). A general equation of MLR model is as follows: where the term δ MLR in equation (10) represents the desired approximation value of the MT thermal error using the MLR model, ϑ expresses the matrix of relevant temperatures, and β represents the vector of calibration coefficients. The calibration coefficient vector may be obtained using the statistical expression as follows.
The MLR model based on a statistical method (11) is a simple approach with no dynamic terms and low physical interpretation. The results are highly dependent on the number of inputs.
All parts of the test piece were manufactured simultaneously. The MT control system switched the compensations on and off automatically during the NC program.
The choice of heat source combinations during the test exercise was based on the calibration measurement results from Fig. 2. A combination (except for the ambient temperature variation) of the spindle rotation and movement in the Z axis caused the greatest impact on thermal errors at the TCP in the Z direction. Therefore, the loading periods consisted of spindle speed spectrum and simultaneous movement in the Z axis (10 h) followed by a cooling phase of the MT (8 h). The progress and setup of the verification test with the thermal test piece are depicted in Fig. 9.
A plan view of the manufactured thermal test piece with an imprint of the measured thermal deformations is shown in Fig. 10.
The measured uncompensated MT state and predicted values along with the approximation model analysis of the contribution of deformational elements during the test piece manufacturing are shown in Fig. 11.
Reduction of thermal deformations in the Z direction with the aid of the thermal error compensation model using TF was as high as 92% compared with the uncompensated state (achieved reduction by the model application from 130 μm to 10 μm). The good results with the thermal test piece confirmed that the calibration time of the TF model was sufficient even for the Z axis movements with the shortest calibration time of 12 h. The result of MLR compensation deteriorated in accuracy compared with the model based on TF over 43%.
The MLR model originally embedded in the control system of the tested MT does not respect the MT's operation range completely. To adequately compare the MLR model approximation quality together with TF model results, the new MLR model should be built based on the available measured data used for TF model calibration. However, a similar analysis has previously been carried out [24] for the compensation of thermally induced displacement at the TCP in the Z direction on a portal milling centre (C axis influence is not considered in the compensation models). The TF model showed much greater accuracy  and reliability, particularly under highly varied operating conditions during a verification test over 5 days in comparison with the MLR model. The approximation quality expressed by the fit values for the model based on the TF is four-fold higher than is fit in case of the MLR model built based on the available measured data used for TF model calibration.

Conclusion
The main objective of the scientific investigation presented in the current paper was enhancement of MT accuracy through minimisation of thermal errors with the help of partial linearisation of the MT thermomechanical issue in terms of solving each heat source impact separately with their subsequent superposition. The tested machine was a gantrytype 5-axis milling centre with a rotary table (no turning operations) and symmetrical structure. The experiments were carried out under specific conditions: no cutting process was involved, and the calibration was performed at one position of the MT axes.
The linearly superposed transfer function approach is extended to include the C-axis of a gantry-type 5-axis milling. The effects of major heat sources (ambient temperature changes, the spindle rotation and X, Y and Z feed drive activities) are also included. Moreover, the model predicts thermally induced displacements at the TCP in all the relevant MT directions i.e. X, Y and Z (linear deformation components of thermal errors, a solution for the angular components of deformations was not considered) based on input parameters provided by the MT control system (indigenous temperature sensors). The developed compensation model based on TFs has an open structure for easy integration of different thermal source impacts (e.g. cutting process as referenced in Ref. [31]). In addition, the compensation model was implemented directly into the MT control system via Python programming.
The approximation quality of the model based on TFs was verified, as follows: � during verification testing (up to 62 h) by combining all calibrated heat sources, � by machining the thermal test piece to visualise thermally induced errors in the Z direction and � comparing the TF model with the MLR compensation model (developed by the MT manufacturer of the tested milling centre) originally embedded in the MT control system.
The presented findings confirm that TFs are a suitable approach for thermal error modelling and lead to promising results. Thermal error reductions of 62%, 56% and 73% (in the X, Y and Z directions, respectively) in thermal error reduction compared with the uncompensated MT state were achieved.
Further research will focus on the impact of the cutting process on the tool-MT-workpiece chain, as already referred to in Ref. [31], and machining tests with different technological conditions (finishing conditions, different workpiece materials, different tool diameters, changes

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.