A Novel Local Smoothing Method for Five-Axis Machining With Time-Synchronization Feedrate Scheduling

This paper presents an analytical continuous smoothing method for the five-axis toolpath by simultaneously scheduling the tool position and tool orientation trajectories. In order to ensure the high-order continuous, the peak-controlled jerk and arclength-parameterized property, a novel curve “airthoid” is proposed for the first time. The biairthoid is involved to smooth the corners of the tool position in the workpiece coordinate system (WCS) and the corners of the tool orientation in the machine coordinate system (MCS), the geometries of which are analytically determined by the user-defined deviation errors. A time synchronization strategy is proposed to extend the duration of the predetermined cubic acceleration profile to a specified time. With the kinematic constraints of the tool position and the tool orientation, the transitional and rotational trajectories are analytically synchronized by sharing the same motion time. To comply with the constraints of the linear feed drives, an optimization strategy is conducted by adjusting the kinematics of the tool position. By doing so, the approximation errors of the tool position and tool orientation in the WCS are strictly satisfied. The analytical arclength expression of the smoothing curves is more suitable for the on-line interpolation. Due to the arclength-parametrized transition curve, the feedrate fluctuation is eliminated. With the proposed time synchronization strategy, the physical limits of the feed drives are all respected. Moreover, the high-order continuous airthoid makes the motion more smoothing-going. Simulations and experiments verify the effectiveness of the proposed algorithm.


