Thermal-Feature System Identification for a Machine Tool Spindle

The internal temperature is an important index for the prevention and maintenance of a spindle. However, the temperature inside the spindle is undetectable directly because there is no space to embed a temperature sensor, and drilling holes will reduce its mechanical stiffness. Therefore, it is worthwhile understanding the thermal-feature of a spindle. This article presents a methodology to identify the thermal-feature model of an externally driven spindle. The methodology contains self-made hardware of the temperature sensing and wireless transmission module (TSWTM) and software for the system identification (SID); the TSWTM acquires the temperature training data, while the SID identifies the parameters of the thermal-feature model of the spindle. Then the resulting thermal-feature model is written into the firmware of the TSWTM to give it the capability of accurately calculating the internal temperature of the spindle from its surface temperature during the operation, or predicting its temperature at various speeds. The thermal-feature of the externally driven spindle is modeled by a linearly time-invariant state-space model whose parameters are identified by the SID, which integrates the command “n4sid” provided by the System ID Toolbox of MATLAB and the k-fold cross-validation that is common in machine learning. The present SID can effectively strike a balance between the bias and variance of the model, such that both under-fitting and over-fitting can be avoided. The resulting thermal-feature model can not only predict the temperature of the spindle rotating at various speeds but can also calculate the internal temperature of the spindle from its surface temperature. Its validation accuracy is higher than 98.5%. This article illustrates the feasibility of accurately calculating the internal temperature (undetectable directly) of the spindle from its surface temperature (detectable directly).


Introduction
High-speed machining is a promising technology to dramatically increase productivity and reduce costs. Most machining errors in machine tools come from thermal problems [1,2]. During high-speed machining of a spindle, a lack of complete information of its temperature variation may lead to an abrupt malfunction due to overheating. Therefore, it is very important to monitor and predict the temperature variation of the spindle during the machining process. In order to accurately monitor the temperature inside the spindle, the authors propose a methodology, which combines hardware and software into an edge-computing module, to identify the thermal-feature of the spindle. One can accurately monitor the temperature inside the spindle from its surface temperature through the proposed module. Research on the thermal analysis of spindles has been going on for quite a long  The firmware used to transform the digitalized voltage signals of RTDs into temperature signals is programmed with Arduino. Figure 3a details the flowchart of the signal transform firmware. First, the pins of A0, A2, and A3 on the Bluno Beetle are defined as the general-purpose output (GPO) pins, they output a high or low potential according to the truth table (Table 1) set up in the register. The analog-to-digital converter (ADC) is initially set up with the external reference, gain of two, excitation current of 210 µA, calibration, 4.17 Hz update rate, and continuous read mode. The channel in connection would be identified according to the examination of the resistance value before the signal process stage. After the initialization stage, the micro control unit (MCU) receives and converts the register number of the analog input voltage into binary code to calculate the resistance of the RTD based on Equation (1) provided by the AD7794 datasheet.
(1)  port. The digitalized voltage signal of the RTD is transferred into a temperature signal by the signal transform firmware. Afterward, the temperature signal is received by the CC2540 through the transceiver (TX) and receiver (RX) pins. Finally, the temperature signal is transmitted to the terminal through Bluetooth or USB.  The firmware used to transform the digitalized voltage signals of RTDs into temperature signals is programmed with Arduino. Figure 3a details the flowchart of the signal transform firmware. First, the pins of A0, A2, and A3 on the Bluno Beetle are defined as the general-purpose output (GPO) pins, they output a high or low potential according to the truth table (Table 1) set up in the register. The analog-to-digital converter (ADC) is initially set up with the external reference, gain of two, excitation current of 210 µA, calibration, 4.17 Hz update rate, and continuous read mode. The channel in connection would be identified according to the examination of the resistance value before the signal process stage. After the initialization stage, the micro control unit (MCU) receives and converts the register number of the analog input voltage into binary code to calculate the resistance of the RTD based on Equation (1) provided by the AD7794 datasheet.
(1) The firmware used to transform the digitalized voltage signals of RTDs into temperature signals is programmed with Arduino. Figure 3a details the flowchart of the signal transform firmware. First, the pins of A0, A2, and A3 on the Bluno Beetle are defined as the general-purpose output (GPO) pins, they output a high or low potential according to the truth table (Table 1) set up in the register. The analog-to-digital converter (ADC) is initially set up with the external reference, gain of two, excitation current of 210 µA, calibration, 4.17 Hz update rate, and continuous read mode. The channel in connection would be identified according to the examination of the resistance value before the signal process stage. After the initialization stage, the micro control unit (MCU) receives and converts the register number of the analog input voltage into binary code to calculate the resistance of the RTD based on Equation (1) provided by the AD7794 datasheet. (1) The temperature can be determined by the relationship of temperature and resistance given by Equation (2) from the DIN EN 60751 standard which defines Pt100 resistance accuracy classes and corresponding tolerance, where The actual temperature could be transmitted to the Bluetooth-MCU (BLE-MCU) (CC2540) through the transceiver (TX) and receiver (RX) pins. The temperature can be determined by the relationship of temperature and resistance given by Equation (2) from the DIN EN 60751 standard which defines Pt100 resistance accuracy classes and corresponding tolerance, The actual temperature could be transmitted to the Bluetooth-MCU (BLE-MCU) (CC2540) through the transceiver (TX) and receiver (RX) pins.

