A Lie Group-Based Iterative Algorithm Framework for Numerically Solving Forward Kinematics of Gough–Stewart Platform

: In this work, we began to take forward kinematics of the Gough–Stewart (G-S) platform as an unconstrained optimization problem on the Lie group-structured manifold SE(3) instead of simply relaxing its intrinsic orthogonal constraint when algorithms are updated on six-dimensional local ﬂat Euclidean space or adding extra unit norm constraint when orientation parts are parametrized by a unit quaternion. With this thought in mind, we construct two kinds of iterative problem-solving algorithms (Gauss–Newton (G-N) and Levenberg–Marquardt (L-M)) with mathematical tools from the Lie group and Lie algebra. Finally, a case study for a general G-S platform was carried out to compare these two kinds of algorithms on SE(3) with corresponding algorithms that updated on six-dimensional ﬂat Euclidean space or seven-dimensional quaternion-based parametrization Euclidean space. Experiment results demonstrate that those algorithms on SE(3) behave better than others in convergence performance especially when the initial guess selection is near to branch solutions.


Introduction
The six degree-of-freedom Gough-Stewart (G-S) type motion platforms have been extensively employed to realize fully operational flight training simulation with high fidelity.This spatial G-S platform had six identical extensible active legs connected with an upper moving platform and a fixed lower platform using six upper passive joints and six lower passive joints.The forward kinematics of the G-S platform allows the platform to configure the pose (position and orientation) of end-effectors on the moving platform when the lengths of the six linear actuators are provided.However, the closed loop kinematic relations between the moving platform and the fixed base platform in 3-D space made forward kinematic problems a challenging issue among all kinematic research of parallel manipulators during the past few decades.
A significant amount of research has been focused on the forward kinematics of the G-S platform, which has lead us to understand the difficulty of this problem.Generally, there exist three main research directions to address the forward kinematics problem, which are the analytic method, numerical iterative method, and auxiliary-sensor-based method.The majority of early research focused on finding all possible closed-form direct solutions for different types of G-S platforms.Since Raghavan [1] and Ronga [2] firstly proved that there existed at most 40 distinct solutions in the complex domain for a general G-S platform, many researchers succeed in determining all 40 direct solutions using algebraic elimination, interval analysis, or numerical continuation methods, etc. Husty [3] obtained a 40-th degree univariate polynomial through complicated algebraic elimination based on a spatial kinematic mapping.For a symmetrical 6-6 G-S platform, the degree of polynomial was ultimately reduced to 14 by Huang [4].Researcher Ji [5] focused on a special type 6-6 G-S platform in which both the fixed base platform and the moving platform are similar hexagons and only univariate quadratic equations are required to derive eight possible symmetry solutions by introducing a quaternion to represent the orientation matrix.As a general analytic form of direction solutions to forward kinematics of a general G-S platform was luxurious, some scholars resorted to the numerical root-finding method.Both Didrit [6] and Merlet [7] preferred to utilize the interval analysis method to numerically solve the forward kinematic problem for the general G-S platform.Wampler [8] implemented a polynomial continuation scheme to deal with this problem.Besides that, Dietmaier [9] proposed a numerical procedure to change the parameters of the G-S platform and ultimately obtained a particular asymmetry for the G-S platform, which possesses 40 real direct solutions.
As there exist multiple direct solutions, research attempting to solve the forward kinematic problem is further expected to determine a unique actual pose (position and orientation) of the G-S platform.Numerical iterative methods and auxiliary-sensor-based methods are two common schemes to pursue this goal.As the forward kinematic problem was easily boiled down to seeking six dimensional variables, satisfying six dimensional kinematic equations, Newton-type iterative algorithms were extensively put into use because of its simplicity and fast convergence speed [10,11].However, it had been well known that the Newton-type algorithm usually converges to a particular solution near to initial guess value.Therefore, some variants of the Newton-type algorithm had been discussed to improve its performance in the following studies.Yang [12] modified the traditional Newton-Raphson algorithm using a monotonic descent operator to achieve global convergence regardless of the initial guess selection.Pratik [13] proposed a neural-network-based hybrid strategy that combined with the standard Newton-Raphson algorithm to yield better performance in dealing with uniqueness problem.Furthermore, Innocenti [14], Parenti-Castelli [15], and Chiu [16] suggested adding one or more redundant sensors to determine a unique real solution.Wang [17] designed a incremental iterative algorithm through a series of continuous small changes in leg length to derive a unique forward kinematic solution.
When iterative algorithms are employed to solve optimization that involves the SE(3) configure pose, the problem of which parametrization scheme is beneficial to design its updating framework must be addressed.As we all know, 3D+YPR and 3D+Quat are two common parametrization schemes, as explained in Reference [18].In the 3D+YPR scheme, the iterative algorithm can be designed on flat Euclidean space R 6 with orthogonal constraints not taken into account.In the 3D+Quat scheme, a unit norm constraint equation is added while state vector is updated on enlarged Euclidean space R 7 .In this work, we address the forward kinematic as an unconstrained optimization problem on a Lie group-structured manifold SE(3).In algorithm construction, the updating direction and step size are computed on the associated tangent space T e SE(3) with differential properties of the Lie group and Lie algebra.The result is then projected back on the manifold SE(3) to update the next configure pose.The advantage in the optimization on manifold SE(3) is that the exponential map that connects the Lie algebra and the Lie group naturally preserves the orthogonal constraint during each iteration.An experimental comparison is expected to show that the ways of simply updating the iteration on the local parameter space make the forward kinematic optimization algorithms more susceptible to being stuck in branch solutions.
The rest of this paper is organized as follows.In Section 2, some preliminary notations and terminology from the theory of Lie group and Lie algebra is reviewed.Section 3 is devoted to formulating the problem of forward kinematics for a general type G-S platform.Two kinds of iterative algorithms framework, the Gauss-Newton (G-N) and Levenberg-Marquardt (L-M) type algorithms are employed in algorithm construction with the mathematical tools from Lie group in Section 4. In Section 5, these two types of Lie group-based iterative algorithms are compared with traditional algorithms on Euclidean space R 6 as well as iterative algorithms using quaternion parametrization through a case study for a general G-S platform.

