The planning of optimal motions of non-holonomic systems

A new method to the planning of optimal motions of the non-holonomic systems is presented. It is based on a non-classical formulation of the Pontryagin Maximum Principle given in variational form, which handles efficiently various control and/or state-dependent constraints. They arise naturally due to both physical limits of the actuators of the non-holonomic systems and potential existence of obstacles in the workspace. The method proposed here provides continuous solutions in infinite-dimensional control space. It seems to be in contrast to majority of known optimization algorithms which project infinite-dimensional control space into finite-dimensional one and then apply techniques of linear and/or nonlinear programming, thus resulting only in near-optimal trajectories. Moreover, the offered control schemes do not require computation of inverse or pseudo-inverse of the Jacobian in the case of classic non-holonomic motion planning what also results in numerical stability of our approach. The performance of the proposed control strategies is illustrated through computer simulations for a chosen class of non-holonomic structures operating in both an obstacle-free workspace and a workspace including obstacles. Numerical comparison of our control scheme with the representative algorithms known from the literature is also given.


Introduction
In recent years, an interest has increased in applying the non-holonomic systems (mobile robots, tractor-trailer vehicles, orbiting satellites, space robot manipulators, etc.) to useful practical positioning tasks such as part transfer or an assembly. In order to accomplish the aforementioned tasks, the non-holonomic motion planning should be involved. It consists in obtaining (generating) a trajectory which joints a given initial posture (state) of the system with some desired final position in a workspace. For many of the repetitious processes, it is economically desirable to obtain the system motions that improve the performance of these processes with respect to some criteria. Reducing the energy required to transfer a part by the mobile robot from one location in the workspace to another constitutes a good example. In such a case, even a small energy decrease in one cycle of the process, multiplied by the total number of cycles in the lifetime of the manufacturing process, may result in essential savings in manufacturing costs. During the system motion in the workspace, some constraints are usually activated. They result from the existence of the control limits reflecting the physical abilities of the actuators and the fact that a system should not collide with obstacles. As is known, [5,57], dynamic systems with non-holonomic constraints possess some peculiar characteristics which make motion planning for them challenging. Namely, non-holonomic velocity constraints are not integrable to equivalent constraints in configuration (posture) space. Consequently, the instantaneous velocity of a system subject to nonholonomic constraints is limited to certain directions. Constraints of this kind occur for any wheeled vehicle under the no-slip condition. Manoeuvres of a system with limited instantaneous velocity directions in a work space containing, e.g. obstacles, are considerably more complex than those accomplished by holonomic system. Moreover, there exists no continuous time-invariant stabilizing feedback of a point-to-point control problem (e.g. a parking task) for drift-less nonholonomic systems [5,57]. Furthermore, the majority of the non-holonomic systems is globally controllable. However, their linearization about an equilibrium point is not controllable. Due to the aforementioned difficulties, motion planning techniques for holonomic systems are not directly applicable to non-holonomic motion planning. On account of the peculiar characteristics which possess systems subject to non-holonomic equality constraints, the solution of the motion planning problem with control and/or state-dependent inequality constraints becomes difficult from both theoretical and computational point of view. In such context, three representative approaches to the non-holonomic motion planning problem may be distinguished. The graph search-based methods including those offered in [13,28,29,33,42] generally involve some types of decomposition of the state space into cells (e.g. probabilistic road maps (PRM), rapidly exploring random trees (RRT) or uniform discrete grid representation). Then, a graph is constructed whose nodes are configurations (postures) and whose arcs are some types of paths (e.g. shortest, collision-free) connecting two configurations. The graph is searched applying the suitable algorithms. According to the discrete grid-based artificial potential field method, a global collision-free path planning of a cooperative cable parallel robot for multiple mobile cranes was presented in [58]. However, the graph search methods based on state space decomposition can be computationally both time-consuming (due to preprocessing of the continuous work space into a suitable graph) and complicated particularly for complex non-holonomic systems such as a tractor pulling the two or more trailers [11]. The second approach applying the Jacobian motion planning algorithms is based on so-called continuation method inspired by Wazewski [54]. It was introduced to robotics by Sussmann [48]. Jacobian algorithms have been developed theoretically in [8,9] and applied to non-holonomic systems in [1,11,12,26,37,[49][50][51]. The augmented motion planning techniques such as the extended Jacobian [2,4] or the configuration control proposed in [45] were adopted to non-holonomic systems in work [50]. Although motion planning algorithms based on the pseudo-inverse of the Jacobian are attractive and further investigated by many researchers, they possess some shortcomings. Namely, generated trajectories provide at most sub-optimal solutions. Moreover, pseudo-inverse motion planning strategies are not, in general, repeatable [43]. (only the extended Jacobian techniques [2,4] may generate repeatable state trajectories provided that the corresponding Jacobian matrix is non-singular.) Consequently, an important class of cyclic technological operations cannot be accomplished by the pseudo-inverse approach. Furthermore, almost all those motion planning algorithms require explicitly pseudo-inverse of the Jacobian which can potentially contain kinematic singularities. In order to tackle the singular configurations, the use of damped least squares has been proposed in works [38] in lieu of the pseudo-inverses. Nevertheless, this technique suffers from the tracking errors due to a long-term numerical integration drift. The non-holonomic planners applying the inverse of the extended Jacobian, can potentially introduce algorithmic singularities [2] which cause the motion planning to stop even though the non-holonomic system is not in a kinematic singular configuration. An iterative learning discontinuous controller using the inverse kinematics has been proposed in [40] to accomplish repetitive tasks. The third approach of motion planning offered in works [7,10,15,16,24,30,32,34,36,41,55] is based on application of optimal control methods or the calculus of variations. Nevertheless, application of Pontryagin's Maximum Principle in its classical form results in a two-point boundary value problem which is hard and sensitive to solve numerically. Moreover, maximizing the Hamiltonian (instead of directly minimizing the performance index) may result in stationary solutions which are not necessarily optimal ones. Using the calculus of variations and different kinds of parameterizations of controls, works [7,15,16,24,30,56] convert the infinite-dimensional optimal control formulation into a finite-dimensional nonlinear programming problem. The result of projection is obtaining only near-optimal solutions. The approach in [10] formulates the planning problem as a problem of the calculus of variations and handles state inequality constraints as equality ones using slack variables which significantly increases the number of unknown functions to be found. Optimal trajectory generation using spline interpolation techniques was offered in works [21,53] only for a simple class of mobile robots. This work also addresses the problem of optimal motion planning for the non-holonomic systems subject to both control and state-dependent constraints. In the case considered, it is very difficult to use Pontryagin's Maximum Principle in its classical form [39], since it presents only a positive form of control. Therefore, it is natural to make an attempt to solve the aforementioned constrained optimal control problem by resorting to other techniques. A new method based on a non-classical formulation of the Pontryagin Maximum Principle [14,18,23] given in variational form, which makes it possible to handle state inequality constraints efficiently, is proposed here to determine optimal motions of the non-holonomic systems. This approach, which consists in transforming the state constraints into control-dependent ones, in contrast to the interior penalty function method, does not require an initial admissible system trajectory (whose determination may be very troublesome in practice). In contrast to optimization algorithms known from the literature [19,20], which project infinite-dimensional control space into a finite-dimensional one and then apply techniques of liner programming problems, thus resulting only in near-optimal trajectories generated by discontinuous controls, the method proposed herein provides continuous solutions in an infinite-dimensional control space. As opposed to the most of the existing control algorithms, which provide only sub-optimal or near-optimal solutions and are based on maximizing the Hamiltonian [44,47], the control schemes offered in our study directly minimize the performance index. As is known, the Hamiltonian-based approaches provide, in fact, only stationary solutions which may not necessarily be optimal ones. Moreover, in the case of classic (constraint-free) non-holonomic motion planning problems analysed and solved in the majority of works by means of the pseudo-inverses of Jacobian matrices [1,11,12,26,37,49,50], our control scheme is significantly simplified in such a way that it does not require computation of inverse or pseudo-inverse of the Jacobian what also results in numerical stability of the approach. The technique of convergence proof for the proposed algorithms is different from those known in the literature and is based on searching for limit trajectory whose existence is theoretically ensured. (Optimal control problems with state inequality constraints may not have limit controls in a class of measurable Lebesgue mappings [17,25]). The paper is organized as follows. Section 2 presents the task of non-holonomic systems in terms of the optimal control problem. Section 3 offers new trajectory generation schemes based on negative formulation of the Pontryagin Maximum Principle to determine the optimal motions of the non-holonomic system. A numerically efficient procedure is proposed in Sect. 4 to determine an initial admissible control which is close (in the L 2 norm) to an optimal one. Section 5 provides computer examples of planning the optimal motions in a workspace without and with obstacles for both a simple benchmark non-holonomic structure and the car pulling the two trailers.

Problem formulation
In order to express the task of a non-holonomic system in terms of an optimal control problem, the state vector χ = χ(t) = (χ 1 (t), . . . , χ n (t)) T ∈ R n is introduced, where t ∈ [0, T ]; T a fixed or final moment of performing the task; n stands for the amount of state variables. Then, the dynamic equations of the non-holonomic system can be written in a general state space form as [3] where v(t) ∈ R m denotes the vector of controls; m stands for the number of controls; f : R n ×R m ×R −→ R n ; f (·, ·, ·) is a smooth mapping with respect to χ , v and t, respectively. Due to non-holonomic system (1) considered herein, the following inequality holds true: m < n. In further analysis, mapping χ(t) is assumed to be absolutely continuous with respect to time. Moreover, control v(t) = (v 1 (t), . . . , v m (t)) T is an integrable Lebesgue mapping, on which some constraints resulting from the abilities of the physical actuators of system (1), are imposed where v l = (v l,1 , . . . , v l,m ) T and v u = (v u,1 , . . . , v u,m ) T are lower and upper limits on control v(t), respectively. If the kinematic non-holonomic constraints are only taken into account when accomplishing the task, then state vector χ in (1) equals χ = q(t), where q(t) is the posture of the non-holonomic system and control vector v(t) corresponds to auxiliary velocities η(t) [6], i.e. v(t) = η(t). In addition, if both kinematic non-holonomic constraints and system dynamic equations (second-order non-holonomic constraints) are analysed by the accomplishment of the task, then state vector χ is equal to χ = (q T (t) η T (t)) T and v(t) in (1) presents torques acting in the wheels of the system. The task is to move system (1) from a given initial state χ(0) = χ 0 to a final desired system position y d in the κ-dimensional task space, where κ ≤ n. Formally, the last condition may be described as follows where g : R n −→ R κ ; g(·) denotes the output equations of system (1) which are assumed to be at least twice continuously differentiable; χ(T ) is an unknown final state to be found. In the particular case, i.e. for g(·) = id(·), where id(·) denotes the identity mapping, χ(t) is assumed to be known. During the system motion, collision avoidance constraints are induced. They result from the fact that non-holonomic system should not collide with the static and/or moving obstacles. The general form of these (inequality) constraints, in the case of moving obstacles, may be written in the following way: where c i = γ i c − r i (χ , t); γ i c is a small positive number (a safety margin around the i-th obstacle); r i (·, ·) denotes either a distance function between the nonholonomic system and the i-th obstacle [22] or an analytic description of the obstacle [31] and O stands for the total number of (moving) obstacles in the work space, in which the system operates. Functions r i are assumed to be continuously differentiable. For O = 0, there are no obstacles in the work space. The motion of the non-holonomic system should be accomplished in such a way as to minimize the following performance index: where φ(·, ·, ·) is a smooth non-negative function of vector variables χ , v and time t. The next section will present an approach that renders it possible to solve optimization problem (1)-(5) making use of the nonclassical formulation of Pontryagin's Maximum Principle.

Generation of optimal trajectory of the non-holonomic system
In what follows, system (1) is assumed to be completely controllable for the pair (χ (·), v(·)) [23]. This assumption means that both matrix ∂g(χ (T )) ∂χ is nonsingular and the rows of are linearly independent as mappings of time t, where Ψ = Ψ (t) is a solution of the fundamental differential ∂χ T Ψ ; Ψ (0) = I n ; I n denotes the n × n identity matrix. Moreover, from the Euler-Lagrange equations of non-holonomic system [3], we can obtain the following upper estimation on the norm of the mapping f (·, ·, ·): || f (χ , v, t)|| ≤ c 1 ||χ || + c 2 ||χ || 2 + c 3 ||v|| + c 4 , where c i denotes some positive coefficients; i = 1:4; c 4 is an upper estimation of the norm of the disturbance signal acting on the non-holonomic system. On account of the obvious inequalities d dt ||χ || ≤ ||χ|| and ||v|| ≤ v u , we see that d(||χ ||) c 5 +c 1 ||χ ||+c 2 ||χ || 2 ≤ dt, where c 5 = c 3 ||v u || + c 4 . Let us assume that Δ = c 2 1 − 4c 2 c 5 < 0 (the case in which Δ ≥ 0 is considered further on). Integrating both sides of the last differential inequality, we obtain the following bound on ||χ ||: ||χ || ≤ A 1 , where Similarly, in the case of Δ ≥ 0, we can obtain the bound on ||χ || in the form ||χ || ≤ A 2 , where Consequently, where A = max{A 1 , A 2 }. In further analysis, penalty functions are introduced to incorporate inequality constraints (2), (4) into our approach. Let η ex v (v i ) be an exterior penalty function defined for the i-th control variable v i , i = 1, 2, . . . , m, in the following way: where w ex,i v denotes a constant positive coefficient (the strength of exterior penalty). Let us observe that application of formula (7) does not guarantee maintenance of control constraints (2). In order to force fulfilment of strong inequalities (2), interior penalty functions should be introduced. They take the following form for the i-th control variable: where γ i is a small positive number (a safety margin for the i-th control variable); w in,i v denotes a constant positive coefficient (the strength of interior penalty). Similarly, let η ex O (c i ) be an exterior penalty function corresponding to η ex v and defined for the i-th obstacle, where w ex,i O > 0 is a constant coefficient (the strength of the exterior penalty). Besides, we additionally introduce for the numerical computation purposes an exterior penalty function η ex O,in corresponding to η in v and defined as where , respectively. In order to maintain control and state-dependent inequality constraints (2), (4) when operating nonholonomic system (1) in the work space, we modify performance index (5) as follows Let us note that penalty terms (7)-(10) can cause a slight increase of the value of performance index (5). Nevertheless, application of interior penalty function (8) leaves a free control space for compensation of modelling uncertainties in the feedback control actions of the system. Moreover, let us note that equality constraints (3) may be equivalently expressed in an integrable form as It should be noted that through Eq. (1), the state vector χ is functionally dependent on control v. In this way, the left-hand side of vector equality (12) also depends on v. For further considerations, it will be reformulated in functional form which is equivalent to the previous one. Thus, constrain (12) assumes the following (functional) form: where (11). To use the negative formulation of Pontryagin's Maximum Principle [14,18,23], we assume that a certain admissible control v 0 = v 0 (t), that is, satisfying relations (2) but not necessarily (3) and (4), is known. Such a control always exists, e.g. v 0 (t) ≡ 0. Moreover, we presume that this control does not minimize performance index (5). Next the increments of functionals given by the left-hand side of relation (13) and the right-hand side of (11) must be determined. Therefore, it is assumed that the is a given small number that guarantees the correctness of the presented method.
According to the theory of small perturbations [39], it results in a small perturbation of state trajectory; that is, the trajectory χ 0 (t) (which corresponds to the admis- with where ||o(δv)|| ||δv|| → 0 as ||δv|| → 0 and δχ(t) should be treated as the exact difference between the perturbed and unperturbed trajectory. For simplicity, in further analysis, we neglect the term o(δv), bearing in mind that in this case, variation δχ differs from the exact difference by o(δv). The value of the functional F 0 (·) for control v 0 + δv can thus be expressed as (its derivation is presented in "Appendix A" section) ∂v , δv dt is the Frechet differential of functional F 0 (·); ·, · means the scalar product of vectors in the Euclidean space; ∂ F 0 ∂v denotes the Frechet derivative of F 0 (·) equal to ; ψ 0 (·) stands for the conjugate mapping computed by solution of the following differential equation: with boundary condition ψ 0 (T ) = 0. Similarly, the value of the vector functional determined by the lefthand side of Eq. (13) for control v 0 + δv becomes (its derivation is also presented in "Appendix A" section) ; ψ 1 (·) stands for the conjugate mapping obtained by solving the following differential equation: with boundary condition ψ 1 (T ) = 0. For properly selected variation δv, Frechet's differentials of functionals (11), (13) can approximate the exact increments of these functionals with any accuracy. Hence, the negative formulation of Pontryagin's Maximum Principle (when neglecting the higher-order terms o(δv)) has the form (17) subject to the vector equality constraint On account of the fact that v 0 does not necessarily satisfy (13), the norm ||F 1 (v 0 )|| ≥ 0 can be large. Hence, in order to ensure required accuracy of the linear approximation, the aforementioned equality constraint is modified as follows is a linear mapping with respect to δv) where 0 < g 0 << 1; δv 0 is now "a small" (algebraically) variation equal to δv 0 = g 0 δv. Minimization of (17) results in variation δv 0 = −g 0 ∂ F 0 ∂v , where g 0 denotes a positive small number, which does not necessarily satisfy equality constrain (18). Therefore, we have to solve the following auxiliary control problem (projection of variation δv 0 onto linear functional constrains (18)): subject to constraint (18) Relations (19), (20) constitute the problem of the relative minima from the calculus of variations whose solution may be expressed as follows Taking into account (21), we obtain the new control v ∂v ; g k , g k are positive coefficients which fulfil the following inequalities: 0 < g k , g k << 1, and the corresponding (according to Eq. (1)) state trajectories χ k+1 are obtained as a result of solving iterative scheme (17)-(22). On account of the obvious rela- Consequently, from the obvious inequality (δv k * is the solution of optimal control problem (19), (20) for sufficiently large k) and the definition of δv k , we have Finally, we obtain that δ F 0 (v k , δv k * ) ≤ 0 what implies decrease (not monotonous for algebraically small values of k) of sequence F 0 (v k ). The natural question arises, under which conditions there exists a limit control for the sequence of v k as k → ∞. Based on the work [17], we can establish the following result for Φ 0 = φ.
where conv{·} denotes the convex hull of the set {·}; then there exists Lebesgue measurable limit control being the solution of successive minimization problems (17)- (22).
Proof The proof of Lemma 1 is similar to that given in [14]. Therefore, it is omitted.
However, in the applications of the control of nonholonomic systems, it is difficult to fulfil equality (23). Moreover, from the practical point of view, it is not essential to know a limit control (which may not exist theoretically for control problem (1)-(5) in a class of Lebesgue measurable mappings), and only the convergence of the trajectories χ k , k = 0, 1, 2, . . . is of great importance. We are now ready to offer the following Theorem.
Proof The proof of Theorem 1 will be limited to the key moments of argumentations. First, observe that inequality (6) expresses the fact of the boundedness of χ . By solving successive minimization problems (17)- and where (24) means that the set {χ k : k = 0, 1, 2, . . .} is uniformly bounded, whereas inequality (25) implies that this set is equi-continuous. Thus, the Arzela-Ascoli theorem [35] holds for χ k . A convergent subsequence χ k n , n = 1, 2, . . ., for which lim n→∞ χ k n = χ * , may be chosen from the sequence χ k . From Filippov's theorem [17], it follows that χ * is an absolutely continuous mapping. Furthermore, the continuity of expression g(χ (T )) − y d with respect to χ yields the following relation: We also know that for sufficiently large n, sequence F 0 (v k n ) monotonous decreases what implies (in general, local) optimality of χ * .
In practice, however, the fact of satisfying Cauchy's condition with a given accuracy by the components of sequence (χ k , F 0 (v k )) may be taken as a solution to optimal control problem (1)- (5). Consequently, iterative scheme (17)-(22) may be stopped provided that the following inequalities are fulfilled: ≤ c , where , c are small numbers (accuracy of iterative solution). A few remarks may be made regarding successive minimization process (17)- (22) and Theorem 1.
-Remark 1 Observe that iterative scheme (17)- (22) provides continuous solutions to variations δv k , k = 0, 1, 2, . . . which result in continuous controls v k . The optimization processes known from the literature [19,20] project infinite-dimensional control space into a finite-dimensional one and then apply techniques of a linear programming problem to find δv k what causes a drastic increase of the amount of numerical computations and results only in nearoptimal solutions. Furthermore, the algorithms from [19,20] search for the solution of control problem (1)-(5) in a class of quasi-constant (discontinuous) controls. Consequently, they are not suitable for non-holonomic systems particularly of the first order, in which the controls (system velocities) have to be continuous. Moreover, by transforming state inequality constraints (4) to equivalent scalar integral form χ , t))dt, control scheme (17)-(22) avoids the troublesome process of computing the Gateaux differentials which were used in works [14,19,20]. To numerically find Gateaux differentials, one has to solve a number of Cauchy's problems. This number depends on amount of time instances activating the state inequality constraints.
-Remark 2 Let us note that in the particular case Φ 0 = 0 and the lack of control and state-dependent inequality constraints (2), (4), we obtain a classic non-holonomic motion planning problem whose solution may be significantly simplified compared with formulas (22). For this purpose, equality constraint (3) is transformed into the scalar functional form as follows where h = 1 2 ||g(χ 0 (T )) − y d || 2 . The value for functional F(v 0 + δv) can be expressed as (its derivation is presented in "Appendix B" section) where where k = 0, 1, 2, . . . is obtained as a result of solving iterative scheme (28). On account of obvious relation δ F(v k , δv k ) ≤ 0, we obtain a sequence of inequalities Observe that expression (28) does not require computation of undesirable Jacobian pseudo-inverse what results also in numerical stability of our control scheme in neighbourhoods of singular configurations. The unconstrained non-holonomic motion planning with Φ 0 = 0 has been also analysed in works [1,11,12,49]. Although studies [26,37] involve the performance index in trajectory planning, the gradient of the Hamiltonian in [26,37] is projected onto the null-subspace of the Jacobian what does not even guarantee an instantaneous sub-optimal solution in the case of the con-flict of the kinematic tasks. The non-holonomic motion planning algorithms proposed in works [1,11,12,26,37,49] require a pseudo-inverse of the Jacobian which can potentially contain kinematic and/or algorithmic singularities (in the case of an extended Jacobian). -Remark 3 Observe that the trajectory from Theorem 1 maximizes (locally, in general) the Hamiltonian for an arbitrary Φ 0 . However, for the specific tasks encountered in practice and determined by state equality constraint (3), convex set (2) and performance index with cost function Φ 0 = ||v|| 2 , the Hamiltonian becomes a quadratic form in control v (strictly concave) provided that non-holonomic system equations (1) are affine with respect to v. Hence, it does not have local maxima in (2), which are different from the global maximum of the Hamiltonian. In other words, the approach proposed here finds the optimal control v * globally maximizing the Hamiltonian. In general, it is difficult to find globally optimal trajectory χ * for the control problem with state inequality constraints (4). This is due to the fact that the optimal control may not exist at all for such problems, or there may be many local minima of functional (11). One way to overcome this difficulty is a suitable choice of the initial control v 0 and then a search for an optimal trajectory whose existence is guaranteed by Theorem 1. The control v 0 could be chosen applying trials and errors technique by e.g. random generation of constant admissible values with condition v 0 = 0. As is easy to see, for many non-holonomic wheeled mobile robots (unicycle, differential drives, car-like wheeled mobile robots, cars with trailers), their corresponding Jacobian matrices become trivially singular for v(t) = 0, i.e. control v(t) = 0 is singular configuration for the aforementioned non-holonomic systems. Nevertheless, calculation by trial and error techniques of a nonzero admissible initial control which is relatively close to an optimal one seems to be both hard and time-consuming. Alternatively, it would be reasonable to take v 0,0 (t) = 0, t ∈ [0, T ], as the initial control (v 0,0 = 0 is admissible) particularly when applying motion planner (17)-(22) in a model predictive control mode and/or in a complex work spaces including many obstacles. Moreover, as the numerical simulations carried out in Sect. 5 show, optimal solutions corresponding to v 0,0 = 0 seem to be more energy-saving than those starting with v 0,0 = 0. As we know, iterative optimization scheme (17)- (22) converges relatively fast to an optimal solution provided that both initial control does not equal zero and it is in a sufficiently small neighbourhood of the optimal solution. The next section offers a numerically efficient algorithm providing an initial nonzero admissible control for iterative scheme (17)- (22). Moreover, this initial control is shown to be close (in the L 2 m [0, T ] norm) to the optimal solution.

