Emmanuel Roussel
email: emmanuel.roussel@isl.eu

Vincent Gassmann

Edouard Laroche

Accuracy-simplicity trade-off for small-scale helicopter models: a comparative study based on flight data

Keywords: Unmanned Aerial Vehicle, Helicopter, Experimental Model Validation, Modeling, System identification, Parameter estimation

published or not. The documents may come

Introduction

A model is a simplified representation of a system, and is necessarily partial and incomplete. Many different models of a particular system can be elaborated, with various structures, complexities and fidelity levels. Physicsbased phenomenological models based on the six-degreeof-freedom (6-DOF) equations are often preferred over black-box models for aircrafts, because their parameters usually have a physical meaning and can be estimated and verified by different means [START_REF] Jategaonkar | Flight vehicle system identification: a time domain methodology[END_REF][START_REF] Klein | Aircraft system identification: theory and practice[END_REF]. Moreover, the literature offers comprehensive mathematical descriptions of the underlying physics and aerodynamics effects and, as the resulting structure comes directly from the physical phenomena, they generally offer large validity domains (see for example [START_REF] Seddon | Basic helicopter aerodynamics[END_REF][START_REF] Leishman | Principles of helicopter aerodynamics[END_REF]). The main drawback of this approach arises in many cases from the complexity of the resulting model which may not be suitable for the expected use. Too simple models cannot reach the required level of fidelity (too many phenomena are neglected), whereas too complex models result in high computational burden, may have identifiability issues (partly due to a higher number of parameters to be estimated), and may not suit for modern controller design.

Consequently, during the modeling work, a trade-off has to be made between complexity and accuracy. Reasona-ble assumptions are required, based on a priori knowledge of the vehicle dynamics and on experimental data. The choice of the model complexity obviously depends on the targeted application (e.g. simulation or control design purposes), but it is often difficult to distinguish a priori 30 between insignificant and prevailing phenomena. In this paper an insight on this particular point is provided, for miniature helicopters, based on experimental data. The 6-DOF equations are quite standard and are used to describe both the translational and rotational motion 35 of a generic rigid-body in space [START_REF] Zipfel | Modeling and simulation of aerospace vehicle dynamics[END_REF]. However, the forces and moments acting on the body are generally much more complex to model, since they entirely depend on the geometry of the vehicle [START_REF] Hurzeler | Modeling and design of unmanned rotorcraft systems for contact based inspection[END_REF]. In the particular case of miniature coaxial helicopters, various levels of complexity are used to 40 model the forces and moments acting on it. The resulting model is sometimes used directly for nonlinear control design, e.g. by Drouot et al. in [START_REF] Drouot | Hierarchical backstepping-based control of a gun launched mav in crosswinds: Theory and experiment[END_REF] with a backstepping approach, but a linearized version is often derived to allow the use of robust control tools, e.g. by Schafroth et al. in [START_REF] Schafroth | Modeling, system identification and robust control of a coaxial micro helicopter[END_REF], where an attitude and heave controller is designed from this linearized model. In these two references, only the main forces and moments are required in the model. In [START_REF] Koehl | Modeling 735 and identification of a launched micro air vehicle: design and experimental results[END_REF], Koehl et al. consider more complex translational dynamics and take into account forces and moments induced by 50 translations in the air (drag on the body) for position and attitude control purposes. A similar model complexity is used by Zarovy and Costello in [START_REF] Zarovy | Extended state observer for helicopter mass and center-of-gravity estimation[END_REF] to design an observer for mass and center-of-gravity estimation, or by Wang et al. in [START_REF] Wang | Flight dynamics modeling of coaxial rotorcraft UAVs[END_REF] and [START_REF] Fankhauser | 745 Modeling and decoupling control of the coax micro helicopter[END_REF], for rotational dynamics identification.

All in all, an essential question is: how far should one go in the description of these forces and moments and in the accuracy of the modeling?

In order to highlight the relevant physical phenomena characterizing miniature helicopters and coaxial rotor configuration in particular, three nonlinear models with increasing levels of complexity and one linear model are considered. The phenomena under study are selected among those most seen in the literature, and they are all based on the 6-DOF equations of motion. The first model is similar to the model described in [START_REF] Drouot | Hierarchical backstepping-based control of a gun launched mav in crosswinds: Theory and experiment[END_REF], [START_REF] Schafroth | Modeling, system identification and robust control of a coaxial micro helicopter[END_REF] and [START_REF] Gassmann | Control of a 750 gun-launched MAV for scene observation[END_REF], and uses a description of the major forces and moment, i.e. caused by the gravity, upper and lower rotor thrust, rotor drag, and swashplate orientation. The second one includes refinements in the expression of the aerodynamic effects acting on the translational dynamics and is similar to the one used in [START_REF] Koehl | Modeling 735 and identification of a launched micro air vehicle: design and experimental results[END_REF], [START_REF] Wang | Flight dynamics modeling of coaxial rotorcraft UAVs[END_REF] and [START_REF] Fankhauser | 745 Modeling and decoupling control of the coax micro helicopter[END_REF]. The third model includes additional effects due to the position of the center of pressure (CoP), and is used for example in [START_REF] Zarovy | Extended state observer for helicopter mass and center-of-gravity estimation[END_REF]. Finally, a linear model is derived from the first model, used by many authors, e.g. in [START_REF] Valavanis | Advances in unmanned aerial vehicles: state 755 of the art and the road to autonomy[END_REF]. The parameters of these four models are then estimated. The models are then compared and their respective fidelity is evaluated against data recorded during a manual flight. A motion capture system provides the "ground-truth" for the position of the helicopter. Performance for both rotational and translational dynamics is evaluated, and allows an informed trade-off between simplicity of the model structure and accuracy of the prediction, depending on the targeted use.

