A Study on Coaxial Quadrotor Model Parameter Estimation: an Application of the Improved Square Root Unscented Kalman Filter

The parametrized model of the Unmanned Aerial Vehicle (UAV) is a crucial part of control algorithms, estimation processes and fault diagnostic systems. Among plenty of available methods for model structure or model parameters estimation, there are a few, which are suitable for nonlinear UAV models. In this work authors propose an estimation method of parameters of the coaxial quadrotor’s orientation model, based on the Square Root Unscented Kalman Filter (SRUKF). The model structure with different aerodynamic aspects is presented. The model is enhanced with various friction types, so it reflects the real quadrotor characteristics more precisely. In order to validate the estimation method, the experiments are conducted in a special hall and essential data is gathered. The research shows that the SRUKF, can provide fast and reliable estimation of the model parameters, however the classic method may lead to serious instabilities. Necessary modifications of the estimation algorithm are included, so the approach is more robust in terms of numerical stability. The resultant method allows for dynamics of selected parameters to be changed and is proved to be adequate for on-line estimation. The studies reveals tracking properties of the algorithm, which makes the method more viable.


Introduction
In recent years Unmanned Aerial Vehicles have been widely used.The most popular UAVs configurations are fixed-wing and multi-rotor [20,21].They are applied in such tasks like: mapping, surveillance and reconnaissance, search and rescue, agriculture and others [30,40].Today one may hear about growing number of projects, where drones play the main role: from the flying cameras, through the personal observer to the sophisticated hardware simulator for the Mars lander [33].Researchers develop improved control mechanisms as well as refined estimation techniques.Currently drone constructors are facing new problems: safety and security.The question is not about how to control the UAV or estimate its parameters, but about reliability and robustness of the proposed systems.
Among many, in literature, one can distinguish three different type of systems: remotely-guided, semi -autonomous, and fully-autonomous.Often these three are distributed into ten-point scale, namely Autonomous Control Level (ACL) [5].In order to achieve higher level of the autonomy it is necessary to incorporate real time self-health diagnostic and mechanisms of adaptation to faults or even failures.As a consequence, specialized control algorithms with adaptive or robust structure are being applied [2,23].In most cases, the control scheme needs prior parametrization, usuallybased on a mathematical model [3,15,32,34].Similar needs may appear during the usage of observers for fault tolerant system in quadrotors (Fig. 1) [17].
Various approaches are used in case of model parameters estimation, namely identification.Since a quadrotor model is highly nonlinear and internally-coupled system which Fig. 1 The Falcon V5 coaxial quadrotor developed at the Institute of Control, Robotics and Information Engineering of Poznan University of Technology has many physical as well as aerodynamic effects, the identification task is challenging.There are two major strategies in system identification: online and offline, both in frequency-domain (FD) and in time-domain (TD).The difference between online and offline [25,36] methods may be understood in terms of needs or availability of onboard computational power, however the choice between time or frequency domain is not yet so trivial.
In frequency domain, e.g. in UAV, the time flight data is converted to frequency domain data.The procedure may by explained based on popular software for identification CIFER ®1 [36,43].Firstly, the flight data is gathered during the flight (CIFER ® is an offline, batch method), then, the frequency response for each input-output pair is calculated together with the coherence coefficients (on this stage the model has already defined structure).Next, the multivariable frequency analysis is carried out; which removes all cross effects in SISO models.Also partial coherence coefficients are calculated.Afterward, the optimization processes are being applied to time window lengths, so the accuracy for low and high frequencies is at a proper level.Finally, the state-space model is identified based on optimization method, driven by frequency response matching.Additionally, the parametrized model is validated at time domain.
Identification in frequency domain is claimed to have several advantages.Firstly, it stem from linear system theory, because of that it gives simple insight in itself and has strong foundations.Also, one has eventuality to choose the range of frequencies in input sequence for the identification (with accordance to an information content and 1 Comprehensive Identification from FrEquency Responses, software developed by the U.S. Army and NASA specifically for rotor-craft applications signal-to-noise ratio).Furthermore, in frequency domain, the model reveals much more information about its characteristics compared to time domain.Therefore it is in general more universal.Models identified in frequency domain are suitable in a case of simulation, control law design or model structure analysis.This approach was successfully used in the case of UAV models in [11] for an insectlike micro aerial vehicle, [10] for fixed-wing UAV or in [19] in a test of flight flutter, namely: vibration, resonance and damping.The most popular methods for an identification in frequency domain are Least Square (LS), recursive Least Squares (RLS) [37,38], Maximum Likelihood (ML) and Recursive Maximum Likelihood (RML) [6,26].
In frequency domain the identification process is able to exploit the MIMO model structure, although each and every sub-model (pair of input-output) must be linear.This is the main feature yet disadvantage of the FD approach.Nonlinear models need to be linearized at an equilibrium point prior to identification.In a case of specialized software like CIFER ® , the identified model will be automatically described by the best linear fit to the response, based on the first harmonic approximation in Fourier series [36].That does not satisfy the requirements for the quadrotor UAV model which is, as already stated, highly nonlinear and has a coupled structure.The model cannot be simply approximated with linear model, and the linearization is only valid for a hoover state.
In some cases that problem can be solved with special techniques.To address this, Adams and Allemang [1] proposed method where model is divided into linear and non-linear parts.In terms of simplicity, the overall method boils down the non-linear state-space model to linear one with the non-linear input function (excitation).The nonlinear input function depends on current input of the prior model and the non-linear function which links the subsequent derivatives of the state.As a result the FNSI gives linear model parameters as well as non-linear model part coefficients.The novel approach was successfully adopted in [37] and [38] under the name: Frequency-domain non Linear Subspace Identification (FNSI).Although the FNSI was proved to be efficient, only one work describes practical use of this method [38].The method can be used to all models with previously mentioned structure, though a conversion from a UAV dynamic model is challenging.The most laborious part of the conversion is to obtain the linear sub-model.In the worst case scenario the linear part is given as an identity matrix, while the nonlinear part forms the overall previously defined model.The FD system identification gives many possibilities, however the cost is disproportionately high.To authors best knowledge, there are not many articles devoted to non-linear quadrotor model identification in frequency domain.In contrast, timedomain methods seem to be more promising.
The most popular and famous method for time-domain identification is Kalman Filter (KF) [29], a recursive algorithm, which works as a two-stages process.In a survey [22] one finds that the KF was successfully applied in various UAVs applications.Primarily, it was designed for linear systems, expressed in a space-state form.Later, its nonlinear version was developed, namely Extended Kalman Filter (EKF).Although the EKF has been used for more than 40 years [12] and has been applied successfully in many works [18,41,44], it is an imprecise method.It suffers from approximations, which are made during the linearization procedure.The EKF approximates mean and covariance up to the first order of the Taylor series.In case of highly non-linear problems it does not reflect the statistical values properly.It is also believed that in extreme cases it exhibits divergence and may be hard to tune when dealing with significant nonlinearities [7].
The EKF was applied for UAV parameters estimation in numerous works.In [28], authors presented results from fixed-wing UAV identification methods, based on EKF and Error Mapping Identification (EMId).In order to verify proposed algorithm, the basic 3-DOF UAV model was given with the assumption that only angular velocities are measured.The experiment was carried out in hardware in the loop (HIL) and in simulations.The obtained results indicate, that the EKF is stable and its uncertainties are at bounded levels (evaluated based on main diagonal elements of the EKF covariance matrix).The EKF proved its effectiveness at presence of partially corrupted sensory data and lower sampling rate.In [8], authors have dealt with online parameters estimation of Vertical Take Off and Landing (VTOL) UAV, by use of the EKF.Identification results were used for Riccatti equation solver in control procedure, providing optimal control even with prominent sensors noise or biases.Common and noticeable for all of the previous works is superiority of the EKF method, however, as already mentioned, the EKF may fail in case of highly non-linear and coupled models.That was partially proven in [7], where authors have proposed two Kalman Filter versions: the Extended and the Unscented.
The Unscented Kalman Filter (UKF) was firstly proposed by Julier and Uhlmann [27] as a solution for EKF problems.The main concept behind the UKF is Unscented Transformation (UT).During the Unscented Transformation a set of carefully chosen points (actually vectors) called sigma points are passed through the nonlinear function.The set of sigma points are chosen in that way, so the mean and covariance of the transformed function can be captured properly.The UT uses the re-sampling technique, so the statistic is calculated based on the set of points.Technically, if the noise distribution is Gaussian, the UKF provides approximation up to the third order, if the input is non-Gaussian the filter still holds the approximation accuracy at least to the second order of the Taylor series [35].The UKF has been applied in different works, focusing mainly on attitude estimation [9], position estimation [42] or parameters estimation [24], though the last one remains practically unexplored.According to [22], the UKF in identification was used only in one work [4].In [4] Brunke and Campbell have introduced the architecture for a model and state (joint) estimation of a fixed-wing plane.It is important to note that, the solution proposed by them was preliminary explained and proven.Firstly the authors have shown that the problem is non-linear and two most popular techniques including Extended Kalman Filtering and Maximum Likelihood are inadequate.Also by citing different authors, arguments for and against the EKF and the proposed UKF were presented.In conclusion Brunke and Campbell have suggested that the most reliable, viable, elegant and numerically stable is the Square Root version of the UKF (SRUKF), proposed by Merwe and Wan [35].In addition, the authors have presented promising results obtained via the SRUKF method.
It is important to remind, that completely different way, which adds an extension to the well-known time and frequency domain horizon also exist.The most popular among others are heuristic techniques.In their scope, two main types of algorithms are placed: evolutionary and swarm.In [45] Yang and Liu used the Particle Swarm Optimization (PSO) method in a quadrotor UAV inertia parameters estimation.The research involved simple model derivation without aerodynamic effects.Also an interesting framework for PSO-based identification was given.The authors emphasized the problem associated with higher dimensions PSO (falling into the local optimum trap).In order to prevent it, two parallel PSOs were used in the same model identification procedure.This ultimate approach seems to be adequate in case of multimodal problems.However, the authors have not mentioned about model's non-stationarity.In different work regarding the PSO in an identification of UAV parameters [16] one finds similar model, enhanced with aerodynamic effects, though limited to one axis.The authors have given comprehensive study on the PSO, with different kind of cost functions and various model and parameters structures.They proved also that only switched models can be used, since the classic models with assumption of parameters constancy cannot properly reflect quadrotor's behavior.Furthermore, when dealing with switched models, it is crucial to introduce few of them as the accuracy needs to be at appropriate level.This involves increased number of parameters to be found and significantly extends the PSO's computational time.In the literature, one will find numerous applications where evolutionary or swarm algorithms are used for identification or simply model parameters estimation, however none of them are suitable for online tracking in case of non-stationary models (numerous parameters, computational burden).
To conclude: in order to identify the model parameters of a quadrotor UAV the time-domain method should be choose.Among many available methods, only a few are suitable for highly non-linear and coupled models.The most popular estimation methods are within the scope of Kalman Filters, where the highest accuracy, robustness and stability are expected to be achieved by the Square Root version of the Unscented Kalman Filter (SRUKF) [35].
In this work, authors have proposed quadrotor model described in [15] (adapted to coaxial quadrotor schame), enhanced with additional effects: flapping, induced, translational and profile drags and friction model.Then the identification method based on Square Root Unscented Kalman Filter (SRUKF) is presented as well as the results obtained via this method.Furthermore, the SRUKF method, proposed by Merwe and Wan [35] is modified, so it can be more robust in terms of numerical stability -regardless of initial settings and conditions.Also the universal framework for non-linear dynamic models parameters identification is given.Studies on different model structure and identifier input-output dependency are also presented.In addition, the comparison of identification results for the same model, based on PSO method [16] and the SRUKF method is included.The novelty of the paper lies in the improved SRUKF method applied to online model parameters estimation, which is assumed to be non-stationary.It also lies in the presented framework, as it gives capabilities in the field of identification, especially in UAV.To our best knowledge there are no works where an issue of online UAV quadrotor non-stationary models parameters estimation was raised.
The article is organized as follows: in the second section the coaxial quadrotor model is derived and described.In the third section, the SRUKF method as well as its modifications are presented.In the fourth section the experiment and all assumptions undertaken prior to the experiment are described.Finally, in the last section; results, conclusions and summary are presented.

