NURBS interpolator with pre-compensation based on discrete inverse transfer function for CNC high-precision machining

In the field of computerized numerical control (CNC) machining, high-speed and high-precision machining has been regarded as the key research by many scholars. In conventional methods, high-speed machining and high-precision machining are contradictory. It is inevitable to reduce the feedrate to improve the processing accuracy. In the paper, a pre-compensation based on discrete inverse transfer function (PDIT) theory is proposed. PDIT is able to improve machining contour accuracy without decreasing feedrate. The proposed PDIT theory is divided into three parts: NURBS interpolator, feedrate scheduling, and interpolator with pre-compensation. The NURBS interpolator has great advantage of directly interpolating the parametric curve. Therefore, the paper adopts the NURBS interpolator to accomplish interpolation. In the feedrate scheduling, S-type flexible acceleration and deceleration are used for path planning, and the maximum starting feedrate is obtained under the feedrate constraint. In the interpolator with pre-compensation, the NURBS interpolator is pre-compensated by PDIT. For the input, the response of the transfer function reaches a steady-state response within periods of time. Before steady-state response, the unsteady-state response exists in the transfer function. The unsteady-state response usually sustains dozens of interpolation periods, which inevitably lead to contour error in machining. Hence, the PDIT theory is employed to compensate the contour error caused by the unsteady-state response of transfer function to NURBS interpolator. The drive system is a transfer function, so the unsteady-state response of drive system can also cause machining errors before the steady-state response. In the paper, the NURBS interpolator is pre-compensated by PDIT theory before drive system to reduce contour errors and improve machining accuracy. Finally, the performance of the proposed PDIT is evaluated by simulations and experiments. The experimental results illustrate that PDIT theory obviously improve the machining accuracy.


Introduction
In the modern manufacturing industry, the demand for highspeed and high-precision machining is increasing, and highspeed and high-precision machining has always been the topic research by many scholars. Many scholars reduce the vibration of the machine tool by reducing the fluctuation of the feedrate to achieve the purpose of improving the machining accuracy. Dong and Wang proposed the target feedrate filter method to achieve smooth feedrate [1,2]. Liu Huan and other scholars proposed a linear programming mathematical model for the optimal feedrate in sensitive areas to obtain the optimal feedrate [3]. However, the disadvantage of this method is that the real-time performance is low. It needs the help of MATLAB linprog. Guo transformed the axial and tangential feedrate constraints into a convex optimization problem to find the optimal feedrate [4]. There are also many scholars to improve the machining accuracy by optimizing the interpolation algorithm. Scholars such as Wolfgang BOHM proposed parameter interpolation algorithms based on Lagrange [5], and Sang et al. proposed interpolation algorithms based on Taylor expansion [6][7][8][9]. Scholars Qin Hu proposed the optimized Taylor parameter interpolation algorithm to compensate the truncation error [10][11][12][13][14][15]. The above scholars mainly focus on the CNC system, and some scholars also focus on the drive system. Scholars who pay attention to drive system mainly study contour errors caused by drive system. Hu proposed an analytical contour compensation algorithm to improve machining accuracy [16,17]. Scholars Li proposed a method based on monocular vision to obtain and compensate contour errors [17]. Whether the feedrate smoothing algorithm and parameter interpolation algorithm in the CNC system, or the contour error compensation algorithm in the drive system, they are all subjectively divided into two parts. However, very little research literature combines the two. Therefore, the paper proposes the PDIT theory and makes full use of the relationship between the two to improve machining accuracy without losing feedrate. The PDIT theory studies the influence of transfer function of drive system on contour error, and directly compensates the contour error to the CNC commanded position.
In order to improve the performance of the NURBS interpolator, PDIT theory is proposed in the paper. PDIT is consist by NURBS interpolator, feedrate scheduling, and interpolator with pre-compensation. NURBS curve has great advantage of expressing free curve and is adopted in the paper. In order to obtain the maximum smooth feedrate and improve the efficiency of planning path, the feedrate scheduling based on S-type flexible acceleration and deceleration algorithm is adopted. This is introduced in Sect. 2. In Sect. 3, the principle of PDIT is introduced in detail. PDIT has ability to resolve the contour error caused by unsteadystate response. In the Sect. 4, simulations and experiments are carried out to evaluate the proposed NURBS interpolator. And conclusions are in Sect. 5.

