Abstract
In micro-line segments machining, transition curves with high harmonic components are more prone to causing vibration issues in the feed drive system, which affects machining efficiency and quality severely. To construct low harmonic trajectories, this paper proposes a corner smoothing algorithm that uses the Trajectory Pattern Method (TPM). The transition curve construction and axial motion scheduling are performed with a specified fundamental frequency in one step, which reduces the smoothing process time and avoids excitation of natural modes of vibration of the system. The synthesized trajectories and axial kinematic profiles are all smooth and only contain the selected fundamental frequency and its first two odd harmonics, which minimizes the number of high harmonic components in the required actuation forces/torques and avoids excitation of the system modes of vibration. Linear programming is used to synthesize the trajectories. The proposed algorithm is shown to achieve near time-optimal trajectories. The provided experimental analysis and comparisons demonstrate that the proposed algorithm achieves smooth axial kinematic profiles with low harmonic contents, which would improve machining efficiency and quality.
Introduction
The most common practice in CNC machining is to express tool paths by small line-segments. Research work into high-speed and high-quality machining has concentrated on methods to improve smoothness of trajectories. To this end, several researchers have investigated methods for smooth tool path generation.1,2 Some researchers have focused on the Local Smooth Method (LSM), in which transition curves are constructed at the junction of small line-segments to improve corner smoothness and velocity since it is simple to implement and is easily controlled. For example, Lv et al. 3 used arcs to blend small line-segments. Zhang et al. 4 constructed transition curves with quintic polynomials and Ernesto and Farouki 5 used rational quadratic Bézier curves. These three algorithms, however, only provide G1 continuity, that is, continuous velocity, at the corner, which would lead to vibration caused due to the present normal acceleration. To obtain G2 continuity, that is, continuous curvature, Du et al. 6 smoothed small line-segments with cubic Bézier curves. The calculation of a proper transition length, however, requires a significant amount of time. To reduce vibration during the machining process, Pythagorean-Hodograph (PH) spline has drawn a lot of attention for corner smoothing studies because it provides analytical solution of arc-length related to spline parameters. Shi et al.7,8 developed corner smoothing algorithms for multi-axis tool paths based on quintic PH splines, which avoided the approximation errors of arc-length calculation. In order to simultaneously achieve real-time interpolation and higher order continuous motion, Hu et al. 9 proposed a C3 continuous analytical corner smoothing algorithm by inserting the specially designed micro PH splines among linear tool path segments. Then, Hu et al. 10 extended the algorithm to construct C3 continuous trajectories for five-axis machine tools. Fan et al. 11 enhanced the continuity further by utilizing two quartic Bézier curves to achieve G3 continuity, that is, smooth curvature. Smooth transition curves and tangential motion profiles can be achieved by above algorithms. However, the methods require two computational steps. In the first step, parametric curves are used at the corners to construct the transition schemes. In the second step, tangent motion profiles are scheduled based on the synthesized trajectories. These algorithms are therefore inefficient. In addition, sudden changes in the axis velocity, acceleration, and/or jerk profiles are encountered since only tangential kinematic limits are considered.
In recent years, one LSM has drawn researchers’ attention, in which in one step, the transition schema construction and motion scheduling are performed. Biagiotti and Melchiorri 12 smoothed tool paths using finite impulse filters. The method, however, leads to large contour errors and unavoidable time delay. Tajima and Sencer 13 developed a direct axis motion profile blending approach to obtain transition curves. Sudden changes in axis jerk may, however, occur at corners due to the use of a method based on jerk-limited acceleration and deceleration process. Zhang and Du 14 used a jounce-limited velocity planning method to achieve transition curves and smooth axis acceleration profiles, but axis jerk profiles are still not smooth.
It is important to note that when expressed in Fourier Series, all above synthesized trajectories, such as those based on polynomial, Akima, or B spline curves contain significant high harmonic components. In addition, the nonlinearity of the machine system dynamics would require actuating forces/torques with even higher harmonic content to follow the prescribed trajectories. The lack of control over the harmonic content of the required actuating forces/torques and their frequencies in currently available and published methods for trajectory synthesis leads to excitation of natural modes of vibration of machining systems. The commonly used actuator input filtering does not alleviate the problem since the filtered components are needed to achieve the prescribed motion and results in path and trajectory errors.
Rastegar 15 introduced the concept of trajectory patterns and their application to synthesis of trajectories with low harmonic content such that the harmonic content of the required actuation forces/torques is controlled and does not contain significant high harmonic content.15–18 The Trajectory Pattern Method (TPM) was then introduced for trajectory synthesis and model-based control of high-speed and precision machinery.15,16 Using TPM, the fundamental frequency of the synthesized trajectories can be selected such that together with all its harmonics that appear in the trajectory and the required actuating forces/torques would not excite the natural modes of vibration of the system. In TPM, trajectory patterns are classes of parametric trajectories that are used to synthesize trajectories with the desired characteristics and derive the system inverse dynamics in an algebraic form.15,16 The actuating forces/torques can then be readily optimized and used in a model-based control schema to achieve high-speed and high-precision operation of the system.
It has been shown that the trajectory patterns that are constructed with sinusoidal functions with an appropriate fundamental frequency can be used to synthesize trajectories with low harmonic content that would satisfy the dynamic response limitations of the system actuators while avoiding excitation of the natural modes of vibration of the system.16,18
In this paper, a corner smoothing algorithm is presented that uses harmonic based trajectory patterns for high-speed and high-quality machining. In Section 2, it is shown that transition curve construction and axis velocity planning with a specified fundamental frequency can be performed in one step using the TPM. The synthesized trajectories and axis kinematic profiles are smooth and only contain the fundamental frequency and its first two harmonics. Linear programming is used to synthesize the trajectories and the proposed algorithm is shown to achieve near time-optimal trajectories. In Section 3, the effectiveness of the proposed algorithm is illustrated through the results of computer simulations studies. Conclusions of the present study are presented in Section 4.
Low harmonic transition schema
In this section, a Trajectory Pattern Method (TPM) based transition schema, which achieves smooth axis motion profiles with chord error and axis kinematic constraints is presented. In this scheme, trajectory patterns that are synthesized with a fundamental frequency and its first two harmonics are used to achieve minimum trajectory harmonic content. 18
Trajectory pattern
The class of trajectory patterns used in the present study is harmonic based and is formed by a fundamental sinusoidal function and (n − 1) number of its harmonics:
where ai and bi are trajectory parameters; t is the time,
In this paper, the trajectory pattern (1) that consists of a fundamental frequency sinusoidal function and its first two odd harmonics is used to construct the kinematic functions of the transition schema. The velocity function is then expressed as,
From equation (2), the acceleration and jerk functions are then derived as,
By integrating equation (2) with respect to the time t, the displacement function is derived as,
Transition schema construction
Figure 1 shows the flow diagram of the developed transition schema. The transition schema construction includes five modules: an initialization module; a boundary conditions definition module; a motion profile constraint definition module; the contour error constraint definition module; and a trajectory parameters calculation module. The above modules are described in more detail in the following sections.

