Pareto Optimization of Cycle Time and Motion Accuracy in Trajectory Planning for Industrial Feed Drive Systems

Manufacturing industries aim to improve product quantity and quality. These trade-off objectives are typically manifested as cycle time reduction and motion accuracy improvement. Corner smoothing approaches that reduce the cycle time of piecewise linear trajectories are proposed in the literature. This study tackles the two objectives by Pareto-optimal corner smoothing with constraints imposed as kinematic limits, continuity conditions and user-specified cornering tolerance. Linear and cornering motions along a contour are respectively described by jerk-limited acceleration profiles and a modified kinematic corner smoothing with interrupted acceleration (KCSIA) approach. A Pareto frontier is generated by the divide and conquer algorithm, where the solution nearest to the utopia point is selected as the best trade-off solution. The effectiveness of the proposed method is validated through experiments, where the best trade-off solution reduces the maximum and average contouring errors by 47.53% and 25.40% while it increases cycle time by 2.53% compared to KCSIA.


I. INTRODUCTION
Manufacturing industries typically use computer numerical control (CNC) machine tools due to their accuracy, repeatability and speed in performing tasks [1]. Feed drive systems (FDSs) actuate CNC machine tools' motion axes [2]. Ongoing demands for higher production quantity and quality drive researches in motion accuracy improvement and cycle time reduction. Several studies have proposed feedback control structures to reduce tracking [3] and contouring [4], [5] errors in feed drive systems. Chen and Sun propose a nonlinear controller for underactuated systems [6]. A feedback control strategy is proposed in [7] for regulating ship yaw and roll perturbations in 5-DOF offshore cranes. Such proposals are limited by the accessibility of in-service feed drive system feedback controllers.
The associate editor coordinating the review of this manuscript and approving it for publication was Jon Atli Benediktsson . Trajectory generation methods have been proposed in the literature to reduce errors and cycle time. Feedforward compensation strategies such as iterative learning control [8], [9] and neural networks and reinforcement learning [10] are proposed for error reduction. A time-optimal trajectory generation approach with consideration to obstacles and dynamic limits is proposed by Uchiyama et al. for robotic manipulators [11]. Sencer et al. propose time-optimal feed scheduling along B-spline tool paths for 5-axis CNC machine tools [12]. Frequency-optimal acceleration profiles are proposed by Sencer and Tajima for vibration suppression [13]. Kucuk proposes minimum time trajectory generation using cubic spline and 7 th order polynomial interpolations [14]. 15-phase sinusoidal jerk profiles are proposed by Fang et al. for cycle time reduction and high-frequency harmonic suppression [15]. Wang et al. propose time-optimal S-curve velocity profile generation for robotic arms [16]. In economic lot scheduling of supply chains, a sub-optimal cycle time of a process lot can be selected to reduce overall production costs [17], [18]. To this end, Jeong et al. present time-optimal and time-fixed jerk-limited velocity profile generation algorithms [19]. Besset and Béarée propose online finite impulse response-based trajectory generation for time-optimal, fixed time and jerk-time fixed cases [20].
Before a contour is fed to feed drive systems, CAD/CAM systems normally discretize it into a set of linear and circular arc segments using G01 and G02/G03 commands, respectively. In order to improve machining quality and speed, corner smoothing methods have been proposed in the literature. Corner smoothing methods can be classified as global and local corner smoothing methods based on the span of a fitted curve. Global corner smoothing methods fit a single curve spanning across all segments while local corner smoothing approaches fit a curve between a pair of adjacent segments [21]. Global corner smoothing is normally used for motion planning along short-segment paths [22]- [24]. However, it is more challenging to evaluate and constrain the smoothing error for global corner smoothing compared to local corner smoothing [25].
Local corner smoothing can be categorized as geometric and kinematic local corner smoothing, where geometric local corner smoothing separately considers geometric and kinematic constraints in two steps while kinematic local corner smoothing plans smooth velocity transitions directly by considering both constraints in one step [26]- [28]. In literature, geometric local corner smoothing using B-spline [29]- [31], Pythagorean Hodograph [32] and Bézier [27] curves are proposed. Regarding kinematic local corner smoothing, kinematic corner smoothing with interrupted acceleration (KCSIA) has been proposed for cycle time optimality [33]. An energy-time trade-off using KCSIA is studied in [34]. Finite impulse response-based kinematic local corner smoothing has been proposed for reducing time and vibrations [35], [36]. Kinematic local corner smoothing approaches using clothoids [37]- [39] and asymmetrical double constant-jerk cornering profiles [21] have been proposed for improving cycle time. Regarding CNC machining, there are several different cases of requirements for surface quality and cycle time: a high-accuracy case, a high-speed case with a certain accuracy level and a time-fixed case. With the aim of reducing cycle time while satisfying accuracy requirements for piecewise toolpaths, existing local corner smoothing algorithms typically maximize cornering velocities, where cornering errors are driven to the upper bounds of the accuracy constraints [27], [40]. Under sub-optimal cycle time scenarios, driving these errors to their upper limits is no longer necessary. Thus, various trade-off solutions for cycle time and cornering error can be selected by a decision-maker.
A Pareto set of a multi-objective optimization problem (MOOP) can be generated by vectorization or scalarization methods [41], [42]. Vectorization methods are stochastic approaches that directly solve the MOOP to produce globaloptimal solutions. Their computation cost and stochastic nature limit their application. Scalarization methods solve a MOOP by parameterizing it into a series of single-objective optimization problems (SOOPs), resulting in a set of locally optimal solutions. Taking a weighted sum of objectives is a commonly used approach, although it has drawbacks in obtaining solutions in non-convex Pareto regions [43]. Approaches such as normal boundary intersection [44], normalized normal constraint [45], [46] methods have been proposed. The normal boundary intersection approach and, to a lesser degree, the normalized normal constraint are prone to generating non-Pareto optimal solutions. Logist and Van Impe propose a criterion for detecting such solutions [47]. The scalarization methods typically distribute the SOOPs evenly, resulting in a waste of computational effort on insignificant Pareto points. Kim and Weck propose a recursive approach for generating a Pareto frontier based on an adaptive selection of objective weights in a bi-objective problem [48]. This approach requires four user-defined parameters to control the Pareto frontier approximation. Hashem et al. propose the divide and conquer algorithm for the recursive exploration of significant trade-off regions on the Pareto frontier, where one user-defined parameter, named the minimum trade-off level, is required to control the Pareto frontier resolution [42].
In order to account for the different cases of surface quality and cycle time requirements, a Pareto-optimal local corner smoothing method that offers a trade-off between cycle time and motion accuracy is proposed. A decision-maker can systematically select Pareto-optimal solutions to address the requirements under consideration. Piecewise linear contours are considered in this study due to their regular occurrence in CAD/CAM systems. The normalized normal constraint method is used to formulate the bi-objective optimization problem since it is known beforehand whether the Pareto frontier does not have non-convex regions. The motion accuracy improvement objective is indirectly represented by the minimization of corner smoothing extent along an entire contour. As part of the optimization problem formulation, it is considered that a smoothed trajectory's linear and cornering segments are respectively described by jerk-limited acceleration profiles (JLAP) [49] and a modified KCSIA. The optimization problem is constrained by kinematic limits, user-defined cornering tolerances and geometric restriction for avoiding path overlaps. The divide and conquer algorithm [42] is used for generating an approximated Pareto frontier since it only requires one user-defined parameter, named a minimum trade-off level, and it is compatible with the normalized normal constraint method. Each Pareto-optimal solution is computed using sequential quadratic programming [50]. In this study, the Pareto-optimal solution that is nearest to the utopia point is selected as the best trade-off solution. In a time-fixed case, a solution can be obtained without the need to generate the Pareto front, where the motion accuracy improvement objective is minimized subject to a userspecified cycle time and the above-mentioned constraints.
Experimental results of the generated Pareto-optimal trajectories demonstrate the effectiveness the proposed approach. The contributions of this work are summarized as follows: • A cycle time and motion accuracy trade-off by Paretooptimal local corner smoothing is proposed.
• The time-optimal solution has a higher motion accuracy and shorter cycle time than KCSIA. The paper is organized as follows: Section II illustrates the design of Pareto-optimal trajectories, Section III presents the optimization and experimental results. A discussion on the findings and concluding remarks are described in Sections IV and V, respectively.

