Novel Feedrate Optimization Method for NURBS Tool Paths Under Various Constraints

A good feedrate profile for NURBS tool paths should be able to consider both machining efficiency and motion limits. In other words, the feedrate profile can lead to the machining time as small as possible and can satisfy various path and axis motion constraints. However, owing to the complicated relationship between trajectory motion and axes, it is difficult and time-consuming to check axis motion constraints (e.g., axis velocity, axis acceleration, axis jerk). In addition, the path length is not easy to calculate, so the machining time cannot be accurately estimated. This paper presents a feedrate profile planning method that can easily include, delete or organize the constraints and minimize the machining time for NURBS tool paths. The key idea is to use the same parameter <inline-formula> <tex-math notation="LaTeX">$u$ </tex-math></inline-formula> to represent both NURBS tool path and feedrate profile. When designing the feedrate profile with <inline-formula> <tex-math notation="LaTeX">$u$ </tex-math></inline-formula>, an analytical form of machining time can be obtained, and the constraints are converted into functions of <inline-formula> <tex-math notation="LaTeX">$u$ </tex-math></inline-formula> (including feedrate, tangential acceleration, tangential jerk, axis velocity, axis acceleration, axis jerk, chord error, and centripetal acceleration), which can be quickly checked. Particle swam optimization algorithm is employed to eliminate the solutions with long machining time or do not meet the constraints, and gradually optimize the feedrate profile. In the simulation, a star-shaped NURBS curve is selected to demonstrate the effectiveness of the proposed method. The results show that the proposed method can achieve not only smooth and high-speed machining under various constraints, but also high accuracy with the minimization of position error at the final interpolation point.