Flow diagram of the developed transition schema.
Initialization module
Figure 2 presents the proposed transition schema. In Figure 2, let the point Pi,s be the origin to establish the x–y coordinate system and lengths of line-segments Pi−1Pi and PiPi+1 be li and li+1, respectively. The two line-segments Pi−1Pi and PiPi+1 form a corner, with θ1 indicating the angle between the vector

Transition schema.

Smooth axial motion profiles: (a) displacement, (b) velocity, (c) acceleration, and (d) jerk.
The user is to provide the required chord error ε and the fundamental frequency f of the trajectory pattern, which must be selected such that together with its first two odd harmonics, the natural modes of vibration of the system are not excited.
Boundary conditions
From the indicated end conditions and equations (2)–(5), the following boundary conditions are obtained,
From equations (6)–(12), it is clear that if the cycle time T is set equal to half the fundamental frequency period, that is, if
Motion profile constrains
To achieve the desired motion profiles shown in Figure 3, the condition a(t) ≥ 0,
where
The functions g1(t), g2(t), and g3(t) with arbitrary amplitudes and the resulting acceleration a(t) are plotted as shown in Figure 4. As expected, and can be seen in Figure 4, the extreme values of a(t) can only occur at one of the extreme points of the fundamental frequency function g1(t), or the functions g2(t) and g3(t). This is easily seen to be the case since the functions g2(t) and g3(t) are the first two odd harmonics of the fundamental sinusoidal function g1(t). Therefore, the condition, a(t) ≥ 0,