I. INTRODUCTION
Five-axis machine tools are widely adopted to machine free-form parts due to the graceful advantages in the cutting efficiency and machining reachability. The parametric curves, such as polynomial splines, B-spline, Bézier or NURBS, are ideal approaches to describe the toolpaths of the free-form parts, which have been experimentally proven to be superior to achieve smooth motion and perfect surface finish [1]. However, most of CNC systems are hard to follow the parametric splines directly, since some bottlenecks still exits, including the poor interpretation ability from CAD to NC [2], The associate editor coordinating the review of this manuscript and approving it for publication was Min Xia . the complex calculation of path lengths [3] and the unwanted feedrate fluctuations [4]. Instead, the parametric curve is usually discretized as the linear segments under the predefined tolerance, and then, segments are programmed as G01 commands [5]. Although the successive G01 commands can be easily followed by machine tools, there are also problems needed to be solved. The tangent vectors at the junction of line segments are different. This enforces the feed drives to stop at each segment junction, owing to the limited motor driving ability. As a result, the machining cycle time increases [6] and the surface finish quality decreases [7].
To overcome these problems, the corners of the multi-axis toolpath are smoothed by inserting high-order micro-splines to achieve the non-stop motion, and then, the feedrate VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ scheduling strategy is involved to generate the feedrate profile with the bounded kinematic constraints. The 5-axis machining toolpath is composed of two independent subpaths, i.e. the tool position and tool orientation [8], [9]. Researchers usually smooth 5-axis linear paths by inserting two different parametric splines to attain the geometric high-order continuity, and at the same time, to guarantee the subspline parameters synchronization. This is achieved by describing the curved 5-axis toolpath with parametric splines. Beudaert et al. [10] developed two cubic B-Spline curves to smooth the shape corners in the workpiece coordinate system (WCS), and then, a third B-Spline curve was designed to synchronize the spline parameters. Tulsyan and Altintas [11] inserted the synchronized quantic and septic curves for the tool tip position and the tool orientation in the WCS. Bi et al. [12] proposed an analytical corner smoothing method with dual cubic Bezier curves in the machine coordinate system (MCS), which synchronously blended the translational and rotational paths by adjusting the equivalent ratios between the transition lengths and the control polygon lengths. Yang and Yuen [13] developed an quintic B-spline to smooth the tool position corner in the WCS and the tool orientation corner in the MCS. The parameter synchronization was analytically guaranteed by ensuring that the third-order derivation of the tool orientation with respect to the tip position displacement was continuous. Huang et al. [14] involved two cubic B-splines to blend the toolpaths in the WCS, the parameters of which were synchronized by converting the remaining linear segments into the parametric form. Yang et al. [15] proposed a smoothing method with PH splines to blend the tool position corner in the WCS and the tool orientation corner in the MCS. The synchronization was analytically reached by converting the remaining linear segments into the B-spline. Following geometric smoothing, the feedrate profile of the tool tip path was determined by the transitional kinematic constraints. The feedrate profile of the orientation path was obtained by synchronizing the parameter of the tool orientation with tool position. It should be noticed that in order to satisfy the real-time requirement of the CNC system, the advanced local smoothing algorithm with the geometric smoothing and feedrate scheduling stages should be implemented within fewer timestamps. Although the parametric curve, such as high-order polynomials [16], [17], B-splines [11], [13], [14], [18]- [20], Bézier splines [10], [12], [21], [22], can be analytically determined in accordance with the specific geometric constraints, there exists two main challenges. First of all, the arclength has no analytical expression with spline parameters, and thus, the estimation of the toolpath length is based on iterative numerical algorithms [18], [23]. It is computationally expensive for real-time interpolation, especially when the splines are inserted to smooth the successive short linear segments or the arclength occupies the vast majority of the toolpath. In addition, the feedrate fluctuation can be easily induced [24], [25], due to inaccurate mapping between the targeted arc displacement and the updated spline parameter. It degrades the machining quality and causes motor torque saturation. The feedrate correction strategies [4], [26], [27] are usually involved to eliminate the fluctuation through remapping the spline parameters to the arclength. Pythagorean-hodograph (PH) curve [28] is a subset of the Bézier spline, but the arclength is the polynomial expression of the spline parameter. It was successfully introduced to blend and interpolate the successive linear toolpaths [6], [15], [29], which did greatly improve the calculation efficiency with the analytical expression of the toolpath length in the lookahead stage. It is noted that, as inherent of non-arclength parameterized splines, the feedrate correction model should be repeatedly executed to revise every interpolation position in the fine interpolation stage [6]. It is still computationally stringent for the real-time interpolation. For these reasons, to realize real-time interpolation and achieve smooth motion, it is important to involve other splines, which possesses the analytical expression of the arclength and the spline parameter directly based on the arclength.
Recently, clothoid draws more attentions for local corner smoothing [30], [31], since it is the arclength-parameterized spline. It provides an analytical expression of the arclength. Moreover, the feedrate fluctuation can be circumvented without any feedrate correction. It is more friendly for the on-line smoothing since more computing resources can be released to handle other tasks [6], [32]. However, the C 2 continuity leads to the discontinues jerk at the junction of the linear and curvilinear sub-paths and the junction of biclothoid [31]. The jerk jump at the junction of the biclothoid can be confined by calculating the critical velocity in the lookahead stage, whereas the jerk jump at the junction of linear and curvilinear splines is hard to control. As a result, it leads to the uncontrollable jerk jump at the junctions, and thus, easily violates the predefined the jerk limitation of the feed drive. In fact, the jerk of the velocity profile of the feed drive is subject to the responsiveness of the motor. Due to the finite responsiveness ability, the jerk peak value should be confined. Otherwise, the exaggerated and discontinuous jerk would excite the structural vibration and decrease the machining surface, as experimentally verified by Yang and Yuen [13] and Fan et al. [22]. Therefore, the arclength parameterized spline with the controllable jerk profile is essential to achieve smoother motion, and thus, implement the real-time interpolation. To the best of the authors' knowledge, the spline has not been analytically applied to on-line blend the corners of the multi-axis machining toolpaths.
In the previous researches [10]- [14], the feedrate profile of the tool position is established on the basis of the transitional kinematic constraints, and the tool orientation path was treated as the slave motion scheduled by the parameter synchronization. Guided by this fact, the transitional kinematic constraints are well bounded, whereas the kinematic constraints of the tool orientation are left out of consideration. As a result, the abrupt change of the rotational motion might occur [33], [34], when synchronizing the parameter of the orientation toolpath to the tool position.
Then, the torque limitation or the maximum speed of the linear and rotary drives are violated. Moreover, due to the nonlinear relationship between the parameter and arclength of the non-arclength parameterized splines, the feedrate fluctuation of the tool orientation may happen evenly with the constant increment of the parameter generated by the tool position [24], [25]. As a consequence, the control stability of the machine tool is destroyed and the visible poor machine quality may be confronted. Therefore, the smoothness of the rotational motion should be carefully conducted.
Some researches attempted to treat the angular (i.e. tool orientation) motion equally with the linear (i.e. tool position) motion, and simultaneously plan the transitional and rotational trajectories respecting the linear and angular motion constraints. Huang et al. [35] proposed a time synchronization method for the tangent jerk limited velocity profile, and then the linear and angular kinematic constraints were both successfully constrained. Liu et al. [9] smoothed the linear and angular motions based on the finite impulse response filters technology, and generated two tangent jerk limited trajectories with the same motion durations. By doing so, the corresponding kinematic constraints were successfully constrained with the smooth transition of local corners. It is noted that the tangent jerk profile in Ref. [9], [35] is bounded but discontinuous. Actually, the tangent jerk reflects the smoothness of the path velocity profile. As experimentally verified by Yang and Yuen [13] and Fan et al. [36], the path velocity profile with the continuous tangent jerk can obtain better machining quality and less machine vibrates than the path velocity profile with the limited but discontinuous jerk. Therefore, it is an important task to develop a path velocity scheduling strategy with the continuous tangent jerk based on the time synchronization for the 5-axis toolpath interpolation.
To simultaneously obtain the real-time high-order continuous interpolation and respect both the kinematic constraints of the tool position and tool orientation paths, this paper proposes a novel continuous smoothing method based on the airthoid curves for the five-axis machining toolpath. First, a C 3 continuous arclength-parameterized curve (airthoid) is designed. Second, the corner of the tool tip position in the WCS and the corner of the tool orientation path in the MCS are analytically blended. The geometries of the smoothing curves are determined by the predefined approximation errors and the linear segment lengths. Third, the transitional and rotational motions with the continuous jerk are analytically synchronized by sharing the same motion time and respecting the linear and angular kinematic constraints. The linear and angular feedrate profiles based on S-curve-type acceleration are separately planned to respect the linear and angular kinematic constraints. Then, the time synchronization strategy extends the duration of the faster sub-path to realize the motion synchronization.
The present article is organized as follows. Characteristics of the proposed airthoid spline are presented in section II. Then, the smoothing strategy for corners of the tool position and the tool orientation is demonstrated in section III.
The time synchronization is presented in section IV to extend the duration of the S-curve-type acceleration profile to any value. Furthermore, a look-ahead corner smoothing method for the tool position and tool orientation trajectories is developed in section V. The simulations and experiments are presented in sections VI and VII, respectively. Finally, the article is summarized and concluded in section VIII.

