A Novel Resolution Scheme of Time-Energy Optimal Trajectory for Precise Acceleration Controlled Industrial Robot Using Neural Networks

: The surging popularity of adopting industrial robots in smart manufacturing has led to an increasing trend in the simultaneous improvement of the energy costs and operational efﬁciency of motion trajectory. Motivated by this, multi-objective trajectory planning subject to kinematic and dynamic constraints at multiple levels has been considered as a promising paradigm to achieve the improvement. However, most existing model-based multi-objective optimization algorithms tend to come out with infeasible solutions, which results in non-zero initial and ﬁnal acceleration. Popular commercial software toolkits applied to solve multi-objective optimization problems in actual situations are mostly based on the fussy conversion of the original objective and constraints into strict convex functions or linear functions, which could induce a failure of duality and obtain results exceeding limits. To address the problem, this paper proposes a time-energy optimization model in a phase plane based on the Riemann approximation method and a solution scheme using an iterative learning algorithm with neural networks. We present forward-substitution interpolation functions as basic functions to calculate indirect kinematic and dynamic expressions introduced in a discrete optimization model with coupled constraints. Moreover, we develop a solution scheme of the complex trajectory optimization problem based on artiﬁcial neural networks to generate candidate solutions for each iteration without any conversion into a strict convex function, until minimum optimization objectives are achieved. Experiments were carried out to verify the effectiveness of the proposed optimization solution scheme by comparing it with state-of-the-art trajectory optimization methods using Yalmip software. The proposed method was observed to improve the acceleration control performance of the solved robot trajectory by reducing accelerations exceeding values of joints 2, 3 and 5 by 3.277 rad/s 2 , 26.674 rad/s 2 , and 7.620 rad/s 2 , respectively.


Introduction
As industrial robots are applied in more and more manufacturing scenarios, high efficiency and low energy cost of the motion trajectory are essential to meet the production needs and reduce resource consumption [1][2][3]. Time-energy cost-optimal trajectory planning for robots given a task path is conducted while surrendering limitations of kinematic and dynamic capabilities for safe operation. Effective trajectory optimization resolution of redundant robots with compliance to single-level physical constraints, such as velocity, is easy to implement. However, existing model-based algorithms for robot path planning may result in joint-angle drift and non-zero final joint velocity and acceleration when systems are subjected to comprehensive velocity, acceleration and torque constraints and controlled at multiple levels, as discussed in [4][5][6].
Recently, many kinds of solver software, combined with various optimization methods, have been utilized in robot planning research to avoid obtaining infeasible solutions. The objective and constraints input into solver software should follow the specified rules of optimization theories, such as convex optimization (CO), linear programming (LP) [7], quadratic programming (QP) [5], second-order cone programming (SOCP) [3], concaveconvex procedure (CCCP) [8], and semi-definite programming (SDP). Each algorithm of the aforementioned software toolkits has its strict solution form and is usually solved by Lagrange duality and the interior point method [9], also known as the barrier method or combined methods. Moreover, the solving process of the above software is implicit, which is not conducive to estimating a global optimum or local minima. Thus, for the trajectory optimization problem of robots with precise multi-level constraints, it is necessary to establish a generalized time-energy consumption optimization model and to develop algorithms with feasible solutions in actual scenes.
To address the issues, this paper proposes a novel multi-objective trajectory optimization approach with the goal of simultaneously minimizing the time and energy consumption of industrial robots. First, we formulate the optimization problem by establishing the weighting function of time and energy consumption in generalized path coordinates. The transformed trajectory variables in the weighting function corresponding to the angle, velocity and acceleration in joint coordinates are considered as generated optimizing variables. Then, we present forward-substitution interpolation functions as basic functions to calculate indirect velocity and acceleration expressions introduced in the discrete optimization model and obtain segmented smooth paths in terms of trajectory variables. Consequently, we impose an iterative metaheuristic scheme to solve trajectory variables in the complex optimization problem based on artificial neural networks, avoiding the fussy conversion of strict convex or linear functions. Moreover, velocity, acceleration and torque constraints controlled in joint performance parameter limits in the overall process of motion are also the focus of the trajectory design and optimization in this paper, according to the actual needs of industrial robots. The proposed trajectory optimization model can adjust the variable values according to the requirement of high-speed motion or low energy consumption. Finally, extensive experiments are conducted to verify the effectiveness of the proposed optimization solution scheme by comparing it with the state-of-the-art trajectory optimization methods using Yalmip software.

