Experimental Parameter Identifications of a Quadrotor by Using an Optimized Trajectory

In this document, the parameter identification of a quadrotor is discussed. More precisely, the aim of this paper is to present results on the application of known methods for estimating the dynamic parameters that capture better the behavior of a quadrotor in comparison with the nominal parameters given by the manufacturer. To take into account the limitations of position, velocity, and acceleration of the quadrotor, an optimized trajectory to excite the quadrotor dynamics adequately is obtained. A proportional-integral-derivative (PID) control scheme is used to implement experimentally the tracking of the optimized trajectory. The obtained data is processed off-line to construct the standard and filtered regression models from which the parameter identification is achieved. Specifically, the least-squares and gradient descent algorithms are applied to the regression models giving four sets of estimated parameters. The four sets of parameters obtained in this work are compared with the parameters provided by the manufacturer by computing the error between simulations and experiments. In addition, the output prediction errors of the regression models are computed, thus providing another validation form. All the comparisons show that the estimated parameters are more precise than the nominal ones. The given results support the functionality of the described methodology.


I. INTRODUCTION
During the last decade, the interest of the scientific and industrial community has been focused on unmanned aerial vehicles (UAVs). These vehicles can be operated by remote control or autonomously. Those propelled by four rotors are called quadrotors. These flying robots are used in many application fields due to their simple structure, small size, high maneuverability, and hovering capability. The main applications of quadrotors are surveillance, deployment, and exploration tasks. In particular, they are used in harvesting and wildlife monitoring, search and rescue, support and reconnaissance on high-risk zones, wind turbine blades and solar panel inspection, supply delivery, and natural disaster zone mapping [1]- [4].
Quadrotors typically consists of a symmetrical lightweight airframe, four brushless DC motors mounted upon it, four The associate editor coordinating the review of this manuscript and approving it for publication was Huaqing Li . fixed-pitch propellers attached to the motors, drivers for each motor, an inertial measurement unit (IMU), a global positioning system (GPS) unit, high-density LiPo batteries and a flight computer loaded with control algorithms [3], [5], [6].
Many works addressing posture regulation and trajectory tracking of quadrotors by using different control schemes and techniques were developed in recent years. The PID control scheme has been studied and implemented to stabilize different kinds of physical systems; it is used in the industry due to its simple structure. The obtained results from different works support the functionality of this control scheme in quadrotors, even when the parameters of the vehicle are inaccurate or unknown.
The PID structure was used in [7] and [8] to control a quadrotor. An intelligent controller scheme was added to the PID scheme in [9]. Also, many robust control schemes have been used to achieve pose regulation and trajectory tracking for quadrotors. Control schemes based on the H ∞ philosophy [10], the feedback linearization methodology [11], VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ sliding mode control [12], [13], and model predictive control [14] can be found in the literature. Despite the existence of robust control schemes, the knowledge of the dynamic model is essential to develop model-based controllers [15]- [20] and to predict the system behavior using numerical simulations. A shortcoming is that the exact knowledge of the parameters of the dynamic model of the quadrotor is not always available. Besides, many schemes make assumptions about the structure and conditions of use of the quadrotors, which give rise to simplified dynamic models. For example, the vehicle weight is uniformly distributed and has symmetric geometry, the movement in the horizontal plane is at low speed, and the roll and pitch angles are small. However, they are not suitable at high-velocity flight, aggressive maneuver, and close to ground flight [3]. Modeling and control techniques neglecting aerodynamic effects at high-velocity are inadequate for quadrotor trajectory tracking tasks [21]. System identification is a useful alternative to precisely obtain some parameters that are difficult to estimate. This procedure is carried out processing data taken from real-time experiments. The effectiveness of the parameter identification for robot manipulators was validated experimentally in multiple works [22]- [25]; those procedures can be extended to other complex systems, such as quadrotors.
A review of the parameter identification of quadrotors is given next. A time-domain system identification software was implemented in [26] to estimate a linearized model of a quadrotor taking the quadrotor flight test data. A parameter identification technique using the state estimation method employing the unscented Kalman filter was described in [27]. A time-domain identification procedure for a pitch, roll, and yaw subsystem of a quadrotor was presented in [28]. Parameters of a quadrotor were obtained in [29] by using the CAD model of the airframe and three different test rigs. The parameter identification results were validated experimentally. The dynamic model of a Parrot AR.Drone obtained by means of a parameter identification process was presented in [30]. The discussion of the obtained results by implementing a continuous-time predictor-based subspace identification approach was presented in [31]. In [32], the least-squares algorithm was used to estimate the parameters of the transfer functions used to represent part of the dynamic model of a quadrotor. In [33], a nonlinear closed-loop multi-variable extremum seeking parameter identification algorithm was proposed to estimate the parameters of a quadrotor. The results of the parameter identification algorithm support the identification method. In [34], a model identification process and a robust controller design for a commercial quadrotor were presented. Comparisons between a PID scheme and an internal model control were given. Parameter identification for a continuous-time black-box model of a quadrotor was achieved in [35].
Real-time experiments are carried out either on laboratorybuilt or commercial platforms aiming to validate the results of several research works. Concerning quadrotors, a known platform in the research area is the Quanser QBall 2. This platform has been used in recent works [36]- [41] due to its programming simplicity, which is developed in MATLAB-Simulink. Another valuable point at its favor is the integration with the motion capture system Optitrack to sense its position. The manufacturer has provided a set of parameters being possible for the user to achieve a simulation phase before implementing in real-time any controller. However, providing parameters that represent the dynamics of the QBall 2 quadrotor more accurately would be attractive for the users.
While literature indicates that much effort has been put into identifying of manipulators and other mechanisms, the theoretical and experimental research for the parameter identification of quadrotors is relatively meager, which suggests an opportunity field. Moreover, many works consider a linearized or simplified model of the quadrotor and assume a symmetrical configuration with uniformly mass distribution. In addition, the experimental tests used for the parameter identification procedures reported in the literature review have not been specially optimized for this purpose. An originality point in our research is that the construction of the regression models considers the inertia tensor as a symmetric matrix with six elements instead of a diagonal matrix, thus representing better the dynamics of quadrotor QBall 2.
Another part of our contribution is the application of optimized trajectories in parameter identification of quadrotors, which has not been reported to the best of our knowledge. It is here that this work presents novel results. This methodology allows taking into account the physical limitations of mechanical systems such as the workspace or velocity limits. Besides, optimized trajectories allow exciting the dynamics of the system so that the identification process provides accurate results. It is noteworthy to mention that optimized trajectories minimize the sensitivity of the system to noise and reduce the variance of the estimated parameters [42]. Examples of works where optimized trajectories have been used for identification are in [43]- [47]. The optimization problem is generally solved using an algorithm developed to minimize the value of a cost function. Recent works using optimization algorithms to address the problem of the cost reduction of electricity dispatch in a smart grid [48], privacy masking on time-varying unbalanced directed networks [49], the collision avoidance [50], and the trajectory optimization in order to reduce flight costs and pollution for commercial aircraft [51] have been found in the literature review.
In this paper, a parameter identification procedure inspired by the methodology used for robot manipulators and industrial robots is applied to the identification of a Qball 2 quadrotor. It is worth mentioning that all the procedures which inspired this research are well-studied and have demonstrated a high degree of accuracy on parameter estimation results. The procedure includes the computation of an optimized trajectory, its experimental implementation with a PID scheme, off-line processing of the obtained data, the construction of a regression model, and the implementation of an identification scheme. The procedure is tested for several variants by using either the standard regression model or the filtered regression model as corresponds. Another tested combination is to use either the least-squares scheme or the gradient descent strategy to estimate the parameters. Thus, four sets of parameters are obtained for the quadrotor QBall 2. It is worthwhile to notice that the given procedure is general and can be applied to quadrotors having access to the same signals used in the method. As mentioned earlier, this experimental quadrotor has been used in recent research [36]- [41], and a new set of parameters with better accuracy may improve the performance of the reported model-based control schemes.
As suggested in [6], [23], [24], [42], [44], [45], [47], [52]- [56], in this work, the nominal parameters provided by the manufacturer and the new estimated parameters are assessed by comparing real-time experiments with simulations results, and by comparing the output prediction error of the regression model, which is also often called identifier error. Our results indicate that the new estimated parameters show more similarity to the experiments and lower output prediction error than the nominal parameter given by the manufacturer.
The rest of this paper is organized as follows. In Section II, an overview of the proposed parameter identification procedure is presented. The application of the procedure to a quadrotor is described in Section III. In Section IV, the discussion of the experimental validation is presented. Finally, the conclusions of this work are given in Section V.