The Coaxial Quadrotor Model
The basic quadrotor's model has six degrees of freedom (DOFs) (zero-order derivatives).First three DOFs describe position, the rest denotes orientation (Euler or Tait-Bryan angles) expressed in global reference system {EF } (Robot local reference system is denoted as {BF }).The global {EF } reference is configured as the NED (North-East-Down).In the paper, authors derive the model of the coaxial quadrotor, which has eight propellers mounted co-axially on a frame of the quadrotor (Fig. 1).Compared to the classic quadrotor, this leads to the similar mechanical structure with increased payload.It is often called octarotor, though it may be confusing since octarotors have in general planar mechanical structure.In order to simplify the vocabulary, authors use the coaxial quadrotor as a name of that sort of robots.In addition, authors would like to point out that planar octarotors are more efficient, but consume more space and are larger.Coaxial quadrotors provide a reasonable trade off between size and payload, especially when applied indoors.In the paper, the discussion is narrowed to the orientation part of the coaxial quadrotor model.The basic equation, which describes the angular acceleration (for convention from the Fig. 2) has the following structure: where the R : {BF } → {EF }, ω BF -angular velocity given in {BF }, [] × -skew symmetric span, m -robot total mass, I -robot inertia matrix, τ -torques acting on the robot, vBF -linear velocity (∈ {BF }), g -gravitational acceleration, f -forces acting on the robot.
It is important to note, that, in order to retrieve the robot's orientation in {EF }, one is forced to integrate angular acceleration ( ωBF ) and transform the result into the Euler rates (4): Although other methods exist, including quaternions, this one remains simple, yet intuitive.Second integration of the Euler rates gives the Euler (or Tait-Bryan) angles.Matrix P EF BF is defined in Eq. 5 (detailed in [15]): The overall definition of P EF BF depends on the Euler or Tait-Bryan representation, here the "3-2-1" convention is used ("flight notation" with roll, pitch and yaw angles).The propellers generate forces, which (in simplest version) can be expressed as product of scalar parameter b i and a second power of rotor's angular velocity ω i , thus the force due to the i-th rotor may be given as in Eq. 6: Thrust is generated only in Z-axis of the local reference system {BF }.The input torque τ is a sum of torques, due to subsequent propeller.The i-th torque is a cross product of the i-th propeller position in respect to the {BF } and the i-th force, Eq. 7.
where l w is equal to 2 l, l stands for arm's length in the quadrotor cross frame and h is a distance between propeller's plane and the center of mass in respect to Z axis of the {BF }.The Eq. 7 covers torques generated in X and Y axis of the {BF }, while the Z axis suffers from the reaction torques, given as: ). (10) In Eq. 10 parameter d i refers to reaction torque gain.Two of them are introduced, which stems from the opposite type of propellers system.Propeller 1, 4, 5 and 8 have clockwise rotation and the rest are counter-wise.This simple notation can be extended to additional effect like different types of drags.The better insight is presented in [3].As it can be found, there are few types of draglike effects: blade flapping, induced drag, translational drag, profile drag and parasitic drag (which for speeds up to 10 m/s may be ignored [3]).
Blade flapping effect is ubiquitous in rotor vehicles.It occurs when the rotor undergoes translational motion.It simply changes the rotor's tip path plane, deflecting it by the flapping angle β f lap .Here, the explanation is narrowed only to the final form of the force generated by the flapping effect, though one may find detailed explanation in [31].Flapping force due to the i-th rotor is given in Eq. 11: where V p is a horizontal velocity in the {BF } (V p = V x V y 0 T ), while the A f lap and B f lap are the matrices with parameters to be identified: where r refers to rotor radius.Another effect is induced drag.It may occur in small quadrotors, where rotors blades are rigid or semi-rigid which cause an unbalanced in flapping forces.This effect can be simply described as a product of planar velocity and a scalar parameter: Similarly to induced drag, it is also possible to distinguish translational drag.This effect occurs, when the air is redirected from any translational movement to the downward of the {BF }.In accordance with [3], this effect depends on a speed in the air.For slow movements, the translational drag is equal to: while for fast maneuvers, it is necessary to incorporate the air velocity v i induced by i-th rotor: where the K T 1 and the K T 2 are scalar parameters.The last two drag effects are profile drag and parasitic drag.Since the parasitic drag comes from non-lifting surfaces of a quadrotor and occurs only for high speeds (> 10m/s) it may be ignored [3].Profile drag is gained with the transverse velocity of the rotor blades in the air: where K P is a lumped parameter.Finally, the total drag force acting on the i-th rotor can be formulated as a sum of individual drags: while the torque due to drag effect, classified for i-th rotor is defined as: The net torque τ D , added to model ( 1) is a sum of all τ D i , where i ∈< 1 : 8 >.In the model, the friction is composed of three main elements: where the h g , b w and c w are friction coefficients, while stands for Hadamarda product.
Final equation for angular acceleration is given as: where τ gyro is a gyroscopic torque due to the rotors and frame velocity.