Time-domain identification techniques are used in this study. Frequency-domain techniques have also been used to successfully identify accurate models of miniature helicopters. For instance, Mettler et al. [START_REF] Mettler | System identification modeling of a small-scale unmanned rotorcraft for flight control design[END_REF] used them to identify linear models of miniature helicopters, specific to a flight condition (around hover, or for slow-speed cruise flight). The paper is organized as follows. In Section 2, the structure for the three nonlinear models is described, while the specific aspects of each model are detailed in Section 3. The helicopter platform used for flight data recording is also presented. The contribution of the paper is then twofold. First, in Section 4, an identification procedure is described and applied on the models. Secondly, the capability of these models to mimic the real system is evaluated using flight data and compared in Section 5, thus allowing the choice of an adequate model for a given aim, i.e. selecting the adequate complexity.

Non-linear dynamic model structure

In this section, the coaxial helicopter platform is first described. Its control and measurement signals are presented. The nonlinear equations governing the motion of an aircraft in a three-dimensional space are then described. They are shared by the nonlinear models considered in this study, and take as input the forces and moments acting on the helicopter, which are modeled in a second 110 part (Section 3).

Description of the coaxial helicopter platform

Description of the structure

Figure 1 shows the platform that has been used to collect flight data for this study. It is a low-cost miniature 115 coaxial helicopter similar in its working principle to the rotor-head and in size to the commercial "Lama v4" or "Big-Lama" helicopters (see also [START_REF] Schafroth | Modeling, system identification and robust control of a coaxial micro helicopter[END_REF] and [START_REF] Fankhauser | 745 Modeling and decoupling control of the coax micro helicopter[END_REF]), except that no stabilizer bar is used: the pitch of the upper rotor blades is fixed (so called "flybarless" structure). This allows 120 faster and more aggressive flight maneuvers but decreases the stability, therefore pitch and yaw axes stabilizing controllers are necessary (Section 2.4).

On the contrary, the cyclic pitch of the lower rotor blades is driven by two servomotors through the swashplate.

125 By this mean, the direction of the thrust generated by the lower rotor can be changed, which allows to control the roll and pitch angles of the helicopter. The two concentric rotors, driven by two brushless motors and rotating in opposite directions, allow yaw-axis control without the need 130 for a tail rotor (by using the motors differential). Finally, as the collective pitch cannot be changed, heave motion is also achieved by varying the speed of the two rotors.

Onboard hardware and software

The total mass is 325 g and the rotor span is 35 cm. The instrumentation system allows flight data recordings on a µSD card (including sensor measurements, instructions from the remote control and state estimates). From these forces and moments, the 6-DOF motion model, in red, provides the states of the helicopter. The sensors dynamics are neglected in this study.

Control signals and measurements Control signals

The inputs of the system, given through servomotor and 160 motor controllers, are fairly common for aerial vehicles:

• u lat : lateral or aileron input (first swashplate angle setpoint), • u lon : longitudinal or elevator input (second swashplate angle setpoint),

165

• u thr : throttle input (total motor speed setpoint),

• u ped : anti-torque or rudder pedal input (motor speed differential setpoint). They are normalized between [-1, 1].

Inputs of the model of forces and moments

Input u ped acts on the motors differential, while u thr acts on the sum of their rotation speed. As a first approximation, the rotation speed of the upper and lower rotors (Ω u and Ω l , respectively) can be expressed as:

Ω u Ω l = K mot u thr + u ped u thr -u ped (1)
Similarly, a gain K ser links the inputs [u lat , u lon ] T to the swashplate angles [δ lat , δ lon ] T :

δ lat δ lon = K ser u lat u lon (2)
These equations (1) and (2) are described by the blue blocks in Figure 2.

Outputs of the system

The output is y

= [ξ I m , v B m , η m , ω B m ]
T , where:

• ξ I = [x, y, z]
T is the position of the center of gravity of the UAV w.r.t. the ground. It is expressed in an inertial frame I = {O, x I , y I , z I }, attached to the ground, and whose axes are oriented North, East and Down. (see Figure 3). • ω B = [p, q, r] T is the body-axes angular rate vector. Figure 3 shows the description of the helicopter attitude by the use of the Euler angles, which links the frames I and B. The equipment used to record the outputs y is 190 described in section 4.3.

Aircraft equations of motion

Considering the helicopter as a rigid body with a constant mass m, the 6-DOF model can be used to describe both its translational and rotational motions in three-195 dimensional space. The equations are composed of a dynamic part based on Newton's and Euler's equations on the one side, and, on the other side, on a kinematic part which allows the determination of the vehicle attitude (Euler angles) and position in the inertial frame. A detailed 200 explanation of this motion model, widely used for aircraft systems, is given by Zipfel [START_REF] Zipfel | Modeling and simulation of aerospace vehicle dynamics[END_REF].

Considering the state vector x given in (3), the 6-DOF model is written as a set of coupled differential equations (4-7),

205 x = [ξ I , v B , η, ω B ] T = [x, y, z, u, v, w, φ, θ, ψ, p, q, r] T (3) ξI =   c θ c ψ s φ s θ c ψ -c φ s ψ c φ s θ c ψ + s φ s ψ c θ s ψ s φ s θ s ψ + c φ c ψ c φ s θ s ψ -s φ c ψ -s θ s φ c θ c φ c θ   v B (4)
  0 -r q r 0 -p -q p 0   v B + 1 m   F x F y F z   +   -s θ s φ c θ c φ c θ   g (5) η =   1 t θ s φ t θ c φ 0 c φ -s φ 0 s φ /c θ c φ /c θ   ω B (6) ωB =        0 I yy I xx .r - I zz I xx .q - I xx I yy .r 0 I zz I yy .p I xx I zz .q - I yy I zz .p 0        ω B +        1 I xx 0 0 0 1 I yy 0 0 0 1 I zz          L M N   (7)
