A Feedrate-Constraint Method for Continuous Small Line Segments in CNC Machining Based on Nominal Acceleration

: When the CNC machining of continuous small line segments is performed, the direction of the machine tool movement will change abruptly at the corner of adjacent line segments. Therefore, a reasonable constraint on the feedrate at the corner is the prerequisite for achieving high-speed and high-precision machining. To achieve this goal, a feedrate-constraint method based on the nominal acceleration was proposed. The proposed method obtains the predicted value of acceleration during the machining process by the machining trajectory prediction and acceleration ﬁltering. Then, the feedrate at the corner is constrained, according to the predicted acceleration. Speciﬁcally, for any corner of adjacent line segments, the proposed method assumes that the CNC machining of a short path centered on the corner is carried out at a constant feedrate. First, the actual machining trajectory is predicted according to the transfer function of the servo system. Then, the nominal acceleration, when the CNC machining is carried out to the corner, is calculated and processed by a low-pass FIR ﬁlter. Last, the feedrate-constraint value at the corner is obtained according to the nominal acceleration and the preset normal acceleration. The advantage of the proposed method is that it can be used for different machining paths consisting of long segments or continuous small segments and it has no special requirement for the accuracy of the machining path. As a result, the feedrate-constraint value obtained is reasonable and the smooth machining process can be ensured. The simulation results in both 2D and 3D machining paths show that the proposed method is insensitive to the length of the line segment and the angle of the corner, and the calculated feedrate-constraint value is close to the theoretical value, which has good stability and versatility. In contrast, the feedrate-constraint values obtained by conventional methods change abruptly along the machining path, especially in the 3D simulation, which will damage the machining quality. The experiment was performed on a three-axis CNC machine tool controlled by a self-developed controller, and a free-form surface workpiece was machined by a conventional feedrate-constraint method and the proposed method, respectively. The experimental results showed that the proposed method can make the feedrate of the machining process higher and more stable. Then, machining defects such as overcutting and undercutting can be avoided and the machining quality can be improved. Therefore, the article proposes a new method to constrain the feedrate at the corner of continuous small line segments, which can improve the machining efﬁciency and quality of the CNC machining.


Introduction
The continuous small line segments are widely used in CNC machining [1][2][3], which are universal and easy to calculate. However, the direction of the machining movement will change abruptly at corners of adjacent line segments, which brings many difficulties to the realization of precise and efficient CNC machining [4][5][6]. One significant difficulty is to constrain the feedrate reasonably at corners [7,8]. Too low a feedrate at the corner will increase the machining time, while too high will cause the machine tool to vibrate and damage the machining quality.
A reasonable constraint on the feedrate at the corner of continuous small line segments is very important for achieving high-speed and high-precision machining. Two commonly used constraint methods are the angle-constraint method and the curvature-constraint method. The angle-constraint method first calculates the angle of the corner of adjacent line segments and then substitutes it into a mathematical model to obtain the corresponding feedrate-constraint value. As in the literature [9,10], scholars assume that the magnitude of the velocity (a vector) of the machine tool remains constant during machining a corner and that the change of the velocity direction is completed within one interpolation cycle, which generates a steering acceleration during the process. Then, according to the relationship among the steering acceleration, the angle, and the velocity magnitude, the feedrateconstraint value can be achieved from the angle, under the constraint that the steering acceleration does not exceed the preset upper limit. The literature [11] assumes that there exists a virtual transition arc at the corner. Firstly, the radius of the arc is obtained according to the corner's angle and the maximum allowable value of transition error. Then, the feedrate-constraint value of the transition arc is calculated according to the radius and the preset upper limit of the normal acceleration, which will be used as the feedrateconstraint value at the corner. In the literature [12], considering that the length of each line segment is usually not exactly the sum of integer interpolation cycles, the interpolation across line segments at the corner will cause the trajectory error in machining. According to the relationship among the trajectory error, the corner's angle, and the feedrate, the feedrate-constraint value can be obtained from the angle under the constraint that the trajectory error does not exceed the preset upper limit. The angle-constraint method above has good performance for the area with long line segments in the machining path (hereafter referred to as "Line Area"), but for the continuous small line segments discretized from the curve (hereafter referred to as "Curve Area"), the calculated angle will be small and unstable due to the short line segments, which will eventually cause the feedrate-constraint value to be overlarge and unstable.
The curvature-constraint method first calculates the radius of curvature of the continuous small line segments and then obtains the feedrate-constraint value based on the radius of curvature. For example, the Syntec CNC [13] uses a virtual arc to approach the path of two adjacent line segments and calculates the chord error. If the chord error is small enough, the line segments will be recognized as Curve Area and its radius of curvature will be calculated. Similarly, the literature [14] estimates the radius of curvature corresponding to the current corner based on the geometric information of several adjacent small line segments around this corner, while the literature [15] first groups small line segments by the length of the line segment and the angle of the corner and then calculates the radius of curvature corresponding to each group. After that, the feedrate-constraint value can be obtained from the radius of curvature and the preset upper limit of normal acceleration. For Curve Area in the machining path, the above curvature-constraint method can achieve good results, but since curvature is a unique property of curves, the feedrate-constraint value obtained by this method is not reasonable for Line Area.
In addition, some scholars obtain the feedrate-constraint value for continuous small line segments based on the limitations of acceleration of each axis [16][17][18] or contour error [19], but the calculation result will be unstable when the length of the line segment is too short. In recent years, some other scholars inserted a transition curve such as PH [20] and B-Spline [21][22][23] at the corner of adjacent line segments. Then, the feedrate can be constrained according to the curvature of the transition curve. At the same time, the machining path will be smoother. However, if the line segment is too short, the transition curve will be hard to insert. Fitting continuous small line segments to a curve such as Akima [24], B-Spline [25], and NURBS [26,27] is another similar method, and the method is also used by Siemens 840D CNC [28]. However, the trajectory error between the fitting curve and the original line segments is inevitable, which will damage machining quality. Additionally, there are many difficulties [29] in the feedrate planning and interpolation of the transition and fitting curve.
To address the above problems, this paper proposes a Feedrate-constraint method for the continuous small line segments based on the Nominal Acceleration (FNA). The FNA method obtains the predicted value of acceleration during machining process by the machining trajectory prediction and acceleration filtering, and the feedrate of the actual machining process is constrained according to the predicted acceleration. Compared with conventional methods, the FNA method is not sensitive to the angle of the corner and the length of the line segment, so it can be used for different machining paths consisting of long segments or continuous small segments, and it has no special requirements for the accuracy of the machining path. Therefore, the feedrate-constraint value obtained by this method is more reasonable and stable, which can improve the quality and efficiency of the CNC machining. In the following section, the FNA method will be introduced in detail and its effectiveness will be verified by simulation and experiment. This paper is organized as follows: Section 2 is the introduction of the FNA method; Section 3 is the simulation and experiment; and Section 4 is the conclusion.

