A contact covariant approach to optimal control with applications to sub-Riemannian geometry

We discuss contact geometry naturally related with optimal control problems (and Pontryagin Maximum Principle). We explore and expand the observations of Ohsawa (Autom J IFAC 55:1–5, 2015), providing simple and elegant characterizations of normal and abnormal sub-Riemannian extremals.


Introduction
A contact interpretation of the Pontryagin Maximum Principle In a recent paper, Ohsawa [17] observed that for normal solutions of the optimal control problem on a manifold Q, the Hamiltonian evolution of the covector Λ t in T * (Q × R) considered in the Pontryagin maximum principle (PMP), projects to a well-defined contact evolution in the projectivization P(T * (Q × R)). Here, Q × R is the extended configuration space (consisting of both the configurations Q and the costs R) and P(T * (Q × R)) is equipped with a natural contact structure. Moreover, Ohsawa observed that the maximized Hamiltonian of the PMP is precisely the generating function of this contact evolution.
The above result was our basic inspiration to undertake this study. Our goal was to understand, from a geometric viewpoint, the role and origins of the above-mentioned contact structure in the PMP and to study possible limitations of the contact approach (does it work alike for abnormal solutions, etc.).
As a result we prove Theorem 3, a version of the PMP, in which the standard Hamiltonian evolution of a covector curve Λ t in T * (Q × R) along an optimal solution q(t) ∈ Q × R is substituted by a contact evolution of a curve of hyperplanes H t in T(Q × R) along this solution. (Note that the space of all hyperplanes in T(Q × R) is actually the manifold of contact elements of Q × R and can be naturally identified with P(T * (Q × R)).) It is worth mentioning that this result is valid regardless of the fact whether the solution is normal or abnormal and, moreover, the contact evolution is given by a natural contact lift of the extremal vector field (regarded as a timedependent vector field on Q × R). Finally, using the well-known relation between contact vector fields and smooth functions we were able to interpret the Pontryagin maximized Hamiltonian as a generating function of the contact evolution of H t .
It seems to us that, apart from the very recent paper of Ohsawa [17], the relation between optimal control and contact geometry has not been explored in the literature. This fact is not difficult to explain as the PMP in its Hamiltonian formulation has been very successful and as symplectic geometry is much better developed and understood than contact geometry. In our opinion, the contact approach to the PMP seems to be a promising direction of studies for at least two reasons. First of all it allows for a unified treatment of normal and abnormal solutions and, second, it seems to be closer to the actual geometric meaning of the PMP (we shall justify this statement below).
About the proof The justification of Theorem 3 is rather trivial. In fact, it is just a matter of interpretation of the classical proof of the PMP [18] (see also [13,15]). Recall that geometrically the PMP says that at each point of the optimal trajectory q(t), the cone K t ⊂ T q(t) (Q × R) approximating the reachable set can be separated, by a hyperplane H t ⊂ T q(t) (Q × R), from the direction of the decreasing cost (cf. Fig. 2). Thus, in its original sense the PMP describes the evolution of a family of hyperplanes H t (i.e., a curve in the manifold of contact elements of Q × R, identified with P(T * (Q × R))) along the optimal solution. This evolution is induced by the flow of the optimal control on Q × R. From this perspective, the only ingredient one needs to prove Theorem 3 is to show that this flow induces a contact evolution (with respect to the natural contact structure) on P(T * (Q × R)). It is worth mentioning that the covector curve Λ t ∈ T * (Q × R) from the standard formulation of the PMP is nothing else than just an alternative description of the above-mentioned curve of hyperplanes, i.e., H t = ker Λ t for each time t. Obviously, there is an ambiguity in choosing such a Λ t , which is defined up to a rescaling.
Applications From the above perspective, it is obvious that the description of the necessary conditions for optimality of the PMP in terms of H t s (the contact approach) is closer to the actual geometric meaning of the PMP as it contains the direct information about the separating hyperplanes. On the contrary, in the Hamiltonian approach this information is translated into the language of covectors (not to forget the nonuniqueness of the choice of Λ t ).
We illustrate the contact approach to the PMP by discussing its applications to the sub-Riemannian (SR) geodesic problem in Sect. 6. Recall that an SR geodesic problem on a manifold Q is an optimal control problem in which the controls parametrize trajectories tangent to a smooth distribution D ⊂ TQ and the cost of a trajectory is its length calculated via a given positively defined bi-linear form g : D × D → R (the SR metric). Actually, due to the Cauchy-Schwartz inequality, the trajectories minimizing the length are exactly those that minimize the kinetic energy and are parametrized by the arc-length. In such a setting, using some elementary geometric considerations, we were able to relate D and g with the separating hyperplanes H t (Lemma 11). In consequence, still using elementary arguments, the following two results about SR extremals were derived: -Theorem 5 completely characterizes abnormal SR extremals. It states that an absolutely continuous curve q(t) ∈ Q tangent to D is an abnormal extremal if and only if the minimal distribution along q(t) which contains D q(t) and is invariant along q(t) under the flow of the extremal vector field is of rank smaller than dim Q. As a special case (for smooth vector fields) we obtain, in Corollary 1, the following result: if the distribution spanned by the iterated Lie brackets of a given D-valued vector field X ∈ Γ (D) with all possible D-valued vector fields, i.e., ad k X (Z ) | Z ∈ Γ (D), k = 0, 1, 2, . . . is of constant rank smaller than dim Q, then the integral curves of X are abnormal SR extremals. -Theorem 6 in a similar manner (yet under an additional assumptions that the controls are normalized with respect to the SR metric g) provides a complete characterization of normal SR extremals. It states that an absolutely continuous curve q(t) ∈ Q, tangent to D, is a normal extremal if and only if it is of class C 1 with an absolutely continuous derivative and if the minimal distribution along q(t) which contains these elements of D q(t) that are g-orthogonal toq(t) and is invariant along q(t) under the flow of the extremal vector field does not contain the direction tangent to q(t) at any point. Again in the smooth case we conclude, in Corollary 2, that if for a given normalized vector field X ∈ Γ (D) the distribution spanned by the iterated Lie brackets of X with all possible D-valued vector fields g-orthogonal to X , i.e., ad k X (Z ) | Z ∈ Γ (D), g(Z , X ) = 0, k = 0, 1, 2, . . .
is of constant rank and does not contain X at any point of q(t), then the integral curves of X are normal SR extremals.
The first of the above results seems not to be present in the literature. Of course a characterization of abnormal extremals in Hamiltonian terms is well known (see, e.g., [20]) and as such has been widely used. In many particular cases (see, e.g., [12]), it allowed to obtain criteria to derive abnormal extremals similar to the smooth version of Theorem 5. The second of the above results appears in its smooth version in [16] for rank-2 distributions, and in the general version in [4] in a formulation equivalent to ours. For these reasons, we do not claim to be the first to obtain the above results (although our simple formulations using the language of flows of the optimal control seem to be new). What we believe, however, to be an added value is the simplicity of derivation of these results in our approach. Indeed, our proofs use only basic geometric tools, actually nothing more sophisticated than the definition of the flow, the derivation of the tangent space to a paraboloid, and the Gram-Schmidt algorithm.
It should be stressed that the language of flows used throughout is much more effective, and in fact simpler, than the language of Lie brackets usually applied in the study of SR extremals. Indeed, the assertions of Theorems 5 and 6 are valid for non-smooth, i.e., absolutely continuous curves and bounded measurable controls do not require any regularity assumptions (contrary to the characterization in terms of Lie brackets) and work for single trajectories (not necessary families of trajectories).
As an illustration of the above results we give a few examples. In particular, in Examples 1 and 8 we were able to provide a surprisingly easy derivation of the Riemannian geodesic equation (obtaining the equation ∇γγ = 0 from the standard Hamiltonian approach is explained in [1,20]). In Examples 3, 7, and 9, we re-discover some results of [16,22] concerning rank-2 distributions.
Organization of the paper We begin our considerations by a technical introduction in Sect. 2. Our main goal in this part is to introduce, in a rigorous way, natural differential geometric tools (Lie brackets, flows of time-dependent vector fields, distributions, etc.) in the non-smooth and time-dependent setting suitable for control theory (in general, we consider controls which are only locally bounded and measurable). Most of the results presented in this section are natural generalizations of the results well known in the smooth case. They are essentially based on the local existence and uniqueness of solutions of ODE in the sense of Caratheodory (Theorem 7). To avoid being too technical, we moved various parts of the exposition of this section (including some proofs and definitions) to the Appendix.
In Sect. 3, we briefly recall basic definitions and constructions of contact geometry. In particular, we show an elegant construction of contact vector fields (infinitesimal symmetries of contact distributions) in terms of equivalence classes of vector fields modulo the contact distribution. This construction is more fundamental than the standard one in terms of generating functions (which requires a particular choice of a contact form). It seems to us that so far it has not been presented in the literature.
In Sect. 4, we discuss in detail a natural contact structure on the projectivization of the cotangent bundle P(T * M). In particular, we construct a natural contact transformation P(F) of P(T * M) induced by a diffeomorphism F of M. Later we study an infinitesimal counterpart of this construction, i.e., a natural lift of a vector field X on M to a contact vector field C X on P(T * M).
In Sect. 5, we introduce the optimal control problem for a control system on a manifold Q and formulate the PMP in its standard version (Theorem 2). Later we sketch the standard proof of the PMP introducing the cones K t and the separating hyperplanes H t . A proper interpretation of these objects, together with our previous considerations about the geometry of P(T * M) from Sect. 4, allows us to conclude Theorems 3 and 4 which are the contact and the covariant versions of the PMP, respectively.
Finally, in the last Sect. 6, we concentrate our attention on the geometry of the cones K t and hyperplanes H t for the Riemannian and sub-Riemannian geodesic problems. The main results of that section, which characterize normal and abnormal SR extremals, were already discussed in detail in the paragraph "Applications" above.

Technical preliminaries
As indicated in the Introduction, in this paper we shall apply the language of differential geometry to optimal control theory. This requires some attention as differential geometry uses tools such as vector fields, their flows, distributions and Lie brackets which are a priori smooth, while in control theory it is natural to work with objects of lower regularity. The main technical difficulty is a rigorous introduction of the notion of the flow of a time-dependent vector field (TDVF) with the time-dependence being, in general, only measurable. A solution of this problem, provided within the framework of chronological calculus, can be found in [2]. The recent monograph [11] with a detailed discussion of regularity aspects is another exhaustive source of information about this topic.
Despite the existence of the above-mentioned excellent references, we decided to present our own explication of the notion of the flow of a TDVF. The reasons for that decision are threefold. First of all, this makes our paper self-contained. Second, we actually do not need the full machinery of [2] or [11], so we can present a simplified approach. Finally, for future purposes we need to concentrate our attention on some specific aspects (such as the transport of a distribution along an integral curve of a TDVF and the relation of this transport with the Lie bracket) which are present in neither [2] nor [11]. Our goal in this section is to give a minimal yet sufficient introduction to the above-mentioned concepts. We move technical details and rigorous proofs to the Appendix.

