Correction to the Euler Lagrange Multirotor Model with Euler Angles Generalized Coordinates

This technical note proves analytically how the exact equivalence of the Newton-Euler and Euler-Lagrange modeling formulations as applied to multirotor UAVs is achieved. This is done by deriving a revised Euler-Lagrange multirotor attitude dynamics model. A review of the published literature reveals that the commonly adopted Euler-Lagrange multiro-tor dynamics model is equivalent to the Newton-Euler model only when it comes to the position dynamics, but not in the attitude dynamics. Step-by-step derivations and calculations are provided to show how modeling equivalence to the Newton-Euler formulation is proven. The modeling equivalence is then verified by obtaining identical results in numerical simulation studies. Simulation results also illustrate that when using the revised model for feedback linearization, controller stability at high gains is improved.


INTRODUCTION
Derivation of an accurate mathematical model is essential and prerequisite to model-based control of complex dynamic systems in general, and of multirotor UAVs in particular.For example, when considering dynamic inversion control, an accurate mathematical model of the system under consideration should allow for compensation of (any) nonlinear effects, and for linearization of the system dynamics.However, as it happens in almost all cases, it is not realistic to expect complete 'capture' of all nonlinear system dynamics effects through a mathematical model, even though the derived model itself does contribute to achieving desirable performance.
When focusing on multirotor UAVs, the Newton-Euler (N-E) and Euler-Lagrange (E-L) formulations are the two main modeling approaches, albeit following different principles; that is, balancing of forces, and the principle of least action, respectively.Regardless, the N-E and E-L formulations are equivalent, and when applied and implemented on multirotor UAVs, they should return identical results.
However, after a thorough literature review, this work reveals that when it comes to the published E-L formulations, when substituting the transformation matrix from the body-fixed frame angular velocity to Euler angle derivatives, the N-E and E-L attitude models are not equivalent -this is shown in detail in Section III.Hence, the motivation and main objective of this paper is twofold: First, derive a revised E-L (r-E-L) attitude dynamics model for multirotors, registering at the same time the main differences with the E-L formulations used in literature.The proposed r-E-L model is based on the one presented in [1].Second, prove analytically the r-E-L model's equivalence to the N-E formulation (position and attitude dynamics).Then, to show the implementation improvements with the adoption of the r-E-L, the modeling equivalence is demonstrated through numerical simulations on quadrotors.
The quadrotor is studied because it is the most widely used configuration of a multirotor UAV.Its mathematical model description may be found in [2][3][4], where both formulations are detailed.To be specific, the first E-L formulations may be found in [5][6][7][8][9][10], while details of the N-E formulation are presented in [11,12].Both formulations have been widely used and have been implemented for model-based control and navigation.
State of the art E-L quadrotor models with global validity [13][14][15] are not affected by the findings in this paper due to the coordinate-free nature of their formulations.Nevertheless, the Euler angle variant of the E-L quadrotor model is still adopted in literature work [16,17], which supports further the correction proposed in this paper.The findings of this work are coherent to the general attitude E-L formulation of [18], which was published at the time of writing this paper.Although following different approaches, both papers arrive, independently, at the same result.
The rest of the paper is organized as follows.Section 2 introduces the required notation, and the N-E and E-L quadrotor dynamics models as found in the related literature.Section 3 details the proposed r-E-L model and proves its equivalence to the existing N-E model found in the literature.Section 4 includes simulation results.Controller performance comparisons between the proposed and existing models verify and illustrate the equivalence between the two modeling approaches as presented in this paper.Lastly, in Section 5, conclusions are offered.
2 Notation and Background Information is the 3 × 3 skew symmetric matrix composed of the elements of a.In addition, consider the unit vector e 3 = [0, 0, 1] T .The notation da dt and ȧ is used interchangeably throughout the paper.In what follows, for clarity purposes, the N-E and E-L quadrotor dynamics are considered, which may be easily generalized to any multirotor UAV dynamics.

Quadrotor Nonlinear Dynamics
The N-E quadrotor dynamics, as presented in [2], are described by the following equations where S(ω) ∈ R 3×3 is the skew symmetric matrix of the angular velocities, ω, defined in the body-fixed frame, J ∈ R 3×3 is the constant diagonal inertia matrix, M is the external torque induced by the quadrotor propellers in the body-fixed frame, v is the quadrotor linear velocity vector in the body-fixed reference frame, m is the total mass of the quadrotor, g is the gravitational acceleration, and T is the total produced thrust.The matrix R(η) ∈ R 3×3 represents the rotation from the body-fixed frame to the inertial frame.Note that the choice of R(η) is not unique since the quadrotor dynamics are invariant to any choice of Euler angles configuration that may be used to represent the attitude of the quadrotor.The E-L formulation, as introduced in [19], is represented by the following two equations where η = [φ , θ , ψ] T is the vector of any choice of Euler angles configuration, and p = [x, y, z] T is the inertial reference frame position vector.J R (η) is the rotated inertia matrix and C(η, η) is the matrix accounting for centrifugal and Coriolis effects.

