Time-optimal path following for non-redundant serial manipulators using an adaptive path-discretization

Abstract The time-optimal path following (OPF) problem is to find a time evolution along a prescribed path in task space with shortest time duration. Numerical solution algorithms rely on an algorithm-specific (usually equidistant) sampling of the path parameter. This does not account for the dynamics in joint space, that is, the actual motion of the robot, however. Moreover, a well-known problem is that large joint velocities are obtained when approaching singularities, even for slow task space motions. This can be avoided by a sampling in joint space, where the path parameter is replaced by the arc length. Such discretization in task space leads to an adaptive refinement according to the nonlinear forward kinematics and guarantees bounded joint velocities. The adaptive refinement is also beneficial for the numerical solution of the problem. It is shown that this yields trajectories with improved continuity compared to an equidistant sampling. The OPF is reformulated as a second-order cone programming and solved numerically. The approach is demonstrated for a 6-DOF industrial robot following various paths in task space.


Introduction
Short duration times are highly requested to maximize the productivity of a robotic system. Therefore, time-optimal motion is desirable. Frequently, the task for a robot is defined by a path that should be followed. Examples are welding and painting robots, CNC machines as well as autonomous vehicles. Thus, several approaches for following paths in a (time-) optimal way are made in the literature. Most of the methods exploit that the motion along a prescribed path is described with a parameter σ and its time derivativeσ . This reduces the multi-dimensional state of a robot to a two-dimensional state-space in the optimization process. A simple, but computationally expensive approach is made in refs. [1][2][3], where dynamic programming is used to solve the path following problem. In refs. [4,5], a less time expensive approach is presented, where the optimal solution is found by moving along extremals in the σ −σ plane. Moreover, the (time-) optimal path following problem can be transformed to a convex optimal control problem through a nonlinear change of variables [6] which is computationally less timeconsuming than the approaches mentioned previously. A further approach is made in ref. [7], where the convex problem is reformulated as second-order cone programming. It is reported that the calculation time compared to the convex optimal control problem is reduced even further. Since latter approach comes with low computation times, it will be used to solve the (time-) optimal path following problem in this paper.
However, the prescribed path is often deduced from process requirement (e.g., welding, polishing) as end effector (EE) path without taking singularities into account. Therefore, the inverse kinematics problem has to be solved in order to receive the corresponding joint paths. Most of the approaches for solving an optimal path following problem discussed in the literature utilize an equidistant discretization C The Author(s), 2023. Published by Cambridge University Press. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited. of σ to solve the optimization problem. Thus, the question of using a non-equidistant sampling strategy in σ which copes with the presence of kinematic singularities arises.
In this paper, an adaptive sampling strategy which copes with kinematic singularities is addressed. The basic idea of such a sampling scheme is depicted in Fig. 1. An arbitrary EE-path z(σ ) with the path parameter σ and uniformly spread samplings over the whole path is assumed to be given. The inverse kinematics problem is denoted with q = f −1 (z), and the solution of this problem leads to the corresponding joint path q(σ ). Due to the nonlinear relation f −1 , the distance of two corresponding joint values is not equidistant anymore. Especially near kinematic singularities latter values tend to vary drastically, due to the ill-posed inverse kinematics problem. This distance can be measured by the joint arc length s. With the knowledge of the nonlinear map s = s(σ ) between s and σ , the joint path and thus the EE-path can be expressed as q(s) and z(s). Prescribing a constant discretization of the joint arc length s, the discrete values of q are spread equally over the new path parameter s in the interval [0, s T ]. Since s is a measure of distance in the joint coordinates and the distance of two consecutive joint values is kept constant ( s = const.), the reverse effect near singularities is applied to the sampling of the EE-path. Thus, the drastically varying joint values with σ as parameter induce an accumulation of samples in the EE-path near kinematic singularities with the parameter s, as illustrated in Fig. 1. The resulting non-equidistant sampling in σ is derived by The adaptive sampling scheme is derived from the joint arc length parameterized inverse kinematics, derived in ref. [8]. Therefore, the inverse kinematics problem for the (1) parameterization of an EE-path z(σ ) in terms of the path parameter σ [9,10] and (2) the parameterization of the joint motion q(s) in terms of the (kinetic) arc length s is discussed briefly. Each of these parameterizations gives rise to a particular inverse kinematics problem. The two parameterizations refer exclusively to the geometry of motion and are independent of the temporal execution of the motion.