II. QUADROTOR PARAMETER IDENTIFICATION
Different parameter identification methods applied to mechatronic systems were presented in the literature review [22]- [25], [43]- [47], [57]. Some of these methods can be adapted and modified to achieve the parameter identification of quadrotors. This Section presents the parameter identification methodology applied to quadrotors by using an optimized trajectory. The steps that define the applied procedure are described as follows.

1) Mathematical model parameterization.
• Model: the linear parameterization property should be satisfied by the quadrotor.
• Optimized trajectory design: to improve the parameter identification results, an optimized trajectory must be obtained.
• Controller implementation: a controller must be implemented to perform the trajectory tracking of the optimized trajectory.
• Since the needed signals to construct the regression models can not be obtained directly by the available hardware, off-line data postprocessing is performed in order to condition the available data and to compute such signals.
• Filter: the sampled signals must be filtered owing to the present noise coming from the data acquisition.
• Numerical differentiation: in many mechatronic systems, there are no sensors to measure certain signals which is motivated by the cost reduction. Therefore, some numerical methods should be implemented to estimate the required signals.
• Signal transformation: the kinematic relation from the inertial reference frame to the body reference frame is used to obtain the needed signals for the parameter identification procedure.

4) Parameter estimation.
• Estimation algorithm: a linear regression method is applied to estimate the system parameters by using the regression matrix constructed with the acquired data and the system output vector.
A graphic representation of the applied procedure is presented in Figure 1. It is worth mentioning that the methodology described in this work and depicted in Figure 1 is based on what is discussed in [22]- [25], [42], [44]- [47], [52]- [63].