Implementation of feedrate scheduling for NURBS interpolator
Non-uniform rational B-spline (NURBS) is defined by the node vector U, a set of control points P, and the weight of each control point. Assuming that the tool path is described by a NURBS curve of order p, the curve can be expressed as [19][20][21]: where P i is the control point, ω i is the corresponding weights of P i , (n + 1) is the number of control points, and p is the order of a NURBS curve. N i,p (u) is the pth-degree of B-spline basis functions defined on the non-uniform knot vector and m + 1 is the number of knots. The de Boor-Cox algorithm can be used to calculate the p-degree B-spline basic function N i,p (u) is [18]: NURBS interpolator has ability to directly interpolate parameter curve. Compared with traditional micro-segmented interpolators, NURBS interpolators are outstanding in improving machining accuracy and efficiency. It is indispensable for high-speed and high-precision machining. Therefore, the paper adopts the NURBS interpolator to accomplish the interpolation, and the relationship between interpolation position u i+1 and interpolation position u i is [20][21][22][23][24][25]: It can be observed that the parameter u i+1 of the next interpolation period can be obtained from Formula (3); the feedrate V(u i ) is indispensable. Therefore, the feedrate scheduling is the highlight in the section. The first part is the S-type flexible acceleration and deceleration algorithm, which is the basis of feedrate scheduling. The second part is the feedrate scheduling based on motion model. The third part is the double-direction look-ahead algorithm. The double-direction look-ahead algorithm is divided into backward planning and forward planning. The purpose of backward planning is to obtain the maximum starting speed under the constraint of sensitive speed area. And the purpose of forward planning is using the maximum starting speed to plan the duration time (τ 1 , τ 2 , τ 3 , τ 4 , τ 5 , τ 6 , τ 7 ) of each segment in the 7-segment S-type acceleration and deceleration for a path with a fixed length of L.

S-type acceleration and deceleration algorithm
Compared with trapezoidal acceleration and deceleration, the S-type acceleration and deceleration algorithm is more stable and has less impact on the drive system. It is suitable for high-speed and high-precision machining.
A complete S-type acceleration and deceleration has a total of 7 segments [26][27][28][29], named acc-acceleration (AA) segment, uniform acceleration (UA) segment, dec-acceleration (DA) segment, uniform speed (US) segment, dec-deceleration (DD) segment, uniform deceleration (UD) segment, and acc-deceleration (AD) segment. Figure 1 is a complete S-type acceleration and deceleration curve. The details of each segment are described in detail in Appendix 1, and the meaning of t i , T i ,τ i , a i (t), a i , v i (t), v i ,s i (t) s i , S total is also introduced in Appendix 1. Jerk, acceleration, velocity, and displacement are respectively obtained from Formulas (4), (5), (6) and (7) for the 7-segment S-type flexible acceleration and deceleration algorithm.
To sum up, in the 7 segments of S-type acceleration and deceleration, the formulas of jerk, acceleration, velocity, and displacement have common characteristics. The unified form of the formula is:

Feedrate scheduling
The process of path planning will not always appear all the 7 segments. Only when the path is long enough, all the 7 segments will exist in the drive system, and it is a complete acceleration and deceleration speed curve. If the path is not long enough, not all the 7 segments but parts of them will appear in the drive system. Therefore, it is important to determine which part appears. The detailed analysis is given below.
In the permutation and combination rules of probability theory, the combination of 7 segments are 2 7 different combinations in total. Some permutation and combination results do not exist in practice, even so, there are still many combinations. It is too complex to analyze for each situation. Therefore, the paper proposes 4 typical motion models to reduce the complexity of combination. For segments 4-7, they can all be mapped to these 4 types. For 3 segments and below, only the acceleration or the deceleration exist, which is easy to determine. The application of the four proposed motion models effectively reduces the complexity of (9) classifications. The content of the section is divided into two parts. The first part is the establishment of the motion model. The second part is the feedrate scheduling based on the proposed motion model.

Establish motion models
For a certain length of L path, 4 special models are analyzed.

Model a 7-segment S-type acceleration and deceleration
In the case, the given path L is long enough, and the speed curve is a complete 7-segment acceleration and deceleration curve. The maximum feedrate can reach V m , and the maximum acceleration can reach A m . The time of each segment satisfies the relational expression T 1 = T 3 = T 5 = T 7 , T 2 = T 6 . And T 1 , T 2 , T 3 , T 4 , T 5 , T 6 , T 7 can be obtained by using Table 1 and Formulas (4), (5), (6) and (7). Formula J m *τ 1 *τ 2 + J m *τ 1 *τ 1 = V m is introduced in Appendix 2. The gray represents the value is 0 or does not exist in the table. Figure 2a is the 7-segment S-type acceleration and deceleration.

