Approximation of Jacobian Inverse Kinematics Algorithms: Differential Geometric vs. Variational Approach

This paper addresses the approximation problem of Jacobian inverse kinematics algorithms for redundant robotic manipulators. Specifically, we focus on the approximation of the Jacobian pseudo inverse by the extended Jacobian algorithm. The algorithms are defined as certain dynamic systems driven by the task space error, and identified with vector field distributions. The distribution corresponding to the Jacobian pseudo inverse is non-integrable, while that associated with the extended Jacobian is integrable. Two methods of devising the approximating extended Jacobian algorithm are examined. The first method is referred to as differential geometric, and relies on the approximation of a non-integrable distribution (in fact: a codistribution) by an integrable one. As an alternative, the approximation problem has been formulated as the minimization of an approximation error functional, and solved using the methods of the calculus of variations. Performance of the obtained extended Jacobian inverse kinematics algorithms has been compared by means of computer simulations involving the kinematics model of the 7 dof industrial manipulator POLYCRANK. It is concluded that the differential geometric method offers a rapid, while the variational method a systematic tool for solving inverse kinematic problems.

tion) by an integrable one. As an alternative, the approximation problem has been formulated as the minimization of an approximation error functional, and solved using the methods of the calculus of variations. Performance of the obtained extended Jacobian inverse kinematics algorithms has been compared by means of computer simulations involving the kinematics model of the 7 dof industrial manipulator POLYCRANK. It is concluded that the differential geometric method offers a rapid, while the variational method a systematic tool for solving inverse kinematic problems.