II. PARETO-OPTIMAL TRAJECTORY DESIGN
A method of generating Pareto-optimal trajectories that provide a compromise between motion accuracy and cycle time for piecewise linear contours is illustrated in this section. The bi-objective optimization problem and constraint formulations are depicted below.

A. PROBLEM FORMULATION
Corner smoothing approaches reduce cycle time by shortening the tool path length and providing non-zero cornering velocities, which consequently deteriorates contouring performance as shown in [22], [27], [40]. Hence, motion accuracy improvement and cycle time reduction are conflicting objectives in trajectory generation for piecewise linear contours.
Corner smoothing methods generate cornering trajectories for given machine tool kinematic limits (ie., jerk, acceleration and velocity limits) and a user-specified cornering tolerance 0 < ε ≤ ε ub as shown in Fig. 1, where p c is the original corner point and the cornering error ε is restricted by upper bound ε ub . The cornering path from position vector p s to p e is symmetrical about line p c p mid , where the distance between the start/end of the path and the original corner p c is the cornering Euclidean distance L c . h s and h e are the direction vectors at the start and end of the path, respectively. The minimization of ε and reduction of L c are non-conflicting objectives. Hence, a SOOP is used to describe the smoothing minimization objective for the m th corner, where µ c,m is parameter vector that consists of variables that describe a cornering trajectory. As detailed in Section II-C2, these variables depend on the corner smoothing method used. g is an equality constraint vector that imposes motion continuity conditions when stitching the corner path with preceding and succeeding linear segments. q is an inequality constraint vector that ensures the cornering motion obeys kinematic limits and the user specified cornering tolerance. The trajectory cycle time is the second objective to be minimized. It is obtained as a sum of the linear and corner segment trajectory durations T l,m s and T c,m s, respectively. µ l,m is a parameter vector that describes the m th linear segment trajectory. n l is the number of linear segments. In this study, it is assumed that the end of a piecewise linear contour is not smoothed, hence n l = n c +1.
The trade-off between corner smoothing minimization and cycle time reduction is represented as the bi-objective optimization problem where L c,tot is the total cornering Euclidean length for corner smoothing. Thus, L c,tot represents the smoothing objective function for all n c corners. µ is the optimization parameter vector consisting of linear and corner segment variables. The g elements describe geometric and C n continuity conditions along a resulting optimal trajectory and q restricts the trajectory within kinematic limits and user-specified cornering tolerances. Each objective extremum is obtained by independently minimizing the components in (4) as: , µ lc = arg min µ L c,tot (µ), 114106 VOLUME 9, 2021 FIGURE 2. An illustration of the reformulated optimization problem (9) with the feasible region reduced by the normalized normal constraint inequality constraint (10). The utopia point is marked by the origin O.
In a normalized objective space, describe the respective individual objectives in (4), where L c,tot = L c,max − L c,min and T cycle = T max − T min are the saving potentials in corner smoothing and cycle time, respectively. In accordance with the normalized normal constraint method [45], (4) is reformulated as: min µT cycle (µ) (9) subject to (5) and an additional normalized normal constraint inequality constraint that limits the feasible region in the normalized objective space ( Fig. 2) with where a point ψ on a Pareto frontier corresponds to the point ρ on the convex hull of individual minima obtained at a weighting factor ζ .
is generated using the divide and conquer recursive structure [42], where (9) is solved by sequence quadratic programming [50] for different ζ . In this work, the Pareto optimal point that is nearest to the utopia point is considered to be the best trade-off solution obtained at the weighting factor ζ * which corresponds to the parameter vector µ * in the decision space. The optimization process described above is demonstrated as Algorithm 1 in the appendix, where it receives a piecewise linear tool path a series of 2D points p c s, cornering tolerance 0 < ε ≤ ε ub and kinematic limits for each axis. Under a time-fixed case, the corresponding Pareto-optimal solution ψ can be directly obtained without the need for generating by solving the single objective problem min µ L c,tot (µ) (13) subject to (5) and a cycle time constraint where T fixed is the fixed cycle time defined by a user.