Forward kinematics
Denote with q ∈ S n the vector of joint variables and with z ∈ SE(3) the EE-pose. In the course of the entire paper, the EE-path is supposed to be given and parameterized in terms of a path parameter σ (t) ∈ [0, 1]. Because the joint motion is derived from the EE-motion, the joint coordinates are also parameterized in terms of the path parameter σ . For motion planning and control, the latter is typically used. It should be mentioned that the relation of path parameter in workspace and joint space is not unique in general. With slight abuse of notation, the same symbol is used throughout the paper. The geometric forward kinematics map determines the EE-pose for prescribed joint coordinates q [11]. The EE-twist, containing the translation and angular velocity, is determined by the velocity forward kinematics where J(q) = ∂V/∂q is the geometric Jacobian and() = d()/dt the total time derivative. Since the EEtwist is linear in the joint rates and asq = q σ , the velocity forward kinematics (2) can be defined in terms ofσ . This then yields to with () = ∂()/∂σ denoting the geometric derivative. The latter results in V(σ ,σ ) = h(σ )σ and the relation where h(σ ) := ∂V/∂σ is the partial path velocity. The relation (4) is referred to as the velocity forward path-kinematics problem and plays a central rule in the singularity-robust adaptive discretization. The EE-path z(σ ) and the partial path velocity h(σ ) originates from a path planning. Further, the path speedσ and thus V(σ ,σ ) are determined by solving an optimal path following problem as described in Section 4.

Inverse path kinematics
The classical inverse kinematics problem is to determine joint coordinates for given an EE-pose z.
For particular robots, this can be solved analytically, for example, wrist-partitioned serial manipulators [12]. However, it can only be solved numerically in general. Various different computational methods were discussed in refs. [12][13][14][15], including nonlinear root-finding, fixed point iteration, or integration of the joint velocity. This paper pursues the approach of integrating the joint velocity numerically. This amounts to solving the relation (2) for q(t) with given EE-twist V(t). The EE-twist is determined by the partial path velocity with (3) when the path is parameterized with σ (t). Then, the curve, and thus q, is uniquely determined by (4), where the partial path velocity h(σ ) describes the motion along the path.
The velocity inverse path-kinematics problem determines q for given partial path velocity h. A solution to the velocity inverse path kinematics is where J(q) + is a generalized pseudo inverse [16,17]. The geometric inverse path-kinematics problem is to determine the joint motion q(σ ) in terms of the path parameter σ . Solving the initial value problem (IVP) (5) with q(0) = q 0 as initial value corresponds to solving the inverse path-kinematics problem. When solving the IVP (5) numerically, numerical drifts are unavoidable. These drifts can be reduced by stabilizing a error system. Thus, the partial path velocity error is introduced as with h d denoting the prescribed partial path velocity. Therewith, the system (5) is amended as where (7b) governs the error dynamics. The latter asymptotically converges to zero for K > 0 [18]. The geometric error e in (7b) comprises the position error and orientation error and can be expressed as where r d is the prescribed EE-position and R d = [n d s d a d ] the prescribed EE-orientation, with columns n, s, a [19]. The skew-symmetric (cross product) matrix of a vector x ∈ R 3 is denoted with x ∈ R 3,3 . The velocity inverse path-kinematics problem can be solved numerically stable using time integration schemes, with the use of latter definitions and appropriate K. Shortly, a parameterization of the joint coordinates q(σ ) in terms of the parameter σ is attained with the solution of (5).

