Prediction and compensation strategy of contour error in multi-axis motion system

The contour error of the machined parts is an important index to evaluate the accuracy of CNC machining. Considering a multi-axis servo control system, a predictive compensation strategy for contour error is presented in this paper. First, the offline identification method is adopted to establish the transfer function of each motion system. In this case, the relationship between the interpolation command and the feedback position of the grating is determined for the selected machine tool. Thus, before machining, the trajectory information of each motion axis can be predicted according to the interpolation instruction and the transfer relationship. It can be converted into the contour trajectory of the part through kinematic analysis, so as to predict the contour error. Finally, the predicted contour error is compensated into the command trajectory to ensure the contour accuracy of the part to be machined. Compared to existing methods, our method can effectively reduce the trial-produced cycle of parts and avoid unnecessary waste of processing materials. Moreover, without increasing the complexity of the control system and greatly reducing the machining efficiency, the dynamic error caused by the dynamic characteristics of the shaft is reduced. Taking the starfish pattern and the contour line of the impeller as the milling processing experiment case, the contour error after compensation is greatly reduced compared with the contour error processed by the original command. Therefore, the predictive compensation method proposed in this paper can significantly improve the machining accuracy. In addition, it also has application value for trial production of complex curved parts.


Introduction
Multi-axis numerical control machining method is a typical representative of the development of CNC technology. It refers to the simultaneous coordinated movement of multiple axes under the control of a computer numerical control system (CNC) to improve the machining accuracy, quality, and efficiency of free-form surfaces, which have a wide range of applications in the field of precision machining and manufacturing [1]. The key performance index of a multi-axis motion control system is the accuracy of contour motion. Especially during the machining of high-speed and large-curvature parts, the dynamic performance mismatch between each axis will be enlarged as the speed and curvature increase. Therefore, the research on the high-precision machining of multi-axis, high-speed, and variable-curvature parts has important theoretical significance and remarkable engineering application value.
According to the existing research, the main reasons that lead to the contour accuracy of the machined parts include the quasi-static error represented by the geometrical size of the machine tool structural components, the material instability, the clamping positioning error and the thermal deformation, the dynamic error caused by the servo lag, and the dynamic adaptation [2]. Among them, the dynamic error is caused in the movement and changes with the trajectory and speed, which often exceeds the influence of quasi-static error on the accuracy of part contour for multi-axis in the high-speed, complex trajectory machining. Therefore, this paper will uniformly consider the dynamic error that affects the contour accuracy of the machined part as the contour error and focus on the prediction and compensation of contour errors.
The contour error is defined as the vector from the actual position point of the coordinated motion of each motion axis to the vertical foot of the expected contour curve planned by the interpolator. The size of its vector modulus is the standard to measure the quality of multi-axis synchronous machining. Scholars from various countries have adopted different solutions to improve the contour accuracy of the multi-axis motion system and ensure processing efficiency. A considerable number of scholars start from the design of the controller and propose control strategies for indirect control and direct control of contour error respectively. The representative indirect control methods are zero phases' error tracking controller [3], disturbance observer [4], model predictive control [5], which are to improve the followability of each axis, so as to broaden the frequency response bandwidth of each motion axis.
However, the load inertia of each motion axis is different, and it is difficult to coordinate the best dynamic characteristics of each axis. Therefore, Koren [6] proposed a cross-coupled control, which could balance the tracking control effect of each axis to achieve the purpose of directly controlling the contour error. In the beginning, most of the cross-coupling controllers (CCC) were based on PID control [7]. With the development of modern control theory, many experts and scholars applied fuzzy control [8,9], adaptive control [10], and robust control [11] methods to achieve feedback control of contour error. In addition, the combination of single-axis controller changes and cross-coupling control has become the current research development trend [12].
At the same time, the structure of the contour error feedback controller has also changed. Yang et al. proposed a twolayer CCC mechanism, which divides the three-dimensional space contour error into two layers of control [13]. At the bottom layer, a two-axis CCC controller was used to control the XY plane contour, and the top layer was designed separately to control the three-dimensional contour error. Each layer could be designed individually, which made the design of the contour error controller more flexible. However, the direct feedback control of the contour error made the whole system have a strong coupling response, which reduced the stability of the system. Secondly, the tool path existed in the form of discrete points after numerical control interpolation, so the servo layer didn't save all the information of the trajectory. At the same time, in order to realize the versatility and stability of the CNC system and the servo system, the CNC manufacturers have carried out relatively independent design and system packaging respectively, so the users are difficult to design the controller according to their needs. Therefore, the method of reducing contour error from the perspective of controller design is difficult to be applied and popularized.
Some scholars who study numerical control interpolation consider contour error as a constraint condition in feed rate planning [14]. Jia et al. established the relationship model of feed rate, trajectory curvature, and tracking error and realized theoretical error constraint by limiting the feed rate [15]. Dong et al. simplified the drive system to a second-order system model and took the contour error as a constraint condition of real-time interpolation to ensure that the error didn't exceed the set constraint limit [16]. In general, the above methods increase the processing accuracy at the cost of efficiency and increase a large amount of computational burden. From the perspective of system parameter correction, Wu et al. sequentially tuned the control parameters in the position loop, speed loop, and current loop to minimize the mismatch of the dynamic characteristics of each axis, so as to obtain the improvement of the contour accuracy [17]. This method required a large number of experiments. As the number of motion axes increases, parameter tuning will become more difficult.
In the multi-axis machining of complex parts, some scholars proposed the tool path optimization method [18,19]. Huang et al. designed an iterative tool path compensation algorithm to reduce machining error and avoid unnecessary interference by modifying the tool path [20]. Hou et al. propose a new blade geometry model, including concaveconvex curve reconstruction based on the secant compensation method [21]. This method could significantly reduce the machining error for the milling of thin-walled parts. Calleja et al. proposed a novel integrated adaptive compensation machining method based on touch-trigger which could reduce the flexible deformation and improve machining accuracy [22]. This paper will propose an offline predictive compensation method. The dynamic characteristics can be obtained by determining the discrete transfer function of each motion axis with the system identification method. Input the interpolation instruction of the part to be machined into the discrete transfer function, convert it into the workpiece coordinate system through kinematics to predict the contour error caused by the dynamic mismatch of each motion axis, and compensate the error to the interpolated discrete trajectory before the actual machining. The method verified by experiments can effectively reduce the contour error of variable curvature parts in the machining process. Especially in the state of lack of processing experience, for the first trialcut parts, it can effectively improve the pass rate of parts machining and avoid the waste of materials. Meanwhile, it can also avoid the fact that the process personnel adopt conservative parameters for machining, and meet the requirements of part accuracy at a very low speed.
The rest of this article is as follows. The second section introduces the identification method of the feed system transfer function of the machine tool. The third section is to evaluate and screen the identification model to determine the best prediction model for each axis position. The fourth section is to solve the predicted contour error through the instruction trajectory of the pre-machined part and the predicted actual trajectory. Predicted error is compensated. The fifth section verifies the effectiveness of the proposed offline prediction and compensation method for contour error through comparative experiments. Our conclusion is in the sixth section.