B. OPTIMIZATION CONSTRAINTS
In order to ensure that a Pareto optimal solution is implementable in a real FDS, kinematic limitations for the k th axis must be incorporated as optimization constraints. Hence, kinematic constraints are defined, where r k is a trajectory position and j lim,k , a lim,k , and v lim,k are the respective jerk, acceleration and velocity limits. (15) is described in a quadratic form so as to guarantee differentiability ∀µ ≥ 0. At each corner, the m th cornering error is bounded by the constraints Geometric constraints (17) are defined in order to avoid overlapping the (m − 1) th and m th cornering paths, where s m is the path length of the linear segment between the smoothed corners. In order to ensure C n smoothness when motion switches from lines to cornering paths and vice versa at time instants t ms , continuity constraints are established, where the motion profiles for linear and cornering paths are respectively r l and r c .

C. TRAJECTORY REPRESENTATION
This section illustrates the motion profiles used for Pareto optimal trajectory design. With no loss in generality, C 2 motion continuity is selected in this study, where JLAP and KCSIA respectively define linear and cornering motions.

1) JERK LIMITED ACCELERATION PROFILE
JLAPs are C 2 continuous motion profiles that connect two points by providing acceleration, constant velocity and deceleration phases while obeying boundary conditions and restrictions on jerk, acceleration and velocity [33], [49]. Fig. 3 shows the JLAP motion phases. A JLAP jerk profile is defined, where j max,k is the k th axis maximum jerk magnitude in the motion. The time intervals t 0,l ≤ t < t 3,a , t 3,a ≤ t < t con and t con ≤ t < t l are respectively the acceleration, constant velocity and deceleration phases. The total motion duration is obtained as a sum of time intervals that are computed according to kinematic constraints (15) and boundary conditions (18). By successive integration of (19)), the k th axis total displacement where the path length s = r l 2 . v s,l,k and a s,l,k are the respective k th axis velocity and acceleration components at the start of the linear segment. Thus, a JLAP can be optimized by describing a parameter vector For JLAPs, kinematic constraints (15) are implemented as: with a maximum k th axis velocity v con,m,k = v s,l,m,k + a s,l,m,k T 1,a,m + T 2,a,m + T 3,a,m along the m th linear segment.