where

f B = [F x , F y , F z ] T and τ B = [L, M, N ]
T are respectively the total force and the total moment acting on the vehicle, expressed in the body frame. I = diag(I xx , I yy , I zz ) is the inertia matrix, considered as diagonal given the symmetries of the UAV. Symbols s α , c α and t α are used to denote sin(α), cos(α) and tan(α), respectively.

Attitude control

The helicopter is unstable in open-loop, especially without a stabilizer bar, and a controller is needed to stabilize the system and allow the pilot to do manual flights for data collection. The controller is kept as simple as possible to ease the identification procedure (detailed in Section 4). Its structure is given in Figure 4. A hierarchical structure is chosen, which is often used for the control of small 220 UAVs (see for example [START_REF] Dong | Modeling and control of a quadrotor UAV with aerodynamic concepts[END_REF] and [START_REF] Drouot | Hierarchical backstepping-based control of a gun launched mav in crosswinds: Theory and experiment[END_REF] that use this control structure): an inner loop acts on the angular rates and an outer loop acts on the attitude. It has been implemented on the embedded electronics for flight tests. The control signals are given as follows:

       u lat = k 2 (k 1 (φ * -φ m ) -p m ) u lon = -k 2 (k 1 (θ * -θ m ) -q m ) u thr = k 3 Ω * u ped = k 4 (r * -r m ) (8)
Scalar parameters k 1 to k 4 allow to tune the closedloop behavior. Notice that this controller is purely static, thus the closed-loop and open-loop systems have the same order. The reference input of the closed-loop system is r = [φ * , θ * , Ω * , r * ] T , and is given by the pilot using a 4-230 channel remote control.

Aerodynamic model equations

The set of equations (4-7) describes the motion of a generic rigid-body in space. Whereas these equations are not specific to a particular flight vehicle, the fidelity of the 235 model entirely depends on the description of forces f B and moments τ B acting on the vehicle.

In this section, three models with increasing complexity are described, called M 0 to M 2 . Equations ( 1), ( 2), (4-7), and ( 8) are common for all three models, but for each 240 successive model, additional forces and moments are taken into account. Then, a linearized model M lin is computed. These models are often used in the literature for different target applications, and they will be compared in the last section (section 5).

Forces and moments generated by the rotors

The first model M 0 includes the major forces and moments acting on the helicopter. The total force f B is divided into three contributions:

f B = f B u + f B l + f B g ( 9
)
where f B u and f B l are the forces generated by the upper and lower rotors, respectively, and f B g is due to the gravity. Similarly, the total moment τ B is expressed as:

τ B = τ B u + τ B l1 + τ B l2 (10)
where moments τ B u and τ B l1 are generated by the upper and lower rotor drag, respectively, and the moment τ B l2 is generated by the component of the lower rotor thrust in the (G, x B , y B ) plane. The latter depends on the swash-250 plate orientation and is null when orientation angles δ lat and δ lon are equal to zero (i.e. the lower rotor thrust is along -z B only).

The upper rotor produces a force along the z B axis, which is directly proportional to the square of the rotor 255 speed Ω u 2 , and to the upper rotor aerodynamic coefficient α u . In the same way, the lower rotor generates a force of magnitude α l Ω l 2 , but its direction depends on the swashplate angles δ lat and δ lon . Coefficients α u and α l are here considered constant, but are in fact function of the Rey-260 nolds number and the air density (see [START_REF] Koehl | Modeling 735 and identification of a launched micro air vehicle: design and experimental results[END_REF]). The drag of the upper rotor creates a moment around z B , which is equal to τ B u = γ u Ω u 2 for the upper rotor and τ B l1 = -γ l Ω l 2 for the lower one. These forces and moments are summarized in the following set of equations:

           f B u = -α u Ω u 2 [0, 0, 1] T f B l = -α l Ω l 2 [s δ lon , c δ lon s δ lat , c δ lon c δ lat ] T f B g = m g [-s θ , s φ c θ , c φ c θ ] T τ B d = (γ u Ω u 2 -γ l Ω l 2 ) [0, 0, 1] T τ B l2 = d l × f B l (11)
where

τ B d = τ B u + τ B l1
, and d l is the coordinate of the lower rotor hub w.r.t. the center of gravity, in the body frame.

In summary, the reference model M 0 consists of the common set of equations ( 1) to ( 8) together with forces and moments described in equations [START_REF] Wang | Flight dynamics modeling of coaxial rotorcraft UAVs[END_REF], with inputs

r = [φ * , θ * , Ω * , r * ] T and state vector x = [ξ I , v B , η, ω B ] T .
Usage of this model in the literature includes near-hover controller design using nonlinear techniques (in [START_REF] Gassmann | Control of a 750 gun-launched MAV for scene observation[END_REF][START_REF] Drouot | An approximate backstepping based trajectory tracking control of a gun launched micro aerial vehicle in crosswind[END_REF]) and near-hover position controller tuning (in [START_REF] Schafroth | Modeling, system identification and robust control of a coaxial micro helicopter[END_REF]). In the next paragraph, additional aerodynamic effects that are specific to the more complex models M 1 and M 2 are described.

Aerodynamic effects due to the motion in the air

Additional aerodynamic forces and moments, neglected in model M 0 , are induced by the airspeed v cp , which is the relative velocity between the air and the helicopter. The airspeed is mainly composed of the wind velocity v wind and the helicopter velocity v B . It is applied at the center of pressure (CoP), and expressed in the body frame as follows:

v cp = [u cp , v cp , w cp ] T = v B -v B wind + ω B × d cp ( 12
)
where d cp is the coordinates of the CoP in the body frame (i.e. the location at which the resultant of aerodynamic forces applies) and × indicates the cross product. In the following, v wind is considered to be a null vector (no wind during the indoor experiments).