Related Work
Generally, to solve a trajectory planning problem, joint angles of an industrial robot, given the desired path described in Cartesian space, are first calculated by an inverse kinematic model. Then, trajectory variables, including velocity and acceleration as functions of angle derivative in terms of time, are optimized to achieve the minimum running time or energy consumption with the demands of multi-level control for safe operation. Feng et al. [10] employed fifth-order polynomial expressions to construct joint trajectories and established a time-energy consumption optimization model by introducing velocity variables obtained from the pre-defined trajectory. Since an industrial robot usually has six joints, the optimization objective of all joint angles is a multidimensional expression. For simplification of the objective function, Palleschi et al. [11] described path coordinates as a normalized trajectory variable based on phase plane theory to reduce the high-dimensional state space and solve the tracking problem of a minimum-time path with smooth and continuous accelerations by the SCIP solver of Opti Toolbox. Although software can generate a set of optimized solutions, the implicit solving process is often difficult to converge when multi-level constraints of the robot are required to be controlled.
For the control of robots with optimal trajectory under the multi-level constraints enforced, a conventional solution is the Jacobian pseudo-inverse method, which can be regarded as an analytical solution with a direct correlation between end-effector and joint trajectory variables. Ramezani and Williams [12] present an optimal redundancy resolution technique by using the augmented Jacobian to provide feasible solutions for the minimum objective function. The traditional Jacobian pseudo-inverse scheme is just applicable to trajectory optimization with joint physical limits expressed as equality equations, while inequality constraints can hardly be introduced to this direct solution scheme.
To achieve efficient redundancy resolution for multi-level control of a robot's trajectory, Verscheure et al. [3] transformed the optimization model with consideration of time, energy cost and smoothness into a convex optimization problem, and presented a transcription method to solve a SOCP with inequality constraints using robust numerical algorithms in Yalmip. Nagy and Vajk [7] converted a profile generation model into linear programming (LP) form and solved it with a sequential optimization method. Steinhauser and Swevers [4] presented a two-step iterative learning algorithm that compensates for inevitable model-plant mismatch of time-optimal motion, which improved tracking performance and ensured feasibility based on a sequential convex log barrier method. Zhang et al. [6] formulated the trajectory resolution problem subject to joint angle, velocity, and acceleration constraints as a QP to overcome computationally intensive pseudo-inverse-based schemes. The optimization models referred to require the objective and constraints to be converted to satisfy the strict solving forms, which sometimes induces a failure of duality and obtains results exceeding limits.
Among optimization approaches, metaheuristic algorithms have shown their capabilities for searching for near-global-optimal solutions to numerical real-valued test problems, such as the genetic algorithm (GA) [13,14], particle swarm optimization (PSO) [15], artificial neural networks (ANNs) [16,17], and so forth. Chen et al. [5] formulated a multi-level simultaneous minimization scheme as a dynamical quadratic program (DQP), which was solved by a piecewise linear projection equation neural network. Abu-Dakka et al. [14] constructed smooth joint trajectories with cubic polynomial functions as basic functions in a segmented path to establish an optimization model of minimum time, energy and jerk. Undetermined parameters in basic functions were solved by a parallel-populations genetic algorithm (PPGA) procedure. Based on the multi-objective genetic optimization algorithm NSGA-II, Shi et al. [18] solved the optimization problem in the multi-objective form.
Although the above state-of-the-art metaheuristic algorithms have good effects on solving general optimization problems, some of these methods applied in trajectory planning cannot guarantee the initial and final zero-velocity and zero-acceleration and joint-angle drift sometimes occurs, so are not suitable for actual conditions of the robot operation. In addition to the solution algorithm, the established trajectory optimization model and the form of internal basic function for iteration calculation, dynamic equation expressions, and other factors will also have an impact on the output trajectories results. Thus, it is necessary to establish a multi-objective optimization model with consideration of the applicable conditions in the actual situation.

