Since in my last post I made a lot of mistakes, I'll try to improve that answer starting from scratch.
First of all: lagrangian coordinates.
It first seemed clear that two coordinates suffice in describing the body orientation. Of course this is false, since in order to describe every possible orientation one needs three coordinates. We'll call them $(\theta, \phi, \psi)$.
Any more coordinates? I don't think so. The only external forces acting of the body are the rope tension and gravity. They are both directed along the vertical, so no oscillations should develop.
Next: rotation description. Refer to this link. You could pick Euler angles or Tait-Bryan, your choice. Given the angles, you then compute the rotation matrix. With the matrix you can compute any vector rotation regarding the body. Will call this $R(\theta (t), \phi (t), \psi (t)) = R(t)$.
Body description: your model depends on many quantities. They essentially are
- Initial body orientation and initial angular velocity
- Initial relative positions of the three gyroscopes
- Initial orientations of the gyroscopes and angular velocities
- Initial center of gravity position
- $\theta(t), \phi(t), \psi(t)$
Given any of these vectors $\vec{x}_0$, you can compute the corresponding vector at time $t$ as $\vec{x}(t) = R(t) \cdot \vec{x}_0 = R(\theta, \phi, \psi) \cdot \vec{x}_0 = \vec{x}(\theta, \phi, \psi)$.
What's next? Compute the kinetic energy of the body. It is given by
$$ T(\theta, \phi, \psi, \dot{\theta}, \dot{\phi}, \dot{\psi}) = \frac{1}{2} I_1 \omega_1^2(\theta, \phi, \psi) + \frac{1}{2} I_2 \omega_2^2(\theta, \phi, \psi) +\frac{1}{2} I_3 \omega_3^2(\theta, \phi, \psi) + \frac{1}{2} I_{TOT} \omega_{body}^2(\theta, \phi, \psi, \dot{\theta}, \dot{\phi}, \dot{\psi}) $$.
Again note that the gyroscopes' angular velocities are constant in magnitude (if you admin no friction and only forces normal to the axis, but they rotate with the body: their relative angle to the body has remain fixed.
Finally compute the potential energy as
$$ U(\theta, \phi, \psi) = mgh(\theta, \phi, \psi)$$
The total lagrangian is
$L = T-U$.
Now you should be able to apply the Euler Lagrange equations and solve numerically. Good luck! You might need to help yourself with mathematica or any other sort of algebra package in order not to damage your brain with the derivatives and scalar products.
There's a constructive proof that can be understood intuitively. I assume z is vertical and y is forward/backward. You rotate the object about the z-axis until the top is somewhere on the xz plane i.e. y=0. This makes the top of the object point perpendicular to the y-axis, so you can rotate about the y-axis until it points up. And now you just need to rotate it on the z-axis until forward points in the right direction.
If you want to know what the original rotation was, just take the opposite of each of those steps and put them in reverse order. For example, if you rotated 10°, -30°, 50°, then to get the original rotation from the correct rotation it's just -50°, 30°, -10°.
This sequence of operations can always be done regardless of how an object is oriented. It's not necessarily unique. If up is already pointing up, then it will be on the correct plane no matter how much you rotate about the z-axis. But the point is that there is some rotation about the z-axis that leaves it on the correct plane, which there is.
Edit:
If you want something strictly more mathematical, suppose you have an orthogonal frame of vectors, $x = (x_1,x_2,x_3), y = (y_1,y_2,y_3),$ and $z = (z_1,z_2,z_3)$.
$z$ is the top of the object, so first we rotate it so $z'_2 = 0$. We just rotate it by $\arctan\frac{z_2}{z_1}$, and it ends up pointing at $z' = \left(\pm\sqrt{z_1^2+z_2^2},0,z_3\right)$. And rotate it another $180^\circ$ if that got a negative instead of a positive, so it's $z' = \left(\sqrt{z_1^2+z_2^2},0,z_3\right)$. $\arctan\frac{z_2}{z_1}$ is only really undefined in the case of $z_1=z_2=0$, in which case don't rotate it at all.
Now we rotate on the $y$-axis by $\arctan\frac{\sqrt{z_1^2+z_2^2}}{z_3}$ and we get $z'' = \left(0,0,\sqrt{z_1^2+z_2^2+z_3^2}\right)$ which must be $(0,0,1)$ since it's a unit vector. Again, if it's undefined, we don't need to rotate it.
Since $y \perp z$ and we're only rotating, $y'' \perp z''$.
From there, we know $y''_3 = 0$, so we have $y'' = (y''_1,y''_2,0)$.
Just rotate it on the $z$-axis by $\arctan\frac{y''_1}{y''_2}$ plus an extra $180^\circ$ if it faces the wrong way, and we get $y''' = \left(0,\sqrt{y''^2_1+y''^2_2},0\right)$. And that's just $(0,1,0)$, since it's a unit vector.
In this case $y''_1$ and $y''_2$ can't both be zero, so we don't have to worry about $\arctan$ being undefined at all.
Since we rotated on the $z$-axis, and $z''$ was already on the $z$-axis, $z''' = z''$.
All we have left is $x$. Since $x = y \times z$, and we're only rotating, $x''' = y''' \times z'''$.
$x''' = (0,1,0) \times (0,0,1)$
$x''' = (1,0,0)$
Best Answer
First Thought (probably not the fastest)
Let us assume you have a vector space in $R^{3}$ with a quaternion defined as:
$$ \mathbf{q} = q_{F}^{*} \ q_{B} \\ = a + b \hat{\mathbf{x}} + c \hat{\mathbf{y}} + d \hat{\mathbf{z}} $$
where $(a, b, c, d)$ are the Euler parameters and $(\hat{\mathbf{x}}, \hat{\mathbf{y}}, \hat{\mathbf{z}})$ defines the reference unit basis set.
If we define the axis of rotation as $\mathbf{n}$ and the angle through which we rotate as $\zeta$, then the Euler parameters are defined as:
$$ a = \cos{\left( \frac{\zeta}{2} \right)} \\ b = n_{x} \ \sin{\left( \frac{\zeta}{2} \right)} \\ c = n_{y} \ \sin{\left( \frac{\zeta}{2} \right)} \\ d = n_{z} \ \sin{\left( \frac{\zeta}{2} \right)} \\ $$
Thus, if you know $\mathbf{q}$, or rather $(a, b, c, d)$, you can find $\mathbf{n}$ and $\zeta$. Once you know the axis of rotation and the angle of rotation, you can determine the Euler angles. First we define the cross product matrix as:
$$ \left[ \mathbf{n} \right]_{x} = \left[ \begin{array}{ c c c } 0 & - n_{z} & n_{y} \\ n_{z} & 0 & - n_{x} \\ - n_{y} & n_{x} & 0 \end{array} \right] $$
and the outer product of $\mathbf{n}$ with itself given by:
$$ \left[ \mathbf{n} \otimes \mathbf{n} \right] = \left[ \begin{array}{ c c c } n_{x} \ n_{x} & n_{x} \ n_{y} & n_{x} \ n_{z} \\ n_{y} \ n_{x} & n_{y} \ n_{y} & n_{y} \ n_{z} \\ n_{z} \ n_{x} & n_{z} \ n_{y} & n_{z} \ n_{z} \end{array} \right] $$
Then we can define the rotation matrix as:
$$ \overleftrightarrow{\mathbf{R}} = \cos{\zeta} \ \overleftrightarrow{\mathbf{I}} + \sin{\zeta} \ \left[ \mathbf{n} \right]_{x} + \left( 1 - \cos{\zeta} \right) \ \left[ \mathbf{n} \otimes \mathbf{n} \right] $$
where $\overleftrightarrow{\mathbf{I}}$ is the unit or identity matrix.
Second Thought (probably faster/easier)
An easier method is to follow the procedure given here. Following that procedure, we define:
$$ \alpha = \frac{ 2 \left( a \ b + c \ d \right) }{ 1 - 2 \left( b^{2} + c^{2} \right) } \\ \beta = 2 \left( a \ c - d \ b \right) \\ \gamma = \frac{ 2 \left( a \ d + b \ c \right) }{ 1 - 2 \left( c^{2} + d^{2} \right) } $$
which gives us the Euler angles:
$$ \phi = \tan^{-1}{ \alpha } \\ \theta = \sin^{-1}{ \beta } \\ \psi = \tan^{-1}{ \gamma } $$
Since you already have $\mathbf{q}$ and you can numerically/analytically determine $(\phi, \theta, \psi)$, then I would just take the time derivative of each of these angles to find $(\dot{\phi}, \dot{\theta}, \dot{\psi})$ rather than using the angular velocities.