Lie Group and Lie Algebra
We assume that readers are familiar with basic mathematical concepts about group and manifold (see e.g., [19,20] and references therein).The configure pose of the G-S platform end-effector belongs to a special Euclidean group SE (3), which is a semi-direct product of R 3 with special orthogonal group SO (3).An element g ∈ SE(3) can be presented using a 4 × 4 homogeneous transformation matrix as follows.
where R ∈ SO(3) is a 3 × 3 orientation matrix, p ∈ R 3 is a 3 × 1 translational vector.The inverse of g can be written in the following formulation, Given a basis {E 1 , • • • , E 6 } for matrix Lie algebra se(3), any arbitrary element X, which is also called a "screw" matrix in the robot community, can be written in the following form, where E 1 , E 2 , • • • , E 6 are listed as Equation (A1) in Appendix A. Given a rigid body motion g(t) ∈ SE(3), its associated Lie algebra element shown as Equation (4) denotes a spatial screw matrix S r in an inertial frame.Meanwhile, Equation (5) denotes a screw matrix S b in body reference frame.
Accordingly, there exist two 6 × 1 column vectors, spatial screw vector s r and body screw vector s b , which can be derived using following operation.
Here the vee operator (•) ∨ : se(3) → R 6 is defined to extract screw coefficients from the screw matrix to form a 6 × 1 column vector.Conversely, a wedge operator (•) ∧ : R 6 → se(3) maps the screw vector back to the matrix Lie algebra se(3), that is s