Robot Kinematics and Dynamics Modeling
In this paper, a serial 6-axis robot is considered to establish the time-energy optimization model, the structure of which is shown in Figure 1. The modified Denavit-Hartenberg (MDH) parameters are presented in Table 1. Forward and inverse kinematic models of the specified robot can be established with these parameters.
Dynamic equations of the robot with several identified parameters are established to express the joint torques NE τ i according to the iterative Newton-Euler formulation [19]. To establish an accurate dynamic identification model, inertia terms r τ i and friction terms f τ i are introduced in the torques expressions [20] as follows: where . q i and .. q i denote the velocity and acceleration of joint i, I ai is the inertia moment for rotor and gears of actuator i, f vi and f ci are the viscous and Coulomb friction coefficients. It should be noted that the inertia tensor C I i related to the center of mass of link in NE τ i should be converted in terms of the inertia tensor A I i related to the origin of the joint coordinate by introducing the expression which contains position coordinate components of the center of mass of link (x Ci , y Ci , z Ci ), inertia moment components of link i (I xxi , I yyi , I zzi , I xyi , I xzi , I yzi ). The differential Equation (1) in terms of 78 dynamic parameters can be linearized and expressed by τ = ω q, where τ = [τ 1 , τ 2 , . . . , τ 6 ] T is the joint torques vector, Ω = [Ω 1 , Ω 2 , . . . , Ω 6 ] T is the standard parameters vector, ω q, q is a 6 × 78 regressor matrix and q is the joint position vector. The above equation can be transformed into a simplified form with a minimal set of 52 identifiable parameters Ω min by regrouping the original parameters Ω, a detailed derivation process for which can be seen in [21,22]. Thus, the torque can be simplified as where ω min q, q . We sample several data of torques and joint angles when the robot moves according to an exciting trajectory at t = t 1 , t 2 , . . . , t N , and over-determined linear equations are obtained as where observation matrix is Ψ = , the torques vector is Γ = τ(t 1 ) T , τ(t 2 ) T , . . . τ(t N ) T , and ε is a vector of sampling error. Based on the least square method, the minimum identification parameters set can be obtained by The calculated minimum set of combined parameters Ω min is listed in Table A1 in Appendix A. The dynamic model with an iterative Newton-Euler formulation is obtained by substituting Ω min into Equation (1).
Then, we rewrite the dynamic model into the state-space equation form, as follows: where M(q) is the mass or inertia matrix, C(q, . q) is the vector of Coriolis and centrifugal terms, G(q) is the vector of gravity terms, and f τ ( q is the vector of friction terms.     I  I  I  I  I  I  I  I Figure 1. Specified robot structure.

Problem Formulation
This section formulates the time-energy optimal trajectory planning problem of a robot with precise multi-level constraints. The optimization objective and complex constraints, including velocity, acceleration and torque, are derived based on the phase plane theory and Riemann approximation method. For effective numerical computing, the optimization model introduces the discrete kinematic and dynamic equation results.

Original Formulation
Minimizing the total motion time t f of an industrial robot is simply expressed as t f 0 1dt. For optimization of the energy cost, the main resistance energy consumption of the servo circuit in each robot joint is modeled by Joule's law, as follows: where i I servo is electric current and i r is the resistance. Due to the relevant relations between torque and current, the above equation can be proportional to the expression in terms of torque: For a single servo motor system, the energy consumption model during the running time should be: which omits the coefficient i r since it is constant for one joint. Then, energy consumption i E is also proportional to the integral of the torque square t f 0 τ 2 i dt. Total energy consumption of a robot can be optimized by minimizing the normalized expression: which divides the square of maximum torque limit (τ i ) 2 for equalized optimization of every joint energy. For simultaneous minimizing of the time and energy of a robot with multi-level constraints, the following synthetical optimization model is established as where β is the weight coefficient; overline and underline of velocity, acceleration and torque represent the maximum and minimum values.

Reformulation Based on Phase Plane Theory
Due to the complexity of the multidimensional-optimal problem for an industrial robot with six joints, the optimization objective is calculated using six different trajectories functions until the minimum operation time value for any of all joints is found. For simplification of the objective function, we describe path coordinates as a parameterized trajectory variable based on phase plane theory to reduce the high-dimensional state space to one dimension.
Based on the Pontryagin principle [23], the optimization independent variable t can be replaced by the trip proportion variable assigned to the path curve s(t) through a generalized coordinate transformation t = s 0 1 . s(t) ds. Then, the reformulated optimization problem is analyzed in the S-phase plane. Joint position q is a function of s, notated as q(s). Joint velocity and acceleration are transformed as where ..
s(t) = ds dt . Moreover, substituting the relations (15) into the dynamic equation (8), yields where the term C(q, . q) is linear to the joint velocity . q i [24], which can be described as C(q(s), . q(s)) = C(q(s), q (s)) . s. Meanwhile, the trajectory velocity . s is positive when the robot runs along a path, such that the sign function of q (s) . s is equal to q (s). We define m(s) = M(q(s))q (s), c(s) = M(q(s))q (s) + C(q(s), q (s))q (s), g(s) = F c sgn(q (s)) + G(q(s)) and f (s) = F v q (s) for simplification of Equation (16).
The normalized time-energy consumption optimization model for the robot trajectory is reformulated by substituting expression (16) and t f = 1 0 1 . s(t) ds as where h(s) represents . s(t) = ds dt . Since the integral optimization model is difficult to be solved, it is converted to discrete sums based on the Riemann approximation method. Then, the corresponding discretization and normalization are carried out for the time-energy optimization objective with multilevel constraints, as follows: where s j is the j-th point proportion of the divided path curve, ∆s j+1 = s j+1 − s j . Simultaneously, velocity and acceleration in joint coordinate and torque constraints are discrete as where u( for j = 1, . . . , n − 1 is assumed as a piecewise function. In the case of the dramatic increase in velocity and acceleration, the corresponding constraints are limited in the range of 1/(n − 1) multiple maximum values at the initial and final interval segments.

Solution Algorithm
In this section, we present a solving scheme of a time-energy optimal trajectory with a forward-substitution interpolation method and metaheuristic algorithm. Firstly, highorder polynomial interpolation functions are applied as basic functions to calculate indirect kinematic and dynamic expressions introduced in a discrete optimization model and other constraints. Then, we develop a solution algorithm based on artificial neural networks to generate candidate solutions for each iteration without any conversion into strict forms, until a minimum objective function value is achieved.

Forward-Substitution Interpolation Basic Functions
A path that a robot usually moves around is specified by a series of space position coordinates P k = (x k , y k , z k ). Geometric curve functions, such as the cubic Bezier or B-spline curve, could be used to fit these points, as where control points 2 Then the interpolated path points P 11 , P 12 , . . . coordinates are calculated by the cubic Bezier curve function, and we add the number of path points to N. For convenience, we notate points as Q 1 , Q 2 , . . . , Q l , . . . , Q N . Although the partitioned Bezier curve function describes the discrete interpolated points space coordinates between the specified points along the path. Each joint positions vector q(Q l ) = [q 1 (Q l ), q 2 (Q l ), q 3 (Q l ), q 4 (Q l ), q 5 (Q l ), q 6 (Q l )] T in the joint coordinate system corresponding to the path point can be resolved by the inverse kinematic model with position coordinate (x, y, z) and specified pose orientation vector, as rotations around the x-axis, y-axis, z-axis (rx, ry, rz).
To express the joint trajectory function, isometric distributed generalized path curve parameters S 1 , . . . , S N are introduced to establish the mapping relation of unevenly distributed trajectory points q(Q 1 ), . . . , q(Q N ). Figure 2 shows the correspondence relations diagram of path points and generalized path curve parameters. Due to the velocity and acceleration requirements, the third-order polynomial function is the lowest order form we choose to be an interpolated trajectory function in terms of the curve parameter S. Moreover, the initial and final trajectory functions with zero-velocity and zero-acceleration are considered by applying the fifth-order polynomial function. The common fitting method of the fixed segmentation function just uses joint angles of two adjacent points, such as q i (Q 1 ) and q i (Q 2 ), to establish the trajectory function. In this case, the higher derivatives of the adjacent trajectory functions have saltation at the connecting point. To solve this problem, we propose a forward-substitution interpolation as a predefined basic trajectory function expressed as . . .

Dynamic Metaheuristic Optimization Algorithm
After the pre-defined value ∆s j , τ( is substituted, the discrete optimization objective v(s j ) is a function in terms of variables h(s j ). The sum terms of a normalized torque square in objective expression (17) will induce high-order nonlinear forms of independent variables h(s j ), which is difficult to be solved by the model-based analytical optimization algorithm or solver software without any conversion. Thus, we develop a solution algorithm based on artificial neural networks (ANN) [17] to handle this issue.
In an n-dimensional optimization problem, a pattern solution representing input data in the ANN, is defined as First, a starting candidate of the pattern solution matrix H is generated, which is randomly generated between the lower and upper bounds of a problem: Cost functions corresponding to pattern solutions are obtained by where f is the objective function ∆s j+1 v(s j+1 ) in expression (17).
The candidate solution with the minimum objective function value for all pattern solutions is selected as the target solution. This target solution will be updated at each iteration. After determining the target solution H Target among the other pattern solutions, the target weight W Target , corresponding to the target solution, must be selected from the population of weight (weight matrix) by the following expression: where W(o) is a matrix generating random numbers uniformly between zero to one during iterations, and o is an iteration index. The weight superscript relates to its pattern solution (e.g., 2 w is related to the second pattern solution) and the weight subscript is shared with the other pattern solutions (e.g., 2 w 3 is shared with the third pattern solution). Every pattern solution has its corresponding weight value which has been involved in generating a new candidate solution. Moreover, the sum of elements in W(o) is 1. After forming W(o), new pattern solutions H New are generated by the following expression: Therefore, the new pattern solution has been updated for iteration o + 1. Based on the best weight value called "target weight", the weight matrix should be updated, as follows: where the i-th new pattern solution . In summary, the optimization problem can be solved by the general behavior of ANN, which can be described by

Results and Discussion
To assess the general applicability and verify the accuracy of the proposed method, optimal time trajectory and time-energy synthesis optimization results were calculated using simulation and actual experiments with the kinematic and dynamic boundary conditions compared with the state-of-the-art algorithms in [3] using Yalmip software. The boundary conditions are listed in Table 2. A circle path of robot motion was applied to test the effectiveness of the proposed method, as shown in    In [3], the minimum objective Error! Reference source not conventional optimization method SOCP, which converts th ( ) into convex form with additiona as follows: The above fussy conversion of strictly convex functions is av optimization model directly solved by the ANN algorithm in S second-order polynomial interpolation functions are applied as b In [3], the minimum objective (17) is solved by the conventional optimization method SOCP, which converts the non-convex function v( convex form with additional constraint conditions, as follows: 2 2 βτ 1 s j+1 /τ 1 . . .
The above fussy conversion of strictly convex functions is avoided by the proposedThe above fussy conversion of strictly convex functions is avoided by the proposed optimization model directly solved by the ANN algorithm in Section 5. Moreover, the second-order polynomial interpolation functions are applied as basic functions to calculate indirect kinematic and dynamic expressions introduced in a discrete optimization model and other constraints in reference [3].
For verification of the accuracy and efficiency of the proposed method, the optimization variable results h(s j ) in objective function solved by SDPT3 optimization toolkit in Yalmip and the ANN algorithm, when the basic functions are the second-order polynomial interpolation functions and the initial and final velocity and acceleration boundaries are neglected as in [3], as shown in Figure 4. We select discrete points when the number of s is 20. As can be seen in Figure 4, optimization variable h(s j ) values at different coordinates solved by the proposed method are similar to that of the commercial software toolkit Yalmip, which indicates the high accuracy of our model. The convergence curve of the optimization objective solved by ANN shows the stable minimum value of the objective function is 1.1454 when the number of iterations is more than 330, which is a little less than the minimum objective value, 1.1484, obtained by Yalmip.
To intuitively display the motion of the robot, the joint angles and joint angular velocities and torques of six joints during the trajectory operation with time-energy optimization weight α = 0 are obtained according to the method in [3] and our proposed method with bounded velocity, acceleration, torque limits, and initial and final boundaries, as shown in Figures 5 and 6.   The joint angular velocity and acceleration curves start with non-zero values calculated by the method in [3], and accelerations of joints 2, 3 and 5 exceed the limits by 3.277 rad/s 2 , 26.674 rad/s 2 , 7.620 rad/s 2 , respectively. The generated trajectories are hardly used in the actual scene. On the contrary, our result showed a smoother optimization trajectory with initial and final zero velocity and acceleration and all kinematic and dynamic constraints were satisfied within the extreme values.
For verification of the effectiveness of the proposed method, the calculated normalized energy consumption 6 ∑ i=1 t f 0 τ 2 i (s)/τ 2 i dt and motion time t f results obtained by the proposed optimization model are listed in Table 3, when weights β = 0, 0.1, 1 and 10 0.4 , of which selected weights β are the same as [3]. In Table 3, the motion time t f rises and the normalized energy decreases simultaneously when the weight increases, which indicates that the proposed trajectory optimization model can adjust the variable values according to the requirement of high-speed motion or low energy consumption. The results of the optimization trajectories with different weights β = 0.1 and 1 can be seen in Figures 7 and 8. As can be seen in Figures 6-8, the trajectory acceleration fluctuation at the corresponding singular point decreases gradually with weight increases. This indicates that the design of the objective function is based on the trade-off analysis of time and energy consumption. Therefore, the overall trend of the joint torque curve of the manipulator still tends to be near the safety limit of the rated torque. However, after the energy consumption modeling and weight distribution of the manipulator servo drive control system, the operation trend of the manipulator joint torque curve tends to decrease with increase in the weight. The time near the safety limit of the rated torque gradually decreases, which also means that the energy consumption of the joint servo drive control system has been effectively improved during the corresponding trajectory operation. From the normalized comparison of energy consumption in Table 3 of the manipulator trajectory optimization, it can be seen that modeling the energy consumption index of the robot servo motor and assigning a corresponding weight design plays a significant role in reducing the energy consumption value of the manipulator trajectory operation.   The joint angular velocity and acceleration curves start with non-zero values calculated by the method in [3], and accelerations of joints 2, 3 and 5 exceed the limits by 3.277 rad/s 2 , 26.674 rad/s 2 , 7.620 rad/s 2 , respectively. The generated trajectories are hardly used in the actual scene. On the contrary, our result showed a smoother optimization trajectory with initial and final zero velocity and acceleration and all kinematic and dynamic constraints were satisfied within the extreme values.
For verification of the effectiveness of the proposed method, the calculated normalized energy consumption ( )  Table 3, when weights β = 0, 0.1, 1 and 10 0.4 , of which selected weights β are the same as [3]. In Table 3, the motion time tf rises and the normalized energy decreases simultaneously when the weight increases, which indicates that the proposed trajectory optimization model can adjust the variable values according to the requirement of high-speed motion or low energy consumption.     To test the actual optimization effect, an experiment of the proposed method was implemented with the robot system, as shown in Figure 9. The path time of the optimal trajectory when 0.4 10  = was 3.2915 s which ensured the joint velocity in the range of the velocity set limits. The energy consumption of the robot running optimal trajectory was 0.357 Wh. The experimental sampling data of six joint angles, torques and velocities To test the actual optimization effect, an experiment of the proposed method was implemented with the robot system, as shown in Figure 9. The path time of the optimal trajectory when α = 10 0.4 was 3.2915 s which ensured the joint velocity in the range of the velocity set limits. The energy consumption of the robot running optimal trajectory was 0.357 Wh. The experimental sampling data of six joint angles, torques and velocities of time-energy consumption optimal trajectory are shown in Figure 10. The torque values of the servo motor of joints were within the torque performance constraints. The safety range of the accelerations of robot joint was reached at several discrete solutions, and the obtained solutions were used in an actual scene, which demonstrated the effectiveness of the above algorithm and the constraint design.

Conclusions
To describe the optimization problem of robot trajectory planning, a time-energy optimization model based on the ANN algorithm was proposed in this paper. The energy consumption model of servo resistance loss was constructed in addition to motion time. For a six degrees of freedom industrial robot, the kinematics and dynamics impact on the coupling constraint condition were considered, including velocity, acceleration and torque. The optimization model was discretized based on the Riemann approximation method. Based on the established kinematic and dynamic model of the robot, a basic discretization model in terms of the generalized path variable mapping was constructed. Forward-substitution interpolation functions were presented as basic functions for the insurance of the initial and final zero-velocity and zero-acceleration of indirect kinematic expressions introduced in the discrete optimization model. Finally, the trajectory optimization parameters and the comprehensive tradeoff time-resistance energy loss index with multi-level performance constraints were solved by a numerical iterative solution strategy based on neural networks. The simulation and actual experiments were implemented with different optimization weights. The proposed method could enhance the acceleration control performance of the solved robot trajectory by reducing accelerations exceeding values of joint 2, 3 and 5 by 3.277 rad/s 2 , 26.674 rad/s 2 , 7.620 rad/s 2 , respectively. Moreover, comparison results between our method and recent optimization methods showed that the improved basic function contributes to the smoothness of the optimization trajectory and guarantees zero velocity and acceleration at the starting and ending points.

Appendix A
The minimum identification parameter set is shown in Table A1.