Joint arc length inverse path kinematics
The inverse path-kinematics problem is ill-posed near singularities. This leads to numerical problems and also causes extremely fast and undesirable joint motions. Following a path with constant speedσ 0 , for instance, leads to an EE-motion with a well-defined bounded EE-velocity. Particularly at singularities of the forward kinematics, the corresponding inverse kinematics solution q(t) will not have well-defined bounds and the joint rateq(t) tends to infinity. As a consequence, when prescribing the EE-motion by an equidistant sampling of the path parameter σ , the 'distance' of joint coordinates q corresponding to two consecutive values of σ tend to vary drastically. A measure of 'distance' between two consecutive joint values is given by the arc length s defined by the differential line element ds = √ dq T Wdq with a metric W (using the generalized mass matrix M of the robot would yield the kinetic line element and ensure homogeneity). The arc length gives rise to a well-defined discretization of the curve in joint space in terms of the step size. For the rest of this paper, W is chosen as identity matrix since the considered robot has revolute joints only.
Expressing the EE-motion z(s) in terms of the arc length s and thus derive parameterizations σ = σ (s) and q = q(s) is referred to as the arc length path-kinematics problem. It is known from the theory of metric curves [20] that there is no such explicit relation in general. Similar to the inverse path kinematics, Section 2.2, latter problem is formulated as an ordinary differential equation (ODE). Integrating the line element ds with given joint path q(σ ), the arc length in joint is [21] s = σ 0 q W dξ (9) with q W := q T Wq denoting the joint path speed. Introducing the jacobian (of above equation) ∂s/∂σ = q W in (4) yields with q = J + h. This ODE system can be solved with given initial conditions σ 0 and q 0 for s = 0 and a given parameter interval s ∈ [0, s T ]. However, this requires the knowledge of the terminal arc length s T , which is unknown a prior. The terminal arc length has to satisfy the condition Latter condition can act as a termination criterion to terminate the integration. System (10) is then solved numerically, until σ (s) = 1 is satisfied. This can be carried out numerically stable using integration schemes with zero-crossing event detection using advanced root-finding methods. In practice, the integration is terminated when σ > 1. It should be noted that the solution σ (s) is only used to determine the terminal arc length s T .
The numerical drift of the arc length parameterized inverse path-kinematics subsystem (10) again has to be stabilized. As previously discussed for the standard geometric inverse path-kinematics problem, this is attained by introducing the path velocity error with ∂h d /∂s = h d / q W . The system (10) is amended as where (13b) governs the error dynamics of the error (8). With K > 0, asymptotically convergence to zero is ensured. For numerical solving, cases where q = 0 occurs have to be taken into account, due to the division by q W . Since q = J + h two cases arise. First, the case when h = 0 is hold. This is avoided by assuming a regular EE-path so that always h = 0. Secondly, the partial path velocity h is an element of the null space h ∈ ker J + . This results in EE movements which do not have any effect on the joint coordinates, which is a contradiction to the forward kinematics.