Time-dependent vector fields and their flows
Let M be a smooth manifold. By a timedependent vector field on M (denoted TDVF) we shall understand a family of vector fields X t ∈ X(M) parametrized by a real parameter t (the time). Every such a field defines the following non-autonomous ODE 1 on Ṁ (2.1) A technical assumption that the map (x, t) → X t (x) is Caratheodory in the sense of Definition 11 below guarantees that solutions of (2.1) (in the sense of Caratheodory) locally exist, are unique and are absolutely continuous with bounded derivatives (ACB, see the Appendix) with respect to the time t. For this reason from now on we shall restrict our attention to TDVFs X t satisfying the above assumption. We will call them Caratheodory TDVFs. In a very similar context, the notion of a Caratheodory section was introduced in the recent monograph [11]. Actually, in the language of the latter work our notion of a Caratheodory TDVF would be called a locally bounded Caratheodory vector field of class C 1 . A solution of (2.1) with the initial condition x(t 0 ) = x 0 will be denoted by x(t; t 0 , x 0 ) and called an integral curve of X t . When speaking about families of such solutions with different initial conditions it will be convenient to introduce (local) maps whenever both sides are defined.
Since X t is Caratheodory, it satisfies locally the assumptions of Theorem 7. Now the justification of Lemma 1 follows directly from the latter result. Properties (2.2) are merely a consequence of the fact that t → A tt 0 (x 0 ) is an integral curve of X t .

Definition 1
The family of local diffeomorphisms A tτ : M → M described in the above lemma will be called the time-dependent flow of X t (TD flow).
Clearly A tt 0 is a natural time-dependent analog of the notion of the flow of a vector field. This justifies the name "TD flow". It is worth noticing that, alike for the standard notion of the flow, there is a natural correspondence between TD flows and Caratheodory TDVFs. Lemma 2 Let A tτ : M → M be a family of local diffeomorphisms satisfying (2.2) and such that for each choice of x 0 ∈ M and t 0 ∈ R the map t → A tt 0 (x 0 ) is ACB. Then A tτ is a TD flow of some Caratheodory TDVF X t .
The natural candidate for such a TDVF is simply X t (x) := ∂ ∂τ τ =t A τ t (x). The remaining details are left to the reader.

Distributions along integral curves of TDVFs
In this paragraph we shall introduce basic definitions and basic properties related with distributions defined along a single ACB integral curve x(t) = x(t; t 0 , x 0 ) (with t ∈ [t 0 , t 1 ]) of a Caratheodory TDVF X t . In particular, for future purposes it will be crucial to understand the behavior of such distributions under the TD flow A tτ of X t .
be an integral curve of a Caratheodory TDVF X t . A distribution B along x(t) is a family of linear subspaces B x(t) ⊂ T x(t) M attached at each point of the considered curve. In general, the dimension of B x(t) may vary from point to point.
By an ACB section of B we will understand a vector field Z along x(t) such that Z (x(t)) ∈ B x(t) for every t ∈ [t 0 , t 1 ] and that the map t → Z (x(t)) is ACB. The space of such sections will be denoted by Γ AC B (B). A distribution B along x(t) shall be called charming if pointwise it is spanned by a finite set of elements of Γ AC B (B).
We shall say that B is A tτ -invariant (or respected by a TD flow This follows from the fact that each map A tτ is a local diffeomorphism. Let us remark that the idea behind the notion of a charming distribution is to provide a natural substitution of the notion of smoothness in the situation where a distribution is considered along a non-smooth curve. Observe namely that a restriction of a smooth vector field on M to an ACB curve x(t; t 0 , x 0 ) is a priori only an ACB vector field along x(t; t 0 , x 0 ).

Proposition 1 Charming distributions appear naturally in the following two situations:
-A restriction of a locally finitely generated smooth distribution on M to an ACB curve x(t) = x(t; t 0 , x 0 ) is charming. -Let A tτ be the TD flow of a Caratheodory TDVF X t and let B be a distribution along an integral curve it is also charming.
The justification of the above result is straightforward. Regarding the first situation it was already observed that a restriction of a smooth vector field to an ACB curve is an ACB vector field. In the second situation, the distribution B is spanned by vector By the results of Lemma 12 these fields are ACB. Given a distribution B along x(t) we can always extend it to the smallest (with respect to inclusion) distribution along x(t) containing B and respected by the TD flow A tτ along x(t). This construction will play a crucial role in geometric characterization of normal and abnormal SR extremals in Sect. 6.

trajectory of a TDVF X t . Let A tτ be the TD flow of X t and let B be a distribution along x(t). Then
is the smallest distribution along x(t) which contains B and is respected by the TD flow A tτ along x(t).
Obviously, any distribution A tτ -invariant along x(t) and containing B x(t) must contain A • (B) x(t) . The fact that the latter is indeed A tτ -invariant along x(t) follows easily from property (2.2).

Lie brackets and distributions
Constructing distributions A tτ -invariant along x(t) introduced in Proposition 2, although conceptually very simple, is not very useful from the practical point of view, as it requires calculating the TD flow A tτ . This difficulty can be overcome by passing to an infinitesimal description in terms of the Lie brackets, however, for a price of loosing some generality. In this paragraph, we shall discuss this and some related problems in detail.

Definition 3
Let X t be a Caratheodory TDVF and x(t) = x(t; t 0 , x 0 ) its integral curve. Given any smooth vector field Z ∈ X(M) we define the Lie bracket of X t and Z along x(t) by the formula i.e., we calculate the standard Lie bracket [X t , Z ] "freezing" the time t and then evaluate it at the point x(t), thus obtaining a well-defined field of vectors along x(t) (the regularity of the map t → [X t , Z ] x(t) is a separate issue that we shall discuss later).
For future purposes, we would like to extend Definition 3 to be able to calculate the bracket [X t , Z ] x(t) also for fields Z of lower regularity. That can be done, but at a price that the bracket [X t , Z ] x(t) would be defined only for almost every (a.e.) t ∈ [t 0 , t 1 ]. The details of this construction are provided below.
As a motivation recall that for M = R n , given two smooth vector fields X, Z ∈ X(R n ) (understood as maps X, Z : R n → R n ) their Lie bracket at a point x 0 equals where t → x(t) is the integral curve of X emerging from x 0 at time 0 (in particular, ∂ ∂t 0 x(t) = X (x 0 )) and s → z(s) is the integral curve of Z emerging from x 0 at time 0 (in particular, ∂ ∂s 0 z(s) = Z (x 0 )). The above formula, actually, allows to define [X, Z ] x 0 on any smooth manifold M, simply by taking it as the definition of the Lie bracket [X, Z ] x 0 in a particular local coordinate system on M. It is an easy exercise to show that [X, Z ] x 0 defined in such a way is a true geometric object (i.e., it does not depend on the particular choice of a local chart). Note that to calculate [X, Z ] x 0 we need only to know X along s → z(s) and Z along t → x(t).
Observe that to use directly the above computational definition to calculate the Lie bracket [X t , Z ] x(t) along x(t) = x(t; t 0 , x 0 ) we should use a separate integral curve of the field X t (with "frozen" time) for every t ∈ [t 0 , t 1 ], i.e., where for each t ∈ [t 0 , t 1 ] the map τ → x t (τ ) is the integral curve of the field X t emerging from the point x(t) at time τ = t and s → z(s, t) = z(s; 0, x(t)) is the integral curve of Z emerging from x(t) at s = 0, i.e., z(t, 0) = x(t) = x(t; t 0 , x 0 ) and ∂ ∂s 0 z(t, s) = Z (x(t)). Observe now that by definitionẋ t (t) = X t (x(t)) =ẋ(t) and thus (2.3) holds for x t (τ ) = x(τ ). What is more, (2.3) is well defined at a given time t ∈ [t 0 , t 1 ] also for any vector field Z on M (not necessarily smooth) such that the map τ → Z (x(τ )) is differentiable at τ = t. This observation justifies the following statement.

Proposition 3
Assuming that t → Z (x(t)) is an ACB map and that X t is a Caratheodory

is completely determined by the values of Z along x(t) and by the values of X t in a neighborhood of x(t).
In other words, formula (2.3) is an extension of Definition 3 which allows to calculate the Lie bracket [X t , Z ] x(t) at almost every point of a given integral curve x(t) of X t , for vector fields Z defined only along x(t) and such that t → Z (x(t)) is ACB. The latter generalization is necessary in control theory, since, as t → x(t) is in general ACB only, even if Z is a smooth vector field, we cannot expect the map t → Z (x(t)) to be of regularity higher than ACB.
The above construction of the Lie bracket [X t , Z ] x(t) allows to introduce the following natural construction.

Definition 4
Let X t be a Caratheodory TDVF, x(t) = x(t; t 0 , x 0 ) (with t ∈ [t 0 , t 1 ]) its integral curve and let B be a distribution along x(t). By [X t , B] we shall understand the distribution along x(t) generated by the Lie brackets of X t and all ACB sections of B: where we consider [X t , Y ] x(t) at all points where it makes sense, i.e., at which the bracket Note that neither [X t , B] nor B +[X t , B] need be charming distributions along x(t) even if so was B as, in general, there is no guarantee that these distributions will be spanned by ACB sections (we can loose regularity when calculating the Lie bracket).
The following result explains the relation between the A tτ -and X t -invariance of distributions along x(t).
, an integral curve of a Caratheodory TDVF X t , and let A tτ be the TD flow of X t . The following conditions are equivalent:

(a) B is respected by the TD flow A tτ of X t along x(t). (b) B is a charming distribution X t -invariant and of constant rank along x(t).
The proof is given in the Appendix. Note that the equivalence between X t -and A tτinvariance is valid only if the considered distribution B along x(t) satisfies regularity conditions: it has to be charming and of constant rank along x(t).
Given a charming distribution B along x(t), it is clear in the light of the above result, that A • (B) x(t) , the smallest distribution A tτ -invariant along x(t) and containing B, should be closed under the operation [X t , ·]. Thus, in the smooth case, it is natural to try to construct A • (B) in the following way. Lemma 3 Let X be a C ∞ -smooth vector field and let B a C ∞ -smooth distribution on M. Assume that along an integral curve the distribution spanned by the iterated Lie brackets of X with all possible B-valued vector fields, i.e., is of constant rank along x(t). Then ad ∞ X (B) x(t) is the smallest distribution along x(t) containing B x(t) and respected by A t , the flow of X , i.e., ad ∞ . Proof The justification of the above result is quite simple. By construction, ad ∞ X (B) x(t) is the smallest distribution along x(t) containing B x(t) and closed under the operation ad X = [X, ·]. It is clear that ad ∞ X (B) is spanned by a finite number of smooth vector fields of the form ad k X (Z ), where Z ∈ Γ (B), and thus it is charming. Since it is also of constant rank along x(t) we can use Theorem 1 (for a time-independent vector field X ) to prove that ad ∞ On the other hand, since A • (B) x(t) is A t -invariant along x(t), again by Theorem 1, it must be closed with respect to the operation [X, ·]. In particular, it must contain the smallest distribution along x(t) containing B x(t) and closed under the operation [X, ·]. Thus, . This ends the proof.
Remark 1 Let us remark that the construction provided by (2.4) would be, in general, not possible in all non-smooth cases. The basic reason is that the Lie bracket defined by (2.3) is of regularity lower than the initial vector fields, i.e., [X t , Z ] may not be ACB along x(t) even if so were X t and Z . Thus, by adding the iterated Lie brackets to the initial distribution B, we may loose the property that it is charming (cf. also a remark following Definition 4) which is essential for Theorem 1 to hold. Also the constant rank condition is important, as otherwise the correspondence between X t -and A tτ -invariance provided by Theorem 1 does not hold. If (2.4) is not of constant rank along x(t) we may only say that ad ∞ Remark 10).
It is worth noticing that this situation resembles the well-known results of Sussmann [19] concerning the integrability of distributions: being closed under the Lie bracket is not sufficient for integrability, as the invariance with respect to the flows of distributionvalued vector fields is also needed. After adding an extra assumption that the rank of the distribution is constant, the latter condition can be relaxed.
By the results of Proposition 3, the property that a distribution B is X t -invariant along x(t) depends not only on B and the values of a Caratheodory TDVF X t along x(t), but also on the values of X t in a neighborhood of that integral curve. It turns out, however, that in a class of natural situations the knowledge of X t along x(t) suffices for checking the X t -invariance.