2) KINEMATIC CORNER SMOOTHING WITH INTERRUPTED ACCELERATION
KCSIA is a 2D method of generating near time-optimal C 2 continuous cornering motions by analytically calculating the cornering velocities, accelerations and durations while enforcing zero acceleration and the same tangential velocities V c at motion boundaries and obeying user-specified cornering tolerances and kinematic constraints [33]. Fig. 4 shows a KCSIA demonstration, where, based on JLAPs, its jerk profile is defined with j c,k is the cornering jerk. By successive integration of (26), a cornering path is defined with where p s,k is the k th axis coordinate of the cornering path starting point. The k th axis total displacement while traversing the cornering path is obtained from (28), where T c = 2T 1,c + T 2,c is the total motion duration. The cornering Euclidean length is derived. The position vector elements at the start, middle and end of the cornering path are respectively, where the cornering error ε = p c − p mid 2 is derived as subject to (15)- (17) and −u c,m ≤ 0. This is followed by the generation of the m th cornering path, where the start/end points are obtained from (31) -(32) and the path is plotted using (28). Afterwards, motions along the line segments that connect the cornering paths are stitched using time optimal JLAPs. In order to solve (9) or (13), (33) is incorporated in the cornering constraints (16). The kinematic constraints (15) are realized in the form of The geometric constraints (17) where a e,l,m and v e,l,m are the acceleration and velocity vectors at the end of the m th linear motion. p s,l,m and p e,l,m are the start and end position vectors of the m th linear motion. All  solutions of (9) or (13) are hereafter referred to as KCSIA* solutions. From each solution, a Pareto-optimal trajectory is generated by plotting linear segment and cornering trajectories in succession according to (19) and (28), where the linear segment and cornering path variables are retrieved from µ.

3) SPACE COMPLEXITY ANALYSIS
Space complexity is described as the memory resource required to execute a computation [51]. For an optimization problem with n var variables and n con constraints, a typical sequence quadratic programming solver requires a O(n 2 var + n var n con ) long double precision working array [52]. KCSIA* consists of n var = n l dim(µ l ) + n c dim(µ c ), (38) variables, where dim(.) denotes vector dimension. According to (23) and (34), n var = 10n l − 3. With reference to (16), (24) and (35)-(37), n con = 32n l − 21. Hence, the KCSIA* memory demand is a double precision working array of length O 420n 2 l − 366n l + 72 ≈ O(n 2 l ). In contrast, KCSIA has a O(1) memory demand since it separately optimizes each corner and linear segment.

III. RESULTS
The validity of the proposed method is tested in this section. The validation process constitutes the generation of Paretooptimal trajectories for a given piecewise linear contour. Subsequently, contouring error performance experiments are conducted for selected optimal trajectories. Experimental results are compared with KCSIA.

A. OPTIMIZATION CONDITIONS AND EXPERIMENTAL SETUP
The industrial biaxial table in Fig. 5 with kinematic limits: j lim = [50000, 50000] mm/s 3 , a lim = [1000, 1000] mm/s 2 and v lim = [80, 80] mm/s is used for verifying the proposed method. The cornering tolerance at each corner is specified as 0 < ε ≤ 200 µm. A sequence quadratic programming implementation in MATLAB R environment is used to solve (9) on a laptop computer having the specifications: core i7 intel processor, 2.50 GHz CPU, 8 GB RAM and Windows 10 operating system. The minimum trade-off level δ min = 0.02 is used in the divide and conquer algorithm.
The experimental setup in Fig. 5 is used to validate the optimization results. Motion along each axis is actuated by computer controlled AC rotary servomotors, where linear motion is acquired via ball screws. 76.29 nm resolution rotary encoders sample table position at a 5 kHz rate. Table velocity is computed as a sampled position numerical difference, where 20 Hz and 60 Hz low-pass filters suppress noise effects for the x 1 and x 2 axes, respectively. Reference trajectories are fed into a desktop computer, having Intel (R) Core i7-3770K CPU, 3.50 GHz, 8 GB RAM and Ubuntu 15.04 64 bit operating system in a Xenomai 3.0 real-time framework, which receives encoder data and sends control signals to the AC rotary servomotors.
In this study, the system dynamics of the industrial biaxial table are considered as a second-order model [53] (39) with

Mẍ(t) + Dẋ(t) + Fsgn {ẋ(t)} = u(t)
where M, D and F are diagonal matrices of axial inertia term, viscous friction and Coulomb friction, respectively. u ∈ R 2 and x ∈ R 2 are the control input and axial position vectors, respectively. sgn {ẋ} is a function that returns a vector with the respective signs of theẋ elements. In a similar manner to [54], the tracking error dynamics    are the proportional and derivative gain diagonal matrices. The controller (41) is used as the control law for conducting experiments with identified plant parameters (see Table 1), k P,k = 7225 s −2 and k D,k = 170 s −1 for ∀k. In order to verify the effectiveness of the proposed method, an experiment is conducted on the industrial biaxial table. KCSIA and ζ = {0, 0.25, ζ * , 0.75, 1} KCSIA* motion profiles are fed into the experimental setup as reference trajectories, where 10 trials are conducted to check for the repeatability of results for each trajectory. In this work, the contouring error estimation method in [40] is used, where contouring errors are computed with respect to the original piecewise linear tool paths. In this study, star-shaped (Section III-B) and complex (Section III-C) tool paths are considered for optimal trajectory generation and experimental verification.

B. CASE I: STAR-SHAPED TOOL PATH
The star-shaped tool path in Fig. 6 is considered for generating Pareto optimal trajectories. This path consists of acute and obtuse corners to check the validity of the proposed method in generating trajectories that trade off cycle time with corner smoothing.

1) OPTIMIZATION RESULTS
The KCSIA* Pareto frontier is shown in Fig. 7, where KCSIA is a dominated solution. The time-optimal KCSIA* solution  has an inferior solution since it maximizes cornering velocities, consequently maximizing the cornering Euclidean lengths as shown in (31) and (34). This results in KCSIA having a reduced cycle time at the cost of a high total cornering Euclidean length (i.e., high cornering errors), while KCSIA* considers both objectives and provides a better performance.
The best trade-off KCSIA* solution is located at ζ * = 0.53125. It offers 2.709 mm less corner smoothing than the KCSIA* time-optimal result while being 0.137 s faster than the minimum corner smoothing result. This corresponds to achieving 53.086 % and 53.006 % of the available cycle time and corner smoothing saving potentials. The optimization results are summarized in Table 2. Fig. 8 shows the KCSIA* jerk, acceleration and velocity profiles for x 1 and x 2 axes at different ζ values. All profiles are within the pre-defined kinematic constraints (15). As ζ → 1, T cycle → T min , where at the time-optimal KCSIA* has a shorter T cycle than KCSIA. The KCSIA and KCSIA* tangential velocities are depicted in Fig. 9(a). At both acute and obtuse corners, KCSIA has higher cornering velocities than KCSIA* for ζ = 1. This is a consequence of (34), where the upper bound cornering constraints are activated to maximize velocity as shown in Fig. 9(b). For KCSIA*, the cornering velocities and ε m depend on ζ , where both parameters increase as ζ → 1 and vice-versa. The correspondence between ε m (Fig. 9(b)) and the generated cornering paths is shown in details 1 and 2 of Fig. 6.
The computation time for generating the Pareto frontier (Fig. 7) at δ min = 0.02 is 188.865 s, where it takes 0.861 s to compute the anchor points (7) and 188.865 s to perform the divide and conquer recursions. The average computation time per Pareto-optimal solution is 6.092 s. The choice of δ min affects computation time and the approximation of the best trade-off solution since it is a measure of Pareto frontier resolution. Table 3 shows influence of δ min on the best tradeoff solution and computation time.

2) EXPERIMENTAL RESULTS
KCSIA* has lower error peaks at acute cornering instances than KCSIA and KCSIA* error tends to decrease as ζ → 0 at obtuse corners (see Fig. 10). The KCSIA* contouring errors are within the pre-defined cornering tolerance (i.e., 0 < ε ≤ 200 µm). The maximum and average contouring errors of each trajectory for all the trials are illustrated in Fig. 11, where it is shown that KCSIA* has a better performance than KCSIA. Contouring error tends to increase as cycle time is reduced (i.e., ζ → 1). This is attributed to increased cornering velocities (see Fig. 9(a)).
The maximum contouring errors for ζ = 0 are nonzero minima ∀ζ even though the cornering errors ε m of the reference trajectory are approximately zero ( Fig. 9(b)). This is the result of vibrations caused by non-linear frictional characteristics such stick-slip and Stribeck effects in pre-sliding regimes and near zero velocity instances, respectively [55].
A 1σ standard deviation is used to validate the contour error results, where Fig. 12 shows the consistency of KCSIA*  in performing better than KCSIA and the ζ * KCSIA* trajectory having the best trade-off between cycle time and contouring error. Table 4 summarizes the experimental results of contouring performance.

C. CASE II: COMPLEX TOOL PATH
A relatively complicated tool path in the shape of a butterfly (Fig. 13) is used to study the effectiveness and limitations of  the proposed approach. The path consists of 50 corner points that are interconnected with linear segments having lengths ranging between 0.3 and 6.1 mm.

1) OPTIMIZATION RESULTS
Pareto-optimal solutions of KCSIA* are represented in Fig. 14. The KCSIA solution is infeasible since it consists of overlapping cornering paths as shown in detail 2 of Fig. 13. This is because KCSIA produces cornering motions in isolation from preceding and succeeding corner points. On other hand, as depicted in detail 1 of Fig. 13, KCSIA* avoids cornering path overlap by considering the geometric constraints (17). The best trade-off solution is located at ζ * = 0.5, where it reduces corner smoothing by 39.856 % while it increases cycle time by 15.022 % relative to the time-optimal KCSIA* solution. A summary of the optimization results is shown in Table 5. At δ min = 0.02, the computation time required to generate the Pareto frontier (Fig. 14) is 80,011.571 s, where the anchor points are computed in 633.501 s and the divide and conquer recursions are done in 79,378.070 s. The average computation time per Pareto-optimal solution is 2963.392 s.
The generated jerk, acceleration and velocity profiles for different ζ values are within the specified kinematic constraints as shown in Fig. 15. The axial velocity constraints are not activated due to the proximity of corner points. Fig. 16  shows the extent of corner smoothing at different ζ values. As ζ → 1, cornering errors tend to approach the upper bound of the tolerance in order to reduce cycle time.

2) EXPERIMENTAL RESULTS
Contouring error profiles for different ζ values are depicted in Fig. 17, where the maximum errors tend to be within the pre-defined cornering tolerance (i.e., 0 < ε ≤ 200 µm). Fig. 18 shows the maximum and average contouring errors of each trial. Similar to Section III-B2, the maximum contouring errors tend to decrease as ζ → 0, where non-zero minima  are obtained at ζ = 0. The impact of vibrations at corner points with ε m ≈ 0, by which the motion must stop once and cause a larger tracking error, is shown in Fig. 18(b), where the average contouring errors increase as ζ decreases beyond ζ * . Pareto-optimal trajectories which are generated at ζ < ζ * have a higher number of corner points with ε m ≈ 0 than those generated at ζ ≥ ζ * (Fig. 16).
The repeatability in the maximum and average contouring error results is shown by the 1σ standard deviations in Fig. 19. The best trade-off KCSIA* trajectory shows the best compromise between cycle time and maximum/average contouring error. A summary of the experimental results is depicted in Table 6.

IV. DISCUSSION
This study aims to solve the bi-objective problem of generating trajectories that minimize cycle time and maximize VOLUME 9, 2021 motion accuracy for piecewise linear contours. The results indicate that the two objectives are contradictory, as illustrated in Figs. 12 and 19, where the best trade-off solutions are selected as the ones that are nearest to the unattainable utopia solutions. This study also depicts a correlation between corner smoothing and contouring error (Figs. 7 and 12 in Section III-B and Figs. 14 and 19 in Section III-C), where this result is in agreement with the findings of [34], [40].
Contrary to KLCS approaches proposed in [21], [33], [35], [36], [38], [39], [56], that achieve near time optimality, where motions at each smooth corner and linear segment are computed separately, the proposed approach generates Pareto optimal motions for the entire smoothed path. This contrast allows the proposed method to have better results in accuracy and cycle time while KCSIA produces a dominated solution as shown in Figs. 7 and 12. The proposed method can also avoid cornering path overlaps while maintaining Paretooptimality (detail 1 in Fig. 13).
The results (Figs.12 and 19) signify that the proposed approach can be used to provide a trade-off between cycle time and motion accuracy within user-defined cornering tolerances. Although Pareto-optimal reference trajectories that favor accuracy can be generated, the minimum achievable contouring error finally depends on the feed drive system controller.
In contrast to [21], [33], [35], [36], [38], [39], [56], which have a O(1) space complexity, the proposed method has a O(n 2 l ) space complexity, where n l is the number of linear segments. This is because, in [21], [33], [35], [36], [38], [39], [56], motions are optimized along each corner and linear segment separately, while the proposed method optimizes motion for the entire trajectory. The computation time of the optimization process is also dependent on the number of linear segments, where it rapidly grows as the number of linear segments increases. Although significant computation time is required for the complex toolpath in Section III-C, there still exist many mechanical parts consisting of simple and moderated number of linear segments. Hence, the proposed method is effective for such parts.

V. CONCLUSION
This paper proposes a method of Pareto-optimal corner smoothing to trade off between cycle time and motion accuracy for industrial feed drive systems with piecewise linear contours. The total cornering Euclidean length is used as a motion accuracy representative in the bi-objective optimization problem formulated by the normalized normal constraint approach, where kinematic limits, continuity conditions and user-specified cornering tolerances are described as constraints. JLAPs describe linear motions, while a modified KCSIA defines smooth corner profiles. A Pareto frontier is generated by the divide and conquer algorithm, where the solution that is nearest to the utopia point is selected as the best trade-off solution. The proposed method's effectiveness is verified via experiments. Relative to KCSIA, the best tradeoff solution reduces the maximum and average contouring errors by 47.53 % and 25.40 % while it increases cycle time by 2.53 %.
The proposed method is limited to symmetrical line-toline corner transitions. Pareto-optimal corner smoothing with asymmetrical line-to-line, line-to-arc and arc-to-arc transitions is left as future works. The proposed method considers tool tip position in trajectory generation. In order to extend the method to 5-or 6-axis CNC machine tool applications, smooth transitions of tool orientation are required at corners. This consideration shall be done in future works. The proposed method requires significant computational resources for complex toolpaths with many corner points. Computational cost reduction is left as future works. The energy consumed by the generated Pareto-optimal was not scrutinized in this study. For future considerations, energy saving can be assimilated by extending the studied bi-objective problem to produce a Pareto surface of optimal solutions. Since the proposed method does not include actual contouring error measurements in the optimization process, the highest achievable cornering accuracy is limited by the feed drive system controller. Feed forward compensation strategies [8]- [10] can be incorporated in the optimization process to overcome controller limitations and improve the proposed method's performance. This research avenue will be considered as future works.

APPENDIX
Algorithm 1 shows a pseudo code of the optimization process described in Section II. A candidate solution ψ m is added into a Pareto set only if it passes removal [47] and significance criteria. For a MOOP with m o objectives, the removal criterion detects non-Pareto regions by checking whether any of the first m o − 1 elements of a permuted Lagrange multiplier vector if ζ l < 1 then 26: ψ u ← the point proceding ψ l in 27: ζ u ← the first ζ > ζ u in Z 28: else 29: break / * Terminate the divide and conquer loop * / µ * ← map ψ * in the decision space 36: return , ψ * , µ * 37: end function ψ m , with neighboring points ψ l = ψ l,1 , ψ l,2 and ψ u = ψ u,1 , ψ u,2 , is considered to be a significant Pareto-optimal point if min |ψ l,1 − ψ m,1 |, |ψ l,2 − ψ m,2 |, |ψ u,1 − ψ m,1 |, |ψ u,2 − ψ m,2 | ≥ δ min , (43) VOLUME 9, 2021 where δ min is a user-specified minimum trade-off level. The choice of δ min affects the resolution of the generated Pareto frontier [42]. Since the removal and significance criteria cannot distinguish between local and global Pareto regions, a Pareto filter [45] is implemented to retain only global Pareto points once the divide and conquer loop is terminated.