Joint Motion Planning of Industrial Robot Based on Modiﬁed Cubic Hermite Interpolation with Velocity Constraint

Featured Application: Industrial robot joint smooth movement via multiple positions with constrained velocities. Abstract: As for industrial robots’ point-to-point joint motion planning with constrained velocity, cubic polynomial planning has the problem of discontinuous acceleration; quintic polynomial planning requires acceleration to be speciﬁed in advance, which will likely cause velocity to ﬂuctuate largely because appropriate acceleration assigned in advance is hardly acquired. Aiming at these problems, a modiﬁed cubic Hermite interpolation for joint motion planning was proposed. In the proposed methodology, knots of cubic Hermite interpolation need to be reconﬁgured according to the initial knots. The formulas for how to build new knots were put forward after derivation. Using the newly-built knots instead of initial knots for cubic Hermite interpolation, joint motion planning was carried out. The purpose was that the joint planning not only satisﬁed the displacement and velocity constraints at the initial knots but also guaranteed C 2 continuity and less velocity ﬂuctuation. A study case was given to verify the rationality and effectiveness of the methodology. Compared with the other two planning methods, it proved that the raised problems can be solved effectively via the proposed methodology, which is beneﬁcial to the working performance and service life of industrial robots.


Introduction
Industrial robots are widely used in handling, painting, welding, assembly, packaging and other fields. With the development of robot technology, industrial robots are gradually applied to polishing, monitoring, automatic manufacture system and so on [1,2]. The manipulator movements become more complex and need to satisfy higher working performances. Therefore, trajectory planning of industrial robots has to meet higher requirements. The most basic task of industrial robot motion planning is to meet the requirements for end-effector displacement and velocity [3]. Moreover, favorable trajectory planning should ensure smooth movements, fewer vibrations, impacts and mechanical wear so that it can improve working performance and extend service life [4].
Motion planning of industrial robots can be divided into two types: one is in the workspace for end-effector, the other is in joint space for joints. For motion planning in the workspace, the inverse kinematics function needs to be solved repeatedly to calculate angular displacements and angular velocities of joints, which requires a large amount of calculation [5]. Furthermore, end-effector singularity possibly appears and will make joints rotate incorrectly so that the end-effector cannot complete work successfully. For motion planning in joint space, only a few inverse kinematics solutions are needed to find out the angular displacements and angular velocities at several target points of joints. Therefore, motion planning in joint space can reduce inverse kinematics calculation largely so that computation efficiency is improved, and also can avoid joints singularity and redundancy [6]. As a result, motion planning in joint space is more suitable for the pointto-point (PTP) motion of industrial robots.
The main methods of joint motion planning are reviewed below. Linear planning with parabolic transition has the advantage that joint speed remains constant in the linear phase, but sufficient transition time must be ensured, and angular acceleration is not continuous. The widely used cubic polynomial planning has the advantages of simple principle and less computation. However, when it is used to multiple points with constrained velocities, the angular acceleration at the connection points is broken, which will lead to vibrations and impacts and make the robot arm tremble. Quintic polynomial planning has continuous angular acceleration. However, quintic polynomial interpolation requires specifying in advance the angular acceleration value of the joint at each target point. If it is not preset properly, the angular velocities will fluctuate back and forth, which is not conducive to joint motion control. Frustratingly, it is difficult to give appropriate angular acceleration values in advance [7]. In addition, quintic multinomial has a higher degree than cubic polynomial so that velocities and accelerations change more intensely, and larger velocity and acceleration peaks appear. B spline interpolation is also used in motion planning [8][9][10]. However, when the cubic B-spline curve is used for multiple points with constrained velocities, acceleration is not continuous at connecting points as well. Some researchers carried out lots of studies about cubic B spline and realized G 2 continuity at the connecting knots [11,12]. However, G 2 continuity refers to curvature continuity, and continuity of second-order derivative C 2 is not fully achieved.
In order to solve the complicated joint planning with multiple points further, Xiangrong Xu proposed a 3-5-3 polynomial splicing method [13]. This method allows to specify angular velocity at the intermediate target points arbitrarily, and angular acceleration is continuous. However, when there are more than 5 knots, angular acceleration at the target knot needs to be given in advance for quintic polynomial. Because it is hard to assign appropriate acceleration in advance, the method will also cause large angular velocity fluctuation. Saeed B. Niku proposed a 4-3-4 polynomial splicing method [14], which guarantees C 2 continuity at connection points. However, only at the beginning and the end target points can velocities be specified. If angular velocities at the intermediate target points are specified, angular acceleration will have broken points and be discontinuous. P. A. Parikh et al. used seventh-order polynomials and ninth-order polynomials to carry out joint planning [15]. The joint rotation trajectory was smooth enough in high order and the joint velocity and acceleration were continuous. However, the computation of high-order polynomial was large, the velocity and acceleration were easy to fluctuate, and the peak values of acceleration and velocity were much higher.
In summary, for PTP joint motion planning via multiple points with constrained velocities, if planning with low-order functions is adopted, angular accelerations at intermediate target points are discontinuous, which will bring additional impacts and vibrations. If planning with high-order functions is adopted [16,17], the problem of acceleration discontinuity can be solved. However, it is necessary to give angular acceleration values at the target points in advance. Since it is hard to give appropriate angular accelerations in advance, angular velocities are likely to fluctuate largely, which is unfavorable to motion control. Moreover, the higher order of planning function, the more intense change of planning curves and the larger the peak value of velocity and acceleration. All these are unfavorable for the steady operation of the robot arm. In view of the problems in the plannings mentioned above, a modified Hermite interpolation is proposed in the paper. In the proposed planning, velocities at knots can be specified at will and C 2 continuity is guaranteed. Moreover, there will be no drastic fluctuations in angular velocity since the maximum order of the planning function is no more than three and acceleration preset is not needed.