for every t. Then the property of B being X t -invariant along x(t) depends only on the values of X t along x(t).
The proof is given in the Appendix.

The basics of contact geometry
Contact manifolds and contact transformations In this section, we shall recall basic facts from contact geometry. A contact structure on a manifold M is a smooth co-rank one distribution C ⊂ TM satisfying a certain maximum non-degeneracy condition.
To formalize that condition we introduce the following geometric construction. From now on we shall assume that the pair (M, C) consists of a smooth manifold M and a smooth co-rank one distribution C on M. Sometimes it will be convenient to treat C as a vector subbundle of TM.
Given (M, C) one can define the bundle normal to C in TM as the quotient NC := TM/C.
Note that NC has a natural structure of a line bundle (since C is of co-rank one) over M. We shall denote this bundle by τ : NC → M.
Let now X and Y be two C-valued vector fields on M. It is easy to check that the class of their Lie bracket [X, Y ] in NC is tensorial with respect to both X and Y . That is, for any pair of smooth functions φ, ψ ∈ C ∞ (M) It follows that the assignment defines an NC-valued 2-form β : Λ 2 C → NC . Now we are ready to state the following Definition 5 A pair (M, C) consisting of a smooth manifold M and a smooth co-rank one distribution C ⊂ TM is called a contact manifold if the associated NC-valued 2-form β is non-degenerate, i.e., if β(X, ·) ≡ 0 implies X ≡ 0.
Sometimes we call C a contact structure or a contact distribution on M.
Observe that C is necessarily of even rank (M is odd-dimensional). This follows from a simple fact from linear algebra that every skew-symmetric 2-form on an odddimensional space has a non-trivial kernel.
where TF stands for the tangent map of F, is called a contact transformation. By a contact vector field (CVF) on M (or an infinitesimal symmetry of (M, C)) we shall understand a smooth vector field X ∈ X(M) preserving the contact distribution C, i.e., Equivalently, X is a CVF if and only if its (local) flow A t consists of contact transformations (cf. Theorem 1).
It is worth mentioning that the above relation between contact vector fields and flows consisting of contact transformations can be generalized to the context of TDVFs and TD flows (cf. Sect. 2). We will need this generalized relation in Sect. 5 after introducing control systems.

and only if the TD flow A tτ consists of contact transformations.
The proof follows directly from Theorem 1 by taking B = C (which is charming, see-cf. Proposition 1).

Characterization of CVFs
It turns out that there is a one-to-one correspondence between CVFs on M and sections of the normal bundle NC.

is well defined and establishes a one-to-one correspondence between sections of NC and CVFs on M.
Remark 2 Throughout we will denote by X ∈ X(M) vector fields on M, by Y ∈ Γ (C) vector fields valued in C and by C (also with variants, like C φ , C [X ] or C φ ) contact vector fields.
Proof Let us begin with introducing the following geometric construction. With every smooth vector field X ∈ X(M) one can associate an NC-valued 1-form α X : C → NC defined by the formula where Y is a C-valued vector field. The correctness of this definition follows from the fact that for every C-valued vector field Y and for any function Using the one-form α X we can prove that [X ] → C [X ] is a well-defined map, i.e., that the value of C [X ] does not depend on the choice of the representative X . Indeed, we can interpret h(X ) as the unique (note that β is non-degenerate) solution of the equation α X (·) = β(h(X ), ·). Now observe that if X and X are two different representatives of [X ], then Y := X − X is a C-valued vector field on M. Thus, we have α Y (·) = −β(Y, ·) and hence, using the obvious linearity of α X with respect to X , we get We conclude that h(X ) = h(X ) − Y and, consequently, Second, observe that C [X ] is a CVF. Indeed, by construction, given any C-valued vector field Y we have Finally, we need to check that every CVF is of the form C [X ] . By construction the class of C [X ] in NC is equal to the class of X in NC (these two vector fields differ by a C-valued vector field h(X )). Thus, the classes of CVFs of the form C [X ] realize every possible section of NC. Now it is enough to observe that the NC-class uniquely determines a CVF. Indeed, if C and C are two CVFs belonging to the same class in NC, then their difference X − X is a C-valued CVF, i.e., [C − C , Y ] ≡ 0 mod C for any C-valued vector field Y . That is, β(C − C , ·) ≡ 0 and from the non-degeneracy of β we conclude that C − C ≡ 0. This ends the proof.

Remark 3
It is natural to call a vector field X ∈ X(M) (or its NC-class [X ]) a generator of the CVF C [X ] . Observe that the NC-class of the CVF C [X ] is the same as the class of its generator X (they differ by a C-valued vector field h(X )).
In the literature, see e.g., [14], a contact distribution C on M is often presented as the kernel of a certain 1-form ω ∈ Λ 1 (M) (such an ω is then called a contact form). In the language of ω, the maximum non-degeneracy condition can be expressed as the non-degeneracy of the 2-form d ω on C. The latter is equivalent to the condition that ω ∧ (d ω) ∧n , where n = 1 Choose a contact form ω such that C = ker ω, then this correspondence is given by an A function φ is usually called the generating function of the corresponding CVF C φ associated with the contact form ω. Notice that given a contact vector field C = C φ and a contact form ω one can recover the generating function simply by evaluating ω on C, i.e., φ = ω(C).
It is interesting to relate the construction φ → C φ with the construction [X ] → C [X ] given above. Namely, the choice of a contact form ω allows to introduce a vector field R ∈ X(M) (known as the Reeb vector field) defined uniquely by the conditions ω(R) = 1 and R d ω = 0. Since R is not contained in C = ker ω, its class [R] establishes a basis of the normal bundle NC. Consequently, we can identify smooth The details are left to the reader.
Note, however, that the description of the contact distribution C in terms of a contact form ω is, in general, non-canonical (as every rescaling of ω by a nowhere-vanishing function gives the same kernel C) and valid only locally (as there clearly exist contact distributions which cannot be globally presented as kernels of single 1-forms). For this reason, the description of a contact manifold (M, C) in terms of C and related objects (e.g., NC, β) is more fundamental and often conceptually simpler (for example, in the description of CVFs) than the one in terms of ω. Not to mention that, for instance, the construction of the CVF C φ does depend on the particular choice of ω, whereas the construction of C [X ] is universal.

Remark 4
In the context of CVFs on a contact manifold (M, C) it is worth noticing an elegant correspondence between CVFs on M and a certain class of control-affine systems on M. A control-affine system on a manifold M is usually understood as a differential equation of the forṁ . . , g m . Note that the distribution G, the linear part of A, is well defined, whereas the drift is defined only relative to G, i.e., f In the light of our considerations about CVFs it is easy to prove the following fact. Indeed, to every CVF C, we attach the affine distribution (control-affine system) A = C +C. Conversely, given an affine distribution (control-affine system) A = X +C, there exists a unique CVF C ∈ Γ (A), namely C = C [X ] , such that A = C + C = C [X ] + C. In other words, on every contact manifold (M, C), there are as many CVF's C as control-affine systems A = X + C, the correspondence being established by the

