EKF-Based Parameter Identification of Multi-Rotor Unmanned Aerial VehiclesModels

This work presents a method for estimating the model parameters of multi-rotor unmanned aerial vehicles by means of an extended Kalman filter. Different from test-bed based identification methods, the proposed approach estimates all the model parameters of a multi-rotor aerial vehicle, using a single online estimation process that integrates measurements that can be obtained directly from onboard sensors commonly available in this kind of UAV. In order to develop the proposed method, the observability property of the system is investigated by means of a nonlinear observability analysis. First, the dynamic models of three classes of multi-rotor aerial vehicles are presented. Then, in order to carry out the observability analysis, the state vector is augmented by considering the parameters to be identified as state variables with zero dynamics. From the analysis, the sets of measurements from which the model parameters can be estimated are derived. Furthermore, the necessary conditions that must be satisfied in order to obtain the observability results are given. An extensive set of computer simulations is carried out in order to validate the proposed method. According to the simulation results, it is feasible to estimate all the model parameters of a multi-rotor aerial vehicle in a single estimation process by means of an extended Kalman filter that is updated with measurements obtained directly from the onboard sensors. Furthermore, in order to better validate the proposed method, the model parameters of a custom-built quadrotor were estimated from actual flight log data. The experimental results show that the proposed method is suitable to be practically applied.


Introduction
Unmanned Aerial Vehicles (UAVs) with rotary wings (multi-rotors) allow great flexibility of movements, which makes them very useful for several tasks and applications. In this case, one of the main objectives of the research community has been the improvement of the autonomy of these systems. In this context, in order to develop autonomous control systems that perform well in terms of robustness, stability, precision, and adaptability, it is fundamental to have mathematical models that represent the actual dynamic behavior of the UAV in a precise manner.
Different methods have been proposed for deriving models of multi-rotor aerial vehicles. Some of them rely only on observed input-output data; for instance in [1], a frequency domain system identification method was used to obtain a linear representation of the dynamics of a quadrotor. Models that are relatively robust to errors associated with processing noise and bias effects can be

Objectives and Contributions
As can be seen from the previous examples, and to our knowledge, all the approaches that can be found in the literature treat the identification of inertia parameters in a separate experimental process from the rotor (thrust and drag coefficient) parameter identification process. Moreover, model parameters are typically identified in a statical manner using custom-built test beds.
In this present research work, the authors investigate the possibility of estimating all the model parameters by means of a unique close-loop estimation process, from measurements that can be obtained directly from the onboard sensors that are commonly available in any multi-rotor setup. Therefore, the objective is to develop an algorithm that is potentially suitable to be used for online identification of the model parameters, thus avoiding the need for laboratory-based identification experiments that make use of custom-built or commercially-available test beds.
The proposed approach is based on the theoretical results obtained from a nonlinear observability analysis. The objective is to find the theoretical conditions needed for estimating the model parameters from different sets of sensor measurements. If a system is observable, then all its internal states can be estimated in a finite time by measuring the system output. Three classes of multi-rotor aerial vehicles are considered in order to carry out the observability analysis. In this case, the state vector of the aerial vehicle is augmented by considering the parameters to be identified as state variables with zero dynamics. Based on this analysis, the sets of measurements and the necessary conditions that must be satisfied will be derived. Those measurements will allow the estimation of the model parameters. An Extended Kalman Filter (EKF) is used for estimating the model parameters. However, since observability properties are inherent to a system, then another state-space-based estimation technique could also be used for the same purpose. In this sense, it is important to note that EKF has been also previously used for parameter estimation of aerial vehicles. For instance, in [20], the aerodynamic parameters of a fixed-wing aircraft and a rotary-wing (helicopter) UAV were estimated. More recently, in [21], a linear Kalman filter-based method was proposed for identifying the dynamic model parameters of quadrotors using flight data. The above method is similar to our work in the sense that both approaches aim to estimate the dynamic model parameters using data available from the onboard sensors. However, in [21], the parameters to be estimated were grouped into coefficients. In our work, the model parameters are explicitly estimated (i.e., inertia parameters). Furthermore, in the previous work, a linear analysis approach was used, while in the present work, a non-linear approach is proposed. Moreover, in [21], only a subset of coefficients was recursively estimated (the least squares method was also used), while in this proposed approach, all the parameters are recursively estimated by means of EKF. Additionally, the proposed approach is developed for three class of multi-rotor vehicles instead of one.
Because the proposed identification method operates in a closed-loop manner and a multi-rotor is open-loop unstable, the identification experiments must be carried out in closed loop, by means of a human operator or automatic control. However, even in the case of custom-built multi-rotors, it is not difficult to use hand-tuned available low-cost autopilots to achieve the maneuvers needed for the identification process. While the proposed technique is suitable for on-line identification, in this work, and for the sake of simplicity, it will be assumed that the data are collected in a flight log during a manually-controlled flight to be processed off-line afterwards.

