On Integrability of the Sub-Riemannian Geodesic Flow for Goursat Distribution

Consider the following optimal control problem:
$$ \dot{q}=u_1 f_1(q)+u_2 f_2(q), \qquad q \in \mathbb{R}^n, \quad u \in \mathbb{R}^2, $$
where $q = (x_1,x_2,\dots,x_n)^T, \ f_1=(1,0,-x_2,-x_3,\dots,-x_{n-1})^T,f_2=(0,1,0,0,\ldots,0)^T,$ boundary conditions:
$$ q(0)=q_0, \ q(t_1)=q_1, $$
quality functional:
$$ l=\int_{0}^{t_1} \sqrt{u^2_1+u^2_2} \ dt \rightarrow min, $$
where the point $q \in \mathbb{R}^n$ determines the state of the system, $u=(u_1,u_2)$ is a control, $t_1$ being fixed.
$ $
Notice that $f_1, f_2$ can be chosen just as they are up to any diffeomorphism. This is how the commutators of $f_1$ and $f_2$ look like:
\begin{align*} f_3&= \frac{\partial{f_2}}{\partial{q}}f_1 - \frac{\partial{f_1}}{\partial{q}}f_2 =[f_1,f_2] =(0,0,1,0,\dots,0)^T, \\ f_4&= \frac{\partial{f_3}}{\partial{q}}f_1 - \frac{\partial{f_1}}{\partial{q}}f_3 =[f_1,f_3]=(0,0,0,1,\dots,0)^T, &\ldots \\ f_n&= \frac{\partial{f_{n-1}}}{\partial{q}}f_1 - \frac{\partial{f_1}}{\partial{q}}f_{n-1} =[f_1,f_{n-1}]=(0,0,\dots,0,1)^T. \end{align*}

Nilpotent Lie algebra is generated by $f_1, f_2$:
$$ Lie(f_1, f_2) = span (f_1, f_2, \dots , f_n), $$
multiplication table being of the form:
$$ [f_1,f_2] = f_3, \ [f_1,f_3] = f_4, \ \dots \ , \ [f_1,f_{n-1}] = f_n. $$
All the others are equal to zero. These relations define the so-called Goursat distribution ([1], [2]).
$ $
Using the Pontryagin maximum principle, one can construct the Hamiltonian
$$ H(q,p,u) = \frac{u^2_1 +u^2_2}{2}=\frac{\dot{x_1}^2+\dot{x_2}^2}{2}=\frac{\langle p,f_1 \rangle ^2 + \langle p,f_2 \rangle ^2}{2}, $$
thus obtaining the following system:
$$ \tag{1} \begin{cases} \dot{x_1}=p_1-x_2 p_3 - \dots -x_{n-1} p_n, \\ \dot{x_2}=p_2, \\ \dot{x_3}=-x_2 \dot{x_1}, \\ \ldots \\ \dot{x_n}=-x_{n-1} \dot{x_1}, \\ \dot{p_1}=0, \\ \dot{p_2}=p_3 \dot{x_1}, \\ \ldots \\ \dot{p_{n-1}}=p_n \dot{x_1}, \\ \dot{p_n}=0. \\ \end{cases} $$
Let us introduce the new coordinates by the following way:
\begin{equation*} \begin{cases} P_1=p_1-x_2 p_3-x_3 p_4-\dots-x_{n-1} p_n, \\ P_n=p_n, \\ P_{n-1}=p_{n-1}-P_n x_1, \\ P_{n-2}=p_{n-2}-P_{n-1} x_1-P_n \frac{x_1^2}{2!}, \\ \dots \\ P_3=p_3-P_4 x_1-P_5 \frac{x_1^2}{2!}-\dots-P_n \frac{x_1^{n-3}}{(n-3)!}, \\ P_2=p_2. \\ \end{cases} \end{equation*}

The following theorem holds.
$ $
Theorem. (1) is the completely integrable system (in the Liouville sense). The whole set of the first integrals is as follows:
\begin{equation*} \begin{cases} F_n=P_n, \\ F_{n-1}=P_{n-1}, \\ \dots \\ F_3=P_3, \\ F_2=P_2-P_3 x_1-\dots-P_n \frac{x_1^{n-2}}{(n-2)!}, \\ F_1=H=\frac{1}{2} (P_1^2+P_2^2), \\ \end{cases} \end{equation*}


Thus one can consider the following “moment map”:
$$\Phi: (x, P) \rightarrow \begin{pmatrix} F_n \\ \ldots \\ F_1 \\ \end{pmatrix}. $$

The primary aim here is to study critical points of this mapping and its properties. That's what we are keep working on.