## Abstract

This paper presents a rigorous analysis of a promising bi-modal multirotor vehicle that can roll and fly. This class of vehicle provides energetic and locomotive advantages over traditional unimodal vehicles. Despite superficial similarities to traditional multirotor vehicles, the dynamics of the vehicle analyzed herein differ substantially. This paper is the first to offer a complete and rigorous derivation, simulation, and validation of the vehicle's terrestrial rolling dynamics. Variational mechanics is used to develop a six degrees-of-freedom dynamic model of the vehicle subject to kinematic rolling constraints and various nonconservative forces. The resulting dynamic system is determined to be differentially flat and the flat outputs of the vehicle are derived. A functional hardware embodiment of the vehicle is constructed, from which empirical motion data are obtained via odometry and inertial sensing. A numerical simulation of the dynamic model is executed, which accurately predicts complex dynamic phenomena observed in the empirical data, such as gravitational and gyroscopic nonlinearities; the comparison of simulation results to empirical data validates the dynamic model.

## 1 Introduction

A bi-modal vehicle leverages complementary modes of transportation by adapting its locomotion to suit the environment. Locomotion adaptability makes bi-modal vehicles well suited for exploratory operations involving unknown terrain [15]. Among bi-modal vehicles, those with aerial capabilities provide unmatched mobility in uncertain and changing environments. As a result, researchers have developed several classes of flying bi-modal vehicles, including rolling-flying [610], walking-flying [11,12], and swimming-flying [1315], each of which utilize various mechanisms for propulsion and lift. Furthermore, vehicles with the ability to hover, such as rotary wing vehicles, demonstrate excellent maneuverability. Such vehicles may fly across impassable terrain, translate in any direction, and be position-controlled with high resolution at zero velocity. However, the penalty for great maneuverability is sustained high power consumption, making rotary wing vehicles particularly inefficient at low speeds. On the other hand, rolling is a supremely efficient means of locomotion [16]. Therefore, a bi-modal vehicle that possesses the complementary capabilities of rolling and rotary wing flight achieves mobility, maneuverability, and efficiency. For example, an energy-efficient rolling vehicle could temporarily fly upon encountering an obstacle that cannot be traversed by rolling; while in flight, the vehicle is trading locomotion efficiency for increased mobility.

A promising embodiment of a rotary wing rolling-flying vehicle (RFV) is illustrated in Fig. 1. This class of vehicle was first conceived by Kalantari and Spenko [6] when they created a hybrid terrestrial/aerial quadrotor (Hy-TAQ), which is a multirotor vehicle that utilizes rotary wings for both propulsion and lift generation. The RFV rolls via two independent, passive wheels on either side of the vehicle. This configuration allows the rotors to be oriented independently of the wheels' angular orientation and combines the efficiency of rolling locomotion with the three-dimensional mobility of rotary wing flight to traverse obstacles or access elevated locations when necessary. The energetic efficiency of the RFV compared to a conventional flying multirotor vehicle (CFMV) has been established along with methods for determining optimal operation on any terrain [17].

Fig. 1
Fig. 1

However, these benefits can only be realized if the RFV's terrestrial motion is controlled in closed-loop. This motivates the development of a rigorous dynamic model of the RFV while rolling. While mathematical models of CFMVs have been studied extensively [1821], the methods and assumptions used to develop models for CFMVs are not valid for an RFV, despite the RFV utilizing the same multirotor mechanism for both rolling and flying. For example, the near-horizontal nominal operating condition of CFMVs permits significant simplification of the dynamics for the purposes of developing control systems [18,19]. However, the dynamic model, flat outputs [22], exogenous forces, and range of motion associated with the RFV's rolling mode differ substantially from those of a CFMV due to the presence of kinematic constraints imposed upon the RFV. To understand these constraints and their impact requires a unique dynamic model.

While other researchers have created rolling-flying vehicles like that shown in Fig. 1, none have developed a dynamic model suitable for terrestrial control system design that also includes the relevant constraints, rotational dynamics, and nonconservative forces. Furthermore, none have validated their dynamic models with experimental data. Kalantari and Spenko [6] invented and patented HyTAQ, consisting of a quadrotor suspended within a single rolling cage that may roll along the ground and fly. The RFV presented here is one of the many embodiments covered by their patent [6]; specifically, the RFV uses two independently passive rotating wheels rather than a single monolithic cage. Takahashi et al. [7] developed the all-round two-wheeled unmanned aerial vehicle (UAV) that contains hemispherical wheels that envelop a UAV suspended along the wheels' axle. While Kalantari and Spenko [6] formulate Lagrange's equations for the terrestrial motion of their single-wheeled design, they do not present their model in a form suitable for control system development, i.e., they do not solve Lagrange's equations for the derivatives of the generalized coordinates. Indeed, their subsequent analysis neglects all dynamic behavior by assuming static operating conditions, and no dynamic simulations are performed. Takahashi et al. [7] create a two-wheeled design but present a simplified terrestrial model that neglects rotational dynamics and nonholonomic constraints. Additionally, the model excludes the wheels entirely. Furthermore, Takahashi et al. [7] make no attempt to measure the velocity either of the wheels or of the vehicle. Mizutani et al. [23] develop a quadrotor UAV with spherical shell capable of rolling; however, they present no model. Many researchers develop dynamic models of generic two-wheeled mobile vehicles for strictly terrestrial locomotion. However, these models typically confine the vehicle's operation to either a horizontal (e.g., Ref. [24]) or vertical (e.g., Ref. [25]) two-dimensional plane. Thus, these models cannot capture the higher dimensional terrestrial motion of the RFV.

This paper extends the work of Kalantari and Spenko [6] and Takahashi et al. [7] and takes the first step toward developing appropriate terrestrial control systems for the RFV by providing an accurate dynamic model of the RFV while rolling. The dynamic model is rigorously derived using the Euler–Lagrange method of variational mechanics. The model takes into consideration the various constraints imposed upon the vehicle while in contact with the ground and includes nonconservative forces. Furthermore, the model does not assume a near-horizontal nominal operating condition, nor constrain the motion of the vehicle to a two-dimensional plane. The result is one first-order and two second-order differential equations that completely govern the RFV's terrestrial motion. The implications of the constraints as they relate to differential flatness are explored. A hardware embodiment of the RFV is developed, which provides empirical motion data, including orientation, linear velocity, and angular velocity. A numerical simulation of the dynamic model is executed, and the results compared to those obtained from the hardware RFV. The comparison reveals that the simulated model captures the prominent nonlinear dynamic behavior of the hardware RFV. Beyond providing a dynamic model suitable for model-based control system development, inspection of the model reveals insights for designing both the control system and the vehicle.

This paper is organized as follows: First, the necessary mathematical nomenclature for deriving a dynamic model of the RFV's terrestrial motion is presented. Second, the terrestrial dynamic model is derived yielding the governing equations of motion for rolling. Third, the flat outputs of the RFV are derived and compared to those of a CFMV. Fourth, a hardware embodiment of the RFV is briefly described. Fifth, simulation and experimental results are presented and discussed. Finally, various ways in which the dynamic model provides guidance for control system and vehicle design are discussed.

## 2 Methods

### 2.1 Mathematical Nomenclature.

This section briefly describes the mathematical conventions and nomenclature used to derive the terrestrial dynamic model of the RFV. Lowercase, uppercase, and double struck sub- and superscripts refer to points, bodies, and coordinate reference frames, respectively. The position vector $r⇀abC$ represents the position of point $b$ relative to point $a$ resolved (i.e., expressed) in the $C$ coordinate frame of reference. $rab,pC$ indexes the element of $r⇀abC$ in the $p$-direction, where $p$ can be $x, y$, or $z$. The velocity vector $d/dt(r⇀abC)D$ is the derivative with respect to (i.e., measured by an observer fixed in) reference frame $D$ of the position of point $b$ relative to point $a$ resolved in coordinate frame $C$. The acceleration vector $d2/dt2(r⇀abC)DE$ is the derivative with respect to reference frame $E$ of the derivative with respect to reference frame $D$ of the position of point $b$ relative to point $a$ resolved in coordinate frame $C$. The angular velocity $ω⇀ABC$ is the angular velocity of frame $B$ with respect to frame $A$ resolved in frame $C$. [ $a⇀$]x represents the skew-symmetric matrix formed from the elements of the vector $a⇀$ such that [ $a⇀$]x$b⇀=a⇀×b⇀$. $CAB$ is a rotation matrix that transforms a vector's frame of resolution from $A$ to $B$ such that $r⇀B=CABr⇀A$. $i,̂ĵ$, and $k̂$ represent the orthonormal basis vectors of $R3$ so that $î=[1 0 0]T, ĵ=[0 1 0]T$, and $k̂=[0 0 1]T$. All other variables are defined in the Nomenclature section.