Contact geometry of P(T * M)
In this section, we shall describe the natural contact structure on P(T * M) and its relation with the canonical symplectic structure on T * M (see, e.g., [5] or [14]). Later it will turn out that this structure for M = Q × R plays the crucial role in optimal control theory.
In other words, C [θ] consists of all vectors in T [θ] P(T * M) which project, under Tπ , to the hyperplane H [θ] = ker θ , see Fig. 1. We shall refer to this contact structure as the canonical contact structure on P(T * M).
The fact that (4.1) defines a contact structure is well known in the literature. The proof is given, for instance, in Appendix 4 of the book of Arnold [5], where the reasoning is based on the properties of the Liouville 1-form M on the cotangent manifold T * M. For convenience of our future considerations in Sect. 5 we shall, however, present a separate proof quite similar to the one of Arnold.
Proof Let R ∈ X(M) be any smooth vector field. We shall now construct a contact form ω R on an open subset of P(T * M) for which R will be the Reeb vector field (cf.
To finish the proof it is enough to check that ω R satisfies the maximum nondegeneracy condition. This can be easily seen by introducing local coordinates Obviously, the pull-back functions x i := (i R ) * x i with i = 0, . . . , n and p j := (i R ) * p j with j = 1, . . . , n form a coordinate system in U R . In these coordinates, the form is a matter of a simple calculation to check that such a one-form satisfies the maximum non-degeneracy condition. We conclude that ω R is, indeed, a contact form on U R for the canonical contact structure on P(T * M).

Contact transformations of P(T * M) induced by diffeomorphisms
In this paragraph, we will define contact transformations of P(T * M) which are natural lifts of diffeomorphisms of the base M.
It is worth noticing that P(F) respects the fiber bundle structure of π : P(T * M) → M: We claim that Lemma 7 P(F) is a contact transformation with respect to the canonical contact structure on P(T * M).
, which ends the proof.
Let us remark that an alternative way to prove the above result is to show that (F −1 ) * maps the contact form ω R to ω TF(R) . To prove that, one uses the fact that the pullback (F −1 ) * preserves the Liouville form.

CVFs on P(T * M) induced by base vector fields
The results of the previous paragraph have their natural infinitesimal version. Definition 8 Let X ∈ X(M) be a smooth vector field. By the contact lift of X we shall understand the contact vector field C X on P(T * M) whose flow is P(A t ), the contact lift of the flow A t of X .
The correctness of the above definition is a consequence of a simple observation that the contact lift preserves the composition of maps, i.e., P(F • G) = P(F) • P(G) for any pair of maps F, G : M → M. It follows that the contact lift of the flow A t is a flow of contact transformations of P(T * M) and as such it must correspond to some contact vector field (cf. Definition 6).
An analogous reasoning shows that given a Caratheodory TDVF X t ∈ X(M) and the related TD flow A tτ : M → M, the contact lift of the latter, i.e., P(A tτ ), will consist of contact transformations and will satisfy all the properties of the TD flow. By the results of Proposition 4, P(A tτ ) is a TD flow associated with some contact TDVF (see also Lemma 2). Obviously, this field is just C X t . The justification of this fact is left for the reader.

Lemma 8
The CVF C X is generated (in the sense of Lemma 5) by the NC-class of X, where X is any smooth vector field on P(T * M) which projects to X under Tπ , i.e., Proof Since P(A t ), the flow of C X , projects under π to A t , the flow of X , we conclude that X = Tπ(C X ). As we already know from the proof of Lemma 5, a CVF is uniquely determined by its class in NC. By (4.1), the NC-class of a field Y ∈ X(P(T * M)) is completely determined by its Tπ -projection. In other words, if two fields Y and Y have the same Tπ -projections, then Y − Y is a C-valued vector field. Thus, the field X has the same NC-class as the CVF C X so, by the results of Lemma 5 (see also Remark 3), it follows C X = C [ X ] .

Remark 5
We shall end our considerations about the contact lift C X by discussing its description in terms of generating functions (cf. our comments following Remark 3). Let us choose a vector field R on M and fix a contact form ω R = (i R ) * M on U R ⊂ P(T * M). Using the results of our previous Sect. 3 and with the help of the contact form ω R , the CVF C X can be presented as C φ for some generating function φ : U R → R. This function is simply the evaluation of ω R at C X . In fact, i.e., θ(R) = 1. In other words, the value of the generating function of C X on the class [θ ] equals the value of the representative θ , defined by θ(R) = 1, on the field X which is being lifted.

Remark 6
It is worth mentioning the following illustrative picture which was pointed to us by Janusz Grabowski. Every contact structure on a manifold N can be viewed as a homogeneous symplectic structure on some principal G L(1, R)-bundle over N . In the case of the canonical contact structure on N = P(T * M) the corresponding bundle is simply T * 0 M, the cotangent bundle of M with the zero section removed, equipped with the natural action of R\{0} = G L(1, R) being the restriction of the multiplication by reals on T * M. The canonical symplectic structure is obviously homogeneous with respect to this action. Now every homogeneous symplectic dynamics on T * 0 M reduces to contact dynamics on P(T * M). For more information on this approach the reader should consult [7,9].

The Pontryagin maximum principle
The Pontryagin Maximum Principle A control system on a manifold Q is constituted by a family of vector fields f : Q × U → TQ parametrized by a topological space U . It can be understood as a parameter-dependent differential equatioṅ For a given measurable and locally bounded u(t) ∈ U , the solution q(t) of (CS) is usually called a trajectory of (CS) associated with the control u(t).
An introduction of a cost function L : Q ×U → R allows to consider the following optimal control problem (OCP) The minimization is performed over u(t)s which are locally bounded and measurable, the time interval [0, T ] is fixed and we are considering fixed-end-points boundary conditions q(0) and q(T ). By a solution of the optimal control problem we shall understand a pair (q(t), u(t)) satisfying (OCP). Let now q(t) ∈ Q be the trajectory of the (CS) associated with a given control u(t) ∈ U . It is convenient to regard the related trajectory q(t) = (q(t), q 0 (t)) in the extended configuration space Q := Q × R, where q 0 (t) := t 0 L(q(s), u(s)) d s is the cost of the trajectory at time t. 2 In fact, q(t) is a trajectory (associated with the same control u(t)) of the following extension of (CS): Here, we treat both f and L as maps from Q × U invariant in the R-direction in Q = Q × R. In other words, we extended (CS) by incorporating the costs q 0 (t) as additional configurations of the system. The evolution of these additional configurations is governed by the cost function L. Note that the total cost of the trajectory q(t) with t ∈ [0, T ] is precisely q 0 (T ). Since the latter is fully determined by the pair (q(t), u(t)), it is natural to regard the extended pair (q(t), u(t)) rather than (q(t), u(t)) as a solution of (OCP). Note that the extended configuration space Q = Q × R q = (q, q 0 ) is equipped with the canonical vector field ∂ q 0 := (0, ∂ q 0 ) ∈ T Q = TQ × TR. We shall denote the distribution spanned by this field by R ⊂ T Q. The ray R − q := {−r ∂ q 0 | r ∈ R + } ⊂ R q ⊂ T q Q contained in this distribution will be called the direction of the decreasing cost at q ∈ Q.
Regarding technical assumptions, following [18], we shall assume that U is a subset of an Euclidean space, f (q, u) and L(q, u) are differentiable with respect to the first variable and, moreover, f (q, u), L(q, u), ∂ f ∂q (q, u) and ∂ L ∂q (q, u) are continuous as functions of (q, u). In the light of Theorem 7 below it is clear that these conditions guarantee that, for any choice of a bounded measurable control u(t) and any initial condition q(0), equation (CS) has a unique (Caratheodory) solution defined in a neighborhood of 0. It will be convenient to denote the TDVFs q → f (q, u(t)) and q → f (q, u(t)) related with such a control u(t) by f u(t) and f u(t) , respectively. In the language of Sect. 2, technical assumptions considered above guarantee that f u(t) and f u(t) are Caratheodory TDVFs. In particular, their TD flows F tτ : Q → Q and F tτ : Q → Q, respectively, are well-defined families of (local) diffeomorphisms. 3 Note that if q(t) with t ∈ [0, T ] is a solution of (CS), then for every t, τ ∈ [0, T ] the map F tτ is well defined in a neighborhood of q(τ ).
In the above setting, necessary conditions for the optimality of (q(t), u(t)) are formulated in the following PMP Theorem 2 ([18]) Let (q(t), u(t)) be a solution of the (OCP). Then for each t ∈ [0, T ] there exists a non-zero covector λ(t) ∈ T * q(t) Q such that the curve Λ t = (q(t), λ(t)) satisfies the time-dependent Hamilton equatioṅ
Proof of the PMP Although the PMP is a commonly known result, for future purposes it will be convenient to sketch its original proof following [18].
Let (q(t), u(t)) be a trajectory of (CS). By F tτ : Q → Q, where 0 ≤ τ ≤ t ≤ T , denote the TD flow on Q of the Caratheodory TDVF f u(t) defined by the control u(t) (cf. Definition 1). In other words, given a point q ∈ Q, the curve t → F t0 (q) is the a trajectory of (CS) associated with the control u(t) with the initial condition q(0) = q.
The crucial step in the proof of the PMP is introducing the, so-called, needle variations and the resulting construction of a family of sets 4 where 0 < t 1 ≤ t 2 ≤ · · · ≤ t k ≤ t < T is any finite sequence of regular points (see the Appendix) of the control u(·), v i are arbitrary elements in U and δt i are arbitrary non-negative numbers. It is easy to see that K t is a closed and convex cone in T q(t) Q, well defined for each regular point t ∈ (0, T ) of the control u(·). What is more, the cones K t are ordered by the TD flow F tτ , i.e., for each pair of regular points 0 < τ ≤ t < T . The above property allows to extend the construction of K t to non-regular points of u(·) (including the end-point T ) by setting Clearly, these sets preserve all important features of K t s, i.e., they are closed and convex cones satisfying (5.5) for any pair of points 0 < τ ≤ t ≤ T . The importance of the construction of the cone K t lies in the fact that it approximates the reachable set of the control system (CS) at the point q(t). In particular, it was proved in [18] that if at any point t ∈ [0, T ], the interior of the cone K t contains the direction of the decreasing cost R − q(t) , then the trajectory t → q(t), t ∈ [0, T ], cannot be optimal.

Lemma 9 ([18])
If, for any 0 < t ≤ T , the ray R − q(t) lies in the interior of K t , then (q(t), u(t)) cannot be a solution of (OCP).
As a direct corollary, using basic facts about separation of convex sets, one obtains the following Proposition 6 ( [18]) Assume that (q(t), u(t)) is a solution of (OCP). Then for each t ∈ (0, T ] there exists a hyperplane H t ⊂ T q(t) Q separating the convex cone K t from the ray R − q(t) . Thus, one can choose a curve of hyperplanes 5 H t ⊂ T q(t) Q separating the cone K t from the ray R − q(t) for each t ∈ (0, T ]. Because of (5.5) and the fact that the canonical vector field ∂ q 0 is invariant under T F tτ (the control does not depend on the cost), we may choose H t in such a way that for each 0 ≤ τ ≤ t ≤ T . Indeed, the basic idea is to choose any H T separating K T from R − q(T ) and define H t := T F −1 T t (H T ) for every t ∈ [0, T ]. We leave the reader to check that such a construction has the desired properties.
The geometry of this situation is depicted in Fig. 2.

Remark 7
Trajectories of (CS) satisfying the above necessary conditions for optimality (i.e., the existence of a curve of separating hyperplanes H t which satisfies (5.6)) can be classified according to the relative position of the hyperplanes H t and the line field R ⊂ T Q. Note that since the hyperplanes H t evolve according to the TD flow F tτ of the TDVF f u(t) , which leaves the distribution R invariant, we conclude that whenever R q(τ ) ⊂ H τ at a particular point τ ∈ [0, T ], then R q(t) ⊂ H t for every t ∈ [0, T ]. We call a trajectory q(t) of (CS) satisfying the above necessary conditions for optimality: normal if R q(t) ⊂ H t for any t ∈ [0, T ]. Note that, in consequence, the ray R −

q(t)
can be strictly separated from the cone K t for each t ∈ [0, T ]; abnormal if R q(t) ⊂ H t for each t ∈ [0, T ]; strictly abnormal if for some t ∈ [0, T ] the ray R − q(t) cannot be strictly separated from the cone K t (and thus R q(t) ⊂ H t for each t ∈ [0, T ]).

Consequently, R q(t) ⊂ H t for each t ∈ [0, T ]
It is worth emphasizing that being normal or abnormal is not a property of a trajectory itself, but of a trajectory together with a particular curve of separating hyperplanes. Thus, a priori, a given trajectory q(t) may admit two different curves of separating hyperplanes, one being normal and the other abnormal. On the other hand, if q(t) is a strictly abnormal trajectory it must be abnormal (and cannot be normal) for any possible choice of the curve of separating hyperplanes. To justify this statement observe that if the ray R − q(t) cannot be strictly separated from the cone K t , then necessarily (since the cones K t are closed) we have −∂ q 0 ∈ K t for some t ∈ (0, T ]. Consequently, also −∂ q 0 ∈ H t , as H t separates −∂ q 0 ∈ K t and −∂ q 0 ∈ R − q(t) (see Fig. 3). Note that since H t is a linear space, the whole line R q(t) spanned by the vector −∂ q 0 is contained in H t in this case.
It is precisely only now when the covector λ(t) of the PMP appears. Namely, one can represent each hyperplane H t ⊂ T q(t) Q as the kernel of a covector λ(t) ∈ T * q(t) Q. Due to (5.6) it is possible to choose these covectors in such a way that for every 0 ≤ τ ≤ t ≤ T the curve Λ t = (q(t), λ(t)) satisfies This reads as the Hamilton Eq. (5.1) in Theorem 2. Finally, the maximum principle (5.3) follows from the fact that Λ t can be chosen in such a way that λ(t), −∂ q 0 ≥ 0 ≥ λ(t), k for any k ∈ K t . The latter inequality for k = f v (q(t)) − f u(t) (q(t)) regarded for each v ∈ U implies (5.3).
Note that, as we have already observed in Remark 7, for abnormal solutions, we have ∂ q 0 ∈ H t = ker λ(t), and thus λ(t), ∂ q 0 ≡ 0. For normal solutions it is possible to choose λ(t) in such a way that λ(t), −∂ q 0 ≡ 1 along the optimal solution.
The contact formulation of the PMP Expressing the essential geometric information of the PMP (see Fig. 2) in terms of hyperplanes H t , instead of covectors λ(t), combined with our considerations about the canonical contact structure on P(T * M) (see Sect. 4) allows to formulate the following contact version of the PMP.

Theorem 3 (the PMP, a contact version) Let (q(t), u(t)) be a solution of the (OCP). Then for each t ∈ [0, T ] there exists a hyperplane H t ∈ P(T * q(t) Q) such that the curve t → H t satisfies the equatioṅ
where C f u(t) denotes the contact TDVF on P(T * Q), being the contact lift of the TDVF f u(t) on Q (see Definition 8 and Lemma 8).
Moreover, each H t separates the convex cone K t defined by (5.4) from the ray R − q(t) .
Proof The family of hyperplanes H t separating the cone K t from the ray R − q(t) and satisfying (5.6) was constructed in the course of the proof of Theorem 2 sketched in the previous paragraph. To end the proof it is enough to check that H t evolves according to (5.7). From (5.6) and Definition 7 of the contact lift we know that H t evolves according to P(F τ t ). By the remark following Definition 8 this is precisely the TD flow induced by the TDVF C f u(t) .
Let us remark that the contact dynamics (5.7) are valid regardless of the fact whether the considered solution is normal or abnormal. We have a unique contact TDVF C f u(t) on P(T * Q) governing the dynamics of the separating hyperplanes H t . The difference between normal and abnormal solutions lies in the relative position of the hyperplanes H t with respect to the canonical vector field −∂ q 0 on Q.
Remark 8 Actually, the fact that the evolution of H t is contact (and at the same time that the evolution of Λ t is Hamiltonian) is in a sense "accidental". Namely, it is merely a natural contact (Hamiltonian) evolution induced on P(T * Q) (on T * Q) by the TD flow on Q defined by means of the extremal vector field. In the Hamiltonian case this was, of course, already observed-see, e.g., Chapter 12 in [2]. Thus, it is perhaps more proper to speak rather about covariant (in terms of hyperplanes) and contravariant (in terms of covectors) formulations of the PMP, than about its contact and Hamiltonian versions. It may seem that the choice between one of these two approaches is a matter of a personal taste, yet obviously the covariant formulation is closer to the original geometric meaning of the PMP, as it contains a direct information about the separating hyperplanes, contrary to the contravariant version where this information is translated to the language of covectors (not to forget the non-uniqueness of the choice of Λ t ). In the next Sect. 6 we shall show a few applications of the covariant approach to the sub-Riemannian geometry. Expressing the optimality in the language of hyperplanes allows to see a direct relation between abnormal extremals and special directions in the constraint distribution. It also provides an elegant geometric characterization of normal extremals.
Although Eq. (5.7) has a very clear geometric interpretation it is more convenient to avoid, in applications, calculating the contact lift. Combining (5.6) with Theorem 1 allows to substitute equation (5.7) by a simple condition involving the Lie bracket.

Theorem 4 (the PMP, a covariant version) Let (q(t), u(t)) be a solution of (OCP). Then for each t ∈ [0, T ] there exists a hyperplane H t ∈ P(T * q(t) Q) such that the curve t → H t satisfies the Eq. (5.6). Equivalently, H t is a curve of hyperplanes which is a charming distribution that is f u(t) -invariant along q(t), i.e.,
Moreover, each H t separates the convex cone K t defined by (5.4) from the ray R − q(t) .
Proof The proof is immediate. The existence of separating hyperplanes H t satisfying (5.6) was already proved in the course of this section. The only part that needs some attention is the justification of equation ( Finally, let us discuss the description of the contact dynamics (5.7) in terms of natural contact forms introduced in the proof of Lemma 6. Recall that a choice of a vector field R ∈ X( Q) allows to define a natural embedding of the open set U R = {[θ] | θ(R) = 0} ⊂ P(T * Q) into T * Q (recall that in the language of hyperplanes, the set U R consists of those hyperplanes in T Q which are transversal to the field R). What is more, the pull-back ω R of the Liouville form Q on T * Q, is a contact form on U R . By the comment of Remark 5, the generating function of the CVF C f u(t) associated with the contact form ω R is simply where λ is a representative of the class [λ] such that λ, R = 1. In other words, the generating function of the contact dynamics (5.7) associated with ω R is precisely the linear Hamiltonian (5.2).
In particular, by choosing R = ∂ q 0 we can easily recover the results of [17]. Note that R = ∂ q 0 is the canonical choice of a vector field transversal to all hyperplanes H t 's in the normal case (note that additionally R = ∂ q 0 is F tτ -invariant). For such a choice of R, the corresponding embedding i R : U R → T * Q is constructed simply by setting λ, ∂ q 0 = 1, which is just the standard normalization of the normal solution. The associated contact form is ω R = π * 2 d q 0 + π * 1 Q , where Q is the Liouville form on T * Q and π 1 : Q × R → Q and π 2 : Q × R → R are natural projections. As observed above, the generating function of the contact dynamics associated with ω R is the linear Hamiltonian (5.2). This stays in a perfect agreement with the results of Sect. 2 in [17].
For the abnormal case there is no canonical choice of the field R transversal to the separating planes. Yet locally such a choice (but not canonical) is possible. The resulting generating function of the contact dynamics (5.7) is again the linear Hamiltonian (5.2).

Applications to the sub-Riemannian geodesic problem
In this section, we shall apply our covariant approach to the PMP (cf. Remark 8) to concrete problems of optimal control. We shall concentrate our attention on the SR geodesic problem on a manifold Q. Our main idea is to extract, from the geometry of the cone K t , as much information as possible about the separating hyperplane H t and then use the contact evolution (in the form (5.6) or (5.8)) to determine the actual extremals of the system.

A sub-Riemannian geodesic problem
To be more precise we are considering a control system constituted by choosing in the tangent space TQ a smooth constant rank distribution D ⊂ TQ. Clearly (locally and non-canonically), by taking u 1 , u 2 , . . . , u d ) and D = f 1 , . . . , f d , we may present D as the image of a map f : Q × U → TQ where U = R d , with d := rank D, is an Euclidean space, i.e., a control system of type (CS). We shall refer to it as to the SR control system. In agreement with our notation from the previous Sect. 5, we will write also f u (q) instead of f (q, u) ∈ D q .
The SR geodesic problem is an optimal control problem of the form (OCP) constituted by considering a cost function L(q, u) : is a given symmetric positively defined bi-linear form (SR metric) on D.
In the considered situation, after passing to the extended configuration space Q = Q × R (q, q 0 ) = q, the extended control function f : Definition 10 By a SR extremal we shall understand a trajectory (q(t), u(t)) of the SR control system satisfying the necessary conditions for optimality given by the PMP (in the form provided by Theorem 2 or, equivalently, Theorems 3 or 4).

The geometry of cones and separating hyperplanes
Observe that the image f (q, U ) Fig. 4). The following fact is a simple consequence of (5.4).