The FNA Method
The principle of the FNA method proposed in this paper is as follows. For a certain corner in the machining path, whose feedrate-constraint value needs to be solved, assume that a constant feedrate F a is used for the CNC machining of the small line segments around this corner. Predict the actual machining trajectory based on the numerical analysis method of discrete systems. Then, obtain the predicted value of acceleration when the machining is carried out to the corner. Finally, calculate the feedrate-constraint value according to the predicted value of acceleration and its maximum allowable value.
Specifically, for a certain corner O of adjacent line segments whose feedrate-constraint value is to be solved, the whole process of the FNA method is shown in Figure 1, including (1) parsing the ISO6983 code and recording the endpoints V k of each line segment; (2) sampling the small line segments in the sampling interval containing the point O in Figure  1a to generate a position sampling sequence P i ; (3) predicting the corresponding machining trajectory of P i after the servo system, so as to generate the position sequence Q i , which is closer to the actual machining trajectory; (4) using the difference method to obtain the corresponding acceleration sequence A i from Q i , and then using a FIR low-pass filter to process A i , so as to obtain the nominal accelerationÂ corresponding to the sampling interval; and (5) obtaining the feedrate-constraint value according toÂ and the preset upper limit of normal acceleration a n . The main parts of the above process will be discussed in detail below.

The Generation of Position Sampling Sequence
In the flow of the FNA method shown in Figure 1, the machining trajectory predicting, the acceleration filtering, and the nominal acceleration solving, is based on discrete digital signals. Therefore, it is necessary to sample the continuous small line segments to generate a discrete position sampling sequence. The sampling process mainly consists of two steps: (1) Starting from the point O, the small line segments are traversed in the backward/forward direction, respectively, until the traversed distance reaches the preset target value l b and l f ; and (2) the traversed region of the machining path is taken as the sampling interval, and Appl. Sci. 2021, 11, 8837 4 of 22 N * points are sampled along the small line segments in this interval, thus generating the position sampling sequence P i : sampling interval, and * points are sampled along the small line segments in th val, thus generating the position sampling sequence : = ( ( ), ( ), ( )), = 0, ⋯ , * − 1 In Figure 1, the sampling interval of each part of the FNA method is a sym interval with point as the midpoint for easy understanding, but, in fact, the sa intervals of these parts are not the same and not all of them are symmetric, which d on the specific needs of each part.