The aerodynamic forces and moments can then be expressed as follows [START_REF] Klein | Aircraft system identification: theory and practice[END_REF][START_REF] Koehl | Modeling 735 and identification of a launched micro air vehicle: design and experimental results[END_REF][START_REF] Wang | Nonlinear modeling of a miniature fixed-pitch coaxial UAV[END_REF]:

               f B aero = - 1 2 ρ S   |u cp | u cp C x |u cp | v cp C y |u cp | w cp C z   τ B aero = - 1 2 ρ S   |u cp | p C lp |v cp | q C mq |w cp | r C nr   + d cp × F B aero ( 13
)
where parameters C x , C y , C z are aerodynamic force coefficients, and C lp , C mq , C nr are aerodynamic damping moment coefficients. These nondimensional coefficients depend on the flight conditions and on a number of quantities in a nonlinear manner, typically the Mach number, altitude, attitude and body rates. Mathematical expressions of these coefficients are given in [START_REF] Klein | Aircraft system identification: theory and practice[END_REF]. However, considering a near-hover or slow cruise flight, a fair assumption is to consider them as constant (as in [START_REF] Schafroth | Modeling, system identification and robust control of a coaxial micro helicopter[END_REF]). Parameter S is a reference surface that depends on the helicopter body, which is used to scale the aerodynamic coefficients: S = 2πR b 2 , with R b = 0.175 m. Finally, ρ denotes the air density at 295 the vehicle altitude.

These additional equations ( 12) and ( 13), together with the previous equations of model M 0 , are included in models M 1 and M 2 . Considering the symmetries of our vehicle, the x B and y B components of the coordinates d cp are 300 set to zero in both models. For model M 1 the additional hypothesis is that the CoP is located at the center of gravity (CoG), i.e. the z B component is also set to zero. Under this assumption, the cross products in equations [START_REF] Fankhauser | 745 Modeling and decoupling control of the coax micro helicopter[END_REF] and ( 13) are null. The differences between the 305 three models are summarized in Table 1.

Linearization of the model

Even if the helicopter is a nonlinear system, it is still worthwhile to have a linearized representation of the system. It is indeed of utmost importance in control design 310 where linear tools are extensive, or in order to have a better understanding of the system properties like the dynamics or the stability using linear stability criteria. Consequently, the last model M lin considered in the study comes from a jacobian linearization of model M 0 . The 315 equilibrium point is the hover flight, where the whole state vector is null (only the position in translation can be chosen arbitrarily). The model can then be written in four sets of equations where state variables are decoupled:

Roll axis dynamics.     ẏ v φ ṗ    =     0 1 0 0 0 0 g 0 0 0 0 1 0 0 0 0         ỹ ṽ φ p    +     0 -(α l Ω l 2 )/m 0 (α l d l,(z) Ω l 2 )/I xx     δlat (14)
Pitch axis dynamics.

    ẋ u θ q    =     0 1 0 0 0 0 -g 0 0 0 0 1 0 0 0 0         x ũ θ q    +     0 -(α l Ω l 2 )/m 0 -(α l d l,(z) Ω l 2 )/I yy     δlon (15)
Yaw axis dynamics.

ψ ṙ = 0 1 0 0 ψ r + 0 0 (2 γ u Ω u )/I zz -(2 γ l Ω l )/I zz Ωu Ωl (16)
Altitude dynamics. Variables Ω u and Ω l are the rotation speed of the rotors at the equilibrium point, where N = 0 (zero yaw angular rate) and F z = m g (gravity is compensated). Their expression is given as:

ż ẇ = 0 1 0 0 z w + 0 0 -(2 α u Ω u )/m -(2 α l Ω l )/m Ωu Ωl (17)
       Ω u = γ l m g a l γ u + a u γ l Ω l = γ u m g a l γ u + a u γ l (18)
The tilde notation denotes deviation variables: with an equilibrium point and input (x, u), deviation variables are given by:

x(t) = x(t) -x ũ(t) = u(t) -u (19)

Model summary 320

The four models considered in the study are summarized in Table 1. In the next section, a parameter estimation methodology is described and applied on these models. The resulting identified models are then compared in Section 5.

System identification

The models presented in the previous section depend on a number of parameters (described in Table 2) that need to be carefully estimated because they have a direct impact on the quality of the model. Some parameters are easily known or can be measured directly. However, others need specific test benches or experiments. In this study, in an attempt to lighten as much as possible the estimation procedure by getting rid of unnecessary testbenches, data recorded during a manual flight are used for both parameter estimation and evaluation of the resulting model.

In the following, the identification methodology is described and applied on the four models described in the previous section. Several steps are followed in order to guaranty the validity of the results: 1/ assessment of the observability of the parameters, 2/ choice of a suitable approach for closed-loop identification, 3/ preprocessing and analysis of exciting trajectories, and finally 4/ evaluation of the results. The latter is conducted by using an independent flight dataset to give an accurate estimate of model performance.

Determining an identifiable model

In practice, the model may not be identifiable: some parameters may have interdependencies (number of para-350 meters to be estimated is too high) and this is problematic for the good convergence of estimation algorithms, leading to incorrect estimated values for parameters. In this case, one may change the model structure or select parameters which are estimated or measured separately and fixed (i.e. 355 considered as constant during the identification process) so that the remaining ones can be distinguished. To this end, structural identifiability analyses are useful to detect collinearities between parameters [START_REF] Walter | Identification of parametric models from Experimental Data[END_REF]. An algorithm developed in [START_REF] Sedoglavic | A probabilistic algorithm to test local algebraic observability in polynomial time[END_REF] allows such an analysis using a symbolic re-360 presentation of the model and has been used in our study: an extended model is constructed, in which the parameters are considered as additional states with zero derivatives, and the structural identifiability problem is then converted into an observability problem. Through this analysis, 365 it appears that all parameters cannot be estimated together from flight data, for the models considered here. For this reason, some of them are considered constant during the estimation procedure. Fortunately, several parameters are easy to measure and are indicated with a cross in Ta-370 ble 3. The inertia matrix is estimated using the trifilar pendulum method (see [START_REF] Harris | Shock and vibration handbook[END_REF]). The aforementioned structural analysis is used to check that the remaining set of free parameters (i.e. to be estimated), shown in Table 3, is observable.

