Here what I have managed to come up with, give it a look. I decided that a Lagrangian approach is best, though I think moving to the Hamiltonian would have some advantages also such as conservation of angular momentum and energy being more explicit. The ellipse' (I worked only in 2D) center of mass can move in the $x$ and $y$ directions and is rotating. Since there is movement in the $y$ direction, there is a need for gravitational potential energy also. My Lagrangian looks like
$$
L = \frac{1}{2}M(\dot{x}^2 + \dot{y}^2) + \frac{1}{2}I\dot{\phi}^2 + Mgy
$$
I take the reference height for potential energy to be $a$ above ground, the semi-major axis. Now there are some serious constraints on this problem, which when combined with this modest Lagrangian make the problem very difficult.
First off, as much as the ellipse rolls (without sliding) so must it cover that much ground. In other words, the arc length traced out by rolling equals the distance in the $x$ direction covered, $f_1 = x - a \, E(\phi,k)=0$ with $E$ the "incomplete elliptic integral of the second kind".
Now motion in the $y$ direction for the center of mass is oscillatory. This is a constraint on it's motion: $f_2 = y + \frac{(a-b)}{2}(1+\cos(\alpha \dot{\phi}t))=0$, with $\alpha$ some possible numerical factor.
Then the equations of motion are
$$
M\ddot{x} = \lambda_1
$$
$$
M\ddot{y} - Mg = \lambda_2
$$
$$
I\ddot{\phi} = -\lambda_1 a \frac{\partial E(\phi,k)}{\partial \phi}
$$
Now
$$
\dot{x} = a\frac{\partial E(\phi,k)}{\partial \phi}\dot{\phi} \implies \ddot{x} = a \left( \frac{\partial^2 E(\phi,k)}{\partial \phi^2}\dot{\phi}^2 + \frac{\partial E(\phi,k)}{\partial \phi}\ddot{\phi} \right)
$$
and
$$
\dot{y} = \frac{\alpha(a-b)}{2} \sin(\alpha \dot{\phi}t) \left( \ddot{\phi}t + \dot{\phi} \right)
$$
so
$$
\ddot{y} = \frac{\alpha (a-b)}{2} \left( \alpha \cos (\alpha \dot{\phi} t) (\dot{\phi} + t \ddot{\phi})^2 + \sin (\alpha \dot{\phi} t)(2\ddot{\phi} + t \dddot{\phi}) \right)
$$
Combining these with the previous equations we get nasty stuff...
$$
M a \left( \frac{\partial^2 E(\phi,k)}{\partial \phi^2}\dot{\phi}^2 + \frac{\partial E(\phi,k)}{\partial \phi}\ddot{\phi} \right) = \lambda_1
$$
$$
\frac{\alpha (a-b)}{2} \left( \alpha \cos (\alpha \dot{\phi} t) (\dot{\phi} + t \ddot{\phi})^2 + \sin (\alpha \dot{\phi} t)(2\ddot{\phi} + t \dddot{\phi}) \right)-g = \frac{\lambda_2}{M}
$$
$$
\ddot{\phi} = -\frac{\lambda_1 a}{I} \frac{\partial E(\phi,k)}{\partial \phi}
$$
This may or may not be right, but I think it is a step in the right direction at least. I hope it helps.