Identification method for the transfer function of machine tool feed system
System identification is to determine a model equivalent to the tested system from a set of given model classes on the basis of input and output data. The parameter model identification method is the most frequently used in system identification. Based on the presumption of the model structure, the parameters of the model are determined by minimizing the error quasi-measurement function between the model and the object. The basic principle is shown in Fig. 1. At the current time, the system output prediction ŷ is calculated according to the estimated parameter ̂( k − 1) at the last moment, the previous input u(k) and output y(k).
Prediction error is as follows: then the prediction error e(k) is fed back to the identification algorithm, and the estimated parameter ̂( k) is calculated at the k moment. Through the iterative loop, until the minimum value of the criterion function is reached, which means that the model estimated parameters ̂( ∞) are obtained, and the model output is also the best approximation to the system output under the criterion. Therefore, in order to achieve the accurate prediction of each axis position information, the following is introduced in turn from the structure of the identification model, the determination of the model order, and the design of the identification signal.

Prediction model structure analysis
In a multi-axis system, each axis is independently controlled, and each motion axis works in a limited frequency range under servo control. Therefore, the high-frequency flexible model part can be ignored in the dynamic mathematical model of the control system. The low-order linear link can be retained. The general structure of the control of the feed axis is shown in Fig. 2. The G code is processed by the CNC system to obtain the interpolation point. Through the three-loop control, the motor is controlled to move, driving the mechanical system. Therefore, this paper utilizes real-time acquisition of the interpolated instruction position signal and grating feedback signal to identify the relationship between them for realizing the prediction of the contour error of the part before machining.
According to the Wold decomposition theorem, this paper will use the prediction-error identification method (PEM) to abstract the identified link into the following discrete transfer function expression: where u k and y k in Eq. (4), respectively, represent the interpolation instruction of the model input and the grating feedback of the model output at the time k . G(z −1 ) is the linear model, H(z −1 ) is the noise model, e k is the Gaussian white noise signal. G(z −1 ) and H(z −1 ) from a different PEM model is equipped with a different structure. According to the characteristics of the feed system, three PEM models are pre-selected in this paper: extended auto-regressive model (ARX), extended auto-regressive moving average model (ARMAX), output-error model (OE), shown in Table 1.