The Prediction of Machining Trajectory
The position sequence shown in Equation (1) has been achieved accordin previous section. However, due to the trajectory characteristics of the small line se and the position error [30], will inevitably contain a random position error ∆ in turn, causes fluctuations in the velocity and acceleration as where s is the sampling period. If ∆ and s are taken to be 1 μm and 1 ms, tively, the acceleration fluctuation ∆ will be as high as 1000 mm/s 2 .
On the other hand, in the actual CNC machining process, will be transm the servo system to control the motion of the machine tool. Since the output of th system will lag behind the input [31], there will be a difference between the act chining trajectory and the command trajectory, which, in turn, will lead to a chang corresponding speed and acceleration.
For the above two reasons, this section will use the analysis method of the system to predict the "actual machining trajectory" corresponding to , and t ing of the nominal acceleration in a later section will be based on to make t feedrate-constraint value be more reasonable.
The function and control structure of the servo system of a typical three-ax machine tool is shown in Figure 2  In Figure 1, the sampling interval of each part of the FNA method is a symmetric interval with point O as the midpoint for easy understanding, but, in fact, the sampling intervals of these parts are not the same and not all of them are symmetric, which depends on the specific needs of each part.

The Prediction of Machining Trajectory
The position sequence P i shown in Equation (1) has been achieved according to the previous section. However, due to the trajectory characteristics of the small line segments and the position error [30], P i will inevitably contain a random position error ∆p, which, in turn, causes fluctuations in the velocity and acceleration as where T s is the sampling period. If ∆p and T s are taken to be 1 µm and 1 ms, respectively, the acceleration fluctuation ∆a will be as high as 1000 mm/s 2 .
On the other hand, in the actual CNC machining process, P i will be transmitted to the servo system to control the motion of the machine tool. Since the output of the servo system will lag behind the input [31], there will be a difference between the actual machining trajectory and the command trajectory, which, in turn, will lead to a change in the corresponding speed and acceleration.
For the above two reasons, this section will use the analysis method of the discrete system to predict the "actual machining trajectory" Q i corresponding to P i , and the solving of the nominal acceleration in a later section will be based on Q i to make the final feedrateconstraint value be more reasonable.
The function and control structure of the servo system of a typical three-axis CNC machine tool is shown in Figure 2, in which H x (s), H y (s), and H z (s) are the transfer functions of each axis. In order to solve the position sequence Q i of the servo system output, the difference equation corresponding to the transfer function of each axis should be obtained first, and then the command position sequence P i can be substituted into it to calculate Q i . Since the servo system of each axis usually has the same control structure, take the X-axis as an example and let its input and output position sequences be p x (k) and q x (k), respectively. Then, when H x (s) is approximated as a second-order system [32,33], its corresponding difference equation is the difference equation corresponding to the transfer function of each tained first, and then the command position sequence can be subs culate . Since the servo system of each axis usually has the same co the X-axis as an example and let its input and output position seque ( ), respectively. Then, when ( ) is approximated as a second-o its corresponding difference equation is The solution of the coefficients and initial conditions of this diffe be detailed below.

The Solution of Difference Equation Coefficients
For the coefficients 0 , 1 , 2 , 0 , and 1 in the difference Equatio output data of the X-axis servo system can be collected first, and the can be obtained through parameter identification. The specific proces Taking the interpolation period s as the sampling period, from the input and output position data of the servo system to obtain and the observation matrix Let the combined observation noise of the input and output sequ The solution of the coefficients and initial conditions of this difference equation will be detailed below.

The Solution of Difference Equation Coefficients
For the coefficients a 0 , a 1 , a 2 , b 0 , and b 1 in the difference Equation (3), the input and output data of the X-axis servo system can be collected first, and then those coefficients can be obtained through parameter identification. The specific process is as follows.
Taking the interpolation period T s as the sampling period, M points are sampled from the input and output position data of the servo system to obtain the output vector and the observation matrix Appl. Sci. 2021, 11, 8837 6 of 22 Let the combined observation noise of the input and output sequence be Then, substituting Equations (4)-(6) into (3), the following equation can be obtained: in which The ξ in the above formula is usually not 0 and unknown. In this case, the least square method can be used to calculate the optimal estimater of the coefficient vector r, which means that the sum of squares of the difference between the observed value q and its estimateq = Φr of the output sequence (Equation (8)) reaches the minimum.
Due to the fact that the derivative of a continuously derivable function reaches 0 at the extreme point, the following equation can be achieved: Thus, the optimal estimate of r iŝ So far, the coefficients of the difference Equation (3) have been obtained.