Model Parameter Estimation -the SRUKF Method
In the Introduction comprehensive study on different methods which are adopted for model parameter estimation was given.In conclusion, authors stressed the role of the time-domain methods based on Kalman Filters by showing their superiority.Among them, probably the most interesting one is the Square Root Unscented Kalman Filter.It was firstly proposed by Merwe and Wan [35] as an extension to well known Unscented Kalman Filter method [27], though it has to be noted that some first attempts in Square-Root filters date back to 1995 [39].As all Kalman Filter methods, the SRUKF is a recursive twostages procedure and in case of parameter estimation it has simplified structure (compared to the state estimation variant).It is important to recall the SRUKF equations since the algorithm is partially changed.The parameter estimator or observer operates on the basic model, given as: where w k are the parameters to be found, r k refers to process noise, G(.) is the non-linear mapping function, dk is a desired output and e k constitutes an error.The notation ( 22)-( 23), clarifies the roles which are assigned to the subsequent vectors and function, however it has to be explained where the non-linear UAV model (discretized by the forward Euler method) is placed.First of all, one needs to assess the measurement system applied in a UAV.Based on that it is possible to choose the observed values.The UAV is equipped with the Inertial Measurement Unit (IMU) or Attitude and Heading Reference System (AHRS) and in case of orientation these quantities are available: angular velocities, angles or quaternions.The model given in Eq.21 has direct impact on the angular acceleration, thus it is highly important to observe or measure this quantity.In the best case scenario, the d k represents the angular acceleration, while the G(.) transfers the previous state vector and parameters to the desired angular acceleration, however dk may also include angular velocities and angles or quaternions.
where λ RLS and γ are internal parameters (to be set based on a process knowledge).In the second stage (measurement phase), the ( 31)-( 37) are used: where L refers to the number of parameters, W (c) i and W (m) i are the method's parameters which need to be properly chosen and R e is a measurement noise -according to [35] and can be set arbitrary, though it has a big influence on convergence speed (to be proven in the next section).Behind the SRUKF, there are few different parameters that rule the algorithm (extensive description in [14,35]).The first one is α, it determines the spread of sigma points around the estimated parameters (eg.value: 0.7).The second is β, which describes the distribution (eg.value: 2.0).Finally, there is the scaling parameter called κ (eg.value: 0 for state-estimation and 3 − L for parameter estimation [35]).The previously listed parameters depend on the α, β and κ in a given form: where i = 1, ..., 2L.The SRUKF algorithm given in Eqs.26-37 is controlled via the λ RLS , which has the same function as in case of Recursive Least Squares Method.It simply introduces the forgetting feature to the estimation process.It is commonly used and set close to 1 (e.g.0.9995).
As it came out during the first investigation, the λ RLS may lead to significant numerical instabilities.In most cases, the safety could be guaranteed by setting the λ RLS to 0.99999, but it was not necessarily the optimal.In order to make the method more robust and less sensitive on the initial setting, the λ RLS is changed in each iteration so it can reach 1.0 in a certain time.Due to that, the algorithm has to be extended with the following formula: where the constant λ 0 change the characteristics of the λ RLS k in the subsequent samples.The influence of the λ 0 on the λ RLS k is depicted in the Fig. 3, for different λ 0 .Further modifications of the SRUKF method address the problem of the parameter convergence speed and confidence for the measurements or reference given in d k .
When the correction phase of the algorithm is invoked, results from the model given in dk are subtracted from the reference d k (35).If the reference introduces errors, the corresponding variance in the covariance matrix R e has to be increased.Although, it has strong influence on a speed of the parameters convergence, Authors of the [35] recommend to set the R e arbitrarily (e.g. to 5I).That makes parameters partially immune to the errors added in the reference.However, high covariance values slow down the convergence speed.Because of this, the fast-changing parameters cannot be tracked well -one will find it difficult to properly choose the values of the R e .In order to find the optimal settings of R e , the optimization is needed.
Although the R e has strong influence on the parameters as well as convergence characteristics of the algorithm, it cannot change parameter characteristic individually.The influence of the R e remains implicit.For example, if d k , λ 0 = 0.9 Fig. 3 The influence of the λ 0 on the λ RLS k dk ∈ R 1×1 , and w k ∈ R L×1 , where L > 1, the R e ∈ R 1×1 changes all the parameters equally.As it can be seen, the user cannot incorporate the prior knowledge about the parameters behavior into the estimation process.If one of the parameters is supposed to change slowly, it would be desirable to prepare the SRUKF algorithm, so it can cut off all sudden changes in the subsequent track.Unfortunately the original form of the algorithm, given in Eqs.26-37 does not provide such tune possibilities.In order to add this, one simple operation is need to be improved.In Eq. 27, the S − w k is enhanced with the square root of diagonal covariance matrix Q, yielding the following: This simple modification stem from the state estimation algorithm, where Q is used as a process noise.In that case, the Q ∈ R L×L .Note, that L parameters have the corresponding constant variance on the main diagonal of the Q.By increasing subsequent variances, the user is able to set the parameters behavior to be volatile.The tracking capabilities are adjusted by Q and R e , however only its ratio explains how fast or how slow the tracking will be carried on.The settings of Q as well as R e and the result achieved by use of the modified SRUKF method are shown and discussed in the next section.

