Velocity and acceleration level inverse kinematic calculation alternatives for redundant manipulators

For both non-redundant and redundant systems, the inverse kinematics (IK) calculation is a fundamental step in the control algorithm of fully actuated serial manipulators. The tool-center-point (TCP) position is given and the joint coordinates are determined by the IK. Depending on the task, robotic manipulators can be kinematically redundant. That is when the desired task possesses lower dimensions than the degrees-of-freedom of a redundant manipulator. The IK calculation can be implemented numerically in several alternative ways not only in case of the redundant but also in the non-redundant case. We study the stability properties and the feasibility of a tracking error feedback and a direct tracking error elimination approach of the numerical implementation of IK calculation both on velocity and acceleration levels. The feedback approach expresses the joint position increment stepwise based on the local velocity or acceleration of the desired TCP trajectory and linear feedback terms. In the direct error elimination concept, the increment of the joint position is directly given by the approximate error between the desired and the realized TCP position, by assuming constant TCP velocity or acceleration. We investigate the possibility of the implementation of the direct method on acceleration level. The investigated IK methods are unified in a framework that utilizes the idea of the auxiliary input. Our closed form results and numerical case study examples show the stability properties, benefits and disadvantages of the assessed IK implementations.

redundant, if there is higher number of DoFs than the dimensions of the task.
In practice, the desired position of the end-effector (tool-center-point) is prescribed and we try to find the corresponding joint coordinate values. This mapping from the task space to the joint space is called inverse kinematics (IK). Several methods exist in the literature to solve the IK problem with different approaches. These methods can be categorized as geometric level approaches, velocity level approaches and acceleration level approaches. Acceleration level resolution methods of the IK problem may improve the performance of redundant robots [13,24] comparing to velocity level methods. There are alternatives for the position error compensation too [9,10]. In this paper, our focus is on the stability, accuracy and efficiency of different combination of approaches used for IK.
Non-redundant cases were considered in our recent paper [22]. The present paper generalizes the theoretical details and the set of case examples for redundant cases.

The idea of auxiliary input in motion control
We formalize and classify the presented IK approaches in a rigorous way. To this end, we borrow the notion of the auxiliary input that is used in the feedback linearization of the dynamic equations of robotic manipulators, e.q. when computed-torque like control or inverse dynamics control is applied for trajectory tracking [26,28]. We summarize here the theoretical background briefly based on the literature. Let us assume the equation of motion of the manipulator in the following general form MðqÞ € q þ Cðq; _ qÞ ¼ HðqÞu : Here the vector of general coordinates is q, the mass matrix is MðqÞ, the inertial forces are represented by the vector Cðq; _ qÞ, the independent control inputs are collected in u and the control input distribution matrix is HðqÞ. One obtains the closed loop system in the following linear form if the control input u is synthesized by using the formula The typical choice of the so-called auxiliary inputṽ, which appears in (3), realizes a proportional-derivative (PD) linear feedback control [26]: The PD controller realized in (4) and the proper choice of the gain valuesP andD provide stable error dynamics (i.e. ðq À q d Þ7 !0), when the desired trajectory q d ðtÞ is tracked. The tuning of the gain values is detailed in [26]. One can conclude that the auxiliary inputṽ is chosen with respect to a particular control goal in general. In analogue, we introduce the auxiliary input for IK.

Synthetization of the auxiliary input, goals
In the subsequent sections, the variety of IK algorithms is realized by means of different synthesis of the auxiliary input: the linear error feedback and an alternative approach which uses the idea of direct error elimination. Both concepts are investigated on velocity and acceleration level, resulting four alternatives: velocity level feedback, velocity level direct elimination, acceleration level feedback and acceleration level direct elimination methods (VF, VD, AF and AD respectively). Three out of the four (VF, VD and AF) variations can be found in the literature, and they are only reformulated in the framework of the auxiliary input. Furthermore, their stability analysis is carried out considering digital control (i.e. time-quantization) [12,17]. Digital effects are quite important and the time-step length affects the stability properties. However, a very few studies conduct discrete-time stability analysis, e.g. [12] and continuous-time stability proofs are typical, e.g. [2].
The question naturally arises: why is not there direct error elimination approach on acceleration level? This fourth alternative (AD) has been proposed in our recent paper [22] for non-redundant cases. In the present paper, we generalize AD for redundancy and investigate its feasibility. The stability analysis is carried out in [22] for all VF, VD, AF and AD variations for non-redundant case. This discrete-time stability analysis is generalized for kinematic redundancy in the present paper.

