Model Formulation Over Lie Groups and Numerical Methods to Simulate the Motion of Gyrostats and Quadrotors

: The present paper recalls a formulation of non-conservative system dynamics through the Lagrange–d’Alembert principle expressed through a generalized Euler–Poincaré form of the system equation on a Lie group. The paper illustrates applications of the generalized Euler–Poincaré equations on the rotation groups to a gyrostat satellite and a quadcopter drone. The numerical solution of the dynamical equations on the rotation groups is tackled via a generalized forward Euler method and an explicit Runge–Kutta integration method tailored to Lie groups.


Introduction
The dynamics of non-conservative systems may be formulated through the Lagrange-d'Alembert principle. This very general principle allows one to deduce the general properties and laws of motion of bodies as well as of non-compressible, compressible and elastic liquids [1]. For continuous systems, the Lagrange-d'Alembert principle is formulated as a variational integral problem, which may be converted into a differential form known as the Euler-Poincaré equation. When the physical problem under investigation is subjected to smooth holonomic constraints, such as in the rigid-body case, the domain of the Euler-Poincaré equation becomes a smooth manifold and, under specific symmetries (i.e., the Lagrangian of the system is left-invariant), even a Lie group. In fact, Lie groups are special instances of differentiable manifolds that also carry on the structure of algebraic groups. In this event, one gets a reduced Euler-Poincaré equation on a Lie algebra paired to a companion differential equation on the associated Lie group. The resulting system of differential equations that describes completely the time-evolution of the system may be regarded as an initial-value problem on a tangent bundle. Such an initial value problem needs to be solved numerically on a computer platform by way of numerical methods that are specifically tailored to upholding the structure of the phase space.
The present paper of archival nature aims at presenting in details, and in a self-contained manner, two interlaced topics, namely: 1) Formulation of Euler-Poincaré equations of motion on Lie groups: This part, covered by Sections 2 and 3, aims at recalling with some mathematical depth the Lagrange-d'Alembert principle expressed through a generalized Euler-Poincaré form. Moreover, this part aims at presenting in a detailed way the derivation of the reduced Lagrangian function for two systems, namely a satellite gyrostat and a quadrotor drone. In addition, this part aims at showing how the general Lagrange-d'Alembert principle boils down to equations of motion through their reduced Euler-Poincaré equations on the special orthogonal group SO(3).
2) Numerical integration of the Euler-Poincaré equations of motion: this part, covered by the Sections 4 and 5, discusses the numerical solution of initial-value problems on smooth manifolds as an extension of two well-known techniques to numerically tackle classical initial-value problems, namely a forward Euler method and an explicit fourth-order Runge-Kutta method.
Section 6 presents some final remarks and concludes the paper.

The Lagrange-d'Alembert Principle and the Forced Euler-Poincaré Equation
We consider a non-conservative dynamical system whose state space is a Lie group G. It is instrumental to recall the following definitions and properties [2,3] (see also [4,5]): Lie group: a smooth manifold G that is also an algebraic group is termed a Lie group. An algebraic group is a structure (G, m, i, e) where m : G × G → G is an associative binary operation termed group multiplication, e ∈ G is an identity element with respect to the multiplication, namely m(g, e) = m(e, g) = g for any g ∈ G, and i : G → G denotes inversion with respect to multiplication, namely m(i(g), g) = m(g, i(g)) = e for any g ∈ G. A left translation L : G × G → G is defined as L g (h) := m(i(g), h). (Our definition is not standard, we adopt it because it simplifies the notation and makes it apparent that left translation and group multiplication are different notions.) In this report, we deal with matrix Lie groups, therefore we may take m(g, h) := gh (matrix multiplication), i(g) := g −1 (matrix inversion) and the identity matrix as identity element. In this setting, L g (h) = g −1 h.
Tangent bundle and metrization: given a state g ∈ G, the tangent space to G at g will be denoted as T g G. The tangent bundle associated with a Lie group G is denoted by TG and plays the role of phase-space for a dynamical system whose state-space is G. The inner product of two tangent vectors ξ, η ∈ T g G is denoted by ξ, η g . A smooth function F : G → G induces a linear map dF g : T g G → T F(g) G termed pushforward. For a matrix Lie group, the pushforward map d(L g ) h : T h G → T g −1 h G associated to a left translation is d(L g ) h (η) := g −1 η.
Lie algebra: the tangent space g := T e G to a Lie group at the identity is termed Lie algebra. The Lie algebra is endowed with Lie brackets, denoted as [·, ·] : g × g → g, and an adjoint endomorphism ad ξ η := [ξ, η]. On a matrix Lie algebra, [ξ, η] := ξη − ηξ. A pushforward map d(L g ) g : T g G → g is denoted as dL g for brevity. Given a smooth function : g → R, for a matrix Lie group one may define the fiber derivative of , ∂ ∂ξ ∈ g, at ξ ∈ g as the unique algebra element such that ∂ ∂ξ , η e = tr (J ξ ) η for any η ∈ g, where J ξ denotes the Jacobian matrix of the function with respect to the matrix ξ.