I. INTRODUCTION
Owing to its smooth properties and capability of generating free-form surfaces, the non-uniform rational B-splines (NURBS) has proven to be superior over linear and circular paths and become the standard format for computeraided design (CAD) systems. To achieve a high-speed and high-precision computer numerical control (CNC) machining, the NURBS interpolator has been extensively studied, and many new interpolation algorithms have been proposed recently [1], [2]. There are two issues to realize a NURBS interpolator. One is to plan the feedrate profile, and the other is to generate accurate interpolation positions for each sampling time. So far, many efficient algorithms have been proposed for acquiring accurate interpolation positions The associate editor coordinating the review of this manuscript and approving it for publication was Zhenzhou Tang .
according to a specified feedrate [3]- [5] or based on finite impulse response (FIR) filtering [6], [7]. Especially the FIR filtering-based algorithm can realize the real-time interpolation. The main concern of the present work is the planning of the feedrate profile, which is an important topic since it directly affects the machining efficiency and accuracy.
In general, it is preferred to have higher feedrate for machining efficiency. However, higher speed may easily lead to larger dynamic errors. In addition, there exist many motion limits on velocity, acceleration, and jerk due to physical constraints. A good feedrate profile should be able to consider both machining efficiency and motion limits. In other words, the feedrate profile can lead to the machining time as small as possible and can satisfy various path and axis motion constraints. However, owing to the complicated relationship between trajectory motion and axes, it is not easy to plan the speed under various constraints.
There have been many literatures on the feedrate scheduling of NURBS paths under various constraints. Gao et al. [8] applied a jerk-varied flexible feedrate scheduling method to slow down the acceleration near the points with local maximum curvature. Jia et al. [9] proposed a feedrate scheduling method with constant speed in the feedrate sensitive regions, which avoids frequent acceleration and deceleration in large curvature areas and guarantees geometric or drive constraints. Zhong et al. [10] developed a parametric curve interpolator that satisfies the chord error, feedrate limit, and axis velocity and acceleration limits. Ni et al. [11] considered the constraints of chord error, centripetal and tangential acceleration, and jerk. In addition, some feedrate profiles composed of S-curve [12] or bell-curve [13], [14] are used to meet acceleration and jerk limits along the path.
There are several problems in the above feedrate scheduling methods. One is the difficulty in changing the considered motion constraints. For example, constraints on feedrate, acceleration, and jerk along the path are considered in [12], but axis motion constraints will be difficult to be incorporated into the method proposed in [12]. Another problem is that to check the constraints (e.g., chord error, axis velocity, axis acceleration, axis jerk), all the interpolation positions along the tool path or on each axis need to be obtained. However, in the practical application, it is time-consuming and not feasible. Wang et al. [15] scanned the tool path to check the points with large curvature, and planned the speed of these points only. However, determining the exact position where the speed should be increased or decreased is still an open problem. Usually, integrals are used to estimate the distance for acceleration and deceleration (acc/dec). However, the numerical integration error may cause feedrate fluctuation [11], especially when the feedrate is high.
In addition to motion constraints, time optimization is another important issue for feedrate planning. For instance, Yuan et al. [16] developed a back and forward check algorithm under an acceleration constraint, which has higher feedrate or shorter operation time. Zhang et al. [17] proposed a greedy algorithm that enables at least one axis to reach its acceleration or jerk limit during machining. Afterwards, some other extended investigations with the linear programming algorithm were presented by Liu et al. [18] and Erkorkmaz et al. [19]. The key issue in time optimization is the difficult estimation of the machining time. There is a complicated relationship between the machining time, the feedrate profile, and path length. Hence, most existing studies did not plan feedrate from a time perspective. That is, the time is not optimized directly; instead, the feedrate is planned as large as possible by considering path geometry or motion constraints, such as curvature, acceleration limits, etc. However, the path length needs to be computed and also complexity increases dramatically when more motion constraints are considered.
To design the feedrate from a time viewpoint, this study presents a novel planning method which designs the feedrate profile in the parameter domain. In general, the feedrate profile is represented in the time domain and specific type of acceleration/deceleration is adopted ( Fig. 1(b)), e.g., trapezoidal, S-curve, etc. It is inconsistent to the NURBS path that is represented in the parameter domain ( Fig. 1(a)). Such inconsistency makes the feedrate panning difficult and require complicated computation in relating the parameter and time. In this paper, the feedrate profile will be represented by the same parameter u used to define the NURBS tool path ( Fig. 1(c)). As a result, the machining time can be easily obtained without the need of computing the path length and can be optimized directly. In the optimization, particle swarm optimization (PSO) [20]- [22] is employed to minimize the machining time, and fitness selection is used to deal with various constraints [22]. When a feedrate profile does not meet the constraints, a smaller fitness value is attached. Since both the NURBS path and feedrate profile are polynomial or rational polynomial functions in u, the check of motion constraints can be efficiently done by utilizing the property of polynomial function. The proposed method has the following advantages: 1) There exists a mapping from tool path to the designed feedrate and thus there is no need to calculate the path length which often causes feedrate fluctuation. 2) The machining time can be accurately estimated in advance. 3) Most constraints can be expressed as the functions of u and checked by evaluating the values from u = 0 to 1. It is more efficient than checking the fine interpolation points at every sampling period. Furthermore, using PSO algorithms, these constraints can be easily included, deleted, or organized. The rest of this paper is organized as follows. Section 2 reviews some nature-inspired algorithms and two common methods to generate NURBS fine interpolation points. Section 3 describes the proposed feedrate planning strategy. Section 4 shows the simulation results of the proposed method applied to a star-shaped NURBS curve. Finally, section 5 summarizes this paper. VOLUME 10, 2022

II. RELATED WORK
In this section, we briefly review some nature-inspired algorithms (Section II.A) and two common methods to generate NURBS fine interpolation points (Section II.B).

A. NATURE-INSPIRED ALGORITHMS
In recent years, various optimization problems were discussed [23], [24]. These optimization problems consist of different type of objective function and difficulty levels such as linear and/or nonlinear constraints, which can be equality and/or inequality, and variety of search space such as discrete and mixed type (continuous, discrete). To handle these problems, a current trend is to use natureinspired algorithms, which work independently of whether the function involved in the problem is differentiable, continuous, or convex, etc. They just evaluate the objective function at given decision variables and consider the optimization problem as a black box. Such algorithms mainly include genetic algorithm (GA) [25], [26], PSO [27]- [32], fruit fly optimization algorithm (FOA) [33], queuing search algorithm (QSA) [34], atom search optimization (ASO) [35], equilibrium optimizer (EO) [36]- [38], evolutionary programming (EP) [39], differential evolution (DE) [40]- [42], etc. Many researchers extended these algorithms to solve the multi-objective optimization problems [43]- [45] or integrated them into multi-population methods to improve the optimization performance [46]- [48]. Among these algorithms, PSO is a considerably popular one which has some attractive features, such as simplicity, less parameters, low computational complexity, and ease to implement. Over the years, several PSO variants have been developed to solve complex industrial and engineering application problems. For example, Liu et al. [49] integrated PSO search strategy and Lévy flight to optimize the deployment of access points for wireless local area networks. Yiyang et al. [50] proposed an inverse kinematics calculation method based on improved PSO algorithm and applicable to general robots. Rahman et al. [51] presented a hybrid GA and PSO algorithm to solve real-time order acceptance and scheduling problems in a flow shop environment. To our best knowledge, no one has used nature-inspired algorithms to solve the NURBS tool path federate planning problems under constraints based on minimum machining time. Since PSO is relatively matured, commercially available, and can be easily accessed by industrial people, it is employed in this study.