An adaptive singularity-consistent sampling scheme
Based on the velocity inverse path-kinematics problem, the adaptive sampling approach is developed. The EE-path with parameter σ is assumed to be equivalent to the curve parameterized with parameter s. As mentioned before, prescribing the EE-motion by an equidistant sampling σ inevitably causes drastically varying joint coordinates q, of two consecutive samples, when solving the inverse path-kinematics problem (5). The arc length in joint space (9) serves as a measure of the 'distance' between two consecutive values of q and is a nonlinear relation of s and σ . For a 6-DOF Comau Racer3 robot following a straight line in workspace Fig. 2(a) with constant orientation and speedṡ 0 , the nonlinear mapping s(σ ) with an equidistant sampling s is shown in Fig. 2(b) for illustration purposes. As the drastically increase of the joint arc length indicates, the manipulator passes near a singularity at The joint coordinate q evolution and the nonlinear relation between s and σ are both preserved by parameterizing the prescribed path using the joint arc length s, Section 2.3. Prescribing a constant arc length speedṡ 0 and equidistant sampling for s for following the trajectory, a varying speedσ due to the nonlinear mapping (9) is induced. As a consequence of the relation, smaller speedsσ are achieved at higher rates ∂s/∂σ leading to more sampling points in σ (locally). Therefore, lower joint velocities are expected in such regions sinceq = q σ . Furthermore, this leads to an adaptive sampling strategy for σ derived from the nonlinear s(σ ). The latter discretization scheme can then be applied for solving optimal path following problems, see Section 5.
In the following, a comparison of the two presented inverse path kinematics is made. Therefore, the functions σ (t) and s(t) are assumed to be given with constant time derivativesσ (t) =σ 0 andṡ(t) =ṡ 0 , respectively. To get comparable results between the two parameterizations, the terminal time t T is fixed. For t T = 5 s, the EE follows the prescribed path with (1) constant speedσ 0 = 0.2 s −1 and equidistant  sampling in σ when using the velocity inverse path kinematics (5), and (2) constant arc length speeḋ s 0 = 1.482 rad/s −1 with equidistant discretization in s for the arc length parameterized velocity inverse path kinematics (10). The resulting joint paths using both inverse path kinematics are depicted in Fig. 3. It can be seen that both strategies lead to the same evolution of the joint coordinates. The arc length parameterized solution also shows an increased number of samples at σ ∈ [0.15, 0.25], which indicates the passing near a singularity. Viewing at the evolution of q 5 , it is observed that the value is going towards zero in this region. This indicates the vicinity to the wrist singularity.
As discussed above, for fixed step size s (which means fixed time step size due to constantṡ 0 ), higher sampling rates of the joint coordinates in this area are induced, and hence, lower joint velocitiesq = q σ are expected, which is illustrated in Fig. 4. This is particularly pronounced forq 4 . The joint velocitẏ q 4 using (10) is more than ten times lower than the solution obtained with the standard inverse path kinematics (5) at σ ∈ [0.15, 0.25]. Moreover, the resulting trajectories are smoother compared to using (5) for the inverse path kinematics. Apparently, this is a consequence of the accumulation of samples at σ ∈ [0.15, 0.25], as indicated previously.

Time-optimal path following problem
In this section, the optimal path following problem is described. The equations of motion (EOM) as M(q)q + g(q,q) = τ (14) describe the dynamic behavior of the manipulator, where M(q) is the mass matrix and τ are the joint torques. The Coriolis, centrifugal, and gravitational terms as well as the Coulomb and viscous friction are described by the vector g(q,q). If not necessary, the dependencies are omitted for simplicity and readability.
The objective of an optimal path following problem is finding a time profile σ (t) for a prescribed path in either the joint variables q(σ ) or the EE z(σ ), while satisfying some constraints. Latter are, for example, the dynamic behavior (14), limiting the joint velocityq, accelerationq and jerk ... q . Therefore, a description depending only on path allocated variables, such as geometrical derivatives () and the path speedσ (σ ), is desired. Given an EE-path z(σ ), the geometrical joint derivatives have to be computed with the inverse path kinematics, as discussed in Section 2.2. The first derivative q is given by (5). The higher derivatives are derived by higher order inverse path kinematics, that is, geometrical derivatives of (4), which yields For the optimization, it is beneficial to introduce the squared path speed z :=σ 2 . By assuming the path is followed only in forward directionσ ≥ 0, a one to one map betweenσ and z exists. Therefore, no case separation has to be taken into account. Introducing the variable z leads to the relations aṡ for higher time derivatives. With these relations and the higher order inverse path kinematics, the joint velocity, acceleration, and jerk can be expressed as followṡ The path parameterization of the EOM (14) is derived by introducing (17a) and (17b) as with a = 1/2Mq and b including the projected Coriolis and centrifugal terms as well as Mq . A pathdepending approximation of the Coulomb friction and gravitational forces are composed in c, and d denotes the path projected viscous friction. For receiving time-optimal solutions, the terminal time t T has to be minimized. The time-optimal path following (OPF) problem is formulated as with the corresponding path description of the objective function [22], the joint constraints (17), and the EOM (18). The lower and upper bounds of the constraints are denoted by (), (). Solving the above problem yields an optimal squared path speed profile z * . The trajectories of the joint velocity, acceleration, jerk, and torque, then can be computed using the relations (17) and (18). The optimization problem can be solved using for example, dynamic programming [1], numerical integration [23], a sequential convex programming (SCP) approach [24], or with second-order cone programming (SOCP), and a SCP approach for non-convex constraints [7].

