Development of Task-Space Nonholonomic Motion Planning Algorithm Based on Lie-Algebraic Method

: In this paper, a Lie-algebraic nonholonomic motion planning technique, originally designed to work in a conﬁguration space, was extended to plan a motion within a task-space resulting from an output function considered. In both planning spaces, a generalized Campbell–Baker–Hausdorff– Dynkin formula was utilized to transform a motion planning into an inverse kinematic task known for serial manipulators. A complete, general-purpose Lie-algebraic algorithm is provided for a local motion planning of nonholonomic systems with or without output functions. Similarities and differences in motion planning within conﬁguration and task spaces were highlighted. It appears that motion planning in a task-space can simplify a planning task and also gives an opportunity to optimize a motion of nonholonomic systems. Unfortunately, in this planning there is no way to avoid working in a conﬁguration space. The auxiliary objective of the paper is to verify, through simulations, an impact of initial parameters on the efﬁciency of the planning algorithm, and to provide some hints on how to set the parameters correctly.


Introduction
A large number of contemporary and practically important robots can be described as nonholonomic systems. Despite different sources of nonholonomy and nonholonomic manipulators, wheel mobile robots [1] and free-floating space robots [2] belong to this class. Even when modeled at the kinematic level, nonholonomic systems are difficult to control because the number of controls is smaller than the dimension of their configuration spaces. Thus, sophisticated maneuvers have to be performed to get a desired location, clearly exemplified by the task of parking a car (especially towing trailers). For some special cases of nonholonomic systems (like Dubbins [3] and Reeds-Shepp cars [4] or systems in a chain form [5]), dedicated methods were proposed to plan their motions. However, there is a constant need to develop general purpose methods applicable to any nonholonomic system. Two main classes of general purpose methods can be distinguished: global [6,7] and local ones [8]. Global methods try to get a whole output trajectory at once (although in an iterative manner), while local ones construct an output trajectory from sub-trajectories derived from local planning. Global methods are rooted in various paradigms. Using a linearization of a nonholonomic system along the trajectory corresponding to current controls, Tchon and coworkers [7] reformulated a nonholonomic motion planning task into an inverse task solved with classical methods. Jakubczyk [9] applied the Volterra series expansion of a nonholonomic system with an output function to get highly oscillatory controls that were able to trace a desired trajectory with an arbitrary precision. Arismendi and coworkers [10] adapted a variant of a wave front technique to design the fast marching square path planning method of nonholonomic motion planning.
Generally, algorithms derived from global methods are computationally involved and very sensitive to geometric constraints due to obstacles. Consequently, they cannot be 1. to extend applicability of the method to plan a motion also within a task-space, and to provide a complete algorithm for implementing the method, 2.
to show through simulations how data parameterizing the algorithm impacts its efficiency, 3. to compare planning within a configuration and a task-space, pointing out similarities and differences.
It is worth noticing that local and global planning tasks belong to a classic triple: geometrical planning, motion planning, and control. Recently, some attempts have been made to combine the last two tasks into one planning-control loop [11,12]. A vast majority of motion planning algorithms are based on discretization of control/configuration spaces to transform continuous tasks, described by differential equations and static output functions, into a graph domain. In this discrete space, graph-searching algorithms, like RRT or its numerous variants [13], can be used to solve the planning task.
The paper is organized as follows: In Section 2 a model of nonholonomic systems is introduced supplemented with an output function, then a motion planning task is defined and, finally, a Lie algebraic algorithm to solve the task is presented. In Section 3 an extensive simulation study is provided to evaluate performance of the algorithm for various data that impact its behavior. Then, some observations are formulated based on the simulation results. Section 4 concludes the paper.