Point-to-Point Joint Motion via Multiple Points with Constrained Velocities
The End-effector of the industrial robot is the final executive component performing operations. The main types of end-effector motion are PTP (point to point) movement and CP (continuous path) movement. The movement of the end-effector is the result of joint motion. So there are corresponding relationships between end-effector movement and joint motion.

Point-to-Point Joint Motion
PTP (point to point) motion of the industrial robot is suitable for spot welding, handling, loading and unloading, circuit board insertion and other operations. In PTP motion mode, the end-effector moves from a start point to an end point. PTP motion only needs to specify positions and orientations of end-effector at the start point and end point, but the trajectory and orientations between the start point and end point are not specified or required.
When the end-effector moves in PTP mode, joints rotates normally in PTP mode as well. As shown in Figure 1, take joint 1 as an example. When the end-effector wants to move from start point P a to end point P b , joint 1, which is one of the six joints, needs to rotate from the start point A to end point B. In general, other joints rotate similarly. Finally, rotations of six joints work together to form the movement of the end-effector. There is a corresponding functional relationship between the position and orientation of the end-effector and the angular displacement of each joint, as shown in formula (1). There is also a corresponding functional relationship between the end-effector velocity and the angular velocities of each joint, as shown in formula (2).
In formula (1), [x y z] and [α β γ] are position and orientation of end-effector, θ i are angular displacements of each joint. In formula (2), [v x v y v z ] and [ω x ω y ω z ] are linear velocity and angular velocity of end-effector, ω i are angular velocities of each joint, and J is the velocity Jacobian matrix. f (θ i ) and J vary in different robots. Formulas (1) and (2) are called forward motion functions, by which the end-effector motion can be obtained from joint motion.
According to the position, orientation and speed of the end-effector, angular displacement and velocity of each joint can be obtained, which is called inverse kinematics, as shown in formulas (3) and (4) below. The speeds of the end-effector at the start and end points are zero in PTP movement so that angular velocities at the start and end points of each joint are also zero. As seen in Figure 1, the angular velocities of joint 1 at the start and end points are all equal to zero, ω a = ω b = 0.
On the basis of angular displacement and velocity constraints at the start point and end point of each joint, motion plannings of each joint can be carried out and these joint motion plannings are finally synthesized into the end-effector motion, as shown in Figure 2.