Introduction
The forward kinematics of a robotic manipulator define the location (i.e. position and orientation) of the end effector as a function of the position of the manipulator's joints. The inversion of the kinematics, i.e. computing the joint position corresponding to a prescribed end effector's location is referred to as the inverse kinematic problem. For redundant serial manipulators, the inverse kinematic problem has an infinite (continuum) number of solutions. In order to provide a unique solution, the inverse kinematics algorithms impose on the solution some additional constraints, like the minimization of the instantaneous joint velocities or other objective functions, avoiding obstacles, joint limits or joint singularities, etc. The state of the art in resolving kinematic redundancy of robotic manipulators has been recently reviewed in [1]. Most frequently, in the case of redundant manipulators the inverse kinematic problem is solved numerically, using Jacobian inverse kinematics algorithms, e.g. the Jacobian pseudo inverse or the extended Jacobian algorithm. By design, the former algorithm minimizes the instantaneous joint velocity; it also distinguishes by quick convergence. The latter algorithm is designed by completion of the manipulator's Jacobian with the derivative of an augmenting kinematics function [2]. Because the extended Jacobian algorithm drives the joints along a level manifold of the augmenting function, this algorithm can solve the inverse problem and simultaneously keep an objective function at its minimum [3]. For the same reason, the algorithm has the property of repeatability, which means that closed paths travelled by the end effector in the task space are converted into closed paths in the joint space [4]. Obviously, repeatability considerably facilitates the control of a manipulator accomplishing cyclic tasks, for after each cycle the joints return to the home position. It is well known that the Jacobian pseudo inverse lacks of repeatability [5]. An alternative to the Jacobian algorithms is the Lagrange multiplier based algorithm [6]. A comparison of the extended Jacobian and the Lagrange multipliers approaches can be found in [7]. The tracking problem of algorithmic singularities in extended Jacobian algorithms has been addressed in [3] and refined in [7].
Within the domain of synthesis of Jacobian inverse kinematics algorithms the idea of shaping the algorithm's performance by combining the advantages of diverse algorithms has been fostered by Roberts and Maciejewski in a series of papers [8][9][10]. These papers addressed the problem of optimal approximation of the Jacobian pseudo inverse by a repeatable Jacobian algorithm. The optimality was defined in terms of an approximation error functional whose minimization was achieved using the tools of the calculus of variations. The resulting necessary optimality conditions assumed the form of a nonlinear partial differential Euler-Lagrange equation. The error functional corresponding to the approximation of the Jacobian pseudo inverse by the extended Jacobian algorithm has been re-defined in [11], in such a way that the Euler-Lagrange equation determining the augmenting kinematics function results in a collection of linear elliptic partial differential equations. Using a geometric interpretation of the inverse Jacobian algorithm as a vector field distribution, an alternative differential geometric approach has been proposed in [12]. This approach is based on the approximation of a non-integrable codistribution by an integrable one, developed for the purpose of approximate feedback linearization of control systems in [13]. Examples of the approximate extended Jacobian inverse kinematics algorithms designed within the differential geometric and the variational approaches have been presented in [14]. Specifically, the variational approach has been applied to academic examples of 3 dof manipulators, for which the corresponding Euler-Lagrange equation can be solved either symbolically or using a standard numeric procedure provided by the MATLAB PDE toolbox. In [14] the differential geometric approach has been used to the simplest nonholonomic system: the unicycle, while in [12] to the POLYCRANK manipulator, also dealt with below.
In an attempt at a further advancement and an assessment of the ideas exposed in the last two references, this paper presents a comparative study of the differential geometric and the variational approaches to the design of the extended Jacobian inverse kinematics algorithms approximating the Jacobian pseudo inverse algorithm, using as a testbed the 7 dof industrial manipulator POLY-CRANK. We begin with a derivation of these algorithms using a continuation method argument that defines the algorithm as a dynamic system driven by the task space error, whose trajectory converges to a solution of the inverse problem. The vector fields constituting this dynamic system span a distribution associated with the algorithm. When the kinematic redundancy is small, it is convenient to replace the distribution by a dual object: a codistribution spanned by differential 1-forms annihilating the associated distribution. In the case of the Jacobian pseudo inverse this codistribution is spanned by 1-forms annihilating the Jacobian transpose (row vectors spanning the Jacobian null space), while for the extended Jacobian algorithm the associated codistribution is spanned by the differentials of components of the augmenting kinematics function. Both the approaches require regularity of the associated distributions/codistribution, so they are applicable outside kinematic and algorithmic singularities. The approximation problem consists in designing a distribution associated with the extended Jacobian inverse that would be in some sense close to the given distribution associated with the Jacobian pseudo inverse. In accordance with the differential geometric approach the singularityfree region of the joint space is converted into a foliation whose leaves, of dimension equal to the redundancy degree of the kinematics, are homotopic to a basic leaf. The homotopy map defines in the joint space a homotopy vector field. The associated codistribution is required to coincide with the given one on the leaves of the foliation and along the homotopy vector field. This requirement translates into a partial differential equation for the augmenting kinematics function that is solved by means of the method of characteristics. Finally, the augmenting kinematics function is computed by numerically integrating backward in time the characteristic (ordinary differential) equation. A central point of the variational approach is the approximation error functional that measures the distance between these two associated distributions/codistributions over a regular region of the joint space. The corresponding Euler-Lagrange necessary optimality condition takes the form of a linear elliptic partial differential equation. For the reason that in the Euler-Lagrange equation the number of independent variables equals the number of dof of the manipulator, this equation cannot be solved using the standard PDE solvers. Instead, the Ritz method of solution has been used allowing to solve the problem efficiently. In order to make the comparison of these two approaches plausible, the definition of approximation error functional has been adapted to the representation of the algorithms by the associated distributions. Using these two approaches the extended Jacobian inverse kinematics algorithms have been designed for the 7 dof POLYCRANK manipulator [15,16]. Their performances have been examined by computer simulations. Thus, the contribution of this paper is two-fold: first, we give a new formula for the approximation error functional of the Jacobian pseudo inverse by the extended Jacobian algorithm, in terms of the associated codistributions (see Section 4), second, we demonstrate that both these approaches are applicable to realistic industrial manipulators, and compare the performance of the corresponding extended Jacobian inverse kinematics algorithms.
This paper is organized in the following way. Section 2 introduces the basic concepts. The differential geometric approach is presented in Section 3. Section 4 focuses on the minimization of the approximation error by means of the variational approach. A comparison of the two approaches, including some computational aspects, is made in Section 5, based on example inverse kinematic problems for the POLYCRANK manipulator. Section 6 sketches a generalization of the presented approach towards higher redundancy of the kinematics. Section 7 concludes the paper.