A. QUADROTOR DYNAMIC MODEL
The dynamic model of a rigid body specifies the relationship between its motion and the applied force and torque. In the case of quadrotors, the dynamics are strongly coupled. Four actuators define the behavior of the position and orientation dynamics of the quadrotor, and a minimal change in just one actuator produces significant changes in the whole dynamics [5]. The control inputs are the total thrust and the torques produced on each rotation axis. The dynamic model of the quadrotor can be obtained under the assumption that it is a symmetric rigid body moving in a 3D space and its center of mass coincides with the geometrical center of the body reference frame. This model is represented with respect to the body reference frame as [5], [64], where m is the quadrotor mass, g is the gravity acceleration constant, v = [u v w] T ∈ R 3 and ω = [p q r] T ∈ R 3 represent the linear and angular velocity vectors in the body frame, respectively. The matrix I ∈ R 3×3 denotes the inertia tensor, e z = [0 0 1] T is the unitary vector along the z-axis in the inertial reference frame, p = [x y z] T ∈ R 3 and η = [φ θ ψ] T ∈ R 3 represent the position and orientation vectors in the inertial frame, respectively. Considering the aeronautical convention ''ZYX'', the vehicle orientation is given by the orthogonal rotation matrix R(η) ∈ R 3×3 expressed as with s x , c x and t x defining the functions sin(x), cos(x) and tan(x), respectively. Besides, W (η) ∈ R 3×3 denotes the transformation matrix F z ∈ R is the total thrust provided by the four actuators along the z axis in the body frame, τ = [τ φ τ θ τ ψ ] T ∈ R 3 denotes the torque vector in the body frame, and S(ω) is the skew symmetric matrix expressed as

1) STANDARD REGRESSION MODEL
The dynamic parameters in the quadrotor model (1) and (2) include the mass m and the inertia tensor I . Specifically, I is a symmetric matrix given by Then, the dynamic model of the quadrotor given in (1) and (2) can be rewritten as [6] u + qw−rv + g bx = 0, and I xxṗ + I xy (q − pr) + I xz (ṙ + pq) with g bx = −g sin (θ), g by = g cos (θ) sin (φ), and g bz = g cos (θ) cos (φ). The system described in (7)- (10) can be rewritten in the form of the linear regression model where (v, ω, η,v,ω) ∈ R 4×7 denotes the regression matrix, ∈ R 7 is the constant parameter vector, and F ∈ R 4 represents the system output vector. Explicitly, the components of the regression model (11) are given as follows: = m I xx I xy I xz I yy I yz I zz T , Notice that in equations (5) and (6), the parameters ∈ R 7 are not present. For simplicity, equation (11) does not take into account the expressions in (5) and (6).

