Linear programming feedrate optimization

This paper focuses on two aspects of feedrate optimization via linear programming methods. Namely, the effect of curve sampling on time optimality of the resultant feedrate profile and a method of feedrate profile adaptation in response to a feedrate override command. A comparison of three distinct curve sampling approaches (uniform in parameter, uniform in arc length and curvature adaptive) is performed on a series of standard tool path curves. Results show that the curvature-adaptive sampling approach leads to substantial machining time reduction for tool path curves displaying high degree of curvature variation. Secondly, a method by which a new feedrate profile can be calculated in response to a feedrate override command is developed. The method formulates a new set of boundary conditions on the control point sequence of the feedrate curve in such a way that the resulting profile is guaranteed to coincide with the currently active profile up to the moment of override command, while minimizing the arc length necessary for transition to the newly commanded feedrate.


L
Total length of toolpath r ′ , r ′′ , r ′′′ First, second and third derivative of toolpath vector with respect to arc length u Toolpath curve parameter v max Tangential velocity limit (feedrate) a max Tangential acceleration limit j max Tangential jerk limit v max Vector of axial velocity limits a max Vector of axial acceleration limits j max Vector of axial jerk limits s Arc length parameteṙ s Tangential velocity (feedrate) s Tangential acceleration ⃛ s Tangential jerk q(s) Squared feedrate expressed as a cubic B-spline d Order of B-spline N i, 2 i-th second order B-spline basis function a Vector of B-spline control points c Uniformly distributed vector of weightŝ q Pseudojerk curvê M Matrix of constraints of LP problem without considering jerk constraints M Matrix of constraints of LP problem with jerk constraints included