Parameter estimation methodology

The principle of the parameter estimation method is illustrated in Figure 5: given the same input u id , one searches for the parameter set that minimizes an error criterion between the model outputs and the measurements.

380

This method is known as the prediction error minimization (PEM) method and the cost function to be minimized is usually chosen as a weighted norm of the prediction

Fixed (×) / Estimated (E) Parameter M 0 M 1 M 2 M lin α l E E E E α u × × × × γ l × × × × ∆ u E E E E I xx , I yy , I zz × × × × m × × × × g × × × × K ser × × × × K mot E E E E d l,(x,y) × × × - d l,(z) × × × × d cp,(x,y) - - × - d cp,(z) - - E - C x , C y , C z - E E - C lp , C mq , C nr - E E - ρ - × × - v B wind - × × - k 1 × × × × k 2 × × × × k 3 × × × × k 4 × × × ×
Nb. of free parameters: error [START_REF] Ljung | System Identification: Theory for the User, PTR Prentice Hall Information and System Sciences Series[END_REF]. This method is well suited for non-linear models based on physical equations and it is one of the most widely used time-domain method to estimate aircraft parameters [START_REF] Jategaonkar | Flight vehicle system identification: a time domain methodology[END_REF]. Other well-known estimation methods are the sub-space methods and the instrumental variable (IV) methods [START_REF] Söderström | Instrumental variable methods 785 for system identification[END_REF], that are not considered here. The main reasons are that the PEM method is more general compared to sub-space methods, and offers simpler mathematical formulation and implementation compared to IV methods. Moreover, the latter is poorly suitable for MIMO systems.

Closed-loop identification approaches

In our case, as stated in section 2.4, the system is unstable in open loop. It is thus necessary to identify it using data from closed-loop operation. Three main approaches may be considered (described in more details in [START_REF] Garnier | A bias-eliminated leastsquares method for continuous-time model identification of closed-loop systems[END_REF]):

• The "direct" identification approach where one makes use of the input u and the output measurements y m of the helicopter system to identify it directly, ignoring the feedback and the reference input r. The main drawback of this method is that when identifying a model of a system in closed-loop from measurements at the input and output of the process, the assump-405 tion that the input and the output noise are uncorrelated is violated: the noise measurements corrupts the control signal. This results in a reduced accuracy on the parameter estimate, unless good noise models are used.

410

• The "indirect" approach where inputs r and outputs y m are used to identify the closed-loop system, from which the open-loop system is obtained by using the knowledge of the controller. In this case, any error in the controller model leads to an error in the identifi-415 cation of the closed-loop system. • Finally, the "joint input-output" method, where a new system with input r and output {u, y m } is considered, allowing identification of both the open-loop system and its controller at the same time. If the con-

420
troller K is known this method is equivalent to the "indirect" method: measurements of u is redundant because it is expressed as u = K(ry m ). In this study the controller is known and linear, so it makes more sense to take advantage of this knowledge by 425 the use of the "indirect" approach. Using this approach, identification methods established for open-loop systems can be applied here as well: the block-diagram in Figure 5 can be used directly, with u id (t) = r(t).

J = 1 N s trace(E T W E) ( 20
)
where E is the matrix of prediction errors between y and y m . It has a size of N s × N o , with N s the number of data The weighting matrix is chosen as diagonal with four gains: a weight on the position (in m -2 ), the linear velocity (in s 2 /m 2 ), the attitude (in rad -2 ) and the angular velocity (in s -2 /rad 2 ). As for the values of the weights, a natural 450 approach is to normalize the signals. They can then be fine-tuned based on the confidence on the measurements and on the model. Here, the position measurement is not used for identification (corresponding weight is set to 0), because it does not improve identifiability: according to equations (4-7), the position states are not used in the equation terms. However, the position measurement is used for validation of the identified model.

Optimization method

The "trust region" optimization method is used to adjust iteratively the values of the parameters in order to minimize the cost function [START_REF] Sedoglavic | A probabilistic algorithm to test local algebraic observability in polynomial time[END_REF]. This "Newton-step"-based method exhibits quadratical speed convergence when initialization is close to the solution [START_REF] Berghen | Levenberg-marquardt algorithms vs trust region algorithms[END_REF]. The "simplex" algorithm [START_REF] Lagarias | Convergence properties of the nelder-mead simplex method in low dimensions[END_REF], which does not require gradient computation but only function evaluations (direct search) has also been tested and leads to similar results, but was found to be slower and encountered convergence issues in our case.

Input-output flight data collection

The structural identifiability analysis (section 4.1) is a necessary condition to the success of the parameter estimation. However, it is not sufficient: the input signal must be an adequate excitation for the experimental data to be informative enough, so that the uncertainty on the parameter estimates is minimized.

Experiment design

Flights of the helicopter are conducted indoor and a motion capture system (MOCAP), provided by VICON, is used to track the helicopter movement: six near-infra-red cameras track the position of five retro-reflective markers attached to the helicopter body. Position and attitude are reconstructed and estimated with a one-millimeter precision on a ground computer. Such systems have been used, for instance, as real-time position and attitude feedback for position control of UAVs by [START_REF] Mellinger | Recent advances in quadrotor capabilities[END_REF], and for identification of a model of coaxial helicopters by [START_REF] Fankhauser | 745 Modeling and decoupling control of the coax micro helicopter[END_REF]. Details on the setup are available in [START_REF] Chellal | A kinematic vision-based position control of a 6-dof cable-driven parallel robot[END_REF].