Acceleration component profiles.
The extreme points of the functions g1(t), g2(t), and g3(t) may be calculated as follows:
As an example, for the function g1(t), the extreme of this sin(j) function occurs at j = (2k + 1)π/2,
Therefore, from equation (18) for the time
Then in general, for all integers k satisfying equation (22), the extreme point times are determined for the function from equation (21). For the case of the g1(t), the extreme point is readily seen to be at,
For the functions g2(t) and g3(t), the extreme point times are similarly determined to be at (see Figure 4),
Now due to symmetry of the functions about t = T/2, to achieve the desired condition of a(t) ≥ 0,
Contour error constraints
To control machining accuracy, the chord error should be bounded. In Figure 2, the transition schema is seen to be symmetric about the bisector of the angle
For the corner motion displacement during the period T, from Figure 2, the corner length Li,c and the displacement sx(T), equation (5), are expressed as,
The corner length Li,c is then determined from equation (26) as,
The projection of corner error ε i on the x-axis is seen in Figure 2 to be,
The x coordinate of the point Pi is also seen to be,
From equation (5) at the time t = T/2, the x coordinate of the point Pi,m is,
As can be seen in Figure 2, the x-axis contour error ε i,x is,
Substitute equations (27)–(30) into equation (31), the maximum chord error ε i is obtained as,
where
To control the machining accuracy, we should guarantee 0 ≤ ε i ≤ ε.
Trajectory parameters calculation
To calculate the trajectory parameters for the x-axis motion profiles, the initial corner velocity vi,c is set at its desired value. From equations (14) and (15), the trajectory parameter a0 is then derived as,
Substitute equations (13) and (36) into equation (27), the corner velocity vi,c in terms of the corner length Li,c is given as,
Now let Li,c be defined as the shorter of half the lengths of the line-segments Pi−1Pi and PiPi+1, that is, li and li+1, respectively, as,
Substituting equation (38) into equation (37), the initial velocity vi,c is then calculated.
It is noted that the boundary conditions, equations (13)–(16); motion profile constraints, equation (25); and the contour error constrain, equation (32), are all linear functions of the trajectory parameters ac, a0, a1, a3, and a5. Thus, their optimal values can be calculated using linear programming. The linear programming problem for determining optimal trajectory parameters to achieve minimum chord error ε i becomes,
If the motion is in a deceleration phase, then the acceleration constraint becomes a(t) ≤ 0,
To synchronize the corner motion, the x-axis and y-axis motion trajectories should be constructed with the same corner velocity vi,c. The process can be as follows:
Simulation validations
In this section, to demonstrate the effectiveness and performance of the proposed corner smoothing algorithm, the results of performed simulations on a typical selected system is presented and compared with results that would be obtained with one of the commonly used trajectory synthesis algorithms. Both algorithms were implemented in Visual Studio 2013 and MATLAB 2015 was used for plotting the results. The following typical settings were used: a fundamental frequency of f = 16 Hz, a sampling period of T = 2 ms, a chord error limit of ε = 0.1 mm, a programmed feedrate of F = 100 mm/s, acceleration limits of the drive axes Ax_max = Ay_max = 2000 mm/s2, and jerk limits of the drive axes Jx_max = Jy_max = 100,000 mm/s3.
Case one
To evaluate the proposed algorithm, a jerk-limited algorithm, 13 which constructs the transition schema and jerk-limited axial kinematic profiles in one step, is used as the comparison algorithm. As shown in Figure 5, a simple pattern with one acute angle is indicated as the desired tool path. The synthesized transition curves by applying the two algorithms are plotted. The cycle time with the proposed algorithm was calculated to be 0.0313 s, which is 9.5% faster than the jerk-limited algorithm of 0.0346 s. Figure 6(a) shows that neither algorithm exceeds the specified chord error limit ε. From Figures 5 and 6(b) to (f), we can see that the tangential and axial displacement and velocity profiles are all smooth. From Figure 6(g) to (i), it is observed that the tangential acceleration profiles of both algorithms are smooth, while the axial acceleration profiles of the comparison algorithm only have G0 continuity. In addition, since the comparison algorithm constructs axial kinematic profiles under axial acceleration and jerk limits, its jerk profiles exhibit sudden changes as shown in Figure 6(j) to (l). Noting that since such unsmooth acceleration and jerk profiles cannot be realized in the feed drive systems, the synthesized trajectories would most likely excite natural modes of vibration of the system. The proposed algorithm, however, has smooth tangential and axial velocity, acceleration, and jerk profiles, as shown in Figure 6, which the feed drive system should be capable of realizing without exciting natural modes of vibration of the system. The proposed algorithm should, therefore, achieve significantly better machining quality.