2) FILTERED REGRESSION MODEL
A filtered regression model may be obtained in order to avoid either measurement or off-line calculation of the linear and angular accelerations [25], [54], [58]- [60]. The system (7)-(10) is rewritten in the following manner where Notice that a (v, ω) ∈ R 4×7 and b (v, ω, η) ∈ R 4×7 do not take into account equations (5) and (6), similarly to the computation of in the linear regression model (11). Therefore, the specific expression of a (v, ω) and b (v, ω, η) are given by A low-pass filter can be defined as where λ represents the cut-off frequency of the filter, and s is the Laplace operator. Therefore, the filtered model is obtained by multiplying equation (13) by the low-pass filter (14), leading to Thus, the filtered regression model is given by with being Details on the discrete computation of F and F F will be given later.
The identification process presented in this work is carried out using the data obtained from the experiments of the tracking optimized trajectories. These signals are selected as finite Fourier Series, similarly to the works [24], [42], [43], [45]- [47], [52], [53]. The desired trajectory and its time derivatives are given by where q di represents the i-th element of the desired position vector q di andq di denote the first and second-order time derivative of q di , respectively, w f = 2π/T is the fundamental frequency, with T being the periodic cycle time, N is the number of harmonics, and c i is the offset of the desired position q di . Each trajectory contains 2N +1 parameters. The parameters a i,l and b i,l determine the amplitude of the sinusoidal functions, and together with c i , can be obtained with an optimization method or randomly.

2) TRAJECTORY OPTIMIZATION
The problem with choosing the parameters a i,l , b i,l , and c i randomly is that the user will not know if the system to identify is able to perform such a trajectory. For this reason, an optimized trajectory is commonly used.
As was discussed in [47], the optimization problem can be established as where J is the cost function. For the case of quadrotor desired trajectories, the cost function is subject to the following constraints, which include initial and final conditions, where x max , y max , and z max are the maximum bound for the position, ψ max is the maximum bound for the angle ψ(t), z min is the minimum bound for the position in the z(t) coordinate, q di stands for any of the elements of the vector q d given in (19),q i max andq i max are the bounds for the velocity and acceleration, respectively, q 0i denotes the initial conditions for position and yaw angle of the quadrotor. According to [62], zero, the identification accuracy and the trajectory tracking may be imprecise.
There are different proposals to optimized the desired trajectory with the cost function J in (21), as can be seen in [24], [42]- [46], [55], [62]. In this work, J = cond(W ), which is the condition number of the matrix W. In particular, W is an array of the samples of the regression matrix (v, ω, η,v,ω) in (12) and is defined in (20), as shown at the bottom of the page, where T is the sampling period and is equal to 0.002 [s], k is the total number of samples, v d , ω d , η d ,v d , andω d are obtained from the kinematic relations in equations (3) and (4) as which are motivated from de position dynamics in the inertial frame, see for example [65].

3) TRAJECTORY TRACKING PID CONTROLLER
Due to the high instability of the system, an open-loop experiment can not be conducted. Therefore, the optimized trajectories q di (t) in (18) should be implemented with a control scheme. The implemented controller is a position and orientation PID control for trajectory tracking defined as where f p = [f x f y f z ] T is the force vector in the inertial reference frame given by K pp ∈ R 3×3 , K ip ∈ R 3×3 , and K dp ∈ R 3×3 are diagonal positive definite matrices for position control actions; K po ∈ R 3×3 , K io ∈ R 3×3 , and K do ∈ R 3×3 are diagonal positive definite matrices for orientation control actions;p = p d − p is the position error, p d is the desired position,η = η c d − η is the orientation error, with (29) and (30) as [65]. Notice that under the assumption that the position errorp(t) is null for all the time t ≥ 0, then θ c d = θ d and φ c d = φ d .

C. DATA PROCESSING
The next step in the identification procedure is to process the obtained signals from the real-time experiment; thus, the regression matrices and F are computed. The regression matrix in (12) is a function of v, ω, η,v, andω, while the regression matrix F in (17) depends only on v, ω and η.
Considering that only measurements of position p(t), orientation η(t) and the angular velocity ω(t) are available, the following process is carried out to obtain the necessary signals.
Notice that due to the controller limitations, a perfect tracking of the optimized trajectory is not achieved, which causes an error between the desired trajectory and the one obtained experimentally. Therefore, only the measured signals are used in the parameter identification procedure, see [24]. As seen later, the process involves the implementation of some filters.