Experiment Setup
The run-in system of the spindle, shown in Figure 3b, consists of a spindle externally driven by a 7.5 kW motor through a belt, a controller accompanied by a personal computer (PC), used to

Spindle Run-In System
The run-in system of the spindle, shown in Figure 3b, consists of a spindle externally driven by a 7.5 kW motor through a belt, a controller accompanied by a personal computer (PC), used to control the motor through a programmable control interface, and two sets of self-made TSWTMs with temperature RTD (PT1000) sensor probes. As shown in Figure 3c, RTDs are located at the outer rings of front and rear bearings, the center point of the spindle's inner housing, the rear and front ends of the spindle's outer housing surface, and the ambient temperature to acquire the training data for the thermal feature identification algorithm. Figure 4 shows the photos of the experimental setup. Each run of experiments lasts for almost a day; therefore, the variation in ambient temperature during the day ( Figure 5) must be considered. Based on the literature survey in Introduction, one knows that bearings are the main heat source of spindles, therefore the bearings are important temperature characteristic points and choosing the outer rings of front and rear bearings to set the RTD is important because they are static rings. As conduction dominates the heat transfer of the metallic inner housing and is very fast, one must choose the center point of the inner housing because it is the same distance away from either end. In fact, spindles do not allow holes drilled in them to set the RTD when in use, therefore one must choose the rear and front ends of the outer housing surface to set the RTD because they are the only positions being capable of detecting temperature directly when in use. control the motor through a programmable control interface, and two sets of self-made TSWTMs with temperature RTD (PT1000) sensor probes. As shown in Figure 3c, RTDs are located at the outer rings of front and rear bearings, the center point of the spindle's inner housing, the rear and front ends of the spindle's outer housing surface, and the ambient temperature to acquire the training data for the thermal feature identification algorithm. Figure 4 shows the photos of the experimental setup. Each run of experiments lasts for almost a day; therefore, the variation in ambient temperature during the day ( Figure 5) must be considered. Based on the literature survey in Introduction, one knows that bearings are the main heat source of spindles, therefore the bearings are important temperature characteristic points and choosing the outer rings of front and rear bearings to set the RTD is important because they are static rings. As conduction dominates the heat transfer of the metallic inner housing and is very fast, one must choose the center point of the inner housing because it is the same distance away from either end. In fact, spindles do not allow holes drilled in them to set the RTD when in use, therefore one must choose the rear and front ends of the outer housing surface to set the RTD because they are the only positions being capable of detecting temperature directly when in use.  The performance of the self-made TSWTM was tested in a precisely temperature-controlled environment provided by the Electronics Testing Center (ETC) in Taiwan. The control the motor through a programmable control interface, and two sets of self-made TSWTMs with temperature RTD (PT1000) sensor probes. As shown in Figure 3c, RTDs are located at the outer rings of front and rear bearings, the center point of the spindle's inner housing, the rear and front ends of the spindle's outer housing surface, and the ambient temperature to acquire the training data for the thermal feature identification algorithm. Figure 4 shows the photos of the experimental setup. Each run of experiments lasts for almost a day; therefore, the variation in ambient temperature during the day ( Figure 5) must be considered. Based on the literature survey in Introduction, one knows that bearings are the main heat source of spindles, therefore the bearings are important temperature characteristic points and choosing the outer rings of front and rear bearings to set the RTD is important because they are static rings. As conduction dominates the heat transfer of the metallic inner housing and is very fast, one must choose the center point of the inner housing because it is the same distance away from either end. In fact, spindles do not allow holes drilled in them to set the RTD when in use, therefore one must choose the rear and front ends of the outer housing surface to set the RTD because they are the only positions being capable of detecting temperature directly when in use.  The performance of the self-made TSWTM was tested in a precisely temperature-controlled environment provided by the Electronics Testing Center (ETC) in Taiwan. The  The performance of the self-made TSWTM was tested in a precisely temperature-controlled environment provided by the Electronics Testing Center (ETC) in Taiwan. The constant-temperature environment is controlled by liquid calibration baths (LCB) and the standard temperature is measured by a high-accuracy (@-100 to 100 • C ± 0.009 • C) platinum RTD thermometer [20]. Given the results of channels 1 and 2 of the TSWTM as examples, Figure 6a,b show their sensing error over the temperature range of 25 to 55 • C. The sensing error at different temperatures can be fitted with a first-order linear equation. The slope of error to temperature of channels 1 to 5 are −0.0056, −0.0058, −0.0061, −0.0051, and −0.0054 respectively, which demonstrate that the five channels of TSWTM are consistent with each other. Therefore, the differences in sensing error between the five channels are negligible. The ordinates of Figure 7a,b are the temperature measured by the TSWTM at the standard temperature of 25 and 55 • C, respectively, while their abscissas are the elapsed times, which demonstrate the high stability and good anti-interference ability of the TSWTM. Tables 2 and 3 list the specifications of the TSWTM and temperature probe, respectively. The bare size of the TSWTM is 43 × 48 × 10 mm 3 and the package size is Ø 80 × 53 mm 3 . The TSWTM measures the temperature with a temperature probe (PT1000, Ø3 probe). Its measurement range is −50 to 300 • C, the accuracy is @ 25 ± 0.05 • C, and the power consumption is 175 mW. Figure 8 shows the components in use and assembly process of the TSWTM. constant-temperature environment is controlled by liquid calibration baths (LCB) and the standard temperature is measured by a high-accuracy (@-100 to 100 °C ± 0.009 °C) platinum RTD thermometer [20]. Given the results of channels 1 and 2 of the TSWTM as examples, Figure 6a,b show their sensing error over the temperature range of 25 to 55 °C. The sensing error at different temperatures can be fitted with a first-order linear equation. The slope of error to temperature of channels 1 to 5 are − 0.0056, − 0.0058, − 0.0061, − 0.0051, and − 0.0054 respectively, which demonstrate that the five channels of TSWTM are consistent with each other. Therefore, the differences in sensing error between the five channels are negligible. The ordinates of Figure 7a,b are the temperature measured by the TSWTM at the standard temperature of 25 and 55 °C, respectively, while their abscissas are the elapsed times, which demonstrate the high stability and good anti-interference ability of the TSWTM. Tables 2 and  3 list the specifications of the TSWTM and temperature probe, respectively. The bare size of the TSWTM is 43 × 48 × 10 mm 3 and the package size is Ø 80 × 53 mm 3 . The TSWTM measures the temperature with a temperature probe (PT1000, Ø3 probe). Its measurement range is −50 to 300 °C, the accuracy is @ 25 ± 0.05 °C, and the power consumption is 175 mW. Figure 8 shows the components in use and assembly process of the TSWTM.   constant-temperature environment is controlled by liquid calibration baths (LCB) and the standard temperature is measured by a high-accuracy (@-100 to 100 °C ± 0.009 °C) platinum RTD thermometer [20]. Given the results of channels 1 and 2 of the TSWTM as examples, Figure 6a,b show their sensing error over the temperature range of 25 to 55 °C. The sensing error at different temperatures can be fitted with a first-order linear equation. The slope of error to temperature of channels 1 to 5 are − 0.0056, − 0.0058, − 0.0061, − 0.0051, and − 0.0054 respectively, which demonstrate that the five channels of TSWTM are consistent with each other. Therefore, the differences in sensing error between the five channels are negligible. The ordinates of Figure 7a,b are the temperature measured by the TSWTM at the standard temperature of 25 and 55 °C, respectively, while their abscissas are the elapsed times, which demonstrate the high stability and good anti-interference ability of the TSWTM. Tables 2 and  3 list the specifications of the TSWTM and temperature probe, respectively. The bare size of the TSWTM is 43 × 48 × 10 mm 3 and the package size is Ø 80 × 53 mm 3 . The TSWTM measures the temperature with a temperature probe (PT1000, Ø3 probe). Its measurement range is −50 to 300 °C, the accuracy is @ 25 ± 0.05 °C, and the power consumption is 175 mW. Figure 8 shows the components in use and assembly process of the TSWTM.