Optimization framework and robotic manipulator
We used an Intel(R) Core(TM) i5-9500 CPU @ 3.00 GHz CPU running on Windows10 for all computations. For a 6-DOF Comau Racer3 robotic manipulator, the OPF problem (19) is solved in SOCP form. Neglecting the jerk and torque constraints temporary, (19) is a convex optimization problem. This convex problem is reformulated as SOCP and solved with a SCP approach for considering the joint jerk and torque constraints, as presented in ref. [7]. In order to reduce the optimization problem (19) to a finite discrete problem, the path parameter is discretized with N values σ i , i ∈ {0, . . . , N − 1}. Also, a linear interpolation profile for the squared path speed z in an interval [σ i , σ i+1 ) is used [6]. For more details, the reader is referred to the corresponding literature, since rewriting the problem is not topic of this paper [6,7,25]. The SOCP form of the optimization problem was implemented in Matlab using YALMIP [26]. As solver, the conic solver MOSEK [27], a tailored solver for SOCPs, was picked. For simplicity, the upper and lower bounds of (19) are chosen to be symmetrical that is, () = −(). These values are listed in Table I.
In order to solve the optimization problem, a discretization strategy for σ has to be chosen. Therefore, (1) an equidistant sampling σ i = const. which acts as benchmark strategy and (2) a non-equidistant sampling σ i = const. derived from the singularity-consistent sampling scheme, Section 3, is used. Thus, the non-equidistant discretization is given by For the rest of the paper, a uniform discretization in s ( s i = const.) is used to compute the non-equidistant sampling σ i = const. Figure 2(b) exemplary shows such a discretization for the straight line in workspace. Clearly, a nonequidistant sampling of the joint arc length s could be used to derive the non-equidistant sampling in σ . However, developing such a discretization in s would not be beneficial, since the latter sampling could be derived directly for σ .
To evaluate the quality of the optimal solution z * , the maximum velocity curve (MVC) [4] is used. This curve describes the highest possible speed z along the path. The MVC is derived by computing the intersection of joint velocity, acceleration, and torque constraints in the z − z plane and choosing the intersection point with the smallest value of z.

Straight line
As the first example, the straight line near a singularity in the task space is selected, Fig. 2(a). The corresponding joint coordinates computed with an inverse path kinematics, Section 2, are shown in Fig. 3 as seen before. Solving the OPF problem (19) for the sampling schemes discussed earlier with N = 100 samples yields the optimal solutions z * in Fig. 5(a). Both strategies lead to nearly the same evolution z * (σ ). Also, the passing near a singularity at σ ∈ [0.15, 0.25] is reflected, since the MVC and thus the optimal solutions become nearly zero in this area. The non-equidistant sampling, marked with the small orange crosses, again shows accumulated sampling points in this region as discussed earlier.
As a consequence, the resulting trajectories are expected to be smoother compared to the equidistant discretization, as discussed previously in Section 3. Moreover, the constraint violation in an interval of two consecutive sampling points, when using the linear interpolation profile for z is much lower. This is pointed out in Fig. 5 The resulting optimal joint velocities in Fig. 6(a) show thatq 4 is limited by the velocity constraint, when close to the singularity. Due to the vicinity to the wrist singularity (q 5 → 0), see Fig. 3 at σ ≈ 0.2, eitherq 4 orq 6 was expected to fulfill the velocity constraint. In the remaining section of passing near the singularity,q 1 is limited by the constraint. Thus, the evolution ofq 1 andq 4 at this region is viewed in more detail, Fig. 6(b). As mentioned before, the joint velocities are much smoother in this section, with the non-equidistant sampling compared to the equidistant sampling of σ . This, besides lower constraint violation when using the linear interpolation profile in z, inevitable leads to smoother higher derivatives that is, joint acceleration and jerk with the adaptive sampling strategy. The impact of the non-equidistant sampling of σ on terminal and computation time is covered in Section 5.2.3.