B. FINE INTERPOLATION POINTS COMPUTATION
NURBS has been popularly used because it can represent both analytic and free-form surface. A NURBS curve can be expressed as follows: where P i is the control point, w i is the weight of P i , n + 1 is the number of the control points, p is the degree of NURBS, N i,p (u) is the B-spline basis function, and u is the interpolation parameter. The pth-degree B-spline basis function can be calculated using the recursive formulas given as follows: where U = {u 0 , u 1 , . . . , u n+p+1 } is the knot vector.
To generate the fine interpolation points, the interpolation parameter u for every sampling time needs to be determined. Taylor approximation is a commonly used method. By employing Taylor series expansion, u at (k +1)th sampling time can be expressed as [52] where T s is the sampling time of the NURBS interpolation, H (t k ) stands for the high order term in the Taylor expansion, v(t k ) is determined according to the specified velocity profile, for example trapezoid or S-curve, and the mth derivatives of a NURBS curve and the B-spline basis function with respect to u are given respectively as [53] To improve the accuracy, higher order Taylor expansion has been applied, but there is always a trade-off between the computing time and the desired accuracy. Another popular method for determining the new interpolation parameter u(t k+1 ) is to compute the incremental interpolation parameter u(t k+1 ) based on the previous interpolation parameter u(t k ) and the incremental path where The new interpolation parameter u(t k+1 ) can then be obtained as follows: After u(t k+1 ) is obtained from (4) or (8), the next interpolation position C (u(t k+1 )) can be determined as well.

III. QUINTIC PROPOSED FEEDRATE PLANNING
A NURBS curve is a curve represented by a vector function in a parameter u denoted as For the machining of a NURBS curve with a CNC machine tool, the command as a time function for each motion axis is required. To this aim, one needs to plan the feedrate profile along the path and then generate the fine interpolation points [3]- [5], [8], [9], [54]. This work presents a novel method to plan the feedrate profile for NURBS tool paths.
The proposed method will optimize the machining time and generate the smooth commands that satisfy all of the prescribed motion constraints, such as axis acceleration and jerk limits, etc. The conventional approach is to perform the planning in the time domain, such as trapezoidal or S-curve feedrate profile [2], [9]- [12]. The problems in such approach include: (i) the total machining time cannot be accurately estimated; (ii) there exists feedrate fluctuation.
In this study, the feedrate planning will be performed in the parameter domain. In other words, both the path function and the feedrate profile will be represented in the same parameter u. The advantage of this approach is that u is always from 0 to 1 and the total machining time can be accurately estimated. To see this, letv (u) be the feedrate profile and s (u) be the arc length function, i.e., where the differentiation convention of Integrating on both sides of Eq. (10), the machining timeT t can be expressed as: Note that (11) can be accurately calculated within any specified error range. Being able to obtain the accurate machining time is important since it allows us to implement the optimization algorithm for minimizing the machining time. Given a desired motion path described by NURBS, i.e., C (u), it is to design a feedrate profilev (u) to minimize the machining timeT t and at the same time satisfy all of the motion constraints. The PSO algorithm is employed. Fig. 2 is a flowchart of the proposed PSO-based feedrate optimization algorithm, which is detailed below.