Lemma 10
Let (q(t), u(t)) be a trajectory of the SR control system and let K t be the associated convex cone defined by formula (5.4). Then K t contains the tangent space of the paraboloid f (q(t), U ) at f u(t) (q(t)), i.e., Proof It follows from (5.4) (after taking k = 1, t 1 = t, and thus F tt 1 = id Q ) that K t contains every secant ray R + · { f v (q(t)) − f u(t) (q(t))} of the paraboloid

Fig. 5
Since the cone K t contains all secant rays R + · f v (q) − f u(t) (q) and is closed, it must contain the tangent space q(t)).
Using these secant rays we may approximate every tangent ray of the paraboloid f (q(t), U ) passing through f u(t) (q(t)) with an arbitrary accuracy. Since K t is closed, it has to contain this tangent ray and, consequently, the whole tangent space of f (q(t), U ) at f u(t) (q(t)) (see Fig. 5). The fact that this tangent space is described by equality (6.1) is an easy exercise.

Remark 9
In general, for an arbitrary control system and an arbitrary cost function, the cone K t contains all secant rays of the image f (q(t), U ) passing through f u(t) (q(t)). Thus, after passing to the limit, the whole tangent cone to f (q(t), U ) at f u(t) (q(t)) is contained in K t . If f (q(t), U ) is a submanifold, as it is the case in the SR geodesic problem, this tangent cone is simply the tangent space at f u(t) (q(t)).
Here is an easy corollary from the above lemma and our previous considerations.

Lemma 11 Let (q(t), u(t)) be an SR extremal and let H t ⊂ T q(t) Q be a curve of separating hyperplanes described in Theorem 4. Then for each t ∈ [0, T ], the hyperplane H t contains a rank D-dimensional linear subspace
If, additionally, (q(t), u(t)) is an abnormal SR extremal, then for each t ∈ [0, T ] there exists a hyperplane H t ⊂ T q(t) Q containing D q(t) , and such that the curve t → H t along q(t) is subject to the evolution equation

