Context
On the following drawing, a platform is hung from the ceiling not exactly from its centre of gravity. Because of this it can't sustain an arbitrary orientation for long; I want to increase its stiffness with the use of 3 orthogonal fast spinning flywheels, but before building those flywheels I need to know how big they need to be by making a dynamic model…
… And I have been unable to model the gyroscopic effect [Edit: end removed for clarity].
Question
How can I model the 3 flywheels reaction to (themselves and) the gravity torque, be it through a rate increment or a torque (to inject in the plant model)?
Implementing Lelesquiz's answer
Top level equations: Euler-Lagrange
The goal is to solve for each timestep for the angles $(\phi, \theta, \psi)$ using the equations, in inertial frame:
$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\phi}}\right)-\frac{\partial L}{\partial \phi}= (\boldsymbol{M_g}-\sum_{k=1}^3 \boldsymbol{M_{wk}})\cdot \frac{\partial \boldsymbol{\omega_b}}{\partial \dot{\phi}}$$
$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\theta}}\right)-\frac{\partial L}{\partial \theta}= (\boldsymbol{M_g}-\sum_{k=1}^3 \boldsymbol{M_{wk}})\cdot \frac{\partial \boldsymbol{\omega_b}}{\partial \dot{\theta}}$$
$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\psi}}\right)-\frac{\partial L}{\partial \psi}= (\boldsymbol{M_g}-\sum_{k=1}^3 \boldsymbol{M_{wk}})\cdot \frac{\partial \boldsymbol{\omega_b}}{\partial \dot{\psi}}$$
Where $M_g$ is the gravity torque, $M_{wk}$ is the motor torque applied on the flywheel k (opposite of what's applied on the platform), and $\omega_b$ is the platform rate vector in inertial frame. 3 other equations based on $\alpha_k$, the angles of the flywheels, should have been added to completely characterize the system but the rates of the wheels are inputs, measured directly. The torques on the wheels are inputs as well (and the command of the closed loop system).
Knowing that the lagrangian $L=T-U$ in the inertial frame.
Kinetic energy T
$$ T(\phi, \theta, \psi, \dot{\phi}, \dot{\theta}, \dot{\psi}) = \frac{1}{2} \sum_{k=1}^3
\left(\boldsymbol{N_{k,Ri}^T}\boldsymbol{[I_{wO,Ri}]}\boldsymbol{N_{k,Ri}} \right) + \frac{1}{2} \left( m_{b} \boldsymbol{v_{b,Ri}}^T\boldsymbol{v_{b,Ri}}+\boldsymbol{\omega_{b,Ri}^T} \boldsymbol{[I_{bG,Ri}]} \boldsymbol{\omega_{b,Ri}} \right)$$
Where $m_b$ is the mass of the system (with wheels on), $\boldsymbol{v_{b}}$ the velocity vector of the center of gravity of the system, $\boldsymbol{\omega_{b}}$ the system's rate vector (wheels not spinning) and $\boldsymbol{[I_{bG}]}$ the inertia matrix of the system about its center of gravity (which needs to be transported from the guimbal centre 0), $[I_{wO}]$ the inertia matrix of the flywheels about their axis, and $\boldsymbol{N_{k}}$ the rate vector of the flywheels in the inertial frame.
$$\begin{array}{lc}\boldsymbol{[I_{wO,Ri}]}=\boldsymbol{[R_{ib}]}\left[\begin{matrix}J_w & 0 & 0\\ 0 & J_w &0\\0 & 0& J_w\end{matrix}\right] \boldsymbol{[R_{bi}]}
\end{array}$$
$$\boldsymbol{v_{b}}=\boldsymbol{\omega_{b,Ri}}\times\boldsymbol{\Delta_{b,Ri}}$$
$$\begin{array}{lc}
\boldsymbol{[I_{bG,Ri}]}= \boldsymbol{[R_{ib}]}\boldsymbol{[I_{bO,Rb}]}\boldsymbol{[R_{bi}]}+m_b& \left(\begin{matrix}
\Delta_{by,Ri}^2+\Delta_{bz,Ri}^2 & -\Delta_{bx,Ri}\Delta_{by,Ri} & -\Delta_{bx,Ri}\Delta_{bz,Ri} \\
-\Delta_{bx,Ri}\Delta_{by,Ri} & \Delta_{bx,Ri}^2+\Delta_{bz,Ri}^2 & -\Delta_{by,Ri}\Delta_{bz,Ri} \\
-\Delta_{bx,Ri}\Delta_{bz,Ri} & -\Delta_{by,Ri}\Delta_{bz,Ri} & \Delta_{bx,Ri}^2+\Delta_{by,Ri}^2
\end{matrix}\right)
\end{array}$$
$\boldsymbol{\Delta_{b,Ri}}$ being the guimbal centre – center of gravity vector of the system expressed in the inertial frame. See below "Euler angle expressions" for the transformation body-inertial $\boldsymbol{<vector>_i}=\boldsymbol{[R_{ib}]}\boldsymbol{<vector>_b}$.
Potential energy U
$$U(\phi, \theta, \psi) = -m\boldsymbol{g_{Ri}}.\boldsymbol{\Delta_{b,Ri}}$$
The minus sign comes from the fact U is minimum with the guimbal-CG vector pointing downward (dot product positive).
Euler angles expressions
The Euler angles 313 have been selected as a coordinate system, which gives the rotation matrix from inertial frame $(0 \boldsymbol{X_i Y_i Z_i})$ to body frame $(0 \boldsymbol{X_b Y_b Z_b})$:
$$ \begin{array}{lc}
\boldsymbol{R_{RbRi}}(\phi, \theta, \psi)= & \left[\begin{matrix}
cos \psi & sin \psi & 0 \\
-sin \psi & cos \psi & 0 \\
0 & 0 & 1
\end{matrix}\right] \left[\begin{matrix}
1 & 0 & 0 \\
0 & cos \theta & sin \theta \\
0 & -sin \theta & cos \theta
\end{matrix}\right] \left[\begin{matrix}
cos \phi & sin \phi & 0 \\
-sin \phi & cos \phi & 0 \\
0 & 0 & 1
\end{matrix}\right]
\end{array}$$
The rate vector is obtained by summing the projections on each axis of the inertial frame the $\boldsymbol{\dot{\phi},\dot{\theta},\dot{\psi}}$ vectors):
$$\begin{array}{lc} \boldsymbol{w_{b,Ri}}= & \left[\begin{matrix}
\dot{\theta}cos\phi+\dot{\psi} sin\theta sin\phi \\
\dot{\theta}sin\phi+\dot{\psi} sin\theta cos\phi \\
\dot{\phi}+\dot{\psi}cos\theta
\end{matrix}\right]\end{array}$$
Problem
The kinetic energy of the wheels is independent on the Euler angles, so the rates of the wheels have no effect, which is wrong! Since it's absolutely normal that no matter the rotation matrix the norm of the wheels rates vector does not change, the kinetic energy expression or even the top level equations must be wrong.
Best Answer
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
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.