### 2.2 Model Formulation.

Figure 2 schematically illustrates the RFV upon a horizontal plane representing the ground. Point $b$ is located at the center of the RFV body $B$, midway between two wheels $W1$ and $W2$, and is the origin of the body reference frame $B$, which is rigidly fixed in $B$. The axes of reference frame $B$ point along the principal directions of $B$. Point $i$ is located at the origin of an inertially fixed reference frame $I$. The rotation matrix from frame $B$ to frame $I$ is
$CBI=[ cos σ−sin σ cos α sin σ sin α sin σ cos σ cos α−cos σ sin α0 sin α cos α]$
(1)
Fig. 2
Fig. 2
Referring to Fig. 2(c), an intermediate reference frame $F$, is defined with respect to $I$ as
$CFI=[ cos σ−sin σ0 sin σ cos σ0001]$
(2)

The origin of frame $F$ is collocated with the origin of frame $B$ at point $b$. As a result, the x-axes of $F$ and $B$ are parallel, and $F$ is free to rotate about the x-axis of $B$. Reference frames $W1$ and $W2$ are located at the wheel center points $w1$ and $w2$ and are rigidly fixed in wheels $W1$ and $W2$, respectively. $w1$ and $w2$ are assumed to coincide with the center of mass of $W1$ and $W2$, respectively. Wheels $W1$ and $W2$ contact the ground at points $c1$ and $c2$, respectively. $W1$ and $W2$ are free to rotate about the x-axis of $B$.

Referring to Fig. 2(a), the RFV is modeled using six generalized coordinates
$q⇀=[xyσαθ1θ2]T$
(3)
$x$ and $y$ are the planar components of the position of point $b$ relative to $i$ and resolved in $I$. That is, $r⇀ib I=[x y R]T$, where $R$ is the radius of $W1$ and $W2$. $σ$ and $α$ are Euler angles specifying the yaw angle and pitch angle, respectively, of the $B$ frame with respect to the $I$ frame. $θ1$ and $θ2$ are the rotational positions of frames $W1$ and $W2$, respectively, with respect to $B$. The generalized velocity is
$q⇀̇=[ẋẏσ̇α̇θ̇1θ̇2]T$
(4)
where the $ẋ$ and $ẏ$ are the time derivatives of $x$ and $y$ with respect to the inertial frame of reference. That is, $d/dt(r⇀ib I)I=[ẋ ẏ 0]T$. The parameters that specify the system are the wheel radius $R$, the distance from point $b$ to the center of each wheel $L$, body mass $mB$, body principal moments of inertia $Ix, Iy, Iz$, wheel mass $mW$, wheel principal moments of inertia $Ixw, Iyw$, rolling resistance coefficient $μrr$, wheel viscous damping coefficient $μvisc$, aerodynamic drag coefficient $CD$, and the location of the center of mass $G$ (which is fixed in $B$) relative to point $b$ resolved in $B$, $r⇀bGB$. Referring to Fig. 2(b), the input to the system is a four-element vector $F⇀B$ containing a force and three moments resolved in the body frame of reference $B$
$F⇀B=[FzBMxBMyBMzB]T$
(5)
where $FzB$ is the force in the $B$ frame z-direction, and $MxB, MyB$, and $MzB$ are moments about the $B$ frame x-, y-, and z-directions, respectively. The angular velocity of frame $B$ with respect to frame $I$ resolved in frame $B$ is related to $σ̇$ and $α̇$ by the kinematic relationship
$ω⇀IBB=[01 sin α0 cos α0][σ̇α̇]=Cψω⇀IBB[σ̇α̇]$
(6)

where $Cψω⇀IBB$ is defined by Eq. (6). Note that $Cψω⇀IBB$ is not an orthonormal rotation matrix, but only relates the angular velocity of $B$ to the Euler angle rates.

### 2.3 Analytical Mechanics.

The Euler–Lagrange method of analytical mechanics is used to determine the equations of motion of the RFV in Fig. 2. Application of Hamilton's principle to nonconservative and nonholonomic systems results in Ref. [26]
$ddt(∂T∂q⇀̇)−∂T∂q⇀=(∂f⇀∂q⇀)Tλ⇀+Q⇀$
(7)
where $T$ is the total kinetic energy of the system, $∂f⇀/∂q⇀$ describes the differential constraints on the system be they holonomic or nonholonomic, $λ⇀$ is a vector of Lagrange multipliers, each element of which is associated with one constraint, and $Q⇀$ is the generalized force containing both monogenic and polygenic exogenous forces. The explicit formulation of the Lagrangian as kinetic energy minus potential energy is omitted because polygenic forces (which are not derivable from a scalar potential function) are involved. Gravity is the only conservative force in the system, and its effect is captured in the formulation of the generalized force $Q⇀$. This leads to a clearer interpretation of Eq. (7) wherein all inertial forces are on the left-hand side and all impressed forces, including forces of constraint, are on the right-hand side. The kinetic energy associated of the RFV body $B$ is
$TB=12(mBddt(r⇀iGB)ITddt(r⇀iGB)I+ω⇀IBBTIGBBω⇀IBB)$
(8)
where $IGBB$ is the inertia tensor of body $B$ about the center of mass $G$ resolved in coordinate frame $B$. The kinetic energy is a scalar quantity, and thus independent of the frame of resolution. Frame $B$ is chosen for convenience because the inertia tensor $IGBB$ is constant in $B$. $IGBB$ can be expressed in terms of the inertia tensor of $B$ about point $b$ resolved in frame $B$ using the parallel axis theorem for tensors
$IGBB=IbBB+mB[r⇀bGB]x2$
(9)
where $IbBB$ is a diagonal matrix containing the principal moments of inertia, $Ix, Iy,$ and $Iz$. Likewise, the velocity of point $G$ can be expressed in terms of the velocity of point $b$ using the transport theorem
$ddt(r⇀iGB)I=ddt(r⇀ibB)I+ddt(r⇀bGB)I=ddt(r⇀ibB)I+ω⇀IBB×r⇀bGB$
(10)
The linear and angular velocities in Eq. (10) can be expressed in terms of the generalized velocity as
$ddt(r⇀ibB)I=Jib,DBq⇀̇ω⇀IBB=JIB,RBq⇀̇$
(11)
where
$Jib,DB=[CIB[îĵ]03x203x2]JIB,RB=[03x2Cψω⇀IBB03x2]$
(12)
are Jacobians whose subscripts $D$ or $R$ indicate their representation of either linear or angular velocity, respectively. The subscript specifying the frame with respect to which differentiation occurs has been omitted in the Jacobian notation for readability, and it is assumed that all linear velocities are measured with respect to the $I$ frame. Postmultiplying $CIB$ by $[î ĵ]$ in Eq. (12) selects the first two columns corresponding to the $x$- and $y$-directions of $I$. Rewriting Eq. (10) in terms of the Jacobians defined in Eq. (12) yields
$ddt(r⇀iGB)I=(Jib,DB−[r⇀bGB]xJIB,RB)q⇀̇=JiG,DBq⇀̇$
(13)
Substituting Eqs. (9) and (13) into Eq. (8) yields
$TB=12q⇀̇TMBq⇀̇$
(14)
where
$MB=mBJiG,DBTJiG,DB+JIB,RBT(IbBB+mB[r⇀bGB]x2)JIB,RB$
(15)
is the symmetric mass matrix of the body $B$. Similarly, the mass matrix of the pth wheel can be obtained as
$MWp=mwJiwp,DBTJiwp,DB+JIWp,RBTIwWBJIWp,RB$
(16)
where
$Jiwp,DB=Jib,DB−[r⇀bWpB]xJIB,RBJIWp,RB=JIB,RB+JBWp,RBJBW1,RB=[03x4î03x1], r⇀bw1B=[L00]TJBW2,RB=[03x403x1î], r⇀bw2B=[−L00]T$
(17)
and $IwWB$ is the inertia tensor of the pth wheel about its center point $wp$ resolved in frame $B$. $IwWB$ is constant despite being resolved in frame $B$ because the wheel is assumed to be symmetric about its $x$-axis (Fig. 2(c)). Therefore $IwWB$ is a diagonal matrix whose elements are ${Ixw, Iyw, Iyw}$. The kinetic energy of each wheel is calculated in the same way as Eq. (14); therefore, the total kinetic energy of the RFV is
$T=12q⇀̇TMq⇀̇$
(18)
where $M=MB+MW1+MW2$. If $r⇀bGB = [0 0 h]T$, indicating that the center of mass of $B$ is centered between the wheels and located a distance $h$ from point $b$ in the $B$ frame z-direction, then
$M=[mB+2mW0mBh cos σ sin αmBh sin σ cos α000mB+2mWmBh sin σ sin α−mBh cos σ cos α00mBh cos σ sin αmBh sin σ sin α2(Iyw+mWL2)+Iy sin2α+Iz cos2α000mBh sin σ cos α−mBh cos σ cos α0Ix+2IxwIxwIxw000IxwIxw0000Ixw0Ixw]$
(19)
Considering the form of Eq. (18) and the fact that $M$ is symmetric, the left-hand side of the Euler–Lagrange equation (Eq. (7)) can be expanded for each of the six generalized coordinates as
$ddt(∂T∂q̇n)−∂T∂qn=∑m=16(Mnmq̈n+∑l=16(dMnmdql−12dMlmdqn)q̇lq̇m)$
(20)