The Solution of Difference Equation Initial Conditions
For the difference Equation (3), solving its initial conditions is to obtain the values of p x (−2), p x (−1), q x (−2), and q x (−1). The values of p x (−2) and p x (−1) can be easily obtained by adding another two sampling points in Figure 1b, and for q x (−2) and q x (−1), according to the error relationship between the input and output of the servo system, the following equation can be obtained: where ε x (−2) and ε x (−1) are the tracking errors at p x (−2) and p x (−1), respectively. According to the control principle of the servo system [34], the tracking error ε x is proportional to the speed f x in the uniform motion state as follows: K x in the above formula is a scale factor, which can be deduced from Appendix A as For the more general non-uniform motion state, the tracking error is also approximately proportional to the speed by the factor K x . Therefore, if the X-axis speeds f x (−2), f x (−1) at p x (−2), and p x (−1) can be obtained, the approximations of ε x (−2) and ε x (−1) can be found; thus, the approximations of q x (−2) and q x (−1) can be obtained.
The values of f x (−2) and f x (−1) can be obtained by decomposing the speedF −2 andF −1 in Figure 3. Although the exact directions ofF −2 andF −1 are difficult to achieve, the direction ofF −2 can be approximately set to be the same as the vector (P −1 − P −3 ), considering the directional asymptotic property of the continuous small line segments. Therefore, the following equation can be obtained: Then, the approximation of (−2) can be achieved as Similarly, the approximation of (−1) can be obtained So far, the approximations of (−2), (−1) in the in ence Equation (3) have been achieved.

The Error Analysis and Control of the Approximations
In the previous section, an approximation to the initial Equation (3) was achieved and the error of the approximatio output sequence ( ) derived from Equation (3) and its tr decrease as increases.
Therefore, if the number of sampling points of the final p acceleration calculation in Figure 1d is set to and the num position sampling sequence ( ) in Equation (1) is set to nominal acceleration ̂ according to the last elements o between ̂ and its true value.
In order to reduce this error as much as possible, the va as possible, but too large a value will take up too much sto time. Therefore, the error of ( ) caused by the approximat tion needs to be analyzed to obtain an appropriate value of Let the errors between the approximations of (−2), be −2 and −1 , respectively. Then, it follows from the sup time-invariant system that the error in ( ) caused by the initial condition can be regarded as the zero input response o at initial conditions of (−2) = −2 , (−1) = −1 , namely, The X-axis component of the above equation iŝ Replacing theF −2 in the above equation with the feedrate F a , the value off x (−2) can be obtained asf Combining the above equation with Equation (12), the approximation of the tracking error at p x (−2) can be obtained as Then, the approximation of q x (−2) can be achieved aŝ Similarly, the approximation of q x (−1) can be obtained aŝ So far, the approximations of q x (−2), q x (−1) in the initial conditions of the difference Equation (3) have been achieved.

The Error Analysis and Control of the Approximations of Initial Conditions
In the previous section, an approximation to the initial conditions of the difference Equation (3) was achieved and the error of the approximation led to an error between the output sequence q x (i) derived from Equation (3) and its true value, which usually will decrease as i increases.
Therefore, if the number of sampling points of the final position sequence for nominal acceleration calculation in Figure 1d is set to N and the number of sampling points of the position sampling sequence q x (i) in Equation (1) is set to N * , when N * > N, solving the nominal accelerationÂ according to the last N elements of q x (i) can reduce the error betweenÂ and its true value.
In order to reduce this error as much as possible, the value of N * should be as large as possible, but too large a value will take up too much storage space and computation time. Therefore, the error of q x (i) caused by the approximation error of the initial condition needs to be analyzed to obtain an appropriate value of N * .
Let the errors between the approximations of q x (−2), q x (−1) and their true values be ε −2 and ε −1 , respectively. Then, it follows from the superposition of the linear and time-invariant system that the error in q x (i) caused by the approximation error of the initial condition can be regarded as the zero input response of the difference Equation (3) To obtain the expression of q x,zero (k), its corresponding characteristic root should be obtained as According to the value of ∆, the expression of q x,zero (k) can be divided into the following three cases: (1) When ∆ > 0, λ 1 , λ 2 are two unequal real numbers and the expression of q x,zero (k) is where C 1 , C 2 are unknown and they can be obtained by substituting q x,zero (0), q x,zero (1) into this equation. Thus, according to Equation (20), the following equation can be obtained: Then, substituting it into Equation (22), the values of C 1 , C 2 can be achieved as Since ε −2 , ε −1 in the formula are unknown, substituting C 1 , C 2 into Equation (22) cannot directly calculate the specific value of q x,zero (k), but its range can be achieved by inequation scaling. Scaling Equation (22) yields Similarly, the range of C 1 , C 2 can be obtained from Equations (24) and (A6) as At this point, let Then, the following inequation can be obtained from Equation (25): Therefore, in order to make the position error |q x,zero (k)| smaller than the preset upper limit ε max , the following inequation can be obtained: In this inequation, k is the only one unknown symbol. Therefore, the minimum value k min of k that satisfies the inequation can be achieved by traversing the integer (the process can be accelerated by using the bisection method).
(2) When ∆ = 0, λ 1 , λ 2 are two equal real numbers and the expression of q x,zero (k) is (3) When ∆ < 0, λ 1 , λ 2 are two unequal complex numbers and the expression of q x,zero (k) is For cases (2) and (3), the corresponding k min can be obtained according to the process of case (1).
So far, the number of sampling points N * of the Sequence (1) can be obtained as where N has been achieved from Equation (42) below. At this point, by substituting the p x (i) in Equation (1) into Equation (3), the output sequence q x (k) for the X-axis can be obtained. Similarly, the output sequence corresponding to the Y-axis and Z-axis can also be obtained. Then, the output position sequence of the servo system of the three-axis CNC machine tool can be achieved as In order to reduce the impact on the subsequent calculation, which is caused by the approximation error of the initial conditions, the last N terms of the above sequence are taken out and renumbered as Compared with P i in Equation (1), Q i achieved here is closer to the actual machining trajectory. At this time, the target values l b and l f of the traversed distance in Section 2.1 can be achieved as where L is the preset length of the Sequence (34).

