Coupled ODEs that model a quad rotor I am working on modeling the vibrations of a quad rotor. The arms that support the rotors are fixed to a center plate; that is, it is pretty much a cantilever beam with an end load. Since this is the case, $m_{eq} = m_{arm} + \frac{m_{motor}}{3}$ and $k_{eq} = \frac{3EI}{\ell^3}$ where $E$ is Young's modulus and $I$ is the moment of inertia of the beam, $I = \frac{bh^3}{12}$. The motors are 880 rpm/Kv and the mass of the whole system is denoted $m_t$.
\begin{align}
m_{eq}\ddot{x}_1 + k_{eq}x_1 + c\dot{x}_1 &= F_1(t)\\
m_{eq}\ddot{x}_2 + k_{eq}x_2 + c\dot{x}_2 &= F_2(t)\\
m_{eq}\ddot{x}_3 + k_{eq}x_3 + c\dot{x}_3 &= F_3(t)\\
m_{eq}\ddot{x}_4 + k_{eq}x_4 + c\dot{x}_4 &= F_4(t)
\end{align}
The equations above are the equations of motion for the four arms and rotors. My intention is to determine a damping constant that will dissipate the most of the damping so a video can be taking without the jell-o affect.


*

*What would be the EOM for the entire system?
$$
m_t\ddot{x}_5 = \mbox{???}
$$

*Is it possible to determine a theoretical forcing function for the given motor?

*How would I go about determining the damping constants? Would they be the same for each arm?


If you need a visual aid, it looks similar to this:


Edit:
I am viewing the arms as cantilever beams and $z$ is being used as the displacement of the body connected to the arms.
\begin{align}                                                                   
  m_{eq}\ddot{y}_1 &= k_{eq}(z - y_1) + c(\dot{z} - \dot{y}_1) - m_{eq}y_1 +    
                      F_1(t)\\                                     
  m_{eq}\ddot{y}_2 &= k_{eq}(z - y_2) + c(\dot{z} - \dot{y}_2) - m_{eq}y_2 +    
                      F_2(t)\\                                     
  m_{eq}\ddot{y}_3 &= k_{eq}(z - y_3) + c(\dot{z} - \dot{y}_3) - m_{eq}y_3 +    
                      F_3(t)\\                                     
  m_{eq}\ddot{y}_4 &= k_{eq}(z - y_4) + c(\dot{z} - \dot{y}_4) - m_{eq}y_4 +    
                      F_4(t)\\                                     
  m_b\ddot{z} &= \sum_i\bigl[F_i(t) + k_{eq}(y_i - z) +                         
                 c(\dot{y}_i - \dot{z})\bigr] - m_bz\\              
  m_t\ddot{y}_5 &= \sum_iF_i(t)\tag{quadrotor displacement}\\           
  COM &= \frac{m_bz + \sum_im_{eq}y_i}{m_t}\tag{constraint equation}
\end{align}
 A: 1) Assuming your coordinates are the motions of the ends of the arms, the equations as you've written them don't allow for any center of mass motion, or at any rate any COM acceleration. Let the center-of-mass coordinates of the arm be located a fraction $a$ of the arm length $l$ from the pivot point. For each arm you have:
$$m_\mathrm{arm}(a\ddot{x}_1 + (1-a)\ddot{x}_5) + \frac{I_{arm1}}{l^2}(\ddot{x}_1 - \ddot{x}_5)+ k_{eq}(x_1-x_5)+c(\dot{x}_1-\dot{x}_5)=F_1(t)$$
where $x_5$ is the coordinates of the body. Good old Newtonian summation relates $x_5$ to net force on the joints:
$$m_\mathrm{center} \ddot{x}_5=-\sum_i k_{eq}(x_i-x_5)+c(\dot{x}_i-\dot{x}_5)$$
2) Outside my area of expertise, someone else will have to answer it. I imagine that the relevant part for you would take the form of a noise, but I can't tell you more. 
3) The usual approach is to Fourier transform. However, since the forces are stochastic and you are really interested in $\langle x_5^2 \rangle$, this is more complicated than the normal "replace derivatives with $i\omega$" operation. This post may give you some useful leads. I doubt that you will find a value or set of values for $c$ that will minimize the response to a flat distribution for $F(\omega)$. If $\langle F(\omega) \rangle$ is peaked at a given value, you may be able to suppress oscillations at that value. For example, if you assume each force is just a simple oscillation at $\omega_{rotor}$ with a random phase offset, Fourier transformation and linear algebra will tell you what $c$ minimizes $x_5(\omega_{rotor})$. 
A: From user27118, I was able to grasp a better understanding of what needed to be done to set up the system of ODEs:
\begin{alignat}{2}                                                              
  m_{eq}\ddot{y}_1 &= k_{eq}(z - y_1) + c(\dot{z} - \dot{y}_1) - m_{eq}gy_1 +   
  F_1(t)\\                                                         
  m_{eq}\ddot{y}_2 &= k_{eq}(z - y_2) + c(\dot{z} - \dot{y}_2) - m_{eq}gy_2 +   
  F_2(t)\\                                                         
  m_{eq}\ddot{y}_3 &= k_{eq}(z - y_3) + c(\dot{z} - \dot{y}_3) - m_{eq}gy_3 +   
  F_3(t)\\                                                         
  m_{eq}\ddot{y}_4 &= k_{eq}(z - y_4) + c(\dot{z} - \dot{y}_4) - m_{eq}gy_4 +   
  F_4(t)\\                                                         
  m_b\ddot{z} &= \sum_i\bigl[F_i(t) + k_{eq}(y_i - z) +                         
  c(\dot{y}_i - \dot{z})\bigr] - m_bgz &&{}= 0\\                                                         
  y_5 &= \frac{m_bz + \sum_im_{eq}y_i}{m_t} &&{}= 0                 
\end{alignat}
where I am using $z$ to denote the displacement of the body of the quad rotor.
Using the last equations, we obtain:
\begin{align}                                                                   
  m_{eq}\sum_i\ddot{y}_i                                                          
  &= \frac{m_b + 4m_{eq}}{m_b}\Bigl[\sum_ic\dot{y}_i + k_{eq}y_i\Bigr] +        
     gm_{eq}\sum_iy_i + \sum_iF_i(t)\\                             
  m_{eq}\ddot{y}_i                                                              
  &= k_{eq}\Bigl[y_i + \frac{m_{eq}}{m_b}\sum_jy_j\Bigr] +                      
     c\Bigl[\dot{y}_i + \frac{m_{eq}}{m_b}\sum_j\dot{y}_j\Bigr] - m_{eq}gy_i +  
     F_i(t)                                                  
\end{align}
which leads to 
\begin{align}                                                                   
  m_{eq}\ddot{y}_1                                                              
  &= k_{eq}\Bigl[y_1 + \frac{m_{eq}}{m_b}\sum_jy_j\Bigr] +                      
     c\Bigl[\dot{y}_1 + \frac{m_{eq}}{m_b}\sum_j\dot{y}_j\Bigr] - m_{eq}gy_1 +  
     F_1(t)\\                                                                   
  m_{eq}\ddot{y}_2                                                              
  &= k_{eq}\Bigl[y_2 + \frac{m_{eq}}{m_b}\sum_jy_j\Bigr] +                      
     c\Bigl[\dot{y}_2 + \frac{m_{eq}}{m_b}\sum_j\dot{y}_j\Bigr] - m_{eq}gy_2 +  
     F_2(t)\\                                                                   
  m_{eq}\ddot{y}_3                                                              
  &= k_{eq}\Bigl[y_3 + \frac{m_{eq}}{m_b}\sum_jy_j\Bigr] +                      
     c\Bigl[\dot{y}_3 + \frac{m_{eq}}{m_b}\sum_j\dot{y}_j\Bigr] - m_{eq}gy_3 +  
     F_3(t)\\                                                                   
  m_{eq}\ddot{y}_4          
&= k_{eq}\Bigl[y_4 + \frac{m_{eq}}{m_b}\sum_jy_j\Bigr] +                      
     c\Bigl[\dot{y}_4 + \frac{m_{eq}}{m_b}\sum_j\dot{y}_j\Bigr] - m_{eq}gy_4 +  
     F_4(t)                                                                     
\end{align}