Exponential Map
For each vector field X in tangent space T e G at group identity e, there exists a smooth one-parameter subgroup γ X (t) of a Lie group G parametrized by a scalar t ∈ R.There also exists an exponential map defined as follows.
This is a unique homomorphism map from tangent space T e G around group identity e to one-parameter subgroups of G.The exponential map is related with a one-parameter subgroup by exp(X) = γ(1) at t = 1.
With these properties of the exponential map, we can define an integral curve of right invariant vector field X r passing through point g ∈ G at time zero as γ r X : R denotes an integral curve of left invariant vector field X l passing through g ∈ G at time zero.
As for special the Euclidean group SE(3), there exists an associated screw matrix as shown in Equation ( 3) at its tangent space T e SE(3).Then, motion around identity I 4 can be obtained by exponentiating these screw matrices.Here, the exponential map exp : se(3) → SE(3) corresponds to matrix exponentiation, which has the closed-form Rodriguez formula as shown in Equation (A2) in the Appendix A. Examples for rotational screw vector x 1 and translational screw vector x 6 are expressed as follows,

Taylor Series Expansion of an Analytic Function on Lie Group
The Taylor series of a function on Euclidean space R n can be extended naturally to Lie group-structured non-Euclidean space.In the same way, there exists similar derivatives for group-valued function f : G → R. Owing to the fact that there exist two invariant vector fields around g ∈ G, The Taylor series expansion around g ∈ G can be formulated in two ways.The following equation provides the right way.
As discussed in [21], the first and second-order term can be expressed in a different form.Analogously to directional derivative, the first-order term are completely equivalent to the following formulation, Here, (X r f )(g) is called right Lie derivative of f (g) with respect to vector filed X.Since E i is a basis for Lie algebra g, which can be shown as X = ∑ n i=1 x i E i , then its associated differential operators will be denoted as E r i f .Thus, the right Lie derivative will be as follows, Finally, the "right" Taylor series expansion around g ∈ G can be written to a second order in the following equation.
A "left" Taylor series expansion can also be derived in an analogous manner.

Problem Formulation of Forward Kinematics for G-S Platform
This section is devoted to formulating the forward kinematic problem of a general type G-S platform as a minimum optimization problem.The ultimate problem formulation depends on which parametrization scheme of configure pose is applied.Here, we take 3D + YPR scheme as an example, where position is in three-dimensional flat Euclidean space R 3 and orientation is parametrized by three angles in Yaw-Pitch-Roll order.
As for a general G-S platform, its geometric structure of a general G-S platform is depicted in the following Figure 1.
where orientation matrix R ∈ SO(3) is defined by successive rotation in Z-Y-X order.Let a i denote the position vector of the lower joint center A i of i-th leg in frame {N}, and let b i denote the position vector of upper joint center B i of i-th leg in frame {B}.
The inverse kinematic function of the G-S platform defined as a six dimensional real-valued function F : SE(3) → R 6 describes a kinematic mapping from configure pose T of upper moving platform to its six leg lengths L 1 ,. . ., L 6 .If T is locally parametrized by 3D+YPR generalized coordinate variable q ∈ R 6 ; thus, the ultimate inverse kinematic function F(q) is defined as follows.
The forward kinematic problem of the G-S platform is formulated as a minimum optimization problem related with kinematic mapping residual error r(T).Each component of r(T) is defined as follows.
where L 1 , L 2 , • • • , L 6 are six scalars for current leg length.Let ϕ denote the objective function in the following formulation.
As shown in Figure 2, local parametrization of objective function ϕ around configure point T is described as follows, Finally, forward kinematic problem of the G-S platform can be formulated as a minimum optimization problem as follows.min where domain D denotes accessible workspace configure pose of the G-S platform.

