Abstract
In spline toolpath interpolation, a crucial point is solving the mapping between the spline parameter (u) and actual arc length (s) accurately, so that the toolpath is traveled without undesirable fluctuations or discontinuities in the feedrate profile. To achieve this, various techniques have been proposed in literature, including Taylor series interpolation, iterative numerical methods, and approximating the mapping between u and s with a feed correction polynomial. This paper presents a new way to parameterize the seventh order feed correction polynomial, which was introduced by Erkorkmaz and Altintas (2005, “Quintic Spline Interpolation With Minimal Feed Fluctuation,” ASME J. Manuf. Sci. Eng., 127(2), pp. 339–349). The proposed technique has a closed-form solution that can be efficiently implemented in real-time, rather than having to construct and solve a linear equation system with 14 unknowns for each spline segment. In this paper, the new solution is derived step by step, and simulation case studies are presented which demonstrate that the new method accurately parameterizes the feed correction polynomial in approximately 43% less computational time, compared to applying the former solution of Erkorkmaz and Altintas (2005, “Quintic Spline Interpolation With Minimal Feed Fluctuation,” ASME J. Manuf. Sci. Eng., 127(2), pp. 339–349). This is because matrix multiplication operations and a dedicated linear equation solver, which are cumbersome to implement inside a real-time computer numerical controller (CNC), are avoided in the new solution.
Introduction
Spline toolpaths offer many advantages in machining freeform shaped components with complex geometry, such as those produced in aerospace, die and mould, automotive, and biomedical applications. In addition to subtractive machining operations [1], they can also be applied to additive manufacturing [2]. Compared to using linear or circular segments, spline toolpaths enable smoother geometry, better surface finish, and also shorter cycle time due to allowing more aggressive feeds and accelerations to be specified. However, spline toolpaths present a challenge in interpolation, due to the inherent discrepancy between the variable (), which they are parameterized with respect to, and the actual arc length () that is traveled, as the spline parameter is incremented.
Figure 1 shows a quintic (fifth degree) spline toolpath, fit in x and y axes, formulated in the form
Earlier research in spline toolpath generation has shown that for polynomial type splines (e.g., cubic, quintic, etc.), there exists no known method that can parameterize to be equivalent to throughout the whole toolpath. However, polynomial splines are preferred in industrial CNCs due to their simplicity in formulation and implementation. Hence, when the spline parameter is interpolated in uniform increments, the discrepancy between and exhibits itself in the form of two detrimental effects [3]:
- (i)
feedrate fluctuations within each segment;
- (ii)
feedrate discontinuity at segment the connections, which leads to acceleration and jerk spikes.
To alleviate these problems, several approaches have been proposed in literature.
One approach has focused on improved parameterization methods for polynomial splines, which has led to the nearly arc length parameterized [4], approximately arc length parameterized C3 [5], and optimally arc length parameterized [3] quintic splines. These methods have shown significant reduction in the feed fluctuation and discontinuity, compared to traditional cubic splines that are currently used in industry. However, they require extra off-line computation.
Another approach has explored different spline structures, such as Pythagorean-hodograph curves [6,7], which allow the arc length to be analytically computed, thereby avoiding the feed inconsistency problem altogether. However, there is still a strong trend in industry to stick with cubic and polynomial splines, which are most easily understood, parameterized, and implemented by practicing engineers. Furthermore, cubic splines provide the necessary degree of freedom to represent toolpaths with C2 continuity, which is sufficient for many of the precision motion control applications in industry.
To improve the feedrate accuracy and continuity during the interpolation of preparameterized toolpaths (such as cubic B-splines), improved interpolation methods have also been proposed. Figure 2 illustrates three common techniques for spline interpolation. The first one is natural interpolation, which is the easiest approach, but it leads to the mentioned feed fluctuation and discontinuity problems. The second one is the approximation of the spline parameter (), by expanding it into a Taylor series as a function of time, around its value used in the previous interpolation cycle (). This approach, which has been proposed in first order [8–10] and second order [11] forms, unfortunately suffers from accumulating round-off errors. These can be problematic when interpolating long toolpaths, unless adequate numerical measures are incorporated into the practical implementation. The third approach, shown in Fig. 2(c), approximates the nonlinear relationship between and with the so-called “feed correction polynomial,” as proposed in Ref. [3]. Typically, seventh degree is used, which allows the six boundary conditions on , (), and () to be imposed at the segment connections; thus avoiding the velocity and acceleration discontinuity problem at these points. The additional two degrees of freedom help establish a reasonably close fit between the feed correction polynomial and the spline parameter () versus arc length () data, which is generated earlier on during the numerical integration of the total arc length () for each spline segment. Fitting of the feed correction polynomial is shown in Fig. 3.
In recent research, the feed correction polynomial has also been applied in an adaptive manner to nonuniform rational B-splines (NURBS) interpolation, in order to constrain the mean squared value of the approximation error for the chord parameter [12]. It has also been extended to the ninth degree, in order to interpolate quintic B-splines with jerk continuity [13]. Other polynomial approximations between the spline parameter and arc length have also been proposed in literature, such as the remapping scheme which employs cubic Hermite splines, as reported by Lei et al. [14].
Aside from the above mentioned solutions, iterative numerical techniques have also been explored in Refs. [3,15]. These techniques yield the most accurate results, but typically require at least 2–3 times additional computational time, compared to natural or feed correction type interpolation. Additional measures also need to be incorporated to ensure the stability of iterative algorithms, especially those relying on the Newton–Raphson method. For brevity, in this paper iterative techniques are not illustrated or benchmarked. It was also demonstrated in Ref. [3] that using the feed correction polynomial can provide a good initial guess for a Newton–Raphson based iterative solution.
Overall, feedrate fluctuation avoidance methods, whether they are based on improved toolpath parameterization, Taylor series approximation of and , polynomial function mapping, or real-time iterative solution, target accurate achievement of the desired arc displacement (), and therefore the correct feedrate at a given point along the fitted toolpath. These methods are meant to be used in conjunction with feedrate planning algorithms (or so-called look-ahead functions), like those reported in Refs. [12,14], which optimize the feed profile along the toolpath typically subject to axis velocity, acceleration, and jerk constraints. While feedrate planning is a much broader subject, outside the scope of this technical brief, without accurate interpolation methods, such as the feed correction polynomial approach adopted in this paper, the full advantages of careful jerk and acceleration limiting cannot be appreciated in practice.
In the remainder of this technical note, brief formulation of the feed correction polynomial is reviewed in Sec. 2. Only the essential information, required to follow the proposed new formulation is presented, and details of secondary importance already available in Ref. [3] are skipped. Section 3 presents the derivation of the new method for computing the feed correction polynomial. Section 4 presents implementation results, which confirm the correct solution and effectiveness of the feed correction polynomial, and provide a computational load and time benchmark with respect to the earlier method in Ref. [3].
Background on Fitting of the Feed Correction Polynomial
Above, ( = 1,…, ) is the numerically integrated arc displacement corresponding to the spline parameter value . is the number of divisions, which are distributed uniformly along [0,]. For each spline parameter value , the corresponding increments in axis positions (,) are calculated using the toolpath definition in Eq. (1) and the resulting arc displacement is calculated as , leading up to . is chosen so that the segment with the shortest parameter range (e.g., chord distance) is split into a minimum number of divisions (e.g., = 100), and segments with longer chord length are split into proportionally larger numbers of divisions.
While minimizing in Eq. (6), it is also important to impose the correct boundary conditions for , , and . This is needed in order to maintain the continuity of the axis velocity and acceleration profiles at the transition points between adjacent toolpath segments and , which are to be interpolated with their respective feed correction polynomials.
where denotes the position vector. and are the derivatives with respect to the spline parameter , which can be obtained by differentiating, e.g., Eq. (1). Subscripts and , in the notation , , designate the kth and k + 1st toolpath segments, respectively. Thus, it can be validated that the toolpaths are second order continuous (C2) with respect to the spline parameter, both within each individual segment, and also at the connection points between the adjacent segments.
Considering Eq. (8), if the commanded feedrate (i.e., tangential velocity ) and tangential acceleration () are continuous throughout the whole toolpath, then the required condition for axis level velocity and acceleration profiles to remain continuous is that the derivatives and also remain continuous.
The coefficients ,…, are computed by substituting the expressions for , , , from Eq. (1) and collecting the terms that multiply different powers of . If the toolpath is cubic, then is only fourth degree.
Examining Eqs. (9) and (10), it can be seen that due to the conditions in Eq. (7) being satisfied, when a toolpath is interpolated by solving the exact value of the spline parameter () which yields the desired arc displacement (), and the commanded and are continuous, the axis velocity and acceleration profiles in Eq. (8) will always remain continuous, even at the segment crossings. However, such an interpolation calls for a very accurate solution technique, such as the iterative solutions reported in Refs. [3,15].
In the practical implementation in Ref. [3], to avoid numerical conditioning problems, the expressions for and were normalized before the computations, and the spline coefficients denormalized afterward. In this paper, the same computer code has been reused for the benchmark study presented in Sec. 4. The normalization steps, for brevity reasons, have not been shown in the preceding review of the earlier method for fitting the feed correction polynomial.
Proposed New Solution
zeroth, first, and second derivative boundary conditions are already known, as explained in Sec. 2. Hence, solving and would enable the feed correction polynomial to be completely parameterized.
Upon solving and , the feed correction polynomial coefficients ,…, are calculated using Eq. (15).
Implementation Results
The proceeding implementation results present the following:
- (1)
correctness of the formulation and implementation, by comparing the result with that of the earlier solution in Ref. [3], using a fan shaped toolpath which had also been used in Ref. [3]
- (2)
justification of using a seventh order feed correction polynomial, as opposed to a lower order (such as fifth)
- (3)
arithmetic computational load and computational time comparison with the earlier solution in Ref. [3]
- (4)
an additional toolpath example (based on an airfoil shape), which further demonstrates the effectiveness of the proposed streamlined solution
Numerical Verification of the Proposed Strategy.
The proposed fitting algorithm for the feed correction polynomial was implemented in matlab and benchmarked against some of the earlier methods that were compared in Ref. [3]. For consistency, the comparison was carried out using the well-known clover leaf profile, originally used by Wang et al. [4,5]. This toolpath, shown in Fig. 4(a), had been used to compare different interpolation algorithms in Ref. [3], and provides a suitable benchmark as it contains a good mixture of low and high curvature regions. Rather than applying feedrate modulation, a constant feed of 100 mm/s was commanded, which makes it very easy to spot and compare feedrate fluctuations consistently throughout the toolpath in terms of percentage (%) values. Implementation results for natural interpolation, first order Taylor interpolation, feed correction based interpolation (using the earlier fitting algorithm in Ref. [3]), and feed correction based interpolation (using the new formulation) are presented in Figs. 4(b)–4(e), respectively. All results were obtained using 1 (ms) as the interpolation period. The improvements obtained with Taylor and feed correction based interpolation are clear. Furthermore, it is seen that the new formulation generates virtually identical results with the earlier method in Ref. [3], which was based on the solution of a 14-unknown linear equation system for each toolpath segment. These were solved using the backslash “\” operator in matlab for solving linear equation systems. Figure 4(e) essentially demonstrates the correctness of the new formulation presented in Sec. 3.
Justification of Seventh Order Feed Correction.
As indicated in Sec. 1, seventh degree feed correction polynomial was chosen to guarantee acceleration level continuity in the interpolated spline toolpath, while also successfully approximating the mapping between the spline parameter and actual arc displacement .
The effect of reducing the feed correction polynomial order to a lower degree (i.e., such as fifth order) is demonstrated in Fig. 5. Figure 5(a) shows the feed consistency achieved with the seventh degree case for the fan profile in Fig. 4(a). Figure 5(b) shows the case in which a fifth degree feed correction polynomial is adopted. Of the six degrees of freedom available, four of them were used to match zeroth and first order boundary conditions at the segment connections, and the remaining two to obtain a least squares fit to the generated versus data. The degradation from the seventh order case in Fig. 5(a) is clear. Figure 5(c) shows another implementation of fifth order feed correction, this time in which all six degrees of freedom are used only to match zeroth, first, and second derivative boundary conditions at both ends of each spline segment, thus leaving no extra degree of freedom available for curve-fitting to the generated versus data. In this case, the degradation in feedrate consistency is even more obvious. These results demonstrate the justification of using a seventh order feed correction polynomial for acceleration-continuous spline interpolation with minimal feedrate fluctuation. If jerk continuity is also desired, a higher order feed correction polynomial, such as ninth, may also be used as proposed in Ref. [13]. The solution methodology for on-the-fly fitting, explained in Sec. 3, can be directly extended to the ninth order feed correction polynomial case as well.
Computational Load Analysis and Benchmarking.
Arithmetic operation breakdowns comparing the earlier feed correction polynomial fitting method [3] and the proposed new technique are provided in Tables 1 and 2. The analysis in Table 1 assumes that the system of 14 linear equations in Eq. (14) is solved using lower-upper (LU) decomposition with partial pivoting (i.e., Crout's method) [16]. This is one of the mainstream methods available in numerical analysis, which suits the problem and can also be implemented efficiently in real-time. As mentioned earlier, is the number of divisions in arc length integration, during which versus data is computed. Using the method in Ref. [3], fitting the feed correction polynomial for one spline segment requires arithmetic operations. With the new formulation, the corresponding number of operations has been reduced to . As is typically is chosen to be larger than 100, the new formulation predicts 41–49% reduction in the number of calculations required. It is important to note that the reasons behind this computational load reduction are: (i) the new streamlined solution does not require the use of a dedicated linear equation solver and (ii) it contains significant simplifications which avoid the tedious matrix multiplication steps, earlier required to construct the and terms in Eq. (14).
Equation no. | Term | Multiplications | Divisions | Additions | Subtractions |
---|---|---|---|---|---|
(15) | … | 6 | 1 | 0 | 0 |
(16) | () | 0 | 0 | 0 | |
(17) | , | 2 | 0 | 0 | 0 |
(17) | … | 6 × () | 0 | 0 | 0 |
(17) | 4 × () | 0 | 3 × () | 1 × () | |
(17) | 3 × () | 0 | 1 × () | 2 × () | |
(17) | 28 × () | 0 | 13 × () | 11 × () | |
(18) | 0 | 0 | 0 | 1 × () | |
(19) | 1 × () | 0 | 1 × () | 0 | |
(19) | 1 × () | 0 | 1 × () | 0 | |
(19) | 1 × () | 0 | 1 × () | 0 | |
(19) | 1 × () | 0 | 1 × () | 0 | |
(19) | 1 × () | 0 | 1 × () | 0 | |
(20) | , | 6 | 1 | 0 | 3 |
(15) | … | 42 | 0 | 18 | 12 |
Total operations | 2 |
Equation no. | Term | Multiplications | Divisions | Additions | Subtractions |
---|---|---|---|---|---|
(15) | … | 6 | 1 | 0 | 0 |
(16) | () | 0 | 0 | 0 | |
(17) | , | 2 | 0 | 0 | 0 |
(17) | … | 6 × () | 0 | 0 | 0 |
(17) | 4 × () | 0 | 3 × () | 1 × () | |
(17) | 3 × () | 0 | 1 × () | 2 × () | |
(17) | 28 × () | 0 | 13 × () | 11 × () | |
(18) | 0 | 0 | 0 | 1 × () | |
(19) | 1 × () | 0 | 1 × () | 0 | |
(19) | 1 × () | 0 | 1 × () | 0 | |
(19) | 1 × () | 0 | 1 × () | 0 | |
(19) | 1 × () | 0 | 1 × () | 0 | |
(19) | 1 × () | 0 | 1 × () | 0 | |
(20) | , | 6 | 1 | 0 | 3 |
(15) | … | 42 | 0 | 18 | 12 |
Total operations | 2 |
While computing the results reported for the fan profile in Fig. 4, the minimum number of divisions, corresponding to the shortest spline segment, was chosen as . As a result, the average value of the number of integration divisions per segment was . Based on this number, the total number of arithmetic operations required to implement the earlier feed correction fitting method in Ref. [3], versus the proposed new technique, has been tabulated in bar graph form in Fig. 6(a). In this figure, the abbreviations “M,” “D,” “A,” and “S” designate the required number of multiplication, division, addition, and subtraction operations, respectively. As can be seen, in fitting the 88 feed correction polynomials for the fan profile, the earlier solution requires on average 65,688 operations per segment, whereas the new solution requires 36,967 operations; thus achieving an overall 43.7% reduction in the computational load.
To validate this prediction, computational time benchmarks were carried out for the cases reported in Figs. 4(d) and 4(e). A “for” loop was incorporated into the code, requiring each segment's feed correction polynomial to be solved using both the old and new methods 10,000 times, in order to allow averaging of the measured computational times. The implementation platform was a 2.26 GHz dual core Intel CPU, running matlab r2011a over windows 7 (64-bit) operating system. While implementing the earlier solution in Ref. [3], two different approaches were tested for solving the system of 14 linear equations in Eq. (14). The first approach was using matlab's built-in implementation of LU decomposition with partial pivoting (i.e., Crout's method), by invoking the command “linsolve.” The arithmetic complexity of this method was analyzed in Table 1 and Fig. 6(a). The second approach was to use matlab's backslash \ operator. The linear equation system in Eq. (14) can be expressed as . According to matlab's documentation [17], when the matrix is Hermitian but not positive definite (which was numerically verified to be the case for Eq. (14)), matlab applies the block lower triangular diagonal matrix (LDL') factorization method for Hermitian indefinite matrices, available from the MA57 library of routines. While an arithmetic operation breakdown of the LDL’ factorization method is beyond the scope of this paper, computational time benchmarks are still presented in order to allow comparison with the results from Crout's method and the proposed on-the-fly fitting method.
The third case that was benchmarked corresponds to using the new feed correction polynomial fitting method, proposed in this paper. The average solution times for all three cases in fitting the feed correction polynomial for a single segment are shown in Fig. 6(b). It is seen that compared to using the earlier method in Ref. [3] with Crout's algorithm, which requires 214 μs per spline segment, the new formulation allows this duration to be reduced to 123 μs, indicating 42.5% reduction in the computational time. It is interesting to note that this reduction is very consistent with the 43.7% reduction in mathematical operations, as predicted in the analysis in Fig. 6(a). It is also seen that while switching to the LDL’ algorithm reduces the computational time from 214 to 195 μs (by 8.8%), this reduction is not as significant as that achieved with the proposed new formulation. This is because both implementations of the earlier method still require matrix multiplications to construct the and terms. Overall, these benchmarks validate the improved suitability of the new feed correction polynomial fitting method, over the earlier approach combined with various linear equation solvers.
Additional Implementation Example.
Implementation of the proposed feed correction polynomial fitting method on a second toolpath example is shown in Fig. 7. In this case, an airfoil profile is constructed by scaling the point data, which can be found in Ref. [18], with a factor of 6 × and fitting an approximately arc length parameterized quintic spline, following the method in Ref. [5]. Two cases are presented: natural interpolation (Fig. 7(b)) and interpolation with a feed correction polynomial, fitted according to the proposed new method (Fig. 7(c)). In natural interpolation, significant feed fluctuation occurs near the highest curvature bend, and reaches 0.5%, whereas the feed correction polynomial contains the fluctuation to within 0.1%; thus, further demonstrating the effectiveness of the proposed scheme.
Conclusions
This technical note has presented a new way to fit the feed correction polynomial, which can be used to interpolate along C2 spline toolpaths with minimal feedrate fluctuation in an acceleration-continuous manner. The new approach has a closed-form solution, which is highly suitable for real-time implementation, as opposed to requiring matrix calculations to construct and solve a system of 14 linear equations for each toolpath segment. The equivalence between the old and new solutions has been demonstrated in terms of achieved feedrate accuracy on a well-known toolpath benchmark. By conducting a thorough computational load analysis, accompanied by execution time measurements, it has been demonstrated that the new solution requires, on average, 43.7% less arithmetic operations and 42.5% less computational time.
Acknowledgment
This research was carried out through the support of NSERC granting agency in Canada.