## 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 [1–5]. 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 [6–10], walking-flying [11,12], and swimming-flying [13–15], 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].

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 [18–21], 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\u21c0abC$ 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\u21c0abC$ in the $p$-direction, where $p$ can be $x,\u2009y$, or $z$. The velocity vector $d/dt(r\u21c0abC)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\u21c0abC)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 $\omega \u21c0ABC$ is the angular velocity of frame $B$ with respect to frame $A$ resolved in frame $C$. [ $a\u21c0$]_{x} represents the skew-symmetric matrix formed from the elements of the vector $a\u21c0$ such that [ $a\u21c0$]_{x}$b\u21c0=a\u21c0\xd7b\u21c0$. $CAB$ is a rotation matrix that transforms a vector's frame of resolution from $A$ to $B$ such that $r\u21c0B=CABr\u21c0A$. $i,\u0302j\u0302$, and $k\u0302$ represent the orthonormal basis vectors of $R3$ so that $i\u0302=[1\u2009\u20090\u2009\u20090]T,\u2009j\u0302=[0\u2009\u20091\u2009\u20090]T$, and $k\u0302=[0\u2009\u20090\u2009\u20091]T$. All other variables are defined in the Nomenclature section.

### 2.2 Model Formulation.

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$.

*z*-direction, and $MxB,\u2009MyB$, 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 $\sigma \u0307$ and $\alpha \u0307$ by the kinematic relationship

where $C\psi \omega \u21c0IBB$ is defined by Eq. (6). Note that $C\psi \omega \u21c0IBB\u2009$ is not an orthonormal rotation matrix, but only relates the angular velocity of $B$ to the Euler angle rates.

### 2.3 Analytical Mechanics.

*p*th wheel can be obtained as

*p*th 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,\u2009Iyw,\u2009Iyw}$. The kinetic energy of each wheel is calculated in the same way as Eq. (14); therefore, the total kinetic energy of the RFV is

*z*-direction, then

where $n=1,2,\u2026,6$.

### 2.4 Constraints.

*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

*p*th wheel and the ground is zero. That is

*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*-direction of the $F$ frame. That is

Therefore, the lateral no slip constraint from Eq. (26) is equivalent to saying that $d/dt(r\u21c0ib,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\u21c0$. Therefore, the nonholonomic constraints do not permit reducing the dimension of $q\u21c0$ 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\u0308=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.

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,\u2009x\u0307,\u2009y$, and $y\u0307$ are all resolved in $I$, and $x\u0307$ and $y\u0307$ 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\u0308$ and $y\u0308$ will be with respect to $I$ as well.

### 2.5 Impressed Forces.

*p*th wheel. If the rolling resistance is proportional to the normal force on the wheel then

*p*th wheel, and $\mu 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)

*p*th wheel with respect to the body $B$. This moment arises due to bearing losses and is

### 2.6 Equations of Motion.

*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)

*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\u22123=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

## 3 Differential Flatness

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 $\alpha $ (except when $\alpha =0,\u2009\pi $) 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 $\sigma ,\u2009\alpha ,\u2009\theta 1,\u2009\theta 2$, and $\omega \u21c0IBB$.

*y*-axis. Measurements of $\sigma ,\u2009\gamma ,\alpha $, and $\omega \u21c0IBB$ are obtained from the INS. $\sigma \u0307$ and $\alpha \u0307$ are computed from $\omega \u21c0IBB$ according to

*p*th wheel

*y*-component of Eq. (59) for each wheel yields

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