Model b 4-segment S-type acceleration and deceleration
In the case, the given path L is not long enough, and the speed curve is a 4-segment acceleration and deceleration curve. This speed curve only includes AA segment (t 0 -t 1 ), DA segment (t 2 -t 3 ), DD segment (t 4 -t 5 ), AD segment (t 6t 7 ). The maximum feedrate can reach V m ' and not reach V m , and the maximum acceleration can reach a m or just reach A m . The time of each segment satisfies the relational expression T 1 = T 3 = T 5 = T 7 , T 2 = T 6 = 0, T 4 = 0. And T 1 , T 3 , T 5 , T 7 can be obtained using Table 2 and Formulas (4), (5), (6) Velocity Accelaration Jerk    Figure 2b is the 4-segment S-type acceleration and deceleration.

Model c 6-segment S-type acceleration and deceleration
In the case, the given path L is not long enough, and the speed curve is a 6-segment acceleration and deceleration curve. This speed curve only includes AA segment (t 0 -t 1 ), UA segment (t 1 -t 2 ), DA segment (t 2 -t 3 ), DD segment (t 4t 5 ), UD segment (t 5 -t 6 ), AD segment (t 6 -t 7 ). The maximum feedrate reach V m /V m ', and the maximum acceleration can reach A m . The time of each segment satisfies the relational expression T 1 = T 3 = T 5 = T 7 , T 2 = T 6 , T 4 = 0. And T 1 , T 2 , T 3 , T 5 , T 6 , T 7 can be obtained by using Table 3 and Formulas (4), (5), (6) and (7). Figure 2c is the 6-segment S-type acceleration and deceleration.

Model d 5-segment S-type acceleration and deceleration
In the case, the given path L is not long enough, and the speed curve is a 5-segment acceleration and deceleration curve. This speed curve only includes AA segment (t 0 -t 1 ), UA segment (t 1 -t 2 ), DA segment (t 2 -t 3 ), DD segment (t 4 -t 5 ), AD segment (t 6 -t 7 ) or DD segment (t 4 -t 5 ), UD segment (t 5t 6 ), AD segment (t 6 -t 7 ), AA segment (t 0 -t 1 ), DA segment (t 2 -t 3 ). Taking the UA segment as an example, the maximum feedrate can reach V m ' which is small than V m , and the maximum acceleration can reach A m . The time of each segment satisfies the relational expression T 1 = T 3 , T 5 = T 7 , T 6 = 0, T 4 = 0. And T 1 , T 2 , T 3 , T 5 , T 7 can be obtained using Table 4 and Formulas (4), (5), (6) and (7). Figure 2d is the 5-segment S-type acceleration and deceleration.

Feedrate scheduling
As shown in Fig. 3, the feedrate scheduling is introduced in the section based on the four proposed motion model in detail.
L is the given path length. Model c is a 6-segment acceleration and deceleration curve, and L6 is defined as the critical displacement when the velocity just reaches the maximum speed V m . Model d is a 5-segment acceleration and deceleration curve, and L5 is defined as the critical displacement when the acceleration can reach the maximum acceleration A m but the velocity cannot reach V m . Model b is a 4-segment acceleration and deceleration curve, and L4 is the critical displacement when the acceleration just reaches the maximum acceleration A m .
it is a complete acceleration and deceleration curve of 7 segments. τ1-τ7 can be solve using Table 1.
is the displacement when the velocity just reaches the maximum speed V m , so the displacement L only reaches the speed V m ' which is smaller than V m . Because the displacement L5 must be able to reach the maximum acceleration A m , the displacement L must be able to reach the maximum acceleration A m . Consequently, L belonging to the interval [L5, L6] is enough long to reach A m but not enough long to reach V m . Therefore, category 2 belongs to model c and (τ 1 , τ 2 , τ 3 , τ 4 , τ 5 , τ 6 , τ 7 ) is solved by the formula in Table 3. For example, τ 1 can be solved by formula A m /J m and τ 2 can be solved by formula s 3 = f(τ 1 , τ 2 ) = L/2.

Category 3
When L4 < L < L5, because the displacement L4 is just length to reach A m , the displacement L must be able to reach the maximum acceleration A m . The displacement L5 is enough long to reach A m but not enough long to reach the maximum speed V m . Consequently, L belonging to the (L4, L5) interval is enough long to reach A m but not enough long to reach V m . Thereby, category 3 belongs to model d and (τ 1 , τ 2 , τ 3 , τ 4 , τ 5 , τ 6 , τ 7 ) can be solved by the formula in Table 4. τ 1 is solved by formula A m /J m , τ 2 is solved by formula J m *τ 1 *τ 2 + J m *τ 1 *τ 1 = V m ', and τ 3 is solved by formula Category 4 When L < L4, because L4 is the displacement just reaching the maximum acceleration A m . L belonging to the (0, L4) interval is only enough to reach a m which is smaller than A m . Therefore, category 4 belongs to model b. (τ 1 , τ 2 , τ 3 , τ 4 , τ 5 , τ 6 , τ 7 ) can be solved by using Formula in Table 2 and τ 1 can be solved by using A m /J m .