Model and Algorithm
Nonholonomic systems modeled at a kinematic level are described by the controlaffine equation [14]q where Q q q q is a configuration, u u u = (u 1 , . . . , u m ) T are controls, and g g g i (q q q) are (smooth) vector fields (generators) that span a null space of nonholonomic constraints. By definition, nonholonomic systems satisfy the Lie algebra rank condition. Thus, according to the Chow's theorem [15], they are small time locally controllable. In order to shorten notations, two-input systems, m = 2, will be considered with g g g 1 = X X X and g g g 2 = Y Y Y. Two input systems are the most challenging as for the systems it is quite easy to generate difficult motion planning tasks by increasing a configuration space dimension.
In many practical applications not all components of a configuration vector are important (for example wheel angles of mobile robots do not impact obstacle collisions), thus an output map is introduced and usually r < n. Now, a motion planning task can be stated as follows: given an initial configuration q q q 0 and a goal point, x x x f find controls u u u(·) which, applied to the system (1), generates an end-point of the resulting trajectory mapped with Equation (2) into the desired location x x x f .
In Lie-algebraic methods, locally, around a current configuration q q q c , possible directions of motion of the system (1) are described by vector fields derived from generators X X X, Y Y Y. To avoid redundancy, it is a common practice to take only a basis of the Lie algebra spanned by the generators. There are at least three such bases, but probably the most frequently used is the Ph. Hall basis. It is composed of an infinite number of Lie monomials, interpreted as vector fields within the scope of the control system (1). The very first elements of the basis are given below The basis is composed of generators (degree 1 Lie monomials) and higher degree monomials derived from generators and their descendants. The system (1) is a nonholonomic one, thus a finite (and usually small) number of initial Ph. Hall basis elements, H H H trunc , can be taken to satisfy the Lie algebra rank condition [15] ∀q q q c ∈ Q rank H H H trunc (q q q c ) = n.
Consequently, the system (1) is a small time locally controllable and it can evolve at any configuration in any direction. The condition (4) is quite easy to check in off-line mode with the assistance of symbolic computation packages. However, a very demanding problem is how to generate a motion in any direction with a small number of controls. Fortunately, the generalized Campbell-Baker-Hausdorff-Dynkin formula [16] answers this question constructively. It states that a configuration shift z z z(t)(q q q c ) at a current configuration q q q c , depends on controls u u u [17] as follows where control-dependent coefficients α α α = (α 1 , α 2 , . . .) T pre-multiplying Lie monomials (vector fields) are expressed as and they become more and more complex as the integration is performed over the k-dimensional simplex, where k is the degree of a Lie monomial the coefficient corresponds to. It is a common practice to express controls as linear combinations of time-dependent functions [7] u i (t) = where φ j (t) are elements of any functional basis (polynomials, harmonic functions) and N i is the number of variables required to describe the i-th control. Thus, a vector p p p collects all variables p j i and uniquely describes controls u u u and also energy of controls. Moreover, the shift (5) can be expressed as a function of p p p and resembles forward kinematics (but at a velocity level when considered on a fixed time horizon) z z z(t)(q q q c , p p p) = F F F(p p p).
The shift (8) from the configuration space is mapped into the task-space, and within this space a desired motion towards the goal x x x f can be computed using the Newton algorithm. A complete algorithm to plan a motion of the system (1) with the output function (2), either in a configuration or a task-space, is presented in Algorithm 1.

Algorithm 1 nonholonomic motion planning within a task-space
Step 1. Read-in initial data: the system (1) together with the output function (2), the initial configuration q q q 0 and the goal point x x x f , desired accuracy of reaching the goal , the initial value of p p p c .
Step 2. Initialize empty resulting trajectory q q q res (·), and set the current configuration q q q c ← q q q 0 .
Step 3. Check the stop condition: iff then stop computations and output resulting trajectory q q q res (·). Otherwise, progress with Step 4.
Step 5. At a current point in the task-space x x x c = k k k(q q q c ), select the planned shift by setting a positive coefficient ξ.
Step 6. Using the Newton algorithm, Ref. [18] find p p p that solves the inverse kinematic task where J J J(q q q) = ∂k k k(q q q)/∂q q q is the Jacobi matrix of the mapping (2).
Step 7. Check whether the solution is acceptable: x where q q q (T) is the end-point of the trajectory q q q (·) initialized at q q q c and generated by the system (1) with controls u u u (t), t ∈ [0, T], corresponding to the vector of parameters p p p , cf. Equation (7).
update the current configuration q q q c ← q q q (T) 3.
and go back to Step 4.
Otherwise, decrease the coefficient ξ and go to Step 5.
Comments on Algorithm 1: • A nonholonomic motion planning in a configuration space is achieved by setting the identity output function (2). • A selected parameterization (7), Step 4, can vary from one iteration to another. Usually, one or at most a few parameterizations can be exploited, and vectors α(u u u(p p p)) from Equation (5) can be computed in an off-line mode. • In Step 7, a minimal requirement was formulated, i.e., a new end-point of trajectory in the task-space should be closer to the goal than it was in the previous iteration. More restrictive requirement assumes that it should also be close enough to the planned sub-goal x x x c + ∆x x x. • Singularities in solving (11) result from a rank deficiency of the Jacobi matrix, derived from the right hand side of (11) Thus, singularities may arise 1. either due to a rank deficiency of J J J(q q q c ) (they are avoidable by a small variation of controls which results in a small shift of q q q c ) 2. or a rank deficiency of the Jacobi matrix ∂F F F(p p p)/∂p p p at the current p p p c . Equivalent to the matrix ∂α(u u u(p p p))/∂p p p, cf. Equations (5), (7), and (8) as the full rank matrix H H H trunc (q q q c ) cannot introduce singularities at all, cf. Equation (4).