The Experiment and Results
The paper concerns an orientation model of the coaxial quadrotor identified by the SRUKF method.In the previous sections the model and the method with improvements were shown.In this section, authors present experiments, restrictions (that are imposed) as well as the results achieved with the proposed method.

Assumptions and Restrictions
It is important to recall that the model given in Eq.21 is in a continuous form, which cannot be applied directly in the SRUKF.The overall discrete form (discretized by the forward Euler method) was given as in Eq. 23, however that does not suits to the general state-space model, where the state vector x k has to be updated in each and every step.In order to make these two forms equally, one has to introduce the measurement system as the reference.As it has been already stated, the reference d k is formed from the state vector, where all angular accelerations, velocities and angles are captured (provided by the AHRS system).The aim is to minimize the error given in Eq. 45: where the dk−1 is equal to x k .If the number of sigma points is 2L + 1, there are 2L + 1 vectors in D i,k|k−1 , where all are calculated based on x k and W k|k−1 .However the x k has to be calculated in the previous iteration.The basic question is how to calculate the x k so, that the state vector will form an optimal trajectory.The solution can be provided based upon the update procedure, which allows to compute the model state vector in k + 1 iteration i.e. x k+1 , when the x k and w k are known: The procedures sequence of model updating and its parameters estimation is shown in the Fig. 4. As it can be seen, the model is updated after the correction phase of the SRUKF, based on the estimated parameters vector ŵk .The reference consists of three main components, i.e. angular accelerations, angular velocities and angles.The most important and required is the first one.Without the angular accelerations it is not possible to perform the estimation process.It is caused by the direct influence of the angular momentum on angular acceleration.The angular velocities and angels are less important.First tests proved that the model where only velocities and angles are used (as a reference) cannot be correctly identified.The SRUKF failed to operate in all cases, where the angular acceleration was not used.The reason is simple, though disputable.The most important part of the SRUKF lies in the calculation of the covariance matrix (actually its decomposition product).If the impact of the parameters to be estimated on the output Fig. 4 The SRUKF procedure and the model update of the model is slender, the covariance matrix is minimize to close to zeros and the resulting Kalman gain is increased.Finally the SRUKF becomes unstable and spoils parameter estimation process.As a matter of fact, low influence of the parameters on the model output is mostly caused by the integration process that occurs when calculating the angular velocities from the acceleration.Another important factor is coupling that exists between the velocities and acceleration.When calculating the integral, the angular velocities are obtained, however in the next iteration those velocities are used in the accelerations equations, where the gyroscopic torques are calculated.This may lead to model instabilities as the higher derivatives are not controlled directly.To handle this, one may set the initial parameters as close to optimal as possible, unfortunately most frequently these are unknown.The second option is to control highest derivatives, i.e. angular accelerations.
The final restriction refers to traceability of the parameters.The model given in Eq.21, cannot be fully-identified as there are infinity sets of parameters.It is why, some of them must be given prior to an identification.The first set includes the inertia tensor parameters (six elements of the symmetric matrix I), the second group consists of the rotorpropellers gains b i and d i .The propellers parameter or the gain will probably change during the flight, while the inertia parameters remain constant, this is way the authors have decided to seek for the gain parameters instead of inertia parameters, which are determined before the identification.