Basic Concepts
We shall study a coordinate representation of the kinematics of a robotic manipulator, with n degrees of freedom and m-dimensional task space. The number r = n − m will be called the degree of redundancy of the kinematics. Let J(x) = Dk(x) denote the manipulator's Jacobian. Given the kinematics (1) and a desirable point y d in the task space, the joint position x d such that k(x d ) = y d is obtained by solving the inverse kinematic problem. This problem is usually solved numerically, by means of a Jacobian inverse kinematics algorithm.
A derivation of these Jacobian algorithms can be accomplished using the following reasoning borrowed from the continuation method [17]. For an initial configuration x 0 ∈ R n , we define a joint space curve x(θ ), θ ∈ R, passing for θ = 0 through x 0 , such that the corresponding task space error e(θ ) = k(x(θ )) − y d decreases exponentially with a prescribed rate γ > 0, so that By the differentiation of the error we obtain an implicit differential equation sometimes referred to as the Ważewski-Davidenko equation. After applying a right inverse J # (x) of the Jacobian, this equation converts to a dynamic system whose trajectory approaches in the limit a solution of the inverse kinematic problem, x d = lim θ→+∞ x(θ ). In this way the dynamic system (2) actually defines an inverse kinematics algorithm. Observe that, by design, any algorithm based on a certain right inverse of the Jacobian will have the same error convergence rate γ . In computations a discrete version of system (2) is employed, i.e.
where θ = 0, 1, . . .. In this paper we shall focus on two right inverses of the Jacobian, namely the Jacobian pseudo inverse and the extended Jacobian inverse. It is well known that at the regular configurations of the manipulator the Jacobian pseudo inverse is defined as The extended Jacobian inverse can be introduced in the following way. We choose an augmenting kinematics map and define the extended kinematics transforming the joint space into itself. Wherever the extended kinematics map is a local diffeomorphism, i.e. rank(J(x) = Dl(x)) = n, we define the extended Jacobian inverse where mf c means "m first columns", by taking the first m columns of the inverse extended Jacobian. By definition, the extended Jacobian is a right Jacobian inverse, that has the annihilation property It is well known that each of these two algorithms has its specific advantages: the Jacobian pseudo inverse algorithms distinguishes by minimizing the joint velocities and fast convergence, while the extended Jacobian algorithm is repeatable, i.e. transforms closed curves in the task space into closed curves in the joint space. In order to design a Jacobian algorithm endowed with properties of both these algorithms, the idea of approximation of the Jacobian pseudo inverse by the extended Jacobian has been invented. Intuitively, it is expected that the manipulator's joint trajectories generated by two inverse kinematics algorithms close to each other will remain close.
In the classical formulation the approximation problem is an optimization problem that can be solved by means of the methods of variational calculus. In a more geometric setting, since on account of Eq. 2 the columns of the Jacobian inverses (4) and (7) can be regarded as vector fields in R n , one associates with these algorithms a pair of distributions At each configuration x ∈ R n these distributions define an m-dimensional linear space of admissible directions of motion. A distribution is referred to as integrable, if through any point x there passes an m-dimensional manifold whose tangent space is defined by the distribution. By the annihilation property, the distribution D E is integrable, while in general the distribution D P is not. A distribution is integrable, if it is involutive, i.e. the Lie bracket of any pair of vector fields from the distribution belongs to it. The integrability of the associated distribution yields directly the repeatability of the Jacobian algorithm.
The differential 1-forms annihilating these distributions define a pair of dual objects called codistributions, where . . , r, and j = 1, . . . , m. By definition of the distribution and Eq. 4, it follows that D ⊥ P is spanned by 1-forms annihilating the Jacobian transpose J T (x). The descriptions of a Jacobian algorithm by means of a distribution or a codistribution are equivalent, however, especially for a low degree of redundancy r, using a codistribution may appear more convenient. For the manipulator used as a testbed of the derived approximate algorithms has 7 dof, in the following developments we shall always assume that the redundancy degree r = 1.