1) Filtering.
A digital low pass filter f p (z) was designed to eliminate the high-frequency noise resulting from the sampling and quantization. Zero-phase forward-backward filtering is implemented in order to avoid distortion of position, orientation, angular velocity, and system output samples [22], [23], [25], [44], [47], [61], [63]. This filtering scheme is easy to implement in MATLAB by using the filtfilt function with a low-pass Butterworth filter with a cutoff frequency of 0.001 [Hz] using a 4-term symmetric Blackman-Harris window with L = 31 that indicates the number of samples contained in the window [25], [61], [63], [66]. This filtering procedure was applied to the position p(hT ), orientation η(hT ), angular velocity ω(hT ), and system output F(hT ) signals, where h = 0, 1, 2, . . . , k − 1, being k the total sample number and T = 0.002 [s] the sampling period. 2) Numerical differentiation. In order to compute the velocities and Euler angles change rate, the central difference algorithm was used to avoid phase shifting [22], [23], [63], [67]. The central difference algorithm is . . .
given byȧ where a f (hT ) ∈ R k is the filtered vector obtained as the output of the digital low pass filter f p (z). More specifically, if a f is either p f (hT ), v f (hT ) or ω f (hT ), then either the signalsṗ f (hT ),v f (hT ) orω f (hT ), respectively, are obtained with (33). 3) Signal transformation. The parameter identification is obtained by using the parameterized model in (11) and the filtered regression model in (16), both are represented in the body reference frame. Due to this, the numerical derivativeṗ f (hT ) obtained by using (33) is transformed into the body reference frame by using the kinematic relation in (3). The result is the linear velocity v f (hT ) in the body reference frame. Subsequently, the derivativesv f (hT ) andω f (hT ) needed in the regression matrix of the parameterized model in (11) are obtained by using the central difference algorithm in (33).

D. PARAMETER ESTIMATION
The parameter estimation procedure presented in this document considers two linear regression models. The first is the so-called standard regression model represented by equation (11) and the second one is the filtered regression model represented by equation (16). All the elements of their respective regression matrix are obtained following the steps described in Section II-C. More specifically, is presented in equation (12) for the standard regression model, and F is defined in (17) for the filtered regression model.

1) STANDARD REGRESSION MODEL CONSTRUCTION
Once the signals F f , v f , ω f , η f ,v f , andω f are computed, the following standard regression model is defined where is the regression matrix evaluated along the off-line processed signals, and ∈ R 7 is the parameter vector to be estimated.

2) FILTERED REGRESSION MODEL CONSTRUCTION
After obtaining the signals v f , ω f and η f needed to construct the matrices a and b , the discretization of the filtered regression model are defined. By defining g(s) = sf (s) from equation (15), the discrete representation of the filters f (s) and g(s) are given by [25], [57], [58], where z is the z-transform operator. By replacing the discrete filters (35) and (36) into equation (15) the discrete filtered regression model (15) evaluated along the off-line processed signals is defined by where a (v f , ω f ) ∈ R 4×7 and b (v f , ω f , η f ) ∈ R 4×7 are the regression matrices of known functions, and ∈ R 7 is the parameter vector. The implementation of the discrete filtering schemes is performed in MATLAB by using the function  filter(b,a,x), where b and a represents row vectors containing the coefficients for the numerator and denominator either of the transfer functions (35) or (36), respectively, and x represents the signal to be processed. In order to reject the high-frequency noise components from the measured signals and avoiding to lose important information from the system dynamics, the cut-off frequency of the filter must be selected considering the highest frequency from the system excitation signals, see [25] and [58]. In this case, the highest frequency corresponds to the value of the fundamental frequency ω f = 36 [ • /s] times the highest number of harmonics quantity, which is N = 5, thus the highest frequency in the excitation trajectory is 180 [ • /s]. Then, the cut-off frequency of the filters was selected as λ = 30 [rad/s] similarly to [25], which provided a good high-frequency rejection without losing any useful information. Equation (37) is finally defined as where

3) ESTIMATION ALGORITHMS
Once the linear-in-the-parameter reconstruction of the system dynamics is obtained either by the standard regression model in (34) or by the filtered regression model in (38), an estimation of the actual parameter vector should be found. The estimated parameter vector will be denoted asˆ . In this paper, two well-known methods are used to computeˆ : the least-squares and gradient descent schemes.

LEAST-SQUARES ALGORITHM
The least-squares algorithm used for the parameter identification through the models (34) and (38) is given by [66], [68] (40) whereˆ (hT ) ∈ R 7 is the estimation of the real parameter vector ∈ R 7 , ∈ R 4×7 represents the regression matrix (either (v f , ω f , η f ,v f ,ω f ) in (34) or dF (v f , ω f , η f ) in (38)) and ϒ ∈ R 4 represents the system output vector (either F f in (34) or F dF in (38)). In other words, to obtain from the standard regression model (34), the least-squares VOLUME 8, 2020 algorithm is applied with (39) when implementing (40) for the filtered regression model in (38).