Paper Outline
The document is structured as follows: Section 2 presents the 12-state model used for representing the kinematics and dynamics of a multi-rotor vehicle. Section 3 presents the specifications and assumptions for the proposed system. Section 4 presents the nonlinear observability analysis from which the theoretical bases of the proposed method are obtained. Section 5 presents the EKF-based scheme used for estimating the model parameters. In Section 6, the results obtained from numerical simulations are shown in order to validate the proposal. Finally, in Section 7, some concluding remarks are given.

Multi-Rotor Modeling
In this section, a mathematical model of a multi-rotor UAV is presented. For the analysis presented in this work, three common classes of multi-rotor vehicles will be considered: (i) a quadrotor in + configuration, (ii) a quadrotor in × configuration, and (iii) a hexarotor. Furthermore, note that it should be straightforward to extend the proposed approach to other multi-rotor UAV configurations. Figure 1 shows the notation used for these three common classes of multi-rotor UAV.
In a multi-rotor UAV, the thrust T i = bω 2 i {i = 1, 2, 3, .., n} is an upward vector generated by each rotor driven by an electric motor. ω i is the rotor speed, and b > 0 is the lift constant that depends on the air density and the characteristics of the blade. T = ∑ T i is the total upward thrust. An aerodynamic drag Q i = kω 2 i is opposed to the torque applied to each propeller by the motor. Pairwise differences in rotor thrust cause the vehicle to rotate. For instant, for the quadrotor in × configuration described in Figure 1b, the torque about the x-axis (rolling torque) is τ x = dT 2 + dT 3 − dT 1 − dT 4 , where d is the distance from center of mass to the rotors generating the thrust.  The total thrust T and torques (τ x , τ y , τ z ) produced by varying each rotor speed ω i can be defined in a compact way using a matrix equation. For the case of the quadrotor in + configuration (see Figure 1a): For the case of the quadrotor in × configuration (see Figure 1b): For the case of the hexarotor (see Figure 1c): The following state vector will be used for representing the motion of the multi-rotor: The north-east-down position (p n , p e , p d ) is defined relative to the inertial frame W . The linear velocities (u, v, w) and the angular velocities (p, q, r) are defined relative to the body frame B . The Euler angles, roll φ, pitch θ, and heading (yaw) ψ represent the orientation of the vehicle.
A six-degree-of-freedom, 12-state model for the UAV kinematics and dynamics is given by the following equations [22]: where: where m is the mass of the UAV and g is the gravitational constant. R WB = (R WB ) T is the rotation matrix that transforms from the inertial frame W to the body frame B . J is the inertia matrix of the vehicle. For multi-rotor vehicles, a common assumption is that the mass distribution is closely symmetrical with respect to the vehicle coordinate frame. In this case, the products of inertia of J are zero, and only the moments of inertia (J x , J y , J z ) will be considered. Since d, d 1 , d 2 , d 3 , and m are parameters that will be easily measured, J x , J y , J z , b, and k are parameters whose values must be identified in order to define a dynamic model of a specific multi-rotor UAV.