Model order determination
According to the Wold decomposition theorem [23], any discrete stationary time series can be decomposed into the sum of two uncorrelated stationary series, one of which is a deterministic series and the other is a random series. At the same time, the deterministic sequence can be expressed as a linear function, which is G(z −1 ) in the Eq. (4), and its expanded form is as follows: where b i and a j are the coefficients of the numerator and denominator of the discrete transfer function, which can be obtained by identification using the least square method [24]. p and q are, respectively, the orders of the numerator and denominator. In previous applications, the commonly used method is to start from the smallest order and keep increasing the order until the accuracy of the model meets the requirements. Some scholars use the BIC criterion of the model to find the optimal order within the specified range, but the result can only be used as a reference for the model order determination [25].
Different from other studies, G(z −1 ) is the dynamic model of the servo feed system with a clear physical meaning in this paper. Therefore, the order of p and q can be determined by establishing the closed-loop transfer function of the theoretical position system.
First, the establishment starts in sequence from the inner ring to the outer ring. The innermost part is the current loop, which usually consists of a PI-controlled regulator, filter link, inverter equivalent link, and armature winding link, as shown in Fig. 3. Where T oi is the current feedback filter time constant, T inv and K inv are the time constant and gain of the inverter equivalent inertia link. L q and R s are the Model structure

Fig. 3
Current closed loop control system diagram armature winding inductance and resistance, K pi and i are the proportional coefficient and integral time constant of the PI controller. After approximating the inertia link in Fig. 3, it is shown in the following Fig. 4. The approximate condition is that the cutoff frequency of the open-loop frequency characteristic of the current loop ci needs to be satisfied Where the small inertia link is combined to obtain the time constant of the inertia link T 1 = T oi + T inv , K 1 = k pi k i nv∕ i , which is the equivalent gain. T 2 = L q ∕R s is the electrical time constant and K 2 = 1∕R s is the electrical gain. Generally, the zero points of the controller is canceled with the large time constant pole of the control object, that is i = T 2 , the transfer function of the final current loop is: The simplified current loop is brought into the velocity loop. Therefore, the final position transfer relationship is shown in Fig. 5.
Where T os is the time constant of the speed feedback filter, J is the equivalent moment of inertia of the controlled object, K s and s are the proportional and integral time constant of the speed loop PI controller. F(s) is the feedforward of the speed loop, which is theoretically composed of a differential link with a coefficient of 1 [26]. K p is the proportional coefficient of the position loop. The closed-loop transfer relation of the simplified position system is shown as Eq. (7).
G p (s) = P r P * r = K I K t K s s s 2 + K I K t K s + K I K t K p K s s s + K t K I K p K s T os s Js 5 + T os K I J s + J s s 4 + K I J s s 3 + K s K t K I s s 2 + K t K I K s + K t K I K p K s s s + K t K I K p K s It can be seen that the order of the numerator and denominator of the transfer function of the continuous system are 2 and 5, respectively. When converted into the discrete domain by the backward difference method, it can be incorporated into the equation s = 1 − z −1 ∕T , where T is the sampling period. The final discrete transfer function does not affect the order of numerator and denominator. Therefore, the values of linear functions p and q of the deterministic sequence are 2 and 5, respectively.