Rectangular path with rounded corners
The following EE-path is inspired by the international norm ISO-9283 for industrial robot calibration and performance testing. The norm involves different geometrical paths in the workspace, such as circles, straight lines, and a rectangle with rounded corners. For performance tests, the various paths are passed in a sequence. In this paper, the rectangle with rounded corners was picked. This shape was preferred of two reasons: First, because the curve has a higher geometrical complexity compared to the circle and the straight line (which is already covered in Section 5.2.1); second, because the length of this path is longer compared to the other shapes. The rectangle is placed in front of the robot with an angle of 45 • to the vertical axes, Fig. 7(a). The path is passed in counterclockwise direction, starting at the point shown in the figure. Again the nonlinear map s(σ ) with an equidistant sampling in s is presented in Fig. 7(b), to illustrate the non-equidistant sampling scheme σ i = σ (s i+1 ) − σ (s i ). Compared to the relation from the straight line, the map is approximately piece-wise linear. Thus, the samples of the non-equidistant scheme are pronounced in regions with high slopes and vice versa. Figure 8 shows the corresponding joint paths space using both inverse path kinematics discussed in Section 2. Viewing the evolution of the joint coordinates and the nonlinear map s(σ ), it is evident  that the robot is not passing near a singularity. Furthermore, the lower straight line of the rectangle at σ ∈ [0.15, 0.35] leads to mainly linear joint paths. This is especially pronounced for q 1 , q 4 and q 6 . As a consequence of the linear joint evolution, the number of samples, using σ i = σ (s i+1 ) − σ (s i ), is reduced in this region. Since following simple geometries, for example, linear paths, does not require much sampling points, the latter behavior is even desirable. At σ ∈ [0.65, 0.85], the robot is passing the upper line of the rectangle, again resulting in mainly linear joint coordinates q 1 , q 4 and q 6 , but with a higher gradient as for the lower line. Thus, more samples are present for this linear path. Solving (19) for the rectangle with both sampling schemes yields to the solution illustrated in Fig. 9(a). At the rounded corners of the rectangle, only lower path speeds can be achieved. This is reflected in the MVC at σ ∈ {0.16, 0.38, 0.65, 0.9}. Both solutions again share nearly the same evolution over the curve parameter σ . The highest deviations are present in the section σ ∈ [0.15, 0.35], which corresponds to the lower straight line of the rectangle mentioned before. Due to the fact that the sampling points are not the same σ i = σ (s i ), the constraints in the optimization problem are evaluated at different values of σ . Therefore, constraints become active or inactive at a different progression of the path, effecting the further evolution. In this case, even higher path speeds z * , for the lower straight line, are achieved with the adaptive sampling strategy   similar for both strategies. Viewing the transition between the rounded corners and the straight line at the beginning and end of this section, it is observed that the non-equidistant sampling scheme leads to a smoother optimal path speed z * . Thus, the joint velocities, accelerations, and jerks are expected to be smoother at this transition areas. As mentioned above, the adaptive sampling results in higher values of z at σ ∈ [0.65, 0.85]. Thus, higher joint velocities are achieved as can be seen in Fig. 10(a), particularly forq 1 andq 4 . This results in a faster traversal of the lower straight line with σ i = σ (s i+1 ) − σ (s i ). In Fig. 10(b), a detailed view of passing the upper straight line is depicted for the joint velocitiesq 1 andq 6 . As previously indicated, the transition between the rounded corners and the straight line at σ ≈ 0.65 and σ ≈ 0.85 is smoother for the non-equidistant discretization.