Point-to-Point Joint Motion via Multiple Points with Constrained Velocities
When there are multiple intermediate target points between the start and end points [18], and the end-effector needs to pass through these intermediate points at specified speeds, we briefly call it PTP motion with constrained speed. When the end-effector moves in speedconstrained PTP mode, joints correspondingly perform velocity-constrained PTP rotation, and both angular displacements and angular velocities at intermediate target points of joints must be specified, so as to meet position, orientation and speed requirements of the end-effector at the target points.
As shown in Figure 3, the end-effector moves from start point P a to end point P d . It needs to pass through two intermediate target points P b and P c at specified speeds, and joints should rotate accordingly. Taking joint 1 as an example, joint 1 correspondingly rotates from the start point A to end point D, passing through the intermediate target points B and C at constrained velocities. Therefore, joint 1 carries out PTP motion via multiple points at constrained velocities. By solving the inverse kinematics function, all joints' angular displacements and angular velocities at each target point can be obtained from the positions, orientations and speeds of the end-effector, so that the rotation trajectory of joint 1 can be piecewise-planned, from point A to point B, from B to C, and from C to D. Other joints can carry out similar rotation plannings. Finally, the movement of the end-effector is synthesized by all joints' rotations together.

Joint Motion Planning Method
The most used joint planning methods are cubic polynomial and quintic polynomial.

Cubic Polynomial Planning
Polynomial interpolation is commonly used to interpolate between discrete points, which is consistent with PTP motion characteristics of industrial robots, so it is often used in motion planning of industrial robots. Among the polynomials, the cubic polynomial is widely used in the PTP motion planning of industrial robots. Cubic polynomial interpolation function is shown as below formula (5): where t is the time variable and a 0 , a 1 , a 2 and a 3 are four coefficients of the cubic polynomial, which need to be found out according to the boundary constraints of angular displacement and angular velocity at two adjacent target points. Put the start time, the end time, the angular displacements and the angular velocities at the two target points into formula (5) and its first derivative, respectively. After sorting, Equation (6) for solving the coefficients can be obtained.
where, t 0 is the start time, and t 1 is the end time. θ 0 and ω 0 are respectively angular displacement and angular velocity at the first target point of a joint, and θ 1 and ω 1 are respectively angular displacement and angular velocity at the next target point. θ 0 , θ 1 , ω 0 , ω 1 can be obtained through inverse kinematics. For general PTP motion, ω 0 and ω 1 are both zeros. For velocity-constrained PTP motion, ω 0 and ω 1 may not be zeros. After finding out the coefficients, the cubic polynomial planning function can be obtained, which will satisfy the boundary constraints of angular displacement and angular velocity at the target points.
Because of simplicity, clarity and small computation, the cubic polynomial is widely used in PTP joint motion planning. However, if a cubic polynomial is applied to PTP motion planning with constrained velocity, angular acceleration is discontinuous at the connecting points between every two adjacent planning segments, as shown in Figure 4c. Broken points of acceleration will cause vibrations and impacts, which is disadvantageous to the working performance and service life of the industrial robot.