Several trajectories have been recorded during manual flights of the helicopter equipped with the attitude controller described in Section 2.4. The flight area was limited by the field of view of the motion capture system to a cube of 3 meters of side. As much as possible in this limited area and with the static controller, the pilot was asked to do fast movements within this cube, in an attempt to excite the system on a wide bandwidth. Several trajectories have been recorded and a frequency analysis has allowed the selection of the most informative ones: good identification requires good exciting trajectories. The focus is on a bandwidth of about 13 rad/s. Two flight datasets of approximately 35 seconds are selected. One is used for parameter estimation, while the other one is used for validation of identified models.

Measurement preprocessing

For identification and evaluation of the models, measurements of the following states are needed: {ξ I , v B , η, ω B }. Position measurements are given by the MOCAP system, from which velocities are computed. These are computed first in the inertial frame by differentiation, then body-axes velocities v B m are obtained from kinematic equations (see equation ( 4)). As measurement noise degrades 510 the derivative results, the use of a noise-robust method is suggested. Here, the Lanczos least-square smoothing derivative is used: the function is approximated by a polynomial near the point of interest, and the derivative of the function is assumed to be equal to the derivative of 515 the polynomial [START_REF] Lanczos | Applied Analysis[END_REF].

For identification, measurements of attitude and body angular rates given by the onboard IMU is used (also used by the attitude controller). The MOCAP system also provides estimates of the attitude, which are used for data 520 integrity verification and to synchronize the MOCAP and IMU measurements by computing the delay that maximizes the cross-correlation between the two attitude estimates.

Finally, the flight measurements datasets are filtered 525 with a second order high-pass filter with a damping ratio of 0.707 and a cut-off frequency of 1 rad/s, in order to remove the central tendency in the data. Otherwise, the estimation algorithm would emphasize the low-frequency range, which may lead to divergences (i.e. wrong estima-530 ted parameter value). The same filter is applied to the model outputs. In this section, the results are analyzed more closely and 550 compared. To this end, a metric is first defined to allow comparisons between the model outputs and the reference data.

Estimation results

Analysis and comparison of the models

Metric for model comparison

The fit value is expressed in % and defined for output i as f it i = 100 (1 -κ i ), where κ i is the Theil's Inequality Coefficient (TIC) [START_REF] Theil | Economic forecasts and policy[END_REF]. The TIC index is a normalized

κ i = 1 Ns Ns k=1 [y i (k) -y m,i (k)]
where y i is a vector of N s sampled values of the model out-555 put i, and y m,i is the vector of measurements for output i.

The TIC index κ i is bounded between zero (perfect fit) and one (worst-case deviation). It is widely used in the aircraft model identification domain and, as a consequence, allows straightforward comparisons with the results presented in 560 the literature. In particular, it is usually considered that when comparing aircraft model outputs to flight data, a value below 0.3 (i.e. a fit above 70%) typically corresponds to an accurate model [START_REF] Dorobantu | Validating uncertain aircraft simulation models using flight test data[END_REF][START_REF] Jategaonkar | Flight vehicle system identification: a time domain methodology[END_REF].

Fit analysis

After identification of the models, the fit between model outputs and measurements is computed and summarized in Figure 10 for the Y axis and in Figure 11 for the Z axis. As the X and Y axes show similar dynamics (see equations (4-7)), the results are close to each other and only the Y axis is shown here. In Figures 10 and11, the fit values obtained with the identification set are shown in dashed lines, the ones obtained with the validation set are in continuous lines. Close fit values with identification and validation datasets means a reliable identification.

Results on the non-linear models

Whereas the four models have similar fits on angular rates (between 54.0% for the Y fit of model M lin and 59.5% for M 2 ), this is not the case for the translational movement. In particular, the position and velocity fits of In fact, taking into account the aerodynamic forces and moments generated by the induced wind (models M 1 and M 2 ) dramatically increases the accuracy on the translational axes. When comparing M 1 to the reference model 590 M 0 , the validation fit on the Y axisgoes from 33.7% to 82.7% for ξ I and from 13.1% to 81.2% for v B (Figure 10). A significant increase is also noticed for the Z fit on the translational axes. These results emphasize the major role played by aerodynamic drag in translational movements.

595

It has a damping effect clearly seen in Figure 7.

Additionally, estimating the z B component of d cp (model M 2 ) further improves the results, especially for the helicopter attitude (from 46.9% to 65.1%). The value of d cp,(z) converges to -2.2 cm, which is credible (see [START_REF] Leishman | Principles of helicopter aerodynamics[END_REF]).

600

Because of its small value, it is often considered negligible in the literature (e.g. [START_REF] Koehl | Modeling 735 and identification of a launched micro air vehicle: design and experimental results[END_REF][START_REF] Wang | Flight dynamics modeling of coaxial rotorcraft UAVs[END_REF][START_REF] Fankhauser | 745 Modeling and decoupling control of the coax micro helicopter[END_REF] consider CoP and CoG at the same point). However, as highlighted in this study, the impact on the fit is quite important for the helicopter platform under study and it adds only moderate complexity 605 to the model and only one additional parameter compared to model M 1 (but 7 additional parameters compared to M 0 ). Indeed, M 1 adds parameters C x , C y , C z , C lp , C mq and C nr , and M 2 adds d cp,(z) on top of that.

Results on the linear model