3)
Proof To justify the first part of the assertion, observe that if, in a linear space V , a hyperplane H ⊂ V supports a cone K ⊂ V which contains a line l ⊂ K (and all these sets contain the zero vector), then necessarily l ⊂ H (each line containing 0 either intersects the hyperplane or is tangent to it). Since, by Lemma 10, K t contains the subspace {Y + g( f u(t) , Y )∂ q 0 | Y ∈ D q(t) }, we conclude that this subspace must lie in H t . Assume now that the considered extremal is abnormal. In this case, as we already observed in Remark 7, H t contains, in addition to the above-mentioned linear subspace, also the line R q(t) and thus we conclude that Since {0 q }⊕R q 0 is the kernel of the natural projection Tπ 1 : T Q → TQ, we conclude that, for every t ∈ [0, T ], the image H t of H t under this projection is a hyperplane in T q(t) Q which contains D q(t) . Obviously, since f u(t) projects to f u(t) , Eq. (5.6) implies (6.2). By Theorem 1, Eq. (6.3) is the infinitesimal form of the latter.
It turns out that in some cases the above basic information, suffices to find SR extremals. Let us study the following two examples. In this case rank D = dim Q so a Riemannian extremal cannot be abnormal from purely dimensional reasons: by Lemma 11 in such a case a (dim Q − 1)-dimensional hyperplane H t ⊂ T q(t) Q would contain a bigger (dim Q)-dimensional space D q(t) = T q(t) Q, which is impossible. Thus, every Riemannian extremal must be normal and, by Lemma 11, necessarily Now any H t -valued vector field along q(t) takes the form

By (5.8) this Lie bracket should be H t -valued, and since
Y ])∂ q 0 belongs to H t , we conclude that the considered bracket belongs to H t if and only if for any Y ∈ T q(t) Q we have In this way, we have expressed the geodesic equation for the metric g in terms of the chosen metric-compatible connection ∇ with torsion T ∇ . In case that ∇ = ∇ LC is the Levi-Civita connection, the torsion vanishes and we recover the standard geodesic equation Example 2 Consider an abnormal SR extremal (q(t), u(t)) in a particular case of the SR geodesic problem where D ⊂ TQ is a co-rank one distribution. By Lemma 11 in such a situation necessarily H t = D q(t) , since the latter space is already of codimension one in T q(t) Q. Now (6.3) gives us for almost every t ∈ [0, T ], i.e., in the considered case any abnormal extremal has to be a characteristic curve of D. The converse statement is also true. Indeed, the reader may check that in this case H t := D q(t) ⊕ R q 0 (t) is the curve of separating hyperplanes containing R q(t) and satisfying the assertion of Theorem 4 (see also the proof of Theorem 5).
In the following two subsections we shall discuss normal and abnormal SR extremals in full generality.

Abnormal SR extremals
Our previous considerations allow us to give the following characterization of SR abnormal extremals.
Theorem 5 For the SR geodesic problem introduced above the following conditions are equivalent: (a) The pair (q(t), u(t)) is an abnormal SR extremal.
The smallest distribution F tτ -invariant along q(t) and containing D q(t) , i.e., is of rank smaller than dim Q. Here F tτ denotes the TD flow (in Q) of the Caratheodory TDVF f u(t) .

Moreover, condition (b) depends only on f u(t) and D along q(t).
Note that Theorem 5 reduces the problem of finding abnormal SR extremals to the study of the minimal distribution F tτ -invariant along q(t) and containing D q(t) . Often, if q(t) is sufficiently regular, this problem can be solved by the methods introduced in Lemma 3, which are more practical from the computational viewpoint.

Corollary 1 Let X be a C ∞ -smooth D-valued vector field and let q(t) with t ∈ [0, T ] be an integral curve of X . Then q(t) is an SR abnormal extremal in the following two (non-exhaustive) situations:
-The distribution spanned by the iterated Lie brackets of X with all possible smooth D-valued vector fields, i.e., is of constant rank r along q(t) and r < dim Q. -There exists a smooth distribution B ⊃ D on Q of constant co-rank at least one, such that The above fact follows directly from Theorem 5, Lemma 3 and Theorem 1. In each of the two cases along q(t) we have a constant rank smooth (and thus charming) distribution which contains D, is X -invariant along q(t) (and thus by Theorem 1 also F tτ -invariant along q(t)) and of co-rank at least one. Clearly such a distribution must contain F • (D) q(t) , which in consequence also is of co-rank at least one.
Remark 10 In Sect. 7.3 in [22] Zhitomirskii considered the following 2-distribution on R 5 where (x, y 1 , y 2 , y 3 , y 4 ) are coordinates on R 5 and smooth functions h 1 and h 2 satisfy the conditions Zhitomirskii proved that the curve (−ε, ε) t → (t, 0, 0, 0, 0) ∈ R 5 (which is obviously an integral curve of X ) is not an abnormal SR extremal, yet, as he claims, the distribution ad ∞ X (D) regarded in the above corollary is of constant rank r = 4 < 5 along this curve. A detailed study of this example reveals, however, that along the investigated curve, r = 4 apart from the point (0, 0, 0, 0, 0), where the rank drops down to 3. Thus, the discussed example does not contradict Corollary 1, as the regularity condition is not matched. In fact, the considered curve consists of two pieces of abnormal SR extremals (for t > 0 and t < 0) which do not concatenate to a single SR abnormal extremal, even though the concatenation is C ∞ -smooth. This example shows that the condition r < dim Q in Corollary 1 is not sufficient (although it is necessary in the smooth case).

Proof of Theorem 5 If (q(t), u(t))
is an abnormal extremal then, by the results of Lemma 11, H t , the TQ-projection of the curve of supporting hyperplanes H t ⊂ T q(t) Q, is a curve of hyperplanes along q(t) (i.e., a distribution of co-rank one along q(t)), it contains D q(t) and is F tτ -invariant along q(t). In particular, it must contain the smallest distribution F tτ -invariant along q(t) and containing D (cf. Proposition 2).
Conversely, assume that rank F • (D) q(t) < dim Q. Now by adding (if necessary) to F • (D) q(t) several vector fields of the form F tt 0 (X ) where X ∈ T q(0) Q, we can extend F • (D) q(t) to H t , a co-rank one distribution F tτ -invariant along q(t). Define now the curve of hyperplanes H t := H t ⊕ R q 0 (t) ⊂ T q(t) Q. We claim that H t is a curve of supporting hyperplanes described in the assertion of Theorem 3. Indeed, the F tτ -invariance of H t should be clear, as on the product Q = Q × R the TD flow F tτ takes the form F tτ (q, q 0 ) = (F tτ (q), B tτ (q 0 )), for some TD flow B tτ on R. Clearly, since H t is F tτ -invariant along q(t), the tangent map of F tτ preserves H t = H t ⊕R q 0 (t) . To prove that H t indeed separates the cone K t from the direction of the decreasing cost R − q(t) observe that any vector of the form f v (q(t)) − f u(t) (q(t)), where f v ∈ D q(t) , lies in D q(t) ⊕ R q 0 (t) ⊂ H t . Moreover, any vector of the form Thus, the whole cone K t is contained in H t (cf. formula (5.4)). Since also R − q(t) ⊂ R q(t) ⊂ H t , we conclude that indeed H t separates K t from R − q(t) (in a trivial way). Finally, to justify the last statement of the assertion we can use Theorem 1 to express the F tτ -invariance of B q(t) := F • (D) q(t) along q(t) as the f u(t) -invariance of the latter distribution, and then use Lemma 4 (for B q(t) ⊃ D q(t) f u(t) (q(t))) to prove that this invariance depends on f u(t) and F • (D) q(t) along q(t) only. Now it is enough to check that F • (D) q(t) itself does not depend on a particular choice of the extension of f u(t) (q(t)) to a neighborhood of q(t). Assume thus that f u (t) is another extension of f u(t) (q(t)), that F tτ is the related TD flow, and that F • (D) q(t) is the minimal distribution F tτ -invariant along q(t) and containing D q(t) . Now repeating the reasoning from the proof of Lemma 4 we would get Yet, for intertwined f u(t) and f u (t) we would get the opposite inclusion in an analogous manner. Thus, We claim that the integral curves of the line bundle φY + ψ Z ⊂ D are SR abnormal extremals (notice that φY + ψ Z ∈ D is a characteristic vector field of D + [D, D]). To prove this we shall use the results of Corollary 1. Indeed, it is easy to check, using (6.4), that for X = φY + ψ Z the smallest distribution ad X -invariant and containing D is simply the 3-distribution Y, Z , [Y, Z ] . This agrees with the results of Prop. 11 in [16] and Sect. 9 of [22].
Example 4 (Zelenko) The following example by Igor Zelenko [21] became known to us thanks to the lecture of Boris Doubrov. The interested reader may consult also [3,8].
Consider a 5-dimensional manifold M with a 2-dimensional distribution B ⊂ TM of type (2,3,5). That is, locally B is spanned by a pair of vector fields X 1 and X 2 such that form a local basis of sections of TM. Consider now the bundle Q := P(B) ⊂ P(TM) → M of lines in B, being a 6-dimensional manifold and a P 1 R-bundle over M. Introduce an affine chart [1 : t] corresponding to the line R · {X 1 + t X 2 } on fibers of Q → M and define a 2-dimensional distribution D := ∂ t , X 1 + t X 2 on Q. Our goal is to find abnormal SR extremals for this distribution. We will use Corollary 1 for this purpose.
First let us show that the integral curves of ∂ t are abnormal extremals. Indeed, it is easy to see that [∂ t , X 1 + t X 2 ] = X 2 and that [∂ t , X 2 ] = 0, i.e., the minimal distribution ∂ t -invariant and containing D is precisely ∂ t , X 1 , X 2 . This distribution is of constant rank smaller than 6 = dim Q, so by Corollary 1, indeed, the integral curves of ∂ t are abnormal extremals.
It is more challenging to find a second family of abnormal extremals of D. Let us look for such a family being the integral curves of the field H = X 1 +t X 2 + F∂ t , where F is some, a priori unknown, function on Q. To calculate the minimal distribution H -invariant and containing D it is enough to consider iterated Lie brackets ad k H (∂ t ). Skipping some simple calculations one can show that the vector fields In such a case, D is a constant rank distribution containing D and closed under ad H (·) (i.e., H -invariant). Since rank D = 5 < dim Q, by Corollary 1 the integral curves of H (for F as above) are abnormal SR extremals related with D.
Example 5 (Strongly bracket generating distributions) Recall that a distribution D ⊂ TQ is called strongly bracket generating (SBG) if for any p ∈ Q and any X ∈ Γ (D) non-vanishing at p we have In the light of Corollary 1 it is clear that an SR geodesic problem related with an SBG distribution does not admit any abnormal SR extremal. In fact, the same conclusion holds for a weaker version of the SBG condition, i.e., it is enough to assume that [X, D]] p + · · · = T p Q for any X ∈ Γ (D) non-vanishing at p. Example 6 (Submanifold) Assume that S ⊂ Q is a submanifold of co-dimension at least one and such that D S ⊂ TS. Then any ACB curve t → q(t) tangent to D and contained in S is an abnormal extremal. Indeed, in this case T q(t) S is obviously a charming distribution F tτ -invariant along q(t) which contains D q(t) and is of co-rank at least one in T q(t) Q. Thus, and, consequently, F • (D) q(t) is of rank smaller than dim Q. By Theorem 5, q(t) is an abnormal extremal.
Example 7 (Zhitomirskii) Let D be a 2-distribution on a manifold Q such that D 2 := D + [D, D] is of rank 3. In [22] Zhitomirskii introduced the following definition.
A distribution Z ⊂ TQ of co-dimension 2 is called nice with respect to D if In this case the intersection L := D ∩Z is a line distribution. We shall show that the integral curves of L are abnormal SR extremals. Indeed, observe that D 2 = D 2 ∩Z +D and thus is a smooth co-rank-one distribution in Q. Clearly D ⊂ H and, what is more, given any section X ∈ Γ (L) we have [X, H] ⊂ H. Indeed, take any H-valued vector field Y . Since H = Z + D we can decompose it (in a non-unique way) as Now it should be clear that the smallest distribution containing D and invariant with respect to the TD flow of Y is contained in H, which is of co-rank one. Thus, by Theorem 5, the integral curves of X are abnormal SR extremals.