The Experiment, Gathered Data
The UAV used in the research (Fig. 5) was equipped with the AHRS (ADIS16488 by Analog Devices) recording samples at a frequency of 400 Hz.It provides the state estimation based on the built-in EKF estimator (400 Hz frequency), however the filter settings are not adjustable.That is why the custom algorithm that provides unbiased angular velocities Fig. 5 The Falcon V5 coaxial quadrotor (start from a safety net) J Intell Robot Syst (2019) 95: -0 491 51 and angles was used instead.The authors have used the sport hall (Fig. 6) as the air disturbances (assuming that the distance from the ceiling, floor and walls is sufficiently high) are lowest and the wind effects can be neglected.
It is important to notice that the research concerns estimation in open-loop system, where in the control the commands were sent with a period of 30ms.A detailed verification of the balance between the control and estimation frequencies (in a closed-loop system with observer) will be provided in next (separate) publication.
As it has already been stated, the reference consists of angular accelerations, velocities and angles (X,Y and Z axis), where accelerations are necessary.Unfortunately, in a great majority of AHRS, the angular acceleration is not available.In order to acquire it, the velocity derivative has to be calculated.This method has two disadvantages, the first one is the undesirably noise gain, and the second one is an introduced delay (at least one sample).In order to minimize the first effect, authors have applied the first-order, low-pass Butterworth FIR filter to velocities measurements (cut-off frequency equal to 0.1 of half of the normalized sampling rate).It strongly reduces the errors, though the same filter must be applied for all derivatives (same delay).As it came out, the introduced delay in angular acceleration has minor influence and could be neglected.
During the experiment the UAV coaxial quadrotor has operated in the closed area in all three axes with low linear speed (Figs.7-8).The date gathered during the experiment consists of set values for BLDC motors (for eight propellers), estimated angles, speed in Z axis (∈ {EF }), estimated height (∈ {EF }), angular velocities (gyroscopebased), linear accelerations, magnetometers readings and the sample time (sample rate equal to 400Hz).
One of the most important aspects of the measurement system, is an unavailability of the propeller speed.These velocities are not measured, however they are needed as they are used in the force calculation.In order to possess the knowledge about the rotors speed, necessary calculation Fig. 6 The sport hall and the UAV coaxial quadrotor Fig. 7 3-D plot of the flight trajectory must be made.The rotor-propeller model can be given as the second-order inertial object [13], where the gain is included in the force (6) as b i .It is important to note that, the second order inertia relation describes how the BLDC controller's set value (duty cycle in %) affects the propeller's speed.This simple procedure allows to estimate the rotors velocities.