Introduction
Interpolator is a component of the computer numerical control (CNC) system, which governs axial servo drives by providing position commands. Vectors of axis positions are being sent to servo drives at specified frequency. In a precision and high-speed machining environment, interpolator is required to optimize axial movement, so that maximum accuracy and shortest machining time are achieved. Considerable effort has been put by researchers into developing algorithms capable of interpolating tool paths which are determined by parametric curves such as NURBS]NURBSNon-Uniform Rational B-Spline. NURBS curves (and surfaces) are a staple in computer assisted design due to their compact representation and efficient evaluation. The necessity of generating feedrate profiles on NURBS tool paths has led to the development of a variety of distinctive interpolation methods.
To create efficient feedrate profiles for NURBS curve tool paths, it is necessary to consider the relationship between the curve's parameter and its arc length. This relationship does not have a closed form solution in general and disregarding it can lead to undesirable feedrate oscillations. This problem was typically solved using Taylor series expansion (see, e.g. [1]). This, of course, introduces errors caused by omitting higher order terms. As a possible solution to this problem, the so-called Pythagorean hodograph curves have been proposed as an alternative to NURBS curve representations (for a comprehensive study of Pythagorean hodograph curves see [2]). Pythagorean hodograph curves enjoy the nice property of polynomial dependence of arc length on curve parameter, making them particularly useful for CNC interpolations. However, due to the widespread use of NURBS curves in CAD systems, Pythagorean hodograph curves still have not found wide application in common practice.
In order to solve the problem of parameter-to-arc length relationship in NURBS curves, it is necessary to devise approximation methods. In [3], the authors develop a recursive approximation algorithm to interpolate quintic spline tool paths with constant distance increment. The feedrate profile is constructed by varying the interpolation period, which is then reconstructed in the control loop period by fitting a fifth-degree polynomial. Another method that approximates the relationship between the spline parameter and its arc length is presented in [4]. The main idea is to use a seventh-degree polynomial (feed correction polynomial) to approximate this relationship. The coefficients of this polynomial are estimated using least squares. This initial approximation is then further refined using the Newton-Raphson iterative method. This method has been further modified in [5], where the parameter-to-arc length relationship is first approximated for a series of discrete points using the adaptive Simpson rule and this discrete representation is then used to fit a series of feedrate polynomials depending on a maximal allowed deviation. Due to the adaptive nature of this algorithm, it can approximate the parameter-to-arc length relationship with arbitrary precision. This particular algorithm was used for the purposes of this paper.
The most common approaches to feedrate optimization found in current literature are: composing the feedrate profile by connecting S-shaped feedrate profiles with piecewise constant jerk (constructed analytically or by using FIR] FIRFinite Impulse Response filters), brute force optimization methods and utilization of linear programming (LP). The last mentioned approach combines computational efficiency with the ability to construct nearly time-optimal feedrate profiles. First presented in [6], this method of feedrate optimization is based on the use of tool path discretization followed by reformulation of the feedrate optimization problem where the feedrate function is represented by a B-spline curve whose control points serve as the free variables. Finding the control points poses a non-linear optimization problem, which requires the introduction of the so-called pseudojerk, that serves as an upper bound to the jerk inequality, thus allowing to linearize the optimization problem. This upper bound is precomputed using only axial velocity and acceleration limits. The proximity of this solution to true optimal feedrate profile is limited by the finite number of control points for feedrate profile curve and the pseudo-jerk relaxation. This method was further expanded on in [7], where a windowing based parallelization method was developed in order to further improve the computational efficiency of the original algorithm.
It is of note, that the methods of linear programming have also been applied to the feedrate optimization problem in an earlier article [8], although in a substantially different way. The feedrate profile on a spline toolpath was defined as a function of time and was obtained by minimizing a square integral of jerk via manipulation of time durations of selected path segments. This solution did not, however, consider the influence of the parameter-to-arc length relationship. For a similar approach using time parametrization, see [9].
Another computationally efficient method to construct the feedrate profile for NURBS tool paths is based on connecting piecewise jerk constant S-curve type feedrate transition segments. Feedrate can be modulated in this way to comply with kinematic limits which are evaluated at a fixed set points (see [10]). Similarly, in [5] feedrate profiles for spline tool paths are generated by connecting constant feedrate segments by S-shaped curve transition segments according to a heuristic algorithm. Axis jerk and acceleration limits are evaluated at knot locations. For a given portion of the spline trajectory, which is bounded by two knots, the maximum feedrate is calculated based on these constraints. The algorithm performs a search for maximum feedrate that is both achievable and respects constraints in all portions of the spline trajectory.
In [11], an axis movement smoothing algorithm for 5-axis milling is developed, which utilizes a heuristic algorithm to conform to contour tolerance limits. In [12], the authors describe a look-ahead algorithm, including analytical expressions, for interpolating linear segments. The article [13] presents a corner smoothing technique for five-axis machining using micro splines inserted between consecutive linear blocks with synchronized position and tool orientation.
A method of bi-directional scan was presented in [14,15], and applied in multiple subsequent works [16][17][18]. As the name suggests, the main idea of the method lies in combination of a forward and backward pass (or scan) of the toolpath. In the first scan, the toolpath is passed in the forward direction while a feedrate profile is iteratively constructed by applying a greedy algorithm that minimizes total cycle time. This algorithm, however, can lead to a feedrate profile with higher order discontinuities. In the second backward scan, these discontinuities are detected and the profile is updated by relaxing its kinematic properties in order to achieve the required order of smoothness. This method can theoretically produce time-optimal trajectories. Its optimality and computation requirements do, however, depend strongly on the length of the parameter step applied in the scanning process.
An alternative approach to feedrate optimization is to adaptively modulate the feedrate based on constraints evaluated by a virtual process model. This approach was applied in [19] to 5-axis flank milling of a jet engine impeller, where feed was adaptively modulated based on multiple constraints such as tool deflection, maximum chip load, torque limits and tool shank bending stress.
Yet another approach to feedrate optimization is to consider the problem as a problem of reparametrization, where the relation between the time and curve parameters is defined by a spline function whose parameters are optimized in order to minimize overall jerk (see [20] and [21]). In [21] this method is specifically applied to the problem of multi-axis flank milling. This type of milling presents its unique set of challenges when compared to end-point milling, such as free-form surface approximation by sweep patches [22] and cutting force modeling [23][24][25]. Flank milling is an interesting alternative to end-point milling because of its higher material removal efficiency and because it is, by definition, scallop free on every sweep patch. For a detailed list of relevant literature on flank milling, see [26].
An alternative approach to generating feedrate profiles relies on the application of sequences of FIR filters. The authors of [27] show that applying a sequence of moving average filters to the feedrate impulse produces a feedrate function that is equivalent to the analytically expressed S-shaped feedrate profile. Their work provides a relationship between feedrate profile and its frequency domain, which can be utilized to suppress vibrations on given resonant frequencies by adjusting filter parameters. In [28], analytically generated acceleration limited feedrate profile is combined with the FIR approach to generate feedrate profile for arbitrary velocity and acceleration conditions. Finally, [29] extends this method by including circular and linear toolpath blending capability with confined contour error. A transition between linear and circular segments can thus be performed at nonzero feedrates by controlling the convolution overlap time of two consecutive feed impulses. This method was recently extended to 5-axis machining in [30].
Iterative method, referred to as VPOp, can be applied to interpolate NURBS toolpaths even in 5-axis milling. This method achieves short (shortest among all the mentioned approaches) machining times while respecting axis acceleration and jerk constraints. However, its high computational demands severely limit its real world application. Articles [31] and [32] use the VPOp algorithm to interpolate directly into the parametric space of the surface of the workpiece CAD model. This original approach removes the issue of CAM tolerance and interpolation tolerance stacking up.
A common method of machining free form curves(and surfaces) is to generate a G-code program by using CAM software. Such programs often consist of a sequence of line commands which approximate the original curve. Although there are methods to optimize CAM block generation with respect to speed and accuracy [33] in this work we focus on methods of interpolating free form curves directly.
To summarize the above comparison, linear programming methods produce results which are close to optimal even on complex curves. Given a sufficient number of discretization points, the method is able to respect all the given kinematic limits while maintaining computational efficiency. The biggest disadvantage is the necessity of linearization of the jerk constraint, which results in loss of optimality of the solution. Nevertheless, when compared to other methods the conciseness of the mathematical formulation, ability to simultaneously incorporate both tangential and axial limits on kinematic variables and the possibility of parallelization all make linear programming methods a favorable choice.
This paper is primarily focused on two research topics: firstly, the effects of several types of spline sampling on the time optimality of a feedrate function obtained via a linear programming feedrate optimization approach are studied and demonstrated with the use of several testing curves. The obtained results demonstrate the potential advantages of adaptive sampling methods on the time optimality of resulting feedrate profiles. Secondly, an original method of feedrate override (again in the context of linear programming feedrate optimization) is presented and its implementation details are discussed. This paper is organized as follows: Section 2 reviews the formulation of the linear programming optimization problem (Section 2.1) as well as knot vector construction and evaluation point selection given a general sampling of the toolpath curve (Section 2.2).
In Section 3, three types of curve sampling (equidistant in curve parameter, equidistant in arc length and adaptive) are presented and the effects of sampling approach on time optimality of the resultant feedrate are discussed based on a comparison for several test curves. The main original result of this section is that an application of a suitable adaptive sampling method can lead to a significant decrease of machining time. The section is concluded with a discussion of limitations and possible future research directions. Section 4 is dedicated to an original method of feedrate override for linear programming feedrate optimization. The method is divided into two distinct cases: override to a higher commanded feedrate (Section 4.2) and override to a lower commanded feedrate (Section 4.3). This method describes how, and under what conditions a new feedrate profile implementing the override command can be obtained given curve parameters and the currently active feedrate profile. The method is suitable for application purposes due to low memory requirements and relative ease of implementation.
The article is concluded with several closing remarks in Section 5 and three appendices: Appendix A (1) providing definitions of the test path curves used in Sections 3 and 4, Appendix B (2) where comparisons of test path curve trajectories interpolated using methods discussed in Section 3 are presented and finally Appendix C (3) in which a definition of an auxiliary function is given.