System Identification for the Thermal-Feature Model of the Spindle
During the operation of the spindle, the temperature variations within it involves complex heat transfer phenomena and unpredictable environmental disturbances. As a result, it is difficult to analyze its heat transfer phenomena by means of theoretical approaches [21]. On the other hand, SID is a method for building the mathematical model of a system from the data measured during its operation [22]. The authors model the thermal-feature of an externally driven spindle in terms of a state-space representation based on the input rotational speed and the output temperature datasets measured during its operation. The process of SID for the thermal-feature model of the spindle can be roughly divided into three stages: data preparation, model structure determination, and model parameter identification, as shown in Figure 9. The three stages are explained in the following: 1. Data preparation: This stage prepares the dataset for SID. More specifically, this stage has to remove the outliers, perform resampling, remove missing data, and justify the starting time for all the time-series sequences.

System Identification for the Thermal-Feature Model of the Spindle
During the operation of the spindle, the temperature variations within it involves complex heat transfer phenomena and unpredictable environmental disturbances. As a result, it is difficult to analyze its heat transfer phenomena by means of theoretical approaches [21]. On the other hand, SID is a method for building the mathematical model of a system from the data measured during its operation [22]. The authors model the thermal-feature of an externally driven spindle in terms of a state-space representation based on the input rotational speed and the output temperature datasets measured during its operation. The process of SID for the thermal-feature model of the spindle can be roughly divided into three stages: data preparation, model structure determination, and model parameter identification, as shown in Figure 9. The three stages are explained in the following:

1.
Data preparation: This stage prepares the dataset for SID. More specifically, this stage has to remove the outliers, perform resampling, remove missing data, and justify the starting time for all the time-series sequences. Model parameter identification: This stage is for identifying the best parameters of the thermal-feature model based on the measured input/output data of the spindle in operation. After the identification, the resulting thermal-feature model is validated by the test data to verify whether the model's prediction is accurate enough. If the accuracy is not good enough, one has to change the method of parameter identification. Furthermore, once the change of the parameter identification method cannot generate satisfactory accuracy, one has to go back to the previous stage to re-determine the model structure. In fact, the stages of model structure determination and parameter identification are often repetitively interleaved until the model with the best structure and parameters is found. It should be mentioned that the test data used to validate the model is always different from the training data used to identify the model. 2. Model structure determination: This stage is for deciding several important model properties, such as linear or nonlinear, and time-invariant or time-varying, and then determine the model structure, namely the types and corresponding structure configurations accordingly. 3. Model parameter identification: This stage is for identifying the best parameters of the thermal-feature model based on the measured input/output data of the spindle in operation. After the identification, the resulting thermal-feature model is validated by the test data to verify whether the model's prediction is accurate enough. If the accuracy is not good enough, one has to change the method of parameter identification. Furthermore, once the change of the parameter identification method cannot generate satisfactory accuracy, one has to go back to the previous stage to re-determine the model structure. In fact, the stages of model structure determination and parameter identification are often repetitively interleaved until the model with the best structure and parameters is found. It should be mentioned that the test data used to validate the model is always different from the training data used to identify the model.