Normal SR extremals
Observe first that the extremal vector field f u(t) is normalized by g( f u(t) , f u(t) ) ≡ 1 along every solution of the SR geodesic problem. Indeed, this follows easily from the standard argument involving the Cauchy-Schwartz inequality. From now on we shall thus assume that the extremal vector field f u(t) is normalized in a neighborhood of a considered trajectory q(t). This assumption allows for an elegant geometric characterization of normal SR extremals in terms of the distribution consisting of those elements of D which are g-orthogonal to f u(t) along q(t). Note that D ⊥ q(t) is a subdistribution of D along q(t).
Then, for the SR geodesic problem introduced above, the following are equivalent: is of class ACB with respect to t, and the smallest distribution F tτ -invariant along q(t) and containing D ⊥ q(t) , i.e., does not contain f u(t) (q(t)) for any t ∈ [0, T ]. Here F tτ denotes the TD flow (in Q) of the Caratheodory TDVF f u(t) .
The condition that the velocity curve t → f u(t) (q(t)) is of class ACB is equivalent to the fact that the velocities are preserved by the flow F tτ , i.e., Theorem 3.1 of [4] contains a formulation of the above result equivalent to ours. Again if q(t) is sufficiently regular we can use the method introduced in Lemma 3 to check condition (b) in the above theorem. The result stated below can be easily derived from Theorem 6 using similar arguments as in the proof of Corollary 1. For the case rank D = 2 it was proved as Theorem 6 in [16].

Corollary 2 Let X be a C ∞ -smooth D-valued vector field and let q(t) with t ∈ [0, T ] be an integral curve of X . Then q(t) is an SR normal extremal in the following two (non-exhaustive) situations:
-The distribution spanned by the iterated Lie brackets of X and all possible smooth D-valued vector fields g-orthogonal to X, i.e., is of constant rank r along q(t) and it does not contain X (q(t)) for any t ∈ [0, T ]. -There exists a smooth distribution B on Q, such that Proof of Theorem 6 The fact that property (6.5) is equivalent to t → f u(t) (q(t)) being of class ACB follows directly from Lemma 12. Indeed, the field f u(t) satisfies [ f u(t) , f u(t) ] q(t) = 0 a.e. Thus, it is ACB along q(t) if and only if it is respected by the flow F tτ .
Assume first that (q(t), u(t)) is a normal SR extremal. Let H t ⊂ T q(t) Q be the related curve of separating hyperplanes given by the PMP. Note that, since H t for each t is a hyperplane transversal to the line R q(t) ⊂ T q(t) Q, it must be of the form where α t : T q(t) Q → R is a linear map. Using the results of Lemma 11 we know that , it is clear that and Y ∈ T q(t) Q. Since T F tτ (H τ ) = H t , the above vector must be of the form That is, α t (TF tτ (X )) = α τ (X ). In particular, t → α t is continuous and, moreover, TF tτ (ker α τ ) = ker α t for every t, τ ∈ [0, T ].
We conclude that ker α t is a distribution along q(t) which is F tτ -invariant, contains D ⊥ q(t) and is transversal to f u(t) (q(t)). Clearly, F • (D ⊥ ) q(t) ⊂ ker α t and thus it is also transversal to f u(t) (q(t)).
To prove that t → f u(t) (q(t)) is ACB, observe first that D ⊥ q(t) = ker α t ∩ D q(t) admits locally a g-orthonormal basis of ACB sections. Indeed, ker α t is charming since it is F tτ -invariant (cf. Proposition 1). Let now {X 1 , . . . , X n−1 } be a local basis of ACB sections of ker α t along q(t). Choose a minimal subset of this basis, say {X 1 , . . . , X s }, such that X 1 , . . . , X s q(t) ⊕D ⊥ q(t) = ker α t for every t in a relatively compact neighborhood of a given point t 0 ∈ [0, T ]. Extend locally the SR metric g to a metric g on ker α t by taking g D ⊥ and by setting vectors X 1 , . . . , X s to be g-orthonormal and g-orthogonal to D ⊥ q(t) . Clearly, this new metric is ACB in the considered neighborhood of t 0 . Now we can apply Lemma 13 to the ACB basis {X 1 , . . . , X n−1 } and obtain an ACB g-orthonormal basis {X 1 , . . . , X s , Y s+1 , . . . , Y n−1 } of ker α t . Clearly, by the construction of the Gram-Schmidt algorithm, {Y s+1 , . . . , Y n−1 } is a g-, and thus also a g-orthonormal basis of D ⊥ q(t) (the relative compactness of the neighborhood is used to assure that the g-lengths of sections X i are separated from zero). Now let us choose any ACB section Y n of D q(t) which is transversal to D ⊥ q(t) . Again using Lemma 13 we modify the ACB local basis {Y s+1 , . . . , q(t)). Now α t ( Y n (q(t))) = ±α t ( f u(t) (q(t))) = ±1. And since both α t and Y n (q(t)) are continuous with respect to t the sign ± must be constant along [0, T ]. We conclude that t → f u(t) (q(t)) is ACB alike t → Y n (q(t)).
Conversely, assume that (b) holds. The crucial step is to build, along the projected trajectory q(t) ∈ Q, a splitting T q(t) Q = B q(t) ⊕ f u(t) (q(t)) , where B q(t) is a co-rank one distribution along q(t), which is F tτ -invariant along q(t) and contains D ⊥ q(t) . Such a B q(t) can be constructed by adding, if necessary, to F • (D ⊥ ) q(t) several vector fields of the form F t0 (X i ), where X i ∈ T q(0) Q together with F u(0) (q(0)) are independent. Clearly, in this way we can build B q(t) which is F tτ -invariant along q(t), of co-rank one and contains D ⊥ q(t) . What is more, since f u(t) (q(t)) is F tτ -invariant, the flow F tτ respects the splitting B q(t) ⊕ f u(t) (q(t)) . Now we can construct the curve of separating hyperplanes H t for q(t) by the formula By construction it is clear that H t is a hyperplane in T q(t) Q which contains the tangent space to the paraboloid (6.1) and does not contain the line R q(t) . From these properties we conclude that H t separates strictly the ray R − q(t) from the elements of q(t)). It is also clear (here we use the normalization of f u(t) , cf. the first part of this proof) that H t is F tτ -invariant along q(t). Thus, it also separates strictly the ray R − Consequently, using the fact that H t and R − q(t) are convex, we can use H t to separate strictly R − q(t) from any finite convex combination of the above-mentioned elements of K t . Since K t is by definition the closure of the set of such finite convex combinations, H t is indeed the separating hyperplane described by the PMP.
A remark on smoothness of normal SR geodesics As was proved above normal SR extremals are C 1 -smooth (and even more: their derivatives are ACB maps). It is worth discussing geometric reasons for this regularity in a less technical manner than in the proof of Theorem 6. Let q(t) be such an extremal and let H t be the corresponding curve of supporting hyperplanes. As we know from Lemma 11 These two facts are enough to exclude, at least in a heuristic way, the existence of singularities of corner-type and of cusp-type along q(t). Indeed, since H t is F tτinvariant it must be continuous. Note that by the continuity of H t , the limit subspaces D ⊥ q(t 0 )± ⊕ 0 · ∂ q 0 coming from both sides of a given point t 0 ∈ [0, T ] must belong to H t 0 . Now if q(t) had a corner-type singularity at t 0 , these limit subspaces would be different and thus they would span together the whole space D q(t 0 ) ⊕ 0 · ∂ q 0 (cf. Fig. 6). In particular, f u(t 0 ) (q(t 0 )) + 0 · ∂ q 0 ∈ D q(t 0 ) ⊕ 0 · ∂ q 0 would belong to H t 0 .
Yet, since f u(t 0 ) (q(t 0 )) + ∂ q 0 ∈ H t 0 , this would mean that also the difference of the latter vectors, 0 + ∂ q 0 , lies in H t 0 , which is impossible since q(t) is normal.
In a similar way one deals with a cusp-type singularity. At a cusp we would have limit vectors ± f u(t 0 ) (q(t 0 )) + ∂ q 0 in H t 0 (see Fig. 6). Now 0 + 2∂ q 0 , the sum of these two vectors, would belong to H t 0 which contradicts the normality of the extremal. Roughly speaking, the existence of singularities of corner-type or cusp-type implies ∂ q 0 ∈ H t 0 , i.e., either a trajectory is not an extremal or it is abnormal. Fig. 6 The existence of singularities of corner-or cusp-type implies abnormality or the lack of optimality Examples Example 8 (Geodesic equation revisited) Theorem 6 provides an alternative way to derive the geodesic equation in the Riemannian case (i.e., when D = TQ). Let (q(t), u(t)) be a trajectory of the SR control system (we shall assume that f u(t) is normalized). Since D = TQ, by the assertion of Theorem 5, in the Riemannian case there are no abnormal extremals. Since for every t, τ ∈ [0, T ]. By the results of Theorem 1 this is equivalent to i.e., g ([ f u(t) , Y ], f u(t) ) q(t) = 0 whenever g(Y, f u(t) ) q(t) = 0. Now for such a Y , after introducing a metric-compatible connection as in Example 1, we have Using the fact that g(Y, f u(t) ) ≡ 0 and that g( f u(t) , f u(t) ) ≡ 1 we get in agreement with the results of Example 1.
Example 9 (Heisenberg system) Consider an SR system on R 3 (x, y, z) constituted by a 2-distribution D (x,y,z) = Y := ∂ x − y∂ z , Z := ∂ y + x∂ z and an SR metric such that the fields Y and Z form an orthonormal basis. Such a system is usually called the Heisenberg system. It is easy to check that the system in question is strongly bracket generating (cf. Example 5) and as such does not admit any abnormal SR extremal. Our goal will thus be to determine the normal SR extremals using the results of Theorem 6. Take now any normalized D-valued vector field X = X t := φ(t)Y + ψ(t)Z , where φ 2 + ψ 2 = 1. We have D ⊥ = X = X t := ψ(t)Y − φ(t)Z and, by the results of Theorem 6, the integral curve q(t) of X is an SR normal extremal if and only if F • (D ⊥ ) q(t) does not contain X at any point q(t). Clearly distribution F • (D ⊥ ) q(t) , being ad X -invariant, contains the fields X , [X, X ], [X, [X, X ]], etc. Skipping some simple calculations one can show that Let us now present vector field [X, X ] as (note that {X, X } is a basis of sections of D) Then Thus, by (6.6), Consequently, which, after substituting φY + ψ Z by X , leads to i.e., the quotient ψ/φ satisfies the equation For α = 0 we get x = const (i.e., φ and ψ are constant along q(t)), and for α = 0 we get x = arctan(αt +γ ) (i.e., φ = cos(αt +γ ) and ψ = sin(αt +γ )). This corresponds to the two well-known families of normal SR extremals of the Heisenberg system (see Sect. 2 of [16]), whose projections to the (x, y)-plane are straight lines and circles, respectively.
Here the limit is taken over open neighborhoods V t and μ(·) denotes the Lebesgue measure on V . By Lebesgue theorem, the set of regular points of a bounded and measurable map f : V → R n is of full measure in V .
A map x : R ⊃ [t 0 , t 1 ] → R n is called absolutely continuous (AC) if it can be presented in the form of an integral for some integrable map v(·). Clearly, an AC map is differentiable at all regular points t of v (and thus, by Lebesgue theorem, a.e.) and the derivative of x(t) at such a point is simply v(t). We will be particularly interested in AC maps x(t) such that the derivative v(t) is locally bounded. In such a case we shall speak about AC maps with bounded derivative (ACB).
Measurable ODE's While speaking about ODE's in the measurable setting we will need to take care of some technical properties of certain functions. In order to simplify the discussion let us introduce the following The notion of a Caratheodory map can be naturally extended to the setting of smooth manifolds, namely we shall call a map F : M × R → N Caratheodory if it is Caratheodory in a (and thus in any) local smooth coordinate chart on M and N . Indeed, it is easy to see that this property does not depend on the particular choice of a chart (cf. the notion of a Caratheodory section in [11]).
Consider now a map G : R n × R → R n and the associated non-autonomous ODE in R nẋ By a (Caratheodory) solution of (A.3) on [t 0 , t 1 ] with the initial condition x 0 at t 0 we shall understand an AC map [t 0 , t 1 ] t → x(t) ∈ R n which satisfies (A.3) a.e. (recall that an AC map is differentiable a.e.), such that x(t 0 ) = x 0 . Note that speaking about Caratheodory solutions makes sense also if the map G is defined only a.e.
The following fact is a straightforward generalization, to the measurable context, of the standard result about the existence and uniqueness of the solutions of ODE's.