Computation of an initial admissible control
Determination of the initial control for iterative procedure (17)- (22) will also be based herein on a suitable modification of the negative formulation of the Pontryagin Maximum Principle. For this purpose, we modify performance index (5) as follows where . Moreover, equality constraint (3) and collision avoidance conditions (4) may be equivalently expressed as a scalar equality constrain of the form 1 2 χ , t))dt = 0, or finally in an equivalent integrable form . For further considerations, the left-hand side of (30) will be reformulated in functional scalar form as follows where Similarly as in the previous section, we take an arbitrary admissible control v 0,0 (t) (including also v 0,0 (t) = 0) as an initial guess for optimal control problem (1)- (5). Let χ 0,0 (·) denote the trajectory corresponding to control v 0,0 with initial condition χ 0,0 (0) = χ 0 . The value of the functional F 3 for control v 0,0 + δv can thus be expressed as (its derivation is presented in "Appendix A" section) Similarly, the value of the functional determined by the left-hand side of Eq. (31) for control v 0,0 +δv becomes (its derivation is also presented in "Appendix A" section) Hence, the modified version of the negative formulation of Pontryagin's Maximum Principle takes the following form: (34) subject to the scalar linear constraint where 0 < g 0,0 << 1; δv 0 = g 0,0 δv. Minimization of (34) results in variation δv 0 = −g 0,0 ∂ F 3 ∂v , where g 0,0 is a small positive coefficient, which does not necessarily satisfy equality constrain (35). Therefore, we have to solve the following auxiliary control problem (projection of variation δv 0 onto linear functional constrain (35)): subject to constraint (35) Relations (36), (37) constitute the problem of the relative minima from the calculus of variations whose solution may be expressed as follows Taking into account (38), we obtain new control v 0, The process of minimization is then repeated for v 0,1 . The sequence of controls v 0,k+1 = v 0,k + δv k , where k = 0, 1, 2, . . .; ∂v ; g 0,k , g 0,k are positive coefficients which fulfil the following inequalities: 0 < g 0,k , g 0,k << 1, and the corresponding (according to Eq. (1)) state trajectories χ 0,k+1 are obtained as a result of solving iterative scheme (34)- (39). On account of the obvious relation δ F 4 (v 0,k , δv k ) = −g 0,k F ( v 0,k ) ≤ 0, we obtain inequalities F 4 (v 0,k ) > F 4 (v 0,k+1 ) ≥ 0 and consequently F 4 (v 0,k ) → 0 as k → ∞. If this is the case, from (39) we have ∂v . Using definition of δv k and the integral Schwartz inequality, it is easy to show implies decrease (not monotonous for algebraically small values of k) of sequence F 3 (v 0,k ). Based on iterative minimization scheme (34)- (39), we can propose the following theorem.
Proof Since functionals F 3 (·) and F 4 (·) are expressed in scalar forms, the proof of Theorem 2 becomes the particular case of the proof of Theorem 1, in which functional F 2 (·) was given in the vector form. Consequently, the proof of Theorem 2 may be omitted.
In practice, it suffices to take control v 0,k for k equal to several hundreds (two or three hundred) as the initial admissible solution v 0 (t) for optimization scheme (17)- (22). Let us note that the scalar form of functional F 4 makes it possible for many non-holonomic structures (unicycle, differential drives, car-like wheeled mobile robots, cars with trailers) to assume v 0,0 = 0 as the initial control for iterative procedure (34)- (39). If this is the case, we have Hence, ∂v v 0,0 =0 dt > 0 and formula (38) may be applied to determine δv 0 . From definition of F 4 (v), it follows that iterative scheme (34)-(39) seem to converge slower than that given by expressions (17)- (22). The reason is that F 4 (v) < ||g(χ (T )) − y d || for small norms ||g(χ (T ))− y d || provided that there are no collisions between system (1) and the obstacles. On the other hand, procedure (34)-(39) of determining the initial admissible solution should converge faster than that given by (17)- (22) for larger norms ||g(χ (T )) − y d || (>2).

