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 (u), which they are parameterized with respect to, and the actual arc length (s) that is traveled, as the spline parameter u is incremented.

Figure 1 shows a quintic (fifth degree) spline toolpath, fit in x and y axes, formulated in the form

Fig. 1
Spline toolpath [3]
(1)
where Axk, Bxk, , Fxk and Ayk, Byk,…, Fyk are the spline coefficients, u is the spline parameter, and lk is the parameter range for the kth segment. As the spline parameter u is incremented, this results in an arc displacement s along the toolpath. The differential relationship between these two variables is
(2)

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 u to be equivalent to s 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 u is interpolated in uniform increments, the discrepancy between u and s exhibits itself in the form of two detrimental effects [3]:

  1. (i)

    feedrate fluctuations within each segment;

  2. (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 (ui), by expanding it into a Taylor series as a function of time, around its value used in the previous interpolation cycle (ui-1). 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 u and s with the so-called “feed correction polynomial,” as proposed in Ref. [3]. Typically, seventh degree is used, which allows the six boundary conditions on u, us (=du/ds), and uss (=d2u/ds2) 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 (u) versus arc length (s) data, which is generated earlier on during the numerical integration of the total arc length (Sk) for each spline segment. Fitting of the feed correction polynomial is shown in Fig. 3.

Fig. 2
Three commonly used spline interpolation methods [3]: (a) natural interpolation,
                            (b) interpolation using Taylor series expansion, and
                            (c) interpolation with feedrate correction
                        polynomial
Fig. 2
Three commonly used spline interpolation methods [3]: (a) natural interpolation,
                            (b) interpolation using Taylor series expansion, and
                            (c) interpolation with feedrate correction
                        polynomial
Close modal
Fig. 3
Parameterization of the feed correction polynomial using data generated from
                        the numerical integration of the arc length (s) over the
                        spline parameter (u) [3]
Fig. 3
Parameterization of the feed correction polynomial using data generated from
                        the numerical integration of the arc length (s) over the
                        spline parameter (u) [3]
Close modal

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 u and s, polynomial function mapping, or real-time iterative solution, target accurate achievement of the desired arc displacement (s), and therefore the correct feedrate ds/dt 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

The seventh degree feed correction polynomial is formulated as
(3)
Data generated during the numerical integration of the segment arc length (Sk) is stored in vectors u and s
(4)

Above, sj (j = 1,…, Mk) is the numerically integrated arc displacement corresponding to the spline parameter value uj. Mk is the number of divisions, which are distributed uniformly along [0,uMk]. For each spline parameter value uj=j·Δu, the corresponding increments in axis positions (Δxj,Δyj) are calculated using the toolpath definition in Eq. (1) and the resulting arc displacement is calculated as sj=sj-1+Δxj2+Δyj2, leading up to sMk=Sk. Mk 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., Mk,min = 100), and segments with longer chord length are split into proportionally larger numbers of divisions.

In the real-time application, to achieve a desired arc displacement (s), the feed correction polynomial is used to find the most suitable value for the spline parameter (u). It is formulated by conducting a least squares fit of u over s. Defining the vector (u) of predicted spline parameter values for individual elements sj within s
(5)
the prediction error vector (e) for the spline parameter becomes e=u-u=u-Φθ. In finding the set of feed correction polynomial parameters θ, the objective function to be minimized is
(6)

While minimizing J in Eq. (6), it is also important to impose the correct boundary conditions for u, us, and uss. This is needed in order to maintain the continuity of the axis velocity and acceleration profiles at the transition points between adjacent toolpath segments k and k+1, which are to be interpolated with their respective feed correction polynomials.

Cubic and quintic spline toolpaths are typically parameterized to satisfy the following continuity conditions at the segment boundaries:
(7)

where r(u)=[x(u),y(u)]T denotes the position vector. ru and ruu are the derivatives with respect to the spline parameter u, which can be obtained by differentiating, e.g., Eq. (1). Subscripts k and k+1, in the notation []k, []k+1, 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.

For the whole toolpath, the axis level velocity and acceleration profiles can be expressed as
(8)