Quintic Polynomial Planning
A quintic polynomial is also commonly used in the joint motion planning of industrial robots [19,20]. Quintic polynomial planning can guarantee angular acceleration continuity at the connection points between every two adjacent planning segments. A quintic polynomial function is shown in formula (7) below.
where t is the time variable. a 0 , a 1 , a 2 , a 3 , a 4 and a 5 are six coefficients of the cubic polynomial, which need to be solved according to the boundary constraints of angular displacement, angular velocity and angular acceleration at two adjacent target points. Similarly, put the start time, the end time, the angular displacements, the velocities, and the preset accelerations at two target points into formula (7) and its first and second derivatives, respectively. After sorting, Equation (8) for solving the coefficients can be obtained.
where, t 0 , t 1 , θ 0 , θ 1 , ω 0 , ω 1 have the same physical meanings as in formula (6) and ε a and ε 1 are the angular accelerations at the target points, which need to be given arbitrarily in advance. After finding out the coefficients, the quintic polynomial planning function can be obtained. As can be seen from Equation (8), as long as the angular accelerations at the connection point of two adjacent plannings are set equal, the angular acceleration continuity is consequent. When a quintic polynomial is used for joint motion planning, it is required to specify in advance the angular acceleration values at target points. If the angular acceleration value is not preset appropriately, the angular velocity will fluctuate greatly. As shown in Figure 5, due to the improper preset of angular acceleration at points A and B, the angular velocity planning curve ω 1 (t) and curve ω 1 (t) between points A and B fluctuate largely. The joint motor needs to accelerate and decelerate back and forth, which is disadvantageous to the control of joint angular velocity.

Function Transformation between Different Intervals
Since the variable interval of cubic Hermite interpolation is [0,1], and the variable interval of industrial robot motion planning is generally not [0,1], function transformation between different intervals is needed. There is an arbitrary function expressed as Given that the function is symbolically independent, this function can also be expressed as , that is to say, the variable interval changes, the function expression should be transformed into the following formulas (9) and (10) in order to keep the function values of y unchanged.
When function interval is changed from [u 0 , u 1 ] to [x 0 , x 1 ], the transformed function y = f (x) can keep the values of the initial function y = f (u) unchanged, but the first-order and second-order derivatives of the function y with respect to u are no longer equal to that of the function y with respect to x, as shown in Equations (11) and (12).
Since the variable interval of cubic Hermite interpolation is [0,1], the following Equations below can be obtained in the interval [x 0 ,x 1 ] instead of the interval [0,1] according to the above eqations. These formulas can be used in the following study.

Cubic Hermite Interpolation
Cubic Hermite interpolation is a commonly used interpolation method. There is n points to be interpolated, (x 0 , y 0 ), (x 1 , y 2 ), . . . , (x n-1 , y n-1 ). Cubic Hermite interpolation function of the i-th curve in the interval [x i-1 , x i ], whose interval is normalized, can be expressed as the following formula (17) in the normalized interval [21].
where n is the number of interpolated points. u is the variable derived from variable x, and the interval of u is [0,1] obtained by normalization of the interval [x i-1 , x i ]. y i-1 and y i are the two interpolated data points of the i-th interpolation curve.
. y i−1 and . y i are the two first-order derivatives of y with respect to x at the points y i are the two first-order derivatives of y with respect to u at the points u = 0 and u = 1.
According to Equations (13) and (14), formula (17) can be rewritten as the following formula (18), and the variable interval is changed from [0,1] to [x i−1 , x i ]. The piecewise cubic Hermite interpolation curve is shown in Figure 6. It can be seen from formula (17) or formula (18), C 0 and C 1 continuity can be guaranteed as long as positions and first-order derivatives at the connection point between every two adjacent piecewise curves are specified equal. Whatever the values of y 0 , y 1 , . . . , y n−1 and y(x). Therefore, in order to guarantee C 2 continuity of cubic Hermite curves at each knot, . y 0 , . y 1 , . . . , . y n−1 must satisfy certain conditions.

Requirement for C 2 Continuity of Cubic Hermite Interpolation Curves
The requirement for C 2 continuity at the connecting point between two adjacent cubic Hermite curves y i (x) and y i+1 (x) is shown as the following Equation (19): ..
By sorting out Equations (19)- (21), formula (22) can be obtained below. First-order derivatives . y i−1 , . y i and . y i+1 at every three adjacent knots must meet formula (22) to guarantee C 2 continuity at connection points.
formula (22) can also be rewritten as formula (23) as below.
It can be seen from formula (22), it is almost impossible that arbitrarily specified y n−1 at knots are constrained arbitrarily, second-order derivatives of cubic Hermite curves at connecting points are discontinuous. In fact, as for cubic Hermtime interpolation curves, only one of them can be satisfied, constraining first-order derivatives willfully at knots or C 2 continuity. In order to solve the matter, a kind of modified cubic Hermite interpolation with willful constraints on first-order derivatives at knots is presented while the interpolation curves are C 2 continuous.