Backward planning
Path planning is divided into backward planning and forward planning [26]. Backward planning is introduced in the part, and forward planning is introduced in next part.
The purpose of backward planning is to obtain the maximum starting speed with a limit of sensitive speed area. The paths are planned in groups not one by one, and how to get the maximum starting speed with the end speed under the feedrate constraint [30][31][32]. This is the meaning of backward planning. The backward planning process is as follows: Table 4 The duration of each segment of 5-segment S-type acceleration and deceleration The relationship between τ 1 -τ 7 and path length L in feedrate scheduling. L is the length of given path, and L6 is defined as the critical displacement just reaching the maximum speed V m in model c.
L5 is defined as the critical displacement that can reach the maximum acceleration A m but cannot reach V m in model d, and L4 is defined as the critical displacement just reaching the maximum acceleration A m in model b Step 1: Getting the end speed V e of last path.
Step 2: For the last path with a given length of L n , V e is regard as the starting speed V s1 . The last path is planned reversely to obtain the end speed V e1 with the feedrate scheduling introduced in previous section.
Step 3: V es2 is the smallest velocity between the end velocity V e1 and the sensitive velocity V sensor . For the penultimate path with a given length of L n-1 , V es2 is the starting speed of L n-1 . L n-1 is planned reversely to obtain the end speed V e2 .
Step 4: Repeat step 3 until the first path with a given length of L 1 is planned reversely and the end speed V en is obtained. This end speed V en is the maximum starting speed V s of the set of paths.

Forward planning
In the previous section, the starting velocity V s of a set of paths is obtained by backward planning. With V s planned by the backward planning and the feedrate scheduling based on Tables 1, 2, 3 and 4, the durationτ 1 , τ 2 , τ 3 , τ 4 , τ 5 , τ 6 , τ 7 of each segment of the 7-segment S-type acceleration and deceleration can be obtained for a given length L [32,33].

Interpolation pre-compensation based on discrete inverse transfer function
In Sect. 2.2, the feedrate scheduling will be generated with the proposed models a-d, and the path will be planned by forward planning and backward planning in Sect. 2.3. After that, CNC system will send the planned path to drive system with Formula (3) in Sect. 2. In Formula (3), u i is the current NURBS interpolation parameter, V(u i ) is the scheduled feedrate. The next interpolation parameter u i+1 will be obtained with u i , current feedrate V(u i ), and interpolation period T s . Each motion axis of CNC system will move to the commanded position with the obtained u i+1 and NURBS curve Formula (1), finally. However, it is impossible to precisely move the axis consisting of the drive system to the commanded position. The reason will be introduced in Sect. 3. Dispersed micro-segment from CAM is input to the CNC system. The CNC system receives micro-segment from CAM and transmits the commanded position P r to the drive system. The output of the drive system is the actual position P a . As shown in Fig. 4.
The commanded position P r output by the CNC system can accurately express the contour of the workpiece, but the actual position P a which is output from the drive system cannot completely follow the commanded position P r . As a result, the actual position P a and the commanded position P r are not exactly same, leading to a reduction in machining accuracy [34]. This machining deviation is called contour error [35][36][37][38][39], as shown in Fig. 5.
The difference is caused by the unsteady-state response of the drive system, which is a transient state before steadystate response of drive system. The paper proposes interpolation pre-compensation theory to avoid contour errors caused by the unsteady-state response. In Sect. 3.1, the reason of leading contour error is analyzed. In Sect. 3.2, before proposing the pre-compensation theory, transfer function discretization theory is induced at first, and pre-compensation theory is introduced in Sect. 3.3. The distance between q and p is the contour err. The distance between q and r is the following error