Design method of identifying signal
In the system identification process, different input signals will affect the identification effectiveness. The identification signal should have enough frequency domain information to completely cover the frequency band of the identified system. The signal should not be too large in power and amplitude, which should avoid the identification object working in the nonlinear region.
Currently, the common identification signals include step signal, square wave signal, sine signal, chirps signal, white noise, etc. But the power of these signals is basically concentrated in a certain frequency band, which cannot be evenly distributed in the whole frequency band. Therefore, this paper adopts the inverse M-sequence signal in the pseudorandom binary sequence signal under the condition that the normal tolerance range of the servo electromechanical drive and the sufficient bandwidth spectrum range can maintain continuous equal power energy input.
In Eqs. (8) and (9), the feedback weight coefficient is . n is the bits number of the shift register and the M-sequence is {m(k)} = m 1, m 2 ⋯ m k . S(k) is a square wave sequence with a period of 2 bits and the element of 0  The inverse M-sequence signal overcomes the static disturbance to the identification system from the DC component in the M-sequence. Besides, it has satisfying correlation characteristics and the power is evenly distributed on each frequency point, as shown in Fig. 6.
The characteristic parameters of inverse M-sequence are the number of signal bits n , signal amplitude a , the minimum time interval of state changes T d .
This paper takes the feed system of the KMC800 machine tool as the research object, and the maximum linear velocity of the feed axis is v max = 48m/min . The sampling frequency of CNC system software is f s = 500Hz . The bandwidth of the position loop of each axis according to the frequency sweep test is B < 100Hz.
In order to facilitate the observation of the frequency characteristics of the signal and eliminate the leakage error, Eq. (10) can be obtained.
where T f is the length of the sampling time, T im is the cycle period of the reverse M-sequence, and the relationship between the number of shift registers n and the cycle period is as follows: then, f s is the sampling frequency and N f is the number of samples. According to the Nyquist sampling theory, there is T d ≥ 2T s . On the other hand, as T d decreases, the spectral spacing decreases. The spectrum is close to the continuous spectrum, which can improve identification accuracy. Then we use q as an adjustment coefficient to get T d = q(1∕f s ) . And N f can set 2q(2 n − 1). Signal amplitude a mainly affects the power spectrum amplitude of the signal, which is mainly subject to physical limitations.
where l max is the maximum range of working. In order to prevent the system from entering the nonlinear region, the value of a is chosen to be 30 mm in this paper. At the same time, in order to realize that the signal bandwidth B im covers the bandwidth of the identified object, therefore: The frequency resolution ratio f im is set to 10 Hz, so a 4-bit register can meet the experiment requirement.

Evaluation and screening of identification models
In this section, the above design signals are input to the feed system of the KMC800 machine tool in the form of G-codes, which will ignore the influence of the CNC interpolation process on signal characteristics. Collecting the actual grating feedback data of each feed axis, the discrete transfer function of the corresponding machine tool X, Y, and Z feeds can be obtained according to three PEM models. Meanwhile, the order of numerator and denominator of a discrete transfer function is derived from Sect. 2.2. Because the accuracy of the identification model directly affects the prediction of the contour accuracy, it is necessary to select according to the evaluation index. The test signal in the time domain is shown in Fig. 7. From the time domain point of view, a time-domain signal is an input into the actual system and the identified model at the same time. By comparing the different identified models, the coefficient of determination and the mean square error between the predicted position value and the position value output by the actual system are obtained. Finally, the identification model which is closest to the actual system is selected. R 2 , MSE calculation expression is as follows: Among them, ŷ i is the predicted value, y i is the actual value, and ȳ is the average of the actual value. The prediction error of the discrete transfer function of the feed system identification is shown in the Table 2.
From the above analysis, it can be seen that the degree of fit of the z-axis OE model is low in the identification process. In other words, the identified model cannot truly reflect the dynamic characteristics of the actual structure. Finally, according to the evaluation index, the discrete transfer functions of X-, Y-, and Z-axis identification are determined as the following equation.
It is necessary to ensure the stability of the identification system while meeting the prediction accuracy. A sufficient and necessary condition for the stability of linear discrete systems is that all the poles of the system are located in the unit circle of the z plane.
From the above pole distribution, it can be concluded that the identified transfer function is stable. In general, the actual system is stable. Our aim is that the identified system should be infinitely close to the real system. If there is an unstable state, it needs to be identified again.