Feedrate optimization -a linear programming problem
This chapter recalls the formulation of feedrate optimization via linear programming (Section 2.1) while briefly discussing knot vector construction and the selection of evaluation points (Section 2.2).

Problem formulation
Assume that the tool path r is defined as a NURBS curve with at least C 2 continuity (in practice this typically means a curve of order three or five) In order to formulate the relations between velocity, acceleration and jerk along the curve, it is beneficial to introduce the arc length parametrization so that one can also express r as where s denotes the arc length parameter and L represents the total length of the curve. Let v max , a max , j max denote the maximal axial limits of velocity, acceleration and jerk, respectively. Furthermore, let v max , ā max , ̄j max denote the maximal limits of tangential velocity, tangential acceleration and tangential jerk, respectively. The axial velocity, acceleration and jerk can be expressed as: where r ′ , r ′′ and r ′′′ denote the derivatives of r with respect to s, while ̇s , s and ⃛ s denote the tangential velocity (feedrate), tangential acceleration and tangential jerk, respectively. The feedrate optimization problem, i.e. derivation of time-optimal feedrate profile, can then be formulated as: In order to linearize the above problem, the authors of [6] (see also [7]) apply the following substitution: i.e. the square of feedrate is expressed as a cubic B-spline where a = [a 1 , … , a K ] is the vector of control points and [N 1,2 , … , N K,2 ] are the basis functions. This substitution satisfies the following properties: Using the properties (1)-(6), the optimization problem can be reformulated as Except for the square root of q appearing in the inequalities (10) and (13), the above optimization problem can be posed as a linear optimization problem in which the control points a i represent the free variables (note that any inequality of the form |h| ≤ c can be equivalently expressed as the combination of inequalities h ≤ c and −h ≤ c ). To overcome the nonlinearities in jerk constraints, the optimization problem (7) is first solved without these constraints, leading to the LP formulation: where the vector c T is a uniformly distributed weighting vector. The matrix M is a constant matrix, whose terms are obtained via evaluation of the constraint Equations (8), (9), (11) and (12) at a set of evaluation points. In this way a solution q (so-called pseudojerk) is obtained which realizes a larger feedrate than any other feasible solution. In the next step, the following constraints are substituted for the original jerk constraints (10) and (13), respectively.
These constraint equations are linear with respect to the control points a i of the squared feedrate function q. Thus, the approximate solution of the original optimization problem (7) can be formulated as an LP problem: where the matrix M is a constant matrix whose terms are obtained via evaluation of the constraint equations (8), (9), (11), (12), (14) and (15) at a set of evaluation points. Note that by virtue of the pseudojerk q realizing higher feedrate than any feasible solution, any solution of (16) is guaranteed to respect the axial and tangential jerk limits ̄j max and j max .