The Solution of the Nominal Acceleration
According to the position sequence Q i obtained above, the nominal acceleration can be solved. Taking the X-axis as an example, the first-order difference calculation of the q x (i) in Q i can be performed to obtain the feedrate sequence as . . .
Similarly, the acceleration sequence a x (i) can be obtained from Equation (36) as The process of calculating the acceleration sequence a x (i) is shown in Figure 4a. The midpoint of the sequence is the X-axis componentÂ x of the acceleration at point O. However, since the above differential calculation only uses the position data of a few adjacent interpolation cycles around point O, the value ofÂ x may fluctuate greatly, which will affect the stability of the final nominal acceleration. For this reason, a FIR low-pass filter will be used, as shown in Figure 4b, to process the acceleration sequence a x (i) of the sampling interval, so as to reduce the fluctuation and obtain a more stable and reasonable nominal acceleration.
The process of calculating the acceleration sequence ( ) is shown in Figure 4a. The midpoint of the sequence is the X-axis component ̂ of the acceleration at point . However, since the above differential calculation only uses the position data of a few adjacent interpolation cycles around point , the value of ̂ may fluctuate greatly, which will affect the stability of the final nominal acceleration. For this reason, a FIR low-pass filter will be used, as shown in Figure 4b, to process the acceleration sequence ( ) of the sampling interval, so as to reduce the fluctuation and obtain a more stable and reasonable nominal acceleration. First, the design of the FIR low-pass filter needs to be carried out, which means obtaining the value of h(i) in Figure 4b. Let the frequency characteristics of the ideal linear-phase FIR low-pass filter be where j is an imaginary unit, ω c is the cutoff frequency, and τ is the group delay. By the Fourier inverse transform of H d e jω , the unit impulse response of the filter can be obtained as Since i in the formula can be taken as any integer, h d (i) is an infinitely long noncausal sequence, which cannot be used in the actual calculation. In this case, a window function [35,36] needs to be used to truncate h d (i). From the data in Table 1, the stopband attenuation δ stop in this paper is −40 dB. Therefore, the Hanning window function with the stopband attenuation −44 dB in Appendix B can meet the requirements, and its expression is N in the above formula can be obtained from the transition bandwidth formula of the Hanning window function, and its value is 6.2π where f s , f pass , and f stop are the sampling frequency, passband cutoff frequency, and stopband start frequency of the FIR low-pass filter, respectively. To ensure that Equation (34) is a position sampling sequence with O as the midpoint, N must be an odd number; therefore, the N obtained from Equation (41) needs to be corrected as The unknown τ in the formula can be obtained from the linear-phase characteristic of the FIR filter as At this point, in order to check whether the amplitude-frequency characteristics of the designed FIR low-pass filter meet the design requirements, the discrete Fourier transform of h(i) can be performed as Then amplitude-frequency response at f stop is obtained as If δ is less than the δ stop in Table 1, the designed filter meets the requirement; otherwise, the value of N should be increased or another type of window function should be used.
So far, the design of the FIR low-pass filter was completed and the filter coefficients h(i) shown in Equation (43) were obtained. By filtering the acceleration Sequence (37), the X-axis acceleration at point O can be obtained aŝ Denoting the filter coefficients, h(i) are written in vector form: The value ofÂ x can be obtained from Equation (37) and (47) aŝ Similarly, the acceleration of the Y-axis and Z-axis at the point O can be achieved as Therefore, the nominal acceleration at point O can be obtained aŝ This is the nominal acceleration of the sampling interval with O as the midpoint in Figure 1b.

The Solution of the Feedrate-Constraint Value
In the process of solving the nominal accelerationÂ, shown in Equation (51), the feedrate remains constant and the tangential acceleration is 0. Therefore, the nominal accelerationÂ is the normal acceleration, so that the nominal radius of curvature of the continuous small line segments in the sampling interval can be calculated as At this point, assume that the small line segments of the sampling interval are machined at another feedrate F O and the normal acceleration generated under the action of the nominal radius of curvature ρ v is just equal to the maximum allowable value of the normal acceleration a n , namely, Then, F O is the feedrate-constraint value at point O and its value is So far, the feedrate-constraint value at any corner of the continuous small line segments can be obtained from the proposed FNA method. Then, combined with the conventional feedrate planning [18,37,38] and interpolation methods, the CNC machining of the continuous small line segments can be realized.