Tool path and synthesized trajectories.

Comparison of the kinematic profiles for the Case one: (a) contour error, (b) x-axis displacement, (c) y-axis displacement, (d) tangential velocity, (e) x-axis velocity, (f) y-axis velocity, (g) tangential acceleration, (h) x-axis acceleration, (i) y-axis acceleration, (j) tangential jerk, (k) x-axis jerk, and (l) y-axis jerk.
To evaluate the proposed algorithm further, frequency spectrum distribution of the synthesized axial kinematic profiles is studied using Fourier integration. In the plots of Figure 7, it can be seen that the x-axis and y-axis displacements, velocities, accelerations, and jerks of the trajectories synthesized using the aforementioned jerk-limited algorithm contain a considerable number of high harmonic components, that is, harmonics with frequencies above the second odd harmonic of the fundamental frequency of the trajectory pattern. As a result, the related actuating forces/torques require producing high frequency components, which for high speed motions would in general be beyond their dynamic response and could excite natural modes of vibration of the system. As expected, Figure 7, the velocity, acceleration, and jerk profiles of the synthesized trajectory with the proposed algorithm only include the fundamental frequency and its first two odd harmonics.

Comparison of the Fourier analysis for the Case one: (a) x-axis displacement of the jerk-limited algorithm, (b) x-axis displacement of the proposed algorithm, (c) y-axis displacement of the jerk-limited algorithm, (d) y-axis displacement of the proposed algorithm, (e) x-axis velocity of the jerk-limited algorithm, (f) x-axis velocity of the proposed algorithm, (g) y-axis velocity of the jerk-limited algorithm, (h) y-axis velocity of the proposed algorithm, (i) x-axis acceleration of the jerk-limited algorithm, (j) x-axis acceleration of the proposed algorithm, (k) y-axis acceleration of the jerk-limited algorithm, (l) y-axis acceleration of the proposed algorithm, (m) x-axis jerk of the jerk-limited algorithm, (n) x-axis jerk of the proposed algorithm, (o) y-axis jerk of the jerk-limited algorithm, and (p) y-axis jerk of the proposed algorithm.
To explain how trajectories with high harmonic components affect vibration of the machine tool, an undamped, 1-DOF, linear spring-mass system modeling the behavior of a vibration isolation system shown in Figure 8(a) is selected. Let the spring have a spring rate and free length of k and l, respectively. The mass of the block is m. The displacements of the mass and the system from their initial position are u and d, respectively. From the Newton’s second law, the system equation of motion is seen to be,