II. AIRTHOID SPLINE
The clothoid spline is an arclength-parameterized spline, where the corresponding curvature increases linearly with the arclength. The C 2 continuity with the uncontrolable and discontinues jerk leads to motion impact at the junction of the linear and curvilinear paths. To achieve smoother and jerk-controlable motion profile, the C 3 continuous arclength-parameterized spline is essential. It can be realized by designing a new arclength-parameterized spline, named airthoid. It delivers the derivation of the curvature increases linearly with the arclength. The curvature of the airthoid spline is expressed as where κ 0 and s are the initial curvature and the arclength, respectively. Moreover, c and θ denote the sharpness of the airthoid spline and the tangent angle, respectively. The tangent angle θ of the spiral can be calculated through the integration of the curvature respect to the arclength s.
where θ 0 is the tangent angle at the start point. It should be indicated that θ 0 and κ 0 are zeros, since the spline is inserted as the smoothing curve to connect the linear segment.
The tangent angle of the airthoid is meaningful to design the geometry of the smoothing curve, because it can be connected with the corner angle. Defining θ = c 3 t 3 , Eq.(3) can be rewritten as where a is the scaling parameter in the form below: and the terms C(θ) and S(θ) are shown as The arclength P 0 P is achieved as It is noted that the arclength of the airthoid can be determined with θ, when the geometry is predefined with a.
The curvature at P is the reciprocal of the curvature radius, obtained as Theorem 1: The path of the biairthoid at the junction of the linear and curvilinear splines is G 3 continuous. The proof of theorem 1 is shown in Appendix A.
Theorem 2: The path of the biairthoid at the junction of the curvilinear splines is G 2 continuous, but with the peak-controlled jerk.
The proof of theorem 2 is shown in Appendix B.

III. FIVE-AXIS TOOLPATH SMOOTHING A. THE ANALYTICAL SMOOTHING FOR THE CORNER OF THE TOOL POSITION
The relative velocity between the tool tip and the workpiece should be strictly controlled in the machining process, since it directly affects the finished surface quality [1]. Therefore, the corners of the tool position are blended and the feedrate profile of the smoothing toolpath is planned in the WCS to guarantee the objective velocity of the tool tip. To achieve the high order continuous motion, a biairthoid spline B 0 p B 1 p is constructed to smooth the adjoining lines, P 1 P 2 and P 2 P 3 , in the WCS, which consists of back to back symmetrical airthoids, as shown in Fig. 2. The airthoid pair is connected at the point Q, which is defined as the critical point. The tangent angle at Q is defined as α, which is equal to 0.5(π −ϕ). ϕ is the corner angle. The length of the transition segments, B 0 p P 2 and B 1 p P 2 , is denoted as l p . The lengths of the smoothing curve, B 0 p Q and QB 1 p , are denoted as l 0 Bp and l 1 Bp , respectively.

1) MAXIMUM DEVIATION ERROR
Due to the symmetrical airthoids, the maximum deviation error ε pw generates at the critical point Q, which can be analytically expressed as It should be indicated that the deviation error ε pw is only determined by the transition length l p under the predefined the corner angle.
With the predefined deviation error, the intersection of the adjoining smoothing curves may occur for the limited length of the linear segments. Therefore, the maximum transition length is restricted to less than half length of the original linear segment, since each linear segment supports the smoothing of two neighboring corners. With the geometric constraints, the transition length should be determined as Based on the constraints of the predefined contour deviation error and the limited length of the linear segments, the geometry of the airthoid is determined through the scaling parameter a p , which is defined as:

2) LENGTH OF THE SMOOTHING CURVE
The length of the smoothing curve can be analytically calculated through combining Eq.(8) and Eq.(12) in the form below:

3) MAXIMUM CURVATURE
The curvature extremum of the smoothing curve occurs at the critical point, which can be analytically calculated as The specific derivation process is provided in the Appendix C.

B. THE ANALYTICAL CORNER SMOOTHING OF THE TOOL ORIENTATION
To plan the rotational trajectory, the length of the orientation toolpath should be obtained. Fig. 3(a) shows that when the smoothing of the orientation toolpath is accomplished in the WCS, the smoothing curve should be normalized to lie on the surface of the unity sphere [11], [14]. The normalization deteriorates the on-line application, since the length of the orientation toolpath should be iteratively calculated. Instead, the problem can be overcome by inserting the biairthoid to smooth the linear segments of rotary commands in the MCS, as shown in Fig 3(b). As a result, the analytical calculation of the smoothing toolpath length can be both realized. The deviation error of the tool orientation ε ow should be strictly constrained to avoid exceeding the predefined machining error caused by the slope of the tool orientation. It can be realized by confining the deviation error ε om between the smoothing curve and the linear segments of rotary commands proposed by Liu et al. [9].
Similar to the transition length of the tool position in the WCS, the needed transition length l omε of the rotational toolpath in the MCS is where β and ε om denote the tangent angle with the linear rotary command and the predefined orientation error tolerance in the MCS, respectively. To eliminate the intersection of the smoothing curves, the maximum transition length l o should be no more than the half length of the linear segment.
It is noted that the geometry of smoothing curve based on the airthoid is confirmed with the scaling parameter a o .
The length of the smoothing curve in the machine coordinate system is analytically calculated as

