How to generalize the Euclidean "unicycle" model?

There is a common system of ODEs known as the unicycle model / Dubin's model which describes the kinematics of an ant-like "unicycle" that can drive forward with some velocity $v(t) \in \mathbb{R}$ and turn in-place with some angular velocity $\omega(t) \in \mathbb{R}$ about the ground normal. The position of the unicycle / ant is described by $\big{(}x(t), y(t)\big{)} \in \mathbb{R}^2$ and its orientation is described by its "heading" angle $\ \theta(t) \in \mathcal{S}^1$.

\begin{align} \dot{x} &= v\cos(\theta)\\[2pt] \dot{y} &= v\sin(\theta)\\[2pt] \dot{\theta} &= \omega \end{align}

These equations rely on the ant traversing a flat 2D plane. I have been thinking a lot about how to generalize this simple model to an arbitrarily curved (but smooth) 2D surface.

We can imagine the infinity of geodesics that intersect $\big{(}x(t_0), y(t_0)\big{)}$. Then $\theta(t_0)$ essentially picks one of those geodesics. While $v > 0$ and $\omega = 0$, the ant moves what it thinks is "straight forward." It will follow the same geodesic until it chooses to turn ($\omega \neq 0$) onto a different geodesic.

Suppose our smooth surface is defined by an injective function that embeds it in 3D Euclidean space: $$ r : \mathbb{R}^2 \to \mathbb{R}^3\ ,\ \ \ r(x,y) = \begin{bmatrix} \bar{x}(x,y) \\ \bar{y}(x,y) \\ \bar{z}(x,y) \end{bmatrix} $$

We have its Jacobian and any necessary higher partial derivatives as well: $$ J(x,y) := \begin{bmatrix}\frac{\partial r}{\partial x}(x,y)\,\ \frac{\partial r}{\partial y}(x,y)\end{bmatrix} \in \mathbb{R}^{3 \times 2} $$

The velocity vector is, $$ \dot{r} = \dot{x}\frac{\partial r}{\partial x} + \dot{y}\frac{\partial r}{\partial y} = J \begin{bmatrix} \dot{x} \\ \dot{y} \end{bmatrix} $$

The scalar velocity $v$ in the flat model would correspond to the magnitude of the velocity vector in the general model, i.e. $v = ||\dot{r}||$. Letting $\hat{\tau}$ be the direction of $\dot{r}$ we can write, $$ \dot{r} = v\hat{\tau} = J \begin{bmatrix} \dot{x} \\ \dot{y} \end{bmatrix} $$

Non-degeneracy of the surface-coordinates implies that $\frac{\partial r}{\partial x}$ is linearly independent of $\frac{\partial r}{\partial y}$. Thus, $J$ is full-column-rank. Left-multiplying the above by $J^\intercal$ and inverting yields, $$ \begin{bmatrix} \dot{x} \\ \dot{y}\end{bmatrix} = v(J^\intercal J)^{-1}J^\intercal\hat{\tau} \tag{1} $$

It's noteworthy that $J^\intercal J$ is the metric tensor in coordinates. This is almost the ODE I'd need for $\big{(}x, y\big{)}$, but $\hat{\tau}$ is undefined when evaluating the right-hand-side. I tried defining $\hat{\tau}$ in terms of some rotation $\theta$ of the surface gradients (with $\dot{\theta} = \omega$) and got cool but erroneous results:

$$ \text{This ^ uses Eq.1 with:}\ \ \ \ \hat{\tau} = \cos(\theta)\frac{\frac{\partial r}{\partial x}}{\big{|}\big{|}\frac{\partial r}{\partial x}\big{|}\big{|}} + \sin(\theta)\frac{\frac{\partial r}{\partial y}}{\big{|}\big{|}\frac{\partial r}{\partial y}\big{|}\big{|}} ,\ \ \ \ \dot{\theta}=\omega$$ The error is that for $\omega = 0$ and $v \neq 0$, holding $\hat{\tau}$ some fixed $\theta$ away from a gradient direction doesn't cause the ant to move along a geodesic. Moreover, it doesn't make sense for $\theta=0$ to yield a $\frac{\partial r}{\partial x}$-follower... $\hat{\tau}$ really needs to be parallel-transported along the surface. Parallel-transport is a statement about how tangent vectors (e.g. $\dot{r}$) change, so we need to examine acceleration.