where $n=1,2,…,6$.

### 2.4 Constraints.

The right-hand side of Eq. (7) consists of the impressed forces acting on the system. Forces of constraint are represented by $(∂f⇀/∂q⇀)Tλ⇀$. Referring to Fig. 2(b), each wheel is subjected to three constraint forces: A normal force that prevents the wheel from penetrating the ground, a lateral force that prevents slipping in the $F$ frame x-direction, and a friction force that permits rolling without slipping in the $F$ frame y-direction. Wheel slip is not modeled because the wheels are towed, rather than driven, reducing the likelihood of slip. All three constraints can be expressed succinctly by stating that the instantaneous velocity (with respect to $I$) of the contact point between the pth wheel and the ground is zero. That is
$ddt(r⇀icpB)I=0⇀$
(21)
Equation (21) is resolved in the $B$ frame simply to permit the use of the previously defined Jacobians in Eqs. (12) and (17), which are also resolved in $B$. However, Eq. (21) is valid when resolved in any coordinate frame because the zero vector $0⇀=[0 0 0]T$ is invariant under rotation. The position of point $cp$ relative to point $i$ is
$r⇀icpB=r⇀ibB+r⇀bwpB+CIBr⇀wpcpI$
(22)
where $r⇀wpcpI=[00−R]T$. Taking the time derivative of Eq. (22) with respect to the $I$ yields
$ddt(r⇀icpB)I=ddt(r⇀ibB)I+ω⇀IBB×r⇀bwpB+(ω⇀IBB+ω⇀BWpB)×CIBr⇀wpcpI=ddt(r⇀ibB)I+ω⇀IBB×(r⇀bwpB+CIBr⇀wpcpI)+ω⇀BWpB×CIBr⇀wpcpI=(Jib,DB−[r⇀bwpB+CIBr⇀wpcpI]xJIB,RB−[CIBr⇀wpcpI]xJBWp,RB)q⇀̇=Jicp,DBq⇀̇$
(23)
By definition Eq. (23) is equal to Eq. (21) and so $Jicp,DBq⇀̇=0⇀$. Also, $Jicp,DB$ can be resolved in the $F$ frame, which will result in a simpler closed-form expression when expanded. Doing so yields
$Jicp,DF=CBFJicp,DB=CBF(Jib,DB−[r⇀bwpB+CIBr⇀wpcpI]xJIB,RB−[CIBr⇀wpcpI]xJBWp,RB)$
(24)
Evaluating Eq. (24) for each wheel yields
$Jic1,DF=[ cos σ sin σ0000−sin σ cos σLRR0000000], Jic2,DF=[ cos σ sin σ0000−sin σ cos σ−LR0R000000]$
(25)
Therefore, the constraints are summarized by $Jic1,DFq⇀̇=0⇀$ and $Jic2,DFq⇀̇=0⇀$. The three rows of $Jicp,DFq⇀̇$ represent the three constraints enforced by the $pth$ wheel's reaction with the ground. These are kinematic constraints in that they are relations between the different elements of the generalized velocity $q⇀̇$, and they are also Pfaffian constraints in that they are linear in the generalized velocity [27]. The first row describes the lateral no-slip constraint in the $F$ frame x-direction. The first row of both $Jic1,DF$ and $Jic2,DF$ are identical, indicating that the conditions under which $W1$ and $W2$ experience lateral slip are identical. From Fig. 2(b), it is apparent that if one wheel were to slip laterally, the other would have to slip also. Expanding this constraint reveals that
$ẋ cos σ+ẏ sin σ=0$
(26)
Equation (26) is recognized to be the velocity of point $b$ relative to point $i$ in the x-direction of the $F$ frame. That is
$ddt(r⇀ibF)I=CIFddt(r⇀ibI)I=CIF[ẋẏ0]T=[ẋ cos σ+ẏ sin σ−ẋ sin σ+ẏ cos σ0]$
(27)

Therefore, the lateral no slip constraint from Eq. (26) is equivalent to saying that $d/dt(r⇀ib,xF)I=0$. This fact simplifies the equations of motion when resolved in the $F$ frame, as will be seen below. The second row of Eq. (25) describes the rolling-without-slipping constraint in the $F$ frame y-direction. The second rows of $Jic1,DF$ and $Jic2,DF$ are different because wheels $W1$ and $W2$ can rotate independently. The constraints represented by the first two rows of Eq. (25) are nonholonomic, indicating that the constraints cannot be integrated and expressed as relations between different elements of the generalized coordinates $q⇀$. Therefore, the nonholonomic constraints do not permit reducing the dimension of $q⇀$ via substitution. Instead, the method of Lagrange multipliers is used to enforce the constraints. The third row of Eq. (25) describes the holonomic normal constraint in the $F$ frame z-direction. This appears to be a trivial constraint because the velocity of point $cp$ (where $p=1$ or $2$) in the $F$ frame z-direction is not described by any combination of the elements of the generalized velocity. The triviality of the constraint is due to the problem statement; by definition, the wheels are constrained to only operate on the horizontal x–y plane of $I$, so it is known a priori that the z-position in the $F$ (or $I$) frame of any point on the system will be holonomically constrained by the other generalized coordinates. Therefore, excluding the z-position from the generalized coordinates is in effect an elimination by substitution, and the zeros in the third row of Eq. (25) reflect that decision. Excluding the z-position from the generalized coordinates is beneficial because it reduces the dimension of the configuration space by one. Alternatively, the generalized coordinates can include a z-position, which ultimately yields an equation of motion of $z̈=0$ while the other equations of motion remain unchanged. Additionally, inclusion of the z-position makes the normal constraint nontrivial and permits calculation of the total normal force on the vehicle. This is useful when calculating the moment due to rolling resistance, described below.