Modified Cubic Hermite Interpolation
Firstly, take three points of cubic Hermite interpolation as an example. Here are three data points to be interpolated A(x 0 , y 0 ), B(x 1 , y 1 ), C(x 2 , y 2 ), and the first-order derivatives at knots A, B, C are assigned willfully to  In order to specify first-order derivatives willfully at points A, B, C and guarantee C 2 continuity, new points are built. As seen in Figure 7b, new point M is built between initial point A and B, and new point N is built between initial point B and C. In order to distinguish the new interpolation function from the old one, the new function is marked as r(t) instead of y(x). Newly-built knots become A(t 0 , r 0 ), M(t 1 , r 1 ), N(t 2 , r 2 ) and C(t 3 , r 3 ), and 1st derivatives at A, M, N, and C are marked respectively as r 2 must be found out.
Because the calculation of new knots' ordinates t 0 , t 1 , . . . , t m-1 will be described in detail in a later section, t 0 , t 1 , . . . , t m-1 can be considered as the knowns. Therefore, firstly explain how to find out r 1 , r 2 , . r 1 , . r 2 . As shown in Figure 7b, the interpolation curve r 2 (t) must pass by the initial knot B. According to formula (18), the second interpolation curve r 2 (t) can be expressed function (24) as below.
where t is the variable, and (t 1 , r 1 ), The initial knot B(x 1 , y 1 ) and its derivative . y 1 must meet the interpolation function r 2 (t) in formula (24), which can be expressed as the following Equation (25).
where, x 1 − t 1 is marked as K 1 , and t 2 − t 1 is marked as H 2 , shown in Figure 7c. As mentioned above, the calculation of t 0 , t 1 , . . . , t m will be described in detail in a later section, so K 1 and H 2 could be considered as the knowns. Besides, y 1 and . y 1 are certainly the knowns of initial knot B. Therefore, Equations (26)  In order to guarantee second-order continuity at connection points, positions and first-order derivatives of new knots A, M, N and C must meet formula (22). Every three adjacent new knots among four new knots need to meet the formula (22). Specifically, knots A, M and N all together should meet formula (22). After sorting out, function (28) can be obtained.
where, t 1 − t 0 and t 2 − t 1 are marked as H 1 and H 2 respectively, shown in Figure 7c. As mentioned above, H 1 and H 2 could be considered as the knowns. Since r 0 and . r 0 are the ordinate and 1st derivative of knot A, they are known constants. Hence, Equation (28) where, t 2 − t 1 and t 3 − t 2 are marked as H 2 and H 3 . H 2 and H 3 could be considered as known constants. Since r 3 and . r 3 are the ordinate and 1st derivative of knot C, they are the knowns. Hence, Equation (29) is also a linear one of the unknowns r 1 , r 2 , r 2 can be obtained. After acquiring the positions and the 1st derivatives at newly-built knots A, M, N, and C, modified cubic Hermite interpolation curves can be obtained according to formula (18). Now let us look at the situation with four initial points interpolated. As seen in Figure 8a, there are four initial knots to be interpolated A(x 0 , y 0 ), B(x 1 , y 1 ), C(x 2 , y 2 ), and D(x 3 , y 3 ). The first-order derivatives at knots A, B, C and D are assigned willfully to . y 0 , . y 1 , . y 2 , and . y 3 . Rebuild new data points. As seen in Figure 8b, new knot M is built between the first and the second knots A and B, and new knot N between the penultimate and the last knot C and D. It must be noted that two new knots P and Q need to be built between knot B and C. If only one new knot is built between knots B and C, set of Equations cannot be solved because the number of Equations is greater than that of the unknowns. Similarly, we need to find out the unknowns r 1 , r 2 , r 3 , r 4 ,  In the same way, the initial knot B(x 1 , y 1 ) and its derivative . y 1 must meet the interpolation function r 2 (t). The initial knot C(x 2 , y 2 ) and its derivative . y 2 must meet the interpolation function r 4 (t), shown in Figure 8b. Then put x 1 , y 1 , . y 1 and x 2 , y 2 , . y 2 into r 2 (t) and r 4 (t) derived from formula (18), and we can get four Equations similar to Equations (26) and (27). When it is extended to more knots, Equations (30) and (31) can be obtained after sorting and inducting.
where, n is the number of the initial knots (n ≥ 3), and m (m = 2n − 2) is the number of the newly-built knots. The subscript i is related to initial knots while subscript j is related to newly-built knots. Array y is the ordinates of initial knots, y 0 , y 1 , . . . , y n-1 , which are all knowns. The array .
y is the first-order derivatives of initial knots,  (32) and (33) below respectively, seen as in Figure 8c for graphic representation. Since all of the initial knots' abscissas x 0 , x 1 , . . . , x n-1 are given and the newly-built knots' abscissas t 0 , t 1 , . . . , t m-1 can be considered as the knowns described in the latter section, array K and array H are known constants. Hence, Equations (30) and (31) are linear ones.
Synthesizing formulas (30) and (31), there are 4(n − 2) unknowns and 2(n − 2) Equations. Hence, the set of Equations has no unique solution unless other Equations are to be added.
Similarly, in order to guarantee C 2 continuity at newly-built knots, the positions and first-order derivatives of new knots A, M, N, P, Q and D in Figure 8b must meet the formula (22). Every three adjacent new knots need to meet the formula (22). When it is extended to more knots, Equation (34) can be obtained after sorting and inducting.
where, the meanings of m(m = 2n − 2), r, After positions and first-order derivatives of newly-built knots are found out, modified piecewise cubic Hermite curves can be obtained according to formula (18). The modified curves would pass by the initial knots, and first-order derivatives at the initial knots can be specified at will while the curves are C 2 continuous.
For more than four initial points to be interpolated, every two new points need to be built between every two initial points from the second and penultimate initial point. The method is similar to the one for four initial points.
The first idea is to average t 0 , t 1 , . . . , t m−1 . Uniform abscissa knots are beneficial to the smoothness of the interpolation curves. However, averaging t 0 , t 1 , . . . , t m−1 may possibly lead to the situation that there are two initial knots existing between two new knots, seen in Figure 9. Since a piecewise curve must pass by the two initial points, too many constraints are imposed on this curve so that the curve function cannot be solved and interpolation functions cannot be found out. To avoid the trouble mentioned above, the following scheme is proposed on the basis of homogenization as much as possible. Take five initial knots for example. As can be seen in Figure 10, two newly-built knots R 2 and R 3 are built between the second and the third initial as described above. Based on the principle of average for new knots, the horizontal distance between the second and the third initial points is divided into three segments, so the abscissas t 2 and t 3 of the new knots R 2 and R 3 can be expressed as Equations (35) and (36) below. The same scheme is applied whenever there are two newly-built knots between every two adjacent initial points.
Between the first and the second initial points, is there only one newly-built knot R 1 inserted, seen in Figure 10. The horizontal distance between the first and the second initial knots is divided by new point R 1 into two segments with unequal length, so the golden ratio is a reasonable method. Consequently, the abscissa t 1 of newly-built knot R 1 can be expressed by Equation (37). The same is as for the last newly-built knot R 6 , and the abscissa t 6 of the newly-built knot R 6 can be expressed by Equation (38).
Extended to the case that the number of the initial points is n (n ≥ 4), formula (39) can be obtained. Abscissas of all newly-built knots t are acquired from the abscissas of initial knots x. If there are only three initial points (n= 3), take the first two Equations and the last two Equations in formula (39).
where n and m (m = 2n − 2) are respectively the numbers of initial and newly-built knots. The subscript i is related to the initial knots while subscript j is related to the newly-built knots.