The Lagrange-d'Alembert Principle
The Lagrange-d'Alembert principle is one of the fundamental principles in mathematical physics to formulate the motion of a physical system on a differentiable manifold and to handle non-conservative forces. Let Λ : TG → R denote a Lagrangian function and F : TG → TG a generalized force field (more precisely, a generalized force field is a smooth map from TG to its dual T G or, for left-invariant force fields, from an algebra g to its dual g ; again, we adopt this non-standard definition because it eases the notation). The Lagrange-d'Alembert principle affirms that a dynamical system generates a trajectory g such that: where the symbol δ denotes variation, which vanish at endpoints and is elsewhere arbitrary and an over-dot (as inġ) denotes (partial) derivative with respect to the parameter t.
We recall that a variational principle consists in taking a continuous family of curves g : U ⊂ R 2 → G, denoted as g(t, ε), where the index ε selects a curve in the family and the index t denotes a point over this curve. All the curves in the family depart from the same initial point and arrive at the same endpoint, namely, g(a, ε) and g(b, ε) are constant with respect to ε. The variations in Equation (1) are defined as δ b a Λ(g,ġ)dt := b a ∂ ∂ε Λ(g(t, ε),ġ(t, ε))dt ε=0 , (2) δg := ∂g(t, ε) ∂ε ε=0 . ( The following lemma affords handling variational principles on Lie groups: Lemma 1 ([6]). Given a smooth function g : U ⊂ R 2 → G on a matrix Lie group, define: where L denotes a left translation and dL its pushforward. Then, it holds that ∂ξ ∂ε =η + ad ξ η.
Assuming that the Lagrangian as well as the generalized force field F are left invariant, we may write Λ(g,ġ) = (dL g (ġ)) and dL g (F(g,ġ)) = f (dL g (ġ)), where : g → R and f : g → g denote a reduced Lagrangian and a reduced force field, respectively. In addition, if the inner product is left-invariant, it holds that F(g,ġ), δg g = f (dL g (ġ)), dL g (δg) e .
Therefore, the integral Lagrange-d'Alembert principle (1) reduces to where it is legitimate to replace dL g (ġ) with ξ and dL g (δg) with η and then set ε to 0.

The Forced Euler-Poincaré Equation
By means of the Lemma 1, the variations in the reduced Lagrange-d'Alembert principle may be worked out to derive equations of motion in differential form, as shown in the following. Theorem 1 ([6]). Let ξ := dL g (ġ) and η := dL g (δg). The solution of the integral Lagrange-d'Alembert Equation (8) under perturbations of the form ∂ξ ∂ε =η + ad ξ η, which vanishes at endpoints, satisfies the Euler-Poincaré equation d dt Proof. For the first term in (8), it holds that where ∂ ∂ξ ∈ g denotes the fiber derivative of the reduced Lagrangian and ξ := dL g (ġ). For the second term in (8), it holds that where η := dL g ∂g ∂ε . Recall that both fields ξ and η depend on t as well as on ε. The Equation (8), thanks to the equivalences (10) and (11) and to the relation Integration by parts yields because the perturbation η vanishes at the endpoints, therefore the Equation (12) or, equivalently, introducing the conjugate ad ξ of the adjoint operator with respect to the metric, b a − d dt Setting ε = 0 and recalling that the perturbation η is arbitrary in the above equality implies the assertion.
The complete dynamical equations then reads The forcing term takes into account a number of external driving phenomena, such as Rayeigh dissipation: this phenomenon is a form of energy dissipation due, e.g., to friction. The forcing term due to Rayleigh dissipation is written as f R = − ∂ρ ∂ξ , where ρ : g → R is a Rayleigh function [6]. For instance, a quadratic function of the state variable ξ leads to a linear dissipation term known as viscous (or aerodynamic) drag.
Other driving causes: other than dissipation, the structure of the forcing term is free and depends on the problem under investigation, as long as it provides a map onto g. This leaves room to incorporate into f control terms, for instance, to stabilize the motion or to drive a dynamical system as discussed, for instance, in [7].
There exist by now quite detailed books, papers, and reviews about these concepts, including the book by Marsden and Ratiu [8] and the two volumes by Holm [9,10], where these notions are explained in great details and with a number of examples going from the single rigid body to more complex multi-body dynamics problems.
Let us elaborate on the notion of dual adjoint operator on the Lie algebra so(3).

Lemma 2.
Let the vector space so(3) be endowed with the canonical inner product ξ, η := −tr(ξη), and define the linear endomorphism ad ξ (η) := ξη − ηξ. The dual of the adjoint operator with respect to the canonical inner product is ad ξ = −ad ξ .
By definition of dual operator with respect to a metric, it must hold that ad ξ (η), ζ = η, ad ξ (ζ) for any η, therefore the assertion follows.
Since, in the state space G := SO (3), it holds that (dL R ) −1 (ξ) = Rξ and ad ξ η = −ad ξ η (by Lemma 2), the Euler-Poincaré equations read where τ is an instance of generalized force and denotes the resultant of all external mechanical torques. In this context, the state variable g := R denotes the attitude of a rigid body and the state-variable ξ denotes its instantaneous angular velocity. It is convenient to define the operator · as: Since any skew-symmetric matrix in so(3) may be written as in Equation (20), it is convenient to define a basis of this vector space as follows: In the present setting, we take the canonical metric ξ, η e := tr(ξ η). With this choice, it not difficult to show that ∂ ∂ξ In order to shorten some relations, it is also convenient to introduce the matrix anti-commutator {A, B} := AB + BA. The anti-commutator is a symmetric form, namely {A, B} = {A, B}, and linear in both arguments. Conversely, the matrix commutator in so (3) is an anti-symmetric bilinear form, namely [ξ, η] + [η, ξ] = 0. (We skip the Jacobi identity that plays no role in this paper.) Moreover, it is convenient to recall a property of the matrix 'trace' operator, namely tr(ABC) = tr(BCA) = tr(CAB) for any square A, B, C.

Satellite Gyrostat
The dynamics of rotating bodies is a classic topic in mechanics. Several aspects of the motion of a rotating rigid body were investigated in the eighteenth and nineteenth centuries, but the study of the dynamics of rotating bodies is still very important for numerous applications as in satellite gyrostat, spacecraft, and robotics [11]. The gyrostats are scientific models or instruments designed to illustrate experimentally the dynamics of a rotating body such as the spinning-top, the hoop and the bicycle, and also the precession of the equinox and the rotation of the earth. The original gyrostat is an instrument designed by Lord Kelvin to illustrate the complicated state of motion of a spinning body when free to wander about on a horizontal plane, like a top spun on the pavement, or a hoop or bicycle on the road [12]. It consists of a massive flywheel concealed in a metal casing. In general, a gyrostat is a system of bodies whose relative motion does not alter the intrinsic mass distribution of the system. A solid gyrostat is an arbitrary rigid body attached to one or more axisymmetric rotors with their axes fixed in the carrier [12].
Let us consider the case of a satellite represented by a gyrostat as in [7,13]. A satellite gyrostat is a suspended rigid body, termed platform, that contains one or more rigid wheels, termed rotors, that are spinning around fixed axes. Since the total angular momentum must stay constant, a rotation of the wheels causes a global counter-rotation of the gyrostat, which makes it possible to stabilize its motion (known as the principle of momentum transfer [13]).
In this paper, the interactions of the translational motion of a gyrostat satellite on its attitude dynamics is completely neglected. Assuming that the center of the gyrostat is fixed in space, one may attach an inertial reference frame F E to the center and denote by x 0 ∈ R 3 the position of each infinitesimal volume of the platform in a platform-fixed reference frame F P , as shown in the Figure 1. At any time t, the position of the infinitesimal volume with respect to F E is x ∈ R 3 and, since the platform is a rigid body, x(t) = R(t)x 0 , where R ∈ SO(3) denotes a rotation matrix that takes the platform-fixed frame F P to coincide with F E . The kinetic energy of the platform P with respect to the inertial reference frame F E may be written as where dm denotes the mass content of the infinitesimal volume. Recalling that ẋ 2 = tr(ẋẋ ), x =Ṙx 0 , and thatṘ = Rξ, with ξ ∈ so(3), we get: where tr(·) denotes matrix trace andĴ p ∈ R 3×3 is defined aŝ and denotes a non-standard inertia tensor [14]. The standard inertia tensor of the platform is defined as These inertia tensors are related as outlined in the following

Lemma 3 ([14]
). The non-standard moment of inertiaĴ of a body is related to its standard moment of inertia J by the relationshipĴ = 1 2 tr(J)I 3 − J.
In addition, assume that the gyrostat contains A rotors spinning on fixed axis inside the platform. Each rotor is indexed by a = 1, 2, . . . , A. The coordinate of the center of each rotor in the platform-fixed reference system F P is denoted as c a ∈ R 3 and the rotation matrix that describes the spinning of each rotor with respect to F P is denoted as R a ∈ SO(3). For a specific rotor R a , the coordinate of each volume element in this rotor, with respect to the reference frame F P , at time t is x(t) = R(t)(c a + R a (t)s), where s ∈ R 3 denotes the position of the volume element in a rotor-fixed reference system. Apparently, it holds thatẋ =Ṙ(c a + R a s) + RṘ a s ⇒ R ẋ = ξc a + (ξR a +Ṙ a )s. (28) Consequently, the kinetic energy per mass element dm of the rotor R a may be written as 1 2 tr(ẋẋ ) = 1 2 tr((R ẋ)(R ẋ) ) = 1 2 tr(ξc a c a ξ ) + 1 2 tr((ξR a +Ṙ a )ss (ξR a +Ṙ a ) ) + tr(ξc a s (ξR a +Ṙ a ) ). (29) The kinetic energy of each spinning rotor R a in the earth frame F E may be written as and straightforward calculations lead to the expression k a = 1 2 M a tr(ξc a c a ξ ) + 1 2 tr((ξR a +Ṙ a )Ĵ a (ξR a +Ṙ a ) ) + M a tr(ξc ac a (ξR a +Ṙ a ) ), where In order to simplify Equation (31), we assume that each rotor is symmetric about its center of mass, which implies thatc a = 0. Moreover, we assume that each rotor may be schematized as a cylinder spinning around its central axis, which implies that R aĴa R a =Ĵ a . The following lemma studies a special case where this result holds true. Lemma 4. Take a cylindrical rotor spinning around the z axis, whose rotation matrix is denoted by R z (t) and its non-standard inertia matrix is denoted byĴ cyl . It holds that R z (t)Ĵ cyl R z (t) =Ĵ cyl .
whereĴ :=Ĵ p + ∑ A a=1 M a c a c a . The spinning of the rotors is assumed to be frictionless. Notice that, upon defining C := [c 1 c 2 · · · c A ], the non-standard inertia tensor may be rewritten aŝ J =Ĵ p + Cdiag(M 1 , M 2 , . . . , M A )C [13].
In the following, we write explicitly the equations of motion of a gyrostat in a number of cases. In all cases, the system is supposed to be subjected to a constant potential energy, which may be assumed to be 0 without loss of generality. The potential energy is non-constant in the case that the gyrostat is subjected, e.g., to a gravitational field, as explained in [15].

Frictionless Gyrostat with Blocked Rotors
This case corresponds to non-spinning rotors. Analytically, this condition corresponds to the case that all ξ a = 0. In addition, no friction is assumed for the rotation of the gyrostat. The reduced Lagrangian of the gyrostat hence reads whereĴ g :=Ĵ + ∑ A a=1Ĵ a denotes the total tensor of inertia of the gyrostat (platform plus non-spinning rotors).
In order to write explicitly the Euler-Poincaré Equation (19), it is necessary to compute the fiber derivative of the reduced Lagrangian: It is not difficult to recognize that the second equation is the classical Euler equation of rotation of an asymmetric rigid body under its own angular momentum unbalance. In fact, taking the platform-fixed reference system F P coincident with the principal axes of inertia and writing it is easy to see that the non-standard tensor of inertia for this configuration iŝ Therefore, by matrix calculations: By equating the homologous entries of the above matrices one gets the classical Euler equations: Notice that for a 'needle' gyrostat (say, J z = 0), the Euler equations degenerate. Furthermore, if the inertia tensor is spherical, namely J x = J y = J z , then the platform spins at constant angular speed around a non-precessing axis.
Besides, by using the notation introduced in this subsection, it is not difficult to show that

Frictionless Gyrostat with Spinning Rotors
We consider now the case of A = 3 rotors positioned along the three axes of inertia of the gyrostat as studied in [7]. This configuration noticeably simplifies the equations and makes it easy to simulate the motion dynamics of the satellite, as well as to design proper control strategies. The rotors along the x (a = 1) and the y (a = 2) axes spin with constant angular velocity, while the rotor along the z axis (a = 3) spins with oscillating angular velocity, namely: where ω 1 , ω 2 , ω 3 > 0, J 11 , J 22 , J 33 > 0, b > 0 and ν > 0 are given constants and the symbol denotes two identical unessential parameters per inertia tensor. The reduced time-dependent Lagrangian reads and its fiber derivatives takes the form whereĴ g :=Ĵ + ∑ 3 a=1Ĵ a . Since everything is constant with respect to the time t except for ξ and ξ 3 , it holds that d dt Let us take again and the non-standard inertia tensor of the gyrostat is as in (41). For a frictionless gyrostat with spinning rotors, the Euler-Poicaré equation reads or, equivalently where we have set The products J aa ω a , for a = 1, 2, 3, are the components of the angular momenta of the three rotors. (For comparison, in [7] these products are denotes as h 1 := J 11 ω 1 , h 2 := J 22 ω 2 and h 3 := J 33 ω 3 .) The following lemma details a case where the formula (56) holds true.

Lemma 5.
Take again a cylindrical rotor spinning around the z axis, as in Lemma 4, with angular speed ω z . It holds that {Ĵ cyl , ω z ξ z } = J zz ω z ξ z .
Proof. The non-standard inertia tensor of the rotor and its orientation matrix take the form Taking the time-derivative of the orientation matrix giveṡ where Computing the matrix products and sum, it is not difficult to prove the assertion.
By equating the homologous entries of the matrices across the two sides of the relation (56), one gets the three differential equations:

Friction Gyrostat with Spinning Rotors
According to [7], we may assume the presence of a significant linear damping along the z axis only. We found that this corresponds to a quadratic Rayleigh damping function where γ > 0 is a given constant friction parameter. Its fiber derivative reads ∂ρ ∂ξ = 1 2 {P, ξ}, hence the generalized forcing term assumes the meaning of a linear viscous drag torque that takes the expression For a friction gyrostat with spinning rotors, the Euler-Poicaré equation reads where all quantities are defined as in the previous Section 3.1.2. By equating the homologous entries of the matrices across the two sides of the relation (66), one gets the system of differential equations: Friction naturally stabilizes the motion of a gyrostat satellite. Stabilization may be achieved by damping effects introduced ad-hoc, such as the internal damping provoked by a cavity filled with viscous fluid [16].

Friction Gyrostat with Spinning Rotors and Feedback Control
According to [7], the equations of the gyrostat may be completed by a state feedback control field 1 2 ϕ : g → g. The complete set of equations of a friction gyrostat with spinning rotors and feedback control reads where all quantities are defined as in the previous subsubsections. The control field suggested in [7] reads where ω r denotes a desired gyrostat spinning velocity and κ 1 , . . . , κ 6 are given control constants. Further terms may be added to the right-hand side of the Equation (70) to control chaos or to achieve synchronization of two gyrostats.

Explicit State-Space Form of Gyrostat Equations
On the basis of the expression (70), we want to express the Euler-Poincaré equation in explicit form, namely, as a first-order abstract systemξ = σ(ξ), with σ ∈ End(so (3)). To this aim, let us define It is straightforward to verify that {Ĵ g ,ξ} = DξD, and that det(D) = J x J y J z , hence the constant matrix D is invertible and the expression (70) may be rewritten in state-space form aṡ This way of expressing the equations of motion, that is equivalent, in principle, to the usual way based on angular coordinates and momenta, is much more compact. We mention that it is also possible to make use of a set of variables introduced by Deprit (see, e.g., [15]), that makes the equations of motion take much simpler forms than the usual ones in terms of the Euler angles and the corresponding angular momenta.
Indeed, the satellite gyrostat equation may be rewritten as This is a first-order system of differential equations on the tangent bundle TSO(3) that may be solved (numerically) upon establishing two initial conditions (typically ξ(0) = 0 and R(0) = I 3 ). The system (74) is also equivalent to the second-order differential equation Recalling that the covariant derivative ∇ṘṘ of the tangent velocity fieldṘ with respect to itself is given byR −ṘR Ṙ , the dynamics in Equation (75) may be rewritten, equivalently, as As a special case, under the hypotheses that (a) the gyrostat is perfectly balanced (its intertia tensor is spherical), (b) the rotors are blocked, (c) the motion is frictionless, and (d) there are no external forcing terms, then σ ≡ 0 and the motion is described by ∇ṘṘ = 0, namely, the gyrostat moves along a geodesic in SO(3), which represents a uniform rotation.

Quadrotor Drone
A quadrotor drone (also referred to as quadcopter [17][18][19]) is a helicopter with four rotors directed upwards and placed in a square formation with equal distance from the center, as in Figure 2. The quadrotor is controlled by adjusting the angular velocities of the rotors, which are responsible for both lift and propulsion. A horizontal movement is achieved by tilting the platform, while a vertical movement is achieved by changing the total thrust of the rotors [20]. In the simplest model, as the one studied here, the only degrees of freedom are the velocity at which the four rotors spin. An unbalance in the spinning velocities of the propellers causes, in general, the platform to tilt (hence to change its attitude and to shift horizontally), to ascend or descend (that is, to shift vertically), and to spin along its vertical axis (which is referred to as yawing).
In general, a quadrotor model is written on the basis of the following assumptions: The structure is supposed to be rigid and axis-symmetrical, with the center of the axis coinciding with the center of mass, thrust and drag are proportional to the square of each rotor's angular speed, and the rotors are supposed to be rigid. Further assumptions may be made to simplify the analytic treatment.
The quadrotor is made of a body B and of four rotors R a , with a = 1, 2, 3, 4, also termed propellers. We assume that one pair of propellers (a = 2 and a = 4) is rotating clockwise, while the other pair of propellers (a = 1 and a = 3) is rotating counterclockwise, as shown in Figure 3. The rotor R 1 is located along the +x axis (hence the rotor R 3 is located along the −x axis), while the rotor R 4 is located along the +y axis (hence the rotor R 2 is located along the −y axis). The coordinates of each volume element of the drone are referred to as an earth reference frame denoted by F E . A reference frame, denoted as F B , is attached to the body of the quadrotor, with the origin in its center of mass, as shown in Figure 3. Figure 3. Quadcopter reference frames.
In the present subsection, as already done for the gyrostat satellite, we proceed on determining a model of quadcopter by a fully geometric approach, namely by means of the Euler-Poincaré equations written in compact form, without resorting to component-wise expressions or classical physics definitions. This procedure makes the system modeling more straightforward and, we believe, notationally easier to manage (compare, e.g., to a more classical procedure used in [21]).

Total Kinetic Energy
The position of the center of mass of the quadrotor with respect to an inertial reference frame F E is denoted by q ∈ R 3 . The orientation of a quadrotor-fixed reference frame F B with respect to the inertial reference frame F E is denoted by a rotation matrix R ∈ SO(3). The matrix R represents the rotation needed to align the orthogonal frame F B to the orthogonal frame F E . A rotation around the x axis is termed roll, a rotation around the y axis is termed pitch and a rotation around the z axis is termed yaw. The matrix R changes over time according to the tilting of the quadrotor. When R = I 3 , the quadcopter is perfectly horizontal and the reference frames are perfectly aligned.
Any volume element of mass dm within the body of the quadcopter is individuated by a coordinate x = q + Rx 0 with respect to the earth reference frame F E , where x 0 denotes the position of the volume element in the body-fixed reference frame F B . The kinetic energy of the body is given by whereẋ =q +Ṙx 0 . The vectorq denotes the linear velocity of the body center in the earth-fixed reference frame F E . Recalling thatṘ = Rξ, it is not difficult to obtain the following expression for the kinetic energy of the body: where The scalar quantity M b denotes the total body mass and the 3 × 3 matrixĴ b denotes the non-standard inertia tensor of the quadcopter body. The vectorc b ∈ R 3 denotes the position of the center of mass of the body in the body-fixed frame. By the assumptions made, it holds thatc b = 0.
Next, we embody in the model the information that the four rotors are located at equal distance from the center of the body, two on the x axis and two on the y axis. A point on the rotor R a has coordinates, in the inertial reference frame, given by where q,q and R are defined as before, p a ∈ R 3 denotes the position of the rotor R a with respect to the body-fixed frame, and R a ∈ SO(3) denotes the instantaneous rotation of a rotor-fixed frame with respect to the body-fixed frame F B ; moreover,Ṙ a = ξ a R a (recall that dp a dt = 0). Then, the kinetic energy of the rotor R a is given by where The last two equalities are due to symmetry of each rotor and to their symmetric disposition within the quadcopter body.
The total kinetic energy of the quadcopter is, thus, given by: where Notice that the rotors have been assumed symmetric with respect to the spin axis z, hence we exploited again the relation R aĴa R a =Ĵ a .

Equations for the Rotational Component of Motion
Since all the rotors are spinning about the z axis, ξ a := (−1) a ω a ξ z . In addition, we assume that all rotors have the same inertia tensor, therefore {Ĵ a , ξ a } = (−1) a ω a J R ξ z , where J R is a given constant. The factor (−1) a was introduced in the equations to take into account the fact that rotors R 1 and R 3 are spinning counterclockwise, while rotors R 2 and R 4 are spinning clockwise. In this manner, their angular velocities ω 1 , ω 2 , ω 3 , ω 4 represent absolute speeds proportional to their RPM (revolutions per minute) rate.
In order to obtain equations that describe the rotational component of motion, we take (ξ) := k(ξ) and compute the fiber derivative: Each rotor exerts a thrust ϕ a ∈ R 3 along the z axis of the body-fixed reference frame (that is, a vertical force in the body-fixed frame F B that tends to push the quadrotor body upward), as shown in the Figure 3. The magnitude of the thrust generated by rotor R a is proportional to ω 2 a via a lift constant denoted by 1 2 b > 0, namely: The lift constant is proportional to the air density and to the rotor cross-section. Each thrust ϕ a is supposed to appear at the tip of each rotor, therefore its moment arm p a produces a mechanical torque on the body of the quadcopter, which has the effect of tilting the quadrotor drone with respect to the earth frame F E during flight. The four moment arms are given by Moreover, each rotor produces a drag torque, which is directed along the z axis and is proportional to the square angular velocity of each rotor through a drag constant 1 2 γ > 0, namely τ a := γ(−1) a ω 2 a ξ z . A drag torque is present due to the spinning of each rotor which, by virtue of the angular momentum conservation law, will produce a counter torque that acts on the whole body of the quadcopter and tends to produce a rotation along its vertical axis.
Therefore, the total external torque acting on the quadcopter is The external mechanical torque τ ∈ so(3) may be written in components with respect to the basis (21) as τ x = ξ z + τ y ξ y + τ z ξ z . By direct computations, it is found that Therefore, the torque due to rotor thrusts has the effect of tilting the quadcopter body along the roll and the pitch axes, while the torque due to drag produces a rotation of the body around the yaw axis. It is worth remarking, taking as the reference Figure 3, that a positive rotation along the yaw axis coincides with the direction of rotation of rotors R 1 and R 3 and is opposite to the direction of rotation of rotors R 2 and R 4 . The expression τ z = γ(ω 2 2 + ω 2 4 ) − γ(ω 2 1 + ω 2 3 ) tells that the drag of rotors R 1 and R 3 tends to rotate the body of the quadcopter in a negative yaw direction, while the drag of rotors R 2 and R 4 tends to rotate the body of the quadcopter in a positive yaw direction. This sign inversion is due to the fact that, when an internal part tends to rotate in a given direction, the whole body tends to rotate in the opposite direction, because of the conservation of the total angular momentum.
In summary, the Euler-Poincaré equations (19) for the quadrotor drone read The second equation (that expresses the dynamics of ξ) is much simpler than the equations of motion derived by using local coordinates (e.g., yaw, pitch and roll angles) because it is written with respect to a body-fixed frame F B and expressed in the Lie algebra so(3) (compare, e.g., with [20]). For completeness and easiness of comparison to other formulations, we report below the differential equations governing the evolution of the components ω x , ω y , ω z of the angular tensor ξ: where diag(J x , J y , J z ) = J q . Some comments on these equations are in order: Rotors spinning velocities: the four rotors spin at almost constant velocity and do not change the spinning direction during operation. Their spinning velocity changes only slightly, namely to ω a ± ∆ω a , in order to create a slight unbalance in the thrusts that produces a slight tilting and affords controlling the attitude/tilt of the body. The quantity Ω r := ω 1 − ω 2 + ω 3 − ω 4 is termed residual rotoric angular velocity and the quantity J RΩr is termed inertial counter-torque [22].
Independence from attitude: the equation {Ĵ q ,ξ} = [Ĵ q , ξ 2 ] + [β, ξ] −β + τ governs the body rotation velocity (in terms of rotation speed and direction). This equation is independent of the current attitude R of the body itself. In fact, the matrix R does not appear in this equation.
In steady state flight (i.e., not takeoff or landing), the term J RΩr is unessential because Ω r is kept constant, since most of the time the rotors will be maintaining an almost constant thrust and will not be accelerating. Moreover, due to the symmetry of the quadrotor, one may take J y ≡ J x . With these assumptions, the model simplifies tȯ The global spinning motion around the z axis is due only to drag. Since the propellers R 2 and R 4 are rotating clockwise, while the propellers R 1 and R 3 are rotating counterclockwise, the Equation (95) tells us that the net difference between the clockwise-directed angular velocities and the counterclockwise-directed angular velocities produces a counterclockwise spinning of the body of the drone.

Equations for the Translational Component of Motion
The translation component of motion obeys the Newton's law written in the inertial (earth) reference frame F E , namely, the linear accelerationq is proportional to the total thrust ϕ T := ∑ 4 a=1 ϕ a , to the gravity acceleration and to an aerodynamic drag term. Notice that the vector ϕ T is directed along the constant direction e z in the body-fixed reference frame F B which, in turn, is rotated of a quantity R with respect to the earth frame F E , therefore, in the earth-fixed frame, the total thrust is directed along the rotated axis R · e z . Formally, the just-recalled Newton's law is written as: wherep denotes the gravitational acceleration scalar and Γ ∈ R 3×3 denotes an aerodynamic drag tensor (typically, diagonal). Whenever the quadcopter is in a horizontal position, namely R = I 3 , its center of mass can only move vertically (the thrust is directed along e z ). Whenever the quadcopter is tilted, namely Re z − (e z Re z )e z = 0, then the quadcopter changes direction.

Explicit State-Space Form of Quadcopter Equations
By adopting again the concept introduced in the Section 3.1.5, that is {Ĵ q ,ξ} = DξD, where the matrix D is defined as in (72), we may rewrite the set of equations governing the flight of a quadcopter drone as For exemplary values of the parameters (OS4 Mini-VTOL quadrotor model), interested readers might want to consult [22]. The model (95) does not include a number of physical effects, such as hub moment and rolling moment due to sideward/forward flight, nor a sub-model of the brushless direct-current motor and of the gearbox.

Numerical Integration of Initial Value Problems
We focus on dynamical systems of the first order, namely initial value problems [23], since it is always possible to reformulate a dynamical system of higher order into a dynamical system of first order by augmenting the number of variables. We start off by recalling two numerical methods for initial value problems on R n . The reason for this choice is twofold: the methods on R n may be generalized to numerical methods on manifolds and Lie groups by appropriate analysis techniques [24], for example to integrate the equations that govern the rotational component of motion; in addition, the equations that govern the translation component of motion need numerical methods on R 3 .

Numerical Schemes in R n
Classical initial-value problems have the Euclidean space R n as domain [25,26]. Well-known numerical integration techniques used in engineering and in applied sciences are the forward Euler method and the explicit fourth-order Runge-Kutta method. The original Euler's method [27] exhibits a good behavior as a numerical integration method and can be used to study a number of important phenomena [28]. The Euler method is traditionally the first numerical technique introduced to solve initial value problems. It is simple to understand and geometrically easy to articulate, but it is known to show limited accuracy for complicated functions. A more robust and intricate family of numerical integration method is the Runge-Kutta scheme. This method is the most widely used one since it gives reliable results and is particularly suitable when the computation of higher derivatives is unpractical. An account for the whole spectrum of both stiff and non-stiff ordinary differential equations solvers is provided by [29,30]. A classical account of Runge-Kutta methods was given in the monograph [31] and in the textbooks on ordinary differential equations [28,32,33].
In addition to Runge-Kutta and linear multi-step methods, Taylor series methods may be considered as a third main direction of development. These three main classes of numerical ordinary differential equations solvers have, respectively, the characteristic of using more derivatives. A variety of new numerical schemes can be constructed which are combinations of methods from the three main classes. An important class of such new schemes consists of the multi-stage multi-step methods known as multi-step Runge-Kutta methods [32]. Another class of numerical schemes is formed by the multi-step-multi-derivative methods, known as Obreshkov methods [34]. The Kutta-Merson method [35] belongs to the family of so-called embedded Runge-Kutta methods. These schemes are useful because they provide a tool by which the accuracy of the calculations can be monitored during the integration process, but they are also quite burdensome to implement on a computing platform.
A generic first-order, non-autonomous, dynamical system whose state variable is denoted by x ∈ R n may be written asẋ where the function σ : R n × R → R n describes completely the dynamical system. In calculus, the above system is referred to as an "initial value problem" which may be solved under mild regularity conditions. The variable x ∈ R n represents the state of the dynamical system, while its derivative with respect to the time,ẋ, denotes the state speed. We notice here that, incidentally, while x andẋ possess a rather different physical meaning (for example, position and speed, respectively, of a particle) they belong to the same space R n . This is, indeed, rather a special case, since systems whose state evolve on curved spaces do not enjoy this property. We also observe that the explicit dependence of the state-transition function σ from the parameter t (which, in most applications, represent time), serves to summarize the possible dependence of σ from a number of time-varying external sources of information (for example, a linear input-state system may be represented by a differential equation like (98) by setting σ(x, t) := Ax + Bu(t), where u(t) denotes a time-varying input signal). In most of the cases of interest in applied sciences, problem (98) is not tractable analytically and needs to be solved numerically.
In the following subsubsections, we recall from the scientific literature two of the above-mentioned numerical integration schemes, namely the forward Euler scheme and the explicit fourth-order Runge-Kutta scheme and illustrate their behavior on the numerical simulation of a non-linear, chaotic system. Although both methods are most likely well-known to the readers, we believe that it is important to briefly recall their essential features in view of their extension from the space R n to curved manifolds, which will be pursued in the next subsection.

Forward Euler Scheme
The forward Euler scheme [36] (hereafter denoted by fEul), is perhaps the simplest (yet effective) numerical integration scheme for an initial value problem known in the scientific literature. It is based on two key ideas, namely, the uniform time-discretization of the state variable and the approximation of its first-order derivative by means of the right-side incremental ratio.
Uniform time-discretization is a procedure that affords switching from a continuous function of a temporal variable to a succession. In other terms, since we are unable to evaluate the state x(t) at any time, we define a uniform succession of time-steps, separated by a reasonably short interval T > 0 termed stepsize, which are denoted by t 0 = 0, t 1 = t 0 + T, t 2 = t 1 + T, and so forth, at which we will try to approximate the values of the state. The fEul scheme is based on the approximatioṅ where ε ∈ R n denotes an error term. Plugging the above approximation into the equation of the dynamical system (98) and applying a uniform time-discretization, we obtain: The fEul numerical integration scheme arises by ignoring the error term and by replacing the actual values of the state by their approximations, which are denoted by x 1 , x 2 , x 3 and so forth. The magnitude ε(t k ) of the error at every step determines the precision of the numerical scheme. The approximations are calculated iteratively by the numerical scheme The error in this method is as small as T 2 . The number of steps depends on the discretization step-size and on the extension of the time-interval that a solution is sought within. Namely, if t ∈ [t 0 , t f ], then the number of steps is the nearest integer to the ratio (t f − t 0 )/T.

Explicit Fourth-Order Runge-Kutta Scheme on R n
The family of Runge-Kutta numerical integration schemes [36] was developed to increase the precision of the solution and is based on the evaluation of the function f at a number of locations in the interval [t k , t k+1 ]. In particular, the explicit fourth-order Runge-Kutta scheme (hereafter denoted as eRK4), is based on four partial increments that, added together, bring to a complete step. The eRK4 approximation of the actual solution to the initial value problem (98) may be expressed by the following equations: for k = 0, 1, 2, 3, . . .. Once again, by ignoring the error term and by denoting the approximations of the state at step k by x k , we get the eRK4 scheme: Notice that x k , H 1,k , H 2,k , H 3,k , H 4,k , x k+1 all belong to the Euclidean space R n . This is not going to be the case when we will deal with curved manifolds rather than a Euclidean space R n . With reference to the interval [t k , t k+1 ], the eRK4 method advances the solution from the step k to the step k + 1 by combining together four estimations of the vector field σ at four specific points, namely, once at the left-hand side of the interval (H 1,k ), twice at the middle of the interval (H 2,k and H 3,k ), and once at the right-hand side of the interval (H 4,k ). The stepping direction is a linear combination of these partial estimates whose coefficients are carefully crafted so as to ensure the error to be as small as T 5 .

Extension of fEul and eRK4 Schemes to Smooth Manifolds
A dynamical system that has a smooth manifold G (which includes Lie groups) as state space is described byġ where σ : G × R → TG. In fact, notice that, by definition of tangent space, it holds thatġ ∈ T g G, therefore σ(g, t) ∈ T g G. This circumstance makes a significant difference in the numerical integration of problem (104) because, now, the system state and its first-order derivative no longer belong to the same space and, in addition, the state-space is no longer flat, which makes most of the operations listed in schemes (101) and (103) invalid. The extension of the integration schemes to smooth manifolds is based on two key concept from differential geometry, namely, manifold retraction and parallel transport, which are surveyed briefly to clarify the notation [2,3] (see also [4,5] for a non-mathematical standpoint): Manifold retraction: a map Ret : TG → G is an operator that maps a pair (g, ξ) ∈ TG into a point in the manifold G. The notation used within the present manuscript is h = Ret g (ξ) which indicates that it takes a point g ∈ G, a tangent vector ξ ∈ T g G and returns a point h ∈ G. A manifold retraction is a generalization of a manifold exponential map and, as such, it shares zeroth-and first-order features, namely, Ret g (0) = g and d dt Ret g (tξ) t=0 = ξ. In the case of a flat space R n endowed with a Euclidean metric, the geodesic connecting any two points is a straight line and the retraction is given by Ret g (ξ) = g + ξ.
Parallel transport operator: a map Par : TG × G → TG maps a tangent vector ξ at a given point g on a manifold into a tangent vector at another given point h on the same manifold. The notation used within the present manuscript is η = Par g h (ξ) where (g, ξ) ∈ TG and (h, η) ∈ TG. Parallel transport is a generalization of rigid translation of free vectors in R n and Par g h is a linear isomorphism between T g G and T h G as well as an isometry. The above notation tacitly subsumes that the manifold G is endowed with a connection, hence geodesics, in such a way the operator Par transports the vector ξ along the geodesic arc connecting points g, h to give the new tangent vector η. In the case of a flat space R n endowed with a Euclidean metric, parallel transport is given simply by η = ξ.

Forward Euler Scheme on a Smooth Manifold G
As a first step, we replace the continuous trajectory t → g(t) ∈ G with the sequence k → g k ∈ G and denote, for simplicity, ξ k := σ(g k , t k ).
It would be very convenient to approximate the solution of the dynamical system (104) by means of a fEul scheme like x k+1 = x k + Tσ(x k , t k ). However, even if we consider G as embedded in a larger Euclidean space, we quickly recognize that, on a curved manifold, the sum operation is not necessarily legitimate, namely, even though g k ∈ G and ξ k ∈ T g k G, their linear combination g k+1 = g k + Tξ k would result in a point laying outside of the manifold (g k+1 / ∈ G). A way to extend a forward Euler stepping method is to replace addition with retraction, which leads to the following forward Euler scheme on manifold (hereafter referred to as fEulM): (Tσ(g k , t k )), k = 0, 1, 2, 3, . . . , K. (105) The above iterative scheme computes the first K points of the trajectory of the dynamical system (104) spaced apart of T in time.

Explicit Fourth-Order Runge-Kutta Scheme on a Smooth Manifold G
Likewise, our extension of the explicit fourth-order Runge-Kutta method eRK4 to a manifold is based on the replacement of linear operations by retraction and parallel transport to give an explicit fourth-order Runge-Kutta method on manifold (hereafter denoted as eRK4M).
As a first step, we recall from scheme (103) that it is necessary to define four interlaced stages ξ 1,k , . . . , ξ 4,k to advance the solution from the point t k to the point t k+1 . The four stages are outlines as follows: • The first stage of the eRK4 scheme is defined as H 1,k := Tσ(x k , t k ). We can extend this calculation directly to a smooth manifold G to get the first stage of the eRK4M scheme by ξ 1,k := Tσ(g k , t k ). Notice that ξ 1,k ∈ T g k G. The tangent vector ξ 1,k represents an estimation of the tangent vector field σ at t k .

•
The second stage of the eRK4 scheme is defined as H 2,k := Tσ(x k + 1 2 H 1,k , t k + 1 2 T). As we already know, the inner sum x k + 1 2 H 1,k is not necessarily legitimate on a manifold but may be extended by means of a retraction operator as Ret g k ( 1 2 ξ 1,k ). Therefore, we define the second stage of a eRK4M scheme as ξ 2,k := Tσ(Ret g k ( 1 2 ξ 1,k ), t k + 1 2 T).
The tangent vector ξ 2,k represents an estimation of the tangent vector field σ at the midpoint t k + T 2 .

•
The third stage of the eRK4 scheme is defined as H 3,k := Tσ(x k + 1 2 H 2,k , t k + 1 2 T). Now, this definition cannot be made legitimate just by replacing the inner sum with a retraction because ξ 2,k / ∈ T g k G. Before applying a retraction, it is necessary to take a preliminary step, namely, to transport the local estimation of the vector field σ to the space T g k G. For the sake of notation conciseness, let us define the intermediate point (the notation k + 1/4 for the subscript is standard in numerical calculus). Notice now that ξ 2,k ∈ T g k+1/4 G. The third stage in the eRK4M scheme is defined as ξ 3,k := Tσ(Ret g k (Par The tangent vector ξ 3,k represents a further estimation of the tangent vector field σ at the midpoint t k + T 2 , as required by the original eRK4 scheme.

•
The fourth stage of the eRK4 scheme is defined as H 4,k := Tσ(x k + h 3,k , t k + T). Again, before applying a retraction, it is necessary to transport the local estimation of the vector field σ to the space T g k G. Let us define the intermediate point g k+1/2 := Ret g k (Par and notice that ξ 3,k ∈ g k+1/2 . The fourth stage in the eRK4M scheme is defined as ξ 4,k := Tσ(Ret g k (Par • In the original eRK4 scheme, the four stages H 1,k , . . . , H 4,k are combined together as ). We already know that, on a curved manifold, the four stages belong to different tangent spaces and cannot, therefore, be combined together as they stand. In particular, the second increment needs to be transported to the tangent space T g k G by means of parallel transport to giveξ 2,k := Par g k+1/4 g k (ξ 2,k ), the third stage needs to be replaced byξ 3,k := Par g k+1/2 g k (ξ 3,k ), and the fourth stage byξ 4,k := Par g k+3/4 g k (ξ 4,k ), where g k+3/4 := Ret g k (Par g k+1/2 g k (ξ 3,k )). The four stages ξ 1,k ,ξ 2,k ,ξ 3,k ,ξ 4,k ∈ T g k G may be combined together with the retraction map to give the next point in the sequence as g k+1 = Ret g k ξ 1,k 6 +ξ 2,k +ξ 3,k 3 +ξ 4,k 6 . (111) In summary, the above stages and stepping equations may be recapitulated as: with k = 0, 1, 2, 3, . . . , K. In the second, third and fourth equation, we encounter a compound function made of the parallel transport operator and of the vector field σ. Notation-wise (if not computation-wise), it is convenient to compose the two operators Par and σ and to define the new operator Parσ g h (t) := TPar g h (σ(g, t)).
Notice that the new operator Parσ encapsulates the step-size T. In summary, the eRK4M method may be expressed in the compact form: with k = 0, 1, 2, 3, . . . , K. A further observation is that the intermediate points g k+1/4 , g k+1/2 and g k+3/4 serve only to compute the tangent-vector stages, therefore it could be convenient to compactify even further the notation by introducing the operator ParσRet g (ξ, t) := TPar Through this compound operator, by adapting slightly the definitions of the ξ i,k to accommodate the new notation, the eRK4M method may be expressed as: with k = 0, 1, 2, 3, . . . , K and g 0 ∈ G given as initial point. It is worth noticing here that the initial value problem might be more complex than Equation (104), as the vector field σ might depend on other variables whose estimation, at midpoints, would then be necessary. This is the case, in fact, of rigid-body equations, which is tackled specifically in Section 5.

Numerical Integration of Rigid-Body Equations
The present section shows application of the numerical recipes to integrate first-order dynamical systems developed in the Section 4 to integrate the rigid-body equations recalled in the Section 3. In particular, we consider the numerical integration of a gyrostat satellite equations as well as of a quadrotor drone equations.
In order to put into effect the numerical schemes of Section 4.2, we choose as metric for the rotation group the canonical inner product, which leads to the following expressions for a retraction map (taken as the manifold exponential map arising from the Levi-Civita connection) and the parallel transport operator [37]: with R, Q ∈ SO(3), ξ ∈ T R SO(3). The symbol 'exp' denotes matrix exponential, and √ · denotes matrix square root. (In a Matlab environment, for example, the matrix exponential may be evaluated through the function expm and the matrix square root through the function sqrtm.) There is by now an extensive literature on the geometric discretization of equations of Euler-Poincaré type. A non-exhaustive list of such contributions is [38][39][40][41][42][43].

Numerical Integration of A Gyrostat Satellite Equations on SO(3) × so(3)
Let us summarize the equations of a gyrostat satellite, from the Section 3.1.5, as In the following numerical experiments, we take t ∈ [0, 600] (min) and initial conditions ξ(0) = 0 and R(0) = I 3 .
The differential equationξ = σ ξ (ξ, t) has support in so (3), which is a vector space isomorphic to R 3 , hence it may be solved numerically by one of the numerical schemes outlined in the Section 4.1.
The differential equationṘ = σ R (R, ξ) has support in SO(3), which is a curved space, hence it needs to be solved numerically by one of the numerical schemes outlined in the Section 4.2.
The numerical values of the parameters are taken from [7], namely: where measurement units have been omitted. The amplitude b of the forcing oscillatory torque along the z axis varies. Moreover, the values of the controller constants κ 1 = . . . = κ 6 may be varied to emphasize or de-emphasize the role of the attitude stabilizer. We chose to employ an explicit Runge-Kutta method on manifold to solve Equation (119). The complete set of equations to simulate the motion of a gyrostat satellite read: where k = 0, 1, 2, . . . , K, ξ 0 = 0 and R 0 = I 3 . The first five equations represent an eRK4 method to advance the solution about the angular velocity matrix ξ, while the following five equations represent an eRK4M method to advance the solution about the attitude variable R. It is worth noticing that the intermediate values ξ k+1/4 , ξ k+1/2 and ξ k+3/4 calculated by the first three equations are reused in the seventh, eighth and ninth equations because the vector field σ R depends on ξ.
The Figure 4a,b show the results obtained by setting the forcing torque parameter b to 0.5. Also, in this simulation we chose κ 1 = . . . = κ 6 = 1. The Figure 4c,d show the results obtained when the forcing torque parameter b was set to 4.5 and again κ 1 = . . . = κ 6 = 1. In both experiments, it was set T = 0.01 and K = 60,000, in such a way that the simulation spans 10 hours of motion. The trajectories look almost periodic, although perfect periodicity was never observed in the model-based numerical simulations. This particular result shows a slightly chaotic behavior in the motion of the modeled gyrostat satellite.
It is interesting to confirm numerically the role of the control torque ϕ in stabilizing the motion of the modeled gyrostat satellite. Figure 5 shows the results obtained by setting the forcing torque parameter b to 0.5 and by zeroing the controller gains κ 1 = . . . = κ 6 . The simulated trajectory looks more chaotic compared to the previous case where the stabilizing control was switched on.   Quadcopter modeling and simulation is a difficult and challenging problem. With six degrees of freedom (three translational and three rotational) and only four independent inputs (rotor speeds), quadcopters are severely underactuated. In order to achieve six degrees of freedom, rotational and translational motion are coupled, which result in coupled differential equations. Unlike ground vehicles, copters have very little friction to prevent their motion, so they must provide their own damping in order to stop moving and remain stable.
In order to implement a numerical simulation of the dynamics of a quadcopter, we recall from the Section 3.2.4 that the model includes a number of physical constants, four input variables, namely the spinning velocities of the four rotors ω a (t) > 0, a = 1, 2, 3, 4, two state variables, ξ(t) ∈ so (3) and v(t) ∈ R 3 , and two output variables, R(t) ∈ SO(3) and q(t) ∈ R 3 . Let us summarize the equations of a quadrotor drone from the Section 3.2.2 (equations for the rotational component of motions) and Section 3.2.3 (equations for the translational component of motion) as for t ∈ [0, 3] (sec). In order to simplify the numerical simulation, we will assume zero residual rotor angular velocity, namely Ω r = 0, in such a way that β = 0 andβ = 0. Therefore, the state transition function σ ξ simplifies to D −1 ([Ĵ q , ξ 2 ] + τ)D −1 . The parameter J R then disappears from the equations.
This differential system contains three equations on vector spaces (the equation for ξ has support so(3) and the equations for q and v have support R 3 ) and an equation on manifold (the equation for R has support SO (3)). The numerical values of the parameters are taken from [22], namely: where measurement units are in the International System and have been omitted. We will also consider an isotropic aerodynamic drag, namely Γ = 1 4 I 3 . Several other aerodynamic effects could be included in the model as, for example, dependence of thrust on angle of attack, blade flapping and airflow disruption [44]. Since some of the aerodynamic effects have significant impact only at high velocities, these were excluded from the present model.
Even though the mathematical model deals with angular velocities ω a as input variables, the input values are chosen to be the RPM n a of the rotors, because this is the more natural way of setting the propeller rotation velocity. The relation between angular velocity and rotation per minute is ω a = π 30 n a . We repeated the four test experiments suggested in ( [20], Table 1) adapted to the present model. The rotors speed was adapted to the present case by means of a steady state design. Namely, we considered the case that the thrust produced by four rotors spinning at equal speed ω ss balances perfectly the gravitational force, formally: The resulting steady-state rotor speed in RPM (rounded to the nearest integer) reads When the quadcopter is tilted, namely R = I 3 , the total thrust only approximately balances the gravitational acceleration, hence the quadcopter is subjected to a vertical lift. Moreover, upon tilting, the quadcopter might be subjected to a horizontal force, whose combined effect with the vertical lift causes its change of direction. An estimation of the net force, in the earth frame, under a common operation condition, is outlined in the following Lemma 6. Define the net force acting on the quadcopter's center of mass as Assume that ω a = ω ss + ∆ω a , with ∑ a (−1) a ∆ω a = 0 and ∑ a ∆ω a = 0. In addition, assume that (∆ω 1 ) 2 + (∆ω 2 ) 2 2ω 2 ss . Then it follows that, to first order, Proof. By the binomial theorem, we get ω 2 a = ω 2 ss + (∆ω a ) 2 + 2ω ss ∆ω a . Summing over the index a we then get ∑ a ω 2 a = 4ω 2 From the conditions ∑ a (−1) a ∆ω a = 0 and ∑ a ∆ω a = 0 it follows that ∆ω 4 = −∆ω 2 and ∆ω 3 = −∆ω 1 , therefore ∑ a (∆ω a ) 2 = 2((∆ω 1 ) 2 + (∆ω 2 ) 2 ) and thus ∑ a ω 2 a = 4ω 2 ss + 2((∆ω 1 ) 2 + (∆ω 2 ) 2 ).
The net force hence reads By way of the assumption that (∆ω 1 ) 2 + (∆ω 2 ) 2 2ω 2 ss , we can neglect the second and third addendum in the parentheses of the relation (130) and write From the definition (124) it follows the assertion.
The condition ∑ a (−1) a ∆ω a = 0 ensures that the residual Ω r stays null. The hypothesis that ∑ a ∆ω a = 0 is a key assumption to simplify the analysis (and to design a controller). The assumption (∆ω 1 ) 2 + (∆ω 2 ) 2 2ω 2 ss reflects the fact that, to produce a manageable unbalance in the thrusts, the RPM of the rotors need to change just slightly. Not surprisingly, the net force f net is proportional to the tilt R − I 3 . The vertical lift, that causes a quadcopter to change its altitude, is estimated by default. It is also important to estimate the yaw effect due to rotors drag under the same conditions of the Lemma 6. An exact estimation is outlined in the following Lemma 7. Define the mechanical torque responsible for the yaw of the drone as in (90) as Under the same hypotheses of Lemma 6, it holds that Proof. By the binomial theorem, we get The first term on the right-hand side is null. From the conditions ∑ a (−1) a ∆ω a = 0 and ∑ a ∆ω a = 0, it follows that the third term on the right-hand side is null and that (∆ω 4 ) 2 = (∆ω 2 ) 2 and (∆ω 3 ) 2 = (∆ω 1 ) 2 . Therefore, ∑ a (−1) a (∆ω a ) 2 = 2((∆ω 2 ) 2 − (∆ω 1 ) 2 ) and thus the assertion.
The rotors RPM used in the experiments, adapted from [20], are outlined in Table 1. Notice that, in all four experiments, the residual angular velocity is null, namely Ω r = 0, and that the angular deviations ∆ω a satisfy all the assumptions made in the Lemma 6. Furthermore, since all |∆ω a | are either 0 or π 30 , and since the drag constant γ is really small in this quadcopter model, according to Lemma 7, we expect the yaw component of the external torque τ z to be practically negligible in all four numerical experiments. We chose to integrate the differential equations that describe the rotational component of motion as well as the translational component of motion by means of a forward Euler scheme. The complete set of equations read: where k = 0, 1, 2, . . . , K, ξ 0 = 0, R 0 = I 3 , q 0 = 3e z and v 0 = 0. In all the following four experiments, it was chosen K = 3000 and T = 0.001 (the experiments will predict the first 3 seconds of dynamics).
In Experiment 1, the rotors RPM is fixed to n ss . In this experiment, the thrusts are perfectly balanced, hence the quadcopter is supposed to maintain the initial position and the initial attitude over time, that means no pitching, rolling nor yawing, nor vertical displacement. The Figure 6 show the results obtained by the numerical scheme (136), which confirm the intuitive prediction.
(a) Pitch, roll and yaw angles over time.
(b) Coordinates of the center of the quadrotor. In the Experiment 2, n 1 = n 3 = n ss , while n 2 = n ss − 1 and n 4 = n ss + 1. The unbalance in the rotation speed of the rotors R 2 and R 4 causes the quadcopter to rotate along the z axis (to yaw) and to tilt along the x axis (to roll) which, in turn, causes the quadrotor to move along the y axis. Since n 4 > n 2 , the rotor R 4 pulls stronger than the rotor R 2 , hence the drone actually shifts in the −y direction. The Figure 7 show the results obtained by the numerical scheme (136). Notice that the model (122) does not include any controller, therefore the slight unbalance in the propellers RPM causes an unbalance between the thrust and the gravitational attraction which, in turn, causes a drop in the quadcopter altitude. (For an example of proportional-derivative controller employed in quadcopter maneuvering, see [45].) As it was predicted by way of Lemma 7, the yaw is negligible.  In the Experiment 3, n 2 = n 4 = n ss , while n 1 = n ss − 1 and n 3 = n ss + 1. The unbalance in the rotation speed of the rotors R 1 and R 3 causes the quadcopter to yaw and to tilt along the y direction (that is, to pitch). The tilting causes the center of the quadrotor drone to shift along the x axis.
The Figure 8 show the results obtained by the numerical scheme (136). Since n 3 > n 1 , the rotor R 3 pulls stronger than the rotor R 1 , hence (readers might want to refer to Figure 3) the drone actually shifts in the +x direction. Again notice that, because the model (122) does not include any feedback controller, the slight unbalance in the propellers RPM causes a drop in the quadcopter altitude, while, again, the yaw effect results to be negligible.
(a) Pitch, roll and yaw angles over time.
(b) Coordinates of the center of the quadrotor. In the Experiment 4, n 2 = n 3 = n ss + 1, while n 1 = n 4 = n ss − 1. The unbalance in the rotation speed of the rotors R 1 and R 3 causes the quadcopter to pitch and to shift toward the +x direction and, at the same time, the unbalance in the rotation speed of the rotors R 2 and R 4 causes the quadcopter to roll and to shift toward the +y direction. Due to the assumed perfect symmetry of the model, equal unbalance in the propellers rotation speed results in equal shifting along both x and y axes. The Figure 9 show the results obtained by the numerical scheme (136), which confirm the intuition. Because of the combined action of all propellers, in this case no drop in the quadcopter altitude nor yawing of the quactopter body are predicted by the model.  Through the numerical scheme, we can compare the values of the vertical net force predicted by the model through the relation (132) in two experiments, for example in Experiment 3 and Experiment 4. The Figure 10 show the results of this comparison. The force f z in Experiment 4 almost doubles the force predicted in Experiment 3, which explains the larger predicted drop in the altitude of the quadcopter in the Figure 9.

Conclusions
A large number of scientific problems are most naturally described by systems of first-order differential and algebraic equations. Numerical methods for solving initial value problems are extremely important, because analytic solutions exist only in very special cases (e.g., for linear equations, exact equations, non-exact equations with integrating factors, and equations that are homogeneous of degree 0). Any first-order initial value problem can be solved (or "integrated") numerically by discretizing the interval of integration into a number of sub-intervals, either with equally or non-equally distributed grid points. As the integration advances, the numerical solutions at the grid points are approximated.
In the recent years, it has come to the attention of researchers in applied sciences and engineering how classical numerical methods, designed to solve initial value problems on R n are unsuitable to solve differential equations formulated on differentiable manifolds. These differential equations represent compactly large systems of mixed non-linear differential/algebraic equations, where the algebraic sub-equations carry on mutual constraints among the variables.
In the present paper, we recalled the formulation of non-conservative system dynamics through the Lagrange-d'Alembert principle and the generalized Euler-Poincaré form of the system equation on Lie groups. The Euler-Poincaré equations for left-invariant Lagrangians were illustrated via derivation of the dynamical equations for a gyrostat satellite and for a quadcopter drone.
Numerical experiments carried out by way of a forward Euler method on a manifold and an explicit fourth-order Runge-Kutta method on manifold illustrated the numerical integration of the Euler-Poincaré equations on the special orthogonal group SO(3). A number of questions are still an open problem and will be the focus of future investigation: • The problem about the order of the numerical schemes, which may be studied either analytically and by means of numerical experiments. The present author has started an analysis of the order of some numerical schemes which rely on ad-hoc Taylor-series-like expansions [46]. It would also be informative to use examples for which the analytic solution is known, and then compare the two numerical solutions to the analytic solution individually. This would also allow for the error order of the Runge-Kutta method to be evaluated numerically. While it is true that closed-form solutions of ODEs on nonlinear spaces cannot be obtained in general, the solutions are available in special cases, as shown in [47].
• The problem of the convergence of the numerical schemes, which is quite involved, being the algorithms considered in this paper of complex non-linear nature. The numerical schemes were stable in the numerical simulations performed in the course of this study. In general, it is expected that they would enjoy/suffer from the same advantages/drawbacks of their precursors.

•
An accurate evaluation of the computational complexity of the forward Euler method and of the Runge-Kutta integration method.

•
It would be very interesting to design more model-based experimental settings to survey the behavior of the studied satellite gyrostat. For example, it would be interesting to study a control strategy to see if it is possible to set the inputs ω 1 , ω 2 , ω 3 so as to make the satellite gain a prescribed attitude with respect to the earth-fixed reference frame. On the same line, it would be interesting to enrich the model (119) with a (possibly differential) equation that describes the relationship between the motors control voltages and the rotors rotation speed, in order to make the control scheme more realistic.
Funding: This research received no external funding.