Linear programming -evaluation points and knot vector construction
The process of construction of the knot vector and evaluation points of the feedrate curve given a sampling of the toolpath r applied in the experiments described in this paper is briefly described below for the sake of completeness and reader's convenience. Given a sampling in the arc-length parameter s the knot vector of the squared-feedrate function q of order d is defined as Next, a sequence of evaluation points needs to be determined such that for every basis function N i,d , there exists at least one point in this sequence lying in its support. A suitable choice is the sequence of Greville points (Greville abscissae) defined as where N is the number of control points of q and d is the order of q (defined as degree of q + 1 ) and The Greville point g i generally lies near the parameter value corresponding to the maximum of the basis function N i,d [34, p. 512]. The main computational advantage of this particular selection of evaluation points is that the matrices M and M (which comprise the evaluations of the respective constraints at points of G) become sparse band matrices with band equal to d, thus increasing stability and efficiency of linear programming optimization methods.
The method described above is not the only method of knot vector construction available. In fact, the topic of knot vector construction continues to draw attention, especially in the context of curve interpolation and approximation and many methods of knot construction have been proposed. In these methods parameter values are typically associated with sampling points and the knot vector is then constructed so as to satisfy the Schoenberg-Whitney condition [35]. The approaches to knot vector construction include sampling point averaging [34], dominant point selection [36,37] application of support vector machines [38], genetic algorithms [39] and others [40,41]. The applications of knot vector construction method are not exclusively limited to curve interpolation or approximation, see e.g. [42,43].
The knot construction method described in this section was selected because the knot vector distribution corresponds to the distribution of sampling points in the arclength parameter. In future research, it would be interesting to focus on the effects of the various knot construction methods in the context of LP feedrate optimization. Further comments on this topic can be found in Section 3.2.

