Time-Optimal Trajectory Planning of 6-DOF Manipulator Based on Fuzzy Control

: Currently, the teaching programming or ofﬂine programming used by an industrial manipulator can manually set the running speed of the manipulator. In this paper, to consider the running speed and stability of the manipulator, the time-optimal trajectory planning (TOTP) of the manipulator is transformed into a nonlinear optimal value search problem under multiple constraints, and a time-search algorithm based on fuzzy control is proposed, so that the end of the manipulator can run along the given path in Cartesian space for the shortest time, and the angular velocity and angular acceleration of each joint is within a limited range. In addition, a simulation model of a 6-DOF manipulator is established in MATLAB, taking a straight-line trajectory of the end of the manipulator in Cartesian space as an example, and the effectiveness and efﬁciency of the algorithm proposed in this paper are proved by comparing the execution time with the bisection algorithm and the traditional gradient descent method.


Introduction
In current industrial production, both the teaching and offline programming can set the running speed of industrial manipulators, but the running speed of the manipulator is still relatively slow in many industrial applications. This is because reducing the running speed of the manipulator can reduce the angular velocity and angular acceleration of the joints of the manipulator, thereby reducing the vibration and jitter during the operation of the manipulator, improving its operation stability, and prolonging its service life [1]. However, reducing the running speed of the manipulator also reduces its production benefits [2,3].
Research into manipulators is divided into several aspects, such as manipulator control algorithms, trajectory planning and servo drive. Trajectory planning is an important part of the design process of manipulator control systems. At present, the mainstream research direction of trajectory planning is to optimize the trajectory of manipulators, including time optimization, jerk optimization [4], energy optimization [5], and multi-objective optimization considering time, jerk and energy [6]. In addition, manipulator obstacle avoidance [7] has become an increasing focus of trajectory planning.
The main goal of this paper is to perform TOTP in Cartesian space, making the planned trajectory time-optimal and smooth. Below, the research background and research methods of the TOTP of the manipulator will be elaborated from joint space and Cartesian space.
For TOTP in joint space, so far, there are some study methods, including limiting the joint torque rate [8], expressing joint torque and joint velocity constraints as functions of path coordinates to generate velocity limit curves [9], and the CPG method based on kinematic constraints [10]. Moreover, some algorithms are used to solve for TOTP, such as the bisection algorithm [11], input-shaping algorithm [12], hybrid-improved whale optimization and particle-swarm optimization (IWOA-PSO) algorithm [13], adaptive cuckoo search (ACS) algorithm [14], genetic algorithm (GA) [3,15], firefly algorithm [16], and simulated annealing (SA) algorithm [17]. Deep learning is also used to plan the trajectory of the grasping movement of the manipulator [18], which greatly shortens the calculation time of trajectory planning.
The above methods can plan a time-optimal and smooth trajectory; however, the TOTP in joint space only allows the manipulator to perform point-to-point (PTP) tasks, such as handling, pick-and-place, and palletizing. If the end of the manipulator moves along straight lines, arcs, or free curves, it is necessary to plan a Cartesian space trajectory.
For the TOTP of Cartesian space, there are two study methods. The first considers the distance and velocity of the end effector along a specified path as the state vector and converts the nonlinear dynamic constraints of the manipulator into state-related constraints of acceleration along the path [2]. The second transforms the time-optimal path tracking problem into a convex optimal control problem of a single state [19]. On the basis of these two study methods, there is a method based on the reachability analysis theory to transform the TOTP problem and achieve efficient solutions through multiple linear programming [20], and the other method that transforms the TOTP problem into a finite-dimensional second-order cone programming problem [21]. The sequential quadratic programming method (SQP) [22] is also used to solve the TOTP of the end of the manipulator along the spline curve, taking into account the continuity of joint acceleration and jerk. However, none of the references [2,[19][20][21] consider acceleration continuity at the end of the manipulator, therefore, during the moving process, the joint torque of the manipulator will change abruptly, resulting in vibration and shaking, which affect the stability and accuracy of the manipulator.
In view of the above research background, to solve the problem of joint space trajectory planning that can only perform PTP tasks, and to solve the problem of manipulator instability caused by the sudden change in joint torque in the TOTP in Cartesian space, this paper proposes a new offline algorithm for TOTP of manipulators based on fuzzy control, which makes the end of the manipulator run with the shortest time along a given path in Cartesian space and avoids sudden changes in the angular velocity and angular acceleration of each joint, thus compensating for the shortcomings of the above research. First, the kinematic and dynamic model of a universal 6-joint industrial robot is established. Subsequently, the TOTP problem of the manipulator is transformed into a nonlinear optimal value search problem under multiple constraints, and an adaptive time search algorithm based on fuzzy control (ATSA-FC) is proposed to calculate the shortest time of Cartesian space trajectory under the constraints of the angular velocity and angular acceleration of each joint of the manipulator. Furthermore, a simulation model of the above-mentioned manipulator is established in MATLAB. Taking a straight-line trajectory of the end of the manipulator in Cartesian space as an example, the method proposed in this paper is used to calculate the shortest time of this trajectory. At the same time, two common nonlinear search algorithms are also selected: the bisection algorithm (BA) [11,23] and the gradient descent method with constant proportional coefficient (GDM-CPC) [24]. The trajectory times and execution times of these two algorithms are compared with ATSA-FC proposed in this paper to verify the efficiency of ATSA-FC.
The remainder of this paper is organized as follows. Section 2 introduces kinematic and dynamic models of the manipulator. Section 3 introduces the trajectory planning method for the end of the manipulator. Section 4 introduces the transformation of TOTP to a nonlinear optimal value search problem and three TOTP algorithms used in this paper, which are BA, GDM-CPC, and ATSA-FC. Section 5 introduces the simulation of TOTP of the manipulator based on MATLAB. Section 6 presents the conclusions of this paper.

Manipulator Kinematics and Dynamics Model
The manipulator used in this paper is a 6-DOF wrist-separated manipulator which satisfies the Piper criterion and has a closed solution [25]. The position-level kinematic model of this type of manipulator is established using the standard D-H method, and the D-H parameters table of the manipulator is shown in Table 1.
Using the above D-H parameters, the schematic diagram of the manipulator in this paper is shown in Figure 1.

Manipulator Kinematics and Dynamics Model
The manipulator used in this paper is a 6-DOF wrist-separated manipulator which satisfies the Piper criterion and has a closed solution [25]. The position-level kinematic  model of this type of manipulator is established using the standard D-H method, and the  D-H parameters table of the manipulator is shown in Table 1.
Using the above D-H parameters, the schematic diagram of the manipulator in this paper is shown in Figure 1. Forward position-level kinematic of the manipulator solves the position and attitude of the end of the manipulator relative to the base by the given joint angles. Let T i i-1 be the homogeneous transformation matrix of the connecting rod coordinate system Σ −1 to Σ . According to the D-H rule, T i i-1 is shown in Equation (1). where Therefore, the homogeneous matrix of the manipulator end coordinate system Σ n relative to the base coordinate system Σ 0 is shown in Equation (2).
This paper uses the axis/angle notation to represent the attitude at the end of the manipulator. For any rotation matrix R, it can be considered as a single rotation around an appropriate axis in space through an appropriate angle, and the axis/angle representation is shown in Equation (3). Forward position-level kinematic of the manipulator solves the position and attitude of the end of the manipulator relative to the base by the given joint angles. Let i−1 i T be the homogeneous transformation matrix of the connecting rod coordinate system Σ i−1 to Σ i . According to the D-H rule, i−1 i T is shown in Equation (1). where Therefore, the homogeneous matrix of the manipulator end coordinate system Σ n relative to the base coordinate system Σ 0 is shown in Equation (2).
This paper uses the axis/angle notation to represent the attitude at the end of the manipulator. For any rotation matrix R, it can be considered as a single rotation around an appropriate axis in space through an appropriate angle, and the axis/angle representation is shown in Equation (3).
where k is the unit vector defining the axis of rotation, φ is the angle rotated around axis k, and the pair (k,φ) is called the axis/angle representation of R [26].
Given any rotation matrix R, whose element is a ij , the corresponding rotation angle φ and axis k are shown in Equations (4) and (5), respectively.
The axis/angle notation for the rotation matrix R is not unique because the rotation angle φ about axis k and the rotation angle −φ about axis −k are equivalent, as shown in Equation (6).
If φ = 0, then R is the identity matrix and axis k is not defined at this time. Because k is a unit vector, the equivalent axis/angle representation can be represented by a single vector r, and the vector r is shown in Equation (7).
where α =φk x , β =φk y , γ =φk z . The length of vector r is the angle φ, and the direction of vector r is the equivalent axis of rotation k.
Therefore, in addition to using a homogeneous transformation matrix to represent the position and attitude of the end of the manipulator, it can also be represented by a 6-dimensional vector X e , where X e is shown in Equation (8).
where y e , z e represent the positions of the end of the manipulator, and α e , β e , γ e represent the attitudes of the end of the manipulator.
The linear velocity υ e and linear acceleration a e at the end of the manipulator are shown in Equations (9) The attitude angular velocity ω e and the attitude angular acceleration . ω e are shown in Equations (11) and (12), respectively. .
According to Equations (9)- (12), the velocity at the end of manipulator . X e and the acceleration at the end of manipulator .. X e are shown in Equations (13) and (14).
The transfer matrix between the joint angular velocity and the end velocity of the manipulator is called the Jacobian matrix J. The Jacobian matrix is a function of joint angle θ, as shown in Equation (15).

of 24
The forward velocity-level kinematic equation of the manipulator is shown in Equation (16).
θ (16) where . θ represents the joint velocity. When J is a reversible square matrix, the inverse velocity-level kinematic equation of the manipulator can be obtained from Equation (16), as shown in Equation (17).
Taking the derivation of both sides of Equation (16), the forward acceleration-level kinematic equation of the manipulator can be obtained, as shown in Equation (18). ..
When J is a reversible square matrix, the inverse acceleration-level kinematic equation of the manipulator can be obtained from Equation (18), as shown in Equation (19). .
J is the derivative of the Jacobian matrix with respect to time, as shown in Equation (20).
The dynamic equation of the manipulator [8] is shown in Equation (21). ..

Trajectory Planning
There are two main types of trajectory planning for manipulators; one is trajectory planning in joint space and the other is trajectory planning in Cartesian space [27]. Given that trajectory planning in joint space is not capable of high-precision work, this paper performs trajectory planning of Cartesian space for the manipulator, and then uses the kinematic model in Section 2 to obtain the corresponding joint-space trajectory.
The traditional trapezoidal velocity curve at the end of the manipulator is shown in Figure 2, and the acceleration curve of the trapezoidal velocity is shown in Figure 3. As shown in Figure 3, the acceleration curve of the trapezoidal velocity changes abruptly. It can be seen from Equation (19) that when the acceleration of the end of the manipulator changes abruptly, the angular acceleration of the joint also changes abruptly.
The forward velocity-level kinematic equation of the manipulator is shown in Equation (16). (16) where θ̇ represents the joint velocity. When J is a reversible square matrix, the inverse velocity-level kinematic equation of the manipulator can be obtained from Equation (16), as shown in Equation (17).
Taking the derivation of both sides of Equation (16), the forward acceleration-level kinematic equation of the manipulator can be obtained, as shown in Equation (18).
When J is a reversible square matrix, the inverse acceleration-level kinematic equation of the manipulator can be obtained from Equation (18), as shown in Equation (19).
J̇ is the derivative of the Jacobian matrix with respect to time, as shown in Equation (20).
The dynamic equation of the manipulator [8] is shown in Equation (21). (21) where M(θ) is the inertia force matrix, C(θ,θ̇) is the cordial force and centrifugal force matrix, G(θ) is the gravity term matrix, and τ is the joint torque vector.

Trajectory Planning
There are two main types of trajectory planning for manipulators; one is trajectory planning in joint space and the other is trajectory planning in Cartesian space [27]. Given that trajectory planning in joint space is not capable of high-precision work, this paper performs trajectory planning of Cartesian space for the manipulator, and then uses the kinematic model in Section 2 to obtain the corresponding joint-space trajectory.
The traditional trapezoidal velocity curve at the end of the manipulator is shown in Figure 2, and the acceleration curve of the trapezoidal velocity is shown in Figure 3. As shown in Figure 3, the acceleration curve of the trapezoidal velocity changes abruptly. It can be seen from Equation (19) that when the acceleration of the end of the manipulator changes abruptly, the angular acceleration of the joint also changes abruptly.    According to Equation (21), it can be shown that an abrupt change in the angular acceleration of the joint indicates an abrupt change in the torque of the joint motor, which causes mechanical vibration, impacting and affecting the accuracy and service life of the manipulator [28]. Conversely, a continuous change in joint angular velocity and angular acceleration causes a continuous change in joint torque. Therefore, in this paper, the Sshaped velocity curve [13] is used to replace the trapezoidal velocity curve shown in Figure 2. The S-shaped velocity curve is shown in Figure 4, and the acceleration curve of the S-shaped velocity is shown in Figure 5. The acceleration and deceleration segments of the S-shaped velocity curve are 5th order polynomials. It can be seen from Figures 4 and 5 that the S-shaped velocity curve and acceleration curve of the S-shaped velocity change continuously. Equations (17) and (19) show that the joint angular velocity and angular acceleration of the manipulator change continuously, so the torque of the manipulator also changes continuously.  In this paper, trajectory planner P was designed to plan a trajectory with continuously changing acceleration in Cartesian space. P is represented by Equation (22). According to Equation (21), it can be shown that an abrupt change in the angular acceleration of the joint indicates an abrupt change in the torque of the joint motor, which causes mechanical vibration, impacting and affecting the accuracy and service life of the manipulator [28]. Conversely, a continuous change in joint angular velocity and angular acceleration causes a continuous change in joint torque. Therefore, in this paper, the S-shaped velocity curve [13] is used to replace the trapezoidal velocity curve shown in Figure 2. The S-shaped velocity curve is shown in Figure 4, and the acceleration curve of the S-shaped velocity is shown in Figure 5. The acceleration and deceleration segments of the S-shaped velocity curve are 5th order polynomials. It can be seen from Figures 4 and 5 that the S-shaped velocity curve and acceleration curve of the S-shaped velocity change continuously. Equations (17) and (19) show that the joint angular velocity and angular acceleration of the manipulator change continuously, so the torque of the manipulator also changes continuously.  According to Equation (21), it can be shown that an abrupt change in the angular acceleration of the joint indicates an abrupt change in the torque of the joint motor, which causes mechanical vibration, impacting and affecting the accuracy and service life of the manipulator [28]. Conversely, a continuous change in joint angular velocity and angular acceleration causes a continuous change in joint torque. Therefore, in this paper, the Sshaped velocity curve [13] is used to replace the trapezoidal velocity curve shown in Figure 2. The S-shaped velocity curve is shown in Figure 4, and the acceleration curve of the S-shaped velocity is shown in Figure 5. The acceleration and deceleration segments of the S-shaped velocity curve are 5th order polynomials. It can be seen from Figures 4 and 5 that the S-shaped velocity curve and acceleration curve of the S-shaped velocity change continuously. Equations (17) and (19) show that the joint angular velocity and angular acceleration of the manipulator change continuously, so the torque of the manipulator also changes continuously.  In this paper, trajectory planner P was designed to plan a trajectory with continuously changing acceleration in Cartesian space. P is represented by Equation (22).  According to Equation (21), it can be shown that an abrupt change in the angular acceleration of the joint indicates an abrupt change in the torque of the joint motor, which causes mechanical vibration, impacting and affecting the accuracy and service life of the manipulator [28]. Conversely, a continuous change in joint angular velocity and angular acceleration causes a continuous change in joint torque. Therefore, in this paper, the Sshaped velocity curve [13] is used to replace the trapezoidal velocity curve shown in Figure 2. The S-shaped velocity curve is shown in Figure 4, and the acceleration curve of the S-shaped velocity is shown in Figure 5. The acceleration and deceleration segments of the S-shaped velocity curve are 5th order polynomials. It can be seen from Figures 4 and 5 that the S-shaped velocity curve and acceleration curve of the S-shaped velocity change continuously. Equations (17) and (19) show that the joint angular velocity and angular acceleration of the manipulator change continuously, so the torque of the manipulator also changes continuously.  In this paper, trajectory planner P was designed to plan a trajectory with continuously changing acceleration in Cartesian space. P is represented by Equation (22). In this paper, trajectory planner P was designed to plan a trajectory with continuously changing acceleration in Cartesian space. P is represented by Equation (22).
.. X e = P(t, kt, p, r, n) where t is the trajectory running time, kt is the trajectory type, p is the constraint points set of the trajectory, which is represented by homogeneous matrices, r is the ratio of acceleration and deceleration time, and n is the number of trajectory points. To plan a straight-line trajectory in Cartesian space, two constraint points are required and the distance between the two constraint points is used as the planning distance. To plan an arc trajectory, three constraint points are required, and the central angle of the arc where these three points are located is used as the planning distance.
Assuming that a straight-line trajectory is planned, the known spatial distance of the two constraint points is l, total running time is t, and ratio of acceleration to deceleration time is r.
The velocity v(x) of this trajectory is shown in Equation (23).
where v a (x) denotes the 5th order polynomial velocity curve of the acceleration segment, v b (x) denotes the velocity of the uniform velocity segment, and v c (x) denotes the 5th order polynomial velocity curve of the deceleration segment.
then v m can be written as Equation (27).
Because the first and second derivatives of v a (x) and v c (x) at 0, rt, (1 − r)t and t are both 0, l a and l c can be written as Equation (28).
Thus, v m can be written as Equation (29).
When a uniform velocity v m is obtained, the 5th order polynomial velocity planning can be performed.
Suppose the time period starting from t s to t e , the velocity of the end of the manipulator is ν(t), and ν(t) is shown in Equation (30).
There are the following six boundary conditions, The first-order derivative ν (t) and the second-order derivative ν (t) of ν(t) are shown in Equations (31) and (32), respectively.
Because A is invertible, Equation (39) can be rewritten as Equation (40).
By substituting the boundary conditions of the acceleration, uniform velocity, and deceleration into Equation (40), the velocity change curve can be obtained, and the acceleration and displacement change curves can be obtained through differentiation and integration, respectively.
Similarly, the angular velocity and angular acceleration of the attitude at the end of the manipulator only need to be planned by changing the spatial distance l to the attitude angle φ.
Through the trajectory constraint points of the trajectory planner P, the attitude matrices R s and R e at the initial and final moments of the end of the manipulator can be determined, and R e is shown in Equation (41).
where R t is the rotation matrix that changes from R s to R e , and R t is shown in Equation (42).
Substituting R t into Equations (4) and (5), attitude rotation angle φ and rotation axis k can be obtained. The rotation matrix R i corresponding to the attitude of each trajectory point, is shown in Equation (43). where φ i represents the rotation angle corresponding to the i-th trajectory point, E 3 is the unit matrix of 3 × 3, k × is the antisymmetric matrix of vector k, and k × is shown in Equation (44).

Time-Optimal Trajectory Planning Algorithm
In this paper, the trajectory running time t is used as the control variable, the joint angular velocity and joint angular acceleration of the manipulator are used as the state variables, and the TOTP problem of the manipulator in Cartesian space is regarded as an optimal-value search problem under multiple constraints. In this paper, the minimummaximum rule is used to solve the problem of multiple constraints, and avoids the situation of local optimal solution when using time-search algorithms to find the trajectory shortest time.

Problem Description of Time-Optimal Trajectory Planning
In the process of trajectory planning, if only the constraint condition of the angular acceleration of joint i is considered, then a time t can be found such that when the manipulator is running along the trajectory, the maximum angular acceleration of joint i reaches the angular acceleration constraint condition of joint i; t at this time is the shortest time that only considers the constraint condition of the angular acceleration of joint i. A block diagram of TOTP is shown in Figure 6. The joint parameters (θ, where ϕ i represents the rotation angle corresponding to the i-th trajectory point, E 3 unit matrix of 3 × 3, k × is the antisymmetric matrix of vector k, and k × is shown in Equ (44).

Time-Optimal Trajectory Planning Algorithm
In this paper, the trajectory running time t is used as the control variable, the angular velocity and joint angular acceleration of the manipulator are used as the variables, and the TOTP problem of the manipulator in Cartesian space is regarded optimal-value search problem under multiple constraints. In this paper, the minim maximum rule is used to solve the problem of multiple constraints, and avoid situation of local optimal solution when using time-search algorithms to find the traje shortest time.

Problem Description of Time-Optimal Trajectory Planning
In the process of trajectory planning, if only the constraint condition of the an acceleration of joint i is considered, then a time t can be found such that whe manipulator is running along the trajectory, the maximum angular acceleration of j reaches the angular acceleration constraint condition of joint i; t at this time is the sh time that only considers the constraint condition of the angular acceleration of join block diagram of TOTP is shown in Figure 6. The joint parameters (θ,θ̇,θ̈) after trajectory planning and inverse kinematic were compared with the joint constraints a time search was performed. The joint constraints of the manipulator are listed in Table 2; θ̇i lim represents the angular velocity constraint, θ̈i lim represents the joint angular acceleration constraint. T were 12 constraints corresponding to the 12 shortest times. Using the minim maximum rule, the maximum value of the 12 shortest times is the shortest time o trajectory. The joint constraints of the manipulator are listed in Table 2; . θ ilim represents the joint angular velocity constraint, .. θ ilim represents the joint angular acceleration constraint. There were 12 constraints corresponding to the 12 shortest times. Using the minimum-maximum rule, the maximum value of the 12 shortest times is the shortest time of the trajectory.  .
Among them, θ imax are the maximum angular velocity and maximum angular acceleration generated by the i-th joint during operation, respectively.
Let min t 1i be the shortest time that only considers the angular velocity constraint of joint i, min t 2i is the shortest time that only considers the angular acceleration constraint of joint i, and minT is the shortest time that the manipulator runs along the Cartesian space trajectory. The problem of TOTP can be described by Equations (47)-(49).
..  (49) must be calculated at least 12 times. Therefore, to reduce the amount of calculation, let minT 1 be the shortest time of the trajectory satisfying the angular velocity constraints of the six joints of the manipulator, and let minT 2 be the shortest time of the trajectory satisfying the angular acceleration constraints of the six joints of the manipulator. The shortest time min T of the Cartesian space trajectory can be expressed as Equations (50)-(52). ..
Solving min T requires calculating Equations (51) and (52) at least once, which reduces the amount of computation to 1/6 compared with using Equations (48) and (49). Considering the 12 constraints of the manipulator joints, the shortest time min T of the Cartesian space trajectory can be further expressed by Equation (53). ..
The min T can be obtained by computing Equation (53) at least once. The condition for determining whether the trajectory is time-optimal is shown in Equation (54).
Equation (53) describes the TOTP problem for the Cartesian spatial. Next, it is necessary to use a nonlinear search algorithm to determine the shortest time t of the trajectory, such that the trajectory of the joint space of the manipulator satisfies Equation (54).

Time-Search Algorithm
This section introduces three kinds of time-search algorithms, which are BA, GDM-CPC, and ATSA-FC. By judging whether the joint trajectory corresponding to the shortest time satisfies Equation (54), the validity of the trajectory is verified. By comparing the execution times of the three algorithms, the efficiency of the ATSA-FC algorithm is verified.
The BA is a widely used search method. Its computational complexity is O(log(n)). Therefore, despite the large amount of data, this search method ensures high computational efficiency [23]. The premise of using BA is that the data must be an ordered sequence, and the time series is exactly an ordered sequence, which makes it suitable for using BA. In this paper, the BA method was used to search for the shortest time of the trajectory in the time interval [t l , t r ]. The input time for the initial trajectory planner is t r . The flowchart of BA is shown in Figure 7. This section introduces three kinds of time-search algorithms, which are BA, GDM-CPC, and ATSA-FC. By judging whether the joint trajectory corresponding to the shortest time satisfies Equation (54), the validity of the trajectory is verified. By comparing the execution times of the three algorithms, the efficiency of the ATSA-FC algorithm is verified.
The BA is a widely used search method. Its computational complexity is O(log(n)). Therefore, despite the large amount of data, this search method ensures high computational efficiency [23]. The premise of using BA is that the data must be an ordered sequence, and the time series is exactly an ordered sequence, which makes it suitable for using BA. In this paper, the BA method was used to search for the shortest time of the trajectory in the time interval [t l , t r ]. The input time for the initial trajectory planner is t r . The flowchart of BA is shown in Figure 7. The algorithm first uses t r as the running time to plan the trajectory. If the planning result does not satisfy the conditions of the time-optimal trajectory, it is necessary to determine whether the maximum angular velocity and maximum angular acceleration of all joints are within the constraints of the angular velocity and angular acceleration of the joints. The judgment condition is shown in Equation (55).
If the joint trajectory satisfies Equation (55), then let t r = t, otherwise, let t l = t. Then, let t = t l + (t r − t l ) / 2 , input it into the trajectory planner as the running time of the trajectory, iterate continuously, and finally determine the shortest time t.
However, these algorithms have limitations. When the shortest time of the trajectory is not in the given time interval, the algorithm will fail and enter an infinite loop, and the time of each planning will be infinitely close to the boundary of the given time interval.
Therefore, this paper uses GDM-CPC to solve this problem. GDM-CPC is a first-order optimization algorithm that can search for a local minimum of the function. Because this paper uses Equation (54) as the judgment condition of time-optimal trajectory, there is only one joint to reach its maximum constraint, and the angular velocity and angular acceleration of the other joints are less than their maximum constraint. Therefore, the use of GDM-CPC here will not fall into the local optimal situation, and must be able to obtain The algorithm first uses t r as the running time to plan the trajectory. If the planning result does not satisfy the conditions of the time-optimal trajectory, it is necessary to determine whether the maximum angular velocity and maximum angular acceleration of all joints are within the constraints of the angular velocity and angular acceleration of the joints. The judgment condition is shown in Equation (55).
If the joint trajectory satisfies Equation (55), then let t r = t, otherwise, let t l = t. Then, let t = t l +(t r − t l ) / 2, input it into the trajectory planner as the running time of the trajectory, iterate continuously, and finally determine the shortest time t.
However, these algorithms have limitations. When the shortest time of the trajectory is not in the given time interval, the algorithm will fail and enter an infinite loop, and the time of each planning will be infinitely close to the boundary of the given time interval.
Therefore, this paper uses GDM-CPC to solve this problem. GDM-CPC is a first-order optimization algorithm that can search for a local minimum of the function. Because this paper uses Equation (54) as the judgment condition of time-optimal trajectory, there is only one joint to reach its maximum constraint, and the angular velocity and angular acceleration of the other joints are less than their maximum constraint. Therefore, the use of GDM-CPC here will not fall into the local optimal situation, and must be able to obtain a shortest time of trajectory that satisfies all joint constraints. GDM-CPC first provides an initial trajectory running time t init , and then determines whether the joint trajectory satisfies Equation (54) after obtaining the joint trajectory of the manipulator through trajectory planning and inverse kinematics. If Equation (54) is not satisfied, then searching for a new trajectory running time, and the shortest running time of the trajectory, will finally be obtained. The advantage of this algorithm is that it only needs to provide a time value greater than 0 to converge to the shortest time of the trajectory, thereby avoiding the limitations of BA.
The flowchart of GDM-CPC is shown in Figure 8.
Actuators 2022, 11, x FOR PEER REVIEW 15 of 27 a shortest time of trajectory that satisfies all joint constraints. GDM-CPC first provides an initial trajectory running time , and then determines whether the joint trajectory satisfies Equation (54) after obtaining the joint trajectory of the manipulator through trajectory planning and inverse kinematics. If Equation (54) is not satisfied, then searching for a new trajectory running time, and the shortest running time of the trajectory, will finally be obtained. The advantage of this algorithm is that it only needs to provide a time value greater than 0 to converge to the shortest time of the trajectory, thereby avoiding the limitations of BA.
The flowchart of GDM-CPC is shown in Figure 8.  where t 1i is the time at which joint i is optimized with k p t∆θ̇i each time and t 2i is the time at which joint i is optimized with k p t∆θ̈i each time.
If ∆θ̇i > 0, the current input time is small, and the maximum angular velocity of joint i exceeds its angular velocity constraint during the trajectory planning process. The time change is k p t∆θ̇i > 0, which increases the input time to reduce the maximum angular velocity of joint i. If ∆θ̇i < 0, the current input time is large, and the maximum angular velocity of joint i is less than its angular velocity constraint during the trajectory planning process. The time change is k p t∆θ̇i < 0, reducing the input time to increase the maximum angular velocity of joint i. The joint angular acceleration has the same adjustment process. Take the maximum value of t 1i and t 2i and assign it to t as the input of the trajectory planner. In the continuous iterative process, the shortest time t of the trajectory will be obtained.
where t 1i is the time at which joint i is optimized with k p t∆ . θ i each time and t 2i is the time at which joint i is optimized with k p t∆ ..
θ i > 0, the current input time is small, and the maximum angular velocity of joint i exceeds its angular velocity constraint during the trajectory planning process. The time change is k p t∆ . θ i > 0, which increases the input time to reduce the maximum angular velocity of joint i. If ∆ . θ i < 0, the current input time is large, and the maximum angular velocity of joint i is less than its angular velocity constraint during the trajectory planning process. The time change is k p t∆ . θ i < 0, reducing the input time to increase the maximum angular velocity of joint i. The joint angular acceleration has the same adjustment process. Take the maximum value of t 1i and t 2i and assign it to t as the input of the trajectory planner. In the continuous iterative process, the shortest time t of the trajectory will be obtained.
Because k p must be adjusted many times, the algorithm will have fewer convergence steps. Therefore, this paper proposes an ATSA-FC. This method adaptively adjusts k p according to ∆  Fuzzy control is a control method that combines an expert system, fuzzy set theory, and control theory, and is very different from traditional control theory based on the mathematical model of the controlled process [29]. The behavior and experience of human experts can be added to fuzzy control. Fuzzy control is practical when establishing a mathematical model for a controlled process is difficult. This paper considers a design for a first-order fuzzy controller to adjust the value of k p . First, the input linguistic variable is fuzzified. Let the input linguistic variable be ∆ϑ, where ∆ϑ is the smallest absolute value between ∆ . θ i and ∆ .. a], and divide it into five fuzzy sets, which are NB, N, ZE, P, and PB, respectively. NB stands for negative big, N for negative, ZE for zero, P for positive, PB for positive big. The membership function corresponding to each fuzzy set is a Gaussian distribution function, as shown in Figure 9. Because k p must be adjusted many times, the algorithm will have fewer convergence steps. Therefore, this paper proposes an ATSA-FC. This method adaptively adjusts k p according to Δθ̇i and Δθ̈i by using fuzzy control.
Fuzzy control is a control method that combines an expert system, fuzzy set theory, and control theory, and is very different from traditional control theory based on the mathematical model of the controlled process [29]. The behavior and experience of human experts can be added to fuzzy control. Fuzzy control is practical when establishing a mathematical model for a controlled process is difficult. This paper considers a design for a first-order fuzzy controller to adjust the value of k p . First, the input linguistic variable is fuzzified. Let the input linguistic variable be Δϑ, where Δϑ is the smallest absolute value between Δθ̇i and Δθ̈i. Let the domain of ∆ϑ be U 1 , U 1 ∈[-a, a], and divide it into five fuzzy sets, which are NB, N, ZE, P, and PB, respectively. NB stands for negative big, N for negative, ZE for zero, P for positive, PB for positive big. The membership function corresponding to each fuzzy set is a Gaussian distribution function, as shown in Figure 9.
where −a < x < a. Second, the output linguistic variable is fuzzified. Let the output linguistic variable be k p , and let the domain of k p be U 2 , U 2 ∈ [b, c], and divided into three fuzzy sets, which are S, M, L. S represents small k p values, M represents medium k p values, and L represents large k p values. The membership function corresponding to each fuzzy set is a Gaussian distribution curve, as shown in Figure 10. The expression of the membership function for each fuzzy set is shown in Equation (58).
Second, the output linguistic variable is fuzzified. Let the output linguistic variable be k p , and let the domain of k p be U 2 , U 2 ∈ [b, c], and divided into three fuzzy sets, which are S, M, L. S represents small k p values, M represents medium k p values, and L represents large k p values. The membership function corresponding to each fuzzy set is a Gaussian distribution curve, as shown in Figure 10. The expression of the function corresponding to each fuzzy set is shown in Equation (59). The expression of the function corresponding to each fuzzy set is shown in Equation (59).
where b < y < c. Fuzzy control rules are then established and fuzzy reasoning is performed. After determining the fuzzy sets of the input and output linguistic variables, fuzzy conditional statements in the form of an IF-THEN are used to establish fuzzy control rules. The fuzzy rules are as follows: When ∆ϑ is NB or PB, it indicates that the difference between the maximum joint angular velocity or maximum angular acceleration and the constraints is large. At this time, a larger k p value should be output and the convergence step should be increased. When ∆ϑ is N or P, it indicates that the difference is medium, and a medium k p value should be output at this time. When ∆ϑ is Z, it indicates that the difference is small. A small k p value should be output to reduce the convergence step and avoid repeated oscillations.
The flowchart of ATSA-FC is shown in Figure 11. The FC is the fuzzy control function, and the FR is the fuzzy control rule, and Δϑ is shown in Equation (60).
Take the maximum value of t 1i and t 2i and assign it to t as the input of the trajectory planner. In the continuous iterative process, the shortest time t of the trajectory will be obtained.

Simulation
The simulation section first sets the parameters of the simulation, and then analyses the simulation results. The FC is the fuzzy control function, and the FR is the fuzzy control rule, and ∆ϑ is shown in Equation (60).
Take the maximum value of t 1i and t 2i and assign it to t as the input of the trajectory planner. In the continuous iterative process, the shortest time t of the trajectory will be obtained.

Simulation
The simulation section first sets the parameters of the simulation, and then analyses the simulation results.

Parameter Setting of the Simulation
This paper simulates TOTP based on the MATLAB environment. The Robotics Toolbox is used to establish the manipulator.
The constraints of each joint angular velocity . θ ilim and angular acceleration .. θ ilim set in the simulation environment are listed in Table 3. The position and attitude of the end of the manipulator is set at the initial and end moments of the straight-line path, which are represented by homogeneous matrices T st and T end , respectively. T st is shown in Equation (61) and T end is shown in Equation (62).
Set the initial time of the trajectory planner P to t = 10 s, n = 1000, the ratio of acceleration and deceleration time to r = 0.3, the trajectory type to be a straight-line, and the constraint points to be T st and T end .
The linear trajectory of the end of the manipulator in Cartesian space is shown in Figure 12. The red triangle represents the starting point and the red circle represents the end point.
Set the initial time of the trajectory planner P to t = 10 s, n = 1000, the ratio of acceleration and deceleration time to r = 0.3, the trajectory type to be a straight-line, and the constraint points to be T st and T end .
The linear trajectory of the end of the manipulator in Cartesian space is shown in Figure 12. The red triangle represents the starting point and the red circle represents the end point. Using Equation (42), the rotation transformation matrix R t of the attitude at the initial and end moments is calculated, and R t is shown in Equation (63). From Equations (4) and (5), the rotation axis/angle representation of R t can be obtained, and the rotation angle φ is shown in Equation (64) Figure 13. For input linguistic variable, x∈[−1, 1], which is due to the normalization of Equations (56) and (57). In order to ensure the completeness of the membership function [30], the membership degree at the intersection of the two membership functions is 0.5, combined with the experience summarized in the simulation debugging process of this study, the σ can be set to be 0.2142. For output language variable, it has three fuzzy sets. In order to make the membership of S and L at 0.675 tend to 0, so that the output has better clarity, the σ can be set to be 0.1. membership function [30], the membership degree at the intersection of the two membership functions is 0.5, combined with the experience summarized in the simulation debugging process of this study, the σ can be set to be 0.2142. For output language variable, it has three fuzzy sets. In order to make the membership of S and L at 0.675 tend to 0, so that the output has better clarity, the σ can be set to be 0.1. Using the fuzzy rules established in Section 4, the mapping curve of the input and output of the fuzzy control is obtained, as shown in Figure 14. Using the fuzzy rules established in Section 4, the mapping curve of the input and output of the fuzzy control is obtained, as shown in Figure 14.

Results of the Simulation
The convergence curves of the controlled time of the three algorithms using BA, GDM-CPC with k p = 0.5, and ATSA-FC are shown in Figure 15.
− − Figure 14. The input and output mapping curve.

Results of the Simulation
The convergence curves of the controlled time of the three algorithms using BA, GDM-CPC with k p = 0.5, and ATSA-FC are shown in Figure 15. As shown in Figure 16, under the same initial conditions, to get the trajectory shortest time, BA required 12 times, GDM-CPC required seven times, and ATSA-FC required five times. It can also be seen that the convergence curves of the controlled time of these three algorithms have oscillation phenomena, in which the convergence step of BA at each iteration is taken as half of the updated time interval at each iteration, since BA simply adjusts the time and does not take into account the difference with the constraint. It has the largest number of convergence steps. GDM-CPC has the smallest convergence step size at the first convergence and produces the smallest oscillation amplitude. The ATSA-FC has the largest convergence step size at the first convergence, but there is only one oscillation phenomenon, and the shortest time to the trajectory is obtained by using the least number of convergence steps, which reflects the superiority of ATSA-FC.
The variation in k p with the number of convergences is shown in Figure 16. The k p value obtained from the current iteration is used to update the input time for the next trajectory planning. Because the input time of the 5th trajectory planning meets the shortest time requirement of the trajectory, k p has only four iterations. As shown in Figure 16, under the same initial conditions, to get the trajectory shortest time, BA required 12 times, GDM-CPC required seven times, and ATSA-FC required five times. It can also be seen that the convergence curves of the controlled time of these three algorithms have oscillation phenomena, in which the convergence step of BA at each iteration is taken as half of the updated time interval at each iteration, since BA simply adjusts the time and does not take into account the difference with the constraint. It has the largest number of convergence steps. GDM-CPC has the smallest convergence step size at the first convergence and produces the smallest oscillation amplitude. The ATSA-FC has the largest convergence step size at the first convergence, but there is only one oscillation phenomenon, and the shortest time to the trajectory is obtained by using the least number of convergence steps, which reflects the superiority of ATSA-FC. oscillation phenomenon, and the shortest time to the trajectory is obtained by using the least number of convergence steps, which reflects the superiority of ATSA-FC.
The variation in k p with the number of convergences is shown in Figure 16. The k p value obtained from the current iteration is used to update the input time for the next trajectory planning. Because the input time of the 5th trajectory planning meets the shortest time requirement of the trajectory, k p has only four iterations. changes with input at each iteration.
The trajectory shortest times obtained by three algorithms are listed in Table 4. The variation in k p with the number of convergences is shown in Figure 16. The k p value obtained from the current iteration is used to update the input time for the next trajectory planning. Because the input time of the 5th trajectory planning meets the shortest time requirement of the trajectory, k p has only four iterations.
The trajectory shortest times obtained by three algorithms are listed in Table 4. It can be seen from Table 4 that the trajectory shortest time planned by BA is the largest, which is 1.6260 s. The trajectory shortest time planned by GDM-CPC is 1.6248 s, which is 1.2 ms less than that of BA. The trajectory shortest time planned by ATSA-FC is 1.6237 s, which is 2.3 ms less than that of BA. Since the judgment condition for reaching the maximum joint parameter specified in Equations (45) and (46) is 99.6-100% of the joint constraints, the trajectory shortest time difference obtained by these three algorithms is very small.
The execution times of the three algorithms are measured using the timing function in MATLAB, as listed in Table 5. As shown in Table 5, the execution time of BA is 8.38 s, and the execution time of GDM-CPC is 5.26 s, which is 37.23% less than that of BA. The execution time of ATSA-FC is 4.24 s, which is 19.39% less than that of GDM-CPC and 49.40% less than that of BA, which proves the efficiency of the ATSA-FC proposed in this paper.
Using the trajectory shortest time obtained by the ATSA-FC, S-shaped velocity planning of the end of the manipulator along the trajectory in Figure 12 is performed. The change curves of the joint angle, angular velocity, and angular acceleration are shown in Figure 17.
GDM-CPC is 5.26 s, which is 37.23% less than that of BA. The execution time of ATSA-FC is 4.24 s, which is 19.39% less than that of GDM-CPC and 49.40% less than that of BA, which proves the efficiency of the ATSA-FC proposed in this paper.
Using the trajectory shortest time obtained by the ATSA-FC, S-shaped velocity planning of the end of the manipulator along the trajectory in Figure 12 is performed. The change curves of the joint angle, angular velocity, and angular acceleration are shown in Figure 17. As shown in Figure 17, the angular velocity and angular acceleration of each joint obtained by using the S-shaped velocity curve change continuously, and it can be inferred that the operation of the manipulator is stable. It can be seen from Figure 17b,c that in the process of TOTP in Cartesian space, the maximum angular acceleration of joint 3 plays a major limiting role, satisfying the maximum angular acceleration judgment condition of the joint. The angular velocity and angular acceleration of other joints do not reach their As shown in Figure 17, the angular velocity and angular acceleration of each joint obtained by using the S-shaped velocity curve change continuously, and it can be inferred that the operation of the manipulator is stable. It can be seen from Figure 17b,c that in the process of TOTP in Cartesian space, the maximum angular acceleration of joint 3 plays a major limiting role, satisfying the maximum angular acceleration judgment condition of the joint. The angular velocity and angular acceleration of other joints do not reach their constraints. This also shows that it is feasible to use the minimum-maximum rule to solve the multi-constraint problem in TOTP.
Because the shortest times of the trajectories obtained by these three algorithms are very close, the difference between the overall angular velocity and angular acceleration cannot be seen in the comparison chart, so only the local enlarged pictures at the maximum angular velocity and maximum angular acceleration of the joint are given here, as shown in Figure 18.
It can be seen from Figure 18a,b that the angular velocity and angular acceleration of the joint are negatively correlated with the trajectory time of the manipulator. Since the ATSA-FC calculates the minimum trajectory shortest time, the corresponding joint trajectory also has the largest peak.
To determine the trajectory planning effect of these three algorithms, beyond comparing the shortest time of the trajectory, it can also be measured by using the degree of TOTP. The degree of TOTP can be described by the ratio of the maximum joint parameters that plays the major limitation role in the joint constraints, and in this simulation, joint 3 s acceleration plays the major role, so the degree of TOTP can be calculated by Equation (66).     the multi-constraint problem in TOTP. Because the shortest times of the trajectories obtained by these three algorithms are very close, the difference between the overall angular velocity and angular acceleration cannot be seen in the comparison chart, so only the local enlarged pictures at the maximum angular velocity and maximum angular acceleration of the joint are given here, as shown in Figure 18. It can be seen from Figure 18a,b that the angular velocity and angular acceleration of the joint are negatively correlated with the trajectory time of the manipulator. Since the ATSA-FC calculates the minimum trajectory shortest time, the corresponding joint trajectory also has the largest peak.   As shown in Tables   Θ 3lim . Therefore, the joint trajectories planned by these three algorithms satisfy Equation (54), which proves that the shortest time of the trajectory obtained by the above three algorithms is effective, and ATSA-FC has the highest degree of TOTP.
In fact, the degree of TOTP is not only related to the algorithm itself, but also to the judgment condition's range set in Equations (45) and (46), which determines the upper and lower limits of the degree of TOTP. In the simulation, the judgment condition's range is 99.6%-100% of the maximum joint parameters, so according to Equation (54), no matter how the algorithm is run, the degree of TOTP will always be between 99.6% and 100%.
Adjust Equations (45) and (46) and simulate a different range of judgment conditions. The trajectory shortest times and algorithm execution times planned by each algorithm are listed in Table 9. Let  θ ilim , where Er is the lower limit of degree of TOTP. As shown in Table 9, under all Er conditions, the ATSA-FC had the shortest execution time. In addition, among the three algorithms, the trajectory shortest time planned by ATSA-FC is also the smallest. Thus, the superiority and effectiveness of the ATSA-FC are verified.

Conclusions
In this paper, the problem of TOTP of the manipulator in Cartesian space is studied, and ATSA-FC is proposed, so that the end of the manipulator can run along the given trajectory of Cartesian space with the shortest running time, while avoiding the sudden change in torque of each joint.
In the simulation, taking a straight-line path of the manipulator in Cartesian space as an example, BA, GDM-CPC, and ATSA-FC are used to calculate the shortest time of this