Simulation and Experimental Results
In order to verify the effectiveness of the FNA method proposed in this paper, the algorithm simulation and machining experiment will be conducted in this section. The algorithm simulation was performed on a personal computer with the CPU i3-9100F, while the machining experiment was performed on a three-axis CNC machine tool DC-6060A (KAIBO, Ningbo, China) controlled by a self-developed controller. The X/Y/Z axes travel of the machine tool were 600/600/260 mm, respectively. The architecture of the experimental platform is shown in Figure 5.
The Embedded Linux controller iTop4412 (TOPEET, Beijing, China) in the figure was mainly used for GUI, while the ARM (STM32H743) controller ran CNC algorithms, including ISO6983 code parsing, nominal acceleration solving, feedrate planning, interpolation, etc. The FPGA (EP4CE22) was mainly used for fine interpolation and generating control signals for servo drives.
The main parameters used in the simulation and experiment and their default values are shown in Table 1.
In the actual CNC machining, the angle-constraint method is easy to realize and simple to calculate and it is the most widely used method to calculate the feedrate-constraint value. Therefore, the simulation and experiment were mainly performed on the angle-constraint method and the proposed FNA method.

2D Machining Simulation
The 2D machining path shown in Figure 6 contains corners, arcs, and reversals, which are typical and complex. Therefore, it can be used for the simulation to compare the effect of the angle-constraint method and FNA method. algorithm simulation and machining experiment will be conducted in this section. T algorithm simulation was performed on a personal computer with the CPU i3-910 while the machining experiment was performed on a three-axis CNC machine tool D 6060A (KAIBO, Ningbo, China) controlled by a self-developed controller. The X/Y/Z ax travel of the machine tool were 600/600/260 mm, respectively. The architecture of the e perimental platform is shown in Figure 5. The Embedded Linux controller iTop4412 (TOPEET, Beijing, China) in th mainly used for GUI, while the ARM (STM32H743) controller ran CNC alg cluding ISO6983 code parsing, nominal acceleration solving, feedrate planni lation, etc. The FPGA (EP4CE22) was mainly used for fine interpolation and control signals for servo drives.
The main parameters used in the simulation and experiment and their de are shown in Table 1.
In the actual CNC machining, the angle-constraint method is easy to real ple to calculate and it is the most widely used method to calculate the feedrat value. Therefore, the simulation and experiment were mainly performed o constraint method and the proposed FNA method.

2D Machining Simulation
The 2D machining path shown in Figure 6 contains corners, arcs, an which are typical and complex. Therefore, it can be used for the simulation the effect of the angle-constraint method and FNA method.  The radius of curvature ρ at Area1 the in above figure is 3 mm, and the curve here is discretized into small continuous line segments with the chord error ε = 3 µm. Therefore, the angle of the corner of adjacent segments is Therefore, the feedrate-constraint value obtained by using the virtual transition-arc constraint method [11] is while, according to the constraint of the radius of curvature and the normal acceleration, the feedrate-constraint value should be Similarly, the feedrate-constraint values at Area2, Area3, and Area4 can also be achieved, as shown in Table 2. The angle-constraint method and the FNA method were used to calculate the feedrateconstraint value of the machining path in Figure 6, respectively. Then, the feedrate planning and interpolation were performed, and the results are shown in Figure 7.
As can be seen from Figure 7a, for Line Area in Figure 6, both the angle-constraint method and the FNA method effectively constrained the feedrate, while for Curve Area, the feedrate-constraint values obtained from these two methods were significantly different. The machining paths in Area1 and Area2 (or Area3 and Area4) had the same radius of curvature, but the angles were different. Therefore, if the angle-constraint method is used to calculate the feedrate-constraint value, the machining paths with the same radius of curvature will be machined by different feedrates, as shown in Figure 7b. Especially in Area1 and Area4, the feedrate exceeded the theoretical value obtained by Equation (58). While in Area1-Area4 of Figure 7c, the feedrates obtained based on the FNA method were all nearly same as the theoretical values in Table 2.
In addition, the small line segments in Area2 and Area3 were relatively longer; so, the velocity direction changed more at the corner. As a result, in Figure 7a, the feedrateconstraint value obtained by the FNA 1 method was obviously lower than the theoretical value, while the feedrate-constraint value obtained by the FNA method was almost equal to the theoretical value, which verified the effectiveness of the machining trajectory prediction in Section 2.2. Therefore, the machining trajectory prediction can suppress the accuracy problem of continuous small line segments and make the calculating results of the proposed method more reasonable and stable. Similarly, the feedrate-constraint values obtained by the FNA 2 method in Area1 and Area2 (or Area3 and Area4) had a difference over 50 mm/min, while the difference of the FNA method was less than 10 mm/min, which verified the effectiveness of the acceleration filtering in Section 2.3.
In summary, the proposed FNA method was not sensitive to the length of the line segment, the angle of the corner, etc. The proposed method was stable and the obtained feedrate-constraint value was more reasonable.