• In
Step 5, the positive coefficient ξ should not be too large, as values of matrices J J J(q q q c ), H H H trunc (q q q c ) T are computed at the current configuration q q q c and neglecting some vector fields, while truncation H H H → H H H trunc , is justified only locally. Its value cannot be too small either, as many iterations might be required to complete the algorithm. • In motion planning within a task-space with r < n, it is impossible to avoid planning within a configuration space, as the mapping (2) is not unique in this case.

Simulations
In order to evaluate Algorithm 1, some tests were performed using two models of mobile robots depicted in Figure 1. The unicycle is described by equationṡ The configuration vector q q q is composed of position (q 1 , q 2 ) and orientation q 3 of the vehicle on a plane, while u 1 , u 2 denote controls interpreted as angular and linear velocities.
The model of a kinematic car is given aṡ where configuration q q q has got an extra component q 4 describing the steering wheel angle. For simplicity, it was assumed that L = 1. Both models are nonholonomic ones because and for the kinematic car, satisfy Condition (4). For both models, the position on the xy plane was selected as the output function If a robot is inscribed in a circle of sufficient radius, its orientation does not influence the possibility of collisions with obstacles. For the kinematic car another output function was used, coupling position on the xy plane and orientation of the vehicle In practice, the orientation of the steering axle is not relevant. In order to check the motion planning Algorithm 1 within a configuration space, the identity output function was selected for both mobile robots.
In all tested cases, an initial configuration was set to the vector composed of zeroes only In the tests, a side-way motion maneuver on the xy-plane was planned for both models and all output functions. The maneuver is challenging as it requires a motion along higher (than generators X X X, Y Y Y) degree vector fields. Depending on a dimension of a task-space, goal points were selected as Any single task was run using three parameterizations of controls. Each control was composed of a few very first items of the orthonormal Fourier basis where ρ = √ 1/T and ω = 2π/T. Later on, it was assumed that T = 1. The total energy expenditure on controls is equal to where #N i is the number of the Fourier basis elements contributing to the i-th control. All tasks, with the data collected in Table 1, were solved using Algorithm 1. In Step 6 of the algorithm, two versions of the Newton algorithm were tested: the basic one and with the energy (25) optimization within the null-space of the Jacobi matrix (13). In order to obtain statistically significant results, Step 6 was repeated 100 times for each task, with initial values of p p p varied (each component of p p p was generated randomly with the uniform distribution on the interval [−1, 1]). Programs used in simulations were written in Wolfram Mathematica, version 11.3, and run on a computer with Intel® Core™ i5-8400 CPU and 24 GiB RAM memory. In Figures 2 and 3, some characteristic paths (projected on the xy-plane and with a target point marked with an asterisk) were visualized for the unicycle and the kinematic car, respectively.     Tables 2 and 3 A path is considered to be accurate if an end-point of generated trajectory x x x f real in the task-space is close enough to the target point x x x f , i.e., where x x x 0 = k k k(q q q 0 )-the starting point and η-an accuracy coefficient. In Tables 2 and 3, data were collected for η = 0.3. The most accurate path satisfies In Tables 4 and 5, statistical data (length and energy ranges for paths satisfying Condition (26)) were collected. The case when no path meets Condition (26) is marked with a dash (or not included at all when both versions of the Newton algorithm failed).      Based on results of the tests, some observations can be formulated: 1. Length and energy of resulting paths strongly depend on initial values of p p p c , especially when the basic Newton algorithm is used. The Newton algorithm, being local and iterative, has to generate trajectory-dependent results (the next iteration output depends on the output of the current iteration) and it tends to get into a desired location as fast as possible. It may happen that the basic algorithm provides better results than the optimized one. However, a smaller variance of results was observed in the energy-optimized version. In this case, the search space is searched more thoroughly.

2.
Trajectories generated with the optimized version of the Newton algorithm usually form a few clusters corresponding to local minima. This feature was not observed in the basic algorithm, as trajectories are not constrained by the energy minimization.