Theorem 7
Assume that the map (x, t) → G(x, t) is Caratheodory. Then, for each choice of (t 0 , x 0 ) ∈ R × R n there exists, in a neighborhood of t 0 , a unique Moreover, x(t; t 0 , x 0 ) is differentiable with respect to x 0 and the derivative ∂ x ∂ x 0 (t; t 0 , x 0 ) is continuous with respect to x 0 and ACB with respect to t. In fact, the derivative ∂ x ∂ x 0 (t; t 0 , x 0 ) is the unique (Caratheodory) solution of the following linear time-dependent ODE, called the variational equation, for a curve of linear maps V (t, x 0 ) : T x 0 R n → T x(t;t 0 ,x 0 ) R n with the initial condition The proof can be found in [6] (Theorem 3.3.2). Also Sect. 3 of [10] may be useful. Note that Eq. (A.4) can be obtained by differentiating the equationẋ(t; t 0 , x 0 ) = G(x(t; t 0 , x 0 ), t) with respect to x 0 and substituting V (t, x 0 ) for ∂ x ∂ x 0 (t; t 0 , x 0 ).
Proof of Theorem 1 In this paragraph we will provide a rigorous proof of Theorem 1. We shall begin with the following lemma which characterizes the tangent map T A tτ of the TD flow of a TDVF X t in terms of the Lie bracket [X t , ·]. Informally speaking, transporting a given vector Z 0 via the map T A tt 0 along an integral curve of X t turns out to be the same as solving the equation [X t , ·] = 0.

Lemma 12
Let X t be a Caratheodory TDVF on the manifold M, x(t) = x(t; t 0 , x 0 ) (with t ∈ [t 0 , t 1 ]) its integral curve, and A tt 0 its TD flow. Let Z 0 ∈ T x 0 M be a tangent vector at x 0 and denote by Z (x(t)) a vector field along x(t) obtained from Z 0 by the action of the TD flow A tt 0 , i.e., Z (x(t)) := T A tt 0 (Z 0 ). Then the assignment t → Z (x(t)) is ACB and, moreover, Conversely, if Z is a vector field along x(t) such that the assignment t → Z (x(t)) is ACB and that Eq. (A.5) holds, then Z (x(t)) = T A tt 0 (Z 0 ), where Z 0 = Z (x(t 0 )).
Proof Consider first the vector field Z (x(t)) := T A tt 0 (Z 0 ) along x(t). The fact that t → Z (x(t)) is ACB follows directly from the second part of the assertion of Theorem 7.
Proof of Theorem 1 Assume first that condition (a) of Theorem 1 holds. Choose a basis {Z 10 , . . . , Z k0 } of B x(t 0 ) , where k is the rank of B, and for i = 1 . . . , k denote Z i (x(t)) := T A tt 0 (Z i0 ). By the results of Lemma 12, the fields Z i are ACB along x(t) and satisfy [X t , Z i ] x(t) ≡ 0 a.e. along x(t). Thanks to condition (a) and the fact that A tt 0 is a local diffeomorphism, Z i 's span B.
Let now Z ∈ Γ AC B (B) be any ACB section of B. We want to present it as a linear combination of fields Z i with ACB coefficients, i.e., Z = i φ i Z i , where φ i are ACB functions along x(t). To prove that such a presentation is possible first take vectors W j0 ∈ T x(t 0 ) M with j = 0, . . . , s such that Z i0 's together with W j0 's form a basis of T x(t 0 ) M. Clearly, the fields W j (x(t)) := T A tt 0 (W j0 ) together with Z i 's span TM along x(t). Since by Lemma 12 these fields are ACB, given any local basis of smooth vector fields U := {U 1 , . . . , U k+s } on M, the transition matrix T U →ZW from the basis U to the basis ZW := {Z 1 , . . . , Z k , W 1 , . . . , W s } is a matrix of ACB functions. As T U →ZW is non-degenerate, the inverse matrix T ZW→U is also a matrix of ACB functions (here we use the fact that if φ is an ACB function separated from 0, then so is 1 φ ). Thus, any vector field with ACB coefficients in basis U (in particular Z ) will have ACB coefficient in basis ZW. As the field Z is B-valued, all W j 's coefficients of Z vanish, i.e., Z = i φ i Z i , where φ i are ACB functions along x(t) as intended. Now by the Leibniz rule 6 we get Thus, (a) implies (b).
Assume now that condition (b) of Theorem 1 holds. Let { Z 1 , . . . , Z k } be any basis of ACB sections of B. The idea is to modify this basis to another basis of ACB sections {Z 1 , . . . , Z k }, such that for every i = 1, . . . , k we have [Y t , Z i ] x(t) ≡ 0 a.e. along x(t). In the light of Lemma 12 this would imply that the new basis is respected by the flow A tt 0 and, consequently, that (a) holds.
Due to (b), [Y t , Z i ] x(t) is a B-valued locally bounded measurable vector field for each i = 1, . . . , k and thus there exists a k × k matrix of locally bounded measurable functions 7 φ j i along x(t) such that . 6 We leave the proof of the fact that the Lie bracket (2.3) satisfies the Leibniz rule as an exercise. 7 The existence of measurable functions φ j i can be justified in a similar manner to the existence of ACB functions φ i above. Now the simple idea is to look for the desired basis {Z 1 , . . . , Z k } in the form Z i = j ψ j i Z j , where ψ j i is an invertible k ×k matrix of function ACB along x(t). Clearly for such Z i 's we have .
As we see [X t , i.e., the matrix ψ j i should be a solution of a linear ODE with locally bounded measurable coefficients. Due to the results of Theorem 7, for a given initial condition, say, ψ j i (x 0 ) = δ j i , this equation has a unique local ACB solution. As a consequence, we prove the local existence of the desired basis {Z 1 , . . . , Z k }, which implies (a).

Proof of Lemma 4
We shall end our considerations by providing the following Proof (of Lemma 4) The idea of the proof is very simple. Consider another D-valued Caratheodory TDVF X t such that X t = X t along x(t). We shall show that B is X tinvariant if and only if it is X t -invariant along x(t). The justification of this statement is just a matter of a calculation. Observe that since X t and X t are both Caratheodory and D-valued, then so is their difference X t − X t . Given any local basis of smooth vector fields {W 1 , . . . , W s } of D we may locally represent X t − X t as where φ i (t, x) are Caratheodory functions (in the sense of Definition 11). Since X t = X t along x(t) and W i 's form a basis of D, we have φ i (t, x(t)) = 0. (A.7) Now for any section Z ∈ Γ AC B (B), using the same notation as in formula (2.3), we have = ∂ ∂s s=0 X t (z(t, s)) − ∂ ∂s s=0 X t (z(t, s)) = ∂ ∂s s=0 X t (z(t, s)) − X t (z(t, s)) = ∂ ∂s s=0 i φ i (t, z(t, s))W i (z(t, s)) = i ∂ ∂s s=0 φ i (t, z(t, s))W i (x(t))) + φ i (t, x(t)) ∂ ∂s s=0 W i (z(t, s)) (7.7) = i ∂ ∂s s=0 φ i (t, z(t, s))W i (x(t)).
Clearly the above expression is D x(t) ⊂ B x(t) -valued. Thus, along x(t) and hence, since Z was an arbitrary ACB section of B, This ends the proof.
A technical result about charming distributions The following result will be needed in the course of Sect. 6.2 to prove that normal SR extremals are of class C 1 . It states that the Gram-Schmidt orthogonalization algorithm works well on charming distributions. Proof Recall that the Gram-Schmidt algorithm maps a set {X 1 , . . . , X s } into a gorthonormal set {X 1 , . . . , X s } constructed in the following way X 1 → X 1 := X 1 , . . .
proj X i X s , where proj X Y := g(X,Y ) g(X,X ) X denotes the g-orthogonal projection of Y on the space spanned by X . Now it is enough to use the following elementary facts concerning ACB functions: -A sum and a difference of two ACB functions is an ACB function.
-A product of two ACB functions is an ACB function.
Clearly, in every step of the Gram-Schmidt algorithm we apply one or more of these elementary operations to ACB sections (we use the fact that g is ACB and that for every ACB non-zero vector X the values of g(X, X ) are separated from zero on [t 0 , t 1 ].) Thus, as a result we also obtain ACB sections X i 's.