Numerical examples
The aim of this section is to illustrate the efficiency of planning strategies (17)- (22), (28), (34)- (39) in the case of existence or non-existence of both control and/or state-dependent inequality constraints on a suitably chosen number of motion planning tasks. The two non-holonomic dynamic systems are assumed in this section to accomplish the corresponding pointto-point tasks. The first one is a simple benchmark non-holonomic system (unicycle or differential drive wheeled mobile robot [46,52]) moving without slip on a horizontal plane. This non-holonomic system is used as a basis for many types of non-holonomic wheeled mobile robots. For this reason, this model has attracted much theoretical attention by nonlinear systems workers. Moreover, this simple non-holonomic structure was chosen herein to compare performance of control strategy (28) with representative motion planning algorithms known from the literature. The second plant, assumed in the simulations, is a car pulling the two trailers without slip which constitutes a complex nonholonomic dynamic system subject to control and state-dependent inequality constraints. As is known, both non-holonomic dynamic systems considered herein are controllable. Nevertheless, a collision-free manoeuvring of the car with trailers in a work space including obstacles does not even seem to be intuitively obvious due to both limited instantaneous velocity directions of the car and a strongly non-convex shape of the plant structure. Consequently, a point-to-point task to be accomplished by the car with two trailers makes a challenge for control strategies (17)- (22), (28), (34)- (39). In all the simulations, SI units are used. The output equation g(·) from (3) is assumed in the whole section to be identity mapping, i.e. g(χ (T )) = χ(T ), (κ = n). Moreover, accuracy of iterative solutions is assumed to be less or equal to 10 −4 , i.e. ≤ 10 −4 and c ≤ 10 −4 , respectively.