System Definition
For the proposed approach, the model parameters (J x , J y , J z , b, k) will be treated as state variables, in order to estimate them from a set of measurements obtained from the sensors commonly available in a multi-rotor UAV. In this case, the state vector of the vehicle x v Equation (4) will be augmented as follows, in order to define the system state x: For the parameters, zero dynamics will be assumed: In practical scenarios, UAVs are commonly equipped with a set of sensors that at least include: an Inertial Measurements Unit (IMU) and a Global Positioning System (GPS). At the same time, IMUs typically include accelerometers, gyroscopes, magnetometers, and barometer. The navigation systems make use of the onboard sensors in order to estimate the state of the vehicle. In this case, the most common approaches for addressing the problem of the UAVs' navigation (position estimation) is the Inertial Navigation Systems (INS), which typically fuses the GPS measurements and inertial measurements from the IMU. Furthermore, the available Attitude and Heading Reference System (AHRS), which commonly relies on the IMU, can be used to provide a reliable estimation of the attitude of the vehicle. Moreover, it is also common that navigation systems now include visual-based sensors like optical flow sensors or cameras in order to improve position estimation. Figure 2 shows a simplified block diagram illustrating a high-level abstraction of the proposed EKF-based model parameter identification approach. At this stage, for the investigation purposes, it will be assumed that the following measurement vector z can be obtained from the navigation system or directly from the sensors mounted on the multi-rotor: z = p n p e p d u v w φ θ ψ p q r T (12) Note that the former assumption is equivalent to saying that the whole state of the vehicle Equation (4) is available. However, later, it will be shown that only subsets of z are required in order to estimate the model parameters (J x , J y , J z , b, k). It is also important to note that there is no theoretical restriction about the manner that measurements are obtained. For instance, in order to measure z, the angular velocities (p, q, r) can be obtained directly from the IMU; the attitude (φ, θ, ψ) can be obtained from the AHRS; and the position (p n , p e , p d ) and linear velocities (u, v, w) can be obtained from the INS. Of course, other configurations may be possible while the required state measurements are available. Even, an external motion capture system (e.g., VICON or OptiTrack) could be used for the same purpose. However, this work is mainly intended to show that the model parameters can be estimated by using only the onboard sensors commonly available in multi-rotor UAVs.
Additionally, it will be assumed that the control input u, which commands the rotors speed ω i , is available from the control system. In many practical scenarios, the vector u will be composed by the values of PWM (Pulse-Width Modulation) signals that are used as reference inputs of the ESC (Electronic Speed Control) units that control the rotational speed of the rotors. Obviously, if tachometers are available, or if the relation between the PWM signal and the rotor speed is known, then the vector u can be composed by the rotor speed (i.e., u = [ω 1 , ω 2 , .., ω n ]).