If the ant never turned or varied speed, then its acceleration would be strictly normal to the surface (acceleration just due to curvature): \begin{gather} \langle \frac{\partial r}{\partial x}, \ddot{r} \rangle \ \ \overset{\dot{v}=\omega=0}{=}\ \ 0\\[3pt] \langle \frac{\partial r}{\partial y}, \ddot{r} \rangle \ \ \overset{\dot{v}=\omega=0}{=}\ \ 0 \end{gather}

The left-hand-side can be expressed more compactly as $J^\intercal \ddot{r}$. Expanding $\ddot{r}$, \begin{align} \ddot{r} &= J \begin{bmatrix} \ddot{x} \\ \ddot{y} \end{bmatrix} + \dot{J} \begin{bmatrix} \dot{x} \\ \dot{y} \end{bmatrix}\\[4pt] J^\intercal \ddot{r} &= J^\intercal J \begin{bmatrix} \ddot{x} \\ \ddot{y} \end{bmatrix} + J^\intercal \dot{J} \begin{bmatrix} \dot{x} \\ \dot{y} \end{bmatrix} \ \ \overset{\dot{v}=\omega=0}{=}\ \ 0 \end{align}

and finally left-multiplying by the inverse metric tensor $(J^\intercal J)^{-1}$, $$ \begin{bmatrix} \ddot{x} \\ \ddot{y} \end{bmatrix} + (J^\intercal J)^{-1} J^\intercal \dot{J} \begin{bmatrix} \dot{x} \\ \dot{y} \end{bmatrix} \ \ \overset{\dot{v}=\omega=0}{=}\ \ 0 \tag{2} $$

This is the geodesic equation. The second term can be simplified with Christoffel symbols, but that isn't really necessary. $\dot{J} = \dot{x}\frac{\partial J}{\partial x} + \dot{y}\frac{\partial J}{\partial y}$ is just another computable function to me. I can use Equation 1 to convert some arbitrary initial $\hat{\tau}$ heading direction into an initial condition on $\big{(} \dot{x}, \dot{y} \big{)}$ and then integrate Equation 2 for the motion along the geodesic as long as $\dot{v} = \omega = 0$.

I suspect that to generalize this to the $\omega \neq 0$ and $\dot{v} \neq 0$ case, the right-hand-side will have to be some function of them rather than $0$. But what function? Edit: And how do we encode "heading" so that the ant can "stand still" when $v = 0$, and even turn "in-place" (about the normal)?

What (if any) ODE system describes the ant / unicycle moving on an arbitrary smooth 2D surface? $$ \overset{?}{f}(x,y,\theta;v,\omega,r) = 0 $$

where $f$ can involve any derivatives up to second-order. Thanks in advance! :)

Addendum

The answer given by Kajelad below seems correct! Though for posterity I would like to add an elaboration of that answer here, using notation consistent with the above and providing more motivation / steps.

First lets define the Frenet-Serret basis: $$ \hat{\tau} := \frac{\dot{r}}{||\dot{r}||},\ \ \ \ \hat{\eta} := \frac{\frac{\partial r}{\partial x} \times \frac{\partial r}{\partial y}}{\big{|}\big{|}\frac{\partial r}{\partial x} \times \frac{\partial r}{\partial y}\big{|}\big{|}},\ \ \ \ \hat{\beta} := \hat{\eta} \times \hat{\tau} $$

These are the standard tangent, normal, and binormal vectors respectively ($\hat{\tau}$ being the same as I had previously defined) and they span $\mathbb{R}^3$ orthonormally at all points on the surface. Thus, the ambient acceleration vector can be expressed as some linear combination of them. $$ \ddot{r} = a\hat{\tau} + b\hat{\beta} + c\hat{\eta} $$

From the previous discussion, we know that if $\dot{v} = \omega = 0$, then $\ddot{r} = c\hat{\eta}$ will yield geodesic motion. Therefore $a$ and $b$ must depend multiplicatively on $\dot{v}$ and $\omega$. We also know that this must reduce to the traditional model in the Euclidean case. As we'll see, the correct relations are (rather intuitively), $$ a = \dot{v},\ \ \ \ b = \omega v $$

For clarity in the following algebra I'll define these shorthands, $$ q := \begin{bmatrix} x \\ y \end{bmatrix},\ \ \ \ \tilde{J} := (J^\intercal J)^{-1}J^\intercal $$