Comparing $\sigma \u0307odom$ and $\sigma \u0307$ 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,\u2009MyF$, 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 $\alpha $ and $\sigma $ signals. The model parameter values appear in Table 1 under the set 1 column.

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 $\alpha ,\u2009\alpha \u0307$, and $\sigma \u0307$ and so excite the nonlinear dynamics of Eqs. (54)–(56). The motion data consist of $\sigma ,\u2009\alpha ,\u2009\theta 1,\u2009\theta 2$, and $\omega \u21c0IBB$ sampled at 2 ms intervals. Values of $\sigma \u0307$ and $\alpha \u0307$ are calculated from the motion data according to $\sigma \u0307=\omega IB,yBsin\alpha +\omega IB,zBcos\alpha $ and $\alpha \u0307=\omega IB,xB$ (see Eq. (6)). $\theta \u03071$ and $\theta \u03072$ are obtained by numerically differentiating measurements of $\theta 1$ and $\theta 2$. Then $y\u0307F$ 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 $\alpha $ and $\alpha \u0307$, respectively, as functions of time, and indicate good agreement between the simulation and empirical data despite uncertainty in the parameter values. $\alpha $ oscillates and decays, as expected, however the response is not a pure sinusoid, particularly at large values of $\alpha $. This is due primarily to the $mBghsin\alpha $ term in the numerator of Eq. (56), and to a lesser extent the $\u22121/2(Iz\u2212Iy)\sigma \u03072sin2\alpha $ term; the remaining terms in Eq. (56) are more than an order of magnitude smaller than the $mBghsin\alpha $ term. That $\alpha $ is not a pure sinusoid is best observed by examining $\alpha \u0307$ in Fig. 4(b), which exhibits sharp rather than smooth peaks when $\alpha $ is large.

Figures 5(a) and 5(b) show $\sigma $ and $\sigma \u0307$, respectively, as functions of time. The empirical $\sigma \u0307$ data are zero-phase low-pass filtered to remove noise. The rate at which $\sigma \u0307$ decays is influenced primarily by $\mu visc$ and $\mu rr$, with the latter causing the plateau in $\sigma $. The ripple in the empirical $\sigma \u0307$ data is caused by the oscillation in $\alpha $, and is predicted by the $\alpha \u0307\sigma \u0307(Iz\u2212Iy)sin2\alpha $ and $mBh\sigma \u0307y\u0307Fsin\alpha $ terms in Eq. (55); the same ripple can be seen in the simulated data.

The empirical and simulated values of the wheel angles $\theta 1$ and $\theta 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 $\theta 1$ and $\theta 2$ coincide with the oscillation in $\alpha $. Figure 6(b) illustrates $y\u0307F$, which is initially responsible for some of the ripple in $\sigma \u0307$ via the $mBh\sigma \u0307y\u0307Fsin\alpha $ term in Eq. (55), though this term decays quickly and is soon overtaken by the $\alpha \u0307\sigma \u0307(Iz\u2212Iy)sin2\alpha $ term. Both the simulated and empirical values of $y\u0307F$ are calculated according to Eq. (61), which requires knowledge of $\theta \u03071$ and $\theta \u03072$, neither of which are acquired directly from the motion data. Therefore, the empirical value of $y\u0307F$ is obtained by using numerically differentiated values of the empirical $\theta 1$ and $\theta 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.

### 5.4 Demonstrating Gyroscopic Nonlinearity.

The primary cross-coupling between the $\alpha $ and $\sigma $ degrees-of-freedom is due to the “gyroscopic” term containing $Iz\u2212Iy$ 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\u2212Iy\u22600$. 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=\u22120.002$ m. This results in the $Iz\u2212Iy$ 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 $\alpha $ and $\alpha \u0307$ appear in Figs. 7(a) and 7(b), respectively. Initially, $\alpha $ is near $\pi $ radians, indicating that the RFV chassis is upside down. Since $h<0$, this configuration is statically unstable due to the gravitational term $mBghsin\alpha $ in Eq. (56). However, rather than immediately diverging toward the stable equilibrium at $\alpha =0$, $\alpha $ instead oscillates about $\pi $ radians for several seconds (see orange oval in Fig. 7(a)) due to the gyroscopic term $\u22121/2(Iz\u2212Iy)\sigma \u03072sin2\alpha $. This “gyroscopically” stabilized behavior continues until of the amplitude of the gyroscopic term decreases (as $\sigma \u0307$ decreases) and the gravitational term $mBghsin\alpha $ begins to dominate. That is, $\alpha =\pi $ behaves like a stable equilibrium of Eq. (56) when $\sigma \u0307$ is viewed as a parameter, with the equilibrium undergoing a bifurcation as $\sigma \u0307$ drops below a critical value. The simulation predicts the same behavior, with $\alpha $ briefly oscillating about $\alpha =\pi $ and then diverging toward $\alpha =2\pi \u2009(=0)$.