Response analysis of drive system
The drive system consists of a closed-loop transfer function G(s) [40][41][42][43]. For a typical drive, system is composed by a position loop and a velocity loop and PMSM mathematical model. The control block diagram is shown in Fig. 6 [44]. The position loop is controlled by the P controller, and the velocity loop is controlled by the PI controller. Therefore, the closed-loop transfer function G(s) of drive system can be calculated as In this situation [45], the drive system can be reasonably approximated as a first-order inertial system (FOIS). The step response of the closed-loop transfer function G(s) is depicted in Fig. 7. From Fig. 7, it is easy to find that, before reaching steady-state, the low bandwidth zone of FOIS with time constant T is about (3 ~ 4)T. The low bandwidth zone can lead to contour error, especially in the field of highprecision machining. T is generally greater than interpolation period T s . Therefore, it is difficult, even impossible, to reach steady-state within T s . To address this issue, the drive system is simplified as a first-order inertial system in Chen et al. [45,46], while the simplify model is not the real mathematical model. Simplifying a model inevitably introduces model errors, even if it is a higher-order simplified system. In the paper, a new theory is proposed to the eliminate the model error.

Transfer function discretization theory
The response of transfer function from unsteady-state to steady-state causes contour errors. Therefore, the article focuses on the unsteady-state process of transfer function to solve the issue. Before discussing the issue, a brief explanation of transfer function discretization theory is induced at first.
The CNC system is a linear time invariant systems (LTI) [46][47][48]. The response of LTI system is:  . x(t) is the input signal of transfer function. y(t) is the output signal of transfer function. In the field of machining, CNC will plan the path before sending the commanded position, which is x(t). The planned path is employed by backward planning and forward planning in Sect. 2.3, and it is continuous. While it is impossible to send the planned position all the time by CNC, CNC has many tasks to perform, and sending position is just one of them. Therefore, CNC will periodically send the commanded position in a time interval, which is the interpolation period T s . Hence, the essence of the CNC system is a discrete time system, not a continuous time system. The smaller the interpolation period, the higher the interpolation accuracy. No matter how small the interpolation period T s is, the series of commanded position are still dispersed. Based on the analysis, x(t) should be replaced by x(n). The form of Formula (13) should be: However, the drive system which is a transfer function is a continuous time system h(t). It can be seen from the above Formula (14), the transfer function h(t) can be degenerated into a discrete-time system h(n). The form of Formula (14) should be:

Theory 1
The transfer function is the ratio of the output function to the input function in Z-transform [46], which is . If the input function is a unit pulse which is impulse signal, the Z-transform of unit pulse is 1. Thereby, the transfer function is the output in Z-transform, Using the Theory 1 and Formula (15), the transfer discrete function h(n) can be obtained.   (

15) y(n) = x(n)*h(n)
Using Theory 2, the serious of commanded position P r (i) can be regarded as a serious of discrete sequence with different amplitudes and different time-shift. The discrete transfer function theory is supported by Theory 1 and Theory 2.

Discrete inverse transfer function algorithm
In Sect. 3.2, the conclusion that the transfer function which is the radio of output and input in Z-transform is the response of the unit pulse is introduced. In Sect. 3.1, the transition process from the unsteady-state response to the steady-state response leads to contour errors. To address the issue, the theory of discrete inverse transfer function is proposed.
According to Theory 1 in Sect. 3.2, the discrete transfer function H(n) is the output of transfer function with a unit pulse δ(n) input. But in machine tools, there are magnetic losses in the actuators and mechanical losses including frictional losses. Therefore, a unit pulse is only enough in simulation but not enough for a machine tool to obtain the discrete transfer function H(n). In order to minimize the effect of magnetic and mechanical losses when obtaining the discrete transfer function by inputting unit pulse, a set of unit pulses are input to the actuator of CNC machine.
From the conclusion of Sect. 3.1, the drive system response can be divided into two parts, the steady-state response and the unsteady-state response. H(m) is defined as the unsteady-state response of discrete transfer function with a unit pulse input and the steady-state response of unit pulse is 0. The length of unsteady-state response is m. Therefore, H(m) is the discrete transfer function of drive system. Besides  Fig. 8 The basic operations of LTI system. a The delay of LTI system. b The addition of LTI system. c The multiplication of LTI system

(a) (b) (c)
In summary, according to Formula (15), if there are no magnetic and mechanical losses, H(m) can be obtained by directly sampling the output of the actuator with a input unit pulse. And H(m) can be obtained by matrix equation (17) with a set of unit pulses in the case of magnetic and mechanical losses.

(16) H(m) = (h(1), h(2), h(3) ⋯ h(m)) T
Because of the last element Δx(0) = 0 in the last line of matrix Δx(i), the last raw of matrix Δx(i) is all zero in Formula (19). Therefore, h(m) will not work and will not affect the discrete transfer function H(m). According to the principle of one-to-one correspondence between the matrix equation and its augmented matrix, the augmented matrix of the matrix equation (19) is,