Figures 10 and 11 also show that the difference between identification and validation fit is small for almost all models and axes, except from the velocity output v B of model M 0 . This shows that the optimization algorithm converges to wrong values of parameters in this case, which is 620 most likely due to the undamped perturbations u q and v p because the effect is not seen for model M 1 . The very close results between identification and validation fits on the other axes denote a high confidence on the identification.

625

It should also be noted that the validation dataset includes large roll and pitch angles of the helicopter (up to ±30 • ), linear speeds up to ±1.5 m/s and angular speeds up to ±60 • /s. Therefore, for variations within this range of speeds and angles, accuracies given in figures 10 and 11 630 can be expected for each model.

Discussion

The identification complexity of a model increases with the number of parameters to be estimated. Consequently, in the standard modeling framework, the number of free 635 parameters is typically used as a simple but fair complexity metrics [START_REF] Spiegelhalter | 810 Bayesian measures of model complexity and fit[END_REF]. An increase in the complexity of the model naturally leads to an increase of goodness-of-fit, and the choice of a model is necessarily a trade-off.

The performance of each model with respect to its com-640 plexity can be conveniently represented in a plane, given in Figure 12, where the mean fit value over the outputs is computed. The points then outline the Pareto front, which gives the maximal accuracy that can be expected with a model of a given complexity (more information on 645 Pareto multi-objective optimization is given in [START_REF] Ngatchou | Pareto multi objective optimization[END_REF]). The trade-off between simplicity and accuracy can then be readily addressed in this form: select the simpler model that reaches the required accuracy for the considered application.

650

For example, considering an application that requires a performance index above 60% on the rotational motion, one would select M lin as a suitable model, based on figure 12. If the same level of performance is required on the translational motion as well, model M 1 would be selected. Lastly, if a performance index above 70% is required on the rotational motion, model M 2 would be of adequate complexity.

Obviously, modeling accuracy can be further improved by taking into account additional phenomena, at the price of a higher complexity. In particular, the flapping motion of the blades has an effect at high frequency, close to the rotation speed of the blades (about 35 Hz, and with harmonics above this frequency). The coupling between the fuselage (rigid-body dynamics) and the rotors (flapping dynamics) is responsible for the 4 Hz vibrations that can be seen on Figure 9. The flapping dynamics is described for example by Hürzeler in [START_REF] Hurzeler | Modeling and design of unmanned rotorcraft systems for contact based inspection[END_REF] where it is modeled as an aerodynamically damped oscillator which is excited through swashplate orientation change and by the helicopter body motion, or in more details by Padfield in [35, p. 93]. This is an inherent behavior of helicopters, and its frequency depends essentially on rotor characteristics (mainly the blade stiffness) and cannot be controlled through the helicopter control inputs. Still, if a higher bandwidth of the model is needed, one should include this phenomena: a model with second order flapping dynamics, and otherwise similar to model M 0 , is given in [START_REF] Tang | Nonlinear system identification of a small-scale unmanned helicopter[END_REF]. This however adds significant complexity and additional states for the flapping dynamics.

The ground effects are also neglected in these models. They could be easily included and would lead to improved results especially on the Z position and velocity (figure 11), again at the price of a higher complexity.

Conclusion 685

The accuracy-simplicity trade-off is a fundamental question in modeling. The best model depends on the target application, and should be as simple as possible while still capturing the relevant phenomena. In this paper, three nonlinear parametric models with increasing complexity 690 and one linearized model for near-hover flight are built and identified. They are compared in terms of complexity (number of parameters, non-linearities) and fits to measurements from a manual flight. If the four models show a good fit on attitude and angular rate, the two simplest fail 695 at reproducing the translational dynamics. Important aerodynamic phenomena are thus highlighted, like the drag from lateral induced wind and the position of the center of pressure. These comparisons allow a sensible choice among these model structures, depending on the needed accuracy 700 in reproducing angular rate, attitude, speed and position of the helicopter. The methodology that has been proposed and applied on a miniature coaxial helicopter is fully generic and can be applied to other types of helicopters and application domains. It should lead to similar conclusions for conventional helicopter MAVs; the differences in accuracy between the models may yet be lower for other types of multi-rotor helicopters (e.g. quadcopters), on which the aerodynamic effects considered here generally have less impact.

Figure 1 :

 1 Figure 1: The low-cost coaxial helicopter platform with reflective markers for tracking.

Figure 2 :

 2 Figure 2: Block diagram of the helicopter model in open-loop

Figure 3 :

 3 Figure 3: Axes definition: the body reference frame B is linked to the inertial frame I by the Euler angles -roll φ, pitch θ and yaw ψ

Figure 4 :

 4 Figure 4: Closed-loop model: controller for manual flight

  245

  325

Figure 5 :

 5 Figure 5: Basic principle of parameter estimation: given the same input u id , find the parameter set θ that minimizes an error between model outputs y and measurements y m .

  reference input r = [φ * , θ * , Ω * , r * ] T , given through the remote control, is recorded during a manual flight together with measurement y m = [ξ I m , v B m , η m , ω B m ] T , and set as input to the closed-loop model. Then, a cost function between y m and the closed-loop 435 model outputs y is computed and minimized to find the optimal parameters. A usual cost function choice is the weighted sum of squared errors:

440

  samples and N o the number of model outputs. The matrix W is a diagonal matrix of size equal to the number of outputs N o , and is used as a weight on the prediction error for each output. The weighting matrix can help improving the convergence of the optimization algorithm and 445 avoiding local minima.

Figures 6 to 9

 9 Figures 6 to 9 show time-domain plots of the model outputs after identification. The estimated values for the 535

Figures 6 to 9

 9 Figures 6 to 9 show that time-domain behavior and performance are uneven depending on the model considered. In this section, the results are analyzed more closely and