Considering Eq. (8), if the commanded feedrate (i.e., tangential velocity s·) and tangential acceleration (s··) are continuous throughout the whole toolpath, then the required condition for axis level velocity and acceleration profiles to remain continuous is that the derivatives us and uss also remain continuous.

From Eqs. (1) and (2), the expressions for us and uss can be obtained as follows:
(9)
For a quintic toolpath, the polynomial P(u) in Eq. (9) takes the form
(10)

The coefficients p0,…, p8 are computed by substituting the expressions for xu=dx/du, yu=dy/du, xuu=d2x/du2, yuu=d2y/du2 from Eq. (1) and collecting the terms that multiply different powers of u. If the toolpath is cubic, then P(u) 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 (u) which yields the desired arc displacement (s), and the commanded s· and s·· 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].

On the other hand, using the output (u) of the feed correction polynomial instead of the true value of u is only an approximate solution. Considering the high order of the feed correction polynomial in Eq. (3), it can be guaranteed that within each toolpath segment, us=du/ds and uss=du/ds will remain continuous. However, at the segment connections, it is vital that second-order continuity conditions be imposed during the fitting of the feed correction polynomial. Hence, for each toolpath segment, the initial and final derivative boundary conditions are computed by evaluating the derivative expressions in Eq. (9) at the beginning and end of that segment
(11)
Compliance with these boundary conditions is expressed in matrix form, together with the parameter range and total arc length constraints (i.e., s=0u=0 and s=Sku=lk), by including the zeroth, first, and second order expressions of Eq. (3) evaluated at u=0 and u=lk
(12)
Hence, the feed correction polynomial is obtained by solving the following optimization problem:
(13)
In Ref. [3], it was shown that this linear equality constrained quadratic minimization problem can be solved using Lagrange multipliers method, which transforms into the problem of solving the following system of 14 linear equations, containing the feed correction polynomial coefficients in θ (8 × 1) and Lagrange multipliers Λ (6 × 1) for the six boundary condition constraints:
(14)

In the practical implementation in Ref. [3], to avoid numerical conditioning problems, the expressions for u and s 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

Denoting the zeroth, first, second, and third derivative boundary conditions at both ends of the seventh order feed correction polynomial as uinit (= 0), ufinal (= lk), usinit, ussinit, usssinit, usfinal, ussfinal, usssfinal, and the segment total arc length as Sk, the feed correction polynomial coefficients in Eq. (3) can be solved by the applying the Gaussian elimination method as
(15)

zeroth, first, and second derivative boundary conditions are already known, as explained in Sec. 2. Hence, solving usssinit and usssfinal would enable the feed correction polynomial to be completely parameterized.

Defining the normalized arc distance as
(16)
substituting the coefficient expressions from Eq. (15) into Eq. (3), and applying the normalization in Eq. (16), the spline parameter prediction u in Eq. (3) can now be expressed as a function of usssinit, usssfinal, and a remainder term (u˜)
(17)
Equations (15) and (17) can be conveniently verified using a Symbolic mathematics package, such as matlab's Symbolic Math Toolbox. The solution of usssinit and usssfinal can now be formulated as a two-parameter Least Squares estimation problem, as shown in the following equation:
(18)
The terms (ΦTΦ)-1 and ΦTY in Eq. (18) can be constructed as follows:
(19)
where the term yj corresponds to the jth element of vector Y, defined in Eq. (18) as yj=uj-u˜(σj). Thus, usssinit and usssfinal are solved per Eq. (20), which is constructed with the “Σ” terms calculated from Eq. (19)
(20)

Upon solving usssinit and usssfinal, the feed correction polynomial coefficients Akf,…, Hkf are calculated using Eq. (15).

Implementation Results