Curve sampling methods
The simplest way to sample the toolpath curve is to sample the curve parameter interval uniformly, i.e.
This sampling technique is denoted as PAR in the following.
Sampling uniform in the curve parameter (PAR) has been considered in literature dealing with feedrate optimization via the LP approach (see e.g. [6,[44][45][46] and [47]). It is of note, however, that in [6] the authors remark that a nonuniform subdivision of the tool path that would take into consideration its local shape could result in a better performance of the algorithm.
An alternative to the uniform parameter sampling is the uniform arc length sampling (denoted as LEN) such that To apply this sampling technique the information about the parameter to arc length relation is required. As this relation cannot be computed analytically in general, it is necessary to find a sufficiently close approximation. To this end, the method described in [5] was used (note that this does not pose an extra requirement, as the arc length evaluation is also used in the formulation of the LP optimization problem).
The LEN sampling was previously applied in feedrate optimization (see [48]) and in the context of LP feedrate optimization specifically (see [49]). In [50] arc length parametrization was also applied, though without specifying a specific sampling method.
Intuitively, one would expect that a sampling technique would produce a set of sampling points with density that is proportional to the curvature of the sampled curve, while maintaining some minimal point density in areas of zero curvature so that the degree of control over the squared feedrate In the following, this adaptive algorithm is referred to as ADA.
For a comparison of the outputs of the three sampling variants see Figs. 1, 2, 3 and 4.

Sampling techniques -comparison of results
In this section, the effects of the PAR, LEN and ADA sampling techniques on the total machining time of the LP optimization output are presented and discussed. The algorithm was programmed in Matlab2017a software in combination with C++ code for the LP optimization and NURBS curve evaluation via the MEX interface. The COIN-OR Linear programming solver [52] has been used to solve the LP optimization task, while the C++ openNURBS Ⓡ library [53] has been used to construct the NURBS curves and evaluate their derivatives. All computations were performed on a computer with an Intel Ⓡ Core TM i7-7700K processor and Windows 10 operating system. Please note that all the results presented in this section have been obtained from an offline simulation.
In order to compare the optimization results several testing curves have been used. These include the Trident curve, the Butterfly curve, the Pentacle curve and the Phobos curve (see Appendix A, Figs. 16, 17, 18, 19 and Tables 2, 3  of curvature and highest rates of its variation. The Phobos curve represents a spline smoothed contour curve of a blade cross section designed in a commercial CAD system. While all of the chosen testing curves are 2D curves, this does not mean that the algorithms discussed throughout the article are in any way limited to 2D curves only. As was mentioned above, the testing curves were chosen mainly on the basis of being previously used in associated research as well as for being relatively simple to visualize. For each curve, the LP feedrate optimization has been performed with the kinematic limits configuration presented in Table 1. This configuration corresponds to one used on a AXA VCC 1200 machining center equipped with a MEFI CNC872 iTQ-E numerical control system (developed in part by the authors). Each curve was then sampled with increasingly larger values of sampling density (defined as number of points per mm of arc length). For every such sampling density the PAR, LEN and ADA sampling methods were used, the respective machining times were recorded and the relative percentage differences (rounded to the nearest percentage point) of the PAR and LEN methods with respect to the ADA method were calculated. The parameters of the ADA method were chosen as s = = 0.5 and The results presented in Figs. 5-8 support the following conclusions: For curves with low curvature variation such as the Trident, Pentacle and Phobos curves, the ADA sampling technique leads to machining times that are comparable ( ±5% ) with the PAR and LEN methods, while typically being a few percent faster. On the other hand, for curves with high curvature variation, such as the Butterfly curve, the ADA sampling technique leads to machining times that are faster (up to 12% ) than the LEN and ADA methods. This behavior is due to the higher density of knot points of the squared feedrate function in areas of high curvature. Thus, the feedrate profile can reflect the local changes of toolpath geometry more closely in these regions, leading to shorter machining times (while still respecting the kinematic limits). Interpolated trajectory, with details of errors in points of high curvature, can be found in Appendix B (Figs. 20,21,22,23). As can be seen from the presented graphs, the contour errors are quite miniscule.   In conclusion, the ADA sampling technique in combination with the LP optimization approach described in Section 2.1 can lead to significantly shorter machining times when interpolating spline toolpaths with both high curvature variation and high maximal curvature. Such toolpaths are typically encountered in practice in machining of injection molds and in side milling (specifically in trimming operations). In case of side milling, however, jerk optimization is a complex topic as jerk needs to be considered along the tool's contact curve with the machined surface (for additional commentary, see Section 3.2)

Discussion and limitations
In the described implementation, the knot vector is defined via one of several methods of sampling of the toolpath curve and then considered as a fixed input during the feedrate optimization process. The locations of knot points could, however, be also considered as optimization variables. The optimization of knot distribution has primarily been dealt with in curve interpolation/approximation research, where the methods applied include particle swarm optimization [54], gaussian-mixture models [55] and multi-objective genetic algorithms [56] (see also commentaries in Section 2.2 for additional references on knot construction techniques). Besides approximation/interpolation, the issue of optimal knot distribution has also been recently studied in the context of isogeometric analysis [42,43]. It seems that not much attention has been given to the topic of knot vector optimization in feedrate optimization related literature so far, however, in [21] the problematic of optimal knot placement was considered in the context of generation of jerk-minimizing toolpaths via reparametrization (see also [20] for additional details). Several possible methods of knot optimization could be considered in future research such as global optimization methods (genetic algorithms, simulated annealing, etc.) or local optimization such as gradient descent-based sequential quadratic programming.
In any case, the application of the above-mentioned optimization methods would incur substantial computational costs unless a method of determination of a suitable initial knot vector distribution which would produce near optimal   results without the necessity of a large number of additional iterations is devised (as in [21]). In future research, it would be interesting to compare the effectiveness of the presented sampling methods in producing a good initial estimate of the knot vector distribution for subsequent optimization.
As has been noted at the end of the previous section, jerk optimization in side milling is a challenging problem. First of all, a 5-axis variation of the presented LP optimization method would have to be applied (such as a modified version of the algorithm presented in [6]). Secondly, it would be necessary to consider the kinematic constraints along the contact curve of the tool with the machined surface. This could theoretically be solved by considering a discrete sampling of the tool axis and incorporating the kinematic constraints at these points into the LP formulation. It would also be interesting to consider the incorporation of the objective function that minimizes the weighted sum of jerks of the respective boundary curves [21] into the LP formulation. This would require some sort of a linear approximation of the jerk integrals, possibly by applying Gaussian quadrature rules (see, e.g. [57]).

Feedrate override
One of the standard commands in machining is the feedrate override command. With this command, the machine tool operator can change the value of commanded feedrate F cmd to a different value expressed as a percentage of the originally programmed value.
The aim of this section is to present techniques by which the squared feedrate function that was obtained using methods described in Section 2 can be appropriately modified in response to an issued feedrate override command.
This chapter is divided into three section: notation is summarized in Section 4.1, method of override acceleration is presented in Section 4.2 and, finally, method of override deceleration is presented in Section 4.3.

Notation
Most of the notation in the following sections follows that introduced in Section 2. Parameters associated with pseudojerk are distinguished by the use of the hat symbol. The knot vector of the squared feedrate curve is denoted by q and the evaluation point vector is denoted by g. The number of control points is denoted by N, while the number of evaluation points is denoted by K (note that K = (N − d) , see (18)). For the sake of brevity the feedrate and acceleration at the i-th evaluation point are denoted by F(i) and a(i), respectively. The currently active commanded feedrate is denoted by F cur , while the commanded feedrate imposed by feedrate override is denoted as F high in case of override acceleration and F low in case of override deceleration.
Subscripts are used to distinguish parameters pertaining to individual feedrate profiles. Specifically, parameters of the currently active feedrate profile are denoted by the subscript cur and parameters related to the feedrate profile realizing the override command are denoted by the subscript ovr. Parameters related to feedrate profiles globally limited by the override command feedrate in cases of override acceleration and deceleration are distinguished by subscripts high and low, respectively.
Specific control point indices defined in the subsequent section are denoted by I xyz , while evaluation point indices are denoted by J xyz . Finally, the phrase "up to" means "up to but not including" in the following text.

Override acceleration
The process of override transition to a new feedrate limit F high = ovr ⋅ F cur > F cur is performed in several steps: first, a feedrate profile with commanded velocity given by F high is found; second, a feedrate profile that realizes the override command is computed using the information from both the currently used profile and the profile limited by Fig. 11 Pentacle curve -feedrate function comparison at 1 p/mm sampling density Fig. 12 Phobos curve -feedrate function comparison at 1 p/mm sampling density F high . The process is described in detail in the rest of this section.
Let s ∈ [0, L] be the displacement where feedrate override command is received and let n ∈ {1, … , N − 1} such that s n ≤ s < s n+1 , where s i ∈ U are knots of the feedrate curve (17). Denote by d the order of the feedrate curve (in all discussed cases d = 4 is assumed). To construct the override profile, it is first necessary to establish the range of control points of the feedrate curve which need to be fixed in order to guarantee that the override profile coincides with the current profile on some neighborhood of s. This range can be identified as The range of control points given by (19) defines the maximal range of basis functions of the feedrate profile whose supports intersect the interval [s n , s n+1 ) and which therefore influence the values of the feedrate profile on this neighborhood of s. In the case where I l ≥ (N − 1) , the computation of the override profile is skipped, since the override command was received at the very end of the tool path curve. If I l < (N − 1) , it is further necessary to check whether the rest of the feedrate profile reaches current commanded feedrate F cur and whether the override command was not issued during, or immediately before, the final deceleration phase of the current feedrate profile. To check this, consider the point g J s where The point g J s is the first evaluation point after the point of override for which the corresponding value of the feedrate profile is unaffected by the choice of {a I f , … , a I l } . If the feedrate at the end of the curve is fixed, then if the rest of the feedrate profile does not reach the current commanded feedrate, i.e. the computation of the override profile can be skipped. Similarly, if the override command has been issued during (or immediately before) the final deceleration phase, i.e. the computation of the override profile is skipped. This follows from the time-optimality of the current feedrate profile (the deceleration starts as late as possible considering the acceleration and jerk limits). Obviously, the previous two arguments do not apply if the end feedrate is allowed to rise above F cur .
As a second step, the feedrate profile with maximal feedrate given by F high is constructed. The override profile (19) (denoted by the subscript ovr) is then constructed in two steps. Firstly, the pseudojerk q ovr of the override profile is constructed by solving the LP optimization problem: where and The feedrate limitation of the right-hand side b ovr is given by the feedrate values of the current feedrate profile up to g J s and by F high from g J s onward.
The override profile q ovr is then found by solving the LP problem: where the velocity-acceleration part of b ovr is equal to b ovr , while the jerk part is constructed using the values of q cur up to g J s and the values of q ovr from g J s onward.
For an example of override acceleration, see Fig. 13.

Override deceleration
The process of override transition to a new feedrate limit F low = ovr ⋅ F cur < F cur requires a more involved approach than the case of override transition to a higher feedrate limit. The main idea is the following: first, as in Section 4.2, the first evaluation point g J s at which the override profile is allowed to deviate from the original profile is found. Then, for each following evaluation point g next , the arc length necessary to transition from the current profile's feedrate at g J s to the feedrate value of the F low -commanded feedrate profile at g next is estimated. The first evaluation point for which this arc length estimate is not higher than the actual arc length between g J s and g next is then used in the LP formulation of the override profile to define the range of control points to (22) maximize c Tâ ovr subject to:Mâ ovr ≤b ovrl b ovr ≤â ovr ≤ûb ovr , maximize c T a ovr subject to: Ma ovr ≤ b ovr lb ovr ≤ a ovr ≤ûb ovr be bounded from above by the respective values of the F low -commanded profile's control points. The rest of this section explains the algorithm in detail. As in Section 4.2, the first order of business is to define the range of control points of the override profile which need to be fixed (19) and the index J s of the first evaluation point at which the current profile and the override profile are allowed to deviate (20).
Secondly, as in Chapter 4.2, if the override command was issued during or immediately before the final deceleration phase (21), the override profile's construction should be skipped. This follows from the time optimality of the current feedrate profile (the deceleration is as fast as possible considering the acceleration and jerk limits).
If the override command was issued before the final deceleration phase, a feedrate profile q low with commanded feedrate given by F low is constructed. The final velocity of q low is defined as Next, the feedrate and tangential acceleration values of q cur and q low are evaluated at the point g J s . These values are denoted as [F cur (J s ), a cur (J s )] and [F low (J s ), a low (J s )] , respectively. The exact form of the override algorithm then depends on whether F cur (J s ) ≤ F low (J s ) (Section 4.3.1), or F cur (J s ) > F low (J s ) (Section 4.3.2).

Case I: F cur
In this case, it is not necessary to compute any kind of estimate of the arc length necessary to transition from the current profile's feedrate to the feedrate of q low . Instead, it is sufficient to take the appropriate range of control points of F end low = min F end cur , F low .
the q low profile as the upper boundary in construction of the override profile. To this end, define an index I t as If I t ≥ (N − 1) , the rest of the trajectory is not long enough to decelerate below the commanded feedrate. Otherwise, q ovr is constructed by solving the LP problem (22), where and The feedrate limitation of the right-hand side b ovr is given by the feedrate values of the current feedrate profile up to g J s and by F new from g J s onward. The override profile q ovr is then found by solving the LP problem: where the velocity-acceleration part of b ovr is equal to b ovr , while the jerk part is constructed using the values of q ovr .
For an example of case I override, see Fig. 14. To formulate the upper and lower boundaries of the LP solution of the override profile, an upper estimate of the arc length necessary to realize the transition from F cur (J s ) to F low (k) (for some as yet unknown index k) is required.
First, define an index J l as (i.e. the index of the first evaluation point after g J s for which q cur is bounded from above by q low . If no such point exists, the index is set as the index of the last evaluation point).
Next, to estimate the minimal arc length necessary for a transition from q cur to q low , the following algorithm is used: (for the definition of the ArcLen function, see Appendix 3, Algorithm 3 and Algorithm 4). Algorithm 1 searches for the closest evaluation point after which the feedrate transition can already be finalized. The search stops at the index J l . As a next step consider the index J e defined as If J e = K , the remaining arc length of the path curve is not sufficient to realize the feedrate transition and the override command should be skipped (or postponed to the start of the following path segment). Otherwise, a control point index I t is defined as Index I t denotes the first control point for which the support of the corresponding basis function lies past the evaluation point g(J e ) . If I t is undefined or I t ≥ K , the remaining arc length of the path curve is not sufficient to realize the feedrate transition and the override command should be skipped (or postponed to the start of the following path segment). Otherwise, the q ovr profile is constructed using the following LP-optimization problem: where the lower and upper boundaries of â ovr are defined as and respectively. The feedrate solution q ovr is then found by solving the LP problem: where b ovr is constructed from b ovr and q ovr in a standard fashion.
For an example of case II override, see Fig. 15.

Case II: Convergence failure
The LP formulations associated with the construction of q ovr or q ovr described in the previous section can occasionally fail to converge to a solution. This is due to the fact that the function ArcLen assumes a so-called "bang-bang"transition 1 feedrate profile, which exhibits maximal absolute values of either acceleration or jerk. When applied to general spline toolpaths this behavior is not always possible, due to limitations posed by local geometry. This convergence failure thus sometimes occurs when the override command was issued during or immediately before a dip in feedrate caused by a significant change in curvature. In such cases the override profile can still be obtained by setting a new point of override as the nearest local extreme of the feedrate curve and repeating the override computation. 1 Specifically, supposing that either of the LP problems failed to converge, the nearest local extreme point of q cur defined as an evaluation point g J ext can be found via Algorithm 3. The index J ext is well defined since the override formulation is skipped whenever the override command is issued during the final deceleration phase.
Note that in the case of repeated override calculation, the feedrate profile q low and its values at evaluation points have already been obtained and do not need to be recalculated.

Conclusions
This paper comprises two original results related to feedrate optimization via linear programming techniques: the effects of curve sampling methods on time optimality of the resulting feedrate profile (Section 3) and a technique by which a feedrate override profile can be implemented (Section 4). This paper demonstrates that the effect of toolpath sampling can be substantial (when combined with linear programming feedrate planning) and suggests a suitable sampling method. This method can be integrated into the feedrate planning process with low computational overhead, thus making it an interesting choice for practical applications. Note that while the adaptive method requires the approximation of the relation between the curve's natural parameter and its arc length as a prerequisite, the existence of this approximation is already required for the formulation of the feedrate LP optimization problem itself. Thus the use of the adaptive method does not substantially increase the computational complexity of the presented feedrate LP optimization technique. In addition, both the approximation and the sampling itself can (should) be performed in the preprocessing stage of feedrate planning and therefore do not add to the computational complexity of the online stage. Given the scope of this paper, sampling methods were compared using only plain axis position data generated by the algorithm. Further comparison could be performed on a CNC machine in the future with the application of a combined speed and error comparison method [59].
In order to successfully apply linear programming techniques to real world feedrate planning, a method capable of recalculating feedrate profile in reaction to a feedrate override command is essential. This topic, however, has not been dealt with in the past. Thus, a possible implementation of such a method is presented in this paper. While this method does require the computation of two additional feedrate profiles (one bounded by the newly issued commanded feedrate and second that realizes the desired feedrate transition), the majority of prerequisites for the computation of these profiles are shared with the original profile. Specifically: the knot vector associated with the feedrate profile, its evaluation points and the constant matrices of both the pseudojerk and feedrate curve optimization problems can be stored in memory the first time they are computed and reused in subsequent optimizations. Main part of the computational complexity is therefore due to the use of linear programming solvers as all the additional computations required are either index manipulations, or simple analytic formulas. The method is therefore suitable for application purposes due to its low memory requirements and simplicity of auxiliary calculations, as long as a suitably fast linear programming solver is available.

Appendix C
This section contains the definition of the ArcLen function used in Section 4.3.2, Algorithm 1. The ArcLen function calculates an estimate of the arc length necessary to decelerate from starting feedrate F s to end feedrate F e with zero acceleration and limits on tangential acceleration and jerk given by a max and j max , respectively. The arc length estimate is obtained Fig. 23 Butterfly test curve -interpolated trajectory comparison Fig. 22 Phobos test curve -interpolated trajectory comparison via a so-called "bang-bang" feedrate profile, i.e. a profile that maximizes the absolute value of either jerk or acceleration at every point. In case the maximization of acceleration would lead to reaching the end feedrate F e with nonzero acceleration, a new maximal deceleration value a mid is calculated, so that the end feedrate can be reached with exactly zero acceleration with a mid replacing a max as a new acceleration limit (see code blocks starting at lines A and B in Algorithms 3 and 4 below).