Patent Publication Number: US-6341245-B1

Title: Machine tool displacement correcting apparatus

Description:
TECHNICAL FIELD 
     The present invention relates to a displacement compensation apparatus for machine tools to compensate displacement resulting from centrifugal force and heat generation caused by spindle rotation in a machine tool. 
     BACKGROUND ART 
     Conventionally, an apparatus to compensate heat displacement of a spindle based on temperature changes on the spindle is known as a heat displacement compensation apparatus to compensate heat displacement resulting from heat generation caused by rotation of the spindle in the machine tool. 
     In the above displacement compensation apparatus, compensation a correlation between the displacement of the spindle and a difference in temperature between places where a heat displacement is most likely to occur and where change in temperature is relatively small is examined in advance, for instance, a difference in temperature between near the spindle and the bed. A compensation value corresponding to the difference in temperature is calculated from the correlation, and the heat displacement is compensated based on the compensation value. But the conventional compensation method of heat displacement has a disadvantage that a difference is often found between an actual heat displacement of a spindle and a compensation value derived from the compensation calculation based on the difference in temperature, when the spindle is rotated at a random rotation frequency. 
     It is clear from FIG. 8 that a large deviation is found between a measured value and a compensation value obtained from the conventional compensation calculation, when the displacement of the spindle is measured after the spindle is rotated at rotation frequencies of, for instance, 1000 rpm, 3000 rpm, 0 (suspension), 2000 rpm, 3000 rpm, 1000 rpm for each 40 minutes and after leaving in suspension for 120 minutes. 
     When a spindle rotates, the spindle is expanded by centrifugal force. There is a disadvantage that an expansion of a spindle by centrifugal force makes the spindle shorten affected by the Poisson&#39;s effect, so that the position of the tool retreats to that extent from the product in process. 
     In the conventional heat displacement compensation method, as the compensation value is determined by the difference in temperature, it is impossible to compensate instant displacement caused by a change in the rotation frequency of the spindle. 
     The object of the present invention is to provide a displacement compensation apparatus for machine tools to precisely compensate displacement resulting from centrifugal force or heat generation caused by rotating spindle rotation even when the spindle is randomly rotated. 
     DISCLOSURE OF THE INVENTION 
     A displacement compensation apparatus for machine tools according to the present invention is to compensate a displacement resulting from centrifugal force and heat generation caused by rotation of a spindle in a machine tool. The apparatus is provided with; a historical data table storing a history of function values of rotation frequency of the spindle in each time period from present to past; a corrective coefficient table in which corrective coefficients are stored in advance relating the coefficients to the corresponding location in the historical data table and a compensation process means. The compensation process means reads out each location and each corrective coefficient which corresponds to the location from the corrective coefficient table and concurrently reads out a function value for the rotation frequency of the spindle at the corresponding location from the historical data table, calculates a compensation value from the function value at each location and the corresponding corrective coefficient, and compensates a commanded position of each shaft in the machine tool based on the compensation values. 
     The corrective coefficient is determined by collecting and storing measured data on displacement of the spindle, shaft rotation time data and rotation frequency data when the spindle is randomly rotated with changing the shaft rotation frequency and by processing a regression calculation based on these data. 
     The function value of the rotation frequency of the spindle is a value proportional to the square of the spindle rotation frequency, or a value based on the proportional value to the square of the spindle rotation frequency, with compensation by other value. 
     For instance, the function value is obtained from dividing a square of the rotation frequency of the spindle by 16384 and subtracting one tenth of the rotation frequency of the spindle. The number thus derived is used as the function value. The equation is expressed as follows. 
      Functional equation=(rotation frequency×rotation frequency/16384) −0.1×rotation frequency  (1) 
     In another example, a function value is calculated as a square of the rotation frequency of the spindle is divided by 16384 and further divided by one plus a number of three times of inverse of the natural logarithm of the rotation frequency divided by 25000. The functional equation is expressed as follows. 
     
       
         Functional equation=(rotation frequency×rotation frequency/16384/(1.0+3.0/exp(rotation frequency/25000))  (2) 
       
     
     In general, as centrifugal force is proportional to a square of a rotation frequency of a spindle, a displacement which is proportional to centrifugal force is also proportional to a square of the rotation frequency of the spindle. And as long as a force acting on the spindle and a friction coefficient stay constant, the friction torque is in proportion to the rotation frequency of the spindle. As friction is a product of friction torque and rotation frequency, it is proportional to a square of the rotation frequency of the spindle. Heat generation is a result of consumed energy by friction, it is also proportional to a square of the rotation frequency of the spindle. Hence, when a value which is proportional to a square of a spindle rotation frequency is used as a function value of the spindle rotation frequency, the accuracy of compensation can be improved. 
     In actuality, as a force acting on the spindle tends to increase with the increase of the rotation frequency, beat generation changes more sharply than a proportion to a square of the rotation frequency of the spindle. Consequently, as a function value for the rotation frequency of the spindle, a number proportional to a square of the spindle rotation frequency as a main component and compensated by another number is used. Then the accuracy of compensation is much improved. 
     In compensating displacement, each stored location and corrective coefficient are read out from the corrective coefficient, concurrently the function value of the spindle rotation frequency in each location is read out from the historical data table. And the compensation value is calculated from the function value in each location and the corresponding corrective coefficient. 
     Now, each location which is read out from the corrective coefficient table is taken as k 00 , k 01 , k 02 , . . . k m , and corresponding each corrective coefficient as b 00 , b 01 , b 02 , . . . b m . And the function value of the spindle rotation frequency at each location read out from the historical data table is taken as R(k 00 ), R(k 01 ), R(k 02 ), . . . R(k m ), then the displacement amounts δ is calculated from the following equation. 
     
       
         δ=R(k 00 )×b 00 +R(k 01 )×b 01 +R(k 02 )×b 02  . . . +R(k m )×b m   (3) 
       
     
     The compensation value ε is obtained by changing the signs of the displacement amount δ. 
     More specifically, as the corrective coefficient value obtained from the equation (3) is very small, 4194304 times of the corrective coefficient value is stored and the displacement amount δ is obtained from dividing the calculated value by 4194304. 
     For instance, when compensating three axes X, Y and Z, the location read out from the corrective coefficient table are expressed as k 00 , k 01 , k 02  . . . k m  the corrective coefficient in X axis corresponding to the above location are b 00x , b 01x , b 02x , . . . b mx , the corrective coefficient in Y axis corresponding to the location are b 00y , b 01y , b 02y , . . . b my , and the corrective coefficient in Z axis are b 00z , b 01z , b 02z , . . . b mz . And the function values of rotation frequency of the spindle at each location read out from the historical data table are expressed as R(k 00 ), R(k 01 ), R(k 02 ), . . . R(k m ). The displacement amounts for each axis δ X , δ y , δ z  can be obtained by the following equations. δ X =(R(k 00 )×b 00x +R(k 01 )×b 01x +R(k 02 )×b 02x  . . . R(k m )×b mx )/4194304 
     
       
         δ y =(R(k 00 )×b 00y +R(k 01 )×b 01y +R(k 02 )×b 02y  . . . R(k m )×b my )/4194304 
       
     
     
       
         δ z =(R(k 00 )×b 00y +R(k 01 )×b 01y +R(k 02 )×b 02y  . . . R(k m )×b my )/4194304  (4) 
       
     
     A compensation value is a value of which absolute number is the same as the displacement amount and which is shown in opposite signs. Therefore, compensation values for the axes (X, Y, and Z axis) ε x , ε y , ε z , are obtained from the following equation. 
     
       
         ε X =−δ x , ε y =−δ y , ε z =−δ z   (5) 
       
     
     Accordingly, the compensation value is calculated from the corrective coefficient considering the rotation time, rotation frequency, elapsed time of the spindle, and the function value of the spindle rotation frequency. As the command position of the machine tool is compensated based on the compensation value, the displacement resulting from centrifugal force caused by the rotation of the spindle or heat generation spindle can be compensated with high accuracy even when the spindle is randomly rotated. 
     That is, in a conventional compensation method, the compensation value is calculated based on a temperature difference at that time. When the spindle is randomly rotated, the compensation value for heat displacement of the spindle varies from small number to big number. And it is impossible to compensate the displacement resulting from changes in rotation frequency of the spindle. But the apparatus according to the present invention gives a compensation value having little difference from the true displacement of the spindle even when the spindle is rotating at a random speed and also compensates a displacement caused by centrifugal force. Thus the apparatus according to the present invention can achieve highly accurate compensation. 
     It is preferable that the function values are stored in the data table at relatively short time intervals while the spindle is rotating and at relatively long time intervals while the spindle is halted. For instance, the function values are stored at intervals from 20 seconds to one minute while the spindle is rotating and stored at intervals from one minute to 3 minutes while the spindle is halted. 
     In this way, a precise compensation of the displacement can be achieved with a single table while the spindle is rotated and also while the spindle is halted. Incidentally, displacement characteristics of the spindle in rotation is different from displacement characteristics of the spindle in suspension. If the function values are renewed and stored still at the same time intervals, the compensation value derived may contain some error. This may be considered as a disadvantage of the method. But the method according to the present invention ensures an accurate compensation of the displacement with one table for both in rotation and in suspension of the rotation of the spindle. 
     It is also preferable that the locations in the data table are arranged in an order reverse to elapsed time. 
     Most of presently available micro processing units for the computer (it will be abbreviated to MPU hereafter) are provided with a batch transfer function of its memory contents. It takes shorter time for carrying out processing with the batch transfer function provided in the MPU than to make and use programs of reading out individual data from a source address and writing in a destination address. But in this function, the processing is sequentially carried out from an address of a lower order to higher order. Accordingly, when the range of address in source address is overlapped with the range of address in destination address and if the destination address is in higher order than the order of the source address, the overlapped range area can not be processed correctly. In the present invention, as the locations in the data table are arranged in an order reverse to elapsed time, the above described problem can be avoided. 
     Moving averages(“ryu-do heikin” in Japanese) of the function values of predetermined plural locations in the past time period may be stored in the data table at the time of data renewal, and one corrective coefficient value in correspondence to the time period above described may be stored in the corrective coefficient table. 
     Then, the size (the number of lines) of the corrective coefficient table can be reduced and the speed for the compensation calculation can be improved. 
     For the value stored in the historical data table, the rotation frequency of the spindle itself can be used instead of the function value of the spindle rotation frequency. 
     In such cases, the compensation procedure is as follows. A location in the historical data table and the corresponding corrective coefficient are read out from the corrective coefficient table and at the same time the rotation frequency of the spindle in the location is read out from the historical data table. The function values are calculated from the rotation frequencies of the spindle in the locations as necessary. And with the corresponding corrective coefficients the calculations of the corrective values are carried out. Then, the commanded position for each shaft of the machine tool is compensated based on the corrective values. 
    
    
     BRIEF DESCRIPTION OF DRAWINGS 
     FIG. 1 is a block diagram of the processing system showing an embodiment of the displacement compensation apparatus relating to the present invention; 
     FIG. 2 is a block diagram of a numeric control system for the above embodiment; 
     FIG. 3 is a corrective coefficients table for the above embodiment; 
     FIG. 4 is a block diagram showing the method of computations for the corrective coefficient in the above embodiment; 
     FIG. 5 is a historical data table for the above embodiment; 
     FIG. 6 is a table after renewal of the historical data table for the above embodiment; 
     FIG. 7 is a graph showing the relations among changes of the rotation frequency of the spindle, the displacement amount obtained through the compensation computations (compensation value) and measured values of the displacement in the above embodiment; and 
     FIG. 8 is a graph showing the relations between the compensation values obtained through the conventional compensation method and the measured values. 
    
    
     BEST MODE FOR CARRYING OUT THE INVENTION 
     An embodiment of the present invention will now be explained in detail with reference to the drawings. 
     FIG. 1 shows a processing system of the displacement compensation apparatus of the present invention. The processing system consists of a horizontal boring machine  1  as a machine tool and a numeric control unit  2  which controls the horizontal boring machine and also carries out the operation of compensation. 
     The horizontal boring machine is constituted by a bed  11 , a table provided on the bed  11  in a movable manner along the X-axis (directions intersecting at right angles with the paper surface of FIG.  1 ), a column  13  mounted and firmly fixed on the bed  11 , a spindle head  14  provided to the column  13  in a movable manner up and down along the Y-axis (the top and the bottom directions on FIG.  1 ), a ram  15  provided to the spindle head  14  in a reciprocally movable manner along the Z-axis directions (right and left directions on FIG. 1) and a spindle  16  accommodated in the ram in a rotatable manner. Incidentally, the designation  17  is a motor driving the spindle  16 . 
     The numerical control unit  2  is provided with a CPU  21 , a memory  22  which stores various programs and data, a keyboard  23 , a display  24 , a printer  25 , and I/O  26 ,  27 ,  28 ,  29  as shown in FIG.  2 . Each I/O  26 ,  27 ,  28  is connected to the actuating systems toward X, Y, Z axis directions of the horizontal boring machine,  18 X,  18 Y,  18 Z and the I/O  29  is connected with the motor  17 . The CPU  21  controls movements of the actuating systems toward X, Y, Z directions  18 X,  18 Y, and  18 Z and performs compensation processing. Thus, the CPU  22  also serves as a compensation processing means. 
     The memory  22  is provided with not only a memory area for storing various processing programs (not shown in figures), but also a coordinate value table  22 A which stores coordinate values of each axis (X, Y, Z), a compensation value table  22 B which stores compensation values ε x , ε y , ε z  of each axis (X, Y, Z), a data table  22 D which stores a history of function values of the rotation frequency of the spindle on the machine tool in each time period from present to past and a corrective coefficient table  22 C which stores, in advance, corrective coefficients of each axes (X, Y, Z) relating the coefficients to the addresses of the data table  22 D. 
     The registration of the function value and the renewal of the history in the data table  22 D are carried out at intervals of 20 seconds to 3 minutes. More specifically, the registration and renewal are carried out at intervals of 20 seconds to 1 minute when the spindle is rotating, and from 1 minute to 3 minutes when the rotation of the spindle is suspended. Thus, the history of the function values stored in the data table  22 D is memorized at intervals of a relatively short period of time while the spindle is rotating and at intervals of a relatively long period of time while the rotation of the spindle is suspended. 
     Here, FIG.  5  and FIG. 6 show the historical data table  22 D. For convenience of explanation, assuming that the intervals of the registration and renewal for the whole period is one minute, the elapsed time is recorded in the elapsed-time column on FIG.  5  and FIG.  6 . As a matter of fact, the elapsed time varies according to rotation or suspension of the spindle, the actual elapsed time is not always the same as shown in the figures. 
     Then, the randomly changeable elapsed time on the historical data table  22 D is each defined as X1, X2, X3, . . . X512. 
     As shown in FIG. 3, in relation to the addresses from k 512  to k 00  (hereafter, they may be sometimes written as “from 512 to 0”) in the data table  22 D, the corrective coefficient of each axis in each time period for every cycle of renewal from present to X4 ago, 
     (b 00x , b 00y , b 00z ) (b 01x , b 01y , b 01z ) . . . (b 04x , b 04y , b 04z ) 
     are stored in the corrective coefficient table  22 C. 
     From X4 through X12 ago, the corrective coefficient of each axis in each time period at intervals of X2, 
     (b 05x , b 05y , b 05z ) (b 06x , b 06y , b 06z ) . . . (b 08x , b 08y , b 08z ) 
     are stored. 
     From X12 through X24 ago, the corrective coefficient of each axis in each time period at intervals of X4, 
     (b 09x , b 09y , b 09z ) (b 10x , b 10y , b 10z ) (b 11x , b 11y , b 11z ) 
     are stored. 
     From X24 through X64 ago, the corrective coefficient of each axis in each time period at intervals of X8, 
     (b 12x , b 12y , b 12z ) (b 13x , b 13y , b 13z ) . . . (b 16x , b 16y , b 16z ) 
     are stored. 
     From X64 through X128 ago, the corrective coefficient of each axis in each time period at intervals of X16, 
     (b 17x , b 17y , b 17z ) (b 18x , b 18y , b 18z ) . . . (b 20x , b 20y , b 20z ) 
     are stored. 
     From X128 through X512 ago, the corrective coefficient of each axis in each time period at intervals of X32, 
     (b 21x , b 21y , b 21z ) (b 22x , b 22y , b 22z ) . . . (b 32x , b 32y , b 32z ) 
     are stored. 
     In this point, the way to obtain the corrective coefficient stored in the corrective coefficient table  22 C is explained. 
     First, measured displacement values of top end portion of the spindle  16 , rotation time and rotation frequencies of the spindle  16  are collected as data through experiments. In this point, in order to make clear the relation between heat generation and heat emission of the machine, the rotation frequency of the spindle  16  is randomly changed (increasing or decreasing the rotation frequency). 
     For instance, the spindle  16  is rotated at a speed of 5000 rpm for 30 minutes, 10000 rpm for an hour, 8000 rpm for 30 minutes, 15000 rpm for an hour, 5000 rpm for30 minutes and 15000 rpm for 30 minutes before it is stopped. The observation begins before the rotation is started and continues even after the rotation is halted until the displacement is substantially recovered without interruption. 
     During the observation, the displacements of the spindle  16  in the X, Y, Z directions are measured at predetermined time intervals (for instance, every 30 seconds) and stored with the spindle rotation frequencies and the time of the measurement. 
     Then, based on the data, each corrective coefficient is determined after a multiple regression analysis. 
     Incidentally, the calculation of the multiple regression analysis is carried out through the procedures shown in FIG.  4 . 
     Initially, the function values of the spindle rotation frequency are set as independent variables. Here, as the amounts of change for one time period in latest several times is large, a function value is set for each time period during the latest several periods as an independent variable. After the elapsed time exceeds X5, as the amount of change for each time period becomes a little slow, an average of the function values of the spindle rotation frequency between X2 is taken as an independent variable. When the elapsed time exceeds X12, the change for each time becomes slower, then an average of the function values of the spindle rotation frequency between X4is set as independent variable. Thus the time period for taking an average is extended according to the change in displacement amount that gets slow with lapse of time. At last, an average between X32 from X481 to X512 is taken as the thirty-third independent variable. 
     The reason for using an average value of the function values as independent variable is; 1. to reduce the number of the independent variables and 2. to prevent confounding between the independent variables. 
     In other words, when the number of the independent variables is large, the simultaneous equations with corresponding number of the independent variables are required to be solved. Under present processing circumstances of using a floating decimal, so-called double precision, the number of independent variables over seventy makes the calculation suddenly difficult. The number of the independent variables of thirty-three is about half of this number seventy and makes no problem. 
     Incidentally, the word “confounding” is used for a situation that a degree of correlation between independent variables becomes stronger than the correlation between the independent variable and the dependent variable. When the confounding occurs, a regression coefficient varies on a large scale with slight change in the data value. In the extreme, solutions to the simultaneous equations may not be obtained. 
     After setting measured values of the spindle displacement as dependent variables, sum of products calculation is carried out. That is, for each independent variable and dependent variable, the sum and average are respectively calculated, and the sum of products are obtained in each relations. The procedure is repeated as many times as the number of the data. 
     Next, sum of products is calculated after data of the sum of products is compensated by data of the sum. Then, coefficients of the simultaneous equations are determined. And after establishing criterion values for convergence, the simultaneous equations are solved through Gauss&#39;s Law of Discharge(sweeping. “Gausu-no-hakidashi-ho” in Japanese) so that the regression coefficients (corrective coefficients) are obtained and printed. The procedure is repeated as many times as the number of the independent variables. And the regression coefficients (corrective coefficients) are printed. 
     As shown in FIG. 5, the history of the function values of the spindle rotation frequency for each time period from present to past in the address is stored from 512 to 0 of the historical data table  22 D. The calculation of the function value of the spindle rotation frequency is based on the equation (1). 
     Here, the method to store the function value of the spindle with renewing the data table  22 D one by one and at the same time calculating the moving averages(“ryu-do heikin” in Japanese) will be explained. FIG. 6 shows the historical data table  22 D after the renewal. 
     When power of the machine is turned on, as the spindle is not rotating before then, the function values for all lapse of time are initialized to zeros. Then, following calculations (1), (2), (3) are repeatedly performed every time. And when the power is turned off, the renewal of the data table  22 D can not be conducted during the suspension of the power. But when the power is turned on again, the number of times corresponding to the suspension are renewed. And if the time of the suspension is longer than the period of the data table  22 D, the function values for all lapse of time are initialized to zeros just as in the case of turning the power on. 
     (1) Calculation of the Sum 
     {circle around (1)} A function value in the address 507 and a function value in the location number 508 are summed up and the total is stored in the memory A 2 . 
     {circle around (2)} A function value in the address 497 and a function value in the location number 499 are summed up and the total is stored in the memory A 4 . 
     {circle around (3)} A function value in the address 481 and a function value in the location number 485 are summed up and the total is stored in the memory A 8 . 
     {circle around (4)} A function value in the address 433 and a function value in the location number 441 are summed up and the total is stored in the memory A 16 . 
     {circle around (5)} A function value in the address 353 and a function value in the location number 369 are summed up and the total is stored in the memory A 32 . 
     (2) Shift 
     Five hundred and twelve function values beginning from the address 1 are shifted to the address 0. 
     (3) Calculation of the Average 
     {circle around (1)} The value in the memory A 2  is divided by 2 and the quotient is written in the address 506. 
     {circle around (2)} The value in the memory A 4  is divided by 2 and the quotient is written in the address 496. 
     {circle around (3)} The value in the memory A 8  is divided by 2 and the quotient is written in the address 480. 
     {circle around (4)} The value in the memory A 16  is divided by 2 and the quotient is written in the address 432. 
     {circle around (5)} The value in the memory A 32  is divided by 2 and the quotient is written in the address 352. 
     Incidentally, many available micro processing units for the computer (hereafter, it will be called MPU for short) are provided with a batch transfer function for its memory contents. It takes shorter time to carry out processing with the batch transfer function provided in the MPU than to make and use programs of reading out individual data from a source address and write in a destination address. But in this function, the processing is sequentially carried out from an address of a lower order to higher order. Accordingly when the range of address in source address is overlapped with the range of address in destination address and if the destination address is in a higher order than the order of the source address, the overlapped range area can not be processed correctly. To avoid this problem, the location numbers are allocated in a sequence reverse to the elapsed time in FIG.  5  and FIG.  6 . 
     The effect of the embodiment of the present invention will be explained next. 
     When power of the horizontal boring machine  1  and the numeric control unit  2  is turned on and the horizontal boring machine  1  starts operating through a command from the numeric control unit  2 , the coordinate values of each shaft driving system  18 X,  18 Y,  18 Z are renewed and stored in the coordinate value table  22 A and the function values of the rotation frequency of the spindle  16  are stored in the historical data table  22 D. 
     All the while, the compensation procedures are carried out at every predetermined interval. The CPU  21  reads out each address and the corrective coefficient from the corrective coefficient table  22 C with reference to the time of performing the compensation, and reads out the function values of each shaft rotation frequency in the address of the data table. After that, the CPU multiplies the function values of the spindle rotation frequency in the address and the corresponding corrective coefficients and obtains through calculation the compensation values based on the total sum of each product. 
     When compensating the three axes X, Y, Z, assuming that the function value of the spindle rotation frequency at the present time X0 related in the addresses from 512 to zero and the corresponding corrective coefficients of each axis are R 512 , b 00x , b 00y , b 00z , the function value of the spindle rotation frequency at the period before X1 and the corresponding corrective coefficients of each axis are R 511 , b 01X , b 01y , b 01Z , the function value of the spindle rotation frequency at the period before X2 and the corresponding corrective coefficients of each axis are R 510 , b 02x , b 02y , b 02z , . . . the function value of the spindle rotation frequency at the period before X5 to X6 and the corresponding corrective coefficients of each axis are R 506 , b 05X , b 05y , b 05 , . . . and the function value of the spindle rotation frequency as of the period before X481 to X512 and the corresponding corrective coefficients of each axis are R 0 , b 32X , b 32y , b 32Z , the amount of displacement for each axis, δ X , δ y, δ   Z  are calculated from the equation (4). 
     That is, the amount of displacement of δ X , δ Y , δ Z  is each calculated from following equations. 
     
       
         δ X =(R 512 ×b 00X +R 511 ×b 01X +R 510 ×b 02X  . . . +R 506 ×b 05x  . . . +R 0 ×b 32X /4194304) 
       
     
     
       
         δ Y =(R 512 ×b 00y +R 511 ×b 01y +R 510 ×b 02y  . . . +R 506 ×b 05y  . . . +R 0 ×b 32y /4194304) 
       
     
     
       
         δ Z =(R 512 ×b 00Z +R 511 ×b 01Z +R 510 ×b 02Z  . . . +R 506 ×b 05Z  . . . +R 0 ×b 32Z /4194304) 
       
     
     The compensation values ε X  (=−δ X ), ε y  (=−δ y ), ε Z  (=−δ Z ) are calculated by substituting these displacement amounts δ X , δ y , and δ Z  into the equation (5). 
     For instance, when R 512 , the function value of the spindle rotation frequency at the present time X0, is 5104, b 00X  the corresponding coefficient of X axis, is −2687, b 00y , the corresponding coefficient of Y axis, is −1080, and b 00z , the corresponding coefficient of Z axis, is 7450. 
     When R 511 , the function value of the spindle rotation frequency at the period before X1, is 5104, b 01X , the corresponding coefficient of X axis, is 2525, b 01y , the corresponding coefficient of Y axis, is 425, and b 01z , the corresponding coefficient of Z axis, is −6467. 
     When the function value of the spindle rotation frequency at the period before X2, R 510  is 5104, the corresponding coefficient of X axis, b 02x  is 1178, the corresponding coefficient of Y axis, b 02y  is −206, and the corresponding coefficient of Z axis, b 02z  is −3456. 
     When R 506 , the function value of the spindle rotation frequency at the period from X5 to X6, is 3065, b 05X , the corresponding coefficient of X axis, is 223, b 05y , the corresponding coefficient of Y axis, is 2000, and b05z, the corresponding coefficient of Z axis, is −3120. 
     When R 0 , the function value of the spindle rotation frequency at the period from X481 to X512, is 0, b 32X , the corresponding coefficient of X axis, is 10, b 32y , the corresponding coefficient of Y axis, is 100, and b32z, the corresponding coefficient of Z axis, is −600. If these situations are assumed, the displacement amounts of the axes δ X , δ y , δ Z  are calculated as follows. 
     
       
         δ X =(5104×(−2687)+5104×2525+5104×1178 . . . +3065×223 . . . +0×10)/4194304 =10185628/4194304=2.428[μm] 
       
     
     
       
         δ y =(5104×(−1080)+5104×425+5104×(−206) . . . +3065×2000 . . . +0×10)/4194304=23201458/4194304=5.532[μm] 
       
     
     
       
         δ z =(5104×7450+5104×(−6467)+5104×(−3456) . . . +3065×(−3120) . . . +0×(−600)/4194304 =65685732/4194304=−15.661[μm] 
       
     
     Accordingly the compensation values ε X , ε y , ε Z  are as follows. 
     ε X =−2.428[μm] 
     ε y =−5.532[μm] 
     ε z =15.661[μm] 
     In such a case that the corrective coefficient for the elapsed time corresponding to the time shown in FIG. 5 is not found in FIG. 3, the corresponding calculation is not given. For instance, the elapsed time corresponding to the elapsed time (from X12 to X13) for address 499 is not found in FIG. 3, the corresponding calculation is not given. But the value not used at that time will be used several times later due to the shift at every X1. 
     Thus, after the displacement amount in each axis δ X , δ y , δ Z  are calculated from the function values of the spindle rotation frequency R 512 , R 511 , R 510 , . . . R 506  . . . R 0  for each time period from present to past and each corrective coefficient corresponding to the above period of time 
     b 00X , b 01X , b 02X , . . . b 05X  . . . b 32x    
     b 00y , b 01y , b 02y , . . . b 05y  . . . b 32y    
     b 00Z , b 01Z , b 02Z , . . . b 05Z  . . . b 32Z    
     the compensation values for each axis ε X , ε y , ε Z  are obtained by reversing the signs of the displacement amount δ X ,δ y , δ Z . These compensation values ε X , ε y , ε Z  are renewed and memorized in the compensation table  22 B. 
     Then, the NC command position commanded after the renewal is adjusted by the amount of compensation values ε X , ε y , ε Z  in the compensation value table  22 B and the displacement amount in each axis is compensated by using the compensation value ε X , ε y , ε Z  memorized in the compensation value table  22 B. 
     In the embodiment of the present invention, a history of function values of the spindle in each time period from present to past are stored in the historical data table  22 D, and the corrective coefficients are stored in advance in the corrective coefficient table  22 C relating the coefficients to the address in the historical data table  22 D. In compensation procedure, the address and the corrective coefficient are read out from the corrective coefficient table and the function value of the spindle rotation frequency corresponding to each address is read out from the data table  22 D. The compensation value is calculated from the function value of the spindle rotation frequency at the corresponding address and corrective coefficient corresponding to this function number. Thus, as the NC command position is compensated based on this compensation value, even when the spindle  16  is randomly rotated, the displacement resulting from centrifugal force and heat generation caused by rotation of the spindle  16  can be compensated accurately. 
     For instance, when the spindle is rotated in a condition shown in FIG.  7 (A), the measured displacement values of three axis X, Y, Z and the displacement amounts of each axis obtained from the compensation calculation in the embodiment of the present invention δ X  (=−ε X ),δ y  (=−ε y ),δ z  (=−ε z ), are shown in FIG. (B), FIG. (C), and FIG. (D). As is clear from these figures, each displacement amounts δ X  (=−ε X ),δ y  (=−ε y ), and δ z  (=−ε z ) are very close to the measured values of displacement with little error, and it is apparent that the compensations of the NC command position based on the compensation values can achieve a high accuracy compensation, even when the spindle is randomly rotated. 
     The calculation of the compensation values ε X , ε Y , ε Z  are obtained from the following procedures. Assuming that the function values of the spindle in each time period from present to past are R 512 , R 511 , R 510 , . . . R 506 . . . R 0 , the corrective coefficients are b 00X , b 01X , b 02X , . . . b 32X , b 00y , b 01y , b 02y , . . . b 05y , b 32y , b 00Z , b 01Z , b 02Z , . . . b 05Z  . . . b 32Z , the displacement amounts δ X , δ Y , δ Z  are calculated from, 
     
       
         δ X =(R 512 ×b 00X +R 511 ×b 01X +R 510 ×b 02X  . . . +R 506 ×b 05X  . . . +R 0 ×b 32X /4194304 
       
     
     
       
         δ y =(R 512 ×b 00y +R 511 ×b 01y+R   510 ×b 02y  . . . +R 506 ×b 05y  . . . +R 0 ×b 32y /4194304 
       
     
     
       
         δ Z =(R 512 ×b 00Z +R 511 ×b 01Z +R 510 ×b 02Z  . . . +R 506 ×b 05Z  . . . +R 0 ×b 32Z /4194304 
       
     
     Then, the compensation values ε X , ε y , εZ Z  are obtained by changing the sign of the displacement amounts δ x , δ y , δ z . Thus the compensation values ε x , ε y , ε z  can be obtained by such a very simple calculation. 
     With regard to X axis, for instance, the corrective coefficients b 00X , b 01X , . . . b 04X  corresponding to each X1 time period from present before X4, corrective coefficients b 05X , b 06X , . . . b 08X , corresponding to every X2 time period from X4before X12, corrective coefficients b 09X , b 10X , b 11X  corresponding to every X4time period from X12 before X24, corrective coefficients b 12X , b 13X , . . . b 16X  corresponding to every X8 time period from X24 before X64, corrective coefficients b 17X , b 18X , . . . b 20X  corresponding to every X16 time period from X64 before X128, corrective coefficients b 21X , b 22X , . . . b 32X  corresponding to every X32 from X128 before X512 are respectively stored in the corrective coefficient table  22 C. In other words, as the period of time becomes farther back, the corrective coefficient is taken with longer intervals, thus the corrective coefficients can be stored with a small memory capacity. 
     And as a value proportional to a square of the main valve rotation frequency, or a compensated value based on a value proportional to a square of the main valve rotation frequency with a compensation by other value, is used as a function value of the spindle rotation frequency to store in the historical data table  22 D, a displacement resulting from centrifugal force and heat generation caused by the rotation of the spindle  16 . 
     In the above described embodiment, the corrective coefficients are stored in every X1 from present till X4, and the corrective coefficients are stored in every X2 from more than X4 till X12, in every X4 from more than X12 till X24, in every X8 from more than X 24 till 64, in every X16 from more than X64 till X128, and in every X32 from more than X128 till X512 are respectively stored. But the time period for the corrective coefficients is not limited within the period written above, and any optional period of time may be used for storing the corrective coefficient. 
     For instance, if a method of storing the corrective coefficients in every one minute from zero minute to 512 minutes ago, is adopted, the corrective coefficient table  22 C needs more than 15 times as large a memory capacity and also needs more than 15 times as large a multiplication and addition circuit to calculate the displacement amount and the compensation value. But the advantage is that the aforementioned calculation of the sum (1) and calculation of the average is unnecessary, and the program for the compensation is simplified. 
     In the above described embodiment, the compensation value is calculated using the function value of the spindle rotation frequency, but direct use of the spindle rotation frequency for the calculation also gives compensation values with improved accuracy than that of the conventional heat displacement compensation method. 
     By the heat displacement compensation apparatus of machine tools according to the present invention, displacement resulting from centrifugal force and heat generation caused by rotation of the spindle, can be compensated with high accuracy, even when the spindle is randomly rotated. 
     Industrial availability 
     The present invention is best suited for compensating the displacement resulting from rotation of the spindle of a machine tool such as a horizontal boring machine and a milling machine which uses rotation of the spindle for processing.