IV. TIME SYNCHRONIZATION ALGORITHM BASED ON THE S-CURVE-TYPE ACCELERATION PROFILE A. VELOCITY PROFILE BASED ON THE S-CURVE-TYPE ACCELERATION
The jerk-continuous feedrate profile performs better than the jerk-limited feedrate profile in reducing the machine tool vibration and increasing the maching accuracy, which has been experimentally verified by Refs. [13], [37]. The S-curve-type acceleration profile is widely utilized to generate the jerk-continuous motion [6], [11], [38]. Fig. 4 illustrates the motion profiles corresponding to the S-curve-type acceleration. It indicates that, with the given start velocity v s , end velocity v e and the displacement L, the trajectory is determined to achieve the shortest motion time with regard to the jerk and acceleration limitations.
where V a and V d are v m − v s and v e − v m , respectively. t a , t d and t c are the durations of the acceleration, deceleration and the constant velocity stages, respectively. Thus, the transition time of the acceleration and deceleration stages FIGURE 4. Motion profiles based on the S-curve-type acceleration [6]. VOLUME 8, 2020 can be determined through integration over the acceleration limitation A max and the jerk limitation J max , which can be expressed as The distance in the acceleration, deceleration and the constant velocity stages denotes as L a , L d and L c , which can be obtained from the following term: Detailed discussion of the above material is available in Ref. [38].

B. TIME SYNCHRONIZATION
The time synchronization method proposed by Huang et.al. [35] focuses on extending the motion duration of the jerk limited acceleration profile to any value. To confront jerk continuous motion, a new time synchronization strategy is proposed to extend the duration of the S-shape-acceleration motion to any specified time. Fig. 5 shows that the five synchronization types A i (i = 1, 2, 3, 4, 5) are involved to extend the trajectory to the goal duration T e . It indicates that synchronization types can be divided into the following categories: (i) the acceleration, uniform and deceleration motions; (ii) the acceleration and uniform motions; (iii) only the acceleration motion; (iv) the deceleration and uniform motions; (v) the deceleration, uniform and deceleration motions, respectively.

1) DETERMINING THE DURATION THRESHOLD
Three duration thresholds T i (i = 1, 2, 3) are involved to determine the revised motion profile types with the given the goal duration T e . The duration threshold T 1 is defined as where T 1a is the duration of the transition velocity from v s to v e with the driving limitations, which is equal to The duration threshold T 2 is the transition interval for the distance L.
The duration threshold T 3 is calculated according to Eq. (25), when only the uniform motion exists.

2) REVISING THE FEEDRATE PROFILE
When T e < T 1 , the motion trajectory is revised as type A 1 , which includes the accelerating and decelerating motions, and the motion with the uniform velocity. Durations of the acceleration motion t a and the deceleration motion t d remain unchanged, and then the duration of the motion with the uniform velocity extends to It should be noted that the kinematic constraints in the acceleration and deceleration stages are respected well due to the reduced v m and invariable t a and t d . The revised velocity function v(t) after synchronization is constructed according to Eq. (19). Without losing generality, the revised velocity function is established with the similar procedures in following synchronization types.
If T 1 < T e < T 2 , the motion feedrate profile is revised as type A 2 , and then the duration of the accelerating phase t a is extended as It is noted that t a should be rounded as ceil(t a /T s ) · T s . The duration of the uniform stage is T e − t a , then v e is modified as If T 2 < T e < T 3 , the feedrate profile is updated as type A 3 so that the duration of the acceleration stage extends to T e . Moreover, v e depends on T e , which is calculated as If T 3 < T e , the profile is revised as type A 4 , and the duration of the deceleration stage t d is obtained as After rounding t d , the duration of the uniform stage is designed as T e − t d . The end velocity v e is revised as negative calculated end velocities. Moreover, the end velocity v e must be zero in some cases, especially in the last segment of the toolpath. To obtain the time synchronization in these cases, type A 5 is considered for designing the motion profile. If the sum of t d1 and t d2 in Eq. (32) is not higher than T e , the duration of the stage with the constant velocity is When t d1 + t d2 > T e , the Newtown-Raphson method is applied to iteratively calculate the constant velocity v m , the termination condition of which is |T ev − T e | ≤ T s . It should be indicated that upper and lower boundaries are set to 0 and v s , respectively.
Hold t d1 and t d2 and set T e as the motion duration of type A 5 by modifying the constant velocity v m according to Eq. (34).
where t d1 and t d2 are obtained from the acceleration and jerk limitations and the additional interpolation cycle T s . Otherwise, the driving limitation in one of the deceleration stages would be violated, since the difference value of the velocity increases caused by the modification the velocity in Eq. (34).

V. FEEDRATE SCHEDULING FOR THE FIVE-AXIS TOOLPATH
Following the geometry smoothing, the velocity profile along the toolpath is scheduled to generate the smooth motion. The toolpath between adjoining critical points and the midpoints of the smoothing curves is regarded as a block unit, which composes two adjacent smoothing curves and the remaining linear segment. The velocity profile of the transitional toolpath is planned in the WCS, which is determined by the normal and tangent constraints. Furthermore, the trajectory of the rotational toolpath is processed in the MCS, which involves the kinematic constraints based on the driving ability of the rotation motion axis to avoid the motor torque saturation. After planning of two trajectories separately, different durations are synchronized by extending the short duration to the longer time.  the initial velocity. It is noted that the initial velocity v s at the point S p is determined by a similar process in the last block. When S p is the start point of the whole path, v s is set to zero. Considering constraints of the predefined tangent acceleration, the velocity v ea at the critical point E p satisfies the following inequality.
where A t p and L p are the maximum acceleration limitation of the tool position and the length between the adjacent critical points, respectively.
Considering the constraint of the tangent jerk, the velocity at the critical point v ej satisfies the following inequality: where J t p is the tangent jerk limitation. The derivations of Eqs. (35) and (36) are provided in Appendix D.
The velocity at the critical point v e is defined as where A n p and J n p are the normal acceleration and jerk constraints of the tool position, respectively. Moreover, v m and κ p denote the maximum feedrate and the curvature at the critical point, respectively. a p and α p are the scaling parameter and the tangent angle at the corner, respectively. It is noted that the last iterm is used to constrain the value of the jerk jump at the junction of the biairthoid. Fig. 7 shows that the orientation toolpath smoothing in the MCS is constrained by the driving ability to avoid exceeding the motor saturation. Such constraints include the tangent/ normal acceleration/jerk limitations and the maximum