Area
Area1 Area2 Area3 Area4 Radius of the curvature (mm) The angle-constraint method and the FNA method were used to calculate the feedrate-constraint value of the machining path in Figure 6, respectively. Then, the feedrate planning and interpolation were performed, and the results are shown in Figure  7. As can be seen from Figure 7a, for Line Area in Figure 6, both the angle-constraint method and the FNA method effectively constrained the feedrate, while for Curve Area, the feedrate-constraint values obtained from these two methods were significantly different. The machining paths in Area1 and Area2 (or Area3 and Area4) had the same radius

3D Machining Simulation
The effectiveness of the FNA method was preliminarily verified by the 2D simulation in the previous section. In this section, the FNA method will be further simulated and analyzed by using the continuous small line segments of a 3D free-form surface.
For the first 3000 lines of the ISO6983 code [39] corresponding to the machining path shown in Figure 8b, the angle-constraint method, the curvature-constraint method, the acceleration-constraint method [17], and the proposed FNA method were used to calculate the feedrate-constraint value, respectively. Then, the feedrate planning and interpolation were performed, and the final feedrate results are shown in Figure 9.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 17 of 23 of curvature, but the angles were different. Therefore, if the angle-constraint method is used to calculate the feedrate-constraint value, the machining paths with the same radius of curvature will be machined by different feedrates, as shown in Figure 7b. Especially in Area1 and Area4, the feedrate exceeded the theoretical value obtained by Equation (58). While in Area1-Area4 of Figure 7c, the feedrates obtained based on the FNA method were all nearly same as the theoretical values in Table 2.
In addition, the small line segments in Area2 and Area3 were relatively longer; so, the velocity direction changed more at the corner. As a result, in Figure 7a, the feedrateconstraint value obtained by the FNA 1 method was obviously lower than the theoretical value, while the feedrate-constraint value obtained by the FNA method was almost equal to the theoretical value, which verified the effectiveness of the machining trajectory prediction in Section 2.2. Therefore, the machining trajectory prediction can suppress the accuracy problem of continuous small line segments and make the calculating results of the proposed method more reasonable and stable. Similarly, the feedrate-constraint values obtained by the FNA 2 method in Area1 and Area2 (or Area3 and Area4) had a difference over 50 mm/min, while the difference of the FNA method was less than 10 mm/min, which verified the effectiveness of the acceleration filtering in Section 2.3.
In summary, the proposed FNA method was not sensitive to the length of the line segment, the angle of the corner, etc. The proposed method was stable and the obtained feedrate-constraint value was more reasonable.

3D Machining Simulation
The effectiveness of the FNA method was preliminarily verified by the 2D simulation in the previous section. In this section, the FNA method will be further simulated and analyzed by using the continuous small line segments of a 3D free-form surface.
For the first 3000 lines of the ISO6983 code [39] corresponding to the machining path shown in Figure 8b, the angle-constraint method, the curvature-constraint method, the acceleration-constraint method [17], and the proposed FNA method were used to calculate the feedrate-constraint value, respectively. Then, the feedrate planning and interpolation were performed, and the final feedrate results are shown in Figure 9.  For area 1 with large curvature in Figure 9a, the distribution of the feedrate in the X-Y plane obtained by the angle-constraint method is shown in Figure 9b, and the feedrate was distributed from 500 mm/min to 3000 mm/min, which is rather chaotic. In a contrast, the feedrate obtained from the FNA method is shown in Figure 9e, which was maintained around 1500 mm/min, and the stability of the feedrate was obviously better than the one of the angle-constraint method and the other two constraint methods.
In order to quantify the stability of the feedrate to make a better comparison between different methods, the feedrate values for each row at X = 43 mm are collected in Figure  9b-e, respectively, and their distribution is shown in Figure 9f. The average, range, relative range (range/average), and standard deviation were calculated (the first three points in the circles of Figure 9b-f are in the reverse region and were not taken into the statistics), and the results are shown in Table 3.  For area 1 with large curvature in Figure 9a, the distribution of the feedrate in the X-Y plane obtained by the angle-constraint method is shown in Figure 9b, and the feedrate was distributed from 500 mm/min to 3000 mm/min, which is rather chaotic. In a contrast, the feedrate obtained from the FNA method is shown in Figure 9e, which was maintained around 1500 mm/min, and the stability of the feedrate was obviously better than the one of the angle-constraint method and the other two constraint methods.
In order to quantify the stability of the feedrate to make a better comparison between different methods, the feedrate values for each row at X = 43 mm are collected in Figure  9b-e, respectively, and their distribution is shown in Figure 9f. The average, range, relative range (range/average), and standard deviation were calculated (the first three points in the circles of Figure 9b-f are in the reverse region and were not taken into the statistics), and the results are shown in Table 3. From the data in this table, it can be seen that the proposed FNA method worked much better than other three methods. Especially compared with the angle-constraint method, the FNA method reduced the relative range of the feedrate from 90.8% to 4.3%, and the standard deviation also decreased from 602 to 17. Therefore, the feedrate obtained by the FNA method had better stability, which is quite important for improving the CNC machining quality of the free-form surface.