A. FEEDRATE PROFILE GENERATION
Most conventional feedrate profiles adopt simple functions (e.g., trapezoidal or S-curve, etc.) in the time domain, which are restricted and difficult to be integrated with NURBS paths. In this study, a more flexible feedrate profile will be taken so that better optimization result can be expected.  (12) where {a i , b i , c i , d i } are constant coefficients to be determined by the smoothness conditions of the profile at FCPs, as shown in Appendix A. For the first and last pieces (i = 1 and M − 1), the profile cannot be designed as a cubic spline (12) since ill condition will occur at the starting and ending points (u = 0 and u = 1), where the federate must be zero. Note that from (10), the tangential accelerationâ (u) can be obtained asâ It is easy to see thatv (u) andâ (u) (and higher order derivatives like jerkĴ (u), etc.) all vanish at u = 0 and u = 1, butv (u) is finite if (12) is assumed. This implies that the path can hardly start at u = 0 and moves slowly at almost zero speed near u = 1. To resolve such situation, a nonzero acceleration at u = 0 and u = 1 is required. From (13), this can be achieved by settinĝ with the condition thatv 1 (0) = 0;v M −1 (1) = 0. Here, p and q are coefficients and again they are to be determined by the smoothness conditions of the profile at FCPs, as discussed in Appendix A. Equation (14) leads tô   axis and F max is the feedrate limit. Note that V 1 = V M = 0. With V i , the feedrate profile given by (16) can be generated.
Next, one needs to check if the generated feedrate profile satisfies all of the motion constraints, such as acceleration limit, jerk limit, etc. To this aim, most conventional methods require to compute all of the interpolation positions along the tool path or on each axis [55]- [57]. In this study, a more efficient way is proposed by representing the motion functions in terms of u. The motion functions considered include tangential velocityv (u), accelerationâ(u), and jerkĴ(u), and axis velocity v(u), acceleration a(u), and jerk J(u), as shown in Table 1

, wherê a(u),Ĵ(u), v(u), a(u), and J(u) are derived in Appendix B.
Note that all of the motion functions are rational polynomial or square root of rational polynomial functions of u. This is because that both the NURBS path function C (u) and the feerate profilev (u) are these types of functions. As a result, all of the motion constraints can be represented as polynomial inequalities. For example, let us consider the limit on the  tangential acceleration: From Eq. (13), it is obvious thatâ 2 (u) is rational polynomial in u. Suppose thatâ where N a (u) and D a (u) are polynomial functions. Thus, inequality (17) can be rewritten as The other motion constraints can be treated similarly. Therefore, there is no need to generate all of the interpolation positions for checking the constraints. All of the motion constraints can easily be checked by solving the polynomial inequality or calculating their values from u = 0 to 1.

C. OPTIMIZATION PROCESS
When one or more constraints are not satisfied, the fitness value is set to zero; otherwise the fitness is set to the inverse of the machining time, i.e., 1/T t , which can be accurately computed by Eq. (11). Thus, an optimal set of V i , i = 2, 3, . . . , M − 1, in FCPs can be obtained for the current iteration. Equivalently, an optimal feedrate profilê v (u) can be obtained. Next, the stop criterion needs to be checked (in this study, a maximum number of iterations is used). If the optimal feedrate profile does not meet the stop criterion, a new population sets of V i will be generated and the PSO optimization process continues. The PSO optimization process is stopped if the stop criterion is met. The PSO optimization algorithm has been employed here. In fact, any other fitness-based optimization algorithm (e.g., GA, QSA, ASO, EO, etc.) can be applied to optimizev (u) (more specifically, V i in FCPs). VOLUME 10, 2022

D. FINAL POSITION ERROR IMPROVEMENT
Since the machining timeT t obtained by (11) may not be an integer multiple of the interpolation period (sampling time) T s , there will be a residual length for the final interpolation period [53]. The problem of residual length can be eliminated by slightly reducing the optimal velocity profile, which is detailed below. Let N = T t T s , where x is the ceiling function denoting the minimum integer greater than x. Then, the optimalv (u) can be modified as:v By Eq. (11), the machining time of the feedrate profilev m (u) will be NT s , an integer multiple of T s . Sincev m (u) is less than (or equal to)v (u) for every parameter u, it is guaranteed that all of the motion constraints satisfied byv (u) will also be satisfied byv m (u). The modification given by (20) is also one of the advantages of planning the feedrate in the parameter domain. If the feedrate had been planned as a time function, it is difficult to adjust the machining time and keep the path length (the integration of the feedrate up to the machining time) unchanged.