Between both wheels, a total of three unique and nontrivial constraints exist: No lateral slipping, wheel $W1$ rolls without slipping, and wheel $W2$ rolls without slipping. These three constraints are combined into a single differential constraint matrix
$∂f⇀∂q⇀=[ cos σ sin σ0000−sin σ cos σLRR0−sin σ cos σ−LR0R]$
(28)
where it is understood that $∂f⇀/∂q⇀$ cannot be integrated to yield $f⇀(q⇀)$ because the constraints are nonholonomic. From Eq. (7), the generalized forces of constraint are furnished through the use of Lagrange multipliers, $λ⇀$
$F⇀constraint=(∂f⇀∂q⇀)Tλ⇀$
(29)
where $λ⇀$ contains three elements associated with the three constraints. The introduction of $λ⇀$ means that the six equations in Eq. (7) have nine unknowns: The six elements of $q⇀̈$, and the three elements of $λ⇀$. Three auxiliary equations are furnished by time-differentiating the equations of constraint. Because $(∂f⇀/∂q⇀)q⇀̇=0⇀$
$ddt(∂f⇀∂q⇀q⇀̇)=∂f⇀∂q⇀q⇀̈+q⇀̇T∂∂q⇀(∂f⇀∂q⇀)q⇀̇=0⇀$
(30)
Substituting Eq. (28) into Eq. (30) and differentiating yields the three auxiliary equations
$ẍ cos σ−ẋσ̇ sin σ+ÿ sin σ+ẏσ̇ cos σ=0−ẍ sin σ−ẋσ̇ cos σ+ÿ cos σ−ẏσ̇ sin σ+Lσ̈+R(α̈+θ̈1)=0−ẍ sin σ−ẋσ̇ cos σ+ÿ cos σ−ẏσ̇ sin σ−Lσ̈+R(α̈+θ̈2)=0$
(31)

Care must be taken when differentiating Eq. (30) because the differentiation of any position or velocity elements will implicitly be done with respect to the frame of resolution [26]. In this case, $x, ẋ, y$, and $ẏ$ are all resolved in $I$, and $ẋ$ and $ẏ$ differentiations are with respect to $I$, (see Eq. (3) and subsequent discussion) so the frame of resolution matches the reference frame and the second derivatives $ẍ$ and $ÿ$ will be with respect to $I$ as well.

### 2.5 Impressed Forces.

The second term on the right-hand side of Eq. (7) is the generalized force, $Q⇀$. $Q⇀$ is the superposition of five groups of exogenous impressed forces and moments acting on the system
$Q⇀=Q⇀grav+Q⇀control+Q⇀rr+Q⇀visc+Q⇀drag$
(32)
where $Q⇀grav$ are gravitational forces, $Q⇀control$ are control forces that can be generated by the systems actuators, $Q⇀rr$ are rolling resistance forces, $Q⇀visc$ are viscous damping forces, and $Q⇀drag$ are aerodynamic drag forces. In general, the virtual work $δU$ performed on a body by an exogenous force $F⇀$ and moment $M⇀$ is $δU=δq⇀T(JDTF⇀+JRTM⇀)$, where $δq⇀$ is the virtual displacement of the body in the configuration space, and $JD$ and $JR$ are the displacement and angular Jacobians, respectively. Therefore, the generalized force associated with $F⇀$ and $M⇀$ is $Q⇀=JDTF⇀+JRTM⇀$. The generalized force due to gravity is
$Q⇀grav=(mB(JiG,DB)T+mW(Jiw1,DB)T+mW(Jiw2,DB)T)CIBg⇀I$
(33)
where $g⇀I=−gk̂$, and $g$ is the gravitational acceleration. The generalized force due to the control forces and moments from Eq. (5) is
$Q⇀control=[Jib,DBJIB,RB]T[00FzBMxBMyBMzB]T$
(34)
The generalized force due to rolling resistance is
$Q⇀rr=(JIW1,RB)TM⇀rr1B+(JIW2,RB)TM⇀rr2B$
(35)
where $M⇀rrpB$ is a rolling resistance moment that acts to oppose to rotation of the pth wheel. If the rolling resistance is proportional to the normal force on the wheel then
$M⇀rrpB=Mrrpî, Mrrp=−μrrRFNp sgn(îTω⇀IWpB)=−μrrRFNp sgn(îTJIWp,RBq⇀̇)$
(36)
where $FNp>0$ represents the normal force on the pth wheel, and $μrr$ is a coefficient of rolling resistance that depends on the terrain. As mentioned previously, the value of $FNp$ is determined by augmenting the generalized coordinates with a $z$-motion component and computing the generalized force of constraint in the $I$ frame z-direction from Eq. (29), which is approximately (assuming equal force on each wheel)
$FNp≈12((mB+2mw)g−FzB cos α)$
(37)
The generalized force due to viscous damping is
$Q⇀visc=(JBW1,RB)TM⇀visc1B+(JBW2,RB)TM⇀visc2B$
(38)
where $M⇀viscpB$ is a viscous damping moment that acts to oppose to rotation of the pth wheel with respect to the body $B$. This moment arises due to bearing losses and is
$M⇀viscpB=Mviscpî, Mviscp=−μvisc (îTω⇀BWpB)=−μvisc (îTJBWp,RBq⇀̇)$
(39)
The generalized force due to aerodynamic drag on the body is
$Q⇀drag=(Jib,DB)TCFBF⇀dragF$
(40)
where $F⇀dragF$ is an aerodynamic drag force that acts on and opposes the direction of motion of point $b$
$F⇀dragF=Fdragĵ, Fdrag=−CDvy|vy|, vy=ĵTCBFJib,DBq⇀̇$
(41)
Substituting the foregoing expressions for the constituent generalized forces into Eq. (32) yields
$Q⇀=[(FzB sin α−Fdrag)sin σ−(FzB sin α−Fdrag)cos σMzB cos α+MyB sin αMrr1+Mrr1+MxB+mBgh sin αMrr1+Mvisc1Mrr2+Mvisc2]$
(42)

### 2.6 Equations of Motion.

The equations of motion of Eq. (7) can now be obtained by combining Eqs. (20), (28), (31), and (42) and solving for $q⇀̈$ and $λ⇀$. As mentioned before, $ẍ$ and $ÿ$ are most easily expressed in the $F$ frame. Differentiating Eq. (27) with respect to the $I$ frame yields $d2/dt2(r⇀ibF)II=CIF[ẍÿ0]T$, which is the acceleration of point $b$ resolved in $F$. If $h=0$ (i.e., $b$ and $G$ are coincident), then the equations of motion are
$d2dt2(rib,xF)II=−σ̇ddt(rib,yF)I$
(43)
$d2dt2(rib,yF)II=−R2FzB sin α+R(Mrr1+Mrr2+Mvisc1+Mvisc2−RFdrag)2Ixw+R2(mB+2mW)$
(44)
$σ̈=R2(MzB cos α+MyB sin α+α̇σ̇(Iz−Iy)sin(2α))+RL(Mrr2−Mrr1+Mvisc2−Mvisc1)2L2Ixw+R2(2(Iyw+mwL2)+Iy sin2α+Iz cos2α)$
(45)
$α̈=MxB−Mvisc2−Mvisc1−12(Iz−Iy)σ̇2 sin(2α)Ix$
(46)
The direct dependence of $α̈, σ̈$, and $d2/dt2(r⇀ibF)II$ on $θ̈1$ and $θ̈2$ has been eliminated by the constraints from Eq. (28), though $θ̇1$ and $θ̇2$ do indirectly appear in the rolling resistance and viscous damping moments, defined in Eqs. (36) and (39), respectively. Despite the lateral no-slip constraint from Eq. (28), which indicates that $d/dt(rib,xF)I=0$, Eq. (43) indicates that the acceleration in the $F$ frame x-direction is nonzero. In fact, Eq. (43) is recognized as the centripetal acceleration. This fact allows the acceleration to be simplified even further by computing the derivative of velocity with respect to the $F$ frame rather than the $I$ frame. To do this, the transport theorem is applied to Eqs. (43) and (44)
$d2dt2(r⇀ibF)IF=d2dt2(r⇀ibF)II+ω⇀FIF×ddt(r⇀ibF)I=d2dt2(r⇀ibF)II−σ̇k̂×(0î+ddt(rib,yF)Iĵ)=d2dt2(r⇀ibF)II+σ̇ddt(rib,yF)Iî=d2dt2(rib,yF)IIĵ$
(47)
Equation (47) indicates that the acceleration of point $b$ with respect to an observer fixed in $F$ is tangent to the velocity of point $b$, and its value is given exactly by Eq. (44). That is, Eq. (47) is the acceleration of $b$along the path on which $b$ moves, which is also constrained to be instantaneously in the y-direction of $F$ because of the lateral no-slip constraint. The change of reference frames performed in Eq. (47) essentially eliminates the x-equation of motion, supporting the intuition that the RFV, characterized by six generalized coordinates and subject to three nonholonomic constraints, will have $6−3=3$ degrees-of-freedom. The nonholonomic constraints do not reduce the number of generalized coordinates because it is still possible to achieve any configuration in the six-dimensional configuration space [27]. However, the motion of the system is constrained to lie on a three-dimensional subspace of the six-dimensional configuration space. The evolution of the RFV on this subspace is given by Eqs. (44)(46). The equations of motion can be further simplified by defining a new set of control inputs that are resolved in the $F$ frame
$M⇀F=[MxFMyFMzF]=CBF[MxBMyBMzB]=[MxBMyB cos α−MzB sin αMyB sin α+MzB cos α]F⇀F=[FxFFyFFzF]=CBF[00FzB]=[0−FzB sin αFzB cos α]$
(48)
Rewriting Eqs. (44), (45), and (46) in terms of $FyF,MxF, MyF$, and $MzF$ and defining the $F$ frame velocity and acceleration as $ẏF≡d/dt(rib,yF)I$ and $ÿF≡d2/dt2(rib,yF)II$, respectively, yields
$ÿF=FyF−Dy2IxwR2+mB+2mW$
(49)
$σ̈=MzF+α̇σ̇(Iz−Iy)sin(2α)−Dσ2(L2IxwR2+Iyw+mwL2)+Iy sin2α+Iz cos2α$
(50)
$α̈=MxF−12(Iz−Iy)σ̇2 sin(2α)−DαIx$
(51)
where
$Dy=1R(Mrr1+Mrr2+Mvisc1+Mvisc2)−FdragDσ=LR(Mrr1−Mrr2+Mvisc1−Mvisc2)Dα=Mvisc1+Mvisc2$
(52)
are dissipation terms, which are all functions of $q⇀$ and $q⇀̇$, in general. Each equation of motion is now grouped into control terms, inertial coupling terms (except for $ÿF$, which has no inertial coupling) and dissipation terms. Finally, the assumption that $h=0$ can be relaxed; if
$|h|≪Ix(mBR2+2(mWR2+Ixw))mBR$
(53)
then the equations of motion become
$ÿF=FyF−Dy+mBhIx(12mBgh sin(2α)−(Dα−MxF)cos α−(Ix(α̇2+σ̇2)+(Iz−Iy)σ̇2 cos2α)sin α)2IxwR2+mB+2mW$
(54)
$σ̈=MzF+α̇σ̇(Iz−Iy)sin(2α)−Dσ+mBhσ̇ẏF sin α2(L2IxwR2+Iyw+mwL2)+Iy sin2α+Iz cos2α$
(55)
$α̈=MxF−12(Iz−Iy)σ̇2 sin(2α)−Dα+mBgh sin αIx−mBhIx(12mBh(α̇2+σ̇2)sin(2α)+(FyF−Dy)cos α)2IxwR2+mB+2mW$
(56)