Iterative Algorithms on SE(3)
As discussed in the introduction section, it will take considerable efforts before we can arrive at an analytical solution to solve the forward kinematic problem in Equation (20).In consideration of this difficulty, many numerical iteration algorithms, e.g., G-N and L-M, have been extensively applied in practical engineering projects.
In order to minimize the quadratic sum of kinematic residual error r(T), the numerically iterative algorithm is updated iteratively by means of a small increment until convergence or the max number of iteration is reached.It should be noticed that all these numerical iteration schemes are originally designed to work on flat Euclidean space, i.e., R n .If the update process was designed on generalized coordinates q ∈ R 6 according to Equation ( 16), there exists a gimbal lock problem as well as three Euler angles that may reach out of their valid ranges.In this section, we begin to build up two elegant iterative algorithms on the Lie group-structured manifold SE(3) with differential properties from theory of the Lie group and Lie algebra.
The core issue in the Newton-type algorithm construction process is to determine the small increments at each step through solving incremental normal equations related with the first-order Taylor series expansion.According to Equation ( 13), the six-dimensional inverse kinematic residual error function r(T) can be approximated in the neighborhood of T ∈ SE(3) with first-order "right" Taylor series expansion as follows.
As for computing the first order term, right Lie derivative of six-dimensional vectorvalued residual function r(T) was defined as X r r = [X r r 1 , X r r 2 , . . ., X r r 6 ] T .Each Lie derivative with respect to body screw matrix S b can be derived using the following formulation.
where tilde on a i and b i means a extension by adding an extra element 1, that is ãi = (a T i , 1) T and bi = (b T i , 1) T .Similarly, a left Lie derivative with respect to spatial screw matrix S r can be used to approximate r(T) in "left" Taylor series expansion.
All these six right Lie derivatives can be formed into a 6 × 6 Jacobian matrix, which is named right Lie derivative matrix J b with respect to body screw matrix S b as follows.
Then Equation ( 21) can be arranged in form of body screw vector s b as follows.)) = 0.This means that the updating direction in Lie algebra vector space at time step t k can be given as follows.
After its updating direction has been chosen, a step factor α k (0 < α < 1) that satisfies the following monotonic descent rule is dedicated to constraint step size at time t k so as to avoid divergence during iterative process.If initial setting of step factor α k cannot meet this rule, then multiply by itself α k = α m k until the descent rule is met.
Once a descent step size has been determined, the next configure pose can be updated by mapping selected descent screw error in Lie algebra to Lie group using the exponential map.
In the end, a Lie group-based G-N algorithm for the forward kinematics of the G-S platform can be constructed as the Algorithm A1 in the Appendix B. Figure 3 provides an intuitive explanation of this algorithm.However, it is well known that G-N type algorithms depend heavily on good local linearization around selected expansion points.As to this issue, L-M type algorithms behave better in adjusting to a descent direction within an appropriate dynamically tuned trust region.An L-M type iterative algorithm under the framework of the Lie group is shown in Algorithm A2 in the Appendix B.

Implementation of Algorithms and Discussions
In this section, we first compare the aforementioned Lie group-based G-N type iterative algorithm (GN-SE(3)) with corresponding G-N type algorithms (GN-Quat and GN-R 6 ).Then, it is followed with three L-M type algorithms (LM-SE(3), LM-Quat, and LM-R 6 ).A general type G-S platform is leveraged to test the performance of these algorithms.The following Table 1 provides those lower joint coordinate vectors on fixed platform {N} and upper joint coordinate vectors on moving platform frame {B} as referred in [22].
Table 1.Coordinate vectors of the G-S platform's upper and lower joints with respect to their own coordinate frame.

