newtonian mechanics - Defy gravity torques with gyroscopes?

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{_i}=\boldsymbol{[R_{ib}]}\boldsymbol{_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.