Numerical comparisons of the proposed control scheme with the Newton algorithm
Due to some analogies between our control strategy (28) which incorporates functional derivatives into iterative process and the Newton algorithms involving the Jacobian pseudo-inverse matrices, it seems natural to numerically compare these two different approaches on a chosen non-holonomic motion planning task. This task is assumed to be accomplished by the simple benchmark non-holonomic system represented herein by either unicycle or differential drive [46,52]. On account of the fact that state constraints of the motion planning problem are directly expressed in task coordinates, the posture χ of this system will also be represented in the same coordinates. As is known, the unicycle in the strict sense is a mobile robot with serious problem of balance in static conditions. (It can fall over.) However, there exist vehicles that are kinematically equivalent in the task coordinates to the unicycle but more stable from a mechanical viewpoint. Among these, the most important is the differential drive [46,52]. This drive consists of two actuated fixed wheels mounted on the left and right side of the vehicle platform. Actuated fixed wheels are driven by motors mounted on fixed positions of the vehicle. Their axis of rotation has a fixed direction with respect to the platform's coordinate frame. The two wheels of radiuses equal to R are independently driven. For the differential drive, denote (χ 1 , χ 2 ) T ∈ R 2 the Cartesian coordinates of the midpoint of the segment of the length d joining the two wheel centres and by χ 3 ∈ R the common orientation of the actuated fixed wheels (hence, of the platform). Then, motion equations of the differential drive expressed in the task coordinates present an analytic drift-less dynamic system of the form [46] where χ = (χ 1 , χ 2 , χ 3 ) T denotes the posture of the differential drive; n = 3; v = (v 1 , v 2 ) T ; v 1 denotes linear velocity of the vehicle; and v 2 is its angular velocity (m = 2). If the wheels rotate at the same angular speed, the differential drive moves straightforward or backward (v 2 = 0). If one wheel is rotating faster than the other, the vehicle follows a curved path (v 1 = 0 and v 2 = 0). If both wheels are rotating at the same velocity in opposite direction, the differential drive turns about the midpoint of the two driving wheels, i.e. v 1 = 0. The linear and angular velocities v 1 , v 2 in (40) are expressed as functions of the physical inputs, i.e. the angular speeds ω 1 , ω 2 of the actuated fixed wheels. As is known, there is a one-to-one correspondence between v in (40) and ω = (ω 1 , ω 2 ) T given by the following expression: From linear Eq. (41), it follows that the upper and lower wheel speed limits v u , v l are determined based on the physical upper and lower angular speed limits ω u , ω l , respectively, which depend on the physical abilities of the motors driving the both wheels of the differential drive. Motion equations (40) of the differential drive are most often used as a basis for many types of non-holonomic wheeled mobile robots. Hence, the numerical comparison of the motion plan-ning algorithms will also be carried out in this section for the differential drive. There are many representative Newton algorithms of non-holonomic motion planning [1,11,12,26,37,49]. Among them, we chose an improved version of the algorithm (originated in [11]) recently proposed in the work [49]. The control law offered in [49] is given by the following iterative formulas: where , 0) T ; = I 2 0 0 0 ; γ = 0.1 and 0 3×3 denotes the 3×3 zero matrix; k = 0, 1, 2, . . .. The first task also analysed in [49] is to move the differential drive from initial posture χ(0) = (0, 0, 0) T to the final location y d = (1, 1, 0) T in the time horizon T equal to T = 2. The initial control v 0 assumed in work [49] equals v 0 = 1, sin( 2π t T ) T . Our iterative scheme (28) with the same initial control and constant gain coefficients g k equal to g k = 0.3, k = 0, 1, 2, . . . , has been applied to solve the aforementioned positioning task. The results of application of iterative schemes (28) and (42) are depicted in Figs. 1 and 2 which show accurate positioning of both schemes. Nevertheless, energy E lost during the differential drive movement equals E = 4.81 for control scheme (42) and E = 4.1 for our iterative procedure, respectively, where E = T 0 ||v|| 2 dt. Inverse-free trajectory generator (28) has also been started with g k = 0.3 and initial control v 0 = 0 which is singular for iterative procedure (42). The convergence result for this simulation is given in Fig. 3.
Moreover, the lost energy by the differential drive equals E = 3.81 in this case. Consequently, iterative scheme (28) seems to be more energy-saving than that  Fig. 3 Convergence to y d corresponding to generator (28) with v 0 = 0 given by formula (42). In the last simulation of this section, we tried to find the most energy-saving solution for the analysed positioning task. For this purpose, iterative optimization scheme (17)- (22) has been applied   (17)- (22) with v 0 = 0 and Φ 0 = ||v|| 2 with performance index (11) equal to Φ 0 = ||v|| 2 . From Remark 3, it follows that our control scheme finds global solution in such a case. Gain coefficients g k , g k , k = 0, 1, 2, . . . , are assumed to be constant for all k and equal to g k = 0.01 and g k = 0.1, respectively. Initial control v 0 is the same as that given in work [49], i.e. v 0 = 1, sin( 2π t T ) T . The results of the simulation are presented in Figs. 4 and 5 which indicate accurate positioning of the differential drive (Fig. 4) and convergence to the global minimum (Fig. 5) which equals E = 3.6.