Experiment
The numerical analysis in the simulations verified the effectiveness of the FNA method, while this section will use the experiment to verify the effectiveness of the FNA method more directly. For the machining path shown in Figure 8b, the feedrate was constrained by the angle-constraint method and FNA method, respectively, and the CNC machining was performed under the same machining conditions. In the machining, a ball milling cutter W2QX4408 (Weidiao, Shanghai, China) was used, and the machining material was AlMg1SiCu aluminum. The command feedrates were both 3000 mm/min, and the machining results are shown in Figure 10.
From the data in this table, it can be seen that the proposed FNA m much better than other three methods. Especially compared with the an method, the FNA method reduced the relative range of the feedrate from and the standard deviation also decreased from 602 to 17. Therefore, the fee by the FNA method had better stability, which is quite important for impro machining quality of the free-form surface.

Experiment
The numerical analysis in the simulations verified the effectivenes method, while this section will use the experiment to verify the effectivene method more directly. For the machining path shown in Figure 8b, the fee strained by the angle-constraint method and FNA method, respectively, and chining was performed under the same machining conditions. In the ma milling cutter W2QX4408 (Weidiao, Shanghai, China) was used, and the m terial was AlMg1SiCu aluminum. The command feedrates were both 3000 the machining results are shown in Figure 10.

By Angle
By FNA Figure 10. The machining results of the angle-constraint method and FNA method When the angle-constraint method was used, the obtained feedrate and unstable and the overcutting occurred at many places, resulting in mach represented by "pits". However, when the FNA method was used, the obt When the angle-constraint method was used, the obtained feedrate was overlarge and unstable and the overcutting occurred at many places, resulting in machining damage represented by "pits". However, when the FNA method was used, the obtained feedrate was reasonable and stable and the machining quality was obviously improved. In addition, the machining times of the two methods were 22 min 45 s and 22 min 8 s, respectively. Therefore, the proposed FNA method can both improve the machining efficiency and the machining quality. The machining experiments and algorithm simulations corroborate each other and they both verify the effectiveness of the FNA method proposed in this paper.

Conclusions
The conventional feedrate-constraint method is sensitive to the angle of the corner and the length of the line segment, and the feedrate-constraint value obtained is unreasonable, which will result in damages to the machining quality. In this paper, the FNA method was proposed, in which the line segment path sampling, the machining trajectory prediction, and the acceleration filtering were used, and they can significantly improve the reasonableness and stability of the feedrate-constraint value.
In the simulation of the 2D machining path, the feedrate-constraint value obtained by the FNA method was more stable and reasonable compared with the angle-constraint method, and was almost equal to the theoretical value. In the simulation of the 3D machining path of the free-form surface, compared with the angle-constraint method, the relative range of the feedrate obtained by the FNA method decreased from 90.8% to 4.3% and the standard deviation decreased from 602 to 17, which means that the stability of the feedrate was significantly improved. Compared with the curvature-constraint method and acceleration-constraint method, the proposed method also worked much better. The machining experiments showed that the proposed FNA method ran stably and can both reduce the machining time and improve the machining quality, thus helping to promote the application of the continuous small line segments in the high-speed and high-precision CNC machining.
In view of the limited time, the FNA method only worked for line segments. It will be improved to be applicable for the machining path consisting of the line segment, the arc, and the parametric curve. Additionally, the current method was only implemented in the three-axis CNC machining and has not yet been applied to the five-axis, which will be realized in the future.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.

FNA
The proposed Feedrate-constraint method based on Nominal Acceleration Line Area Area with long line segments in the machining path Curve Area Area with continuous small line segments discretized from the curve in the machining path F o The feedrate-constraint value at the corner of adjacent line segments F a The const feedrate used for machining the small line segments in FNA method P i Position sampling sequence of the sampling interval N * The length of the sequence P i Q i Predicted position sequence corresponding to P i N The length of the sequence Q î A Nominal acceleration corresponding to a corner T s Interpolation period and sampling period a n Normal acceleration in which 1+b 0 +b 1 T s It follows that when the X-axis moves at a constant speed, the tracking error ε x is proportional to the speed f x and the scale factor is K x .
Appendix B Table A1. Several commonly used window functions in the design of the FIR filter [35,36].