GRADIENT DESCENT ALGORITHM
As described in [68], the gradient descent algorithm is a steepest descent approach to minimize the square of the identification error e 2 (t), with the identification error defined as where e ∈ R 4 is the identification error vector,ˆ ∈ R 7 is the estimation of the real parameter vector ∈ R 7 , ∈ R 4×7 represents the regression matrix and ϒ ∈ R 4 represents the system output vector. Since Therefore, update law is given bẏ where σ > 0 is an adaptation gain. Thus, the estimation of the parameter vectorˆ is obtained by the discrete integration of the update law in (41) aŝ For the application of (42) to the standard regression model (34) for the filtered regression model (38), = dF (v f , ω f , η f ), and ϒ = F dF . In this work, the adaptation gain was selected as σ = 30.

III. APPLICATION OF THE IDENTIFICATION PROCESS IN THE QBall 2 QUADROTOR A. EXPERIMENTAL PLATFORM DESCRIPTION
In this document, the proposed identification methodology is applied to the QBall 2 quadrotor from Quanser shown in Figure 2. The controller implementation was performed in MATLAB-Simulink and Quarc software, which are required to develop, compile, and upload the executable code in the on-board computer. The IMU of the QBall 2 consists of a 3-axis accelerometer and a 3-axis gyroscope. The OptiTrack motion capture system available at the Laboratory of Control of the Instituto Politécnico Nacional-CITEDI consists of six synchronized Flex 3 cameras connected to a ground station.
The experiments were carried out using a sampling rate of 500 [Hz] for the inertial measurement unit of the quadrotor and a sampling rate of 30 [Hz] for the motion capture system. The quadrotor signals φ(t) and θ(t) are obtained by the  IMU of the quadrotor. The position p(t) and the angle ψ(t) are obtained by the motion capture system OptiTrack. The nominal parameters of the experimental platform given by the manufacturer are presented in Table 1. Let us notice that in the model (1)-(4), the inertia products I xy , I xz and I yz have been taken into account, in contrast to the nominal parameters in Table 1. As was discussed in [24], [42], [43], [45]- [47], [52], [53] the number of harmonics is directly related to the highest frequency on the excitation signal. By using too many will produce noise, reducing the signal-to-noise ratio, which also increases the number of zero velocity passings and impairs the excitation of the system, and in consequence, the parameter estimation accuracy. Moreover, by increasing the number of harmonics, the optimization of the trajectory is complicated since it increases the number of coefficients to be computed. Thus, the number of harmonics chosen in this work was N = 5. The optimization problem is solved using the fmincon function in MATLAB and the Optimization Toolbox with random initial conditions for a i,l , b i,l , and c i . Figure 3 shows the resulting desired position and orientation trajectories. The obtained coefficients from the optimization process that describes the desired position and orientation trajectories are shown in Table 2.

C. PID CONTROL IMPLEMENTATION
The gains of the PID controller (31)-(32) selected to track the optimized trajectories (18) were obtained by a trial and error procedure until an acceptable tracking was achieved, and are defined as The implementation of the PID controller (31)- (32) ensures that the optimized trajectory (18) is tracked in real-time.

D. PARAMETER IDENTIFICATION RESULTS
The parameter estimation results were obtained using the procedure described in Section II. As suggested in the literature [24], [25], the initial data are not used in the identification procedure. Therefore, the identification procedure was implemented for 10 [s] ≤ t ≤ 50 [s]. This period excludes the samples of the transitory period. The time evolution of the estimated parametersˆ (t) are shown in Figure 4, and Figure 5, where label LS and label GD correspond to the standard regression model in (34) with the least-squares algorithm (40) and the gradient descent algorithm (42), respectively. Labels FMLS and FMGD denotes the least-squares algorithm (40) and the gradient descent algorithm (42), respectively, by using the filtered regression model (38). Table 3 shows the four sets of parameters obtained in the last time sample of the applied methods. It is worth noting that all the estimated parameters remain bounded and seem to converge to some value as time increases. Considering the geometrical characteristics of the quadrotor, and the fact that its mass is concentrated at the origin of the body reference frame, the estimated values of the inertial parameters are consistent.