Figure 6 :

 6 Figure 6: Validation flight dataset: Position ξ I of the UAV along x I , y I and z I . Measurements are in thick black line, output of model M 0 in red dashed line, M 1 in blue dotted line, M 2 in thin magenta line, M lin in green dotted-dashed line.

580 model M 0

 0 are relatively low. Moreover, a large gap can be seen between the identification and validation fit values on the velocity. Hence, this model should not be trusted to predict translational behavior of the helicopter. It lacks important phenomena which are included in models M 1 585 and M 2 .

Figure 7 :

 7 Figure 7: Validation flight dataset: Velocity v B of the UAV along x B , y B and z B . Measurements are in thick black line, output of model M 0 in red dashed line, M 1 in blue dotted line, M 2 in thin magenta line, M lin in green dotted-dashed line.

Figure 8 :

 8 Figure 8: Validation flight dataset: Euler angles η: roll, pitch and yaw axes. Measurements are in thick black line, output of model M 0 in red dashed line, M 1 in blue dotted line, M 2 in thin magenta line, M lin in green dash-dotted line.

Figure 9 :

 9 Figure 9: Validation flight dataset: Angular rate ω B of the UAV along x B , y B and z B . Measurements are in thick black line, output of model M 0 in red dashed line, M 1 in blue dotted line, M 2 in thin magenta line, M lin in green dotted-dashed line.

Figure 10 :

 10 Figure 10: TIC-based fit value in % on Y axis, with the identification set (dashed lines) and the validation set (continuous lines)

Figure 12 :

 12 Figure 12: Accuracy (mean fit value of the outputs, in %) versus complexity for the four models considered. The red dashed line represents the performance for the rotational dynamics, the blue line represents the mean value over all outputs (rotational and translational dynamics).

Table 1 :

 1 Summary of the four models considered in the study

Table 2 :

 2 Description of model parameters: a × indicates that the parameter is present in the model I xx , I yy , I zz Inertia along x B -, y B -and z B -axis, resp.

	Symbol	Description	M 0 M 1 M 2 M lin
	α l	Thrust coefficient of lower rotor	×	×	×	×
	α u	Thrust coefficient of upper rotor	×	×	×	×
	γ l	Drag coefficient of lower rotor	×	×	×	×
	∆ u	Ratio between γ u and γ l	×	×	×	×
			×	×	×	×
	m	Mass of the UAV	×	×	×	×
	g	Gravitationnal acceleration	×	×	×	×
	K ser	Servomotor constant	×	×	×	×
	K mot	Motor and ESC constant	×	×	×	×
	d l,(x,y)	x B and y B coord. of lower rotor hub	×	×	×	-
	d l,(z)	z B coord. of lower rotor hub	×	×	×	×
	d cp	Coord. of center of pressure in body axes	-	-	×	-
	C x , C y , C z	Aerodynamic force coefficients	-	×	×	-
	C lp , C mq , C nr Aerodynamic damping moment coef.	-	×	×	-
	ρ	Air density	-	×	×	-
	v B wind k 1	External wind velocity in body axes Roll/pitch gain	-×	× ×	× ×	-×
	k 2	Roll/pitch rate gain	×	×	×	×
	k 3	Thrust gain	×	×	×	×
	k 4	Yaw rate gain	×	×	×	×

Table 3 :

 3 Parameters estimated for each model

Table 4 :

 4 Estimated parameter values

		Value of parameter for model:
	Parameter	M 0	M 1	M 2	M lin Unit
	α l	36.2	64.4	34.7	79.3 µN/s 2
	α u		36.2			µN/s 2
	γ l		1.62			µN m/s 2
	∆ u	1.067	1.134	1.076	1.000 n.d.
	I xx		1.22			g m 2
	I yy		1.23			g m 2
	I zz		0.77			g m 2
	m		325			g
	g		9.81			m/s 2
	K ser		0.43			rad
	K mot	52.5	83.0	129.4	82.8 rad/s
	d l,(z)		-76			mm
	d cp,(z)	-	-	-22	-	mm
	C x	-	4.5	1.0	-	n.d.
	C y	-	2.1	0.6	-	n.d.
	C z	-	0.8	1.0	-	n.d.
	C lp	-	1.8 • 10 -3 1.6 • 10 -2	-	n.d.
	C mq	-	1.5 • 10 -3 8.0 • 10 -3	-	n.d.
	C nr	-	1.7 • 10 -2 1.0 • 10 -2	-	n.d.
	ρ		1.226			kg/m 3
	k 1		6.5			s -1
	k 2		0.13			s
	k 3		1.0			n.d.
	k 4		0.2			s
	version of the root mean square error, computed for each
	output i, and defined by the following equation:

ing Practice 25 (2014) 1-15.

Another noteworthy result is the surprisingly good accuracy of the linear model M lin on the velocity v B , especially on the Z axis. Compared to model M 0 from which it is derived, the Y fit goes from 13.1% to 36.8%, and the Z fit goes from 14.6% to 71.3%. The results with identification and validation sets are close, which means that the identification is reliable. To analyze this result, let us write the differential equation associated to ẇ for M 1 . It can be developed from equation ( 5) and [START_REF] Wang | Flight dynamics modeling of coaxial rotorcraft UAVs[END_REF] as follows:

where the first term accounts for the rate of change of v B , and the second term accounts for the axis system rotation and contains coupled velocity terms u q and v p (see also [START_REF] Klein | Aircraft system identification: theory and practice[END_REF]). The coupling terms are close to zero if changes in orientation are small. They disappear in the linearized model M lin , and, around the equilibrium point, equation ( 22) becomes (see equation ( 17)):

Without the damping effect -K w 2 brought by the aerodynamic drag (included in models M 1 and M 2 ), the coupling terms add undamped errors and perturbations to the velocity output of M 0 and thus leads to a lower fit to measurements than its linearized version M lin .