And correspondingly, the unknown vector H(m) to be solved for is
Using the elementary transformation principle of the matrix and Δx(0) = Δy(0) = 0, augmented matrix (20) is simplified to Δx(i)-Δx(i-1) is the impulse pulse in the ith interpolation period, and Δy(i)-Δy(i-1) is represented by the y inc (i). Augmented matrix (22) can be simplified as matrix (23). According to the principle of elementary transformation of matrices, matrix (23) can be further simplified as matrix (24), However, a set of unit pulses cannot be directly input to the actuator, nor are they regular inputs in CNC machine. Assuming that the regular inputs are x(n) = {x(0),…., x(n-1)}, which are the commanded position. x(0) is usually regarded as the initial value, named init x , and the increments Δx(i) which is impulse is x(i) relative to x(0), that is x(i)-x(0), i = 0,…,n-1. The outputs y(n) is same to x(n). Obviously, Δx(0) = Δy(0) = 0.
On the basic of Formula (15), H(m) can be obtained by Formula (19) with the known Δx(i) and Δy(i).

And obviously, the discrete transfer function H(m) is
When the motor is in enabled state, a slight vibration will be kept all the time. Considering the long-term existence of jitter, the average of the discrete transfer function H(m) is adopted and shown in Formula (26), For the law of conservation of energy, the sum of the output should be a unit pulse if the input is a unit pulse. However, the sum of the output is not a unit pulse duo to the data accuracy error, the H(m) should be normalized. P r , P a , and P c represent the commanded position, the actual position, and the pre-compensated position. Generally, the commanded position P r is transmitted to the drive , P a is not in accordance with P r . If the consistency is achieved between P a and P r after processed by H(m), P r must be compensated and the P c is the pre-compensated position. That is, the commanded position P r is replaced by the precompensated position P c , and actual position P a is replaced by the commanded position P r . The pre-compensation principle for the commanded position will be introduced in detail. If a constant value is continuously input to the transfer function, it is always a steady-state response after a while. Only when the input changes, will the unsteady-state response appear. Therefore, the concept of input increment W(n) which is named weight is introduced.
The input increment W(n) is the increment of the current commanded position P r (i) relative to the first commanded position P r (0) which is the initial value. Based on the above analysis of pre-compensation principle, W(n) is the output of the transfer function. Similarly, PW(n) is the increment of pre-compensation position and is the input.

Then with the discrete inverse transfer function H(m), the relationship between input PW(n) and output W(n) is
Then the pre-compensation position P c (i) which is inputted to the drive system is: system which is a transfer function, and the actual position P a is output. Based on the analysis of the response of transfer

Simulation and experiment
To validate the effectiveness and functionality of the proposed algorithms, the simulations and experiments are conducted in this section. Butterfly and spiral curves are applied to show the details of the PDIT and compared to the other two methods.

Simulation
In the section, simulation is used to evaluate the performance of the proposed method. FOIS is a method to improve the contour accuracy and is used for comparison, as is the contour error estimation and compensation (CEEC) method proposed by Hu et al. [49]. FOIS is a contour compensation method based on a mathematically simplified pattern of the first-order inertia of the drive system. The pattern consists of position loop, velocity loop and the mathematical model of PMSM. The contour error is deduced by the FOIS theory and compensated to the commanded position. In CEEC, Newton's method is used to analyze the contour error, and the iterative learning control (ILC) [50] method is used to compensate the contour error. In order to avoid the interference of mechanical errors such as bearing clearance, the (31) performance of the proposed PDIT is first evaluated by simulation, and a simplified drive system model based on Simulink module is adopted [51][52][53][54]. Figure 9 is the simplified drive system. Since the same method applies to the XYZ axis, the X axis is used as an example to introduce. In X axis, the K p , K i , and K d parameters are 1200, 300, and 5 respectively.
A butterfly and spiral graphics are used for simulation to evaluate the PDIT proposed in this paper. The butterfly is used in the first simulation, and the spiral is used in the second simulation. In Sect. 2, the parameter interpolation based on Taylor expansion and the feedrate scheduling are used to obtain the position C(u i+1 ) of next interpolation period. The PDIT proposed in Sect. 3 is used to pre-compensate the position C(u i+1 ) of the next interpolation period.
In Fig. 10, simulations are carried out to verify the improvement of the butterfly angle among the six corners marked 1-6. In the figure, the red trajectory represents the compensation effect of PDIT theory proposed in the paper, the blue trajectory is the CEEC theory, and the black trajectory is the FOIS method. The black point is the ideal commanded position, which is used as a reference for comparing the PDIT, FOIS, and the CEEC theory. In Fig. 10, at the selected corner of the butterfly trajectory, it can be clearly seen that the proposed PDIT is closer to the ideal interpolation curve than the CEEC. The CEEC trajectory is closer than the FOIS trajectory.
At the same time, the contour error of all the six corners is given in Fig. 11. The contour error of compensation by PDIT theory is almost a horizontal curve, while the CEEC and FOIS theory is a curve that fluctuates up and down. The contour error of PDIT theory is significant small than that of CEEC theory, and the contour error of corner in CEEC is small than the FOIS. The contour error data analysis of the six corner is given in Table 5 with PDIT FIOS and CEEC. The maximum contour error of PDIT theory is 0.4um, the maximum contour error of CEEC theory is 2.9um, and the FOIS is 15um. The interpolation Fig. 9 The mathematical simplified model of drive system pre-compensation theory based on the discrete inverse transfer function proposed in this paper has a significant improvement effect in terms of contour.
In the second simulation, a spiral curve is used as an example. The spiral curve [49] Formula (32) is as follows. In Fig. 12, the red trajectory is the actual trajectory compensated by PDIT theory proposed in this paper, the blue trajectory is the CEEC theory, and the black point is the ideal commanded position. Figure 12 shows that the red trajectory of the PDIT theory is closer to the ideal curve than the blue curve of the CEEC theory. Figure 13 and Table 6 are the contour error of PDIT theory and CEEC theory. The contour error of PDIT is smaller than CEEC. For the contour error of the spiral curve, the contour error of PDIT is 0.000849-0.147475 mm, and the contour error of CEEC is 0.001189-0.496462 mm.