Data Preparation Stage
Recall that the input and output data of the thermal-feature model of the spindle are its rotational speeds and temperatures, respectively. The rotational speed of the spindle is controlled by a programmable control interface, while the temperatures are measured by the self-made TSWTM. The spindle was set to rotate at a certain constant speed, meanwhile its temperature variation was recorded continuously. The spindle was turned down to stationary if the temperature variation was less than 1 ℃ for 3 consecutive hours. The temperature was recorded continuously until the spindle naturally cooled to ambient temperature. The whole process lasted a full day for each specific speed of the spindle, the records are detailed in Table 4. Figure 10 shows the raw data of 4000 rpm as an example.
As the temperatures at the different locations in and on the spindle are acquired by different channels of the TSWTM, the temperature sampling timings of different locations are not synchronous. The center symbols shown in the zoom-in window of Figure 10 demonstrate that the temperature raw data acquired by different channels are in different sampling times. For system identification, the raw data has to be cleaned and synchronized through the following procedure: 1. Plot the raw data for visual inspection. If there some extremely erroneous data exists, which may be due to device failure or some unknown reasons, then consider such data to be an outlier and remove it all directly. 2. As the sampling timings of all channels are nonsynchronous, it is necessary to resample the all-time series of output temperatures and input speeds to reach the same sample rate of 1 Hz.

Data Preparation Stage
Recall that the input and output data of the thermal-feature model of the spindle are its rotational speeds and temperatures, respectively. The rotational speed of the spindle is controlled by a programmable control interface, while the temperatures are measured by the self-made TSWTM. The spindle was set to rotate at a certain constant speed, meanwhile its temperature variation was recorded continuously. The spindle was turned down to stationary if the temperature variation was less than 1 • C for 3 consecutive hours. The temperature was recorded continuously until the spindle naturally cooled to ambient temperature. The whole process lasted a full day for each specific speed of the spindle, the records are detailed in Table 4. Figure 10 shows the raw data of 4000 rpm as an example.
As the temperatures at the different locations in and on the spindle are acquired by different channels of the TSWTM, the temperature sampling timings of different locations are not synchronous. The center symbols shown in the zoom-in window of Figure 10 demonstrate that the temperature raw data acquired by different channels are in different sampling times. For system identification, the raw data has to be cleaned and synchronized through the following procedure: 1.
Plot the raw data for visual inspection. If there some extremely erroneous data exists, which may be due to device failure or some unknown reasons, then consider such data to be an outlier and remove it all directly.

2.
As the sampling timings of all channels are nonsynchronous, it is necessary to resample the all-time series of output temperatures and input speeds to reach the same sample rate of 1 Hz. 3.
The data of not a number (NaN) has to be removed from the time series because it represents missing-data. Furthermore, the all-time series should start simultaneously with the first time sample that has all non-NaN temperatures.
Given the spindle speed of 4000 rpm as an example, Figure 10 shows the measured raw data, wherein each temperature curve has a different and possibly uneven sampling rate. Figure 11 demonstrates that the data resampled from the raw data have the same sample rate of 1 Hz. For system identification, one still needs to remove the NaN data and subtract the ambient temperature from the measured temperature data, and adjust the starting time. The final data ready for system identification is shown in Figure 12, wherein the rotational speed and the temperature serve as the input and output data, respectively, for system identification. 3. The data of not a number (NaN) has to be removed from the time series because it represents missing-data. Furthermore, the all-time series should start simultaneously with the first time sample that has all non-NaN temperatures.
Given the spindle speed of 4000 rpm as an example, Figure 10 shows the measured raw data, wherein each temperature curve has a different and possibly uneven sampling rate. Figure 11 demonstrates that the data resampled from the raw data have the same sample rate of 1 Hz. For system identification, one still needs to remove the NaN data and subtract the ambient temperature from the measured temperature data, and adjust the starting time. The final data ready for system identification is shown in Figure 12, wherein the rotational speed and the temperature serve as the input and output data, respectively, for system identification.

Structure Determination Stage
When a system is too complicated to be modeled based on a theoretical approach, one can model it mathematically based on the correlation between input and output data. For the thermal-feature model of the spindle, the input is its rotational speed while the output is its temperature variation measured by the TSWTM. Conceptually, a model can be classified into linear or nonlinear and time-invariant or time-varying. Figure 13 shows the temperature variations of seven spindle's characteristic points at the same rotational speeds to that of Figure 10b, which were measured on different dates. These figures reveal that the temperature progressions do not change much on different dates and thereby demonstrate the time-invariance of the thermal-feature model of the spindle. According to modern control theory, a dynamic system can be represented by a

Structure Determination Stage
When a system is too complicated to be modeled based on a theoretical approach, one can model it mathematically based on the correlation between input and output data. For the thermal-feature model of the spindle, the input is its rotational speed while the output is its temperature variation measured by the TSWTM. Conceptually, a model can be classified into linear or nonlinear and time-invariant or time-varying. Figure 13 shows the temperature variations of seven spindle's characteristic points at the same rotational speeds to that of Figure 10b, which were measured on different dates. These figures reveal that the temperature progressions do not change much on different dates and thereby demonstrate the time-invariance of the thermal-feature model of the spindle. According to modern control theory, a dynamic system can be represented by a

