I am working on sub-Riemannian geometry and try to understand what are the tools to find the equations of a sub-Riemannian problem. Here is an example:

Let us consider the system defined by a lagrangian:

\begin{equation}
L=\frac{1}{2}m(\overset{\cdot}{x}^2+\overset{\cdot}{y}^2)+\frac{1}{2}I\overset{\cdot}{\theta}^2
\end{equation}

And take the following non-holonomic constraint 

\begin{equation}
\mathscr{D}=span\left\{\cos\theta\frac{\partial}{\partial x}+\sin\theta\frac{\partial}{\partial y},\frac{\partial}{\partial\theta}\right\}
\end{equation}
So we have,
\begin{equation}
\mathscr{D}^\perp=span\{\omega\}
\end{equation}
with $\omega=-\sin\theta dx+\cos\theta dy$

If I understand clearly the sub-Riemannian problem is to minimize 
\begin{equation}
\int_0^TLdt
\end{equation}
under the constaint $\dot{z}(t)\in\mathscr{D}(z(t))$. Now, how do I find these equations? So far I am not asking how to solve these equations but just to understand how to find them in the first place.

Thanks in advance for your help.

Edit:

I quite understand your answer, thank you for the clear explanation. I am though not really familiar with canonical coordinates. I agree with the expression of the Hamiltonian.

Now, I am right if I say that the Hamilton's equations are given by:

\begin{equation}
	\begin{cases}
		\overset{\cdot}{x}(t)&=\frac{\partial H}{\partial p_x}=p_x\cos^2\theta +p_y\cos\theta\cdot\sin\theta \\
		\overset{\cdot}{y}(t)&=\frac{\partial H}{\partial p_y}=p_y\sin^2\theta +p_x\cos\theta\cdot\sin\theta \\
		\overset{\cdot}{\theta}(t)&=\frac{\partial H}{\partial p_\theta}=p_\theta\\
		\overset{\cdot}{p_x}(t)&=-\frac{\partial H}{\partial x}= 0 \\
		\overset{\cdot}{p_y}(t)&=-\frac{\partial H}{\partial y}= 0 \\
		\overset{\cdot}{p_\theta}(t)&=-\frac{\partial H}{\partial \theta}= (p_x^2-p_y^2)\sin(\theta)cos(\theta)+p_xp_y(\sin^2(\theta)-\cos^2(\theta)) \\
	\end{cases}
\end{equation}

I'm not sure to understand how to interpret them.
Are these the equations of the sub-Riemannian problem?