Lower Joint Coordinates on Frame {N} (cm)
Upper Joint Coordinates on Frame {B} (cm) Take a group of current leg lengths as our destination L d , which corresponds to the true value of platform pose q C = [0, 0, 50, 20  ] T that correspond to the same lengths group L d .Although the two branch solutions q S1 and q S2 are represented with different Euler angles, their orientation parts actually correspond to the same orientation matrix.Therefore, both branch solutions q S1 and q S2 should be regarded as one configure pose T s as follows.The following Figure 5 gives its platform state in 3D Cartesian space.In algorithm initialization phase, five configure poses as shown in Table 2 are provided as initial guess states to testify to three kinds of G-N type forward kinematic algorithms.All these initial states are chosen to be far away from the true value q C .The choice of these five configure poses are based on the performance that G-N and L-M algorithms on R 6 (GN-R 6 and LM-R 6 ) act at these initial guess states.State 1, state 2, and state 4 represent a group of initial guess values that GN-R 6 and LM-R 6 deteriorated quickly, while state 3 and state 5 represent the opposite type of initial guess value, these two algorithms converged to the true value normally.
The other algorithms that were constructed on quaternion-based Euclidean space R 7 and SE(3) are expected to improve the convergence performance when starting from state 1, state 2, and state 4 because both of them have kept the orthogonal constraint at each iteration process.Meanwhile, state 3 and state 5 are necessary moderate choices to verify the reliability of the other two algorithms.In order to evaluate the influence of the step factor α on algorithm performance, the step factor α is selected from 0.50 to 0.99 with an equal interval of 0.01.The reason to choose a relatively large value for the step factor is to make sure at least half of the derived updating direction could be exploited at each iteration.All tested G-N type algorithms share the same accuracy threshold ( 1 = 1.0 × 10 −14 , 2 = 1.0 × 10 −14 , 3 = 1.0 × 10 −14 ), and max iteration time k max = 200 to end their iteration.Figure 4b-f corresponds to these five initial guess states of the G-S platform, respectively.In this experiment, two criteria were employed to evaluate algorithm performance, one is the range of step factor that converges to true value q C , the other is convergence speed.The algorithm performance reflects significant difference among these five initial states.Their performance difference comparisons are shown in Figures 6 and 7, which take step factor α and iteration time as their coordinate axes.Let the positive direction of longitudinal axis indicate the iteration time that an algorithm reaches final true value q C within given accuracy threshold while negative direction indicates the iteration time that it takes to converge to other branch solutions.This comparison analyst yields the following conclusions.

1.
For GN-R 6 , it achieves fast convergence when it starts from state 2, state 3, and state 5.However, it is more likely to reach one solution q S2 when state 2 is selected.For state 4, the algorithm is inclined to another branch solution q S1 within some isolated range.As to state 1, due its divergence under the large range of the step factor, it means that state 1 is not an appropriate for the GN-R 6 guess value.

2.
For GN-Quat, it ends up with the true value q C among some isolated range of step factor when starting from state 1.For state 2, it always converges to q S2 .For state 3 and state 5, although it achieves convergence to the true value q C , but it could possibly behave worse if the convergence speed is slowed down at some specific larger step factors.The state 4 is definitely no longer a suitable choice for GN-Quat algorithm initialization.

3.
For GN-SE(3), when it starts from state 3 and state 5, it always converges to the true value and the larger step factor would be helpful to reduce iteration time.When state 1 is provided as the initial guess state, it converges to the true value q C with a probability of more than fifty percent.For state 4, it converges to the true value with a consecutive range from 0.60 to 0.84.
In the second comparison experiment, these five configure poses were provided as input to L-M type algorithms (LM-R 6 , LM-Quat and LM-SE(3)).Their performance comparison with respect to the damping ratio τ is given in Table 3.The logarithmic of the damping ratio τ that was taken as horizontal coordinate is assigned to a range between −9 and −3.12 with an equal interval of 0.12.Its iteration time is taken as its longitudinal coordinate.All these tested L-M type algorithms share the same accuracy threshold ( 1 = 1.0 × 10 −14 , 2 = 1.0 × 10 −144 ) and max iteration time k max = 200 to terminate their iteration.
As seen in Table 3 and Figure 8, if state 3 is employed for algorithm initialization, all three kinds of L-M type algorithms converge to the true value q C at an acceptable speed in the whole range of the damping ratio τ.However, there exists obvious performance differences between these L-M type algorithms when the four other states are used to initialize these algorithms.Their performance comparisons are arranged in Figures 8 and 9. From these comparisons, we can draw the following conclusions, 1.
For LM-R 6 , when its is initialized from state 5, it achieves fast convergence toward true value q C in the whole range of damping factor.However, it is more likely to converge to branch solution q S2 when state 2 is employed.If state 4 is taken into consideration, it reaches the true value q C within some isolated regions for damping factor selection.As to state 1, it has a higher chance to converge to the branch solution q S1 within some isolated regions.

2.
For LM-Quat, it would converge to the true value q C when it starts form state 5.
From state 2, it also achieves almost total convergence towards branch solution q S2 .For state 1, it is definitely not suitable for initialization as it cannot converge to any solution within the whole selection range.Finally, when the algorithm initialized from state 4, it can achieve some convergence to the true value q C within some isolated regions for damping factor τ.