3.
Although the null-space optimization Newton algorithm does not necessarily guarantee the best accuracy, still much more often it promptly reaches the target. 4.
The failure ratio, cf. Equation (26), is much bigger for non-redundant parameterizations. What is interesting, in some tasks, parameterizations with the same number of parameters, cf. (23) and (24), generate significantly different results. Sometimes a task is solvable in one parameterization, while it is not in the other. 5.
The optimization within the null space significantly increases computational costs. In order to decrease the costs, this version of the Newton algorithm might be run after the target point was reached using the basic algorithm. 6.
For the unicycle robot and the identity output function, the quality of the solutions was generally better than that for tasks with smaller dimensional task-spaces, likely due to simplicity of the model and a small number of local minima. However, the length of a path and an energy of motion is smaller for tasks with low-dimensional output functions as it is less restricted. 7.
Obviously, the higher dimensional output space is, the more difficult task is generated. Additionally, a high redundancy in a parameter space is more computationally costly, but gives more flexibility as well. Thus, a compromise between flexibility and costs has to be found. 8.
Computational costs are smaller for the identity output function than for other output functions, as there is no need to transfer data between configuration and task spaces. 9.
Short length steps, characterized by small values of δ, are easier to control. 10. There are a few sources of inaccuracies in motion planning using Lie-algebraic methods: the truncation of the Lie series (5), evaluation of vector fields at q c although a trajectory surrounds the point (quite well visible for long-range motions) numeric properties of the Newton algorithm. All of them can be controlled by the right choice of parameters. 11. A planning accuracy varies substantially with the current configuration q c , as vector fields in the Lie-algebraic methods are evaluated at the point, cf. Equation (5), and vector fields exhibit different variability around different q q q c . 12. In almost all xy-projection plots one can notice characteristic cusps, which are typical in Lie-algebraic motion planning. At a close vicinity of a cusp a tangent line to the trajectory remains almost the same, while direction of motion is changed. It means that due to particular values of controls and generators, some components of velocitẏ q in Equation (1) change their signs while passing through zero. An illustrative example of such behavior can be constructed by setting: X X X long, Y Y Y short, u 2 small, and u 1 control passing though zero. Therefore, one can expect that the number of cusps strongly and positively correlates with frequencies of controls. For mobile robots tested the explanation is even simpler, as the cusps result from linear velocity sign changes. 13. Many trajectory cusps result in a zigzag motion, which is especially useful in obstaclecluttered environments as the volume of maneuvers can be controlled and kept small. However, it is also visible that generating a relatively short motion along high degree vector fields requires a relatively long path to follow. 14. Results for the kinematic car are the most interesting, especially for the identity output function. In this case, a hard to meet condition q 4 = 0 forces extensive movement. This condition is particularly difficult to satisfy for non-redundant parameterizations. Therefore, in some sub-figures of Task 4, cf. Figure 3, scales of axes were extended to present the whole paths. Nevertheless, the point x x x f real (cf. Figure 3 b,c,e,f) is close to the desired point x x x f as required. 15. It should be noted that planning a motion of the kinematic car is qualitatively more difficult than of the unicycle, as a high degree (equal to 3) of vector fields have to be involved, comparing to degree 2 for the unicycle, cf. (16), (17). 16. It may look strange that in Task 6, cf. Figure 3e, the projected path with the biggest energy is relatively short. In this case, the major contribution to the energy consumption was caused by control u u u 2 , responsible for reorientation of the robot.

Conclusions
In this paper, the Lie-algebraic method of motion planning for nonholonomic systems, originally designed to plan within the configuration space, has been extended to also work in a task-space. There are some advantages of planning within the task-space: • a smaller than the configuration space dimension either increases redundancy and facilitates optimizing a motion or allows to decrease the number of parameters of controls, • it allows to define and solve many practical tasks, especially in obstacle-cluttered environments when some components of the configuration vector are not important.
However, there are also some disadvantages: • a new source of singularities introduced by the output mapping, • the impossibility to plan only within the task-space, without referring explicitly to the configuration space (usually the task-space dimension is smaller than the configuration space dimension, thus no unique mapping between the two spaces exists).
Some general observations based on simulation results: • Depending on the planning purpose, the Newton algorithm, extensively used in the proposed algorithm, can be applied with or without optimization. The first version (low computational costs) should be used for tasks with no extra requirements or as a first try solution. • A small redundancy in setting the search space size (= dim(p p p)) is advised not to cause a significant increase of computational costs, while preserving flexibility and solvability of planning tasks.
• An initial value of the parameter vector p p p highly impacts solvability of the planning task and also the resulting trajectory shape. Therefore, the multi-start technique is recommended, especially when planning is performed in the off-line mode.
Our future research will concentrate on applications of the Lie-algebraic method in difficult environments (narrow passages, non-convex obstacles) in the task-space. In such environments, the planning should be performed through selection of sub-goals on the way to the target position. The main difficulty is the sub-goal selection is caused by the fact that nonholonomic systems are not governed by the Euclidean geometry, but rather the sub-Riemannian one [19]. The Lie-algebraic method is local and general-purpose, thus it is particularly well suited for solving tasks in partially known or dynamic environments.

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