(a) Mass-spring system and (b) its vibrations.
Let
Then, equation (43) can be rewritten as,
Now as an example, take the y-axis motion profiles calculated using the aforementioned two algorithms as input to the dynamic system described by equation (45). The input motion d(t) is then y-axis acceleration profile, that is,
Substituting
As can be seen in Figure 8(b), the displacement z, which indicated the tool displacement, of the proposed algorithm is smoother than that of the jerk-limited algorithm. In addition, when the y-axis motion ends, the mass m comes to a stop with the proposed algorithm, while the mass of the jerk-limited algorithm continues to vibrate with its natural mode. It is noted that the trajectories synthesized using the jerk-limited algorithm contain a significant number of high harmonic content, which would in general cause excitation of the system modes of vibration irrespective of their frequencies. The synthesized trajectories of the proposed algorithm, however, only contain the selected fundamental frequency and two of its odd harmonics in its velocity, acceleration, and jerk terms (and even in their higher order derivatives), which would minimize the number of high harmonic components in the required actuation forces/torques. The fundamental frequency of the trajectory pattern can also be selected to avoid natural modes of vibration of the system. As a result, the system should be capable of operating with minimal vibration and achieve higher quality machining.
Case two
In this Case, a polynomial algorithm, 4 a jerk-limited algorithm, 13 and a jounce-limited algorithm 14 are implemented as the basis for comparison with the proposed Trajectory Pattern Method (TPM) based algorithm. The polynomial algorithm constructs the transition schema with the quintic polynomial in arc length domain, then achieves tangential kinematic profiles by using a traditional tangential trigonometric velocity scheduling method. The jerk-limited algorithm uses a jerk-limited velocity scheduling method to obtain the transition schema and jerk-limited axial kinematic profiles in one step in time domain. The jounce-limited algorithm utilized a jounce-limited velocity scheduling method to construct smoother acceleration transition profiles, which can reach G3 continuity. The proposed algorithm adopts the TPM to get the transition schema and jerk-smooth axial kinematic profiles in one step in time domain as described in Section 2.
The simulation conditions of this Case are the same as those of the Case one. Figure 9(a) shows the trident-shaped pattern that is selected as the tool path. The synthesized corner tool paths with the above three algorithms are shown in Figure 9(b).

(a) Trident-shaped tool path and (b) enlarged corners.
The comparison of the simulation results is demonstrated in the plots of Figure 10. Figure 10(a) shows that tangential velocities of the four algorithms, which are all smooth and do not exceed the programmed feedrate F. As shown in Figure 10(b) and (c), axial velocity profiles of the jerk-limited, jounce-limited, and proposed algorithms are smooth and similar, while those of the polynomial algorithm is not smooth at the corners. In addition, Figure 10(d) to (g) show that axial acceleration and jerk profiles of the polynomial algorithm exceed specified limitations and fluctuate frequently around corners since this algorithm does not use an axial velocity planning method, but a traditional tangential velocity scheduling method. In the case of jerk-limited algorithm, its axial acceleration profiles are also not smooth and its axial jerk profiles exhibits sudden changes at the beginning and end of acceleration or deceleration periods around the corners. These fluctuations and sudden changes cannot be realized in the feed drive system, which makes the machining process prone to vibration.

Comparison of kinematic profiles for Case two: (a) tangential velocity, (b) x-axis velocity, (c) y-axis velocity, (d) x-axis acceleration, (e) y-axis acceleration, (f) x-axis jerk, and (g) y-axis jerk.
As can be seen in the plots of Figure 10, in contrast to polynomial, jerk-limited, and jounce-limited algorithms, the proposed algorithm has smooth axial velocity, acceleration, and jerk profiles, which should make it achievable by the drive systems of the machinery system.
Fourier analysis results of the axial acceleration profiles of the four algorithms showing their frequency contents are shown in Figure 11. As can be seen in the plots of Figure 11(a) to (c) and (e) to (g) for the polynomial, jerk-limited and jounce-limited algorithms, respectively, the axial acceleration profiles of these three algorithms contain a large number of high frequency components since their trajectories are synthesized by polynomial basis curves. The trajectories of the proposed algorithm, however, are constructed with only three harmonic functions, a fundamental frequency sinusoidal function and its first two odd harmonics, in the present simulations with a fundamental frequency of 16 Hz and the first two odd harmonic frequencies of 48 and 80 Hz. The plots of (c) and (f) of Figure 11 also shows that the maximum axial acceleration frequency for synthesized trajectory with the proposed algorithm is 80 Hz, which should minimize the number of high harmonic components in the required actuation forces/torques and minimize the excitation of the system modes of vibrations. The proposed algorithm should, therefore, yield a significantly better machining quality and flexibility.