Observability Analysis
The observability properties of the nonlinear system defined by Equations (5)- (8) and (11) will be studied by means of a nonlinear observability analysis. The observability is a measure about how well internal states of a system can be inferred from knowledge of their external outputs (measurements). The objective is to find the theoretical conditions needed for estimating the state vector x Equation (10) and in particular the state variables (J x , J y , J z , b, k) from the measured vector z Equation (12) or a subset of it.
A system is defined as observable if the initial state x 0 , at any initial time t 0 , can be determined given the state transition model of the systemẋ = f(x), the measurement model of the system y = h(x), and the observations z[t 0 , t], from time t 0 to a finite time t. In [23], it was proven that a non-linear system is locally weakly observable if the observability rank condition rank( where L n f h m is the nth-order Lie derivative [24], of the mth scalar field h m with respect to the vector field f. Independent of the sensory source from which each component of the measurement vector z is obtained (see Section 3), it will be assumed that the vector z can be predicted directly from the state vector x. In this case, the measurement models h i (x), where i is the index for the predicted component of z, will simply be defined as follows: There are no rules regarding the order of Lie derivatives used for constructing the observability matrix O O O. In this work, only {0, 1, 2}-order Lie derivatives are tested for each measurement model h i (x).
The observability matrices were computed numerically by means of the symbolic MATLAB R toolbox. Table 1 shows the sets of measurements from which the full rank matrix condition That means that if the state vector x v is available to be measured, then the parameters (J x , J y , J z , b, k) can be indirectly estimated when treated as state variables (see Configuration (a)). Moreover, this means that the state vector x Equation (10), and thus the model parameters, can be observable even if the linear velocities (u, v, w) are not measured (see Configuration (b)). Table 2 shows the sets of measurements from which the full rank matrix condition (rank(O O O) = dim(x)) is not accomplished, but from which the unobservable modes are related only to the vehicle state variables of x v , meaning that the model parameters (J x , J y , J z , b, k) are still observable. That means that for instance, and at least theoretically, the model parameters can be estimated using only inertial measurements (see Configuration (e)). Moreover, the minimal set of sensors needed for estimating the model parameters appears to be the Configuration (f). In this case, only angular measurements and measurement of the vertical position p d are considered.
Observing Tables 1 and 2, it can be seen that the measurement of angular velocities (p, q, r) represents a necessary condition for the observability of the model parameters. Furthermore, in order to obtain the above observability results, a necessary condition that must be satisfied is p = 0 and q = 0. That means that the aerial vehicle must be varying its attitude during the data capture. The same observability results were found for the three multi-rotor configurations shown in Figure 1.

Kalman Filtering
In Section 4, the observability properties of the multi-rotor UAV model were investigated. From that analysis, some conditions were found that theoretically should allow estimating the multi-rotor model parameters from data obtained with the on-board sensors. In this work, an Extended Kalman Filter (EKF) is used as the estimation technique for this purpose. On the other hand, it is important to note that the observability properties are inherent to a system, and therefore are independent of the estimation technique. Therefore, the theoretical results presented in this work should also be valid for implementing another state-space-based estimation techniques (e.g., particle filtering, unscented Kalman filtering, etc.) for identifying the model parameters.
In order to apply the EKF, and since the sensor data were obtained in a discrete manner, a discrete-stochastic system model is characterized by the continuous dynamicsẋ = f(x, u) defined by Equations (5)-(8) and (11). In this work, the Euler method is used for this purpose, and therefore: The system measurement prediction model is defined by: Let n k ∼ N N N (0, Q k ) and r k ∼ N N N (0, R k ) be the noise vectors that affect the state and the measurements, which are assumed to be mutually uncorrelated. Let ∆t be the differential of time and k the sample step.
The prediction equations of the EKF are: The update equations of the EKF are: with: and: Let z k be the vector of actual measurements at step k, P be the covariance matrix of the system state, and K be the Kalman gain.
Let A k be the Jacobian defined by the derivatives ∂f ∂x of the prediction model f with respect to the system state x. Let C k be the Jacobian defined by the derivatives ∂h ∂x of the prediction model h with respect to the system state x. For the sake of simplicity, the expressions for derivatives are not included.
The update Equations (19)- (21) are used according to the available set of measurements. For instance, Figure 3 shows a diagram illustrating the data flow for the EKF in/Configuration (b). In this case, the linear velocities (u, v, w) are not used for updating the system statex (see Section 4). It is well worth noting that, as indicated in Figure 3, the filter was updated in a separate manner just after each type of measurement is available.

Results
In this section, experimental results are presented in order to show the performance of the proposed method. In this case, computer simulations were carried out to validate in a numerical way the proposed approach and the theoretical findings. Furthermore, experiments with real data obtained from a quadrotor were performed to have a better insight into the performance of the proposed method.

Simulations
In order to test the proposed EKF-based model parameter identification approach, under all the observability configurations defined in Section 4, a multi-rotor system as described in Figure 2 was simulated using Simulink R from MATLAB R . In this case, since all observability results were the same for the three multi-rotor configurations investigated, only one configuration was tested: a quadrotor in × configuration (Figure 1b). Table 3 shows the actual parameter values of the simulated quadrotor. % To implement the control system, a hand-tuned altitude and orientation PI control scheme was used. In this case, the control reference signals were: (p d r , φ r , θ r , ψ r ). Note that by controlling only these four variables, the remaining states can also be excited. The control vector u was obtained from the control system.
In order to obtain the measurement vector z Equation (12), the following sensors scheme was used. Angular velocities (p, q, r) were obtained by emulating the noise specifications of an Invensense MPU-6000 inertial measurement unit. The attitude (φ, θ, ψ) of the multi-rotor was provided by an AHRS, as described in [25], which in turn used the same IMU measurements. Position measurements (p n , p e , p d ) were obtained from smoothing GPS measurements with a linear Kalman filter. For the altitude measurement (p d ), the filter also fused measurements obtained from an altitude sensor configured to emulate the noise specifications of the Adafruit BMP183 barometric low-cost sensor. To model the transient behavior of the GPS error, the approach of [26] was followed. This Kalman filter also included, as state variables, the derivatives of the position variables, in order to provide the linear velocities (u, v, w) of the vehicle. Figure 4 shows a comparison between actual states of the quadrotor and the measured states used for estimating its model parameters. In this case, note that for instance the position and linear velocities measurements were not very precise. Table 4 shows the diagonal values of the matrices P ini (initial values), Q, and R used in simulations.  configurations (a-f), from whose parameters (J x , J y , J z , b, k) were observable (see Section 4), were tested. Figure 4 shows the control reference signals (p d r , φ r , θ r , ψ r ) used for exciting the quadrotor state (green signals). Note that the quadrotor was commanded to maintain a stable altitude p d , while its orientation (φ, θ, ψ) was excited by sinusoidal signals. Furthermore, the actual values for all the vehicle states in x v are shown (blue signals). The measured vehicle states, used for estimating the model parameters, are illustrated by the red signals. Figure 5 shows the evolution over time of the estimated model parameter values obtained from each sensor configuration (a-f). For instance, in this case, the first row of plots (Configuration (a)) shows the results obtained when the measurements of all the vehicle states are used for updating the EKF. In the same manner, the last row of plots (configuration (f)) shows the results obtained when only measurements of altitude and angular velocities were used for updating the filter. Observing Figure 5, it can be seen that, with all the sensors configurations, it was possible to estimate in a good manner the model parameters values (J x , J y , J z , b, k). Table 5 shows the average values of the parameters estimated and their respective mean errors obtained for each sensor configuration. The above results validated the theoretical findings obtained from the observability analysis regarding the measurements needed for estimating the parameters.

Observability Conditions Test
From the observability results (Section 4), it was found that a necessary condition, which must be satisfied in order to make observable the full set of model parameters, had to fulfill that p = 0 and q = 0. This result had important implications for the success of the method because it implied that the aerial vehicle must be varying its attitude in order to estimate the parameters' values. In order to validate the observability necessary conditions, two tests were carried out. For these tests, only the results obtained from the sensors in Configuration (d) are shown, but similar results were obtained with the other configurations. The following cases were tested: (i) Case 1. For the first case, the quadrotor was commanded to ascend and maintain a stable altitude (hovering position). Therefore, from time t = 15 s to the end of the simulation, the yaw angle ψ of the quadrotor was varied.
(ii) Case 2. For the second case, the quadrotor was commanded to ascend and maintain a stable altitude (hovering position). Then, from time t = 10 s to the end of the simulation, the roll angle φ of the quadrotor was varied. From time t = 25 s to the end of the simulation, the pitch angle θ of the quadrotor was varied. Finally, from time t = 40 s to the end of the simulation, the yaw angle ψ of the quadrotor was varied. Figure 6 shows the results obtained Case 1 (left column plots) and Case 2 (right column plots). The upper two rows of plots show the movements (position and attitude) of the quadrotor for each case. The lower plots show the evolution over time of the values estimated for each model parameter. Analyzing the plots, the following remarks can be stated. The estimation of parameter b did not depend on the movement of the multi-rotor vehicle; even if the aerial vehicle was maintained in a hover position. Observing Case 1, it was clear that only moving the multi-rotor vehicle over its z-axis was not enough to estimate the parameters k and J z . In this case, a variation in the yaw ψ excited the states of these parameters (k and J z ), but their estimates clearly diverged from the actual values. Observing Case 2, it can be observed that just moving the multi-rotor over its x-axis (roll ψ), the parameter J x became observable, but not the parameters J z , J y , and k. It is not illustrated in the plots, but the same effect occurred if the multi-rotor was rotated only over its y-axis (in addition to b, only the parameter J y was observable). On the other hand, the fact of accomplishing the sufficient conditions (p = 0 and q = 0) made all the parameters observable. In this case, it is important to note that varying roll and pitch angles produced a small perturbation over the yaw of the vehicle (see t the plots from time t = 25 s to t = 40 s), thus exciting the states of the related z-axis parameters (i.e., k and J z ). However, in the above case, the convergence was very slow. Note that, when the yaw ψ of the multi-rotor vehicle varied abruptly (which happened in t = 40 s), the convergence of parameters k and J z became faster. The results obtained from these simulation tests validated the theoretical findings derived from the observability analysis.

Experimental Case
In order to validate experimentally the proposed method with real data, the model parameters of a custom-built quadrotor in × configuration were estimated (see Figure 7). The quadrotor was equipped with a standard Pixhawk-PX4 autopilot unit. The on-board sensors included in the Pixhawk were: Accel/Gyro ICM-20689, Accel/Gyro BMI055, Magnetometer IST8310, and Barometer MS5611. For performing the experiments, the quadrotor was flown manually by means of a radio transmitter. In order to accomplish the observability conditions, during the flights, the quadrotor was commanded to vary as much as possible its attitude. During each flight, the Pixhawk-PX4 created a log file in which the sensors' data and other several data, like the attitude, position, and velocities estimates, were stored. Therefore, for estimating the model parameters of the quadrotor, the flight log data were used as input to a MATLAB R implementation of the proposed method. In this case, the following data were used: (i) raw angular velocity from gyros {p, q, r}, (ii) estimated attitude {φ, θ, ψ}, (iii) estimated altitude {p d }, and (iv) PWM signals used for controlling the rotors speed, (see Figure 8).   For comparison purposes, the thrust coefficient b and the drag coefficient k were measured using a custom-built test bench in a similar manner as in [14]. Figure 10 shows the measured data and the linear approximation that best fit the data. According to the linear approximation, the following values were measured: b = 0.0064 (N/PWM) and k = 1.766 × 10 4 (N·m/PWM). Inertia matrix parameters (J x ,J y ,J z ) were experimentally measured using a test bench (see [8]). Therefore, in order to have at least an insight about the values of these parameters, a geometric approximation method as that proposed in [27] was used. Observing Figure 9 (lower plot), a small offset between the estimated values and the values obtained from the geometric approximation can be seen.   Figure 10. Linear approximation models obtained for the thrust coefficient b and the drag coefficient k using a custom-built test bench.

Conclusions
In this work, a method for identifying the model parameters of a multi-rotor vehicle was presented, using a single online EKF-based estimation process that integrated measurements that can be obtained directly from onboard sensors commonly available in this class of UAVs. The objective was to provide a practical and reliable estimate of the model parameters, but without the need for static test-bed-based identification experiments.
In order to investigate the theoretical conditions that were necessary for estimating the model parameters from different sets of sensor measurements, a nonlinear observability analysis was carried out. As an outcome of this analysis, several sensor configurations along with sufficient conditions of observability were obtained. In this work, an extended Kalman filter was used as the estimation technique for identifying the model parameters. However, in future works, it could be of interest to test other state-space-based estimation techniques, since observability properties are inherent to a system, and they do not depend on a particular estimation technique. Since the proposed method worked in a closed-loop manner, the requirement of operation was by means of a human operator or automatic control of the multi-rotor vehicle that carried out the flight maneuvers needed for accomplishing the observability conditions of the system. In order to validate the proposed method, an extensive set of computer simulations was carried out. The results obtained from simulation supported the theoretical findings. Hence, it was feasible to estimate all the model parameters in a single estimation process from measurements that could be obtained directly from the onboard sensors of the multi-rotor. Moreover, the model parameters of a custom-built quadrotor were estimated from actual flight log data. The experimental results showed that the proposed method could also be applied in a practical manner.