Inverse kinematics calculation alternatives
Transferring the motion to the workspace from the joint space and vice versa is an essential component of motion control algorithms. In IK problems, the desired trajectory r d ðtÞ : R7 !R m of the tool-center-point (TCP) position (or end-effector position) is prescribed in the workspace in the case of position trajectory tracking task [4,14,16,26]. The tracking error e is introduced as with rðqÞ : R n 7 !R m . The term rðqÞ in (5) maps the joint coordinates to the actual TCP position and it involves the geometric properties of the manipulator. The desired position r d ðtÞ is purely time dependent. The IK calculation aims to find the joint variable vector qðtÞ 2 R n that satisfies the error equation Note that we do not need to restrict the IK calculations to position tracking. In general case, (5) may interpret any virtual geometric constraint equations, which are referred as servo-constraints or control-constraints in the literature [6-8, 18, 30, 32]. Essentially, we distinguish the kinematically not redundant case, when n ¼ m, and the redundant case, when n [ m. In the case of redundancy there is no unique IK solution; rather, motion optimization is carried out [3,4,14,16,20,26,29,31].
For the geometric level approach, (6) is solved directly as a nonlinear algebraic problem, e.g. [25,27]. Yet, geometric solution is not typical in practice due to the possibly high computational demand. Typically, the IK is solved on velocity [1,20,29] or acceleration level [13,20,24]. The numerical methods presented in this paper are realized in velocity and acceleration level.

Auxiliary input synthesis on velocity level
The velocity level implementation has the advantage that the unknown joint speed _ q appears in linear form. Hence, it's calculation demands low computational capacity. After differentiating (6) with respect to time, is obtained. Equation (7) must be satisfied by determining the solution for _ q. The linear connection of the TCP velocity _ r and the joint speed is given by where the Jacobian matrix J 2 R mÂn is obtained as JðqÞ ¼ r q r. The commanded joint speed is expressed from (9) as where J y denotes the Moore-Penrose generalized inverse [5,19,20,23] of the non-square Jacobian. The term v v in (10) is the auxiliary input, i.e. it has the same role asṽ in (3). Considering Eq. (8), seems to be a natural choice for the auxiliary input.
The problem is that this choice does not provide position error elimination e, although lim t7 !1 is a requirement of IK methods. Therefore a variety of numerical implementations is established in the literature. Providing the possibility of the elimination of position errors, two alternatives for the choice of v v are detailed and studied in the following sections. Note that (10) might be extended depending on the control goal [3,4,14,16,20,26,29,31]. On one hand, a scalar cost function u ¼ _ q | W _ q is minimized when the generalized inverse incorporates the weight matrix W: We assume W to be identity I, which results the simplified formula J y ¼ J | ðJ J | Þ À1 . On the other hand, the null-space N J ¼ I À J y J of the Jacobian is often used in order to minimize a scalar cost function cðqÞ. Then the joint speed is obtained as _ q ¼ J y v v þ N J _ n with n ¼ r q c. In this paper, cðqÞ ¼ 0 is assumed.

Velocity level linear feedback approach (VF)
The present IK calculation alternatives differ in the synthesis of the auxiliary input v v that appears in (10). The TCP position error is eliminated by using v v :¼ _ r d À jðr À r d Þ, [29] with which we obtain the commanded joint speed as By substituting (14) into (9), the TCP velocity arises in the form Equation (15) shows that the error dynamics is stable for the scalar gain parameter j [ 0. By neglecting the digital effects [17], the following first order ordinary differential equation describes the error dynamics of with the exponentially decaying solution is eðtÞ ¼ Be Àjt . We detail the stability analysis and the proper choice of j in Sect. 3.1 in the presence of timequantization.

Velocity level direct error elimination approach (VD)
Comparing to the VF method, the auxiliary input v v (see: (10)) is synthesized in an alternative way. The fundamental goal is the direct elimination of the tracking error by a proper choice of the v v [9,10].
The actual values of the joint variables and the TCP position in the current time instance t i are introduced as: These variables and the time-step h :¼ t iþ1 À t i are shown in Fig. 1 left panel. Let us take the difference d i ¼ r d iþ1 À r i of the actual and the desired TCP position respectively in the current time instant t i and upcoming time instant t iþ1 . The goal is to express the joint position increment directly with d i . Since the time-step h\\1 is small, constant TCP velocity is assumed _ r i . We assume that considering the yet unknown velocity _ r i , the upcoming desired TCP position r d iþ1 is reached: The TCP velocity _ r i is expressed from (17). Indeed, the error Substituting the auxiliary input into (10), the commanded joint speed is obtained as The commanded joint position in the upcoming time is directly expressed by applying the single explicit Euler [11] The TCP position error is shown in the right panel of Fig. 1. The TCP positioning error e iþ1 occures, because the relation of the joint coordinates and the TCP position is nonlinear. However, if rðqÞ is linear, then e iþ1 ¼ 0 is achieved. In general, the error e iþ1 converges to zero for the time-step h7 !0. The convergence analysis is detailed in Sect. 3.1.
The advantage of the direct error elimination method over the linear feedback approach is that the joint position increment (see: (19)) is obtained directly, and there is no need for considering the joint or TCP speed. Furthermore, the joint position increment is directly related to the TCP positioning error.

Auxiliary input synthesis on acceleration level
The IK problem is maintained on acceleration level in some cases, e.g. [13,24]. We gain the acceleration level error equation after time differentiating (6) twice: Equation (21) needs to be satisfied in the acceleration level approaches. Similarly to (9), we express the TCP acceleration Since the joint acceleration is linear in (22), the commanded joint acceleration is expressed as: where v a is again the auxiliary input of the IK control. Considering Eq. (21), one could synthesize the auxiliary input as The stable position error dynamics is not guaranteed by (24). To avoid this problem, a variety of the synthesis of v a is applied in the literature. Similarly to the velocity level approaches, two variations for the choice of v a are detailed in the following sections, which are capable of the positioning error elimination.

Acceleration level linear feedback approach (AF)
By chosing the proper gain parameters j P and j D , the auxiliary input v a :¼ € r d À j D ð _ r À _ r d Þ À j P ðr À r d Þ provides the elimination of the TCP tracking error [13]. The auxiliary input is substituted in Eq. (23) and we obtain the commanded joint acceleration as By substituting (25) into (22), the TCP acceleration can be expressed as: which leads to a stable error dynamics with the scalar gain parameters j P [ 0 and j D [ 0 governed by the second order differential equation The error dynamics is asymptotically stable, if j P [ 0 and j D [ 0 and we assume continuous dynamics. However, the digital effects [12,17] are important. The stability analysis and the choice of j P and j D are demonstrated in Sect. 3.1 in the presence of the discrete dynamical effects due to digital control.

Acceleration level direct error elimination approach (AD)
Although the description of the IK methods in Sects. 2.1.1, 2.1.2 and 2.2.1 can be found the literature, they typically appear without using the framework of the auxiliary inputs. The method explained in the present Section was proposed first in [22] and here it is generalized for redundant cases. The goal is the direct elimination of the position error by choosing the auxiliary input appearing in (23) properly. The method explained here is the acceleration level analogy of the velocity level direct error elimination (VD) method. The deviation d i ¼ r d iþ1 À r i between the realized and the prescribed TCP position is considered, in analog to the VD approach, as it is shown in Fig. 2. By the proper choice of the constant TCP acceleration € r i , the TCP position prescribed in the upcoming time instant r d iþ1 is reached: The error d i is approximately eliminated if the auxiliary input is chosen as v a;i ¼ 2ðr d iþ1 À r i Þ=h 2 À 2 _ r i =h. This formula is obtained by the solution (28) for the unknown acceleration € r i . By substituting v a;i into (23) and applying (9) for the actual time instant ( _ r i ¼ J i _ q i ), and introducing the joint speed approximation _ q i % J y ðr d iþ1 À r i Þ=h, the commanded joint acceleration Fig. 2 Approximation of the constant TCP acceleration can be synthesized. In order to obtain explicit formula for the joint position, one-step time integration is carried out. First, the velocity and after that the joint position is expressed with the Adams-Moulton family of numerical integrators [11]: We directly obtain the upcoming commanded joint position. Again, small TCP position error e iþ1 arises because of the nonlinearity of rðqÞ; however, e iþ1 converges to zero when h7 !0. The results of the convergence analysis are presented in Sect. 3.1.
In analogy with the VD approach, the increment of the joint position vector is obtained without including the joint or TCP acceleration in the formulae.

Convergence analysis
We couldn't find the rigorous stability analysis and comparison of the presented four IK methods in the presence of digital effects. The stability analysis of the four methods (VF, VD, AF and AD) is carried out by the time integration (if necessary) and linearization (for nonlinear case examples) of the commanded joint position/speed/acceleration formulae (14), (19), (25) and (29). The linearization is carried out around an arbitrarily chosen configuration q 0 . In each case, we obtain the discrete mappings [15] that map the commanded joint variables from the time instant t i to t iþ1 . In the case of the VD approach, the joint positions are mapped only. Although the mapped variables x ¼ ½q | ; _ q | | include the joint speed too in the case of the VF and AD methods. For the AF approach, the joint acceleration is also included in z ¼ ½q | ; _ q | ; € q | | . The coefficient matrix A and the vector b change depending on the formulae for the auxiliary input v v or v a and on the integration scheme. The eigen-values of A define the convergence: jk i j\1; 8i guarantees asymptotic stability [15], defined in (12). The eigen-value analysis of each particular method is summarized in Sect. 3.1.5. It was shown in our recent paper [22] that the matrix A for models B and C are identical. It was also shown that the matrix A for models B and C have the same pattern as for model A, but the pattern is repeated twice. It is because the workspace is two-dimensional instead of one. As a consequence, the mapping matrix is shown only for the models A and D with onedimensional workspace.

Convergence analysis of the VF approach
In the VF approach, the numerical integration scheme was the second order Adams-Bashforth scheme [11]: This integration scheme is applied on the commanded joint velocity defined in (14), with which the mapping matrix is obtained in the form below for the model A and configuration q 0 ¼ 0: Note that the matrix A VF repeats the same pattern for models B and C. Indeed, the eigen-values are also the same, but multiple. The real eigen-values of A VF are In redundant case (model D) with q 0 ¼ 0, the mapping matrix is obtained as: Àj 0 0 guarantees the decay of the tracking error e. Note that the stable region has infinitely long but monotonically narrowing tails both in the direction of the time-step duration h and in the direcion of the gain parameter j.

Convergence analysis of the VD approach
The numerical integration is already incorporated into the commanded joint position formula defined in (19).
The mapping matrix for model A with q 0 ¼ 0 reads: of which the eigen-value is always zero, indeed. For models B and C with two-dimensional workspace, the matrix is 2 by 2 zero matrix, i.e. the pattern is repeating itself again. For the redundant case (model D) with q 0 ¼ 0, the mapping matrix is: The redundancy introduces one extra row and column for each extra DoF, and as a consequence, a new spurious eigen-value 1 appears. The exponential convergence of the tracking error to zero is always guaranteed.

Convergence analysis of the AF approach
Model A is considered again with q 0 ¼ 0. In case of the acceleration level feedback IK approach, the integration scheme was the second order Adams- Bashforth scheme combined with the second order Adams-Moulton scheme [11]: This integration scheme is applied on the commanded joint acceleration defined in (25) with which, the mapping matrix reads The real eigen-value k a1 ðj P ; j D ; hÞ and the complex eigen-values k a2 ðj P ; j D ; hÞ, k a3 ðj P ; j D ; hÞ of A AF can be obtained in closed form; however, we do not show the resulting formula because of its high complexity.
In the redundant case, when model D is investigated, the mapping matrix arises in the form: Besides the eigen-values k a1 , k a2 and k a3 , additional three eigen-values 0, 1, 1 appear for each extra DoF. The stable region in the parameter space is shown in the right panel of Fig. 4. The choice of parameters inside the stable region (yellow) guarantees the decay of the tracking error e. There is no stable region for the gain value j D ¼ 0; although increasing j D modifies the shape of the stable region. The most important difference comparing to the stability chart of the VF method is that the maximum time interval h is limited.
A very important observation about the linear feedback approaches is that Fig. 4 clearly shows that the stable region for the time-step h and the proportional gains j and j P is much larger for the acceleration level approach than for the velocity level one.
The practically important benefit is that larger timesteps are allowed without loosing stability.

Convergence analysis of the AD approach
Model A and q 0 ¼ 0 is investigated first. The commanded acceleration is defined by (29), and the numerical integration scheme is defined in (30) and (31), with which the mapping matrix is: of which the eigen-values are 0 and À1 for any timestep h. This two-element eigen-value set is multiple for each dimension of the workspace. For the redundant case (model D with q 0 ¼ 0), the mapping matrix is of which the eigen-value set extends with two spurious 1 eigen-value for each extra DoF. The new spurious eigen-values do not change the fact that the zero error solution is marginally stable because of the À1 eigenvalue. In practice À1 eigen-value causes a nondecaying timestep-by-timestep numerical oscillation.

Summary on the convergence analysis
For all the A, B, C, D, E, and F case examples, the same stability chart is valid for the velocity and acceleration level linear feedback approaches (VF and AF) (see: Fig. 4). Two parameters, namely h and j determine the convergence for the VF method, while the convergence of the AF method is determined by h, j P and j D . The convergence of method VD is guaranteed, since the non-spuriois eigen-values are always zeros. Method AD is proven to be marginally stable, since there is always an eigen-value, which is À1. The eigen-value sets are summarized for each IK method in Table 1. In redundant case, the number of eigen-values is 2m, m, 3m and 2m respectively for the VF, VD, AF and AD methods. The number of eigenvalues introduced by the redundant extra DoFs is 2ðn À mÞ, n À m, 3ðn À mÞ, 2ðn À mÞ respectively.
The case examples in the subsequent section show that the eigen-value analysis and the numerical simulations are in correspondence. In each case, both stable parameter sets and parameter sets from the vicinity of the stability border were tested.

Case examples A and D
Case example models A and D (see Fig. 3) are the simplest possible ones in the sense that the workspace is m ¼ 1 dimensional. If there is no redundancy, n ¼ 1. However, for redundant case n [ m must be true, which yields n ¼ 2 DoF in the simplest case. We consider a nonlinear mapping from the R 1 or R 2 joint space to the R 1 workspace without any particular physical meaning.

Kinematically non-redundant case: model A
The nonlinear function r(q) and the desired TCP position are given. The desired TCP trajectory r d ðtÞ and the function r(q), which maps from the joint space to the workspace, were r d ðtÞ ¼a 0 þ a 1 sin xt ; ð48Þ in the numerical simulations and also in the analytical convergence analysis. The parameter values were a 0 ¼ 1 m, a 1 ¼ 1 m, x ¼ 1 rad=s and c ¼ 0:2 m. The simulations were performed with an initial error as it is shown in the top panels of Fig 5. The gain parameters and the time-step were j ¼ 19 1=s, j P ¼ 250 1=s 2 , j D ¼ 8 1=s and h ¼ 0:05 s respectively for the simulations that are carried out in the inner vicinity of the stability border shown in Fig. 4 (close to the marginal stability). A parameter set deep inside the stable region was j ¼ 5 1=s and j P ¼ 50 1=s without changing any other parameter. Note that every other simulation is carried out with the same gain parameters and timestep. The initial conditions are q 0 ¼ 1, _ q 0 ¼ 0.

Kinematically redundant case: model D
The joint variables are mapped to the workspace by Every other parameter is the same as in the nonredundant case example A. The initial conditions were q 0 ¼ ½1; 1 | and _ q 0 ¼ 0. The simulation results are shown in the bottom panels of Fig 5. As we can see in the graphs for VF, and AF methods, the almost unstable parameter sets provide stability; however the error decays very slowly with high oscillations. The convergence speeds up when choosing the parameters from deep inside the stable range. Method VD suppresses the error quickly, but method AD causes abrupt peaks in the position error. We expected this kind of stability problem based on the eigen-value analysis.

Case examples B and E
The non-redundant and the redundant cases were tested on the linear models B and E. Models B and E are a bit more general than A and D, since the workspace is not one-dimensional. The R 2 7 !R 2 and R 3 7 !R 2 linear mapping functions for the planar PP and PPP manipulators, which consist of two and three perpendicular prismatic drive (see: Fig. 3), are respectively Table 1 Eigen-values of the discrete mapping for each IK approach Eigen-values in non-redundant case Extra eigen-values in redundant case VF m sets of k v1 , k v2 n À m sets of 0, 1 VD m sets of 0 n À m sets of 1 AF m sets of k a1 , k a2 , k a3 n À m sets of 0, 1, 1 AD m sets of 0, À1 n À m sets of 1, 1 The eigen-values for the non-redundant cases are collected in the middle row. The number of eigen-value sets is equal to the n ¼ m dimensions of the workspace. The rightmost column shows the extra spurious eigen-values which are additionally introduced for the extra n À m number DoFs r E ðqÞ ¼ xðqÞ where l ¼ 1 m is a geometric parameter. The desired path was a circle defined by with parameters a 0 ¼ 1 m, a 1 ¼ 0:5 m and x ¼ 2 rad=s. For model B, the initial conditions were q 0 ¼ ½0:5; 0:5 | and _ q 0 ¼ 0. For model E, the initial conditions were set to q 0 ¼ ½0; 0; 0 | and _ q 0 ¼ 0. First, we carried out numerical simulations with parameters intentionally chosen from the border of the stable range; i.e. the norm of the eigen-values of A is about 1 (for the IK parameter values see: Sect. 3.2.1). The Euclidean norm of the TCP positioning error is shown in Fig. 6. The behaviour of the IK methods is similar to the previous case examples A and D. However, models B and E are linear. Because of the linearity, methods VD and AD are able to track the given trajectory without error. Fig. 7 depicts the desired and the realized path of the TCP. The vibrations in the vicinity of the stability border are clear for VF and AF (left column). However stable parameters provide smooth decaying of the tracking error (middle column). Although the direct methods practically jump to the desired path, in the physical reality, the motion controller and the actuators limit the speed of this motion.

Case examples C and F
Nonlinear non-redundant and redundant models (C and F respectively) are taken as case examples. Models C and F are planar RR and RRR manipulators. The joint variables are mapped to the TCP position by the following R 2 7 !R 2 and R 3 7 !R 2 functions r C ðqÞ ¼ xðqÞ yðqÞ r F ðqÞ ¼ xðqÞ yðqÞ where l ¼ 1 m is the length of the bars and the abbreviations are c i ¼ cosðq i Þ, s i ¼ sinðq i Þ. For model C, the initial conditions were q 0 ¼ ½22:5; 81 | and _ q 0 ¼ 0 for methods VF, VD, AF and q 0 ¼ ½30:5; 81 | and _ q 0 ¼ 0 for method AD. For model F, the initial conditions were q 0 ¼ ½À10; 48; 132 | and _ q 0 ¼ 0 for methods VF, VD, AF and q 0 ¼ ½À1; 48; 132 | for method AD. Method AD is very sensitive to any large initial error, that is the reason for choosing different initial conditions.
Simulations with close-to-unstable and stable parameters were carried out. For the parameters j, j P , In the case of these nonlinear example models, the behavior is a bit more intricate then for cases B and E. As we can see in Fig. 8, there is a small positioning error remaining for IK methods VF and AF because of the nonlinearity. Again, method AD develops oscillations after settling down (right, bottom panel). This phenomena questions the practical applicability. Fig. 9 shows the desired and the realized trajectories.

Conclusion
Four alternatives for the numerical implementation of the inverse kinematics calculation were studied: the feedback approach and the direct error elimination approach were implemented on velocity and acceleration level in the framework of auxiliary inputs. The comparison of the stability properties were carried out and the stable parameter regions were presented, which are essential in the practical application of these Fig. 6 The Euclidean error norm e k k ¼ rðqÞ À r d ðtÞ for the case example models B (top row) and E (bottom row). Almost unstable and stable cases for VF and AF approaches and the results for VD and AD methods are shown in separate columns Fig. 7 The desired r d ðtÞ and the realized trajectory rðqÞ of the TCP for the case example models B (top row) and E (bottom row). Close-to-unstable and stable cases for VF and AF approaches and the results for VD and AD methods are shown in separate columns. The initial positions are shown by larger markers algorithms. The acceleration level feedback method is proven to be beneficial over the velocity level type in the sense that the stable parameter region is larger and therefore larger time-step is allowed in the digital implementation. In order to test the possible benefits of acceleration level methods, we introduced the acceleration level direct error elimination method (AD). The analytical convergence analysis showed that VF, VD, AF methods are convergent and therefore feasible in practice. The stable parameter range was defined for VF and AF methods. It was however proven that the AD method is marginally stable, which means that its practical application is not feasible in the present form. Some research is worth being carried out on the possibilities of stabilization, because the convergence of the AD method is much faster than the AF feedback method when the underlying system is linear. The analytically obtained stability charts were tested by numerical case example simulations, and the behavior was as one could expect based on the stability charts. Fig. 8 The Euclidean error norm e k k ¼ rðqÞ À r d ðtÞ for the case example models C (top row) and F (bottom row). Almost unstable and stable cases for VF and AF approaches and the results for VD and AD methods are shown in separate columns Fig. 9 The desired r d ðtÞ and the realized trajectory rðqÞ of the TCP for the case example models C (top row) and F (bottom row). Close-to-unstable and stable cases for VF and AF approaches and the results for VD and AD methods are shown in separate columns. The initial positions are shown by larger markers