Case Study
In order to verify the rationality and effectiveness of the proposed methodology, the following study case is employed. The usage of the proposed methodology is further described in detail during the study case analysis.
Case: In order to accomplish a task, joint 1 of the industrial robot needs to carry out velocity-constrained PTP motion. As shown in Figure 11a Table 1. Using the proposed planning of modified cubic Hermite interpolation presented above, carry out velocity-constrained PTP motion planning for joint 1.  There are four initial knots A, B, C, and D to be interpolated. As seen in Figure 11b, newly-built knot M needs to be inserted between initial A and initial B, and N between C and D, and P, Q between B and D, according to the proposed method. Therefore, the initial knots A, B, C, D are replaced by the newly-built knots A, M, P, Q, N, and D. The initial time knots become the newly-built time knots correspondingly. In order to be unified with the previous formula, the initial time knot is denoted as x, and the newly-built time knots is denoted as t. According to formula (39), the newly-built time knots are arranged from the initial time knots x = 0, 5, 15, 25 (s) and the newly-built time knots are t = 0, 3, 8, 13, 19, 25 (s) respectively with proper rounding.
We need to find out the angular displacements and angular velocities at the newlybuilt knots M, P, Q and N from the initial knots A, B, C, and D. According to formula (30), 31, 32, 33, and 34, we can obtain Equation (40) to solve out the angular displacements and angular velocities of the rebuilt knots, The values of θ M , θ P , θ Q , θ N , ω M , ω P , ω Q , ω N can be obtained by solving the linear Equations. In addition to the previously calculated amount of time, the information of the newly-built knots for joint motion planning can be listed in Table 2. According to the motion conditions at newly-built target points in Table 2, the piecewise motion planning of joint 1 can be carried out by formula (18).
The planning function of AM segment is shown in Equation (41) below.
The planning function of the MP segment is shown in Equation (42) below.
The planning function of the PQ segment is shown in Equation (43) below.
The planning function of the QN segment is shown in Equation (44) below.
The planning function of the ND segment is shown in Equation (45) below.
Validate whether the motion planning functions meet the requirements of angular displacements and angular velocities at the initial target points, as shown in Table 3. Validation results show the planning functions meet the requirements.
Validate the continuities of the motion planning functions at the connecting points between every two adjacent planning segments, as shown in Table 4. The validation results show that angular displacements, angular velocities and angular accelerations at connecting points are all continuous.
According to the motion planning functions above, the motion planning curves can be obtained, shown in Figure 12. As can be seen from Figure 12a, the angular displacement curve is very smooth and meets the angular displacement requirements at the initial target points A, B, C, D. As can be seen from Figure 12b, the angular velocity curves are also smooth, the velocity fluctuation is gentle, and meets the angular velocity requirements at the initial target points A, B, C, D. As can be seen from Figure 12c, the angular acceleration curves are continuous at the connecting points M, P, Q, N, and the whole acceleration curves does not fluctuate much.   The motion planning curves once more verify that the motion planning functions can not only meet the constrained requirements of angular displacements and angular velocities at the initial target points but also guarantee the continuity of angular acceleration and angular velocity fluctuates gently. It proves that the proposed method can be applied well to joint motion planning with velocity constraints, and avoid the problems of discontinuous angular acceleration in cubic polynomial planning and large fluctuation of angular velocity in the quintic polynomial.