The proceeding implementation results present the following:

  1. (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. (2)

    justification of using a seventh order feed correction polynomial, as opposed to a lower order (such as fifth)

  3. (3)

    arithmetic computational load and computational time comparison with the earlier solution in Ref. [3]

  4. (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.

Fig. 4
Implementation comparing the new parameterization method for the feed
                            correction polynomial with the earlier results reported in Ref. [3]. (a) 88-Segment
                            toolpath, (b) natural interpolation,
                                (c) first order Taylor, (d) feed
                            correction polynomial (earlier solution, using matrix calculations
                                [3]), and (e)
                            feed correction polynomial (computed with the proposed solution).
Fig. 4
Implementation comparing the new parameterization method for the feed
                            correction polynomial with the earlier results reported in Ref. [3]. (a) 88-Segment
                            toolpath, (b) natural interpolation,
                                (c) first order Taylor, (d) feed
                            correction polynomial (earlier solution, using matrix calculations
                                [3]), and (e)
                            feed correction polynomial (computed with the proposed solution).
Close modal

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 u and actual arc displacement s.

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 u versus s 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 u versus s 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.

Fig. 5
Comparison between seventh order feed correction polynomial and fifth
                            order variations (i.e., justification of using seventh order):
                                (a) seventh order feed correction polynomial
                            parameterized considering u–s data, as well as zeroth,
                            first, and second order boundary conditions (proposed),
                                (b) fifth order feed correction polynomial
                            parameterized considering u–s data and only zeroth and
                            first order boundary conditions, and (c) fifth order
                            feed correction polynomial parameterized considering only zeroth, first,
                            and second order boundary conditions
Fig. 5
Comparison between seventh order feed correction polynomial and fifth
                            order variations (i.e., justification of using seventh order):
                                (a) seventh order feed correction polynomial
                            parameterized considering u–s data, as well as zeroth,
                            first, and second order boundary conditions (proposed),
                                (b) fifth order feed correction polynomial
                            parameterized considering u–s data and only zeroth and
                            first order boundary conditions, and (c) fifth order
                            feed correction polynomial parameterized considering only zeroth, first,
                            and second order boundary conditions
Close modal

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, Mk is the number of divisions in arc length integration, during which u versus s data is computed. Using the method in Ref. [3], fitting the feed correction polynomial for one spline segment requires 144Mk+2328 arithmetic operations. With the new formulation, the corresponding number of operations has been reduced to 84Mk+7. As Mk 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 ΦTΦ and ΦTu terms in Eq. (14).

Table 1

Breakdown of calculations required for the earlier proposed feed correction method [3], using LU decomposition with partial pivoting (i.e., Crout's method) [16] to solve the system of 14 linear equations

Equation no.TermMultiplicationsDivisionsAdditionsSubtractions
(5) and (14)ΦTΦ64 × (Mk+1)064 × Mk0
(12)Sk2Sk76000
(12)L or LT11000
(5) and (14)ΦTu8 × (Mk+1)08 × Mk0
(14)[ΦTΦLTL0][θΛ]=[ΦTuξ]11974101001
Total operations72Mk+12864172Mk1001
Equation no.TermMultiplicationsDivisionsAdditionsSubtractions
(5) and (14)ΦTΦ64 × (Mk+1)064 × Mk0
(12)Sk2Sk76000
(12)L or LT11000
(5) and (14)ΦTu8 × (Mk+1)08 × Mk0
(14)[ΦTΦLTL0][θΛ]=[ΦTuξ]11974101001
Total operations72Mk+12864172Mk1001
Table 2

Breakdown of calculations required to implement the proposed on-the-fly feed correction method

Equation no.TermMultiplicationsDivisionsAdditionsSubtractions
(15)1/Sk1/Sk76100
(16)σj(Mk-1)000
(17)Sk2,Sk32000
(17)σj2σj76 × (Mk-1)000
(17)φi(σj)4 × (Mk-1)03 × (Mk-1)1 × (Mk-1)
(17)φf(σj)3 × (Mk-1)01 × (Mk-1)2 × (Mk-1)
(17)u˜(σj)28 × (Mk-1)013 × (Mk-1)11 × (Mk-1)
(18)yj0001 × (Mk-1)
(19)Σii1 × (Mk-1)01 × (Mk-1)0
(19)Σif1 × (Mk-1)01 × (Mk-1)0
(19)Σff1 × (Mk-1)01 × (Mk-1)0
(19)Σiy1 × (Mk-1)01 × (Mk-1)0
(19)Σfy1 × (Mk-1)01 × (Mk-1)0
(20)usssinit,usssfinal6103
(15)AkfHkf4201812
Total operations47Mk+9222Mk-415Mk
Equation no.TermMultiplicationsDivisionsAdditionsSubtractions
(15)1/Sk1/Sk76100
(16)σj(Mk-1)000
(17)Sk2,Sk32000
(17)σj2σj76 × (Mk-1)000
(17)φi(σj)4 × (Mk-1)03 × (Mk-1)1 × (Mk-1)
(17)φf(σj)3 × (Mk-1)01 × (Mk-1)2 × (Mk-1)
(17)u˜(σj)28 × (Mk-1)013 × (Mk-1)11 × (Mk-1)
(18)yj0001 × (Mk-1)
(19)Σii1 × (Mk-1)01 × (Mk-1)0
(19)Σif1 × (Mk-1)01 × (Mk-1)0
(19)Σff1 × (Mk-1)01 × (Mk-1)0
(19)Σiy1 × (Mk-1)01 × (Mk-1)0
(19)Σfy1 × (Mk-1)01 × (Mk-1)0
(20)usssinit,usssfinal6103
(15)AkfHkf4201812
Total operations47Mk+9222Mk-415Mk

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 Mk,min=100. As a result, the average value of the number of integration divisions per segment was Mk440. 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.

Fig. 6
(a) Computational load and (b) average
                            computational time comparison between the earlier solution [3] and proposed new method for
                            computing the feed correction polynomial for one segment. The reported
                            results are based on the 88-segment toolpath in Fig. 4(a), which was
                            split, on average, to 440 divisions per segment
                                    (= Mk)
                            during the computation of the segment arc-lengths and feed correction
                            polynomials.
Fig. 6
(a) Computational load and (b) average
                            computational time comparison between the earlier solution [3] and proposed new method for
                            computing the feed correction polynomial for one segment. The reported
                            results are based on the 88-segment toolpath in Fig. 4(a), which was
                            split, on average, to 440 divisions per segment
                                    (= Mk)
                            during the computation of the segment arc-lengths and feed correction
                            polynomials.
Close modal

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 Ax=b. According to matlab's documentation [17], when the matrix A 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 ΦTΦ and ΦTu 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 C3 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.

Fig. 7
Implementation of the proposed feed correction method on an airfoil
                            profile. (a) 48-Segment toolpath, (b)
                            natural interpolation, and (c) feed correction
                            polynomial (computed with the proposed solution).
Fig. 7
Implementation of the proposed feed correction method on an airfoil
                            profile. (a) 48-Segment toolpath, (b)
                            natural interpolation, and (c) feed correction
                            polynomial (computed with the proposed solution).
Close modal

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.

References

1.
Li
,
Z.-L.
, and
Zhu
,
L.-M.
,
2014
, “
Envelope Surface Modeling and Tool Path Optimization for Five-Axis Flank Milling Considering Cutter Runout
,”
ASME J. Manuf. Sci. Eng.
,
136
(
4
), p.
041021
.10.1115/1.4027415
2.
Pan
,
Y.
,
Zhou
,
C.
,
Chen
,
Y.
, and
Partanen
,
J.
,
2014
, “
Multitool and Multi-Axis Computer Numerically Controlled Accumulation for Fabricating Conformal Features on Curved Surfaces
,”
ASME J. Manuf. Sci. Eng.
,
136
(
3
), p.
031007
.10.1115/1.4026898
3.
Erkorkmaz
,
K.
, and
Altintas
,
Y.
,
2005
, “
Quintic Spline Interpolation With Minimal Feed Fluctuation
,”
ASME J. Manuf. Sci. Eng.
,
127
(
2
), pp.
339
349
.10.1115/1.1830493
4.
Wang
,
F.-C.
, and
Yang
,
D. C. H.
,
1993
, “
Nearly Arc-Length Parameterized Quintic-Spline Interpolation for Precision Machining
,”
Comput. Aided Des.
,
25
(
5
), pp.
281
288
.10.1016/0010-4485(93)90085-3
5.
Wang
,
F.-C.
,
Wright
,
P. K.
,
Barsky
,
B. A.
, and
Yang
,
D. C. H.
,
1999
, “
Approximately Arc-Length Parameterized C3 Quintic Interpolatory Splines
,”
ASME J. Mech. Des.
,
121
(
3
), pp.
430
439
.10.1115/1.2829479
6.
Farouki
,
R. T.
, and
Shah
,
S.
,
1996
, “
Real-Time CNC Interpolators for Pythagorean-Hodograph Curves
,”
Comput. Aided Geom. Des.
,
13
(
7
), pp.
583
600
.10.1016/0167-8396(95)00047-X
7.
Farouki
,
R. T.
,
Al-Kandari
,
M.
, and
Sakkalis
,
T.
,
2002
, “
Hermite Interpolation by Rotation-Invariant Spatial Pythagorean-Hodograph Curves
,”
Adv. Comput. Math.
,
17
(
4
), pp.
369
383
.10.1023/A:1016280811626
8.
Huang
,
J.-T.
, and
Yang
,
D. C. H.
,
1992
, “
Precision Command Generation for Computer Controlled Machines
,”
Precision Machining: Technology and Machine Development and Improvement
, Vol. PED 58,
ASME
,
New York
, pp.
89
104
.
9.
Shpitalni
,
M.
,
Koren
,
Y.
, and
Lo
,
C.-C.
,
1994
, “
Realtime Curve Interpolators
,”
Comput. Aided Des.
,
26
(
11
), pp.
832
838
.10.1016/0010-4485(94)90097-3
10.
Otsuki
,
T.
,
Kozai
,
H.
, and
Wakinotani
,
Y.
,
1998
, “
Free-Form Curve Interpolation Method and Apparatus
,” U.S. Patent No. 5,815,401.
11.
Lin
,
R.-S.
,
2000
, “
Real-Time Surface Interpolator for 3-D Parametric Surface Machining on 3-Axis Machine Tools
,”
Int. J. Mach. Tools Manuf.
,
40
(
10
), pp.
1513
1526
.10.1016/S0890-6955(00)00002-X
12.
Heng
,
M.
, and
Erkorkmaz
,
K.
,
2010
, “
Design of a NURBS Interpolator With Minimal Feed Fluctuation and Continuous Feed Modulation Capability
,”
Int. J. Mach. Tools Manuf.
,
50
(
3
), pp.
281
293
.10.1016/j.ijmachtools.2009.11.005
13.
Yuen
,
A.
,
Zhang
,
K.
, and
Altintas
,
Y.
,
2013
, “
Smooth Trajectory Generation for Five-Axis Machine Tools
,”
Int. J. Mach. Tools Manuf.
,
71
, pp.
11
19
.10.1016/j.ijmachtools.2013.04.002
14.
Lei
,
W. T.
,
Sung
,
M. P.
,
Lin
,
L. Y.
, and
Huang
,
J. J.
,
2007
, “
Fast Real-Time NURBS Path Interpolation for CNC Machine Tools
,”
Int. J. Mach. Tools Manuf.
,
47
(
10
), pp.
1530
1541
.10.1016/j.ijmachtools.2006.11.011
15.
Cheng
,
C. W.
,
Tsai
,
M. C.
, and
Maciejowski
,
J.
,
2006
, “
Feedrate Control for Non-Uniform Rational B-Spline Motion Command Generation
,”
Proc. Inst. Mech. Eng., Part B
,
220
(
11
), pp.
1855
1861
.10.1243/09544054JEM401
16.
Press
,
W. H.
,
Flannery
,
B. P.
,
Teukolsky
,
S. A.
, and
Vetterling
,
W. T.
,
1992
,
Numerical Recipes in C: The Art of Scientific Computing
, 2nd ed.,
Cambridge University Press
,
New York
.
17.
MathWorks
, “
MATLAB R2014b Documentation
,” MathWorks, Inc., Natick, MA, http://www.mathworks.com/help/matlab/ref/mldivide.html
18.
Altintas
,
Y.
,
2000
,
Manufacturing Automation: Metal Cutting Mechanics, Machine Tool Vibrations, and CNC Design
,
Cambridge University Press
,
Cambridge, UK
.