Differential Geometric Approach
The differential geometric approach to the approximation problem of the Jacobian algorithms relies on designing an integrable codistribution that coincides with the given codistribution in a certain region of the joint space. More specifically, we need to find a codistribution D ⊥ E , associated with the extended Jacobian inverse, that approximates the codistribution D ⊥ P corresponding to the Jacobian pseudo inverse. In accordance with the assumption made at the end of the previous section, the presentation of the differential geometric approach will be adopted to r = 1; a general theory has been exposed in [12]. To begin with, we define in R n a foliation with 1-dimensional leaves E α ∼ = R, parametrized by α = (α 1 , . . . , α m ) ∈ R m , with a distinguished zero-leaf E 0 . Then, we associate with this foliation a homotopy map t : R n → R n , where t ∈ [0, 1], such that 1 = id R n , Using this definition we deduce that On the other hand what means that i.e. the homotopy map t (x) may be regarded as a flow of the time-dependent vector field X(x)/t. To proceed further, it is advantageous to introduce into R n specific coordinates x = (y, z) such that y ∈ R m is constant on every leaf, and z ∈ R varies along the leaf. In these coordinates the zero-leaf E 0 = {(0, z)|z ∈ R}, so 0 (y, z) = (0, z). Now we choose a 1-form ω annihilating J T (x) such that ω| Eα = dz, so the codistribution D ⊥ P = span C ∞ (R n ) {ω(x)}. Denoting the augmenting kinematics map by h(x), we get the codistribution In this setting the approximation problem of the Jacobian pseudo inverse by an extended Jacobian inverse can be given the following formulation: find an augmenting kinematics map h(x) satisfying the condition h(0, z) = z, such that the codistributions D ⊥ P and D ⊥ E coincide on the leaves of the foliation {E α } and along the vector field X(x). Taking into account the form of these codistributions this means that for some smooth function b (x) = 0 and ωX = bdhX.
To proceed further, let us recall that two 1-forms ω,ω coincide on a submanifold M, if ω(x)v = ω(x)v, for any x ∈ M and any vector v ∈ T x M tangent to M at the point x. Now, for x = (y, z) we have Next, setting F(x) = ω(x)X(x) and substituting x = t (y, z), we derive from Eq. 14 the identity y, z)).
Taking into account Eq. 12, this identity is transformed further to and finally, since it follows that the condition (14) is equivalent to With the notation H(t, y, z) = h( t (y, z)) and using Eq. 15, the condition (16) can be expressed in the form of a partial differential equation parametrized by the y-coordinates. The Eq. 17 can be solved by means of the method of characteristics. Indeed, it is easily seen that H(t, y, z(t)) = const along the solution z(t) of the characteristic equation initialized at z(0) = z 0 . Using the identity H(1, y, z(1)) = H(0, y, z(0)) and the properties of the homotopy map 1 (y, z) = (y, z), 0 (y, z) = (0, z), we compute
Observe that the conditions (13) and (14) imply that the codistributions D ⊥ P and D ⊥ E coincide on a leaf E α and along the trajectories t (x) of the vector field (11) joining this leaf with the zero-leaf E 0 . So obtained subsets of the joint space are called pages. It follows that the 1-forms ω and dh coincide on the pages.