Non-holonomic motion planning for a car pulling the two trailers
The dynamic structure considered in this section is a car pulling the two trailers. The whole system is subject to the rolling of the wheels without sliding of both the Fig. 6 Kinematic scheme of the car pulling the two trailers and the positioning task to be accomplished car and the trailers. The car will be represented by two driving wheels connected by an axle. The posture of the system is given by two position coordinates of the last (second) trailer and three angles with respect to an absolute coordinate system. There are only two inputs, namely one linear velocity and one angular velocity which represent the actions of the car. The posture χ of the whole non-holonomic system (shown schematically in Fig. 6) is expressed by five task coordinates (n = 5) χ = (χ 1 , . . . , χ 5 ) T , where (χ 1 , χ 2 ) are the coordinates of the centre of the axle between the two wheels of the second trailer; χ 5 is the orientation angle with respect to the Oχ 1 axis of the pulling car; χ i for 3 ≤ i ≤ 4 denotes the orientation angle of the trailer 5 − i with respect to the Oχ 1 axis. The motion equations of a car pulling the two trailers present dynamic system of the form [27] where l 1 is the distance from the wheels of the last trailer to the wheels of the first trailer; l 2 denotes the distance from the wheels of the first trailer to the wheels of the pulling car, 2w stands for the distance between the wheels of the car; l is the length of the car; v = (v 1 , v 2 ) T ; v 1 denotes the linear velocity of the pulling car, and v 2 is its angular velocity (m = 2). Distances between the two wheels of the first and second trailers are assumed also to be equal to 2w. Let us note that the pulling car represents physically a differential drive with one or two passive castor wheels which are attached in the front of the car for balance and stability. Castor wheels are not actuated but they can also rotate freely about an axis perpendicular to their axis of rotation. The linear and angular velocities v 1 , v 2 in Eq. (43) of the pulling car may be expressed as linear functions of the physical velocity inputs, i.e. angular speeds ω 1 , ω 2 of the both actuated fixed wheels of the differential drive. The one-to-one correspondence between v in (43) and ω = (ω 1 , ω 2 ) T is given by the following formula: The drawbar of the first trailer is for simplicity attached to the midpoint of the segment joining the two wheel centres of the car. Similarly, the drawbar of the second trailer is attached to the midpoint of the segment joining the two wheel centres of the first trailer. All the wheels of the trailers are passive (not actuated). Moreover, the drawbars of the trailers are rigidly fixed to the trailer platforms.