This makes our previous Eq.1 appear nicely as $\dot{q} = v\tilde{J}\hat{\tau}$, and we also have $\ddot{r} = J\ddot{q} + \dot{J}\dot{q}$. Thus, $$ J\ddot{q} = \dot{v}\hat{\tau} + \omega v \hat{\beta} + c\hat{\eta} - \dot{J}\dot{q} $$

Left-multiplying by $\tilde{J}$ not only cancels the leftmost $J$, but also kills the $c\hat{\eta}$ term because $J^\intercal \hat{\eta}=0$. $$ \ddot{q} = \dot{v}\tilde{J}\hat{\tau} + \omega v \tilde{J}\hat{\beta} - \tilde{J}\dot{J}\dot{q} \tag{3} $$

Next, lets define our generalization of "heading" as the unique vector $h(t) \in \mathbb{R}^2$ satisfying, \begin{align} \hat{\tau} &= Jh\\[4pt] h &= \tilde{J}\hat{\tau} = \tfrac{1}{v} \dot{q} \end{align}

This vector exists because $\hat{\tau}$ lives in the tangent space (the span of $J$). Furthermore, $\hat{\beta}$ also lives in the tangent space as a $90^\circ$ rotation about the normal from $\hat{\tau}$. Calling that transform $R^{\hat{\eta}}_{90^\circ} \in \mathbb{SO}3$, we have, \begin{align} \hat{\beta} &= R^{\hat{\eta}}_{90^\circ} \hat{\tau}\\[4pt] \tilde{J} \hat{\beta} &= (J^\intercal J)^{-1} J^T R^{\hat{\eta}}_{90^\circ} J h \tag{$\hat{\tau} \to Jh$}\\[4pt] &= \sqrt{|J^\intercal J|}(J^\intercal J)^{-1}\begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix}h\\[4pt] &=: Rh \end{align}

We can derive an ODE for $h$ by using the above relations to put Eq.3 in terms of $h$. \begin{align} \ddot{q} &= \dot{v}h + \omega v Rh - v\tilde{J}\dot{J}h\ \ \tag{$\tilde{J}\hat{\tau} \to h,\ \tilde{J}\hat{\beta} \to Rh,\ \dot{q} \to v h$} \\[4pt] \dot{v}h + v\dot{h} &= \dot{v}h + \omega v Rh - v\tilde{J}\dot{J}h \tag{$\ddot{q} \to \dot{v}h + v\dot{h}$}\\[4pt] \dot{h} &= \omega Rh - \tilde{J}\dot{J}h \tag{cancellations} \end{align}

The $\dot{J}$ is still hiding a $\dot{q}$ dependence, so lets expand it in terms of the Hessian $H := \frac{dJ}{dq}$. $$ \dot{J} = \dot{x}\frac{\partial J}{\partial x} + \dot{y}\frac{\partial J}{\partial y} =: H \odot \dot{q} = v H \odot h $$

Finally, we arrive at, $$ \dot{h} = \omega Rh - v\tilde{J}(H \odot h)h $$

The rightmost-term can be simplified using Christoffel symbols, but that isn't critical as it is already something I can compute. This equation mirrors Kajelad's result: $$ \dot{T}^i = \omega R^i{}_j T^j-v\Gamma^i{}_{jk}T^jT^k $$

So in summary, the generalized model is: \begin{gather} \begin{bmatrix} \dot{x} \\ \dot{y} \end{bmatrix} = vh \tag{4}\\[4pt] \dot{h} = \omega Rh - v(J^\intercal J)^{-1}J^\intercal(H \odot h)h \tag{5}\\[4pt] R = \sqrt{|J^\intercal J|}(J^\intercal J)^{-1}\begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix} \tag{6} \end{gather} where heading $h(t) \in \mathbb{R}^2$, $J$ and $H$ are the embedded surface's Jacobian and Hessian respectively, and $v$ and $\omega$ are the specified ant/unicycle movements. This could be (more commonly) written in terms of the metric tensor ($J^\intercal J$) and Christoffel symbols (gradients of the metric tensor). The appeal of that form is that it's completely independent of any embedding of our surface $r$; only the metric tensor field is needed.

Note that the units are consistent ($H$ has units of inverse-space). In the Euclidean case, the Hessian vanishes and the ODE reduces to the classic model. Lastly, if $\omega = 0$, the $\dot{h}$ equation becomes a parallel-transport / geodesic equation. Now to test in simulation!