## 3 Differential Flatness

A differentially flat system is one wherein the system's states and inputs can be written in terms of the system's outputs and their derivatives [22]. If these conditions are satisfied, then the outputs are referred to as flat outputs. A consequence of being differentially flat is that any trajectory specified in terms of the flat outputs can be directly mapped to the inputs, and thus used to specify feasible motion trajectories. It can be shown (see, e.g., Refs. [28] and [29]) that the flat outputs of a CFMV are the three components of linear velocity in the $I$ frame and $σ$; $α$ is absent, indicating that it cannot be used to specify motion trajectories for the CFMV. On the other hand, if the RFV's state and input vectors are defined as $[σ σ̇ α α̇ ẏF]T$ and $[MzF MxF FyF ]T$, respectively, then inspection of the equations of motion Eqs. (54)(56) reveals that $z⇀=[ σ α ẏF ]T$ are flat outputs because
$FyF=FyF(α,α̇,σ̇,ÿF)=FyF(z⇀,z⇀̇,z⇀̈)MxF=MxF(α,α̇,σ̇,α̈,FyF(α,α̇,σ̇,ÿF))=MxF(α,α̇,σ̇,α̈,ÿF)=MxF(z⇀,z⇀̇,z⇀̈)MzF=MzF(α,α̇,σ̇,ẏF,σ̈)=MzF(z⇀,z⇀̇,z⇀̈)$
(57)

The differences in flat outputs are due to the RFV's lateral no-slip constraint (Eq. (26)) and the fact that the normal constraints on the wheels result in normal forces that partially support the vehicle's weight. The result is that a given horizontal force $FyF$ can be prescribed independently of $α$ (except when $α=0, π$) because the vertical component of the net thrust need not completely offset the vehicle weight. This contrasts with a CFMV, whose vertical component of the net thrust must always equal the vehicle's weight to maintain constant altitude.

## 4 Hardware Embodiment

Figure 3 illustrates a hardware embodiment of the RFV. Attached to each wheel is a quadrature encoder (U.S. Digital E3 series, Vancouver, WA) used to resolve the angular position and angular rate of the wheel. A strapdown inertial navigation system (INS) (VectorNav Technologies VN-100, Dallas, TX) measures the orientation of the RFV. Data from the encoders and INS are obtained and processed by an on-board ARM cortex-M7 microcontroller (ST Microelectronics STM32F767ZIT6U, Geneva, Switzerland). A 2.4 GHz radio transceiver (Digi XBee Pro-S2C, Hopkins, MN) transmits motion data to a laptop computer every 2 ms. The data transmitted include $σ, α, θ1, θ2$, and $ω⇀IBB$.

Fig. 3
Fig. 3
Although the RFV model derived in Sec. 2 assumes operation on a horizontal plane, the hardware embodiment of the RFV must account for nonhorizontal ground in order to accurately compute $σ̇$ and $α̇$. The angle of the ground is reflected in the lateral rolling angle $γ$ about the $F$ frame y-axis. Measurements of $σ, γ,α$, and $ω⇀IBB$ are obtained from the INS. $σ̇$ and $α̇$ are computed from $ω⇀IBB$ according to
$[σ̇γ̇α̇]=[0 sin α sec γ cos α sec γ0 cos α−sin α1 sin α tan γ cos α tan γ]ω⇀IBB$
(58)
$ẏF$ is computed from $θ̇1$ and $θ̇2$ (measured by the encoders) and $ω⇀IBB$ as follows: From Eq. (23), the velocity of the body can be expressed as a function of the angular velocity of pth wheel
$ddt(r⇀ibF)I=CBF(−ω⇀IBB×r⇀bwpB−(ω⇀IBB+ω⇀BWpB)×CFBr⇀wpcpF)$
(59)
where $r⇀wpcpF=[00−R]T$. This results in two equivalent expressions for $d/dt(r⇀ibF)I$ (one for each wheel), which can be used to solve for $ẏF$ as a function of $θ̇1$ and $θ̇2$. Solving Eq. (58) for $ω⇀IBB$, substituting the result into Eq. (59), and then expanding the y-component of Eq. (59) for each wheel yields
$y˙F=ddt(rib,yF)I=−R(θ˙1+α˙)+σ˙(R sin γ−L cos γ)=−R(θ˙2+α˙)+σ˙(R sin γ+L cos γ)$
(60)
Adding the two expressions in Eq. (60), substituting $σ̇$ and $α̇$ from Eq. (58) and solving for $ẏF$ yields
$ẏF=−R(12(θ̇1+θ̇2)+ωIB,xB)$
(61)

The result is independent of $γ$, which makes sense considering that $F$ frame is defined with respect to the ground, i.e., the $F$ frame rotates with $γ$.

Additionally, an estimate of $σ̇$ is obtained using odometry data from the wheel encoders; equating the two expressions from Eq. (60) and substituting in $σ̇$ from Eq. (58) yields
$σ̇odom=R2L(θ̇2−θ̇1)sec γ$
(62)

Comparing $σ̇odom$ and $σ̇$ from Eq. (58) provides a method for determining if the RFV's wheels are in contact with the ground, assuming the ground is horizontal.

## 5 Model Validation

To validate the dynamic model developed in Sec. 2, empirical motion data from the hardware RFV are collected and compared to motion data from an RFV simulation. These experiments are carried out with the control inputs $FyF,MxF, MyF$, and $MzF$ all set to zero, i.e., they are open loop. The goal is to validate the nonlinear dynamics of the system in response to nonzero initial conditions. For the empirical experiments, this is accomplished by manually providing an initial excitation with no control input, and then recording the subsequent motion. Then, the initial conditions of the simulations are set to match those of the empirical data, and the two sets of data are compared.