Constraint-free singular motion of the car with two trailers
The second task is to move the car pulling the two trailers from its initial posture χ(0) = (−4, −2, 0, 0, π/2) T to desired location y d = (−4, −2, 0, π/2, π) T . The final moment T of accomplishing the task equals T = 12. As is known [27], postures χ(0) and y d are singular for non-holonomic system (43). The initial control v 0 equals v 0 = 0. Our control scheme (28) has been applied to solve the third task with gain coefficients g k equal to g k = 0.05, k = 0, 1, . . . The results of the simulation are presented in Figs. 7, 8, 9 and 10 which indicate that accurate positioning of the car with trailers is achieved (Fig. 9). Moreover, as is seen from  Fig. 8 Algebraic difference of controls v 1 − v 2 corresponding to generator (28) with v 0 = 0 and non-trivial singular motion Figs. 7 and 8, controls generating the system movement schematically shown in Fig. 10, are both continuous and non-trivially singular (v 1 = v 2 ) according to the results of the work [27]. Consequently, iterative scheme (28) is also able to tackle non-trivial singular controls.

Constraint optimal motions of the car pulling the trailers
The third task is to move non-holonomic system (43) Fig. 10 The singular motion of the car with two trailers corresponding to generator (28) for v 0 = 0 imize performance index (11). There are two obstacles in the work space schematically shown in Fig. 6.
The threshold values γ 1 c and γ 2 c taken for computations are equal to γ 1 c = 0.14 and γ 2 c = 0.2, respectively. Let us note that both initial and desired postures of non-holonomic system (43) are located in neighbourhoods of obstacles to make difficult the collision-free manoeuvring of the car pulling two trailers. In order to calculate the values of r i which are the distance functions, both the length l of the car and its width (represented by the length 2w of the axle connecting the two driving wheels) have been discretized into 8 points. Analogous discretization process was carried out for the trailers. For simplicity of computations, all the wheels of non-holonomic system (43) are not taken into account in the collision avoidance. Depending on the initial control (nonzero or zero) and performance index Φ 0 , we have carried out two simulations with , where φ = ||v|| 2 , respectively. In this subsection, exterior penalty functions (7), (9) are assumed to be incorporated in iterative optimization schemes (17)-(22), (34)- (39). The following values for the exterior penalty function parameters w ex,i v , w ex,i O are taken in the first two simulations: w ex,i v = 1.5, w ex,i O = 7.5, i = 1, 2, respectively. The first two simulations depend on initial control. In the first simulation, initial control v 0,0 is equal to zero, i.e. v 0,0 = 0, t ∈ [0, 16]. Consequently, iterative scheme (34)-(39) from Sect. 4 has to be applied to determine nonzero initial admissible control (close to the optimal one) for optimization procedure (17)- (22). Gain coefficients g 0,k , g 0,k , g k g k , k = 0, 1, 2, . . . are assumed to be constant for all k, i.e. g 0,k = g 0,k = 0.015 and g k = g k = 0.15, respectively. In order to find nonzero initial control for optimization process (17)- (22), five hundred iterations of the algorithm given by expressions (34)-(39) were executed. Then, proper optimization process (17)- (22) has been started. The results of the first simulation are presented in Figs. 11, 12, 13, 14 and 15 which show accurate positioning of dynamic system (43) (Fig.  13).
As is seen from Figs. 11 and 12, control constraints (2) are not violated. Figure 14 presents convergence of performance index F 0 (v) to the (global) minimum which guarantees collision-free movement of the car with two trailers (Fig. 15). The energy E lost dur-   (7), (9) ing the movement equals E = 8.76. By applying the trial and error technique, we have found a nonzero initial admissible control for the second simulation. It  (7), (9) equals v 0 = 0.4, −0.2 sin( 2π t T ) T . Figure 16 shows the movement of non-holonomic system (43) corresponding to control v 0 . Observe that v 0 violates both state equality and inequality constraints (3), (4). Gain coefficients g k , g k taken in this (second) simulation equal g k = g k = 0.04, k = 0, 1, 2, . . . Iterative optimization scheme (17)- (22) has been applied to find optimal solution. Figures 17, 18, 19, 20 and 21 present the results of the second simulation. As is seen from Figs. 17, 18, 20 and 21, control and state-dependent inequality constraints are maintained in the whole time horizon of the system movement. Furthermore, Fig. 19 presents convergence of system (43) to desired location y d . The energy E required to move the car in the second simulation equals E = 9.47. Similarly as in the case of the differ-    (7), (9) and v 0 = 0 ential drive, nonzero initial admissible control provides more time-consuming optimal motion than that equal to zero.  Fig. 19 Convergence to y d corresponding to generator (17)- (22) with exterior penalty functions (7), (9) and v 0 = 0  (7), (9) and v 0 = 0 In the third simulation, exterior penalty functions (7), (9) are applied with  (7), (9), v 0 = 0 and φ = 0 is equal to that from the previous simulation. The following values for the exterior penalty function parameters w ex,i v , w ex,i O were taken in the current simulation: w ex,i v = 150, w ex,i O = 1500, respectively. Gain coefficients g k , g k are assumed to be constant for all k, k = 0, 1, . . . , and equal to g k = 0.0005, g k = 0.2, respectively. The results of the current simulation are presented in Figs. 22, 23 and 24 which indicate accurate positioning of the car with trailers (Fig. 22). As is seen from Fig. 23, performance index F 0 (v) − J (v) reflecting fulfilment of control and state inequality constraints (2), (4), converges to the (global) minimum what guarantees the collision-free movement of the car pulling the two trailers. Figure 24 presents the convergence of J (v) to the (global) minimum equal to J (v) = 7.2. Observe that energy E required to move the car with trailers in the first simulation for initial guess v 0,0 = 0  (7), (9) and v 0 = 0 differs slightly from that obtained in the current simulation and equals to E = 7.2.