Computation time vs. terminal time
In this section, the terminal time t T and computation time t comp. are compared. The corresponding time values are depicted in Table II, for different number of samples N and four EE-paths. The viewed paths include the straight line and rectangle as discussed above as well as a meander-shaped EE-path and an arbitrary spatial curve. The latter two paths are illustrated in Fig. 11, with their corresponding nonlinear map s(σ ). In addition, the relative difference p t T and p tcomp. compared to the equidistant discretization strategy in σ is also listed in the table. Thus, negative numbers reflect a loss in time (higher values of t T ) and vice versa. First, the results A) of the straight line are viewed. The adaptive sampling has an overall positive impact on the computation time for this path. This is especially pronounced for N = 100 with a reduction of 27%, compared to the equidistant sampling. Viewing the terminal time a loss in time is observed. This can be explained due to the fact that the integral of (19) is approximated numerically in the optimization process. Therefore, the size of a sampling interval σ i noticeably influences the quality of the approximation. Since the number of samples N is fixed for the whole path, a raise of samples in one section reduces the number of samples for the remaining curve in return. For the straight line, most of the samples are accumulated in the region near the singularity, see for example, Fig. 5. Therefore, the sampling intervals σ i of the remaining path are larger on average, for the adaptive sampling strategy. This results in a poor approximation of the time interval in sections not near the kinematic singularity which leads to an increase or decrease of the terminal time t T . In the case of the straight line, the poor approximation induces an increased terminal time, compared to the equidistant sampling scheme. Comparing the terminal and computation time of (B) the rectangle with rounded corners, the adaptive sampling strategy barely impacts both. For (C) the meander-shaped EE-path, Fig. 11(a), the behavior is quite the same as for (A) the straight line. The computation time of the optimal solution is reduced for all three numbers of samples, but on the contrary the terminal time is increased when using the nonequidistant sampling scheme. The strategy σ i = σ (s i+1 ) − σ (s i ) comes with a slight improvement of both the terminal t T and computation time t comp. , for D) the arbitrary curve, Fig. 11(c).

Conclusion
The inverse kinematics problem was expressed in terms of the arc length s in joint space when the EEmotion is parameterized by a path parameter σ . This leads to an ODE system for q(s) and σ (s) which is solved for s ∈ [0, s T ]. It is shown that using the arc length as independent parameter (instead of σ ) allows for coping with kinematic singularities. Moreover, the arc length parameterized velocity inverse path-kinematics solution (10) gives rise to an adaptive sampling scheme.
This derived discretization strategy was used to solve the optimal path following problem (19) with SOCP. The solution was discussed in detail for two exemplary EE-paths. The optimization results, with the non-equidistant sampling, were compared with an equidistant sampling in σ . It is shown that the adaptive sampling strategy leads to smoother trajectories near singularities and geometrically challenging paths (in joint space). The non-equidistant discretization also copes well with the presence of simple geometries, where the sampling rate is reduced, see Section 5.2.2.
Furthermore, the impact on the terminal and computation time was discussed in detail. It was observed that non-equidistant sampling leads to an overall reduction of the computation time. In return, the terminal time shows slightly increased values.
Finally, it should be mentioned that the derived adaptive sampling strategy has to be investigated in detail for every prescribed EE-path since the inverse path kinematics, and thus, the non-equidistant sampling is highly depending on the shape of the EE-path and position in the workspace.