A geometric formulation of multirotor aerial vehicle dynamics

A geometric dynamic modeling framework for generic multirotor aerial vehicles (MAV), based on a modern Lie group formulation of classical screw theory, is presented. Our framework allows for a broad range of rotor-wing conﬁgurations: any number of rotors can be attached in arbitrary conﬁgurations to either the body or wings, with the rotors and wings also tiltable. Our framework takes into account all masses and inertias of the MAV body and rotors, and accounts for both rotor thrust forces and moments as well as external aerodynamic and other forces. Compared to existing methods, our Lie group framework possesses several practical advantages useful for applications ranging from design optimization to model identiﬁcation and trajectory optimization: (i) the dynamic equations can be easily transformed to coordinates of any reference frame; (ii) kinematic and mass-inertial parameters can be easily factored from the dynamic equations; (iii) exact, closed-form analytic derivatives of the dynamics with respect to the conﬁguration variables are easily derived. We demonstrate our systematic modeling procedure on examples of ﬁxed-tilt, variable-tilt, and hybrid MAVs with wings.


Nomenclature {a}
Reference frame a Multirotor aerial vehicles (MAV) come in a wide variety of designs and configurations, from the standard quadrotor-type unmanned aerial vehicle (UAV) to those with rotor axes aligned in different directions, some actively tiltable, possibly with fixed or moving wings attached to various parts of the vehicle (e.g., vertical takeoff and landing aircraft designs). With the growing diversity of designs, together with the increasingly dynamic and complex maneuvers that are being performed by MAVs, accurate mathematical models are becoming more and more critical to their control and design. Such models are essential not only for accurate physics-based simulation and design optimization, but also for applying advanced model-based nonlinear control laws and trajectory planning algorithms, which with advances in computational power are now becoming more and more amenable to real-time computation.
Most of the existing work on MAV dynamic models have been limited to the most basic designs, in which the vehicle is modeled as a single rigid body subject to rotor thrust, aerodynamic, and other external forces. For example, the dynamics of standard quadrotors are addressed in [2,11,12,16,24]. The work of [2,12] takes into account the dynamics of the rotor actuators. The work of [16] also considers the aerodynamic forces acting on the vehicle. All of the previous works assume the standard quadrotor configuration, i.e., the four rotor axes are parallel and symmetrically arranged about the vehicle center of mass. The dynamics of six rotor MAVs are addressed in [1], and eight-rotor MAVs in [25]; here the rotor axes are all assumed parallel and arranged symmetrically about the vehicle center of mass. More recently, [4,14,21,28] address the dynamics of simple omnidirectional aerial vehicles with rotors aligned in arbitrary directions, but only the mass-inertia of the main body is considered, and the dynamics is formulated in the conventional way with respect the center of mass frame only.
More recently, in an effort to overcome the inherent power limitations of wingless rotorcraft, MAV designs with both fixed and moving wings attached have been proposed in the literature. The authors of [7,10] derive the dynamics of a MAV with two fixed wings, while the authors of [6] derive the kinematics (but not the dynamics) of a quad tilt-wing MAV with four rotors. A systematic analysis of the dynamics of generic winged MAVs has yet to be addressed in the literature.
In practice many of the underlying design assumptions made in the previous works will usually not hold entirely: manufacturing tolerances and errors will result in, e.g., rotor axes that may not be exactly parallel to one another, or rotors may not be arranged in a perfectly symmetric fashion about the vehicle's center of mass, or some of the rotor axes may fail to pass directly through the rotor body's center of mass. In some cases deliberately designing a MAV in such an asymmetric or skew way, with different types of wings arranged in atypical configurations, may even enhance some aspect of the MAV's performance. For a recent survey, the reader is referred to [23].
In this paper we present a systematic modeling framework for generic MAVs that encompasses a wide spectrum of possible designs. We assume the MAVs are driven by a set of rotors, possibly with tiltable axes, and that an arbitrary number of fixed and adjustable wings may be attached to the vehicle body. A subset of the rotors may be attached to the wings, and some to the main MAV body. Our modeling framework considers the rigid multibody dynamics of both the vehicle body's mass and inertia as well as the masses and inertias of the rotor bodies; other articulated bodies that may be attached to the MAV-for example, an open chain manipulator-can also be treated straightforwardly within our framework. In addition to the thrust forces and moments generated by the rotors, our framework allows for the consideration of external aerodynamic forces that are applied to the MAV.
One of the distinguishing features of our approach is the adoption of modern screw-theoretic methods for multibody dynamic modeling, or more precisely, using notation and operations associated with the Lie group of rigid body motions, also commonly known as the Special Euclidean group SE(3). The use of such methods is by now quite standard in the robot dynamics literature, e.g., [15,[18][19][20], as the algorithms are expressed in a way that is invariant with respect to the placement of reference frames on the bodies, and closed-form expressions for derivatives of the dynamics are available for optimization and sensitivity analysis, among other applications. Here we apply the same concepts and transformations to MAV dynamic modeling, and derive a reference frame-invariant set of dynamic equations for generic winged and wingless MAV models. The dynamics of MAVs with articulated structures such as a multi-dof robot arm attached to the body can also be systematically included within this framework.
By adopting the Lie group framework and the ensuing rules for the transformation of six-dimensional velocities, forces, and inertias, the dynamics of complex MAVs can be derived in a straightforward and systematic way, without regard to the way in which reference frames are attached to the MAV. Moreover, the resulting equations are shown to be easily factored with respect to any of the model parameters, e.g., inertias of the body and rotors, rotor axes, aerodynamic constants, etc, and also easily differentiated. These features make our geometric formulation particularly attractive for applications such as design, calibration and identification, and also dynamics-based trajectory op-timization where analytic gradients involving the dynamics are required.
The paper is organized as follows. Section 2 reviews the formulation of the dynamic equations for a single rigid body using Lie group concepts and notation. Section 3 presents the main results of our dynamic modeling framework for a generic MAV that incorporates fixed-tilt, variable-tilt, and hybrid MAVs. In Section 4, we apply our generic framework, together with complete solutions, for representative MAV structures like the standard quadrotor and also more complex structures like quadrotors with tilting rotors and multirotor tilt-wing MAVs. Finally, we conclude the paper in Section 5.

Mathematical preliminaries
We begin with some basic definitions and remarks about the notation, mostly following that used in [20] with some variations. Let {a} denote a right-handed orthonormal coordinate frame defined as {a} := (o a ,î a , a ,k a ), with o a denoting the origin of the frame and (î a , a ,k a ) are linearly independent free vectors.
Given two frames {a} and {b} situated in space, we denote by T b a ∈ SE(3) the relative displacement of frame {a} with respect to frame {b}, where SE(3) is the matrix Lie group of rigid body transformations. An element T b a has the explicit form where R b a ∈ SO(3) is a 3 × 3 rotation matrix representing the orientation of frame {a} with respect to frame {b}, while p b a ∈ R 3 denotes the position of o a in {b}. A general rigid body motion on SE(3) is a combination of rotational and translational motion. The geometric treatment followed in this work allows us to describe these two aspects using a single mathematical object, called a twist, which is an element of the Lie algebra of SE(3). We denote by se(3) the Lie algebra of the Lie group SE(3), consisting of 4 × 4 matrices of the form where v ∈ R 3 represents linear velocity, and [ω] ∈ so(3) is a 3 × 3 skew-symmetric matrix representing angular velocity.
We can associate to any [ω] ∈ so(3) a vector ω ∈ R 3 using the isomorphism Consequently, we can always identify the twist [V] ∈ se(3) with a six-dimensional vector V ∈ R 6 of the form For convenience V will also sometimes be written V = (ω, v), and with an abuse of notation we will call both V and [V] a twist. If either {a} or {b} is a moving frame, then the twist of {a} with respect to frame {b}, expressed in {a} where the over-dot denotes differentiation with respect to time, and ω a,b a , v a,b a ∈ R 3 are respectively the angular and linear velocity components of the twist V a,b a . Another important connection between SE(3) and se(3) is the Chasles-Mozzi Theorem 1 , which states that every rigid body motion can be expressed as a screw motion. Mathematically, every T ∈ SE(3) can be expressed as a matrix exponential of the form where the Plucker coordinates S = (ω, v) correspond to an element of se(3), and θ is a scalar. Closed-form expressions for the matrix exponential and its inverse (i.e. the matrix logarithm) can be found in, e.g., [15,29]. We denote the dual space of se(3) by se * (3) which is the space of wrenches (i.e. generalized forces) that are dual entities to twists. The pairing between a wrench [W] ∈ se * (3) and a twist [V] ∈ se(3) is a scalar that corresponds to the total power that is supplied to generate the rigid body motion described by [V]. By duality, we can associate to every wrench [W] ∈ se * (3) a 6-dimensional vector 2 W ∈ (R 6 ) * given by where m, f ∈ (R 3 ) * are the moment and force components of the wrench, respectively. Similar to twists, we abusively call both [W] and W wrenches. The duality pairing of a wrench W and a twist V yields 1 Most textbooks contribute this theorem only to the work of M. Chasles in 1830, however it dates back to G. Mozzi in 1763.
2 Although one has that (R 6 ) * ∼ = R 6 , we denote in this work the space of wrenches by (R 6 ) * to stress their covector nature which is crucial to note as wrenches and twists change coordinates differently. which corresponds to the total mechanical power. Now consider a three-dimensional rigid body of mass m with a body-fixed frame {b} attached to the body's center of mass. Let I b ∈ R 3×3 be the 3 × 3 rotational inertia matrix of the rigid body expressed in {b}. We assume that {b} is chosen to be aligned with the principle directions of the inertia ellipsoid of the corresponding rigid body. Consequently, I b is a diagonal matrix.
Consider an inertial fixed frame {0} fixed to the ground. Let V b,0 b be the twist of the rigid body with respect to this ground frame expressed in {b}. Further, let W b,b ext be the net external wrench applied to the body expressed in {b}. The equations of motion of the rigid body motion can then be expressed as where G b b ∈ R 6×6 is the generalized inertia tensor of the rigid body defined by with 1 denoting the 3 × 3 identity matrix. The 6 × 6 matrix [ad V ] represents the adjoint action of the Lie algebra se(3) given by the following matrix: Note that the total kinetic energy of the rigid body can be compactly represented as Along solutions of the dynamical equation (9), the kinetic energy satisfies the power balance: The compact form of the dynamic equations (9) is equivalent to the well-known rigid body equations ext .
An important advantage of the Lie group formulation of rigid body motion presented so far and used extensively in this work is the ease of changing coordinates. Let {a} be any arbitrary frame one wishes to represent the dynamics of the rigid body associated to {b} in. Then, the relative spatial displacement T b a ∈ SE(3) (and its inverse T a b ) is all that is needed to express the dynamic equations in the new frame {a}. For that purpose, one uses the adjoint representation of the Lie group SE(3) given by the 6 × 6 matrix: where [p] denotes the skew-symmetric matrix representation of p given by (3). The coordinate transformation rules for twists, wrenches, and generalized inertias are as follows: In the special case that {a} is another arbitrary body-fixed coordinate frame attached to a different point on the rigid body, then one has that V a,b a = 0 and thus The rigid body dynamic equations in terms of the new body-fixed frame {a} then become G a aV a,0 which is exactly of the same form as (9). This forminvariance property under change of coordinates is an advantage of the differential geometric framework utilized to describe body motion. Note that in the new coordinate frame {a}, the generalized inertia G a a is a symmetric positive definite matrix but does not have the diagonal form in (10).

Dynamics of generic MAVs
In this section, we present the systematic modeling procedure for a generic MAV. As for many other types of robots, it is common to treat a MAV as a multibody system consisting of a set of rigid bodies linked together via joints that are usually activated by electric motors. In what follows, we describe the dynamic model of two rigid bodies connected by a 1 degree-of-freedom revolute joint, excluding any other externally applied wrench. Then, we show how these basic atomic units can be generalized systematically to model a complete generic MAV with different external sources.

Parent-Child Relations
Consider the parent rigid body, with body-fixed frame {p}, connected through a revolute joint to a child rigid body, with body-fixed frame {c}, as depicted in Fig. 1.
The revolute joint constraints the relative motion of the two rigid bodies with a unique twist V p,p c ∈ R 6 given by whereθ c ∈ R is the joint's velocity and S p,p c ∈ R 6 is a constant unit twist fixed in {p} and given by S p,p c = (n c , q c ×n c ), withn c ∈ S 2 denoting the unit vector that describes the rotation axis of the joint, with positive rotation defined in the sense of the right-hand rule, and q c ∈ R 3 is any vector from the origin of {p} to any point on the axisn c (bothn c and q c are expressed in coordinates of {p}).
Let T p c ∈ SE(3) denote the relative displacement between the two frames {c} and {p}, then one has that where θ c is the joint's angular displacement and T p c (0) is the relative displacement of the two frames when θ c is zero.
A very important distinction between the parent and child rigid bodies is that the twist of the child body V c,0 c ∈ R 6 is completely determined by the parent's motion and the joint's velocity: On the other hand, the twist of the parent body V p,0 p ∈ R 6 is determined by the differential equation: where W p,p c ∈ (R 6 ) * is the reaction wrench exerted by the child body on the parent body, expressed in {p}. This reaction wrench corresponds to the inertial forces and moments associated with the relative motion of the child's body (as well as other external wrenches applied to it which were not considered so far).
The reaction wrench W p,p c is computed by with W c,c p ∈ (R 6 ) * denoting the reaction wrench that the parent body exerts on the child body (expressed in {c}) which is given by with V c,0 c given by (22) and its rate of change given by the closed-form expressioṅ where S c,p c is the joint's unit twist expressed in {c} and related to S p,p By substituting equations (22,25,26) into (24) and expressing all quantities related to the child body in {p} using (15)(16)(17)(18), one can express the wrench W p,p c applied by the child on the parent as where G p c is the child's generalized inertia expressed in {p}, and X p c ∈ (R 6 ) * is given by Finally using (28), one can re-express the parent's dynamic equations of motion (23) as with the total (or locked) inertia tensor G p T ∈ R 6×6 given by We conclude this section by some remarks on the three basic units mentioned in this section, with reference to Fig. 1. First, while both the parent and child bodies are governed by the rigid body dynamics (9), they differ in the causality. In particular, the wrench is an input to the parent's dynamic model (23) while the twist is an output which is in contrast to the child's dynamic model (25). Thus, we refer to the parent's dynamics to be in integral causality while the child's dynamics to be in differential causality. The reader is referred to [8] for more details on this terminology and its relation to energy-based modeling using bond graphs. Second, the differential geometric approach employed above allows a compact representation of the motion of two rigid bodies in space. Furthermore, one can compact (22) and (24) into the following relation with the additional output τ c ∈ R denoting the torque applied by the joint to generate the motion, while J p c (θ c ) is the matrix defined as with T p c given by (21). Therefore, J p c (θ c ) is determined completely by θ c and is constructed by the constant matrices S p,p c and T p c (0) which are known by definition of the joint's rotation axes and the choice of the reference frames {p} and {c}.

Generic MAV with wings
A generic MAV model with n rotors and m wings is depicted in Fig. 2. The MAV's rotors could be attached to the main body or to one of the wings. Such generic model encompasses fixed-tilt MAVs, variable-tilt MAVs, and hybrid MAVs with wings. In our framework, we treat the tilting mechanisms of a variable-tilt MAV as wings with zero surface area.
We denote by n α the number of rotors connected to the α-th wing, and by n b the number of rotors directly connected to the body, so that n = n b + m α=1 n α . Let N α be the set of indices for the rotors connected to wing α, and N b be the set of indices for the rotors directly connected to the body. In what follows, we will use the index α = {1, · · · , m} for all the physical quantities associated with the wings. While the indices i ∈ N b and j ∈ N α will be used for the corresponding quantities associated with the body-rotors and the α-th wing-rotors, respectively.
Let frame {0} denote the inertial (ground) frame, while {b} is a moving frame (the "body" frame) attached to the MAV's main body. Frame {w α } is a moving frame attached to wing α, for α ∈ {1, · · · , m}, whereas frames {r i } and {r j } are moving frames attached respectively to rotors i, for i ∈ N b , and rotors j, for j ∈ N α . All of the aforementioned frames can be arbitrarily oriented and attached to their respective bodies. Furthermore, the tilting axes of the wings and the rotation axes of the rotors can be arbitrary. The angles for wings are denoted by φ α , for α = 1, ..., m, whereas the angles for rotors are denoted by θ i for i ∈ N b and θ j for j ∈ N α . Now we show how the parent-child dynamics described in the previous section can be employed to systematically construct the dynamic model for the generic MAV described above. Possible combinations of parentchild frames (p,c) include (b,α) for α = 1, . . . , m, (α, j) for j ∈ N α , and (b, i) for i ∈ N b . The overall dynamic model is depicted in Fig. 3. Next, we explicate the individual submodels corresponding to the MAV's components.

Main Body Dynamics
First we consider the dynamic modeling of the MAV's main body. Similar to the parent dynamics previously discussed in (23), the dynamic equations for the MAV's body (expressed in {b}) are in the integral causality form given by where W b,b α , W b,b i ∈ (R 6 ) * denote the reaction wrenches exerted on the main body by the wing α and body-rotor i, respectively. Further, W b,b ext ∈ (R 6 ) * denotes the external wrench (with respect to the MAV as a whole) applied to the main body which is given by the sum of the gravity wrench W b,b grv and the aerodynamic wrench aer . The wrench W b,b grv ∈ (R 6 ) * exerted by gravity on the body, expressed in {b}, is given by where A 0 grv = (0, a 0 ) ∈ R 6 , with a 0 ∈ R 3 denoting the gravitational acceleration vector expressed in the ground frame {0}. Note that in the special case that {b} was chosen to coincide with the center of mass of the body, the gravity wrench takes the more usual form However, the expression in (35) accounts for the additional moments due to the arbitrary choice of {b}.
The aerodynamic wrench W b,b aer in turn does not solely depend on the body's ground speed, corresponding to V b,0 b , but rather its air-speed, corresponding to V b,a b which denotes the relative twist of {b} with respect to the air given by where V 0,0 a = (0, v 0,0 a ) denotes the local velocity of the wind with respect to the inertial frame in a region upstream of the MAV's body. In case the wind is neglected (i.e. V 0,0 a = 0), the ground speed and air speed coincide with each other.
The aerodynamic wrench W b,b aer acting on the main body is usually due to drag that is modeled as a nonlinear function of the twist V b,a b that depends on properties of the body like its shape and material, and whose precise form is obtained from aerodynamic experiments. Depending on V b,a b , the drag force can be dominated by either the viscous force term (which is linearly dependent on the speed) or the inertial drag term (which is quadratically dependent on the speed).

Rotor Dynamics
Next, we consider the dynamic modeling of the rotors attached directly to the MAV's body through revolute joints.
Let T b i ∈ SE(3) be the relative displacement between the MAV's body frame {b} and the rotor frame {r i }, for i ∈ N b . Denoting by θ i the angle of rotation for rotor i, one can express T b i in the following exponential form: where (3) is the displacement of rotor frame {r i } relative to the body frame {b} when θ i is zero and S b,b i ∈ R 6 is the constant screw vector corresponding to rotor's rotation axis as seen from {b}.
In our presented framework, one has to provide the constant matrices T b i (0) and S b,b i based on the choice of the reference frames and the definition of the rotor's axis. In applications, it is more convenient to define the unit twist S i,b i corresponding to the rotor's axis expressed in {r i }. Thus, one can defined S i,b i = (n i , q i × n i ) , wheren i ∈ R 3 is a unit vector in the direction of the rotor axis and q i ∈ R 3 is any vector from the {r i } frame origin to any point on rotor axis i. Then, similar to (27), one could compute S b,b i by Ideally the rotor frame would be placed such that its origin lies on the rotor axis, in which case q i will be zero. However, with this more general formulation one can reflect, e.g., possible manufacturing tolerances and other errors introduced in the rotor assembly modeling that are inevitable in practice. Using i and the joint angle, one can construct the joint matrix J b i (θ i ) given by (33), replacing the scripts c and p by i and b, respectively. The joint matrix J b i (θ i ) will be essential for relating the main body's ground twist V b,0 b , the rotor's ground twist V i,0 i , and the joint's velocityθ i , as well as the reaction wrench on the main body W b,b i , the reaction wrench on the rotor W i,i b , and the joint's torque τ i . Based on our earlier geometric treatment of the dynamics of a child rigid body in (25), the dynamic equations for rotor i (expressed in {r i }) are in the differential causality form where W i,i b is the reaction wrench applied by the MAV's main body to the rotor and W i,i ext is the wrench applied due to external sources given by the sum of the gravity wrench W i,i grv and the aerodynamic wrench W i,i aer . The gravity wrench W i,i grv has the exact form as (35) with the script b replaced by i. On the other hand, the aerodynamic wrench W i,i aer depends on the rotor's relative air twist V i,a i which in turn depends on the rotational velocityθ i and the main body's relative air In practice, MAVs are usually operated in nearhovering regimes and thus it is common to neglect the main body's relative air twist and experimentally fit the aerodynamic wrench W i,i aer as a quadratic function with respect to the rotor's rotational velocityθ i . As outlined in [3,11,13,17], for MAVs W i,i aer can be written as where K i ∈ R 6×6 is the following constant matrix: where i = ±1 depends on whether the rotor rotates counterclockwise (+1) or clockwise (−1) with respect to the rotor axis directionn i , and k di and k li are constants whose values are specified by the propeller geometry, air density, air velocity, etc. [13]. Note that in the special case in which the coordinate frame {r i } is placed along the rotation axisn i , then one has that q i is zero and the aerodynamic wrench W i,i aer takes the usual form In practice however, one could have that q i = 0 as mentioned above and consequently an extra moment appears in (42) compared to the ideal case (43). Finally, all the previous constructions in (40-43) hold also for each rotor j ∈ N α attached to the α-th wing with the indices i and b replaced by j and α, respectively. Following the same line of thought presented above, one can construct the joint matrix J α j (θ j ) using the constant matrices T α j (0) and S α,α j based on the choice of the reference frames and the definition of the rotor's axis. Then, the dynamic model for rotor j is given by

Wing and Tilt-Mechanism Dynamics
Now we consider the dynamic modeling of the wings and/or tilting mechanisms of the MAV. In our framework, the tilting mechanism of a variable-tilt rotor is considered as a wing with zero surface area. Thus, the only distinction in its dynamic model compared to a general wing is the absence of the aerodynamic wrench. Let φ α ∈ R denote the tilting angle of wing α, S b,b α ∈ R 6 denote the constant screw vector corresponding to the wing's tilting axis, and T b α (0) ∈ SE(3) denote the relative displacement of the wing frame {w α } and the body frame {b} when φ α = 0. Using (33), one can construct the joint matrix J b α (φ α ) which relates the main body's ground twist V b,0 b , the wing's ground twist V α,0 α , and the joint's velocityφ α , as well as the reaction wrench on the main body W b,b α , the reaction wrench on the wing W α,α b , and the joint's torque τ α .
The dynamic equations for wing α are given in differential causality form by where W α,α j is the reaction wrench applied on the wing by rotor j attached to it, which is related to its counterpart W j,j α in (44) by the joint matrix J α j (θ j ). The wrench W α,α ext applied to the wing due to external sources is given by the sum of the gravity wrench W α,α grv , similar to (35), and the aerodynamic wrench W α,α aer . In case the wing α represents the tilting mechanism of a rotor, then one has that W α,α aer = 0. Next we discuss how the aerodynamic wrench W α,α aer is calculated for general wings. The primary difference between winged MAV models considered in this work and typical fixed-wing aerial vehicles is that in the former, aerodynamic forces need to be considered in the dynamic equations for each wing. In typical fixed-wing designs where identical wings may be multiply attached in symmetric configurations, aerodynamic forces can usually be derived directly with respect to the main body frame. This is no longer the case for MAVs with wings of different shapes and sizes that in general are actively tilted; aerodynamic forces need to be considered wing-by-wing for such generic MAV model.
The aerodynamic wrench applied on each wing will typically be calculated from the dimensionless aerodynamic coefficients of the finite-wing. However, a crucial issue that should be considered is the coordinate-frame in which these aerodynamic coefficients are represented. This issue is sometimes a source of confusion since the terminology used in the aerodynamic literature to describe the components of the aerodynamic wrench implicitly defines the coordinate-frame.
Note that in this section so far, our treatment allowed an arbitrary location and orientation of the wing frame {w α }. However, it is standard to define the aerodynamic coefficients with respect to a frame located at the wing's center of mass. For that purpose, we introduce another body-fixed reference frame attached to wing α, denoted by {w α }, that is located at the wing's center of the mass and oriented such that theî α axis is aligned with the wing's chord, thek α axis is normal to the wing's chord, andĵ α completes the coordinateframe as depicted in Fig. 4.
In terms of the coordinate-frame {w α }, the wing's ground twist and aerodynamic wrench applied to it are related to their counterparts in the arbitrary frame {w α } by The aerodynamic wrench W α ,α aer acting on the wing depends on two factors; the relative motion of the wing with respect to the air as well as on the orientation of the wing with respect to the air flow. The first factor is represented by the wing's twist relative to the air V α ,a α given by similar to (37). The magnitude of v α ,a α , the linear velocity component of the twist V α ,a α , is known as the free stream speed V ∞ ∈ R, i.e.
The second factor is represented by two aerodynamic angles; the angle of attack β aa and the side-slip angle β ss . If we denote the individual components of v α ,a α by (v α x , v α y , v α z ), then the aerodynamic angles are defined by Therefore, the aerodynamic wrench W α ,α aer can be expressed in the form [7,9,26] W α ,α where ρ ∞ is the free-stream air density, S α is the planform area of wing α (shown shaded in Fig. 4), V ∞ is the free-stream speed, and C α ∈ (R 6 ) * is the sixdimensional vector of aerodynamic coefficients: wherec α , b α are the mean aerodynamic chord and the span of wing α, respectively, while the rest of the coefficients are defined in Table 1. An alternative and more common way to define the aerodynamic coefficients is in another frame that is aligned with the direction of the free-stream velocity instead of being aligned with the wing's geometry. Using the aerodynamic angles (50), it is common in the literature to introduce another frame {wᾱ} that is also located at the center of mass of the wing but oriented relative to {w α } by: where Rot(z, β), Rot(y, β) ∈ SO(3) denote, respectively, the standard rotation matrices about the z-axis and yaxis by β. Let Vᾱ ,ā α denote the wing's twist relative to the air and Wᾱ ,ᾱ aer denote the aerodynamic wrench applied to the wing, both expressed in the new frame {wᾱ} and related to their counterparts in {w α } by In the frame {wᾱ}, the linear velocity component of Vᾱ ,ā α takes the form vᾱ ,ā α = (V ∞ , 0, 0) ∈ R 3 , while the aerodynamic wrench Wᾱ ,ᾱ aer is expressed by where Cᾱ ∈ (R 6 ) * is the six-dimensional vector of aerodynamic coefficients: where the individual coefficients are defined in Table 1. Note that it is quite common in the aerodynamic literature that the same notation is used for the moment coefficients as in (52) and (57), although they are represented in two different frames [26]. In summary, the calculation of the wrench W α ,α aer depends on the coordinate-frame chosen to represent the dimensionless coefficients. The representation of the aerodynamic coefficients in terms of the normal, axial, and side-force components implies that one should use (51, 52) to calculate the aerodynamic wrench applied to the wing. On the other hand, the representation of the drag, cross-wind, and lift components implies that one should use (55, 56,57) to calculate the aerodynamic wrench.
In general, the aerodynamic coefficients are obtained from experimental wind-tunnel measurements or using computational techniques. Furthermore, they depend on the attack angle β aa , the sideslip angle β ss , as well as on the Reynolds Mach numbers. Another important point is that the aerodynamic coefficients (57) or (52) are in fact static coefficients that do not include the effects of the MAV maneuvering or the wings being actively tilted. More accurate models that address such unsteady aerodynamic effects are thus needed for a generic MAV with wings. We refer the reader to [9,26] for further discussion on static and dynamic aerodynamic coefficients.

Compact MAV body model
With the constructions presented so far, one can develop the dynamic model for a generic MAV in the framework shown in Fig. 3. The framework's compactness and modularity is beneficial for conceptual analysis and simulation purposes.
In case it is desired to compact the MAV's dynamic model, e.g. for control or identification purposes, then the dynamic equation (40) for each rotor i connected to the main body should be represented in {b} and substituted in (34). Then, the dynamic equation (44) for each rotor j, connected to wing α, should be represented in {w α } and substituted in (45). Next, one would represent the result of (45) in {b} and then substitute it in the dynamic equation (34). These steps are summarized in Algorithm 1.
To demonstrate Algorithm 1, consider the special case of a generic fixed-tilt MAV depicted in Fig.5. The MAV has n rotors all attached to the main body with no wings or tilting mechanisms i.e. i ∈ N b = {1, · · · , n} and m = 0.
By applying lines (1)(2)(3)(4)(5)(6) in Algorithm 1, the dynamic equations for the MAV body including all rotor wrenches Algorithm 1 Dynamics of a generic MAV model i into (34); 6: end for 7: for α = 1 : m do 8: for j ∈ N α do 9: Construct J α j (θ j ) using T α j (0) and S α,α j ; 10: Evaluate W j,j α from (44); 11: Compute W α,α j and τ j using J α j (θ j ); 12: Substitute W α,α j into (45); 13: end for 14: Construct J b α (φ α ) using T b α (0) and S b,b α ; 15: Evaluate W α,α b from (45); 16: Compute W b,b α and τ α using J b α (φ α ); 17: Substitute W b,b α into (34); 18: end for can now be written in frame {b} coordinates as follows: Following the same line of thought as in (28), in addition to using the transformation rules (15)(16)(17)(18) and the expressions for the external wrenches W b,b ext and W i,i ext , the above dynamic equation can be compactly rewritten as follows: where the inertia tensors G b T and G b i are given by whereas the matrixK i ∈ R 6×6 is given bȳ A straightforward calculation shows in fact thatK i = K i . Furthermore, the screw vectors S b,b i and S i,b i are related to each other by (39), and finally X b i ∈ (R 6 ) * is equivalent to (29) with c and p replaced by i and b, respectively.
For the small-scale fixed MAVs usually used in research, the aerodynamic wrench W b,b aer applied on the body due to air drag is usually neglected in near-hovering situations. Also, the inertial-coupling effects of the rotors are typically quite small compared to other wrenches acting the body, so that they can be treated as zero in many situations. In this case the dynamic equation (58) reduces to whereθ 2 := [θ 2 1 , · · · ,θ 2 n ] ∈ R n is a vector of the squared rotor speeds and K ∈ R 6×n is the constant control allocation matrix with its i-th column given by i . The above simplified formulation (62) can be adopted to many practical MAV control applications. Often to further simplify matters, it is common practice to formulate the translational dynamics in terms of a bodyfixed frame attached at the center of mass, with the further assumption that the rotor axes are placed symmetrically about the center of mass. One of the advantages of Lie group formulation presented in our work is that the dynamics can be formulated independent of the choice of body frame, and without assuming the rotors are symmetrically placed or aligned in the same direction. Note also that all the model parameters in (62), i.e., the screw parameters of the rotor axes S b,b i , matrices of aerodynamic coefficients K i , and the massinertial parameters G b T , appear in explicit and linear form, independently of each other. This feature is particularly important for model identification and adaptive control purposes.

Examples
Based on the dynamic modeling framework formulated above, we now derive the dynamics for four MAV designs: (1) a standard quadrotor with symmetric parallel rotor axes; (2) a standard quadrotor with tilting rotors; (3) a quadrotor with two wings; and (4) a quadrotor with four tilting wings. The four examples highlight the applicability of our framework to fixed-tilt, variable-tilt, and hybrid MAVs.

Standard quadrotor
A standard quadrotor together with its assigned body and rotor frames is depicted in Fig. 6. Assume that orientations of all rotor frames are identical to that of {b}, and each of the rotor and body frames are attached to its corresponding center of mass. The first and third rotors rotate clockwise, while the second and fourth rotors rotate counterclockwise, so that i = (−1) i . The four rotors are parallel and symmetrically located about its center of mass. Denote by d the distance between the body's center of mass and each of the rotor axes. If each rotor is assumed to be ideally balanced about its center of mass, each q i will be zero, and the screw axis for each rotor is given by for i ∈ N b = {1, . . . , 4}. The initial relative configuration T b i (0) of each frame {r i } with respect to {b} is defined by the pair where ψ i ∈ {π, −π/2, 0, π/2}. Substituting (63) and (64) into (38), the orientation and position of each rotor frame for a given rotor angle θ i is given by Furthermore, using (63-65), one has that the gyroscopic wrench X b i exerted by each rotor on the body, defined in (29), simplifies to Under the assumption that the frames {b} and {r i } are aligned with the principal axes of inertia and all rotors have identical mass-distributions, one has that the inertia tensors G b b and G i i are diagonal, i.e., where I bx , I by , and I bz denote the moments of inertia of the body aboutî b , b ,k b , respectively, and I rx , I ry , and I rz denote the moments of inertia of the rotor aboutî i , i ,k i , respectively. Further, m b denotes the main body mass and m r the rotor mass. The dynamic equations for this fixed-tilt MAV are given in (58-60), which can be rewritten as where where m T = m b + 4m r , I T x = I bx + 4I rx , and I T z = I bz + 4I rz .

Quadrotor with tilting rotors
The variable-tilt quadrotor shown in Fig. 7 has four rotors whose axes can be tilted [24]. We assume that four identical rotors are symmetrically attached about the body center of mass, with the first and third rotors rotating clockwise, and the second and fourth rotors rotating counterclockwise, so that as in the previous example we have j = (−1) j . All frames including the body, rotors, and tilting assemblies are assumed attached to the corresponding center of mass. If the four tilting assemblies including each rotor are considered as a type of wing with zero platform area, then Algorithm 1 can be used to derive the dynamics for this quadrotor. For the rotor frames {r j } and wing frames {w α } shown in Fig. 7, the configuration of each of the wings and rotors for α ∈ {1, · · · , 4} and j = α is given by Using Lines 8-15 of Algorithm 1 to derive the dynamic equations for each tilting assembly, we multiply −[Ad T j α ] to both sides of (44) and substitute the result into (45), so that all terms are now expressed in {w α }. Then, we have that where the aerodynamic wrench W α,α aer was neglected and the last term on the right-hand-side corresponds to the gyroscopic wrench similar to (66).
Applying Lines 16-17 of Algorithm 1 to (77), the final form of the MAV dynamic equations expressed in the main body's frame {b} becomes where the last three wrenches are given by while the wings' tilting axes and rotors' rotation axes are given by and finally the wings', rotor's and total inertia tensors are given by If the rotor masses and inertias are much smaller then those of the body such that one can neglect the inertia coupling wrench, the tilting angles φ α vary in time slow enough such thatφ α andφ α are negligible, and the drag wrench on the main body is insignificant, the dynamics (78) further simplifies to

Quadrotor with two fixed wings
For this example we consider a quadrotor with two wings and four rotors as shown in Fig. 8. Two of the four rotors are connected to one wing, while the other two are connected to the other wing. Assume that all rotors are ideal, rotating in the way described in the preceding examples. All frames are attached to the corresponding center of mass, and their orientations are the same as that of the body frame. Under these assumptions, the screw axis for each rotor S i,b i is given by (63) and the current configuration T b i is given by (65);for example, p b i=1 is (−d, h, w) as shown in Fig.  8.
Since all the wings are fixed, we have for α = 1, 2, thatφ α = 0 and consequently V b,0 b and V α,0 α are the same, while T b α are constant transformations. Therefore, the MAV dynamics are given by where while Y b and X b are the same as given in (71) and (72), respectively, and Kθ A fundamental difference between the quadrotor with two fixed wings of this example and the standard quadrotor in Sec. 4.1 is that the rotors are mostly used in the latter to counteract gravity in near-hovering regimes. On the other hand, the quadrotor with two fixed wings uses the rotors to provide thrust in thek b direction only, as indicated in (80), and is usually operated at high forward-speeds. Consequently, using a static hovering model as in (42)  for estimating the rotors' aerodynamic wrench for such MAVs. Instead, the aerodynamic wrench W i,i aer should take into account also the main body's relative air twist in (41) and the aerodynamic drag wrench on the body aer should be also be considered.

Quadrotor with four tilting wings
We now consider the quadrotor design shown in Fig. 9 and first presented in [6]. This particular quadrotor consists of four tilting wings and four rotors. Each rotor is attached to one of the wings, with the four wings identical and symmetrically arranged. Assume that all rotors are ideal, rotating in the way described in the preceding examples, and the rotating axis of each wing passes through the wing's center of mass along the dashed lines indicated in Fig. 9. All frames including the body, rotors, and wings are assumed attached to the corresponding center of mass. This MAV with four tilting platforms and four rotors is quite similar in structure to the quadrotor with tilting rotors addressed in Sect. 4.3. The dynamic equations are in fact structurally identical to (78), since the tilting assemblies in Sect. 4.3 correspond to the wings of this quadrotor, with each of the tilting assemblies (and wings) possessing only one rotor. The only remaining difference is the presence of the aerodynamic wrenches W α,α aer in (77) and (78), since unlike the tilting rotor case, the planform area of each wing is not zero. Therefore, using the geometric characteristics and aerodynamic coefficients of the wings, one compute the wings' aerodynamic wrench as described in Sec. 3.2.3.

Conclusion
With the growing diversity of designs for multirotor aerial vehicles (MAVs), together with the increasingly dynamic and complex maneuvers that are being performed by MAVs, accurate dynamical models are becoming critical to MAV trajectory planning, control, Fig. 10 Bond graph representation of a generic MAV equivalent to the block-diagram representation shown in Fig.3. The bond graph notation follows [8].
and design. In this paper we have presented a systematic dynamic modeling framework for generic MAVs, ranging from the standard fixed-tilt quadrotor to wingless MAVs with tiltable rotors, and to hybrid MAVs with tiltable wings. Our approach considers the dynamics of both the vehicle body's mass and inertia as well as the masses and inertias of the rotors. Our framework does not rely on the standard design assumptions, and allows for the possibility that, e.g., rotor axes may not be exactly parallel to one another or arranged in a perfectly symmetric fashion about the vehicle's center of mass, or some of the rotor axes may fail to pass directly through the rotor body's center of mass. Both thrust forces and moments generated by the rotors, as well as external aerodynamic and other forces applied to the MAV, are taken into account. Furthermore, our generic mathematical representations highlight the hidden modeling assumptions behind widely used dynamic models in the literature.
A key feature of our approach is the adoption of notation and operations associated with the Lie group of rigid body motions SE(3). The Lie group framework allows for, among other things, a reference frameinvariant set of dynamic equations that can be easily transformed to coordinates of any arbitrary reference frame by using the tensorial transformation rules for six-dimensional velocities, forces, and inertias. The kinematic and mass-inertial parameters can also be transparently factored in linear form, making our model useful for model identification and calibration purposes.
A second important advantage of our geometric framework is the ease with which the dynamics can be differentiated in exact fashion; this feature allows for exact analytic gradients to be computed for optimal control and trajectory optimization problems, which is critical for fast and numerically robust computation of optimal trajectories. An important point that was highlighted in our work is the geometric dual structure of twists and wrenches. The exploitation of this structure for modeling of dynamical systems in general is the main reason for the modularity that is clearly shown in Fig. 3. This duality between twists and wrenches is essential also for an energy-based interpretation of the MAV's dynamic model. The pairing between every twist and its corresponding wrench yields the power-flow network of the system, shown using the bond graphs notation in Fig.10. This alternative view of dynamical systems have proven very effective for modeling, analysis and control purposes within the framework of port-Hamiltonian systems. The interested reader in this topic can refer to [8,22,27].
Finally, an interesting extension of our presented framework is to include flexible flapping-wings. While the kinematics and inertial effects of rigid flappingwings can be handled by our framework, the true challenge of such complex dynamical systems lies in the phenomenological relation between the aerodynamic wrench on the wing and the wing's relative air twist in addition to modeling the wing's deformation. Such additions require the coupling of the finite-dimensional dynamics presented in this work to the infinite-dimensional dynamics of air flow and flexible structures. The interested reader can refer to [5] for a view of how energy-based modeling is well-suited for such complex MAVs.