The Results
The experiment involved coaxial quadrotor with main board equipped with STM microprocessor and internal memory.Data collected during flights was stored in the memory.After a flight the data was downloaded on a stationary computer and was used in further processing.The SRUKF estimation algorithm was used together with the discrete version of the model given in Eq.21, fed with gathered and filtered data.The discrete version of the model was obtained via explicit forward Euler method.Similarly, the discrete integration was realized with the same approximation.
The discrete model of the coaxial quadrotor has sixteen gains b i in X and Y axes (eight per axis, two per propeller), eight gains d i in Z axis (one per propeller), nine friction coefficients hg x , bw x , cw x , hg y , bw y , cw z , hg z , bw z , cw z (three per axis, for first, second and third power of the velocity in the subsequent axis), two coefficients of the flapping effect B 1 , B 2 and six offsets for velocities and angles in X, Y and Z axis.The quadrotor has also constant parameters, which are calculated or measured prior to the experiment (Table 1).
After the first estimation trial it was observed that some of the parameters are redundant.The coaxial propellers have the same gain, which reduces the number of b i to eight.Similar applies to d i , where only four of them are necessary.The final equation for τ has a following form: As it can be seen, the coaxial propellers pairs in X and Y axis differ with the gains.This is crucial, as the model with further simplification cannot be identified properly.
The most important question is why some of the naturally or theoretically equal parameters like b 1 , b 5 or b 2 , b 6 cannot be simplified?The answer may lie in the side effects when maneuver in the air.Some of the drag-like effects (induced, translational) are not accounted as the linear velocity is not available.Also the propellers may deflect, which results in different torques in X and Y axis.In the Z axis, that problem did not occur.The d i parameters cover the doublets of upper and lower propellers with the same rotary direction.
During the first stage of the research, data gathered in the experiment was filtrated and used in the standard SRUKF algorithm.The results are shown in Figs. 9, 10, 11 and 12.The angular acceleration is well reconstructed.The difference between the model and the reference can be captured only during the first 20 ms.However these 20 ms has a huge impact on angular velocity which results in the offset. 2Similar situation occurred in case of the lowest derivations, namely Tait-Bryan angles.Although the offsets were estimated, the SRUKF could not follow its changes sufficiently fast.Authors have tested the algorithm with different settings for R e without positive results.Moreover, the algorithm was stable only for R e ii ∈< 0.01; 9 >.In order to intentionally change the convergence for offsets and to boost the robustness of the proposed method, modification of the SRUKF was added.The √ Q was modified, so its main diagonal elements responsible for the noise of the offset parameters were set arbitrary to 0.01.The results are shown in Figs. 13, 14, 15 and 16.The applied modifications have allowed to follow the offset properly.The estimation process takes less than 0.8s to fully reconstruct the reference.It also reveals high stability, as the R e was changed from small to large values (tested within the range: < 0.001; 50 >).As the R e was decreased, the method started to be sensitive for small changes in the reference signal.The advantage of the low values in R e can be visible as fast convergence, however this may lead to errors while some sudden changes occur.For example, when the reference is distorted by state estimation or measurement noise, the SRUKF apace follows this noise by appropriate changes in parameters, which definitely corrupts the estimation.Although reliable method for assessment of the hyperactive behavior of parameters estimator do not exist, one is able to intentionally introduced noise which should not have influence on parameters.Authors have added the noise as a peak value in the An adverse effect is clearly observed for R e ii = 0.005.The most interesting case is depicted in the Figs.19 and 20, where in 150-th sample of data, sudden reference noise occurred.It shows how simple data corruption changes the overall characteristics of the estimation.As it can be deduced the process of the algorithm tuning requires optimization.Authors have decided to use the PSO algorithm with only two parameters to be tuned-optimized.The first one is the R e ii , the second is Q ii .An optimization must be performed in presence of noise, that have significant influence on parameters estimation for small R e ii .A trade off is made between the SRUKF's hyperactivity and the convergence speed.The PSO's settings are described in [16].In order to introduce fair evaluating factor, the cost function based on sum of the squared errors (SSE) for acceleration, velocity and angles is adopted.One of the PSO simplification is a length of input data, which was shortened to 200 samples.This not only reduces an optimization time, but also gives possibility to choose the most interesting scope of flight data.The obtained trajectory of particles in a chase for the global optimum is shown in the Figs.21 and 22.
The SRUKF settings found via PSO satisfy the estimation performance, though the optimal performance depends on many variables, which cannot be provided  2. The best performance is achieved for R e ii = 5.0658, √ Q ii = 0.0978 and R e ii = 6.4941, √ Q ii = 0.0077 (data disrupted).Interestingly, the best settings for disrupted reference data, which were found via PSO, have low sum of errors (22.884) for undistorted data.Therefore settings: R e ii = 6.4941, √ Q ii = 0.0077 become optimal for SRUKF parameters estimator.Note that the cost function given in Figs.21 and 22 are significantly different as the PSO optimization operated only on first two hundred samples, while the err sum was calculated for one thousand input sequences.

Discussion and Conclusions
In the article the method of coaxial quadrotor model parameter estimation was presented.The method is based on Square Root Unscented Kalman Filter which has better stability compared to Unscented Kalman Filter (UKF).

Fig. 2
Fig. 2 The quadrotor skeleton and kinematic model (propellers marked by light gray circles have clockwise rotation)

Fig. 10
Fig. 10 Motor-Propellers gains b i , d i

Table 1
Parameters of the coaxial quadrotor