IV. EXPERIMENTAL VALIDATION OF THE IDENTIFICATION PROCEDURE
Two methods have been used to validate the obtained results [25], [57], [63]. First, simulations have been performed with the parameters provided by the manufacturer shown in Table 1, and with the new estimated parameters  TABLE 3. Estimated parameters of the QBall 2 quadrotor by using the identification procedure of Section II. shown in Table 3. Then, a comparison was made between the simulations and the experiments to determine which set of parameters represents closer the behavior of the quadrotor during the experimental test. The second method computes and compares the output prediction accuracy of the regression models (34) and (38), which is defined as the regression matrix times the estimated parameter vector minus the corresponding output vector. Real-time experiments corresponding to different desired position and orientation trajectories are employed in order to assess the quality of the new estimated parameters, which are described later. In order to save space, only the simulation comparison figures corresponding to the identification experiment are presented. • Signal sampling: The sampling of the inertial sensors and the motion capture system was added by using the ''zero-order hold'' block with the corresponding sampling frequency for each subsystem. The sampling frequency of the position p, the yaw angle ψ, and the yaw angle derivativeψ is 30 [Hz]. The sampling frequency of the roll φ and pitch θ angles and its derivativesφ anḋ θ is 500 [Hz].

B. VALIDATION USING THE IDENTIFICATION EXPERIMENT
In the remaining of this document, the signals obtained experimentally are represented with black lines and the label ''Exp''. The signals obtained in simulation using nominal parameters (provided by the manufacturer) are depicted with red lines and the label ''Sim( nom )''. The signals obtained in simulation using the estimated parameters obtained with the standard regression model (34) and the least-squares algorithm (40) are presented with blue lines and the label ''Sim(ˆ LS )''. The signals obtained in simulation using the estimated parameters obtained with the standard regression model (34) and the gradient descent algorithm (42) are presented with purple lines and the label ''Sim(ˆ GD )''. The signals obtained in simulation using the estimated parameters obtained with the filtered regression model (38) and the least-squares algorithm (40) are presented with green lines and the label ''Sim(ˆ FMLS )''. Finally, the signals obtained in simulation using the estimated parameters obtained with the filtered regression model (38) and the gradient descent algorithm (42) are presented with orange lines and the label ''Sim(ˆ FMGD )''.  Figure 6 shows the path drawn by the quadrotor in the real-time experiment using the PID controller to track the optimized trajectory. Figure 6 also shows the simulation results in the same conditions as the experiment by using the nominal parameters and the parameters obtained by the identification procedure of Section II. Figure 7 presents the time evolution of the position signals p(t) obtained experimentally and by numerical simulations tracking the optimized trajectory. As can be seen, the differences between the experiment and the simulations are relatively small. Figure 8 shows the time evolution of the orientation signals η(t) obtained experimentally and using the simulations when the tracking of the optimized trajectories q d (t) is performed. The similarity of the Euler angles obtained by experiment to the given by simulation is remarkable considering the complexity of the system dynamics and the strong coupling of them with the position signals.
The time evolution of control signals obtained experimentally and the control signals from the simulations recreating the experiment with the optimized trajectory are depicted in Figure 9.
In order to provide a quantitative comparison of the experiments with the simulations, the root mean square value (RMS) of the signals e αnom = α Exp − α Sim( nom ) and e αβ = α Exp − α Sim(ˆ β )  are computed, where α represents any of the signals x, y, z, φ, θ, ψ, F z , τ φ , τ θ , τ ψ , and β represents one of the four sets of parameters given by LS, GD, FMLS, and FMGD.
To quickly identify the improvement of the simulation results using the estimated parameters with respect to the nominal parameters, the relative percentage of improvement (P imp %) was calculated as P imp %(α) = RMS(e αnom ) − RMS(e αβ ) RMS(e αnom ) × 100%, VOLUME 8, 2020  where e αnom is the error between experiment and simulation corresponding to the signal α and the simulation uses the nominal parameters, and e αβ is the error between experiment and simulation for the signal α and the simulation uses one of the sets of the new parameters with β indicating either LS, GD, FMLS or FMGD. The RMS value of each one of the signals e αnom and e αβ and their respective improvement percentages are presented in Table 4. The positive values of the improvement percent (in blue) indicate a favorable result of using the estimated parameters instead of nominal parameters. The negative values of the improvement percent (in red) indicate a reduction in similarity with respect to the experiment for the signals obtained in simulation using the estimated parameter instead of nominal parameters. Notice that the improvement is negative mainly for two signals, specifically for e z and e ψ , while the improvement for the remaining signals is positive, meaning that any of the new sets of identified parameters represents better the behavior of the quadrotor.
The another validation method was described in many parameter identification works [22], [24], [25], [44], [46], [47], [54]. The validation consists of computing the difference of the predicted output of the regression models (34) and (38) obtained in the experimental test. The predicted output is obtained as the product of the regression matrix from the experimental test times by the estimated parameters. Therefore, the output prediction error is given by where ι indicates the set of parameters used for the calculation, that is, LS, GD, FMLS, and FMGD. The matrix (hT ) ∈ R 4×7 represents the regression matrix (either (38)), ϒ(hT ) ∈ R 4 represents the system output vector (either F f in (34) or F dF in (38)), and e yι (hT ) is the output prediction error computed for each sample.
Since different regression models are considered in this work, in order to provide a fair comparison, two sets of the output prediction error are computed, each one related to the implemented regression model. The first set is associated with the standard regression model in (34) and is computed only by using the regression matrix (v f , ω f , η f ,v f ,ω f ) and the corresponding system output vector F f . The set of parameters considered for this calculation are the nominal ones and the obtained through this regression model, which are presented in Table 3, denoted as LS and GD. The second set is related to the filtered regression model in (38) and is computed only by using the regression matrix dF (v f , ω f , η f ) and the corresponding system output vector F dF . The parameters considered are those obtained through the filtered regression model denoted as FMLS and FMGD in Table 3 and the nominal parameters.
The RMS value of each component of the vector e yι defined in (43) is computed for each set of parameters. The percentages of improvement with respect to the results obtained with the nominal values were also computed. The results concerning the standard regression model are presented in Table 5 and the results concerning the filtered regression model are given in Table 6.