Comparison and Discussion
In order to compare with the proposed planning methodology, the most commonly used PTP joint planning method, cubic polynomial and quintic polynomial, are applied to the same case study above as well. The results of the study and analysis below reflect there are more or fewer deficiencies for these two methods to be used for joint planning with constrained velocity. The detailed comparison and analysis process are as follows. Figure 13 shows the angular displacement curves obtained by the proposed methodology, cubic polynomial, and quintic polynomial respectively. As can be seen from Figure 13a,b, the angular displacement curve obtained by the proposed method is similar to the curves obtained by the cubic polynomial. Both of the curves are smooth and the fluctuation is very small. In Figure 13c, the angular displacement curves obtained by quintic polynomial fluctuate a little bit. In fact, a curve with a higher order is more likely to fluctuate than a curve with a lower order. Therefore, we should try to avoid using higher-order functions for planning.  Figure 14 is the angular velocity curves obtained by these three methods. In Figure 14a, the velocity curves obtained by the proposed method are smooth, the velocity fluctuation is gentle, and the variation trend is roughly similar to that in Figure 14b by the cubic polynomial. In Figure 14b, the overall velocity waveform by a cubic polynomial is simple and appropriate. However, there are turning points at junctions B and C, which will lead to discontinuous angular acceleration at the connection points. In Figure 14c, the overall fluctuation of the velocity curves by a quintic polynomial is obviously large. In particular, because of inappropriate preset angular acceleration at point C, the velocity near point C fluctuates back and forth sharply, which is not conducive to the control of joint rotation speed. This shows that quintic polynomial planning is likely to cause angular velocity to fluctuate too much. Moreover, the higher the order of the polynomial is, the higher the velocity peak will be, which makes greater volatility. These are the weaknesses of quintic polynomial planning.  Figure 15 is the angular acceleration curves obtained by these three methods. As can be seen from Figure 15a,c, the angular acceleration curves, obtained by proposed methodology and quintic polynomial respectively, are both continuous at connection points. Because the order of the quintic polynomial is high, the angular acceleration in Figure 15c is also easy to fluctuate. In Figure 15b, the angular acceleration curves obtained by cubic polynomial are obviously broken at the connecting points B and C. The angular acceleration discontinuity of the joint will cause vibrations and impacts, which will degrade the working performance and decrease the service life of industrial robots. It can be seen from the above comparative analysis of this case, as for joint PTP motion planning with constrained velocity, the angular acceleration of cubic polynomial planning is discontinuous, which will lead to vibrations and impacts. Angular velocity of quintic polynomial planning fluctuates greatly back and forth, which is not conducive to velocity control. The proposed motion planning methodology not only meets the constraints of angular displacements and angular velocities at target points but also guarantees the continuity of angular acceleration. The maximum function order of proposed planning is three, which effectively avoids the fluctuation caused by high-order functions.