The numerical simulation results demonstrate the same qualitative behavior in $\sigma ,\u2009\alpha ,\u2009\theta 1,\u2009\theta 2$, and $y\u0307F$ 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 $\sigma $, $\alpha $, and $y\u0307F$ (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., $\u22121/2(Iz\u2212Iy)\sigma \u03072sin2\alpha $ 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 $\alpha $ and $\sigma $ degrees-of-freedom is the “gyroscopic” term $Iz\u2212Iy$ appearing in both Eqs. (55) and (56). This suggests that the coupling could be nearly eliminated by designing the RFV such that $Iz\u2212Iy\u22480$. Decoupling the degrees-of-freedom and removing nonlinearities would simplify the design of closed-loop controllers for $\alpha $ and $\sigma $.

### 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\psi \omega \u21c0IBB$ =
Euler angle rates-to-angular velocity transformation matrix

- $d/dt(\u22c5)A$ =
time derivative of $(\u22c5)$ with respect to the $A$ reference frame

- $d2/dt2(\u22c5)AB$ =
time derivative of $d/dt(\u22c5)A$ with respect to the $B$ reference frame

- $Dy,D\alpha ,D\sigma $ =
dissipation terms

- $F\u21c0constraint$ =
generalized constraint force

- $Fdrag$ =
aerodynamic drag force (N)

- $FNp$ =
*p*th wheel normal force - $FzB$ =
body force (N)

- $g$ =
gravitational acceleration (m/s

^{2}) - $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 m

^{2}) - $Ixw,Iyw$ =
wheel principal moments of inertia of (kg m

^{2}) - $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$ =
*p*th wheel rolling resistance moment (N·m) - $Mvisc,p$ =
*p*th wheel viscous damping moment (N·m) - $MxB,MyB,MzB$ =
body moments (N·m)

- $q\u21c0$ =
generalized coordinates

- $Q\u21c0,Q\u21c0grav,Q\u21c0control,Q\u21c0rr,Q\u21c0visc,Q\u21c0drag$ =
generalized forces

- $R$ =
wheel radius (m)

- $r\u21c0ab$ =
position of point $b$ relative to point $a$

- $T,TB,TW$ =
kinetic energies

- $x$ =
position coordinate (m)

- $y$ =
position coordinate (m)

- $y\u0307F,y\u0308F$ =
velocity and acceleration in the $F$ frame

*y*-direction - $z\u21c0,z\u21c0\u0307,z\u21c0\u0308$ =
flat outputs and derivatives

- $\alpha $ =
pitch Euler angle (radians)

- $\gamma $ =
tilt Euler angle

- $\theta 1,\theta 2$ =
wheels angles (radians)

- $\lambda \u21c0$ =
Lagrange multipliers

- $\mu visc$ =
viscous damping coefficient (N·m s)

- $\mu rr$ =
rolling resistance coefficient

- $\sigma $ =
yaw Euler angle (radians)

- $\omega \u21c0AB$ =
angular velocity of reference frame $B$ with respect to reference frame $A$

- $\u2202f\u21c0/\u2202q\u21c0$ =
differential constraint matrix

### Superscripts

## References

**172**, p. 03008. 10.1051/matecconf/201817203008

**13**(5).

*Kinematics and Dynamics of Multi-Body Systems, CISM International Centre for Mechanical Sciences*(Courses and Lectures, Vol. 360), Springer, Vienna, Austria.