Variational Approach
The differential geometric approach has provided an approximate extended Jacobian inverse whose associated codistribution agrees with that of the Jacobian pseudo inverse on pages. Outside the pages the relationship between these two inverse Jacobians is not controlled. An alternative formulation of the approximation problem involves a suitably defined approximation error. By the minimization of this error over a certain region in the joint space we can design an extended Jacobian inverse that is as close as possible to the Jacobian pseudo inverse over this region. A formula for the approximation error will be derived taking into account the codistributions representing the Jacobian inverses.
Assuming that the redundancy degree r = 1, we can represent the inverse by the differential dh(x) of the augmenting kinematics map. This being so, the approximation error can be defined in the following way (ω 1 (x), . . . , ω n (x)) annihilates the Jacobian transpose and restricts to dz on the leaves of the foliation {E α }, and D denotes a subset of R n over which the approximation is made. The error functional (21) should be minimized with respect to augmenting maps, so the approximation problem becomes a problem of the variational calculus. It is easily checked that the integrand in Eq. 21 hence the corresponding necessary optimality condition takes the form of the Euler-Lagrange equation where denotes the Laplace operator, = n i=1 ∂ω(x) ∂ x i . The partial differential equation (22) is known as the Poisson equation; in the case when its right hand side equals zero, the optimal augmenting kinematics map appears to be a harmonic function. For realistic robotic manipulators the Poisson equation cannot be solved, even numerically, therefore a direct minimization method of the error functional is preferred [18]. Applying the Ritz method we assume that h( (ϕ 1 (x), . . . , ϕ p (x)) T , for some basic functions {ϕ j (x)}, and f ∈ R p being a vector of parameters. After a substitution into Eq. 21, we get a quadratic error form for Obviously, the error (23) reaches a global minimum at f * = M −1 N.

Comparison
The approximate extended Jacobian algorithms devised within the differential geometric and the variational approaches and operating in accordance with Eq. 3 have been summarized in Fig. 1. POLYCRANK a unique design of a 7 d.o.f manipulator with revolute joints, distinguishing by practically unlimited range of joint motion as well as kinematic and dynamic isotropy [15,16]. Using the Cartesian position coordinates and Z X Z Euler angles, the kinematics of POLYCRANK can be represented as As usual, we have adopted the notations s i = sin x i and c i = cos x i . The geometric parameters of the manipulator are equal to l 1 = 0.2975m, l 2 = 0.18m, l 3 = 1.552m, l 4 = l 5 = 0.16m, l 6 = 0.2562m. The POLYCRANK's Jacobian is computed as It is easily checked that singular configurations of the POLYCRANK lie at x 5 = jπ or x 2 − x 1 = lπ , j, l = 0, ±1, . . .. The remaining configurations are regular.

Differential Geometric Approach
In the manipulator's joint space we choose a regular region D, and define its foliation by straight lines parallel to x 4 coordinate parametrized by α = (α 1 , α 2 , α 3 , α 5 , α 6 , α 7 ) ∈ R 6 , and dependent on the numbers a, b , c = 0, a = b whose role is to place the zero-leaf On its basis we compute the vector field (11) The 1-form annihilating the Jacobian transpose, and complying with the foliation is the following [12]  Finally, we need to introduce coordinates x = (y, z) such that z = x 4 changes along the leaves, and y i = x i , i = 1, 2, 3 and y i+3 = x i+4 , i = 1, 2, 3 are transverse to the leaves. Following the procedure described in Section 3, we arrive at the characteristic equation for the z coordinate, parametrized by y ∈ R m , Let ϕ(t, y, z 0 ) denote the flow of Eq. 27 initialized at z 0 . Then, the augmenting kinematics map satisfies the identity y, h(y, z)).

Variational Approach
In the variational approach the 1-form appearing in the error functional (21) is the same as in Eq. 26, while the integration extends over the region The error is minimized using the Ritz method.
Having chosen a quadratic augmenting kinematics function and performed the computations in accordance with Eq. 23, we obtain the following result The coefficients of h(x) have been computed by MATHEMATICA with a quite high precision. However, this precision is not crucial in the practical implementation.

Computations
The solutions of problems 1 and 2 provided by the approximate extended Jacobian algorithms and the original, Jacobian pseudo inverse algorithm, are displayed in Figs. 3 and 4. The first 7 plots show trajectories of each joint. The last plot, at the bottom right side, presents the norm of the instantaneous joint velocities ẋ(θ ) 2 for the 3 algorithms compared with the norm produced by a non-optimal extended Jacobian algorithm with a randomly selected augmenting kinematics map. It follows that both the geometric as well as the variational approach have solved the problems in a satisfactory way in the sense that the solutions stay acceptably close to the solution provided by the Jacobian pseudo inverse algorithm. Specifically, in Problem 2 the approximation provided by the variational approach has been more accurate than the differential geometric. To explain this difference, Fig. 5 shows two pages in the POLYCRANK's joint space, passing through the leaves E α of the initial and E β of the desired joint space positions, for Problems 1 and 2. It can be seen that the joint trajectories computed by the extended Jacobian algorithms do not lie on the pages passing through E α . If the trajectory belonged to the leaf, the approximation would be perfect, so intuitively, the closer to the leaf the trajectory stays, the better approximation is achieved. Taking the variability of ω as a measure of the closeness, we have computed the increment ω = ω(x(θ max )) − ω(x 0 ) when moving from the initial to the final joint space position. Assuming that θ max = 10, it turns out that in Problem 1 || ω|| = 0.490342, whereas in Problem 2 || ω|| = 0.62882. By design, a distinctive feature of the control strategy based on the Jacobian pseudo inverse is the minimization of the joint velocities of the manipulator during the motion. The bottom-right plots in Figs. 3 and 4 have confirmed that the instantaneous joint velocities obtained by the approximate extended Jacobian algorithms are close to those resulting from the Jacobian pseudo inverse algorithm. They are not close in the case of a non-optimal extended Jacobian algorithm. Furthermore, Table 1 collects the values of the norm θ max 0 ẋ(θ ) 2 dθ corresponding to the averaged joint velocity. An inspection shows that again the velocity norms for the approximate extended Jacobian and the Jacobian pseudo inverse algorithms are comparable, in contrast to the non-optimal extended Jacobian algorithm. Finally, let us make some comments on computational aspects of the presented synthesis methods of Jacobian algorithms. The structure of underlying computations has been shown in Fig. 1. Given an inverse kinematics problem, it follows that the extended Jacobian algorithm based on the differential geometric approach contents itself with computing numerically the augmenting function only at the configurations generated by the algorithm: in the case of redundancy r = 1, in order to determine the value of J E# (x θ ) in the θth step, one needs to integrate the characteristic equation (18) 2n times to find ∂h( . . . , n, where n equals the dimension of the joint space. On a PC equipped with AMD Phenom II X6 1055T 2.8 MHz processor with 8 Gb operational memory, working under Ubuntu Linux 10.04 x86 64 and using Matlab R2010a 64bit, the one step of computation for POLYCRANK takes in average 19.6 ms for the method of finite differences and 7.8 ms, when the forward sensitivity analysis method is applied. This means that, potentially, the differential geometric approach can be implemented in real time. The computations involved in the variational approach are organized differently: the augmenting function can be found offline in the closed form, once over the whole region D ⊂ R n of the joint space, plugged into the extended Jacobian inverse, and then used to solve inverse kinematic problems in real time. Computationally, in solving a single inverse kinematic problem for POLYCRANK the variational approach has been from about 5 (linear augmenting function) to about 60 (quadratic augmenting function) times more time-consuming than the differential geometric approach. Therefore, on condition that approximation accuracies are similar, the differential geometric approach could be recommended as a rapid solution method of a single or a few inverse kinematic problems. Contrary to that, the variational approach results in an extended Jacobian algorithm that can be used systematically to long series of inverse kinematic problems, that justifies the investment into the determination of the augmenting function. Once it has been determined, the computations can be accomplished in real time. Observe that both the approaches involve only numerical, not symbolic, computations.

Extension to r > 1
The presented approaches generalize in a natural way to the case of redundancy r > 1. As-sume as a starting point assume the codistributions (10). Then, in the differential geometric approach the leaf E α will be r dimensional, with coordinates z 1 , . . . , z r chosen in such a way that ω i | E α = dz i . Consequently, Eqs. 13 and 14 will take the form ω i | E α = dz i = r j B ij dh j and F = ω i X = r j B ij dh j X for an r × r matrix B with entries B ij smoothly depending on x. After suitable mathematical developments we arrive at a vector version of Eq. 18, dz(t) dt = − F( t (y,z)) t , where x = (y, z) and F(x) = (F 1 (x), . . . , F r (x)). The augmenting kinematics maps h = (h 1 (x), . . . , h r (x)) are obtained from the analog of Eq. 19. The details can be found in [12] along with a computational example for r = 2.
Within the variational approach the approximation error (21) for r > 1 will be defined as the integral of the squared Frobenius matrix norm Assuming that h i (x) = a T i ϕ(x) we compute the optimal parameter vectors a * i = M −1 N i where N i = D Dϕ(x)ω i (x)dx, i = 1, . . . , r. The decomposition of the minimization problem into independent problems for i = 1, 2, . . . , r unveils an additional advantage of the error formula (21).

Conclusion
We have examined two approaches to devising extended Jacobian inverse kinematics algorithms that approximate the Jacobian pseudo inverse algorithm. To enable a comparison, a codistribution framework has been employed, and a novel approximation error functional introduced. These approaches have been applied to the kinematics model of the 7 dof industrial manipulator POLY-CRANK. The computations show that, basically, both the methods offer satisfactory approximations. Specifically, the quality of approximation resulting from the differential geometric approach depends on the variability of the 1-form (26) around the pages in the joint space, that depends on the choice of the joint space foliation and, ultimately, on the displacement of the kinematic singularities. It is an open problem how to increase the quality of approximation by a proper definition of this foliation. In the case of the variational approach the quality of approximation is directly measured by the error functional, and is dependent on the choice of the basic functions involved in the Ritz method as well as on the accuracy of computation of the integrals appearing in the error formula (23). As has been said above, potentially, both the approaches can be used in real time. Specifically, if the approximation accuracy is not the primary issue, the differential geometric approach is preferred as a rapid, whereas the variational approach as a systematic tool of solving inverse kinematic problems for robotic manipulators.