Conclusions
For the industrial robots' PTP joint motion planning with velocity constraint, cubic polynomial planning has the problem of angular acceleration discontinuity, which will cause vibrations and impacts. As for quintic polynomial planning, angular accelerations at target points need to be given in advance and because it is hard to give appropriate angular accelerations in advance, it will lead to large velocity fluctuations which are not conducive to velocity control.
In order to solve these problems in joint planning, a kind of methodology based on modified Hermite interpolation is proposed. On the basis of Hermite interpolation, new target knots are reconstructed from the initial target knots. Then, newly-built knots instead of old initial knots are used to carry out joint motion planning. Using the proposed methodology, not only does the motion planning meet the requirements for angular displacement and angular velocity at the initial knots, but also the angular acceleration of the whole planning is continuous without broken points. Moreover, the maximum order of the planning functions is three, which avoids the velocity fluctuation caused by the higher-order functions such as quintic polynomial. Therefore, the problems existing in cubic polynomial planning and quintic polynomial planning can be solved well.
A study case is given to verify the rationality and effectiveness of the proposed method. The case study results show that, compared with the other two planning methods, the proposed methodology is applied well to PTP joint motion planning with velocity constraints, and is conducive to the working performance and service life of industrial robots.
In the future, efforts will be made on optimization investigation as continued research in the area. The basic motion planning of industrial robots, including improved motion planning, focuses on the continuity and smoothness of the trajectory. However, efficiency, energy, impact and other problems in the actual operation requirements are not considered. Thus, optimization of single or multiple objectives among efficiency, energy, the impact will be the further potential studies on the basis of the motion planning. Moreover, motion planning needs to be realized through motion control. Due to strong nonlinear coupling between each joint and end-effector, there will be disturbances, errors and other uncertain factors, so optimal control will also be considered as a further investigation to guarantee motion planning realized precisely.