Experiment
In addition to the simulation studies applied in the previous section, experiments are performed to verify the contour accuracy improvements of the proposed method. As shown in Fig. 14, the experimental platform is a CNC machine tool consisted of two axes with build-in position loops, and the motion parameter, which is K p , K i and forward parameter of position, velocity, and current loop in servo, are in Table 7. The position loop is a proportional controller, while the velocity loop and the current loop are both proportional and integral controllers. The axes are controlled by a PC-based CNC system with a real-time Linux operating system and EtherCAT network. EtherCAT is the real-time industrial  In order to minimize the effects of magnetic and mechanical losses, a total of 50,000 inputs are used to obtain the discrete transfer functions H(m) and named x(0), x(1), …, x (49,999). The 50,000 inputs are transmitted to the actuators by EtherCAT network with object 607A, and the duration is 50 s in total if the interpolation period is 1 ms. Correspondingly, the output of drive are y(0), Fig. 13 Contour error of the entire path of the spiral curve  (27). The method proposed in the paper focuses on reducing the effect of unsteady-state response caused by the increment between adjacent positions. The position increment changes significantly at the corners, resulting in an overshoot with unsteady-state response. There are many types of corners in the butterfly, so the butterfly is still used for processing experiments. The butterfly is machined by the experimental bi-axial table in Fig. 14. All the selected six corners used to evaluate the performance of the proposed PDIT theory are same to the simulation in Sect. 4.1. The contour of all the six corners of the experiment is shown in Fig. 15, which is comparison among PDIT, FOIS, and CEEC methods. The red line is actual axes trajectory which is compensated by PDIT, the blue line represents the output of commanded position in servo with FOIS, the purple line is the actual trajectory compensated by CEEC, and the black line is the commanded position which is the ideal position without contour error. As we can see, in Fig. 15, the red line is closer to the black line than the blue line and purple line. It is obvious that the proposed PDIT theory can significantly improve the contour accuracy in corners. The contour error is shown in Fig. 16, and the contour error data analysis of the six corner is given in Table 8. The maximum contour error of PDIT theory is 0.038 mm, the maximum contour error of CEEC theory is 0.09 mm, and the maximum contour error is 0.071 mm in FOIS. In terms of mean contour error, the proposed PDIT improves the contour accuracy of the selected six corners by an average of 28.712% compared to FOIS, and 57.35% compared to CEEC. The interpolation pre-compensation  15 Contour comparison of the selected six corners with PDIT CEEC and FOIS methods theory based on the discrete inverse transfer function proposed in this paper has a significant improvement effect in terms of contour. It can be seen that the compensation effect of CEEC is slightly worse in experiment than in simulation, because the pre-compensated position P c (i) in the experiment is approximated. The position in servo is based on encoder and is represented by the number of pulses, so the command position P r (i) and the actual position P a (i) are both integers in servo. However, the pre-compensated position P c (i) obtained by CEEC from the actual position P a (i) is float. Therefore, in experiment, the float P c (i) needs to be rounded to the nearest integer and then sent to servo, while  In the second experiment, in order to maintain consistency between the simulation and experiment, the spiral curve is still taken as the example. In Figs. 17 and 18, the red line is about the PDIT theory proposed in this article, the purple line is about the CEEC theory, and the blue line is the FOIS theory. The black point is the ideal commanded position in Fig. 17. Figure 17 shows that the red trajectory of  the PDIT theory is closer to the ideal curve than the purple curve of the CEEC theory and the blue curve of the FOIS theory. Figure 18 is the contour error after compensation by PDIT CEEC and FOIS methods for the full path of a spiral curve. As we can see in Fig. 18, The red line representing the compensation effect of PDIT theory proposed by the paper is minimal. And Table 9 is the contour error data analysis of the entire spiral curve. The contour error of compensation by PDIT is smaller than FOIS and CEEC. For the contour error of the spiral curve, the maximum contour error of PDIT is 0.0739 mm, the CEEC is 0.138 mm, and FOIS is 0.0772 mm. In terms of mean contour error, the proposed PDIT improves the contour accuracy of the entire spiral curve by an average of 47.372% compared to CEEC, and 5.5% compared to FOIS.