Contour error prediction and compensation
Contour error refers to the shortest distance from the current actual position to the desired trajectory. Compared with the tracking error from the current position to the desired tracking point, the contour error is coupled with the tracking errors of every single axis, which is more practical in the processing of high-precision complex curved surfaces. Contour error can directly reflect the contour accuracy of processing. As shown in Fig. 9, Pe1 has a larger tracking error than Pe2 , but has a smaller contour error.
In the process of NC machining, the curve information has been completely discretized by the interpolator. After interpolation, the output data points are discrete according to the interpolation period. Because the interval of discrete interpolation points is small, most of them are in 2 ms. And the error of bow height is almost negligible. So, the discrete data points according to the interpolation period are regarded as the expected trajectory.
In the previous section, by identifying the mathematical model of the selected machine tool feed system, the actual position of the machine tool feed axis during the machining process can be predicted. Our objective is to get the vector from the predicted actual position point to the vertical foot of the expected contour curve. Therefore, determining the foot point is the key. Next, a method for solving contour errors based on interpolation point sequence will be introduced. This paper will use two search strategies to determine the footsteps. The first step is to establish a storage area with the length of m, and then search for the instruction position of the closest point to the actual point in the area. The definition of the buffer area is shown in Fig. 10.
In Fig. 10, Pe i is the actual position point at time i . Pc i is the expected interpolation position at time i , Pc i−1 , and Pc i+1 represent the expected positions of the previous interpolation cycle and the next interpolation cycle, where n = (m − 1)∕2 , which fully include leading and lagging interpolation points. The length of the storage area m is user-defined. As the command trajectory speed increases, m needs to be expanded. After the first search, it is deteris the same distance from both sides, as shown in Fig. 11(c). Our goal focuses on solving the contour error, so the third case can be merged into the two cases. Finally, the calculation process of the above case is unified as  Therefore, in order to reduce the contour error in the actual machining process, the contour error value can be predicted according to the interpolation command. Before processing, the error has been compensated to the interpolation point at the corresponding position. The value of the interpolation after compensation Pn a is shown in Eq. (23).

Experimental verification and analysis of contour error prediction and compensation
KMC600 UMT turning-milling vertical machining center was selected as the machine tool. The starfish pattern with a large curvature variation and the profile of integral impeller with a wide range of applications were selected as the parts to be machined, as shown in Fig. 12.
In the starfish pattern machining, we adopted the B-spline curves interpolation method [27]. After feed rate planning, interpolation points are formed and input into the CNC system. The actual cutting maximum feed speed of the starfish pattern is 6000 mm/min. In the impeller machining, we adopted high-performance micro line block interpolation. G-code is generated by Ultra CAM. In addition, the CNC system has opened the interpolation point collection function, so that we can also get interpolation information. The spindle speed is 8000 rad/min, and the feed speed increases from 600 to 2000 mm/min.
In the two experiments, the processing material of the parts is aluminum alloy 2A70. The tool is made of a hard alloy ball end tool. The cylindrical radius, the radius of the ball, and the taper of the tool are 12 mm, 4 mm, and 4°, respectively.
According to the identification model of the feed system selected in Sect. 3 and the contour error solution method introduced in Sect. 4, the predicted contour error is shown in Fig. 13.
From Fig. 13, the predicted contour error value is consistent with the actual true contour error value. It can be seen from the data that most of the predicted contour error values can accurately express the contour error of the part in the actual machining process, except for a few points. Then the predicted contour error is used as Fig. 11 Three cases of searching the location of foot point Fig. 12 Milling experiment of contour error prediction and compensation the compensation vector to compensate to the interpolation point of the response position. To make the discussion clear, the compensation position diagram is partially enlarged, and the enlarged position is consistent with Fig. 13, as shown in Fig. 14.
Input the compensated interpolation point into the CNC system for machining. The contour error value after machining is shown in Fig. 15. The contour error defined at this time is the distance between the compensated machining tool point and the theoretical contour of the original design. In order to illustrate the effectiveness of compensation, this paper compared the contour errors of parts machined with the original interpolation points.
It can be seen from Fig. 15 that the contour error after compensation is obviously reduced. After calculation, the maximum value of the contour error before compensation is 0.0135 mm and the maximum value of the contour error after compensation is 0.0078 mm. The average error value is reduced by 61.76% after compensation. Therefore, the proposed method can effectively reduce the contour error. By predicting and compensating for the contour error before machining, the machining accuracy of the part can be obviously improved.