### 5.1 Numerical Simulation Methods.

While precise knowledge of the model parameters is not required to compare qualitative characteristics of the RFV's dynamic behavior, the simulation model parameters are based on those of the hardware RFV. Hardware RFV parameters that can be directly measured (e.g., mass, length, center of mass) are so obtained. Moments of inertia are calculated in a piecewise manner; inertia values for components designed using cad software are obtained therein, while other components are modeled with simple geometries and the moments of inertia are calculated accordingly. Dissipation parameters (e.g., rolling resistance, viscous damping), are estimated based on the decay of the measured $α$ and $σ$ signals. The model parameter values appear in Table 1 under the set 1 column.

Table 1

Model parameters for the hardware RFV and simulations

Value
ParameterSet 1Set 2
$L$ (m)0.1650.165
$R$ (m)0.22540.2254
$mB$ (kg)1.10.7
$mW$ (kg)0.089750.08975
$Ix$ (kg m2)0.003710.001829
$Iy$ (kg m2)0.003060.002191
$Iz$ (kg m2)0.003420.002567
$Ixw$(kg m2)0.0012010.001201
$Iyw$ (kg m2)0.0006030.000603
$g$ (m/s2)9.819.81
$h$ (m)−0.006−0.002
$μrr$0.0050.005
$μvisc$ (N m s)0.0010.001
$CD$ (kg/m)00
Value
ParameterSet 1Set 2
$L$ (m)0.1650.165
$R$ (m)0.22540.2254
$mB$ (kg)1.10.7
$mW$ (kg)0.089750.08975
$Ix$ (kg m2)0.003710.001829
$Iy$ (kg m2)0.003060.002191
$Iz$ (kg m2)0.003420.002567
$Ixw$(kg m2)0.0012010.001201
$Iyw$ (kg m2)0.0006030.000603
$g$ (m/s2)9.819.81
$h$ (m)−0.006−0.002
$μrr$0.0050.005
$μvisc$ (N m s)0.0010.001
$CD$ (kg/m)00
The simulation entails numerically integrating Eqs. (43), (54)(56) along with expressions for $θ̇1$ and $θ̇2,$ which are obtained by rearranging the constraint equations in the second and third rows of Eq. (28)
$θ̇1=1R(ẋ sin σ−ẏ cos σ−Lσ̇−Rα̇)=1R(−ẏF−Lσ̇−Rα̇)θ̇2=1R(ẋ sin σ−ẏ cos σ+Lσ̇−Rα̇)=1R(−ẏF+Lσ̇−Rα̇)$
(63)

Initial conditions are specified to match those of the empirical data. The simulation is executed using matlab's ode45 function, which is a fourth-order Runge–Kutta numerical integration function with a variable time-step.

### 5.2 Empirical Data Collection.

Empirical motion data are collected from the hardware RFV after manually applying an external moment about the $F$ frame z-axis and then letting the RFV come to rest. The external moment is intended to provide a nonzero initial $α, α̇$, and $σ̇$ and so excite the nonlinear dynamics of Eqs. (54)(56). The motion data consist of $σ, α, θ1, θ2$, and $ω⇀IBB$ sampled at 2 ms intervals. Values of $σ̇$ and $α̇$ are calculated from the motion data according to $σ̇=ωIB,yBsinα+ωIB,zBcosα$ and $α̇=ωIB,xB$ (see Eq. (6)). $θ̇1$ and $θ̇2$ are obtained by numerically differentiating measurements of $θ1$ and $θ2$. Then $ẏF$ is determined from Eq. (61).

### 5.3 Comparing Empirical and Simulated Data.

In the following time series plots, the solid and dashed lines represent the simulated and empirical data, respectively. Figures 4(a) and 4(b) show $α$ and $α̇$, respectively, as functions of time, and indicate good agreement between the simulation and empirical data despite uncertainty in the parameter values. $α$ oscillates and decays, as expected, however the response is not a pure sinusoid, particularly at large values of $α$. This is due primarily to the $mBghsinα$ term in the numerator of Eq. (56), and to a lesser extent the $−1/2(Iz−Iy)σ̇2sin2α$ term; the remaining terms in Eq. (56) are more than an order of magnitude smaller than the $mBghsinα$ term. That $α$ is not a pure sinusoid is best observed by examining $α̇$ in Fig. 4(b), which exhibits sharp rather than smooth peaks when $α$ is large.

Fig. 4
Fig. 4

Figures 5(a) and 5(b) show $σ$ and $σ̇$, respectively, as functions of time. The empirical $σ̇$ data are zero-phase low-pass filtered to remove noise. The rate at which $σ̇$ decays is influenced primarily by $μvisc$ and $μrr$, with the latter causing the plateau in $σ$. The ripple in the empirical $σ̇$ data is caused by the oscillation in $α$, and is predicted by the $α̇σ̇(Iz−Iy)sin2α$ and $mBhσ̇ẏFsinα$ terms in Eq. (55); the same ripple can be seen in the simulated data.

Fig. 5
Fig. 5

The empirical and simulated values of the wheel angles $θ1$ and $θ2$ appear in Fig. 6(a). The simulated signals exhibit the same trends as the empirical signals; they oscillate while asymptotically approaching a terminal value. The oscillations in $θ1$ and $θ2$ coincide with the oscillation in $α$. Figure 6(b) illustrates $ẏF$, which is initially responsible for some of the ripple in $σ̇$ via the $mBhσ̇ẏFsinα$ term in Eq. (55), though this term decays quickly and is soon overtaken by the $α̇σ̇(Iz−Iy)sin2α$ term. Both the simulated and empirical values of $ẏF$ are calculated according to Eq. (61), which requires knowledge of $θ̇1$ and $θ̇2$, neither of which are acquired directly from the motion data. Therefore, the empirical value of $ẏF$ is obtained by using numerically differentiated values of the empirical $θ1$ and $θ2$ and zero-phase low-pass filtering the result. Despite the filtering, the numerical differentiation results in poor signal quality. However, the overall oscillation and decay of the signal closely match the simulated result.

Fig. 6
Fig. 6

### 5.4 Demonstrating Gyroscopic Nonlinearity.

The primary cross-coupling between the $α$ and $σ$ degrees-of-freedom is due to the “gyroscopic” term containing $Iz−Iy$ appearing in both Eqs. (55) and (56). For the purposes of demonstrating the nonlinear behavior of the RFV and validating the nonlinear model of Sec. 2.6, the RFV is explicitly configured such that $Iz−Iy≠0$. To further emphasize the gyroscopic nonlinearity, a second set of tests are executed wherein the center of mass is moved nearer to point $b$ so that $h=−0.002$ m. This results in the $Iz−Iy$ terms having a greater influence on the dynamic behavior of the system because the other terms are attenuated when $|h|$ is small. The model parameters for this test appear in the set 2 column of Table 1. Empirical and simulated data for $α$ and $α̇$ appear in Figs. 7(a) and 7(b), respectively. Initially, $α$ is near $π$ radians, indicating that the RFV chassis is upside down. Since $h<0$, this configuration is statically unstable due to the gravitational term $mBghsinα$ in Eq. (56). However, rather than immediately diverging toward the stable equilibrium at $α=0$, $α$ instead oscillates about $π$ radians for several seconds (see orange oval in Fig. 7(a)) due to the gyroscopic term $−1/2(Iz−Iy)σ̇2sin2α$. This “gyroscopically” stabilized behavior continues until of the amplitude of the gyroscopic term decreases (as $σ̇$ decreases) and the gravitational term $mBghsinα$ begins to dominate. That is, $α=π$ behaves like a stable equilibrium of Eq. (56) when $σ̇$ is viewed as a parameter, with the equilibrium undergoing a bifurcation as $σ̇$ drops below a critical value. The simulation predicts the same behavior, with $α$ briefly oscillating about $α=π$ and then diverging toward $α=2π (=0)$.

Fig. 7
Fig. 7