Conclusions
A new method based on the non-classical (negative) formulation of the Pontryagin Maximum Principle to find optimal trajectories of non-holonomic systems is proposed in the paper. An important factor which affects the speed and accuracy of determining the minimizing sequence of controls, compared with works [19,20] , is elimination of both the time-consuming linear programming problem and discontinuity of the controls. Moreover, proposed iterative control scheme (28) does not require explicit computing any inverse or pseudo-inverse of the (general) Jacobian matrix. This important feature of the offered generator for the classic (constraint-free) non-holonomic motion planning tasks does not cause to stop the algorithm if singular configurations appear when operating the non-holonomic system in the workspace, as the computer simulations carried out in Sect. 5 have shown. A case of statedependent control constraints can be tackled by introducing a new control vector ω such thatv = ω for v l (χ ) ≤ v ≤ v u (χ ), where v l (χ ) and v u (χ ) are statedependent lower and upper limits imposed on vector v, respectively. Thus, a modified control problem is obtained with a new trajectory (χ , v) and control ω. It is important to note that the method presented here does not call for the knowledge of an initial solution satisfying constraints (3) and (4). It is only required that admissible control should satisfy relations (2). Numeri-cal simulations confirm the theoretical results obtained in Sects. 3 and 4. The problem formulation and given approach may be directly applicable to complex multiple non-holonomic systems interacting in a workspace with obstacles. This method is also useful in finding admissible trajectories (in the sense of satisfying state constraints (3) and (4)), which is often encountered in practice.

Open Access
∂v . Let us note that for the scalar form of F j , we obtain ∂v is the m-dimensional column vector).