2) CRITICAL VELOCITY IN THE ORIENTATION TOOLPATH
where A t o , ω s and L o are angular acceleration limitation determined by the driving ability of the rotation motion axis, the initial velocity and the length of the smoothing segment between two adjoining critical points of the orientation toolpath, respectively.
Considering the predefined jerk limitation, the angular velocity ω ej at the critical point E o should satisfy the following inequality: where J t o is the jerk limitation of the rotational toolpath in the MCS, which is determined by the driving ability of the motion axis.
The angular velocity ω e is determined by where ω m is the maximum rotational velocity of the orientation toolpath. Moreover, A n o , J n o and κ o denote the normal acceleration and jerk constraints of the tool orientation, and the curvature at the critical point, respectively. a o and α o are the scaling parameter and the tangent angle at the corner, respectively.

B. KINEMATIC CONSTRAINTS OF LINEAR FEED DRIVES
Previous section discussed how to respect the kinematic constraints of the tool position in WCS and the tool orientation in MCS. By doing so, the physical constraints of the rotary axis are well bounded. However, due to the nonlinear transformation between the WCS and the MCS, the fair feedrate scheduling of the tool position in the WCS might lead to violate the linear feed drive constraints. It results in the oversize position error or instability of the motion control. To overcome the problem, a real-time optimizing strategy is developed for the feedrate profiles in this sub-section.
Firstly, designate some observation points of the five-axis toolpath, which includes the two critical points and some equally spaced points of the current block. Second, the positions of those points in the MCS are obtained through the inverse kinematics transformation (IKT) of the five-axis machine tool [39], [40]. Then, by differentiating the forward kinematics transformation (FKT) with respect to time, the relationship between the velocity of the tool tip in the WCS and axial velocities is conducted as where J T and J R are the first order differential of the tool position relative to the linear drive displacements and the rotational drive displacements, respectively. Denote (Ẋ ,Ẏ ,Ż ) T and (Ȧ,Ċ) T as v and w, respectively. Eq. (41) is rewritten as Take the i-th observation point hereby to illustrate the proposed optimizing strategy, and similar procedures hold for other points. To respect the velocity constraints of the linear and rotational axes, the path velocity f i r at the i-th observation point in the WCS is determined by Repeating the process, the allowable path velocity f m at the block is determined by considering the axis kinematic Differentiating Eq.(42) with respect to time, the acceleration of the tool tip in the WCS can be determined as follows.
where J TTk and J RRk (k = 1, 2, 3) are the differential of J Tk and J Rk to the drive displacements, respectively.v andẇ are (Ẍ ,Ÿ ,Z ) T and (Ä,C) T , respectively. Different from the path velocity, the path acceleration of the toolpath is the comprehensive expression of the tangent and normal accelerations. Still, take the i-th observation point as the example. The path acceleraion a i r at the i-th observation point in the WCS should satisfy the following inequality. The allowable tangent acceleration a m of the tool position in the WCS is the intersection of the axis kinematic accelerations and the scheduled tangent acceleration a lm . a m = min a 1 rt , a 2 rt , · · · , a N rt , a lm (47) Thus, the prefedined normal acceleration a ln is corrected as the allowable normal acceleration a n , which is equal to a m a lm ·a ln . Following similar procedures to determine the path acceleration, the path tangent and normal jerk constraints can be obtained analogously.

C. REAL-TIME LOOKAHEAD SCHEDULING BASED ON THE TIME SYNCHRONIZATION 1) BACKWARD SCANNING
The bidirectional scanning, consisted of the backward and forward scanning, is combination of the traditional scanning method and the time synchronization strategy. The backward scanning stage executes the feedrate planning from the final block to the start block based on the length of the smoothing unit and the curvature extremum at the critical point. First, velocities at the critical points, motion durations of the tool position and the tool orientation, are separately obtained. Second, referring the different durations, the time synchronization strategy is involved to eliminate the unmatched motion and then updates the backward critical velocities. Finally, the geometry massages and the obtained backward critical velocities are stored into the memory buffer. This process is repeated until the start point. Details of the backward scanning are demonstrated as the following:

2) FORWARD SCANNING
Forward scanning is performed from the start point to the end. In this scheme, the forward critical velocity is obtained based on the toolpath length and the tangent kinematic constraints. The critical velocity is the minimum of the backward and forward velocities, which is repeatedly implemented for the tool position and tool orientation sub-paths. Then, the final critical velocities of two paths are determined by the time synchronization strategy. Finally, the maximum velocity and acceleration of linear feed drives are evaluated and optimized to conform the physical limits. The critical velocities, the maximum velocities and the durations of the acceleration, uniform and deceleration stages, in the two paths are stored. The aforementioned procedures are repeated for all blocks, from the start block to the end block. The details are provided as following.