C. VALIDATION USING DIFFERENT EXPERIMENTS
Additional to the identification experiment, three different experiments consisting in the tracking of a circular path at different speeds were performed to evaluate the quality of the parameter identification results. The RMS value of the output  prediction errors was computed, including their respective percentage of improvement. The desired position p d (t) and orientation ψ d (t) were defined to produce a circular path in the Cartesian space and designed to be completed in 5 [s] while different radius were specified. The gains of the PID controller (31) and (32) used for the additional experiments are the same as the identification experiment.
The desired trajectories that were used in the experiments of the three circular paths are defined by where ϑ is the radius of the circular path and is a constant used to define the altitude trajectory.   as well, where the values in blue mean more similarity with the output of the regression models using the estimated parameters than the nominal parameters, while the values in red represent the opposite. Only two negative values were obtained in this comparison. In general, all the results represent an improvement in the output prediction; this indicates that the new estimated parameter values capture better the real behavior of the quadrotor.

2) CIRCULAR PATH EXPERIMENT 2
Similarly, another experiment by using the trajectory in (44) was carried out. The radius and constant altitude values selected for the circular path experiment 2 are given by The RMS values of the output prediction errors for the regression models (34) and (38) and as well as their corresponding percentages of improvement are shown in Table 9 and 10, respectively. All the obtained values in this experiment show improvements when using the estimated parameters instead of the nominal parameters. For this experiment, the sets of estimated parameters represent in a better way the real behavior of the quadrotor.

3) CIRCULAR PATH EXPERIMENT 3
The radius and constant altitude values selected for the circular path experiment 3 are given by  Similar to the data presented in Tables 7-10, the corresponding results for the implementation of circular path experiment 3 are presented in Tables 11 and 12. All the percentages of improvement related to the standard regression model are positive, while only one negative value is obtained with the filtered regression model. Therefore, these results indicate that the output prediction results more similar to the experimental ones were obtained using the new estimated parameters instead of the nominal parameters.

V. CONCLUSION
This document describes a procedure for identifying the parameters of quadrotors. The procedure was inspired by manipulator identification literature. The standard and filtered regression models were used. Besides, the least-square algorithm and the gradient descent algorithm were applied. The procedure identified the mass and the inertia moments of the quadrotor. An optimized trajectory to take into account the experimental platform limitations, and to excite the quadrotor dynamics was crucial in the given procedure. The quadrotor parameterized dynamic model in the body reference frame was proposed considering the inertia tensor as a symmetric matrix. Off-line data processing operations to reduce the effect of the noise in the sampled signals and to construct the regression models were described. The presented procedure was validated with different experimental tests and two validation criteria. In all the cases, the results using the estimated parameters turned out to be better than the nominal ones.
In conclusion, the obtained results of the comparisons validate the new estimated parameters, which supports the parameter identification procedure described in this work. Furthermore, the results demonstrate that this method can be used in quadrotors whose dynamic parameters are unknown or have been modified with payloads. In future work, a comparative study with respect to other parameter estimation approaches will provide a better idea of the advantages or shortcomings of the given procedure. Investigation of different excitation signals will be considered, and other criteria for the trajectory optimization will be tested to improve the accuracy of the parameter estimation results.