The numerical simulation results demonstrate the same qualitative behavior in $σ, α, θ1, θ2$, and $ẏF$ as the empirical data obtained from the hardware embodiment of the RFV. The quantitative differences between the two (e.g., decay time, number of oscillations, absolute magnitude) are due to uncertain parameter values. In particular, the moments of inertia and the coefficients of viscous and rolling friction are uncertain.

## 6 Conclusions and Future Work

The class of bi-modal rolling-flying vehicle presented in this paper offers potential energetic benefits compared to CFMVs. However, the control algorithms developed for CFMVs are not ideal for controlling the terrestrial motion of the RFV because of fundamental differences in the system dynamics (Sec. 2) and flat outputs (Sec. 3).

### 6.1 Model Use in Control System Design.

This paper presents an accurate dynamic model for RFV terrestrial motion, providing the foundation for nonlinear, model-based control system development. Furthermore, inspection of the model and simulation results provides some insights for the control system designer. First, the equations of motion (54)(56) reveal the appropriate outputs for which to specify motion trajectories, i.e., the flat outputs $σ$, $α$, and $ẏF$ (Sec. 3). Second, the quantitative differences between the experimental and simulated data indicate sensitivity to parameter uncertainty, including difficult-to-measure inertial and dissipative parameters. That is, while the structures of the model's nonlinearities are known (e.g., $−1/2(Iz−Iy)σ̇2sin2α$ from Eq. (56)), the precise values of their coefficients are uncertain. This motivates selecting a nonlinear control technique that is robust to parameter uncertainty, such as sliding mode control. Such a technique requires uncertainty bounds, which the model developed herein can provide; the uncertainty in each parameter can be estimated; therefore, the coefficient for each term in the equations of motion can be bounded. Furthermore, the model reveals which terms of the equations of motion are dominant, providing the opportunity for informed model simplification while also providing bounds on the model uncertainty resulting from simplification. That is, dominant terms can be compensated for directly within a simplified model, while nondominant terms (e.g., the final terms in Eqs. (54)(56)) can be aggregated together as a bounded modeling uncertainty.

The dynamic model also motivates vehicle design considerations. For example, locating the center of mass near to the wheel axle (small $h$) reduces the impact of gravitational nonlinearities, thereby reducing control effort and simplifying control algorithms. The model also reveals that the primary coupling mechanism between the $α$ and $σ$ degrees-of-freedom is the “gyroscopic” term $Iz−Iy$ appearing in both Eqs. (55) and (56). This suggests that the coupling could be nearly eliminated by designing the RFV such that $Iz−Iy≈0$. Decoupling the degrees-of-freedom and removing nonlinearities would simplify the design of closed-loop controllers for $α$ and $σ$.

### 6.2 Future Work.

While the numerical simulation accurately predicts key qualitative behaviors captured in the empirical data, such as gravitational and gyroscopic nonlinearities, more precise parameter knowledge would improve the quantitative predictive power of the model. To that end, future work entails performing a rigorous system identification of the hardware RFV and developing robust, closed-loop control algorithms based on the dynamic model presented herein. Furthermore, the dynamic model could be improved by including aerodynamic effects, such as rotor thrust models [17] or ground effect [30].

## Acknowledgment

The authors would like to thank Tyler Jenkins, Ryan Lynch, and Ellis Stevens for their work in creating the hardware embodiment of the RFV.

## Funding Data

• U.S. Army Research Office and the U.S. Army Special Operations Command (Contract No. W911-NF-13-C-0045; Funder ID: 10.13039/100000183).

## Nomenclature

• $CD$ =

aerodynamic drag coefficient

•
• $Cψω⇀IBB$ =

Euler angle rates-to-angular velocity transformation matrix

•
• $d/dt(⋅)A$ =

time derivative of $(⋅)$ with respect to the $A$ reference frame

•
• $d2/dt2(⋅)AB$ =

time derivative of $d/dt(⋅)A$ with respect to the $B$ reference frame

•
• $Dy,Dα,Dσ$ =

dissipation terms

•
• $F⇀constraint$ =

generalized constraint force

•
• $Fdrag$ =

aerodynamic drag force (N)

•
• $FNp$ =

pth wheel normal force

•
• $FzB$ =

body force (N)

•
• $g$ =

gravitational acceleration (m/s2)

•
• $h$ =

location of RFV body center of mass (m)

•
• $IaBC$ =

inertia tensor of body $B$ about point $a$ resolved in reference frame $C$

•
• $Ix,Iy,Iz$ =

RFV body principal moments of inertia of (kg m2)

•
• $Ixw,Iyw$ =

wheel principal moments of inertia of (kg m2)

•
• $Jab,DC$ =

translational motion Jacobian

•
• $JAB,RC$ =

rotational motion Jacobian

•
• $L$ =

RFV characteristic length (m)

•
• $M,MB,Mwp$ =

mass matrices

•
• $mB,mW$ =

masses (kg)

•
• $Mrr,p$ =

pth wheel rolling resistance moment (N·m)

•
• $Mvisc,p$ =

pth wheel viscous damping moment (N·m)

•
• $MxB,MyB,MzB$ =

body moments (N·m)

•
• $q⇀$ =

generalized coordinates

•
• $Q⇀,Q⇀grav,Q⇀control,Q⇀rr,Q⇀visc,Q⇀drag$ =

generalized forces

•
• $R$ =

•
• $r⇀ab$ =

position of point $b$ relative to point $a$

•
• $T,TB,TW$ =

kinetic energies

•
• $x$ =

position coordinate (m)

•
• $y$ =

position coordinate (m)

•
• $ẏF,ÿF$ =

velocity and acceleration in the $F$ frame y-direction

•
• $z⇀,z⇀̇,z⇀̈$ =

flat outputs and derivatives

•
• $α$ =

•
• $γ$ =

tilt Euler angle

•
• $θ1,θ2$ =

•
• $λ⇀$ =

Lagrange multipliers

•
• $μvisc$ =

viscous damping coefficient (N·m s)

•
• $μrr$ =

rolling resistance coefficient

•
• $σ$ =

•
• $ω⇀AB$ =

angular velocity of reference frame $B$ with respect to reference frame $A$

•
• $∂f⇀/∂q⇀$ =

differential constraint matrix

### Superscripts

Superscripts

• $B$ =

variable resolved in the $B$ (body) reference frame

•
• $F$ =

variable resolved in the $F$ reference frame

•
• $I$ =

variable resolved in the $I$ (inertial) reference frame

## References