VI. SIMULATION
The proposed smoothing algorithm is employed to round the corners of the five-axis toolpath. To evaluate the advantages of the proposed method over the conventional schemes, it is intended to perform appropriate simulations in this section. The approximation errors for the linear and angular paths in the WCS are set to 80 µm and 60 µrad, respectively. The linear constraints are specified as the following: the desired feedrate, normal/tangent accelerations, normal/ tangent jerks are F = 12 mm/s, A n p /A t p = 100 mm/s 2 , and J n p /J t p = 1000 mm/s 3 , respectively. The maximum velocity, acceleration and jerk of linear feed drives are V l = 24 mm/s, A l = 150 mm/s 2 , and J l = 2000 mm/s 3 , respectively. Angular constraints are represented as the following: the desired feedrate, tangent/normal accelerations and tangent/normal jerks are ω = 0.1 rad/s, A n o /A t o = 1 rad/s 2 , and J n o /J t o = 10 rad/s 3 , respectively. The maximum velocity, acceleration and jerk of the rotational feed drives are V r = 0.1 rad/s, A r = 1.5 rad/s 2 , and J r = 20 rad/s 3 . The interpolation period is set as 1 ms. The smoothing strategy is performed in the MATLAB environment on a computer with 3.1GHz CPU.
In the five-axis spline interpolation, the forward and inverse kinematic transformations between the machine coordinate system and the workpiece coordinate system are dedicated to a specific machine tool model [39], [40]. The simulation in the study is performed on a table-tilting 5-axis machine tool model. The forward kinematics transformation can be expressed as where P x , P y ,  Solving Eq. (48), the inverse kinematic transformation from the WCS to the MCS is obtained as To verify the proposed algorithm for corner smoothing, a 5-axis toolpath composed of 4 linear segments and 3 corners is employed, as shown in Fig. 8. Fig. 9(a) shows the smoothing path of the tool position in the WCS. Fig. 9(b) shows the original and smoothing paths of the tool orientation in the MCS. The original linear segments of the rotary command are achieved from the cutter location data through the inverse kinematic transformation. It can be seen that, after performing the local smoothing, the corners are interpolated smoothly. The approximation errors of the tool position and tool orientation in the WCS are shown in Fig. 10. It is found that the corresponding errors are constrained well within the predefined tolerances.
To demonstrate the time synchronization handling ability, Fig. 11 shows the linear velocity of the tool tip position in the WCS and the angular velocity of the tool orientation in the MCS. When planed separately, the velocity profile of the tool tip position considers only the normal and tangent kinematic constraints of the tool position and ignores the angular kinematic constraints. The velocity profile of the tool orientation performs the same action. As a result, the durations of the linear and rotational trajectories in each block are different, as shown in Table 1. It is noted that the durations of the linear and rotational motions without the time synchronization are 7.886 s and 9.289 s. To address the problem, the proposed synchronization strategy shares the motion times through stretching the shorter duration to the longer one. Five synchronization types, as previously mentioned in section IV-B, are employed to accomplish the goal. It is found that the durations of the linear motion for the first and third blocks are longer than those of the angular motion, as presented in Table 1. Fig. 11(b) shows that the   angular trajectory is extended to share the same duration by using the synchronization type A1. For the second and fourth blocks, the durations of the angular motion are longer. The synchronization type A2 and type A5 are involved to stretch the linear trajectory to synchronize with the angular motion, as can be observed from Fig. 11(a). After synchronization, the durations of the linear and angular trajectories in each block are the same, as illustrated in Table 1. The total times are equal to 9.518 s, which are longer than any trajectory planed separately. What is more, Fig. 12(a) shows the acceleration and jerk profiles of the tool position in the WCS. Fig. 12(b) illustrates the acceleration and jerk profiles of the tool orientation in the MCS. It should be indicated that although the motion profiles of the subpaths are revised, the accelerations and jerks are constrained well under the predefined limitations.
To evaluate the advantages of the proposed method to conform the linear and rotational kinematic constraints, the traditional smoothing method with the quintic B-spline is employed. The feedrate profile of the tool position in the comparison method involves only the tool position kinematic constraints, including the normal/tangent accelerations and    jerks, desired velocity and chord error. Moreover, the tool orientation interpolation is considered as the slave motion, which is synchronized to the tool position motion by sharing the curve parameter. Fig. 13(a) shows the velocity and acceleration profiles of the tool position in the WCS, and Fig. 13(b) presents the velocity and acceleration profiles of the tool orientation in the MCS. It is found that the tool position kinematic constraints are well bounded. However, the situation is quite different for the tool orientation. The tool orientation kinematic constraints are aggressively violated in some blocks. The reason is caused by the curve-parameter synchronization procedure and regardless the rotary constraints. Fig. 14(a) shows the axial velocities and accelerations with the proposed method. Due to the tool orientation kinematic constraints and the time synchronization, the velocity and acceleration profiles of the rotational axis follow the specified limitations of the rotational motor. With the proposed strategy to confine the linear axial kinematic constraints, the velocity and acceleration of the linear axis feed drives are bounded. The axial velocities and accelerations with the traditional feedrate scheduling method are shown in Fig. 14(b). It is observed that the axis velocities and accelerations of A and C axes violate the limitations, which may result in exceeding the specified limitations of the rotational motor. Thus, the control stability of the feed drive will be destroyed and the visible mark on the workpiece may be confronted. Obviously, it is adverse to the machining process.