IV. SIMULATIONS
To verify the proposed method, a star-shaped NURBS curve as shown in Fig. 4(a) is used as an example. This curve has four high-curvature regions, as marked with c1, c2, c3, and c4, respectively, in Fig. 4(b). The high-curvature regions are critical for machining and thus the feedrate should be slowed down in these zones. So, this curve is helpful to validate the effectiveness of feedrate planning. In this example, 30 FCPs (i.e., V 1 , . . . , V 30 ) are created uniformly on u ∈ [0, 1] to construct a feedrate profile, where V 1 and V 30 are set to zero. PSO is employed to optimize V 2 , . . . , V 29 . The population size is 15, and the maximum number of iterations is 500. Other relevant parameters are listed in Table 2.   In the simulation, various constraints are analyzed. First, only the axis velocity limits are considered. In this case, we are going to design a high feedrate, but each axis speed does not exceed 20 mm/s. Figure 5(a) shows the result of the optimization of the feedrate profile, which produces a machining time of 10.50843 seconds. Clearly, it is not an integer multiple of interpolation period (i.e., 0.5 ms). This means that there will be a little time left when executing the last interpolation. To improve this, the feedrate in Fig. 5(a) is multiplied by a scale factor ''10.50843/10.5085'', where 10.5085 is a modified machining time that is slightly larger than 10.50843 but integer multiple of 0.5. Figure 5   value when approaching these high-curvature zones and then increase gradually to a permissible value. Figure 6(a) shows a new profile of the feedrate. It can be found that there exist four slow-down zones in the feedrate profile corresponding to the four high-curvature zones on the starshaped path curve. The total machining time is obtained as 11.3550, which is slightly higher than that for the case that only considered the axis velocity constraint. Nevertheless, Figs. 6(b) and 6(c) show that both axis velocities and axis accelerations are constrained within the given limits of [−20 mm/s, 20 mm/s] and [−50 mm/s, 50 mm/s], respectively.
Next, consider the different constraints listed in Table 3.  are considered as constraints for smoother motion. The results are shown in Fig. 8. As can be seen, the proposed method VOLUME 10, 2022   achieves a maximum velocity along the x-axis of 20 mm/s and/or the y-axis of 15 mm/s and a maximum acceleration along the x-axis of 25 mm/s 2 and/or the y-axis of 20 mm/s 2 without exceeding the jerk constraints (x-axis: 65 mm/s 3 ; yaxis: 70 mm/s 2 ). The corresponding total machining times for these cases are shown in Table 3.
For PSO, increasing the population size will increase the diversity of the solution, but also increase the calculation time under the same number of iterations. Table 4 shows the optimization results and computational time under population sizes of 20, 30, and 50. As can be seen, PSO can always find fairly consistent solutions (machining time of 15.10837 seconds). It can also be seen from Fig. 9 that convergence to a stable solution is achieved after 200 iterations for all three choices. In addition, we also observed the convergence behavior with respect to the choice of (c 1 , c 2 ) and w [58]. We found that these parameters have little effect on final results and have quick convergence (Figs. 10 and 11), thus saving a lot of parameter design work.
For the star-shaped NURBS curve path, the number of FCPs is determined by trial and error through simulations. Although in theory, using more FCPs helps to optimize the feedrate profile, it will inevitably take more optimization time. In future work, we will propose a systematic design strategy for the number and location of FCPs.

V. CONCLUSION
This study proposed a feederate optimization method for NURBS tool paths under various constraint combinations. The conclusions include the following: 1) The feedrate profile is designed with parameter u instead of commonly used time t, there is no need to calculate the path length which often causes feedrate fluctuation.
2) The machining time can be expressed as an analytical form. This facilitates optimizing the machining time and alleviate the residual length after executing the last interpolation. 3) Most constraints can be expressed as the functions of u (including feedrate, tangential acceleration, tangential jerk, axis velocity, axis acceleration, axis jerk, chord error, and centripetal acceleration) and easily checked. In addition, using the fitness-based optimization algorithms, these constraints can be easily included, deleted, or organized. The proposed optimization method is also possible to be applied to a five-axis machine tools with two additional rotary axes. However, this requires further research.

APPENDIX A
In the proposed feedrate profile, M FCPs will produce M − 1 cubic splines (as shown in Fig. 2). Except for the first and last splines, which have only one polynomial coefficient (i.e., p and q, respectively), all other splines have x (u) 2 + y (u) 2 + z (u) 2 four polynomial coefficients (i.e., a j , b j , c j , and d j ). So, there are a total of 1 + 4(M − 3) + 1 = 4M − 10 unknown coefficients. Accordingly, 4M − 10 equations are needed to determine these coefficients. First, we require that each spline connects adjacent FCPs, that is: where i = 2, . . . , M − 2. Hence, there are 2(M − 3) of these conditions. Then, to makev (u) as smooth as possible (first and second derivatives are continuous) at everyũ i , we require: where s (u), as shown at the top of the page. According tȯ u,ü, and u, the axis velocity v (u), axis acceleration a (u), and axis jerk J (u) can be expressed with u, respectively, as: