A Kinematically Constrained Reparameterization Approach to Optimal Time and Jerk Motion of Industrial Machines

The increasing demands of a modern industry require the usage of industrial machines to perform various tasks with complex profiles. Tasks are usually given by a set of motion trajectory via-points; thus the geometric path representation and the proper reparameterization of the trajectory with respect to time are necessary. This study presents a kinematically constrained reparameterization approach for generating the trajectory by formulating an optimal control problem for time and jerk. Both geometric path generation and trajectory reparameterization are adopted by B-splines to ensure the trajectory smoothness. Inequality constraints are proposed to satisfy the axial kinematic limits, and equality constraints are included to ensure zero velocity and acceleration at the start and end of the trajectory. The trade-off between time and jerk is formulated as a bi-objective optimization problem. Moreover, the normalized normal constraint method with divide and conquer algorithm is applied to generate a Pareto front with significant trade-offs. Subsequently, the best trade-off solution is chosen. The proposed method is investigated by simulation with several geometric profiles, and the effectiveness is validated with an experimental result.


I. INTRODUCTION
The development of science and technology enables the usage of robotic manipulators and computer numerical control (CNC) machines for increased production rates and precision in given tasks. Therefore, research has focused on strategies for generating a smooth and optimal trajectory that can be incorporated with advanced controllers. The geometric path of a robot usually consists of a set of points, called viapoints, in the Cartesian space, where the tasks to perform and obstacles to avoid are determined. Trajectory planning takes an input of the geometric path, kinematic, and dynamic limitations to industrial machines and generates an output as a reference trajectory either in a joint or Cartesian space. Therefore, the optimal trajectory generation of a geometric path is generally translated into an optimal control problem (OCP) while considering the equality and inequality constraints of industrial machines.
The associate editor coordinating the review of this manuscript and approving it for publication was Kai Li .
There are several approaches to solve the OCP in the literature [1]. Time-optimal trajectories were generated with multiple switching curves along a specific path [2]. Korayem et al. provided optimal trajectory planning strategies for mobile manipulators based on the Pontryagin's minimum principle [3]- [5]. Similarly, time-energy optimal trajectory generation with small tracking errors was proposed in [6]. A method for determining the maximum allowable load for mobile manipulators with redundancy constraints was presented in [7]. Posa et al. applied a direct optimization method to generate the optimal trajectory of rigid body systems with environmental contact [8]. The pseudospectral method [9] and fast direct multiple shooting algorithm [10] were successfully used to solve the OCP in robotics. However, to simultaneously satisfy geometric positions of given via-points and desired smoothness of a contour profile, a smooth interpolation or approximation function with sufficient continuity is required to parameterize the motion trajectory of industrial machines.
Splines are piecewise polynomial functions widely used to generate smooth trajectory primitives to ensure acceleration or jerk continuity [11], [12]. The spline representation in terms of a compact basis form known as B-spline gives a simpler computation and a precise representation of smooth free-form curves and surfaces [13]. The interpolation or approximation of via-points is independent of the B-spline order. For general usage of splines, the required geometric path is usually generated in a parametric form, and the specific reparameterization of the curve is applied to modify the motion trajectories [14]- [17]. However, these methods initially require computing the arc length of the parametric curve to determine the required time and velocity profile. The curve parameter is then updated in each time step. As a drawback, the excitation of the machine vibration may occur due to the inaccurate mapping between the curve parameter and the arc displacement. Therefore, the spline-based trajectory generation problem is formulated as an OCP using the spline properties for constraint satisfaction.
An optimal criterion for increasing the productivity of industrial needs is the generation of time-optimal trajectories [18]- [20]. Wang and Horng [21] generated near minimum-time trajectory planning for robot manipulators by adopting a cubic B-spline with velocity, acceleration, and jerk constraints on every intermediate point. Duan and Okwudire [22] suggested a trajectory that smoothly transversed along the corner segments with a required tolerance in an optimal time. Mercy et al. [23] proposed a method for minimizing the motion time of the CNC machine's tool while considering the generated workpiece accuracy. In [24], the combination of trajectory planning with cubic splines in the Cartesian space and B-splines in the joint space produced time-optimal and jerk-continuous trajectories.
The generation of trajectories with minimum jerk has been investigated [25]- [27]. High-jerk trajectories induce significant vibrations in actuators, which affect the tracking performance of control algorithms [28]. The jerk minimization problem is contradictory to the time minimization because the trajectory reaches the extremities of the kinematic limits when time is kept minimum. Based on this concept, time-jerk optimal trajectory generation was investigated [29]. Gasparetto and Zanotto [30], [31] formulated the bi-objective optimization problem for the trade-off between the total execution time and the jerk of the trajectory, where the jerk square integral was considered as an objective function. In [32], Gasparetto et al. compared the results with the global jerk minimization algorithm and stated the advantage of a faster convergence time while considering the kinematic constraints along the trajectory.
As regards the abovementioned methods, optimal trajectories were generated based on changing the geometric shape to achieve the desired smoothness or time criteria. One exception is the work of Hashemian et al. [33], which presented a method for minimizing the total jerk of the trajectory while maintaining the geometric shape at a predefined time. In this method, the geometric path was reparameterized in time by a B-spline, and the nonlinear relationship between the curve parameter and time was formulated by jerk min-imization. However, time must be considered as a priori, and the kinematic constraints along the trajectory were not addressed.
This study proposes a kinematically constrained splinebased reparameterization approach for the trajectory generation of industrial machines while maintaining the geometric shape. It is a direct approach to parameterizing the trajectory in terms of splines and formulates an OCP as a bi-objective optimization problem. The properties of splines are exploited to ensure the kinematic constraint satisfaction for all times. Two objectives are considered herein: total time and jerk square integral of the trajectory. A parametric curve with jerk continuity is used to represent the via-points in the Cartesian space; therefore, the predefined geometric path is obtained. As in [33], the reparameterization function in terms of the B-spline is adopted to have a nonlinear relationship between the curve parameter and time. Fig. 1 depicts an overview of the B-spline reparameterization and bi-objective optimization method. Considering the total time as an unknown parameter, the computation for the jerk square integral of the reparameterized trajectory is very expensive; hence the trapezoidal method is used to estimate the jerk value at each partitioning interval.
A finite set of inequality constraints are determined based on the convex hull of the reparameterization function derivatives and the maximum parametric curve derivatives to satisfy the kinematic limits along the trajectory. Moreover, adding equality constraints at the start and end of the reparameterization function allows smooth start-and end-transitions (i.e., zero initial and final velocities and accelerations). For the bi-objective optimization problem, the process of revealing the Pareto front comprising trade-off solutions between time and jerk is implemented by applying the divide and conquer algorithm (D&C) [34]- [36] with normalized normal constraint (NNC) method [37], [38], where each solution is computed using the sequential quadratic programming (SQP) [39]. This achieves an efficient representation of the Pareto front and determines the constrained reparameterization of the trajectory with a significant trade-off. In contrast with the previous studies, the main contributions of this paper are summarized as follows: 1. A kinematically constrained B-spline function is adopted for the reparameterization of the predefined geometric path (parametric curve) with respect to time to generate time-jerk optimal trajectories. Thus, compared with the related studies in [29]- [32], the geometric path executed from interpolation or approximation of via-points remains unchanged; therefore, it is effective for industrial product design that requires a precise presentation of geometric profiles, e.g., computer-aided design (CAD). 2. Motion time is not considered as a priori, and the trade-off solutions between time and jerk are formulated as a bi-objective optimization problem. A finite set of equality and inequality constraints are proposed using the spline properties in the reparameterization function. Therefore, compared to the B-spline reparameterization implemented in [33], optimal motion time is achieved, and the semi-infinite kinematic constraints along the trajectory are satisfied for all times. Moreover, the proposed method ensures zero start-and end-velocities and accelerations of the trajectory.
The rest of the paper is organized as follows: Section II briefly describes the representation of the given via-points as a predefined geometric path by B-spline interpolation or approximation; Section III adopts the reparameterization function as a B-spline; Section IV proposes the method for limiting the kinematic constraints for the reparameterized trajectory; Section V states the bi-objective optimization scheme for time and jerk, followed by the application of the NNC method and the D&C algorithm; Section VI presents the effectiveness of the proposed method as validated by simulation and experimental results; and Section VII draws the conclusions.

II. REPRESENTATION OF VIA-POINTS AS A GEOMETRIC PATH BY B-SPLINES
B-splines are widely used for the geometric modeling of curves and surfaces due to their flexibility. The construction and properties of the B-splines are discussed in detail [40]. Given via-points, the following two types of curve fitting methods are considered to represent the geometric path: interpolation, where the curve passes through all via-points, and approximation, where the curve passes near the via-points by minimizing the error between them. The parametric B-spline of order k consisting of piecewise polynomial functions of degree k −1 and n+1 control points is determined as follows: with a non-decreasing knot vector where, s(u) represents the position of the geometric path with parameter u; p j denotes the three-dimensional position coefficient; and the knot vector U consists of n + k + 1 knots that are k-times clamped at both ends. The first and last position coefficients coincide with the first and last viapoints, respectively. The distribution of the inner knots in (2) reflects the interpolation or approximation curve to be fitted [40]. B j,k (u) is the basis function defined by the Cox-de Boor recursion formula as follows: The r th derivative of the parametric curve is determined as follows: with its r th derivative of the basis function given by . (5) The unknown position coefficients for the parametric curve can be solved in (1) by knowing the basic functions and the via-points to be interpolated or approximated. Thereafter, the parametric curve derivatives can be determined using (4). In this study, the sixth-order parametric B-splines are used to represent the geometric paths to generate trajectories with a continuous jerk.

III. TRAJECTORY REPARAMETERIZATION A. REPRESENTATION OF REPARAMETERIZATION FUNCTION BY B-SPLINE
After the geometric path is defined as a parametric curve, a particular reparameterization method is generally adopted for obtaining motion trajectories with respect to time. Reparameterization refers to velocity, acceleration, and jerk vectors modification without changing the geometric shape. More specifically, the curve parameter u is expressed in terms of time t by a strictly monotonic increasing function u(t). Therefore, the curves s(u) ands(t) are geometrically the same in position, but are kinematically different functions.
The conventional reparameterization method of the parametric curve is constant scaling, u(t) = λ t, where λ is a constant. In this case, the velocity, acceleration, and jerk vectors are obtained by multiplying the respective r th derivatives of the parametric curve with λ r . However, this method does not guarantee a smooth trajectory, which has the lowest possible jerk value, especially at the start and end [13]. The reparameterization function for the parametric curve is adopted herein in terms of a piecewise continuous function by a B-spline [33].
The q th order B-spline reparameterizaion function with m + 1 control points is given as follows: with its knot vector where, c  (7) are assumed to be uniformly distributed between the interval [0, t f ]. Therefore, the unknown total time and position coefficients of the reparameterization function are determined through the bi-objective optimization procedure presented in Section V.

B. CONVEX HULL PROPERTY OF B-SPLINE
The convex hull property of B-spline states that the piecewise segment of the curve lies within the convex hull of its coefficients. Therefore, the segment of the reparameterization function u(t), t ∈ [t i , t i+1 ] must exist within the convex hull formed by the position coefficients c pos i−q+1 , . . . , c pos i and the change in these points locally affects the function shape. The reparameterization function u(t) is a B-spline; thus, its velocity, acceleration, and jerk are also B-splines defined as The velocity, acceleration, and jerk coefficients of the reparameterization function must be determined separately to define the kinematic constraints for the trajectory. The velocity of the reparameterization function can be reconstructed as follows by discarding the first and last values of the knot vector in (7):u with Similarly, the acceleration of the reparameterization function is given byü with The jerk function is defined as follows: ... with where, c vel i , c acc i , and c and represent the velocity, acceleration, and jerk coefficients of the reparameterization function. Fig. 2 depicts the formation of those velocity, acceleration and jerk coefficients, which are the convex hull vertices of (9), (11), and (13), respectively.

IV. CONSTRAINED KINEMATIC REPARAMETERIZATION
After the adoption of the B-spline reparameterization of the parametric curve s(u), the velocity, acceleration, and jerk of the reparameterized trajectorys(t) can be expressed in terms of time by the chain rule as follows: where,ṽ(t),ã(t), andj(t) are the respective velocity, acceleration, and jerk of the reparameterized trajectory required to satisfy the kinematic limits of industrial machines. The initial and final velocity and acceleration of the trajectory must be zero to obtain a smooth transition at the start and end of the trajectory. Hence, the boundary conditions for velocity and acceleration are added as equality constraints of the reparameterized trajectory in (15) and (16) Moreover, the inequality constraints of the velocity, acceleration, and jerk must be determined such that the reparameterized trajectory satisfies the kinematic limits. Here, the change in the velocity, acceleration, and jerk of the reparameterization function affects the kinematic values of the trajectory. Therefore, the inequality constraints are proposed based on the convex hull of the reparameterization function derivatives and the maximum parametric curve derivatives. For instance, the semi-infinite velocity constraints of the reparameterized trajectory where, v min and v max denote the minimum and maximum symmetric velocity limits, are satisfied by proposing the following finite set of constraints for (15): where, |ṡ(u)| max is the maximum absolute first derivative of parametric curve, and v lim denotes the absolute velocity limit for the trajectory. Velocity and acceleration coefficients of the reparameterization function are required to define the acceleration constraints for (16). Since the acceleration coefficient c acc i is the rate of change of two velocity coefficients, c vel i and c vel i+1 , there are two combinations for each i th acceleration coefficient of the reparameterization function, which are required to satisfy the acceleration limit for (16) as follows: where, |s(u)| max is the maximum absolute second derivative of parametric curve, and a lim denotes the absolute acceleration limit. Similarly, the jerk coefficient c jerk i is the rate of change of c acc i and c acc i+1 , which again consist of three velocity coefficients, namely c vel i , c vel i+1 , and c vel i+2 , respectively. Therefore, four combinations for each i th jerk coefficient of the reparameterization function are required to satisfy the jerk limit for (17) as follows: where, | ... s (u)| max and j lim represent the maximum absolute third derivative of the parametric curve and the absolute jerk limit of the trajectory, respectively.

V. BI-OBJECTIVE OPTIMIZATION APPROACH A. OBJECTIVE FUNCTIONS AND CONSTRAINTS
This section formulates an OCP as the bi-objective optimization problem for time and jerk of the reparameterized trajectory. The objective functions are adopted with two contradictory terms: the total time and the jerk square integral of the reparameterized trajectory. Moreover, the kinematic constraints of the trajectory described in Section IV are considered as equality and inequality constraints for the optimization. Therefore, the problem is formulated as follows: where  (28) where, . 2 is the Euclidean norm, and t max is the maximum time limit defined by the user. Jerk is the third derivative of position; thus, the orders for the parametric curve and the reparameterization function must be chosen as sixth or higher to have a continuous jerk profile. In other words, the F 2 computation is very expensive. Here, the jerk square integral of the reparameterized trajectory is estimated by the trapezoidal method as follows: where,j(t h ) is the jerk value at each time step t h , and N is the number of partitions of the function.

B. GENERATION OF PARETO OPTIMAL SOLUTIONS
For the bi-objective optimization, optimal solutions are presented as the trade-off solutions between total time and jerk square integral by the Pareto front. Pareto optimality conditions can be found in [41], [42]. The extrema of the objective set in (26) are formulated by minimizing each objective independently as follows: subject to (28). The solutions in (30) are mapped into the normalized objective space formulated bỹ with The anchor pointsF * 1 andF * 2 , which are the extreme coordinates in a normalized plane, are represented as (0, 1) and (1, 0), respectively by solving (31) and (32). Therefore, the optimization problem in (26) is reformulated as follows by the NNC method: subject to (28), and: with n =F * where, ω is the weighting factor adjusted to give the solutions with a specified trade-off. Thereafter, the D&C algorithm is used to find a set of Pareto-relevant solutions with a minimum significant trade-off tolerance between consecutive trade-off solutions. The best trade-off solution ϑ * is chosen by the following condition corresponding to the weighting factor ω * .

VI. EXECUTION OF ALGORITHM A. CALCULATION CONDITIONS
In this section, the proposed kinematically constrained reparameterization with a bi-objective optimization for time and jerk is investigated with two geometric paths: S-shaped and GEMINI airfoil profiles. Thereafter, the benefit of the B-spline reparameterization over the linear reparameterization is discussed with an S-shaped profile. The optimization problems in Section V are solved using SQP (''fmincon'' function) in a MATLAB R environment of Core i7-7500U CPU and 8 GB RAM laptop computer with a Windows 10 64-bit operating system.

B. S-SHAPED PROFILE
Given nine via-points in an x-y plane, the geometric path is interpolated by the sixth-order parametric B-spline with nine control points. Fig. 3 illustrates the positions of the predefined geometric path as a parametric curve of the S-shaped profile. For this example, the velocity, acceleration, and jerk limits of optimization are set as v lim = [80, 80] mm/s, a lim = [500, 500] mm/s 2 , and j lim = [20000, 20000] mm/s 3 , respectively. The maximum time limit for the trajectory is defined as t max = 6 s. The reparameterization function is chosen as the sixth-order B-spline function with 16 control points. The number of partitions N = 1000 is used to estimate the jerk square integral. The minimum trade-off tolerance for the D&C algorithm is set as 0.01. Fig. 4 represents the Pareto front consisting of significant solutions between total time and jerk square integral. The solution for the optimal time is  achieved with the highest possible jerk value that satisfies the kinematic limits of the trajectory. Similarly, the solution with the lowest jerk value is found at t f = t max . The best trade-off solution occurs at time t f = 3.63 s with the weighting factor ω * = 0.51, and approximately 75% of the highest jerk value can be reduced by using approximately 60% of t max .
Three conditions of the relation between parameter u and motion time are provided to investigate the effect of the B-spline reparameterization on bi-objective optimization; time-optimal, jerk-optimal, and best trade-off shown in Fig. 5. The generated motion can be faster or slower depending on the reparameterization function to optimize the trajectory. The parameter u is transversed slowly at the start and end of the motion for all cases to avoid abrupt changes in velocity and acceleration values, which achieves smooth start-and end-transitions. For the rest of the motion in the time-optimal case, the reparameterization function finds the fastest way for traversing parameter u from 0 to 1. For the jerk-optimal case, the function gives more flexible reparameterization of the trajectory by modifying its position coefficients to provide the lowest jerk value at t max . The reparameterization function finds an optimal solution, which considers time and jerk equally as the best trade-off case. Fig. 6 presents a comparison of the velocity, acceleration, and jerk in x-and y-axis for time-optimal, best tradeoff, and jerk-optimal trajectories. With the proposed method, the initial and final velocity and acceleration of the trajectory become zero. Moreover, the velocity, acceleration, and jerk limits are satisfied in each axis along the trajectory of the S-shaped profile.

C. GEMINI PROFILE
For the second example, the GEMINI profile with a chord length of 100 mm, which consists of 79 points, are downloaded from the UIUC airfoil database [43]. These points are adopted as the via-points. The sixth-order parametric B-spline with 12 control points is used to approximate the geometric path (Fig. 7). In this case, the curve does not exactly pass through all the via-points, but the geometric path is approximated in the sense of least squares with the total error of [123. 34, 8.17] mm 2 . For this trajectory, the axial velocity, acceleration, and jerk limits are determined as v lim = [150, 150] mm/s, a lim = [1000, 1000] mm/s 2 , and j lim = [30000, 30000] mm/s 3 , respectively. The maximum time limit is set as t max = 8 s.
The reparameterization function, number of partitions, and minimum trade-off tolerance for the D&C algorithm are same  to those in the previous example (Section VI-B). Fig. 8 shows the Pareto front representation between total time and jerk square integral of the GEMINI profile consisting of significant trade-off solutions. The best trade-off solution is found at t f = 4.52 s with the weighting factor of ω * = 0.52 and has approximately 43% of the time saving potential compared to the jerk-optimal case and approximately 81% of the jerk saving potential compared to the time-optimal case. Fig. 9(a) and 9(b) illustrate the verification of the velocity, acceleration, and jerk limits of the GEMINI airfoil profile in x-and y-axis, respectively. Zero start-and end-velocities and accelerations are also satisfied by the proposed method.

D. COMPARISON
The simulation results were compared with an S-shaped profile to investigate the effect of the B-spline reparameterization over the linear reparameterization on the bi-objective optimization. In linear reparameterization, the curve parameter is expressed as the constant scaling function of time u = λ t [13], [30]. To achieve zero starting and ending velocities and accelerations, the additional position coefficients and knots are required to impose the first and second derivatives of the geometric path to be zero at both ends such that the geometric path representation deviates from the proposed one, as shown in Fig. 10. Fig. 11 shows the Pareto front representation for total time and jerk square integral of the proposed and linear reparameterization. All the trade-off solutions of the B-spline reparameterization are superior to those of the linear reparameterization. Compared with the best trade-off cases, the proposed method is approximately 3% faster, and the value of the jerk square integral is approximately 75% lesser than the linear reparameterization. Fig. 12 depicts a comparison of motion time, velocity, acceleration, and jerk for x-and y-axis of the respective best trade-off trajectories. By the proposed method, the peak values of jerk are reduced, especially at the start and the end of the trajectory, mainly due to the flexibility of the B-spline reparameterization function obtained by the optimization of the position coefficients and the total time. Besides, the boundary conditions for the zero velocities and accelerations are satisfied without changing the geometric shape, consequently leading to a smoother trajectory with a lower jerk value and faster time.

E. EXPERIMENTAL VALIDATION
To validate the performance of the proposed method, experiments were conducted with an industrial biaxial feed drive system shown in Fig. 13. It consists of two AC servomotors which transmit the motion in x-and y-directions via ball screws. The axial limits of the feed drive system are set    The experiments were conducted seven times to guarantee the repeatability of the results. Velocity estimation was conducted by numerical differentiation of position measure-  ments. Fig. 14 depicts the velocities for x-and y-axis of the linear and proposed method compared to their reference values. It was observed that the reference trajectories were trackable by the control system, and the time of the proposed method was faster than the linear reparameterization. Since the acceleration and jerk estimations were noisy, a comparison of tracking error was conducted for both cases instead because the lower jerk trajectory produced the smaller tracking errors [28]. Fig. 15 provides the mean absolute tracking errors in x-and y-axis of linear and proposed reparameterization for each iterative result. Due to the lower jerk trajectory of the proposed method, the mean absolute tracking errors  in the x-and y-axis were reduced approximately 17.62% and 25.52%, respectively. Moreover, variations of the x-and y-axial tracking errors from the mean were approximately 17.37% and 25.66% smaller than the linear reparameteri-zation. Therefore, these results validated that the proposed trajectory generation with kinematically constrained B-spline reparameterization provided better motion conditions in both time and jerk of industrial machines.

VII. CONCLUSION
This study proposed a kinematically constrained reparameterization method for optimal trajectory generation. The sixth-order parametric B-spline was used to represent the given via-points as a predefined geometric path, and the sixth-order B-spline reparameterization function was adopted for the optimization without changing the geometric shape. This technique considered both time and jerk minimization of the trajectory by the reparameterization function modification. The velocity and the acceleration at the start and end of the trajectory were imposed as zero by considering the equality constraints; therefore, smooth starting and ending were achieved. The inequality constraints were proposed herein based on the maximum parametric curve derivatives and the convex hull formed by the velocity, acceleration, and jerk coefficients of the reparameterization function. Therefore, the optimal solution provided a guaranteed satisfaction of the kinematic limits for all times.
The Pareto front comprising significant trade-off solutions between total time and jerk square integral was revealed by the NNC method and the D&C algorithm, and each solution was computed by SQP. The best trade-off solution was chosen subsequently. The proposed method was investigated herein with different geometric paths and limits, and simulations were implemented. The comparison with the linear reparameterization proved that the proposed method generated a smoother trajectory with a lower jerk value and faster time without changing the geometric shape. Experiments were conducted to verify the effectiveness of the proposed method in the real system. Therefore, the proposed optimal trajectory generation considering time and jerk is practically applicable to industrial machines.