Structure Determination Stage
When a system is too complicated to be modeled based on a theoretical approach, one can model it mathematically based on the correlation between input and output data. For the thermal-feature model of the spindle, the input is its rotational speed while the output is its temperature variation measured by the TSWTM. Conceptually, a model can be classified into linear or nonlinear and time-invariant or time-varying. Figure 13 shows the temperature variations of seven spindle's characteristic points at the same rotational speeds to that of Figure 10b, which were measured on different dates. These figures reveal that the temperature progressions do not change much on different dates and thereby demonstrate the time-invariance of the thermal-feature model of the spindle. According to modern control theory, a dynamic system can be represented by a state-space model [23], which is usually sufficient to describe a dynamic system accurately [24]. Therefore, the authors chose a linearly time-invariant state-space model provided by the System ID Toolbox of MATLAB [24] to model the thermal-feature of the spindle in the following parameter identification stage. state-space model [23], which is usually sufficient to describe a dynamic system accurately [24]. Therefore, the authors chose a linearly time-invariant state-space model provided by the System ID Toolbox of MATLAB [24] to model the thermal-feature of the spindle in the following parameter identification stage.

Parameter Identification Stage
According to modern control theory [23], the linear time-invariant model for model physical system can be expressed in terms of the state-space equation where u(t) is a vector of input data whose components are the time-series of all input data, y( vector of output data whose components are the time-series of all output data, e(t) is the vec noise signals that cannot be reduced during the system identification process, x(t) is the vec state variables whose dimension is dependent on the number of order (n) of the model, an entries of matrices, A, B, C, D, and K, are the parameters of the model to be identified to fit the training data during the identification process and whose dimensions are dependent on the nu of order of the model as well as the number of input/output data. For the case of the spindle i article, the input data is only its time-series of rotational speed (single input), while the output are the time-series of temperature variations of its seven characteristic points (7 outputs) show Figure 3c. Therefore, u(t) is a 1-dimensional vector whose component is the time-series of rotat speed, y(t) is a 7-dimensional vector whose components are the time-series of the temper variations of the 7 characteristic points, and e(t) is a 7-dimensional vector. There are several commands, such as n4sid, ssest, and ssregest, in the System ID Toolb MATLAB for system identification [24]. As the first one, n4sid, is the fastest compared wit others, the authors use it throughout this article. The training data, i.e., the fold 1 of Table 4, wa into the command n4sid to conduct the system identification process. The command n4sid a one to set a range of the model order (n) and determine the best model order within this ran the Hankel singular value. The Hankel singular values measure the contribution of each state t input/output behavior, namely the states with small Hankel singular values can be discard simplify the model. As a result, n4sid determines the best model order to be 8, as shown in F 14. Therefore, with the 8-order state space model, the state space equation of the spind expressed as

Parameter Identification Stage
According to modern control theory [23], the linear time-invariant model for modeling a physical system can be expressed in terms of the state-space equation where u(t) is a vector of input data whose components are the time-series of all input data, y(t) is a vector of output data whose components are the time-series of all output data, e(t) is the vector of noise signals that cannot be reduced during the system identification process, x(t) is the vector of state variables whose dimension is dependent on the number of order (n) of the model, and the entries of matrices, A, B, C, D, and K, are the parameters of the model to be identified to fit the given training data during the identification process and whose dimensions are dependent on the number of order of the model as well as the number of input/output data. For the case of the spindle in this article, the input data is only its time-series of rotational speed (single input), while the output data are the time-series of temperature variations of its seven characteristic points (7 outputs) shown in Figure 3c. Therefore, u(t) is a 1-dimensional vector whose component is the time-series of rotational speed, y(t) is a 7-dimensional vector whose components are the time-series of the temperature variations of the 7 characteristic points, and e(t) is a 7-dimensional vector. There are several commands, such as n4sid, ssest, and ssregest, in the System ID Toolbox of MATLAB for system identification [24]. As the first one, n4sid, is the fastest compared with the others, the authors use it throughout this article. The training data, i.e., the fold 1 of Table 4, was put into the command n4sid to conduct the system identification process. The command n4sid allows one to set a range of the model order (n) and determine the best model order within this range by the Hankel singular value. The Hankel singular values measure the contribution of each state to the input/output behavior, namely the states with small Hankel singular values can be discarded to simplify the model.
As a result, n4sid determines the best model order to be 8, as shown in Figure 14. Therefore, with the 8-order state space model, the state space equation of the spindle is expressed as . . .
After determination of the best model order and the parameters identification of the corresponding model, the predictive temperature and goodness of fit of the model can be evaluated by invoking the command "compare," which will be discussed in the following section, Results and Discussion. ( ( ) ( ) , After determination of the best model order and the parameters identification of the corresponding model, the predictive temperature and goodness of fit of the model can be evaluated by invoking the command "compare," which will be discussed in the following section, Results and Discussion. Instead of using the variation in Hankel singular value (as used in "n4sid" of MATLAB), the authors propose another way to determine the best model order by means of k-fold cross-validation, which is commonly used in machine learning [25]. For simplicity, the authors adopt 2-fold cross-validation with varying model order to determine what order is the best one that leads to the highest validation accuracy (goodness of fit). The criterion of goodness of fit of each output is defined by the normalized root mean square error (NRMSE),  Instead of using the variation in Hankel singular value (as used in "n4sid" of MATLAB), the authors propose another way to determine the best model order by means of k-fold cross-validation, which is commonly used in machine learning [25]. For simplicity, the authors adopt 2-fold cross-validation with varying model order to determine what order is the best one that leads to the highest validation accuracy (goodness of fit). The criterion of goodness of fit of each output is defined by the normalized root mean square error (NRMSE), where X obs (:,i) and X(:,i) are respectively the observed and predicted time-series of i-th output and mean( ) is the mean value of a time-series of data. Then the performance of SID is evaluated by the mean value of NRMSE i , that is where N is the number of outputs and N = 7 for this article because there are 7 characteristic temperature points in/on the spindle. Let If the given data D f old rpm is the same as the training data of M f old rpm,n , then the accuracy is called the training accuracy (TA) otherwise its called the validation accuracy (VA). For a given order n, one can evaluate the overall training accuracy (TA) by Similarly, the overall validation accuracy (VA) is evaluated by Based on the concept of cross-validation, which is commonly used in machine learning for classification and regression [26], it is necessary to vary the model complexity, namely the model order n, and search for a n* that maximizes the validation accuracy, Using cross-validation to determine the optimal complexity of a model can effectively strike a balance between bias and variance, such that both under-fitting and over-fitting can be avoided. Table 5 lists all the results of 2-fold cross-validation for the model at given rotational speeds as mentioned above, which shows that the model order of 25 results in the best MNRMSE in validation accuracy. As an example, Figure 15 shows the goodness of fit (MNRMSE) of the model, which is trained by the data D 1 5000 and tested with the data D 2 5000 .
Based on the concept of cross-validation, which is commonly used in machine learning for classification and regression [26], it is necessary to vary the model complexity, namely the model order n, and search for a n* that maximizes the validation accuracy, * arg max ( , ).
Using cross-validation to determine the optimal complexity of a model can effectively strike a balance between bias and variance, such that both under-fitting and over-fitting can be avoided. Table 5 lists all the results of 2-fold cross-validation for the model at given rotational speeds as mentioned above, which shows that the model order of 25 results in the best MNRMSE in validation accuracy. As an example, Figure 15 shows the goodness of fit (MNRMSE) of the model, which is trained by the data 1 5000 D and tested with the data 2 5000 D .

Speed-Dependence of the Thermal-Feature Model of the Spindle
As the thermal-feature model of the spindle is identified by the correlation between the input data, i.e., the time-series of spindle speed, and the output data, i.e., the time-series of spindle temperature at the input spindle speed, then, for different input spindle speeds, the identified parameters of thermal-feature models are slightly different. Hence it is worthwhile investigating the influence of the spindle speed on the variance of the parameters of the thermal-feature model. Here the authors use the data in the fold 1 of Table 4 as the training data to identify the model and those in the fold 2 as the test data to test the validation accuracies of each model. The results are displayed in both a matrix form and a bar chart as shown in Figures 16 and 17 . Each row entry of the matrix forms are also displayed in bar charts, Figures 16b and   17b, corresponding to a specific spindle speed. If the model is highly dependent on the spindle speed, then the VAs of the diagonal entries in the matrix forms should be much higher than the other entries. However, as shown in the bar charts, this is not true. The VAs of all models are very close, for 8-and 25-order models, they are about 98.5% and 99%, respectively. Therefore, the influence of the spindle speed on the variance of the model is approximately negligible. This means that one may use the thermal-feature model identified at a certain spindle speed to predict the spindle temperature at other speeds. It is rational to assume the thermal-feature model of the spindle is linear and time-invariant.

Speed-Dependence of the Thermal-Feature Model of the Spindle
As the thermal-feature model of the spindle is identified by the correlation between the input data, i.e., the time-series of spindle speed, and the output data, i.e., the time-series of spindle temperature at the input spindle speed, then, for different input spindle speeds, the identified parameters of thermal-feature models are slightly different. Hence it is worthwhile investigating the influence of the spindle speed on the variance of the parameters of the thermal-feature model. Here the authors use the data in the fold 1 of Table 4 as the training data to identify the model and those in the fold 2 as the test data to test the validation accuracies of each model. The results are displayed in both a matrix form and a bar chart as shown in Figures 16 and 17, wherein the former is identified by an 8-order model and the latter by a 25-order model. In the matrix forms, Figures 16a and 17a, the entry (i, j) represents the percentage VA of the model using the training data D 1 rpm i and test data D 2 rpm j . Each row entry of the matrix forms are also displayed in bar charts, Figures 16b and 17b, corresponding to a specific spindle speed. If the model is highly dependent on the spindle speed, then the VAs of the diagonal entries in the matrix forms should be much higher than the other entries. However, as shown in the bar charts, this is not true. The VAs of all models are very close, for 8-and 25-order models, they are about 98.5% and 99%, respectively. Therefore, the influence of the spindle speed on the variance of the model is approximately negligible. This means that one may use the thermal-feature model identified at a certain spindle speed to predict the spindle temperature at other speeds. It is rational to assume the thermal-feature model of the spindle is linear and time-invariant. Appl. Sci. 2018, 8,

Using the Thermal-Feature Model to Predict the Temperature Variation of the Spindle at Various Speeds
The previous section concluded that the influence of spindle speed on the variance of the thermal-feature model is approximately negligible, namely the parameters of the thermal-feature model are constant. This section tries to use the thermal-feature model  Table 4 to predict the temperature variation of the spindle rotating at various speeds. The spindle was set to rotate at a certain constant speed, meanwhile its temperature variation was recorded continuously. The spindle was turned down to stationary if its temperature variation was less than 1 ℃ for 3 consecutive hours. The temperature was recorded continuously until the spindle naturally cooled to ambient temperature. Figure 18 99

Using the Thermal-Feature Model to Predict the Temperature Variation of the Spindle at Various Speeds
The previous section concluded that the influence of spindle speed on the variance of the thermal-feature model is approximately negligible, namely the parameters of the thermal-feature model are constant. This section tries to use the thermal-feature model  Table 4 to predict the temperature variation of the spindle rotating at various speeds. The spindle was set to rotate at a certain constant speed, meanwhile its temperature variation was recorded continuously. The spindle was turned down to stationary if its temperature variation was less than 1 ℃ for 3 consecutive hours. The temperature was recorded continuously until the spindle naturally cooled to ambient temperature. Figure 18 99

Using the Thermal-Feature Model to Predict the Temperature Variation of the Spindle at Various Speeds
The previous section concluded that the influence of spindle speed on the variance of the thermal-feature model is approximately negligible, namely the parameters of the thermal-feature model are constant. This section tries to use the thermal-feature model M 1 6000,8 , namely the 8-order model identified from the training data D 1 6000 listed in Table 4 to predict the temperature variation of the spindle rotating at various speeds. The spindle was set to rotate at a certain constant speed, meanwhile its temperature variation was recorded continuously. The spindle was turned down to stationary if its temperature variation was less than 1 • C for 3 consecutive hours. The temperature was recorded continuously until the spindle naturally cooled to ambient temperature. Figure 18 shows the test results of the rotational speeds of 4500, 5300 and 7800 rpm, respectively, where the zoom-in window clearly illustrates that the prediction is very close to the measured data, the goodness of fit (MNRMSE) is over 98.94% for all locations. These results confirm that the thermal-feature model of the spindle can be represented by a linearly time-invariant state-space model with constant parameters.
shows the test results of the rotational speeds of 4500, 5300 and 7800 rpm, respectively, where the zoom-in window clearly illustrates that the prediction is very close to the measured data, the goodness of fit (MNRMSE) is over 98.94% for all locations. These results confirm that the thermal-feature model of the spindle can be represented by a linearly time-invariant state-space model with constant parameters.

Predicting Internal Temperature of the Spindle from Its Surface Temperature
In fact, the temperature inside the spindle is undetectable because there is no space ins spindle to embed the temperature sensor. In addition, the spindle manufacturer or user allows drilling on the spindle to place the temperature sensor inside because these hole reduce the mechanical stiffness of the spindle. However, the internal temperature may prov important index for prevention and maintenance of the spindle. Sections 3.1 and 3.2 conclu the speed-dependence of the thermal-feature model of the spindle is approximately negligib it can be represented by a linearly time-invariant state-space model with constant para Therefore, this section tries to use the temperature on the surface of the outer housing of the which is detectable directly, as the input data to identify the thermal-feature model of the s and then, by the resulting thermal-feature model, to figure out the undetectable temperature the spindle from its surface temperature. We set the spindle to rotate at a constant speed rpm, meanwhile we recorded its temperature variation continuously. The spindle was turned to stationary if its temperature variation was less than 1 ℃ for 3 consecutive hour temperature was recorded continuously until the spindle naturally cooled to ambient tempe The time-series of temperature of seven points inside/outside the spindle, as shown in Figure  recorded and serve as input or output data for the identification of the thermal-feature described in Table 6. According to the number of input (I) and output (O) data, four thermalmodels of the spindle, namely 1I7O, 3I5O, 2I5O, and 1I5O, are respectively identified by the l time-invariant 8-order state-space model. Then the resulting thermal-feature models were u predict the temperature of the spindle at various speeds, say 4500, 5300, and 7800 rpm.  Figures 19-21 show the results predic the models of 3I5O, 2I5O, and 1I5O, respectively. The prediction results agree very well w measured data; as shown in Figure 22, whose accuracies are over 99%. This illustrates the fea of predicting the internal temperature of the spindle from its surface temperature. Temperature of rear bearing B 3.
Temperature of front bearing C 4.
Temperature of front bearing D 5.
Temperature of inner housing 6.
Temperature of front outer housing

Predicting Internal Temperature of the Spindle from Its Surface Temperature
In fact, the temperature inside the spindle is undetectable because there is no space inside the spindle to embed the temperature sensor. In addition, the spindle manufacturer or user never allows drilling on the spindle to place the temperature sensor inside because these holes may reduce the mechanical stiffness of the spindle. However, the internal temperature may provide an important index for prevention and maintenance of the spindle. Sections 3.1 and 3.2 conclude that the speed-dependence of the thermal-feature model of the spindle is approximately negligible and it can be represented by a linearly time-invariant state-space model with constant parameters. Therefore, this section tries to use the temperature on the surface of the outer housing of the spindle, which is detectable directly, as the input data to identify the thermal-feature model of the spindle, and then, by the resulting thermal-feature model, to figure out the undetectable temperature inside the spindle from its surface temperature. We set the spindle to rotate at a constant speed of 6000 rpm, meanwhile we recorded its temperature variation continuously. The spindle was turned down to stationary if its temperature variation was less than 1 • C for 3 consecutive hours. The temperature was recorded continuously until the spindle naturally cooled to ambient temperature. The time-series of temperature of seven points inside/outside the spindle, as shown in Figure 3c, are recorded and serve as input or output data for the identification of the thermal-feature model, described in Table 6. According to the number of input (I) and output (O) data, four thermal-feature models of the spindle, namely 1I7O, 3I5O, 2I5O, and 1I5O, are respectively identified by the linearly time-invariant 8-order state-space model. Then the resulting thermal-feature models were used to predict the temperature of the spindle at various speeds, say 4500, 5300, and 7800 rpm. Figure 16 shows the results predicted by the 1I7O model, while Figures 19-21 show the results predicted by the models of 3I5O, 2I5O, and 1I5O, respectively. The prediction results agree very well with the measured data; as shown in Figure 22, whose accuracies are over 99%. This illustrates the feasibility of predicting the internal temperature of the spindle from its surface temperature. (e) Figure 19. Using the 3I5O thermal-feature model ( Table 6) Table 6.

Conclusions
This article proposes a methodology for the thermal-feature system identification of a mac tool spindle. A self-made temperature sensing and wireless transmission module was develope acquire the training data required for system identification, with an accuracy of 0.05℃ at 25 measurement range −50 to 300 ℃, 5 support channels, sampling rate 4.17 Hz, power consump 175 mW, and protection level IP65. The thermal-feature of the spindle was modeled by a line time-invariant state-space model. The parameters of the state-space model are identified by command "n4sid" provided by the System ID Toolbox of MATLAB. The best model-order doubly checked by the Hankel singular value (used in "n4sid" of MATLAB) and k-fold c validation (common in machine learning), wherein the accuracy of an 8-order model is about 98   Table 6.

Conclusions
This article proposes a methodology for the thermal-feature system identification of a machine tool spindle. A self-made temperature sensing and wireless transmission module was developed to acquire the training data required for system identification, with an accuracy of 0.05℃ at 25 ℃, measurement range −50 to 300 ℃, 5 support channels, sampling rate 4.17 Hz, power consumption 175 mW, and protection level IP65. The thermal-feature of the spindle was modeled by a linearly time-invariant state-space model. The parameters of the state-space model are identified by the command "n4sid" provided by the System ID Toolbox of MATLAB. The best model-order was doubly checked by the Hankel singular value (used in "n4sid" of MATLAB) and k-fold cross validation (common in machine learning), wherein the accuracy of an 8-order model is about 98.5%, Figure 22. The accuracies of temperature prediction by the thermal-feature models listed in Table 6.

Conclusions
This article proposes a methodology for the thermal-feature system identification of a machine tool spindle. A self-made temperature sensing and wireless transmission module was developed to acquire the training data required for system identification, with an accuracy of ±0.05 • C at 25 • C, measurement range −50 to 300 • C, 5 support channels, sampling rate 4.17 Hz, power consumption 175 mW, and protection level IP65. The thermal-feature of the spindle was modeled by a linearly time-invariant state-space model. The parameters of the state-space model are identified by the command "n4sid" provided by the System ID Toolbox of MATLAB. The best model-order was doubly checked by the Hankel singular value (used in "n4sid" of MATLAB) and k-fold cross validation (common in machine learning), wherein the accuracy of an 8-order model is about 98.5%, while that of a 25-order model is about 99%. The 25-order model is the optimal model for the present case. However, for the calculation efficiency, the 8-order model is a good model for the present case. The proposed method to determine the optimal complexity of a model, namely model order, can effectively strike a balance between bias and variance, such that both under-fitting and over-fitting can be avoided. Section 3.1 illustrates that the influence of spindle speed on the variance of the thermal-feature model of the spindle is approximately negligible. It is rational to assume the thermal-feature model of the spindle is linearly time-invariant. Section 3.2 illustrates that the thermal-feature model of the spindle can be represented by a linearly time-invariant state-space model with constant parameters. This means that one may use the thermal-feature model identified at a certain spindle speed to predict the spindle temperature at other speeds. Based on the assumption that the spindle speed has a negligible effect on the variance of the thermal-feature model, Section 3.3 uses the temperature on the surface of the outer housing of the spindle as input data and its internal temperature as output data to identify its thermal-feature model. Section 3.3 illustrates the feasibility of predicting the internal temperature (undetectable directly) of the spindle from its surface temperature (detectable directly). The internal temperature may provide an important index for prevention and maintenance of the spindle. It should be mentioned here that the research objective of this article is an externally driven spindle. However, for a motorized spindle, its thermal-feature model might be highly dependent on the spindle speed, which is worthy of studying in the near future. The applicability of the present methodology is not limited to externally-driven spindles, although it is validated by the run-in test of our available externally-driven spindle. Based on the results of the run-in test, it is indeed accurate to simulate the externally-driven spindle by a state-space model with constant coefficients. However, based on the literature, the state-space model of motorized spindles ought to use non-constant coefficients. For example, Lin and Tu [27] had mentioned that high-speed rotation can cause substantial changes in the dynamic and thermal behaviors of motorized spindles, and Liu and Zhang [28] concluded that the electromagnetic power loss of the built-in motor with motor slip at high-speed rotation will change the thermal-mechanical properties of motorized spindles.