Conclusions
The NURBS interpolator with pre-compensation based on discrete inverse transfer function is important to guarantee high-efficiency and high-accuracy machining. In conventional methods, high-speed machining and high-precision machining are contradictory, and the proposed PDIT breaks this point. Without decreasing feedrate, the proposed PDIT can effectively improve the machining contour accuracy, especially at the corners. With the proposed method, the NURBS interpolator can pre-compensate for contour errors due to the unsteady-state response of the servo system. As compared to the other two algorithms, the proposed method can significantly improve machining accuracy. In simulation, the commanded position with PDIT can completely avoid the influence of unsteady-state response, which is helpful to realize high-speed and high-precision machining in CNC machining field. Even in the experiment, the contour accuracy is improved by 28.163% and 57.35% compared to the other two methods, respectively. If the interpolation period is smaller, the contour accuracy is further improved. Moreover, PDIT only includes basic addition, subtraction, multiplication and division operations, without complex operations such as integral in FOIS. Therefore, PDIT does not tax the CPU and hardly affects the real-time performance of the NURBS interpolator.

Appendix 1
In the acc-acceleration segment, the jerk is a positive constant value, and the acceleration increases at a constant value. In the uniform acceleration segment, the value of the is 0, and the acceleration reach a maximum value A m . In the dec-acceleration segment, the jerk is a negative constant value, and the acceleration decreases at a constant value. In the constant speed segment, the values of jerk and acceleration are both 0, and the velocity reach the maximum value V m . In the dec-deceleration segment, the jerk is a negative constant value, and the acceleration decreases at a constant value. In the uniform deceleration segment, the value of jerk is 0, and the acceleration reach a minimum value − A m . During the acc-deceleration segment, the jerk is a positive constant value, and the acceleration increases at a constant value until it reaches zero. t i represents the critical moment in segments, so t 1 -t 7 are the end moments of the first to seventh segments. T i represents the duration of each segment, so T 1 = t 1 -t 0 is the duration of the first segment, and T 7 = t 7 -t 6 is the duration of the seventh segment. Others are the same. τ i represents the time difference between the any moment in one segment and the starting moment of the same segment. For example, in the interval t 6 < t < t 7 , τ 7 = t-t 6 . The acceleration in each segment is represented by a i (t), so a 1 (t)-a 7 (t) are the acceleration of the first to seventh segments. The acceleration at the end of each segment is represented by a i , so a 1 -a 7 are the final accelerations of the first to seventh segments. The speed in each segment is represented by v i (t), so v 1 (t)-v 7 (t) are the speed of the first to seventh segments. The speed at the end of each segment is represented by v i , so v 1 -v 7 are the final speeds of the first to seventh segments. The displacement in each segment is represented by s i (t), so s 1 (t)-s 7 (t) are the displacement in the first to seventh segments. The displacement at the end of each segment is represented by s i . The total displacement of the entire 7 segments is represented by S total .

Appendix 2
In model a, after the first three segments, which are AA, UA, and DA, the speed increases from 0 to the maximum speed V m . The speed increment in AA segment is Δv 1 . The speed increments in UA and AD segment are Δv 2 and Δv 3 , individually. And: t 1 = T 1 , t 2 = T 1 + T 2 , t 3 = T 1 + T 2 +T 3 , T 1 = T 3 = T. 0-t 1 : t 1 -t 2 : (33) ▵ 1 = ∫

Declarations
Ethics approval Approval.

Consent to participate Consent.
Consent for publication Consent.

Competing interests
The authors declare no competing interests.