1.
Choi
,
S. H.
, and
Zhu
,
W. K.
,
2012
, “
Performance Optimisation of Mobile Robots for Search-and-Rescue
,”
Appl. Mech. Mater.
,
232
, pp.
403
407
.10.4028/www.scientific.net/AMM.232.403
2.
Murphy
,
R. R.
,
,
S.
, and
Kleiner
,
A.
,
2016
,
Disaster Robotics
, B. Siciliano and O. Khatib O, eds.,
Springer Handbook of Robotics, Springer
, Cham, Switzerland.
3.
Bishop
,
B.
,
Crabbe
,
F.
, and
Hudock
,
B.
,
2005
, “
Design of a Low-Cost, Highly Mobile Urban Search and Rescue Robot
,”
,
19
(
8
), pp.
1
27
.
4.
Neumann
,
M.
,
Predki
,
T.
,
Heckes
,
L.
, and
Labenda
,
P.
,
2013
, “
Snake-Like, Tracked, Mobile Robot With Active Flippers for Urban Search-and-Rescue Tasks
,”
Ind. Rob.
,
40
(
3
), pp.
246
250
.10.1108/01439911311309942
5.
Santos
,
J. M.
,
Krajník
,
T.
, and
Duckett
,
T.
,
2017
, “
Spatio-Temporal Exploration Strategies for Long-Term Autonomy of Mobile Robots
,”
Rob. Auton. Syst.
,
88
, pp.
116
126
.10.1016/j.robot.2016.11.016
6.
Kalantari
,
A.
, and
Spenko
,
M.
,
2014
, “
Modeling and Performance Assessment of the HyTAQ, a Hybrid Terrestrial/Aerial Quadrotor
,”
IEEE Trans. Rob.
,
30
(
5
), pp.
1278
1285
.10.1109/TRO.2014.2337555
7.
Takahashi
,
N.
,
Yamashita
,
S.
,
Sato
,
Y.
,
Kutsuna
,
Y.
, and
,
M.
,
2015
, “
All-Round Two-Wheeled Quadrotor Helicopters With Protect-Frames for Air-Land-Sea Vehicle (Controller Design and Automatic Charging Equipment
,”
,
29
(
1
), pp.
69
87
.10.1080/01691864.2014.991754
8.
Morton
,
S.
, and
Papanikolopoulos
,
N.
,
2017
, “
A Small Hybrid Ground-Air Vehicle Concept
,”
IEEE International Conference on Intelligent Robots and Systems
, Vancouver, BC, Sept. 24–28, pp. 5149–5154. 10.1109/IROS.2017.8206402
9.
Kossett
,
A.
, and
Papanikolopoulos
,
N.
,
2011
, “
A Robust Miniature Robot Design for Land/Air Hybrid Locomotion
,”
Proceedings of the IEEE International Conference on Robotics and Automation
, Shanghai, China, May 9–13, pp.
4595
4600
.10.1109/ICRA.2011.5979845
10.
Kawasaki
,
K.
,
Zhao
,
M.
,
,
K.
, and
Inaba
,
M.
,
2013
, “
MUWA: Multi-Field Universal Wheel for Air-Land Vehicle With Quad Variable-Pitch Propellers
,”
IEEE International Conference on Intelligent Robots and Systems
, Tokyo, Japan, Nov. 3–7, pp.
1880
1885
11.
Boria
,
F. J.
,
Bachmann
,
R. J.
,
Ifju
,
P. G.
,
Quinn
,
R. D.
,
Vaidyanathan
,
R.
,
Perry
,
C.
, and
Wagener
,
J.
,
2005
, “
A Sensor Platform Capable of Aerial and Terrestrial Locomotion
,”
IEEE/RSJ International Conference on Intelligent Robots and Systems
, Edmonton, AB, Aug. 2–6, pp.
3959
3964
.10.1109/IROS.2005.1545597
12.
Sreevishnu
,
S.
,
Koshy
,
M.
,
Krishnan
,
A.
, and
Das
,
G. P.
,
2018
, “
Kinematic Design, Analysis and Simulation of a Hybrid Robot With Terrain and Aerial Locomotion Capability
,”
MATEC Web of Conf.
, 172, p. 03008. 10.1051/matecconf/201817203008
13.
Stewart
,
W.
,
Weisler
,
W.
,
Macleod
,
M.
,
Powers
,
T.
,
Defreitas
,
A.
,
Gritter
,
R.
,
Anderson
,
M.
,
Peters
,
K.
,
Gopalarathnam
,
A.
, and
Bryant
,
M.
,
2018
, “
Design and Demonstration of a Seabird-Inspired Fixed-Wing Hybrid UAV-UUV System
,”
Bioinspiration Biomimetics
, 13(5).
14.
Siddall
,
R.
, and
Kovač
,
M.
,
2014
, “
Launching the AquaMAV: Bioinspired Design for Aerial-Aquatic Robotic Platforms
,”
Bioinspiration Biomimetics
,
9
(
3
).
15.
Maia
,
M. M.
,
,
D. A.
, and
Diez
,
F. J.
,
2017
, “
Design and Implementation of Multirotor Aerial-Underwater Vehicles With Experimental Results
,”
IEEE International Conference on Intelligent Robots and Systems (IROS)
, Vancouver, BC, Sept. 24–28, pp.
961
966
.10.1109/IROS.2017.8202261
16.
Kuo
,
A. D.
,
2007
, “
,”
IEEE Robot. Autom. Mag.
,
14
(
2
), pp.
18
29
.10.1109/MRA.2007.380653
17.
Atay
,
S.
,
Jenkins
,
T.
,
Buckner
,
G.
, and
Bryant
,
M.
,
2020
, “
Energetic Analysis and Optimization of a Bi-Modal Rolling-Flying Vehicle
,”
Int. J. Intell. Robot. Appl.
,
4
, pp.
3
20
.10.1007/s41315-020-00119-2
18.
Stebler
,
S.
,
Mackunis
,
W.
, and
Reyhanoglu
,
M.
,
2016
, “
Nonlinear Output Feedback Tracking Control of a Quadrotor UAV in the Presence of Uncertainty
,” 14th International Conference on Control, Automation, Robotics and Vision (
ICARCV 2016
), Phuket, Thailand, Nov. 13–15, pp.
1
6
.10.1109/ICARCV.2016.7838569
19.
Mahony
,
R.
,
Kumar
,
V.
, and
Corke
,
P.
,
2012
, “
Multirotor Aerial Vehicles: Modeling, Estimation, and Control of Quadrotor
,”
IEEE Rob. Autom. Mag
,
19
(
3
), pp.
20
32
.10.1109/MRA.2012.2206474
20.
L'Afflitto
,
A.
, and
,
K.
,
2017
, “
Equations of Motion of Rotary-Wing Unmanned Aerial System With Time-Varying Inertial Properties
,”
J. Guid. Control. Dyn.
,
41
(
2
).10.2514/1.G003015
21.
Lee
,
T.
,
Leok
,
M.
, and
Mcclamroch
,
N. H.
,
2012
, “
Nonlinear Robust Tracking Control of a Quadrotor UAV on SE(3)
,” American Control Conference (
ACC
), Montreal, QC, Canada, June 27–29, pp. 4649–4654.10.1109/ACC.2012.6315143
22.
Murray
,
R. M.
,
Sluis
,
W.
,
Rathinam
,
M.
, and
Sluis
,
W.
,
1995
, “
Differential Flatness of Mechanical Control Systems: A Catalog of Prototype Systems
,” Proceedings of the 1995 ASME International Congress and Exposition.
23.
Mizutani
,
S.
,
,
Y.
,
Salaan
,
C. J.
,
Ishii
,
T.
,
Ohno
,
K.
, and
,
S.
,
2015
, “
Proposal and Experimental Validation of a Design Strategy for a UAV With a Passive Rotating Spherical Shell
,” IEEE International Conference on Intelligent Robots and Systems (
IROS
), Hamburg, Germany, Sept. 28–Oct. 2, pp.
1271
1278
.10.1109/IROS.2015.7353532
24.
Whitman
,
A.
,
Clayton
,
G.
, and
Ashrafiuon
,
H.
,
2018
, “
Prediction of Wheel Slipping Limits for Mobile Robots
,”
ASME J. Dyn. Syst. Meas. Control
,
141
(
4
), p.
041002
.10.1115/DSCC2018-9053
25.
Yin
,
H.
,
Chen
,
Y.-H.
, and
Yu
,
D.
,
2019
, “
Controlling an Underactuated Two-Wheeled Mobile Robot: A Constraint-Following Approach
,”
ASME J. Dyn. Syst. Meas. Control
,
141
(
7
), p.
071002
.10.1115/1.4043112
26.
Goldstein
,
H.
,
1950
,
Classical Mechanics
,
,
.
27.
Luca
,
A. D.
, and
Oriolo
,
G.
,
1995
, “
Modelling and Control of Nonholonomic Mechanical Systems
,” Kinematics and Dynamics of Multi-Body Systems, CISM International Centre for Mechanical Sciences (Courses and Lectures, Vol. 360), Springer, Vienna, Austria.
28.
Powers
,
C.
,
Mellinger
,
D.
, and
Kumar
,
V.
,
2015
, “
,”
Handbook of Unmanned Aerial Vehicles
, k. Valavanis and G. Vachtsevanos, eds., Springer, Dordrecht, The Netherlands.
29.
Lu
,
H.
,
Liu
,
C.
,
Coombes
,
M.
,
Guo
,
L.
, and
Chen
,
W.-H.
,
2016
, “
Online Optimisation-Based Backstepping Control Design With Application to Quadrotor
,”
IET Control Theory Appl.
,
10
(
14
), pp.
1601
1611
.10.1049/iet-cta.2015.0976
30.
He
,
X.
,
Kou
,
G.
,
Calaf
,
M.
, and
Leang
,
K. K.
,
2019
, “
In-Ground-Effect Modeling and Nonlinear-Disturbance Observer for Multirotor Unmanned Aerial Vehicle Control
,”
ASME J. Dyn. Syst. Meas. Control
,
141
(
7
), p.
071013
.10.1115/1.4043221