Comparison of frequency distribution for Case two: (a) x-axis acceleration profile of polynomial algorithm, (b) x-axis acceleration profile of jerk-limited algorithm, (c) x-axis acceleration profile of jounce-limited algorithm, (d) x-axis acceleration profile of proposed algorithm, (e) y-axis acceleration profile of polynomial algorithm, (f) y-axis acceleration profile of jerk-limited algorithm, (g) y-axis acceleration profile of jounce-limited algorithm, and (h) y-axis acceleration profile of proposed algorithm.
The calculated computing and machining times and contour errors for the polynomial, jerk-limited, jounce-limited, and proposed algorithms are presented in Table 1. The results show that the computing time cost for the polynomial algorithm is almost twice that of the jerk-limited, jounce-limited, and proposed algorithms. It can also be seen that the machining time spent by the proposed algorithm is the shortest (1.7203 s), which is 4.7%, 1.8%, and 5.3% shorter than the jounce-limited (1.8058 s), jerk-limited (1.7515 s), and polynomial (1.8165 s) algorithms, respectively. The reason for the proposed algorithm providing the shortest machining time is due to the use of an axial velocity planning method and since no limit on the axial jerk maximum. It can, therefore, be concluded that the proposed algorithm is near time-optimal and satisfy the requirement of real-time machining, while yielding a significantly better machining quality.
Comparison of computing and machining times and contour error.
Figure 12 shows that none of the four algorithms exceeds the prescribed chord error ε, and that the maximum always happens around the corners. As can be seen in Table 1, the contour error maximum for the proposed algorithm is 0.0740 mm, which is 0.9%, 1.2%, and 14.9% lower than those of the jounce-limited, jerk-limited and polynomial algorithms, respectively. In addition, compared with the other three algorithms, the average contour error for the proposed algorithm is also smaller, that is, it is 0.0017 mm as compared to 0.0023 mm for the polynomial, 0.0019 mm for the jerk-limited, and 0.0020 mm for the jounce-limited algorithms. The proposed algorithm can, therefore, better satisfy machining accuracy requirements, while yielding a significantly better machining quality.

Contour errors.
Conclusions
In this paper, a corner smoothing algorithm using the Trajectory Pattern Method (TPM) for high-speed and high-quality machining is proposed. In this algorithm, the transition schema and axial kinematic profiles are constructed simultaneously with a fundamental frequency sinusoidal function and its first two odd harmonics in time domain.
In comparison with previously published algorithms, the proposed algorithm has the following advantages: (I) the transition schema and axial kinematic profiles are constructed with a sinusoidal function with a selected fundamental frequency in one step, which reduces the smoothing process time, while avoiding excitation of the natural modes of vibration of the machinery system; (II) all velocity, acceleration, and jerk profiles are expressed by a sinusoidal function consisting of a selected frequency harmonic and its first two odd harmonics, which minimizes the number of high frequency harmonic components in the required actuation forces/torques, reducing their dynamic response requirements, and by proper selection of the trajectory fundamental frequency, the excitation of the system modes of vibration is minimized; (III) the proposed algorithm can realize smooth axial kinematic control action under chord error and dynamic response limitations of the axial drives, noting that for high-speed and high-quality machining, smooth axial kinematic profiles are more important than smooth tangential kinematic profiles; and (IV) the proposed algorithm is near time-optimal and can significantly improve machining efficiency and quality.
Since the proposed algorithm is a general corner smoothing method, it could be extended to five-axis machines and its machining efficiency and quality would be improved more significantly. However, compared to the three-axis machine, there are mainly two challenges for a five-axis corner smooth algorithm. First, the tool tip position and tool orientation need to be smoothed and synchronized. Second, due to the non-linear relationship between tool orientation vectors and rotary axis coordinates, the smooth motion of the rotary axis is not capable of being guaranteed by the smooth process of tool orientation vectors. Especially when tool orientation vectors approach singular points, fierce vibrations of the rotary axis would cause the servo alarm, even damage machine parts. Efforts should be, therefore, made to obtain smooth motion of rotary axes under the tool orientation deviation limitation in the five-axis tool center point interpolation processing.
Future work will focus on two aspects. First, we are in the process of implementing the proposed algorithm and more comparison algorithms on an open-source CNC system to analyze the performance further. Second, more limitation conditions will be considered in the proposed algorithm. We plan to publish the results shortly after.
Footnotes
Handling Editor: James Baldwin
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Scientific Research Startup Foundation for High-level Talents of Nanjing Institute of Technology (Grant No. YKJ2019112), the Natural Science Foundation of Jiangsu Provincial of China (Grant No. BK20181024), and the Shuangchuang Doctoral Project of Jiangsu Province.
