\documentclass[a4paper, 11pt, pdftex]{article}
\usepackage[pdftex]{graphicx}
\usepackage{ae}
\usepackage{latexsym}
\usepackage{rotating}
\usepackage[intlimits, sumlimits]{amsmath}
\usepackage[pdftex, bookmarks]{hyperref}
\pdfcompresslevel=9

\newcommand{\unit}[1]{\textrm{#1}}
\newcommand{\vi}[1]{_{\textrm{#1}}}
\newcommand{\EtaProp}{\eta_{\textrm{prop}}}
\newcommand{\EtaSteller}{\eta_{\textrm{steller}}}

\def\jwwTitle{How to adjust multicopter parameters, multicopter model version 1}
\def\jwwSubject{}
\def\jwwKeywords{multicopter, quadrocopte, CRRCSim}
\def\jwwAuthor{Jens Wilhelm Wulf}

\hypersetup{pdftitle={\jwwTitle},%
    pdfsubject={\jwwSubject},%
    pdfkeywords={\jwwKeywords},%
    pdfauthor={\jwwAuthor}}

\author{\jwwAuthor}
\title{\jwwTitle}
\date{06.12.2009}

\begin{document}

  \numberwithin{equation}{section}
  \maketitle
  \tableofcontents
  \clearpage

  \setlength{\parindent}{0pt}
  \setlength{\parskip}{5pt plus 2pt minus 1pt}
  
  \section{General file format}

    Take a look at {\tt qc01.xml} as an example.

    Hard points/wheels, mass, inertia are just like with fixed wing models
    (see \cite{l:C_FILE}),
    the only difference is about aerodynamic data and controllers.  

    You need to define where propellers are mounted. Describe x- and
    y-offsets to center of gravity and whether props are spinning
    counter-clockwise ({\tt dir="-1"}).

    All motors/propellers are assumed to be of the same type (except for
    direction of rotation) and defined in the power-section as it is done
    for fixed-wing models. Note that the simulation uses one battery for
    every motor, so its capacity defined in the file is smaller than you
    might think. {\tt qc01.xml} does not use the automagic type of description in
    order to allow a more detailed fit to the real model.

    Coordinate system: x positive forward, y positive right, z positive
    down. Rotation rates around these axes are $p$, $q$, $r$.

    In order to adjust the CRRCSim model to my real one, I used the
    following steps. Beware: this is NOT easy, and some parts even plain
    dangerous. I guess that people who are able to do this are able to
    do the math themselves and much more.

  \section{Statics}
  
    Hover data measured on the real copter and used for some of the calculations here:
    
    \begin{tabular}{r|r|r|r|r}
      PWM-ratio &  battery voltage &      rpm & mass  &  $M$ (calculated, see below)\\
      \hline
      80/255    &  11,2 V          & 2768 rpm & 405 g &  0.02 N\\
      96/255    &  10,8 V          & 3100 rpm & 505 g &  0.02494 N\\
      110/255   &  10,6 V          & 3397 rpm & 612 g &  0.03022 N\\
    \end{tabular}

    First of all, prop positions, mass and inertia should be quite accurate.
    All this can simply be measured and calculated.
    
    Describe your power plant.
    
    Set {\tt yaw.dist} and {\tt roll.dist} to zero.
    
    \subsection{Propeller}
    
    Weight is 400\,g, so every motor/prop needs to create
    $F = 0.1\,\unit{kg} \cdot 9.81\,\unit{m}/\unit{s}^2=0.981\,\unit{N}$.
    
    Take a look at the source code to see how propeller force is calculated:
    {\tt src/mod\_fdm/power/propeller.cpp}. When hovering, this is
    \begin{subequations}
    \begin{eqnarray}
      \EtaProp &=& 0.65\\
      \varrho &=& 1.225\,\unit{kg}/\unit{m}^3\\
      F &=& \pi \cdot 0.25 \cdot D^2 \varrho \left(H n\right)^2 \cdot 1/2\, \EtaProp \label{f:PROP_F}
    \end{eqnarray}
    \end{subequations}
    I know that my setup needs around 2800\,rpm to hover, $D$=0.235\,m.
    Using these values and equation~\eqref{f:PROP_F} results in $H$=0.162\,m.
    
    I recalculated using all the numbers from the table above, the mean
    value is $H$=0.1644\,m.
    
    You may ask about why I didn't just use the propeller's $H$ as defined by
    the manufacturer. First, my prop is a modified standard one, so usual
    rules don't apply. Secondly, I often read that manufacturers $H$ is not
    accurate.
    
    Set {\tt engine.log} to 1 and hover at fixed altitude using an attitude
    controlled setup. CRRCSim will write lots of data to stdout: one line
    per timestep and motor. Redirect this output to a file (is automatically
    done on windows) in order to not slow your PC down.
    
    Testing shows that the copter hovers at 2800\,rpm now. Fine.
    
    Motor parameters need to be adjusted. Also see \cite{l:C_POWER}. However,
    in this case I have some data at three operating points, but it does not
    include motor current. Going back to the propeller code, lets calculate
    the torque needed. Equations found in the code are:

    \begin{subequations}
    \begin{eqnarray}  
      P &=& F_X \, V_p / (2 \, \EtaProp) \\
      P &=& F_X \, H \, n / (2 \, \EtaProp) \\
      M &=& P / \omega \\
      M &=& P / (2 \, \pi \, n) \\
      M &=& F_X \, H \, n / (2 \, \EtaProp \, 2 \, \pi \, n)\\
      M &=& F_X \, H / (4 \, \pi \, \EtaProp)
    \end{eqnarray}
    \end{subequations}
      
    With $F_X = m/4 \cdot 9.81\,\unit{m}/\unit{s}^2$ this results in $M$ as
    written into the table above.
    
    \subsection{Engine}
    
    So let's take a look at the engine's equations.
    
    {\tt src/mod\_fdm/power/engine\_dcm.cpp} reveals
    \begin{subequations}
    \begin{eqnarray}  \label{f:MOTOR}
      \EtaSteller &=& 0.95 \\
      U\vi{K} &=& \textrm{throttle} \cdot U \cdot \EtaSteller \\
      U\vi{Gen} &=& \omega \cdot k\vi{M} \\
      I\vi{M} &=& (U\vi{K} - U\vi{Gen}) / R\vi{I} \\
      M\vi{M} &=& k\vi{M} \cdot I\vi{M} - k\vi{r} \cdot \omega - M\vi{r} \\
      M\vi{r} &=& k\vi{M} \cdot I_0 \\
      k\vi{r} &=& 0 \\
    \end{eqnarray}
    \end{subequations}
    Which boils down to
    \begin{equation}
      M\vi{M} = k\vi{M} \cdot (\textrm{throttle} \cdot U \cdot 0.95 - \omega \cdot k\vi{M}) / R\vi{I} - I_0 \cdot k\vi{M}
    \end{equation}
    for us. There are three unknown parameters here and luckily I have three sets of data, 
    which solves to $k_M=0.0131$ Vs, $R_I=0.1839\,\Omega$, $I_0=-4.013$\,A. Of course a negative 
    $I_0$ is bullshit. Something wrong with the simulation equations? However, 
    I continue and force $I_0$ to be zero and solve again. Due to three data sets, I get 
    three solutions:
    
    $k\vi{M}$ = 0.0081483226538728\,Vs, $R\vi{I}$ = 0.39817956295345\,$\Omega$ \\
    $k\vi{M}$ = 0.0085130006948808\,Vs, $R\vi{I}$ = 0.37102864543948\,$\Omega$ \\    
    $k\vi{M}$ = 0.0089455381882278\,Vs, $R\vi{I}$ = 0.34435633250623\,$\Omega$ \\
    
    Now motor current and voltage used are right. Don't be surprised because current 
    is high according to log output: this is the current flowing through the motor, 
    not the one which has to be delivered by the battery. The latter is lower 
    (we use PWM here).

  \section{Dynamics}
  
    I measured the real one's behaviour using step responses. This is
    dangerous and you might need special firmware! However, I did it and
    will use the results to tweak simulation parameters.

    I used a firmware which records raw sensor data of every control cycle.
    In this case motor voltage (as commanded to the controller) and rate
    gyro output is important. After turning a switch, it turns of the
    controller of the axis to be examined for a defined time (about 500\,ms in
    my case). Afterwards, commands a predefined voltage to this axis for a
    defined time (some hundred milliseconds). By turning the switch back,
    the process can be interrupted at any time.

    Of course you can only do this on a calm day. If you're out of luck, the
    copter already leaves its attitude while the commanded voltage is zero.
    However, it will do this once the predefined voltage is applied. React
    quickly, it will go down and crash otherwise.

    I tried getting the parameters from data gathered while the control loop
    was closed. This is not dangerous at all, but only gave sane parameters
    for rotation about the z-axis.

    For x-axis, my controller commands voltage like this:
    \begin{subequations}  \label{f:U}
    \begin{eqnarray}
      u\vi{MotLeft}  &=& u_0 + u\vi{ControllerX}\\
      u\vi{MotRight} &=& u_0 - u\vi{ControllerX}
    \end{eqnarray}
    \end{subequations}
    y-axis is the same.

    So the controller is turned of for some time, $u\vi{ControllerX}$=0.
    Afterwards $u\vi{ControllerX}$=const for some time. Rotation rate around
    x-axis is recorded, let's call it $p$. For shortness, $u\vi{ControllerX}$=$u$.

    I searched for parameters to make the following equation fit measurements (in laplace domain):
    \begin{equation}    
      s \, p = a \, u/(\tau\,s+1)^2 - b \, p
    \end{equation}    
    For x and y I found $a$=20\,rad/s$^2$/V, $b$=0, $\tau$=15\,ms.
    
    For z I found $a$=5.3\,rad/s$^2$/V, $b$=0.5/s, $\tau$=5\,ms. $\tau$ is not that important and may be wrong here.

  \subsection{Motor}

    There is something undefined up to now: inertia of motor and propeller.
    I approximate that from the x axis transfer function. The motors
    acceleration is $M\vi{M} / J\vi{M}$; using equations~\eqref{f:MOTOR}:
    \begin{equation}    
      \frac{d\,\omega}{dt} = k\vi{M}\,\left( U - \omega \, k\vi{M} \right) / (R\vi{I}\,J\vi{M})
    \end{equation}
    
    Bringing that from time domain to laplace and solving for $\omega/U$ gives
    
    \begin{equation}    
      \omega/U = 1/k\vi{M}^2 \, \frac{1}{(R\vi{I} \, J\vi{M} / k\vi{M}^2) \, s + 1}
    \end{equation}    
    
    So there is a PT1 (I guess its called first order time lag in 
    english) with $\tau=R\vi{I} \, J\vi{M} / k\vi{M}^2$.

    I assume the motor speed as proportional to the copter's rotational
    acceleration. But the motor's transfer function contains a PT1 while my
    copter's rotation transfer function (see above) contains a PT2. Well, I
    just say that a PT2 with $\tau\vi{pt2}$ is similar to a PT1 with
    $\tau\vi{pt1}=\tau\vi{pt2} \cdot 1.7$. Therefore, 
    $15\,\unit{ms} \cdot 1.7=\tau \cdot 1.7=R\vi{I} \, J\vi{M}/k\vi{M}^2$.
    Solving this for $J\vi{M}$ gives $J\vi{M}=5\cdot 10^{-6} \, \unit{kg} \, \unit{m}^2$.

    Of course, all this assumed that only the motor's inertia is active
    here, but there is the propeller's inertia, too. However, for the
    simulation, only their sum is important, therefore I set {\tt
    propeller.J=0}
    and {\tt engine.J\_M=5E-6}.

    The sum of both inertias could of course also be measured by applying a
    fixed voltage to the motor and recording its rotational speed. However,
    when getting to much details from the real system very accurately, keep
    in mind that CRRCSim doesn't model everything that accurately -- so
    using 'wrong' values might result in a more realistic simulation. In the
    case of a multicopter this is very important, because for example the
    time lag from measuring a rotation rate, calculating the controller,
    transferring a new setpoint to the motor's speed controller, until this
    speed controller applies the new voltage to the motor does depend on
    software heavily. And it depends on software of at least two different
    embedded electronics usually (main flight control and brushless
    controller). Furthermore, different copter firmwares use different
    filters here and there and so on...
  
    So in this regard CRRCSim's multicopter model can't realistically model
    every real one. It would be possible to do, but I guess noone would be
    willing to find all those parameters...

  \subsection{Aerodynamics}
  
    For yaw I found $b$=0.5/s.
    
    The simulator uses
    \begin{subequations} 
    \begin{eqnarray}
      \frac{d\,r}{dt} &=& M\vi{z} \, \frac{I\vi{xx}}{I\vi{xx}\,I\vi{zz} - I\vi{xz}^2}\\
      M\vi{z}         &=& - \texttt{yaw.damp1} \cdot r - \texttt{yaw.damp2} \cdot |r| \, r
    \end{eqnarray}
    \end{subequations}
    as can be seen in {\tt EOM01::ls\_accel()} and 
    {\tt CRRC\_AirplaneSim\_MCopter01::aero()}.
    
    Therefore,
    \begin{equation}    
      \texttt{yaw.damp1} = b \, \frac{I\vi{xx}\,I\vi{zz} - I\vi{xz}^2}{I\vi{xx}}
    \end{equation}            
    which results in {\tt yaw.damp1=0.00341, yaw.damp2=0}.
    
    Although I found $b$=0 for nick and roll, {\tt roll.damp1=0.0009}, {\tt roll.damp2=0}
    because the simulation likes at least some damping.        
    
    I have not thought about something to measure {\tt speed.damp}, so this
    is done by try and error.

  \section{Controllers}
  
    Well, we're not through yet. Physics are defined now, but the control
    loop of the copter is not\dots
    
    Looking at the source, any controllers output means how much voltage is
    applied in the sense of equation~\eqref{f:U} with 
    \begin{equation}
      u\vi{ControllerX} = \textrm{URel} \cdot U\vi{Rel} = \textrm{URel} \cdot U\vi{Battery0} \cdot 0.7
    \end{equation}
    So if the controller output is one, the voltage $u\vi{ControllerX}$ that
    will be applied is 0.7 times the voltage of the fully charged battery.
    The controllers action will be independent of battery voltage as
    long as it is at least 0.7 times the full charge voltage and the
    controller output is smaller than one.

    \subsection{Controller type "Omega"}

    It gives rate controlled behaviour. Because it does not know about the
    scaling with $U\vi{Battery0} \cdot 0.7$, its parameters depend on this
    value. Using variable names for roll axis, with $p\vi{sp}$ being the
    setpoint and $p\vi{pv}$ being the process value, its output is
    \begin{subequations}
    \begin{eqnarray}
      p\vi{d}           &=& p\vi{sp} - p\vi{pv}\\
      u\vi{ControllerX} &=& U\vi{Battery0} \cdot 0.7 \left( K\vi{P} \, p\vi{d} + K\vi{D} \, \frac{d\,p\vi{pv}}{dt} + K\vi{I} \, \int p\vi{d} \right)
    \end{eqnarray}
    \end{subequations}
    
    Compared to the settings of my real copter I had to halve $K\vi{D}$,
    otherwise the simulation got unstable.
    
    To the controller, maximum stick input is 0.5 (just like every other
    signed stick input in CRRCSim) and $p\vi{sp} = {\tt stick} \cdot ${\tt
    scale}. If {\tt scale\_roll\_pitch.exp} $>$ 1, a special scaling is used for
    roll and nick to allow flips and loopings:    
    \begin{equation}
      p\vi{sp} = \texttt{scale} \left( \texttt{a} \cdot \texttt{stick}
                 + \texttt{b} \cdot \texttt{stick}^\texttt{exp}
                 \right)
    \end{equation}
    Here {\tt a}, {\tt b} and {\tt exp} are from the {\tt
    scale\_roll\_pitch} section of the input file.
        
              
  \clearpage
    
  \begin{thebibliography}{99}
  
    \bibitem{l:C_FILE}
      CRRCSim fixed wing file format description,
      CRRCSim installation: {\tt documentation/file\_format/index.html}
      
    \bibitem{l:C_POWER}
      CRRCSim power and propulsion system description,
      CRRCSim installation: {\tt documentation/power\_propulsion/power\_propulsion.html}
    
  \end{thebibliography}

\end{document}