3.
For LM-SE(3), it can converge to the true value q C with almost the whole range of damping factor selection except some large values.For state 2, it behaved better in convergence to the true value q C when the logarithm value of the damping factor is selected from −5.5 to −3.12.For state 1 and state 4, its convergence probability towards the true value q C is quite small, which means that the two states cannot be appropriate for LM-SE(3) algorithm initialization.

Iteration time
Solution 1 (a) LM-R   Finally, we consider two evaluation indicators for the initial guess state selection, convergence probability toward the true value q C and consecutive region for influence factors' selection.For GN type algorithms, state 1, state 4, and state 5 are totally or conditionally suitable for LM-SE(3).State 5 is the only suitable initial guess state for GN-R 6 and conditionally acceptable for GN-Quat.For LM type algorithms, state 5 is totally suitable initial guess state choice for all three algorithms.Meanwhile, state 2 is conditionally acceptable for LM-SE(3) algorithm.
In summary, we can arrive at some conclusions through the above comparative analysis.Firstly, the iterative algorithms constructed on SE(3) would leave more space for initial guess state selection in the workspace of the G-S platform.Secondly, the algorithm on SE(3) has few chances to be stuck in other branch solutions even though the initial guess value is near to the branch solution.

Conclusions
In this paper, we discussed how to construct G-N and L-M type iterative algorithms (GN-SE(3) and LM-SE(3)) to solve the forward kinematics problem of a G-S platform with mathematical tools from the Lie group.The key part of the group-based iterative algorithm was to determine an updating direction in Lie se(3) and was then projected back to the Lie group with exponential map, which differentiates the algorithm from traditional iterative algorithms based on locally Euclidean space R 6 or quaternion-based parametrization space R 7 .Five different initial configure poses of a general G-S platform were used to compare Lie group-based algorithms (GN-SE(3) and LM-SE(3)) with two other kinds of iterative algorithms.Experiment results demonstrated that Lie group-based iterative algorithms behave better in converging to the true solution while the two other algorithms leave more chance in converging to another branch solutions.In the future, more work is needed to study how geometric parameters of the G-S platform can influence the choice of the iterative forward kinematic algorithm.
Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Appendix A.1 The following E 1 , E 2 , • • • , E 6 are six basis components for matrix Lie algebra se(3).

Figure 1 .
Figure 1.Geometric structure of a general Gough-Stewart (G-S) platform.It consists of one lower fixed platform and one upper moving platform connected through six identical extensible legs.Both two platforms are irregular hexagon structure with their six vertices located at one circumcircle.In order to describe the relative configure pose between two platform, one coordinate frame {O N − X N Y N Z N } is fixed at geometric center of lower fixed platform O N , another frame {O B − X B Y B Z B } is fixed at geometric center of upper moving platform O B .Let T ∈ SE(3) denote the configure pose (position and orientation) of geometric center O B of upper moving platform represented in lower fixed frame {N}.
r(T • exp(tS b )) ≈ r(T) + J b (T)ts b (25) Thus we can obtain body screw error s b = −J −1 b (T)r(T) by assuming r(T • exp(tS b

Figure 3 .
Figure 3. Updating process of G-N type algorithm on Lie group SE(3).

Figure 4 .
Figure 4. True state and five initial guess states of the G-S platform.

Figure 5 .
Figure 5. Platform state T s corresponding to branch solutions q S1 and q S2 .

Figure 6 .
Figure 6.G-N algorithm performance comparison for initial guess state 1, state 2, and state 3.

Figure 7 .
Figure 7. G-N type algorithm performance comparison for initial guess state 4 and state 5.

Figure 9 .
Figure 9. L-M algorithm performance comparison for initial guess state 4 and state 5.

Table 2 .
Five initial configure poses of the G-S platform and performance comparison of G-N type algorithms.

Table 3 .
Figure 8. L-M type algorithm performance comparison for initial guess state 1, state 2, and state 3. Five initial configure poses of the G-S platform and performance comparison for L-M algorithms.