## 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
Fig. 1
Close modal
$x(u)=Axku5+Bxku4+Cxku3+Dxku2+Exku+Fxky(u)=Ayku5+Byku4+Cyku3+Dyku2+Eyku+Fyk} for 0≤u≤lk$
(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
$ds=(dx/du)2+(dy/du)2·du=xu2+yu2·du$
(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
Fig. 2
Close modal
Fig. 3
Fig. 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
$u∧=Akfs7+Bkfs6+…+Hkfwhere: 0≤s≤Sk}$
(3)
Data generated during the numerical integration of the segment arc length ($Sk$) is stored in vectors $u$ and $s$
$u=[0u1u2…uMk]T, s=[0s1s2…sMk]T$
(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$
$u∧=[0u∧1u∧2⋮u∧Mk]=[000…1s17s16s15…1s27s26s25…1⋮⋮⋮⋱⋮sMk7sMk6sMk5…1]︸Φ[AkfBkfCkf⋮Hkf]︸θ=Φθ$
(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
$J=12eTe=12(u-Φθ)T(u-Φθ)$
(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:
$[r(lk)]k=[r(0)]k+1, [ru(lk)]k=[ru(0)]k+1, [ruu(lk)]k=[ruu(0)]k+1$
(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
$r·(t)=drdt=rss·=ruuss·r··(t)=dr·dt=rsss·2+rss··=dds(ruus)︸rss·2+ruuss··=(ruuus2+ruuss)s·2+ruuss··}$
(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:
$us=duds=[P(u)]-1/2uss=d2uds2=-12[P(u)]-2·dPdu(u)$
(9)
For a quintic toolpath, the polynomial $P(u)$ in Eq. (9) takes the form
$P(u)=xu2+yu2=p0u8+p1u7+…+p8Pu(u)=dPdu=2xuxuu+2yuyuu=8p0u7+7p1u6+…+p7$
(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, $u∧s=du∧/ds$ and $u∧ss=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
$usinit=[[P(u)]-1/2]u=0ussinit=[-12[P(u)]-2·dPdu(u)]u=0usfinal=[[P(u)]-1/2]u=lkussfinal=[-12[P(u)]-2·dPdu(u)]u=lk}$
(11)
Compliance with these boundary conditions is expressed in matrix form, together with the parameter range and total arc length constraints (i.e., $s=0⇔u=0$ and $s=Sk⇔u=lk$), by including the zeroth, first, and second order expressions of Eq. (3) evaluated at $u=0$ and $u=lk$
$[000000010000001000000200Sk7Sk6Sk5Sk4Sk3Sk2Sk17Sk66Sk55Sk44Sk33Sk22Sk1042Sk530Sk420Sk312Sk26Sk200]︸L[AkfBkfCkfDkfEkfFkfGkfHkf]︸θ=[0usinitussinitlkusfinalussfinal]︸ξ(12)$
(12)
Hence, the feed correction polynomial is obtained by solving the following optimization problem:
$minθ12(u-Φθ)T(u-Φθ) subject to:Lθ=ξ$
(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:
$ΦTΦθ+LTΛ=ΦTuLθ=ξ}⇒[ΦTΦLTL0][θΛ]=[ΦTuξ]$
(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
$Akf=16Sk4(usssinit+usssfinal)+2Sk5(ussinit-ussfinal)+10Sk6(usinit+usfinal)+20Sk7(uinit-ufinal)Bkf=-16Sk3(4usssinit+3usssfinal)+12Sk4(13ussfinal-15ussinit)-1Sk5(36usinit+34usfinal)+70Sk6(ufinal-uinit)Ckf=1Sk2(usssinit+12usssfinal)+1Sk3(10ussinit-7ussfinal)+1Sk4(45usinit+39usfinal)+84Sk5(uinit-ufinal)Dkf=-16Sk(4usssinit+usssfinal)+52Sk2(ussfinal-2ussinit)-1Sk3(20usinit+15usfinal)+35Sk4(ufinal-uinit)Ekf=16usssinit Fkf=12ussinit Gkf=usinit Hkf=uinit$
(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
$σj=sj/Sk, such that: 0≤σj≤1, 1≤j≤Mk-1$
(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:
$[φi(σ1)φf(σ1)φi(σ2)φf(σ2)⋮⋮φi(σMk-1)φf(σMk-1)]︸Φ[usssinitusssfinal]︸θ=[u1-u˜(σ1)u2-u˜(σ2)⋮uMk-1-u˜(σMk-1)]︸Y⇒θ=(ΦTΦ)-1ΦTY(18)$
(18)
The terms $(ΦTΦ)-1$ and $ΦTY$ in Eq. (18) can be constructed as follows:
$ΦTΦ=[∑j=1Mk-1φi(σj)φi(σj)∑j=1Mk-1φi(σj)φf(σj)∑j=1Mk-1φf(σj)φi(σj)∑j=1Mk-1φf(σj)φf(σj)]=[ΣiiΣifΣifΣff], (ΦTΦ)-1=1ΣiiΣff-Σif2[Σff-Σif-ΣifΣii],ΦTY=[∑j=1Mk-1φi(σj)yj∑j=1Mk-1φf(σj)yj]=[ΣiyΣfy](19)$
(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)
$[usssinitusssfinal]=1∑ii∑ff-∑if2[ΣffΣiy-ΣifΣfy-ΣifΣiy+ΣiiΣfy]$
(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
Fig. 4
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
Fig. 5
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

(5) and (14)$ΦTΦ$64 × $(Mk+1)$064 × $Mk$0
(12)$Sk2 …$$Sk7$6000
(12)$L$ or $LT$11000
(5) and (14)$ΦTu$8 × $(Mk+1)$08 × $Mk$0
(14)$[ΦTΦLTL0][θΛ]=[ΦTuξ]$11974101001
Total operations$72Mk+1286$41$72Mk$1001
(5) and (14)$ΦTΦ$64 × $(Mk+1)$064 × $Mk$0
(12)$Sk2 …$$Sk7$6000
(12)$L$ or $LT$11000
(5) and (14)$ΦTu$8 × $(Mk+1)$08 × $Mk$0
(14)$[ΦTΦLTL0][θΛ]=[ΦTuξ]$11974101001
Total operations$72Mk+1286$41$72Mk$1001
Table 2

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

(15)$1/Sk$$1/Sk7$6100
(16)$σj$($Mk-1$)000
(17)$Sk2$,$Sk3$2000
(17)$σj2$$σj7$6 × ($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)$yj$0001 × ($Mk-1$)
(19)$Σii$1 × ($Mk-1$)01 × ($Mk-1$)0
(19)$Σif$1 × ($Mk-1$)01 × ($Mk-1$)0
(19)$Σff$1 × ($Mk-1$)01 × ($Mk-1$)0
(19)$Σiy$1 × ($Mk-1$)01 × ($Mk-1$)0
(19)$Σfy$1 × ($Mk-1$)01 × ($Mk-1$)0
(20)$usssinit$,$usssfinal$6103
(15)$Akf$$Hkf$4201812
Total operations$47Mk+9$2$22Mk-4$$15Mk$
(15)$1/Sk$$1/Sk7$6100
(16)$σj$($Mk-1$)000
(17)$Sk2$,$Sk3$2000
(17)$σj2$$σj7$6 × ($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)$yj$0001 × ($Mk-1$)
(19)$Σii$1 × ($Mk-1$)01 × ($Mk-1$)0
(19)$Σif$1 × ($Mk-1$)01 × ($Mk-1$)0
(19)$Σff$1 × ($Mk-1$)01 × ($Mk-1$)0
(19)$Σiy$1 × ($Mk-1$)01 × ($Mk-1$)0
(19)$Σfy$1 × ($Mk-1$)01 × ($Mk-1$)0
(20)$usssinit$,$usssfinal$6103
(15)$Akf$$Hkf$4201812
Total operations$47Mk+9$2$22Mk-4$$15Mk$

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 $Mk≅440$. 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
Fig. 6
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.

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
Fig. 7
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
,”
,
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
.