Equivalence of the N-E and E-L Modeling Formulations
To verify if the E-L and N-E models are equivalent, a coordinate transformation should lead from one formulation to the other.Starting from the position dynamics, the linear velocity Thus, by substituting ( 6) in ( 3), the following equation is derived (with all steps shown in detail) @ @ @ @ @ @ @ S(ω) T R T (η) ṗ + R T (η) p = @ @ @ @ @ @ @ −S(ω)R T (η The resulting equation is the same as (5), thus, the equivalence between the N-E and E-L position dynamics formulations (as found in the literature) is proven.However, the relationship between the Euler angles derivatives η and the angular velocities ω is less intuitive.This relationship may be derived by considering the rotation from the inertial reference frame (IF) to the body-fixed frame (BF) as the composition of three elementary rotations.To achieve this, two intermediate reference frames are defined, F1 and F2, respectively.Following [20], without loss of generality, to show how this relationship is derived, the Tait-Bryan 321 sequence of rotations is employed as shown in Fig. 1, and it is detailed next.

• IF−→F1
F1 is obtained from the elementary rotation of an angle ψ with respect to the z-axis of the IF, defined as R 3 (ψ) T = R 3 (−ψ).The related angular velocity in the IF is ω ψ,I = [0, 0, ψ] T .It follows that the angular velocity in F1 is given by ω F2 is obtained from the elementary rotation of an angle θ with respect to the y 1 -axis of F1, defined as R 2 (θ ) T = R 2 (−θ ).The related angular velocity in F1 is ω θ ,1 = 0, θ , 0 T .It follows that the angular velocity in F2 is given by ω θ ,2 = R 2 (−θ )ω θ ,1 = ω θ ,1 .Moreover, in F2, the angular velocity related to the angle BF is obtained from the elementary rotation of an angle φ with respect to the x 2 -axis of F2, defined as R 1 (φ follows that the angular velocity in BF is given by ω Moreover, in BF, the angular velocity related to the angle Finally, in BF, the angular velocity related to the angle Therefore, the angular velocity as expressed in the BF is where Thus, in more compact form where W depends on the choice of the Euler angle sequence.One example, for clarification purposes, may be found in [21].Kinematic relations ( 9) and ( 10) are invariant to the choice of the Euler angles, but the configuration needs to be consistent with the choice of R.However, different Euler angles configurations will result in a different W , and so, given (10), it is recommended to choose a configuration in which W is invertible for the whole flight envelope, see [19], except when the pitch angle is equal to π (an uncommon state outside acrobatic manoeuvres).
When substituting (9) in the N-E attitude model (2), this leads to where S(W η) is the skew symmetric matrix of W (η) η.
When compared to the literature E-L model shown in (4), which can be rewritten as it leads to the following inequalities and, in general, contrary to what should be expected, the multirotor attitude dynamics following the E-L and N-E formulations found in the literature do not produce an equivalent result given that However, this shortcoming may be rectified based on the approach introduced in [1], in which a proof is provided for the equivalence of the projective N-E equations and the E-L equations of second kind for spatial rigid multibody systems.
Considering the proof in [1], it is shown that for multirotor UAVs (including quadrotors), the attitude dynamics E-L equations should be written in the following form instead of the one presented in [2], that is with the Lagrangian, L = 1 2 ω T Jω, which, in the case of attitude dynamics, is equivalent to the rotational kinetic energy.
An intuitive reason for writing the r-E-L form as ( 15) is given by noticing that a premultiplication of W T would result in ( 12) and ( 13) being equalities.Moreover, this approach is similar to the quaternion variant of the E-L formulation as shown in [22].Building on the above observation, it is now shown that the same steps presented in [1] may be followed to model multirotor UAVs / quadrotors.Since it is considered that the forces are applied to the center of mass of the quadrotor (multirotor), the position and attitude dynamics can be analyzed independently.Hence, only the E-L equations for the angular dynamics are derived.
To prove the equivalence of the r-E-L model, to the N-E model, the following relations are defined first.For better readability, W is written without explicit dependence on η.
Relation 1.Consider W (η) −1 , the rows of which are given by w inv,1 , w inv,2 , w inv,3 , as follows Next, define the matrix This matrix, for any Euler angle sequence, is composed of the skew symmetric matrices of S(w inv,i ) with i = 1, 2, 3.
Relation 2. Note that the time derivative of W −1 can be expressed as Relation 3. The following is obvious Relation 4. From Rel. 3, the following holds Relation 5.In addition, using Rel.3, it is true that that leads to the following Relation 6. Rearranging Rel. 5 and substituting ( 19), (21), and ( 18), the following equality is obtained leading to Relation 7. Considering (9) it is easy to show that Given the relationships above, it is shown next that the r-E-L formulation leads to an equivalent result with the N-E attitude model.
Proof.The multirotor r-E-L formulation (15) can be rewritten as Since J is a constant symmetric matrix, this is equivalent to From (28), by using Rel.6 and Rel.7 one obtains And, for W having full rank, (29) is rewritten as This corresponds to the N-E quadrotor model formulation.This equivalence does not hold when using the E-L formulation from the literature ( 16).This concludes the proof.
Therefore, the proposed attitude r-E-L model is which rectifies the literature E-L model (11) and this result is reconfirms and agrees with findings in [18] for the general attitude dynamics.

Comparison with Previous Formulations
Given the proof of the mathematical equivalence of the r-E-L and N-E, it is straightforward to state that the presented formulation leads to improved performance with respect to the E-L quadrotor model found in literature.This section shows the difference in performance considering the model employed in section 4.1 and 4.2, as simulation platform, and in section 4.3, for model-based control.

Implementation Comparison between the different models
In this section, it is shown that the error between the r-E-L and N-E formulations is lower compared to the one between the E-L and N-E formulations.To do so, the same rotor input velocity, u = [475.9+ 0.1 sint, 476.2 + 0.1 sint, 476, 476.1] is applied to the (N-E, E-L, r-E-L) quadrotor models for a period of 60s, and the root mean square error (RMSE) is computed for the generalized coordinates (p, η, ṗ, η).The choice of these rotor inputs is such that the quadrotor will cover a long range in x, y, z, while having non constant attitude and not reaching any singular configuration.From Table 1, which illustrates results with an integration step of 10ms, it is evident that the r-E-L has a RMSE of 3 to 6 orders of magnitude smaller than the literature E-L.Moreover, Table 2 shows that by further decreasing the integration step to 1ms, the error of the revised model decreases even more, while this is not true for the E-L model found in literature E-L, leading to an error of 6 to 9 order of magnitude smaller.Note that, even though they are studied independently, the position dynamics are affected by the attitude dynamics due to the influence of the rotation matrix.For this reason it is worth considering the RMSE for p and ṗ. as well.

Comparison with Multibody Dynamic Simulator
Now, the comparison of these mathematical models with respect to a multibody dynamic simulator such as Mathworks's Simacape Multibody is considered.Note that Simscape Multibody uses its own dynamic engine to solve the equation of motions, while only the structure and the inertia of the quadrotor is provided by the user.As shown in Table 3, providing the same input as before to all models, the N-E and the r-E-L state dynamics are almost identical to the one of a multibody dynamic simulator, while the literature E-L model has a much larger error.Again, as shown in Table 4, this becomes even more evident when decreasing the simulation step to 1ms.
It is essential to state that in order to achieve these error magnitudes, all mathematical models (N-E, E-L, r-E-L) need to account for the gyroscopic effect, which is automatically computed in the dynamic simulator.The gyroscopic effect is modeled as in [19], and since it depends on the external input u, for the r-E-L model, it is pre multiplied by W T .

Effect on Model Based Control
By implementing a PID controller with feedback linearization on the multibody dynamic simulator quadrotor the different effects of the E-L and r-E-L are analyzed.The helix trajectory to track is displayed in Fig. 2, along the resulting in the attitude trajectory.The attitude errors resulting from the numerical simulations with different PID gain configurations are shown in Fig. 3.While at some gain values the controller performance is identical, it is shown that by increasing the controller gains, the controller that uses the r-E-L for dynamic compensation has bigger margin of stability.The feedback linearization using the literature E-L model leads to unstable behaviour with smaller gain compared to the r-E-L.

Conclusions
A revised E-L attitude dynamics model of quadrotors / multirotor UAVs has been presented and derived.The equivalence to the N-E formulation has been demonstrated through analytical steps and numerical simulations.Compared to existing E-L formulations found in Building on this work, previously presented feedback linearization controllers using E-L formulation (including the ones from the authors) could be improved.

2. 1
NotationLet I n denote the n × n identity matrix.Moreover, given vectors a, b ∈ R 3 , denote S(a)b = a × b, where

Fig. 1
Fig. 1 Tait-Bryan 321 Sequence of Elementary Rotations from IF to BF

Table 2
Comparison with E-L Model from