VII. EXPERIMENT
To verify the geometric smoothing method and the velocity planning strategy further, a 5-axis motion platform is designed, as shown in Fig. 15. The system is equipped with a computer, a motion controller dSPACE 1005 and an in-house five-axis motion table. The computer interprets NC codes and executes the interactive operation. The dSPACE hardware system is used to smooth the toolpath, generate axial commands and position control. For the in-house-developed five-axis machine tool, the X-axis feed system carries the Y-axis and a rotary table, and the C-axis is mounted on the A-axis. The general schematic diagram of the linear or rotational feed system in the five-axis machine tool is shown in Fig. 16. The feed drive is controlled by the traditional proportional, integral, and derivative (PID) controller,    which is composed of three independent feedback loops, namely position, velocity and current loops. Besides, an additional velocity feedforward controller is involved to improve the response performance. K p , K vp and K vi are the gain of position controller, the proportional and integral gains of velocity controllers, respectively. K fv is the gain of velocity forward controller. K t is the torque constant of the motor. J and B are the equivalent inertia and viscous damping of the feed drive, respectively. R g is the transmission ratio from the motor to the worktable. The tracking error caused by the friction is eliminated with the friction forward compensation strategy [41], [42]. Unless special stated, the main parameters are presented in Table 2.
The C 3 continuous arclength-parameterized curve is involved to smooth the tool position and tool orientation paths in this study. The arclength-parameterized curve can analytically calculate the toolpath length and effectively eliminate the feedrate fluctuation. Moreover, the smoothing path with airthoid curve generates smoother jerk profiles than the traditional C 2 continuous trajectory. To verify these, the cutter location data generated by CAM is conducted as the testing toolpath, as shown in Fig. 17. The approximation errors for the linear and angular paths in the WCS are set to 80 µm and 60 µrad, respectively. The adopted kinematical constraints are listed as follows. The linear constraints are designed as the following: F = 20 mm/s, A n p /A t p = 100 mm/s 2 , and J n p /J t p = 5000 mm/s 3 , respectively. Angular constraints are set as the  With the proposed smoothing method, the actual approximation errors of the tool position and tool orientation smoothing toolpath with airthoids in the WCS are shown in Fig. 18. It is found that the errors of the tool position and tool orientation are constrained within the responding tolerances. Table 3 presents the calculation time for the smoothing arclength in the first three corners of the tool position and tool orientation paths. It is found that the computation time of the airthoid is much shorter than that of the B-spline's. The reason is that the arclength of the airthoid is analytically obtained, while the value of B-splines should be iteratively calculated with the numerical method until the predefined approximation accuracy is satisfied. The specified tolerance of the B-spline in the case is set as 10 −14 mm. It illustrates that the proposed airthoid spline is more efficient to on-line smooth the corners, when the original toolpath consists of a large number of short linear segments or the length of the inserted spline occupies the vast majority of the smoothing toolpath. Fig. 19 shows the path linear velocity of the tool position in the WCS and the path angular velocity of the tool orientation in the MCS with three smoothing methods. As can be observed from Fig. 19 (b), the rate of angular velocity change with the M1 method is higher than the other methods. This phenomenon can be explained as follows. For the parameter synchronization method, the movement of the tool orientation is designed as the slave motion of the tool position. The abrupt motion change of the tool orientation leads to the large acceleration and jerk of the rotational axis. As a result, the kinematics of the tool orientation may be violated. Fig. 20 presents the feedrate fluctuation profiles based on the non-arclength parameterized splines with two traditional interpolation methods. One interpolation method is the second-order Taylor's expansion to approximately calculate the interpolation points for the given arclength [25]. The other one is a feedback interpolation method accurately mapping the interpolation points with the object arclength [27]. The feedrate fluctuation is the ratio of the difference between the desired feedrate and the actual feedrate to the desired feedrate. Fig. 20 (a) shows the feedrate fluctuation profiles of the tool position and tool orientation with the M1 method, and the updated parameter is obtained based on the second-order Taylor's expansion and the arclength of the tool position. Fig. 20 (b) shows the profiles with the M1 method, and the updated spline parameter is determined by a feedback correction scheme to map the arclength of the tool position. As can be seen, the feedrate fluctuation of the tool position is    evidently eliminated, whereas the fluctuation of the tool orientation has no significant change. It can be understood by the following reason. Though the updated parameter is revised on the basis of the arclength of the tool position, the incremental arclength of the tool orientaion is still unsmooth, which leads to the feedrate fluctuation of the tool orientation. In addition, Fig. 20 (c) presents the feedrate fluctuation profiles of the tool position and tool orientation with the M2 method, and the updated parameters of two sub-paths are severally obtained based on the individual arclength owing to the second-order Taylor's expansion. The feedrate fluctuation profiles with the feedback correction scheme is shown in Fig. 20 (d). It is noticed that the M2 method with the feedback correction scheme exhibits a better feedrate accuracy. However, the price to be paid for the achievement is computationally stringent, as shown in Table 4. It is noted that Table 4 illustrates the mean calculation time for the non-arclength parameterized splines with the interpolation methods, and the proposed smoothing strategy with the airthoid. The proposed method achieves shorter calculation time since it directly conducts the commands without the parameter correction module, as mentioned in Table 4. In addition, the feedrate fluctuation of the proposed method is zero, since the spline parameter is based on the arclength. It indicates that the proposed method with the arclength-parameterized airthoid can achieves smoother motion and release more computing resources to handle other tasks in the fine interpolation stage, which is friendlier to the CNC for the on-line executing smooth motion.    21 shows the axis acceleration and jerk profiles with three smoothing methods. As can be seen from Fig. 21(a), the axis accelerations with the M1 method in some blocks are much larger than the values generated with the time synchronization, which is more severe in X-, Y-and C-axis profiles. As a result, the axis jerk profiles of the M1 method exceed the scale range. This phenomenon can be explained as follows. For the parameter synchronization method, the tangent angular acceleration and jerk of the M1 method is extremely high, which leads to abrupt change of the rotational drives. Due to the kinematic transformation between the MCS and the WCS, the linear axis of the machine tool also be influneced. In contrast, M2 method generates lower amplitude accelerations and jerks in all the axes, as shown in Fig. 21(a) and (b). This result conforms to the fact that the kinematic constraints of the tool orientation are involved.
In addition, the proposed method generates the similar acceleration and jerk profiles with the M2 method due to the same consideration of the constraints of the tool orientation, as observed from Fig. 21(a) and (b). Moreover, the proposed airthoid achieves the C 3 continuity at the junction of the linear and curvilinear path and the jump-controlled jerk at the junction of the biairthoids than the M2 method. What is more, the airthoid spline can avoid the feedrate fluctuation, since it is an arclength-parameterized curve. As it can be seen in the Fig. 21(c), it is found that the proposed method with airthoid generates lower jerk profiles. Verification test illustrates that the proposed method with the airthoid spline and the jerk-continuous feedrate scheduling conducts smoother and bounded axial kinematics. Fig. 22 shows the mean values of the absolute tracking errors for the all motion axes. It is observed that the proposed  method with the airthoid curve performs better in reducing the position errors. This confirms that the proposed method provides smoother motion and higher tracking accuracy.