(Concluding remark: While very mechanical, the above derivation is a reflection of some rather fundamental and cool concepts in differential geometry, which Kajelad's answer used explicitly and may have gone beyond my background with. The key idea here, I think, is the Levi-Civita connection, or more specifically, the "covariant-derivative." We took a covariant-derivative in the above when we took a regular derivative of $\dot{r}$ and then projected it onto the tangent space by left-multiplying $\tilde{J}$. Instead, I could have used the abstract parlance and just defined $\nabla_{\dot{r}} \dot{r} = \ddot{r} - c\hat{\eta}$. The ugliness in my derivation above is in how tied it is to coordinate representations of all the geometric quantities involved. To those unfamiliar with differential geometry, however, it is far more grounded, and can even serve as a gateway to the beautiful abstraction.)

Let $R$ be a consistently chosen rotation by $90^\circ$. Such an operator always exists locally, e.g. for embedded surfaces in $\mathbb{R}^3$ one can choose a unit normal vector field $n$ and define $R(v):=n\times v$. Let $\gamma:\mathbb{R}\to S$ be the parameterized curve of interest. The "most general" second order equations of motion for such a curve can be written in terms of the vector acceleration $A$, which describes how $\gamma$ departs from geodesic motion: $$ \nabla_{\dot{\gamma}}\dot{\gamma}=A $$ Which in local coordinates can be written in terms of the Christoffel symbols $\Gamma$: $$ \ddot{\gamma}^i+\Gamma^i{}_{jk}\dot{\gamma}^j\dot{\gamma}^k=A^i $$ One can split $A$ into a parallel and perpendicular components, which can be parameterized by the scalar acceleration $a$ and the angular velocity $\omega$. $$ \nabla_{\dot{\gamma}}\dot{\gamma}=a\frac{\dot{\gamma}}{\|\dot{\gamma}\|}+\omega R(\dot{\gamma}) \\ $$ Not that $a$ is indeed the acceleration as typically defined by $a=\frac{d}{dt}\|\dot{\gamma}\|$. The angular velocity $\omega$ behaves exactly as you would expect in the flat case, but on a curved manifold there is not a geometrically meaningful notion of absolute angle $\theta$, since we cannot compare angles at different points. In coordinates, the resulting system can be written as $$ \ddot{\gamma}^i=a\frac{\dot{\gamma}^i}{\|\dot{\gamma}\|}+\omega R^i{}_j\dot{\gamma}^j-\Gamma^i{}_{jk}\dot{\gamma}^j\dot{\gamma}^k $$ If you're interested in allowing $\dot{\gamma}=0$, you can split $\dot{\gamma}$ into two degrees of freedom $T=\dot{\gamma}/\|\dot{\gamma}\|$ and $v=\|\dot{\gamma}\|$ (with $T$ constrained by $g_{ij}T^iT^j=1$). With a bit of algebra, this simplifies to a system of ODEs with degrees of freedom $(\gamma,T,v)$, where the $v$ equation conveniently decouples from the rest. \begin{align*} \dot{\gamma}^i=&vT^i \\ \dot{T}^i=&\omega R^i{}_j T^j-v\Gamma^i{}_{jk}T^jT^k \\ \dot{v}=&a \end{align*}

There are, of course, other ways to parameterize curves like this. The best known are the Frenet-Serret equations, which use auxiliary vector-valued degrees of freedom rather than a rotation, but this is probably unnecessary unless you are interested in higher dimensions.

Note on computing $R$ in terms of $g$

Throughout all matrices are 2x2. Let $\epsilon=\begin{bmatrix}0&-1\\1&0\end{bmatrix}$ and $g=J^TJ$. Note that for any matrix $A$, we have $A^T\epsilon A=|A|\epsilon$. Let $U$ be a matrix whose columns form an oriented, $g$-orthonormal basis, i.e. $U^TgU=I_2$, $|U|>0$. This implies $U=g^{-1}(U^{-1})^T$ and $|U|=|g|^{-1/2}$. In an oriented orthonormal basis, a $90^\circ$ rotation is just $\epsilon$, so we have \begin{align} R=&U\epsilon U^{-1} \\ =&g^{-1}(U^{-1})^T\epsilon U^{-1} \\ =&g^{-1}|U^{-1}|\epsilon \\ =&\sqrt{|g|}\ g^{-1}\epsilon \end{align}