VIII. CONCLUSION
This paper develops a high-order continuous smoothing method based on the airthoid for the five-axis toolpaths. The high-order continuity, the peak-controlled jerk, the constraints of the geometry deviation, the motion synchronization of the tool position and the tool orientation, the kinematic constraints of the feed drives are guaranteed. On the one hand, a high-order continuous spline, airthoid, is proposed, the parameter of which is based on the arclength. The derivation of curvature increases linearly with the arclength. When inserted as the smoothing curve, it is analytically determined respecting the user-defined deviation errors of the tool position and the tool orientation. Compared with the numerical calculation in previous works, the arclength of the toolpath is analytically obtained. Moreover, without any iterative correction of the updated spline parameter at each interpolation timestamp, the feedrate fluctuation is eliminated for the arclength-parameterized expression of the airthoid. The high calculation efficiency and smooth motion make it more suitful to the real-time interpolation. On the other hand, instead of regarding as the slave motion, the rotational motion is synchronously scheduled with the transitional motion by sharing the motion time. Thus, the angular kinematic constraints are strictly respected, which avoids exceeding the physical limits of the rotary feed drives. To comply with the constraints of the linear feed drive, an on-line optimization strategy is involved. As a result, the proposed method successfully plans the velocity profiles under the physical limits of the machine tool. Furthermore, the high-order continuous airthoid with the peak-constrained jerk improves the smoothness of the toolpath. It generates lower and smoother jerk profiles of the drives than the traditional C 2 interpolation. The experiments verify that the proposed smoothing algorithm achieves higher tracking accuracy.

APPENDIX A G 3 CONTINUITY AT THE CONJUNCTION OF THE LINEAR AND CURVILINEAR SPLINES
This is an analytical discussion on the jerk continuity at the conjunction of the original linear segment and the airthoid spline. As shown in Figure. 23, the corner, constructed by lines P 1 P 2 and P 2 P 3 , is smoothed by a biairthoid spline B 1 B 2 , which comprises two symmetric airthoid splines. A local coordinate system {t − n} is built at the conjuction B 1 Substituting Eq.(51) into Eq.(50), the differentials of the airthoid with respect to the time at the conjunction arė P s=0 + =ṡt;P s=0 + =st; ...
The differentials of the linear segment with respect to the time at the conjunction arė P s=0 − =ṡt;P s=0 − =st; ...
The result indicates that the G 3 continuity is proved with the continues velocity, acceleration and jerk at the joint point.

APPENDIX B G 2 CONTINUITY AT THE CONJUNCTION OF THE BIAIRTHOID
Based on Figure. 23, an another local coordinate system t − n is built at the conjuction B 2 , as shown in Figure. 24.
The arclengths between the conjuction C and B 0 and B 1 are s + j and s − j , respectively. The length of B 0 B 1 is equal to s total . The relationship between them is The first, second and third order time derivatives of the equation are expressed aṡ s − j = −ṡ + j ;s − j = −s + j ; ...
The first, second and third order differentials of the airthoid with respect to the arclength at the conjunction in the coordinate frame {t − n} are Similarly, combined with Eq. (54), the first, second and third order differentials of the airthoid with respect to the arclength at the conjunction in the coordinate frame t − n Combined with Eq. (50), (58) and (59), the difference between the first-order time derivative of the curve at the conjunction C isṖ Similarly, the difference between the second-order time derivative of the curve at the conjunction C is However, the difference between the third-order time derivative of the curve at the conjunction C is ...
To constrain the peak jerk, the jerk jump is no more than two times the normal jerk J n . Therefore, the velocity at the conjunction should be restricted as v ≤ 9 J 3 n a 6 24π 2 α (65)

APPENDIX C THE DETERMINATION OF THE SMOOTHING SPLINE
Set P i is the original point of the planar coordinate system. As mention before, the curvature maximum generates at the critical point Q. The radius of curvature is defined as r α . The coordinates of the center point are x c = x(α) − r α sin α = aC(α) − r α sin α y c = y(α) + r α cos α = aS(α) + r α cos α