Steering Control of an Autonomous Unicycle

1 Steering Control of an Autonomous Unic ycle M ´ at ´ e B. V izi, G ´ abor Orosz, D ´ enes T ak ´ acs and G ´ abor St ´ ep ´ an Abstract —The steering control of an autonomous unicycle is considered. The underlying dynamical model of a single r olling wheel is discussed regarding the steady state motions and their stability . The unicycle model is intr oduced as the simplest possible extension of the rolling wheel wher e the location of the center of gravity is contr olled. With the help of the Appellian approach, a state space repr esentation of the controlled nonholonomic system is built in a way that the most compact nonlinear equations of motions are constructed. Based on controllability analysis, feedback controllers are designed which successfully carry out lane changing and turning maneuvers. The behavior of the closed- loop system is demonstrated by numerical simulations. Index T erms —Unicycle, Nonholonomic dynamics, Stability , Feedback control, Maneuvering I . I N T RO D U C T I O N M ICR O-MOBILITY solutions are spreading rapidly in urban en vironments [ 1 ]. Among these, human-ridden electric unicycles (EUCs) become more and more popular transportation de vices; see Figure 1 (a). These micro-mobility vehicles can match the speed of automobiles in urban traf fic while their compact size make them appealing for commute in congested en vironments. Due to the three dimensional spatial rolling of the wheel and the stabilization of an unstable equilibrium, the unique dynamics of the unicycle combines agility and maneuverability . T o exploit these properties, one may consider making EUCs autonomous (see Figure 1 (b)) which opens up a challenging a venue for modeling, dynamics and control. During the last few decades, several autonomous unicycle designs have appeared in the literature which differ on various aspects such as the number and/or types of actuators that can be used for control. The first publication related to autonomous unicycles known to the authors is [ 2 ] in which the longi- tudinal/pitch motion is controlled by balancing an in verted pendulum, and the lateral/tilt motion is controlled by mo ving a mass perpendicular to the wheel. T wo other approaches are presented in [ 3 ]. In the first case, the longitudinal/pitch motion The research reported in this paper has been supported by the Hungarian National Science Foundation under Grant No. NKFI K 132477, NKFI KKP 133846 and by the HUN-REN Hungarian Research Network; and also by Project no. TKP-6-6/P AL Y -2021 provided by the Ministry of Culture and Inno vation of Hungary from the National Research, De velopment and Innov ation Fund, financed under the TKP2021-NV A funding scheme. The research of D. T . was supported by a J ´ anos Bolyai Research Scholarship of the Hungarian Academy of Sciences. M ´ at ´ e B. V izi and D ´ enes T ak ´ acs are with the Department of Applied Mechanics, Budapest Univ ersity of T echnology and Economics, Budapest, Hungary and HUN-REN–BME Dynamics of Machines Research Group, Budapest, Hungary (e-mails: vizi@mm.bme.hu, takacs@mm.bme.hu) G ´ abor Orosz is with the Department of Mechanical Engineering and with the Department of Civil and Environmental Engineering, University of Michigan, Ann Arbor , MI 48109, USA (e-mail: orosz@umich.edu). G ´ abor St ´ ep ´ an is with the Department of Applied Mechanics, Budapest Univ ersity of T echnology and Economics, Budapest, Hungary (e-mail: stepan@mm.bme.hu). Fig. 1. Human riding a unicycle (a), and a simple autonomous EUC (b). is also controlled by balancing an inv erted pendulum, while the turning/yaw motion of the unicycle is controlled by an ov erhead flywheel. In the second case, the tilt is controlled by adding a second pendulum swinging in the lateral plane. The o verhead flywheel approach was further explored in [ 4 ], [ 5 ]; the lateral pendulum approach can be found in [ 6 ]. The tilt motion of the unicycle can also be controlled by a lateral flywheel, see, for example, [ 7 ]–[ 9 ], while the combination of overhead and lateral flywheels can be found in [ 10 ], [ 11 ]. Furthermore, the application of gyroscopes for lateral stabilization and steering is presented in [ 12 ]–[ 15 ]. Humanoid- type autonomous unicycles are introduced and analyzed in [ 16 ], [ 17 ], see Figure 1 (a) as illustration. Most of the above approaches provide satisfactory dynamic behavior , but their complexity prohibits closed-form analysis of the controlled system. Our goal here is to develop a simple autonomous unicycle model that is capable of carrying out a v ariety of maneuvers while it can still be inv estigated analytically . Thus, we consider the simple mechanical model shown in Figure 1 (b), which consists of a rolling wheel and added mass that can be mov ed along the axle to balance the lateral motion. W e derive the equations of motion for this simple system, analyze the stability properties, and design controllers for maneuvering the autonomous unic ycle. The rolling of the wheel can be described using kinematic constraints [ 18 ]–[ 21 ]. Thus, the unicycle is considered as a nonholonomic mechanical system. Such systems are often described by the generalized Lagrangian equations of the second kind (or Routh–V oss equations) [ 22 ], [ 23 ]. This method yields a differential-algebraic system of equations. Howe ver , eliminating the algebraic variables to obtain a system of ordinary differential equations that is appropriate for control design is a challenging task. Alternativ ely , for conservati ve 2 nonholonomic systems, a Routhian-like model reduction tech- nique, the so-called Lagrange–d’Alembert–Poincar ´ e equations can be used to identify conserv ed quantities and also to e xplore symmetries [ 24 ]. The Appellian approach [ 21 ], [ 25 ], [ 26 ], which is used in this study , results in a system of first order ordinary dif ferential equations as a compact and simple representation of the underlying nonholonomic system. Moreov er , an innov ativ e definition of the pseudo velocities can significantly reduce the algebraic complexity of the resulting equations of motion while describing the same dynamical system, which simplifies the subsequent analysis. This enables one to deploy a plethora of control techniques. The cost of the Appellian formalism is that accelerations must be calculated, howe ver , the mentioned benefits justify this cost in the examples considered here. Further details about nonholonomic systems may be found in [ 27 ]–[ 34 ]. In this study , we first omit the mass moving along the axle in Figure 1 (b) and explore the dynamics of the uncontrolled rolling wheel. W e categorize dif ferent steady state motions (e.g., straight rolling, turning) of interest. A strong sufficient condition of stable rolling was given by Bloch in [ 34 ], while here the local necessary and suf ficient condition is derived. This simplified case also enables us to explain the self stabi- lizing effects in the tilt direction abov e a critical speed. By adding the moving mass and an internal force between the wheel and the mass, we create a control system and study how the steady states are af fected. This enables us to design feedback controllers which stabilize the steady states at any speed and enables the unicycle to perform maneuvers such as lane changes and sharp turns. Our control design exploits the inherent instabilities of the system in order to demonstrate high lev el of maneuv erability . The article is structured as follo ws. In Section II , the modeling frame work and steady state analysis of the rolling wheel example are presented. Section III introduces a nov el autonomous unicycle model and analyzes the steady states of the open-loop system. Section IV proposes controller designs that successfully perform lane changing and turning maneu- vers. The performance of these controllers are demonstrated by numerical simulations. W e conclude our results in Section V and provide future research directions. I I . D Y NA M I C S O F T H E R O L L I N G W H E E L W e introduce the modeling framework and notation on the rolling wheel example which represents the uncontrolled be- havior of the unic ycle. W e re veal the dynamical characteristics, including the self-stabilization phenomenon, which can be exploited for control design. A. Governing equations A wheel has 𝑁 = 6 de grees of freedom (DoF) in the three dimensional space. That is, the spatial position and orientation are described by six v ariables: the position of the center of gravity r G = [ 𝑥 G 𝑦 G 𝑧 G ] T and the yaw ( 𝜓 ), tilt ( 𝜗 ) and pitch ( 𝜑 ) angles; see Figure 2 . Note that tilt is often called roll in the vehicle dynamics literature but here we do not use this Fig. 2. Mechanical model of the rolling wheel con vention to av oid confusion with the fact that disc rolls on the horizontal plane. T o describe the motion of the wheel, three coordinate frames are introduced; see Figure 2 . The axes 𝑥 0 and 𝑦 0 of the ground fixed frame F 0 span the horizontal plane and its 𝑧 0 axis gi ves the v ertical direction. The frame F 1 is mo ving with the wheel such that its origin is the wheel-ground contact point P . This frame is rotated with respect to F 0 around the 𝑧 0 axis with yaw angle 𝜓 so that the 𝑥 1 axis is tangential to the path of P (the fixed polode). The frame F 2 is rotated with respect to F 1 around the 𝑥 1 axis with the tilt angle 𝜗 so the 𝑥 1 and 𝑧 1 axes span the plane of the wheel and the 𝑦 2 axis is aligned with the wheel axle. The origin of frame F 2 is placed at the wheel center point G . Assume that the wheel rolls without slipping; the kinematic condition of rolling is that the instantaneous center of rotation coincides with the contact point P : v P = 0 . (1) This yields the kinematic constraints: ¤ 𝑥 G = ¤ 𝜓 𝑅 cos 𝜓 sin 𝜗 + ¤ 𝜗 𝑅 sin 𝜓 cos 𝜗 + ¤ 𝜑 𝑅 cos 𝜓 , ¤ 𝑦 G = ¤ 𝜓 𝑅 sin 𝜓 sin 𝜗 − ¤ 𝜗 𝑅 cos 𝜓 cos 𝜗 + ¤ 𝜑 𝑅 sin 𝜓 , (2) and the geometric constraint ¤ 𝑧 G = − 𝑅 ¤ 𝜗 sin 𝜗 ⇒ 𝑧 G = 𝑅 cos 𝜗 , (3) where the vertical position 𝑧 G depends on the tilt angle 𝜗 ; see Figure 2 . Therefore the rolling wheel is a nonholonomic mechanical system with 𝑛 g = 1 geometric constraint ( 3 ) and 𝑛 k = 2 kinematic constraints ( 2 ). Note that the geometric and kinematic constraints are also referred to as holonomic and nonholonomic constraints, respecti vely . The equations of mo- tion are deri ved in Appendix A using the Appellian approach [ 21 ], [ 25 ], [ 26 ] to provide the most compact algebraic form. According to the number of geometric and kinematic con- straints, 𝑛 𝑞 = 6 − 𝑛 g = 5 generalized coordinates hav e to be chosen to describe the system unambiguously; let these be:  𝑥 G , 𝑦 G , 𝜓 , 𝜗 , 𝜑  . (4) 3 Moreov er , 𝑛 𝜎 = 𝑛 𝑞 − 𝑛 k = 3 pseudovelocities ha ve to be cho- sen; let these be defined by the components of the angular velocity 𝝎 resolved in frame F 2 (cf. ( 73 ) in Appendix A ): 𝜔 1 : = ¤ 𝜗, 𝜔 2 : = ¤ 𝜓 sin 𝜗 + ¤ 𝜑 , 𝜔 3 : = ¤ 𝜓 cos 𝜗 . (5) Then, the Appellian approach yields the equations of motion:                  ¤ 𝜔 1 = 6 5 𝜔 2 𝜔 3 − 1 5 𝜔 2 3 tan 𝜗 + 4 𝑔 5 𝑅 sin 𝜗 , ¤ 𝜔 2 = − 2 3 𝜔 1 𝜔 3 , ¤ 𝜔 3 = − 2 𝜔 1 𝜔 2 + 𝜔 1 𝜔 3 tan 𝜗 , ¤ 𝜗 = 𝜔 1 ,              ¤ 𝜓 = 𝜔 3 1 cos 𝜗 , ¤ 𝜑 = 𝜔 2 − 𝜔 3 tan 𝜗 , ¤ 𝑥 G = 𝜔 1 𝑅 sin 𝜓 cos 𝜗 + 𝜔 2 𝑅 cos 𝜓 , ¤ 𝑦 G = − 𝜔 1 𝑅 cos 𝜓 cos 𝜗 + 𝜔 2 𝑅 sin 𝜓 , (6) that is, the rolling wheel is an 𝑛 = 6 − 𝑛 g − 𝑛 k / 2 = 4 DoF nonholonomic mechanical system. The equations in ( 6 ) are ordered such that the system can be separated into essential dynamics (the first four equations) and hidden dynamics (the second four equations) where the essential dynamics is independent of the hidden dynamics [ 22 ]. The equations are in the form ¤ x = 𝑓 ( x ) where the state is defined as x =  𝜔 1 𝜔 2 𝜔 3 𝜗 𝜓 𝜑 𝑥 G 𝑦 G  T , (7) where the dashed line separates the essential states from the cyclic coordinates describing hidden motion. B. Steady state motions The rolling wheel exhibits a steady state motion when the essential dynamics (first four equations in ( 6 )) possess an equilibrium. That is, the pseudovelocities and the tilt angle are constants: 𝜔 1 ( 𝑡 ) ≡ 𝜔 1 ∗ , 𝜔 2 ( 𝑡 ) ≡ 𝜔 2 ∗ , 𝜔 3 ( 𝑡 ) ≡ 𝜔 3 ∗ , 𝜗 ( 𝑡 ) ≡ 𝜗 ∗ . (8) Then, according to the fourth equation in ( 6 ), the tilt rate must be zero, i.e., 𝜔 1 ∗ = 0 . Substituting this into ( 6 ) the first equation yields 6 5 𝜔 2 ∗ 𝜔 3 ∗ − 1 5 𝜔 2 3 ∗ tan 𝜗 ∗ + 4 𝑔 5 𝑅 sin 𝜗 ∗ = 0 , (9) while the hidden motion can be expressed as ¤ 𝜓 ( 𝑡 ) ≡ ¤ 𝜓 ∗ = 𝜔 3 ∗ 1 cos 𝜗 ∗ , ¤ 𝑥 G ( 𝑡 ) = 𝜔 2 ∗ 𝑅 cos 𝜓 ( 𝑡 ) , ¤ 𝜑 ( 𝑡 ) ≡ ¤ 𝜑 ∗ = 𝜔 2 ∗ − 𝜔 3 ∗ tan 𝜗 ∗ , ¤ 𝑦 G ( 𝑡 ) = 𝜔 2 ∗ 𝑅 sin 𝜓 ( 𝑡 ) . (10) That is, the yaw rate ¤ 𝜓 and the pitch rate ¤ 𝜑 are constants, while the horizontal velocity components ¤ 𝑥 G , ¤ 𝑦 G of the center of gravity v ary with time through the yaw angle 𝜓 . Integrating ( 10 ), the generalized coordinates become 𝜓 ∗ ( 𝑡 ) = ¤ 𝜓 ∗ 𝑡 + 𝜓 0 , 𝜑 ∗ ( 𝑡 ) = ¤ 𝜑 ∗ 𝑡 + 𝜑 0 , 𝑥 G ∗ ( 𝑡 ) =         ¤ 𝜑 ∗ ¤ 𝜓 ∗ + sin 𝜗 ∗  𝑅 sin ( ¤ 𝜓 ∗ 𝑡 + 𝜓 0 ) + 𝑥 0 if ¤ 𝜓 ∗ ≠ 0 , ¤ 𝜑 ∗ 𝑡 𝑅 cos 𝜓 0 + 𝑥 0 if ¤ 𝜓 ∗ = 0 , 𝑦 G ∗ ( 𝑡 ) =        −  ¤ 𝜑 ∗ ¤ 𝜓 ∗ + sin 𝜗 ∗  𝑅 cos ( ¤ 𝜓 ∗ 𝑡 + 𝜓 0 ) + 𝑦 0 if ¤ 𝜓 ∗ ≠ 0 , ¤ 𝜑 ∗ 𝑡 𝑅 sin 𝜓 0 + 𝑦 0 if ¤ 𝜓 ∗ = 0 , (11) where the 𝜓 0 and 𝜑 0 denote the initial yaw and pitch angles while 𝑥 0 and 𝑦 0 originate in the initial position of point G . The center of gravity G follows a circular path of radius 𝜌 G = | ¤ 𝜑 ∗ / ¤ 𝜓 ∗ + sin 𝜗 ∗ | 𝑅 if the yaw rate is not zero ( ¤ 𝜓 ∗ ≠ 0 ). Correspondingly , the contact point P draws a circle of radius 𝜌 P =   ¤ 𝜑 ∗ / ¤ 𝜓 ∗   𝑅 on the ground plane; see Figure 2 . W e refer to this motion as turning-rolling in the rest of the paper . For zero yaw rate ( ¤ 𝜓 ∗ = 0 ), the center or gravity moves along a straight path. W e refer to this as straight rolling in the rest of the paper . Keep in mind that the steady state tilt angle, yaw rate and pitch rate are not independent of each other . T o establish a relationship between these quantities, one must replace the steady state pseudovelocities ( 𝜔 1 ∗ , 𝜔 2 ∗ , 𝜔 3 ∗ ) in ( 9 ) with the generalized velocities ( ¤ 𝜓 ∗ , ¤ 𝜑 ∗ ) according to ( 5 ). This yields ¤ 𝜓 2 ∗ sin 𝜗 ∗ cos 𝜗 ∗ + 6 5 ¤ 𝜓 ∗ ¤ 𝜑 ∗ cos 𝜗 ∗ + 4 𝑔 5 𝑅 sin 𝜗 ∗ = 0 , (12) which is visualized in Figure 3 (a). Since ( 12 ) is linear in the pitch rate ¤ 𝜑 ∗ , one may express ¤ 𝜑 ∗ as a function of the the tilt angle 𝜗 ∗ and the yaw rate ¤ 𝜓 ∗ for turning-rolling ( ¤ 𝜓 ∗ ≠ 0 ). For straight rolling ( ¤ 𝜓 ∗ = 0 ), ( 12 ) reduces to sin 𝜗 ∗ = 0 , (13) yielding 𝜗 ∗ = 0 , that is, the wheel must be non-tilted. Observ e that ( 13 ) is independent of the pitch rate ¤ 𝜑 ∗ , so the straight rolling is feasible for arbitrary pitch rates. Another special case is when the pitch rate is zero ( ¤ 𝜑 ∗ = 0 ). Then ( 12 ) results in ¤ 𝜓 2 ∗ sin 𝜗 ∗ cos 𝜗 ∗ + 4 𝑔 5 𝑅 sin 𝜗 ∗ = 0 , (14) which can only hold when the tilt angle is zero, that is, 𝜗 ∗ = 0 . W e refer to this solution as spinning on the spot in the rest of the paper . Finally , the very special case ¤ 𝜓 ∗ = ¤ 𝜑 ∗ = 0 corresponds to the static equilibrium of the standing disc, which is, indeed, unstable. The stability of the steady states (turning-rolling, straight rolling and spinning) are analyzed consecutiv ely below . C. Stability of the steady state motions The steady state motions of the rolling wheel are defined with respect to the essential motion independently of the cyclic states of the system. Accordingly , when we talk about the stability of the steady state motions, it is only about the essential dynamics, while the hidden motions may present instability in the L yapunov sense. 4 Let us consider two examples when the wheel initially rolls straight along the 𝑥 axis. In the first case, a small perturbation may increase the speed of rolling. If the essential states 𝜔 1 , 𝜔 2 , 𝜔 3 , 𝜗 remain in the small vicinity of the original v alues then the essential motion is stable. Howe ver , the increased velocity 𝜔 2 causes the cyclic coordinates 𝜑, 𝑥 G mov e away linearly in time from the original unperturbed trajectories; therefore, the hidden motion is unstable. In the second case, let the perturbation cause a small tilt angle 𝜗 0 , tilt rate ¤ 𝜗 0 = 𝜔 10 and/or a small yaw rate ¤ 𝜓 0 = 𝜔 30 . If the essential states remain in the vicinity of the original unperturbed values then the essential motion is stable. Ho wev er , this small perturbation yields a circular motion of the wheel with lar ge radius instead of the straight motion, so the hidden motion is unstable. For simplicity assume the initial values 𝜓 0 = 0 , 𝜑 0 = 0 , 𝑥 0 = 0 , and 𝑦 0 = 0 . According to ( 5 ), ( 8 ), and ( 11 ), the steady state motion is gi ven as x ∗ =                 𝜔 1 ∗ 𝜔 2 ∗ 𝜔 3 ∗ 𝜗 ∗ 𝜓 ∗ ( 𝑡 ) 𝜑 ∗ ( 𝑡 ) 𝑥 G ∗ ( 𝑡 ) 𝑦 G ∗ ( 𝑡 )                 =                 0 ¤ 𝜓 ∗ sin 𝜗 ∗ + ¤ 𝜑 ∗ ¤ 𝜓 ∗ cos 𝜗 ∗ 𝜗 ∗ ¤ 𝜓 ∗ 𝑡 ¤ 𝜑 ∗ 𝑡 . . . . . .                 , (15) where the dots in 𝑥 G ∗ ( 𝑡 ) and 𝑦 G ∗ ( 𝑡 ) refer to ( 11 ). Introducing the state perturbations ˜ x : = x − x ∗ , ( 6 ) leads to the linearized equations of the form ¤ ˜ x = A ˜ x where the state matrix is obtained as A = 𝑓 ′ ( x ∗ ) . This yields A =                 0 𝐴 12 𝐴 13 𝐴 14 0 0 0 0 𝐴 21 0 0 0 0 0 0 0 𝐴 31 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 𝐴 53 𝐴 54 0 0 0 0 0 1 𝐴 63 𝐴 64 0 0 0 0 𝐴 71 𝐴 72 0 0 𝐴 75 0 0 0 𝐴 81 𝐴 82 0 0 𝐴 85 0 0 0                 , (16) where 𝐴 12 = 6 5 ¤ 𝜓 ∗ cos 𝜗 ∗ , 𝐴 13 = 4 5 ¤ 𝜓 ∗ sin 𝜗 ∗ + 6 5 ¤ 𝜑 ∗ , 𝐴 14 = − 1 5 ¤ 𝜓 2 ∗ + 4 𝑔 5 𝑅 cos 𝜗 ∗ , 𝐴 21 = − 2 3 ¤ 𝜓 ∗ cos 𝜗 ∗ 𝐴 31 = − ¤ 𝜓 ∗ sin 𝜗 ∗ − 2 ¤ 𝜑 ∗ , 𝐴 53 = 1 cos 𝜗 ∗ , 𝐴 54 = ¤ 𝜓 ∗ tan 𝜗 ∗ , 𝐴 63 = − tan 𝜗 ∗ , 𝐴 64 = − ¤ 𝜓 ∗ cos 𝜗 ∗ , 𝐴 71 = 𝑅 sin 𝜓 ∗ cos 𝜗 ∗ , 𝐴 72 = 𝑅 cos 𝜓 ∗ , 𝐴 75 = − 𝑅 sin 𝜓 ∗  ¤ 𝜓 ∗ sin 𝜗 ∗ + ¤ 𝜑 ∗  , 𝐴 81 = − 𝑅 cos 𝜓 ∗ cos 𝜗 ∗ , 𝐴 82 = 𝑅 sin 𝜓 ∗ , 𝐴 85 = 𝑅 cos 𝜓 ∗  ¤ 𝜓 ∗ sin 𝜗 ∗ + ¤ 𝜑 ∗  , (17) and the pitch rate ¤ 𝜑 ∗ must be substituted according to ( 12 ). The dashed lines distinguish the parts of A related to the essential and to the hidden dynamics. The linear stability of the steady state (essential) motion can be determined by the characteristic equation det ( 𝜆 I − ˆ A ) = 0 , where ˆ A is the top quadrant of A in ( 16 ). This yields  𝜆 2 − 𝐴 12 𝐴 21 − 𝐴 13 𝐴 31 − 𝐴 14  𝜆 2 = 0 , (18) and the characteristic roots become 𝜆 1 , 2 = ±  4 𝑔 5 𝑅 cos 𝜗 ∗ −  ¤ 𝜓 2 ∗ + 14 5 ¤ 𝜓 ∗ ¤ 𝜑 ∗ sin 𝜗 ∗ + 12 5 ¤ 𝜑 2 ∗  , 𝜆 3 , 4 = 0 . (19) The corresponding steady state is unstable if the characteristic roots are real, 𝜆 1 , 2 ∈ R and 𝜆 1 = − 𝜆 2 . Otherwise, if the first two characteristic roots constitute a purely imaginary complex conjugate pair 𝜆 1 , 2 ∈ C and 𝜆 1 = ¯ 𝜆 2 , a ‘doubtful L yapunov’ case arises due to the zero eigen value of multiplicity two. One may show , howe ver , that both the algebraic and geometric multiplicity of the zero eigenv alue are two, and in Appendix B the L yapunov stability of the essential dynamics is proven. Note that the full motion is not stable in the L yapunov sense due to the zero eigen value of algebraic multiplicity four and geometric multiplicity one belonging to the hidden dynamics (see the lower right quadrant of state matrix ( 16 )). T o sum up the analysis of the linear model, the necessary and sufficient condition for the stable rolling of the wheel is ¤ 𝜓 2 ∗ + 14 5 ¤ 𝜓 ∗ ¤ 𝜑 ∗ sin 𝜗 ∗ + 12 5 ¤ 𝜑 2 ∗ − 4 𝑔 5 𝑅 cos 𝜗 ∗ > 0 , (20) where ( 12 ) holds. When varying the tilt angle 𝜗 ∗ , the yaw rate ¤ 𝜓 ∗ and pitch rate ¤ 𝜑 ∗ , one may change the stability of the steady-state motion. 1) T urning-r olling: In case of a general turning-rolling type steady state motion, neither the tilt angle nor the yaw and pitch rates are zero, i.e., 𝜗 ∗ ≠ 0 , ¤ 𝜓 ∗ ≠ 0 and ¤ 𝜑 ∗ ≠ 0 . Then stability of the steady state changes at the critical ya w rates ¤ 𝜓 crit , 1 , 2 =  2 𝑔 5 𝑅 v t 3 − 6 cos 2 𝜗 ∗ ±  76 sin 4 𝜗 ∗ − 96 sin 2 𝜗 ∗ + 9 ( 2 sin 2 𝜗 ∗ − 3 ) cos 𝜗 ∗ (21) for | 𝜗 ∗ | ≤ V = arcsin   12 19 − 9 √ 5 38  ≈ 18 . 62 ◦ . (22) For a gi ven tilt angle | 𝜗 ∗ | ≤ V , if | ¤ 𝜓 ∗ | < ¤ 𝜓 crit , 1 or ¤ 𝜓 crit , 2 < | ¤ 𝜓 ∗ | then the characteristic roots in ( 19 ) are purely imaginary 𝜆 1 , 2 ∈ C and 𝜆 1 = ¯ 𝜆 2 and the corresponding motion is sta- ble. Otherwise, if ¤ 𝜓 crit , 1 < | ¤ 𝜓 ∗ | < ¤ 𝜓 crit , 2 then the characteristic roots are real 𝜆 1 , 2 ∈ R and 𝜆 1 = − 𝜆 2 , and the corresponding motion is unstable. Notice that, if the tilt angle is large enough | 𝜗 ∗ | > V then all corresponding steady state motion are stable, since V is independent of any physical parameter . This is a fundamental physical constant related to the steady state stability of any rolling wheel. 2) Straight r olling: As discussed abov e the straight rolling steady state motion ( ¤ 𝜓 ∗ = 0 ) can only occur for zero tilt angle 𝜗 ∗ = 0 . That is, the steady state ( 15 ) simplifies to 𝜔 2 ∗ = ¤ 𝜑 ∗ , 𝜑 ∗ ( 𝑡 ) = ¤ 𝜑 ∗ 𝑡 , 𝑥 G ∗ ( 𝑡 ) = 𝑅 ¤ 𝜑 ∗ 𝑡 , (23) 5 Fig. 3. Lev el sets and stability of steady state motion ( 12 ) of the rolling wheel: stable and unstable steady states are denoted by green and red, respectively , the blue line is the stability boundary ( 20 ) and the dotted blue line is the folding of the surface. In panel (d), white denotes that no steady state exists, while yellow corresponds to the cases when there are both stable and unstable steady states due to the folding of the surface. and all the other states become zero. The state matrix ( 16 ) is A =                 0 0 6 5 ¤ 𝜑 ∗ 4 𝑔 5 𝑅 0 0 0 0 0 0 0 0 0 0 0 0 − 2 ¤ 𝜑 ∗ 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 𝑅 0 0 0 0 0 0 − 𝑅 0 0 0 𝑅 ¤ 𝜑 ∗ 0 0 0                 , (24) and the characteristic roots ( 19 ) simplify to 𝜆 1 , 2 = ±  4 𝑔 5 𝑅 − 12 5 ¤ 𝜑 2 ∗ . (25) This leads to the critical pitch rate ¤ 𝜑 crit =  𝑔 3 𝑅 , (26) that is, the straight rolling steady state is stable if the wheel rolls fast enough, | ¤ 𝜑 ∗ | > ¤ 𝜑 crit . Otherwise, it is unstable. W e remark that instead of the pitch rate ¤ 𝜑 ∗ one may introduce the parameter 𝑣 ∗ = ¤ 𝜑 ∗ 𝑅 which is the steady state velocity of the center of gravity . This way the critical pitch rate ( 26 ) can be con verted to a critical velocity . 3) Spinning on the spot: The spinning steady state motion ( ¤ 𝜑 ∗ = 0 ) can also occur for zero tilt angle 𝜗 ∗ = 0 . Then the steady state becomes 𝜔 3 ∗ = ¤ 𝜓 ∗ , 𝜓 ∗ = ¤ 𝜓 ∗ 𝑡 , (27) and all the other states are zero. The state matrix ( 16 ) becomes A =                 0 6 5 ¤ 𝜓 ∗ 0 4 𝑔 5 𝑅 − 1 5 ¤ 𝜓 2 ∗ 0 0 0 0 − 2 3 ¤ 𝜓 ∗ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 − ¤ 𝜓 ∗ 0 0 0 0 𝑅 sin 𝜓 ∗ 𝑅 cos 𝜓 ∗ 0 0 0 0 0 0 − 𝑅 cos 𝜓 ∗ 𝑅 sin 𝜓 ∗ 0 0 0 0 0 0                 , (28) while the characteristic roots ( 19 ) and the critical yaw rate ( 21 ) simplify to 𝜆 1 , 2 = ±  4 𝑔 5 𝑅 − ¤ 𝜓 2 ∗ , ¤ 𝜓 crit =  4 𝑔 5 𝑅 . (29) That is, the spinning steady state is stable if the wheel spins fast enough | ¤ 𝜓 ∗ | > ¤ 𝜓 crit . Otherwise, it is unstable. All results about the steady state motions and their stability are sho wn in Figure 3 . P anel (a) sho ws the steady state motions 6 T ABLE I N U ME R I C AL P A R A ME T E RS F O R T H E RO L L IN G WH E E L A N D U N I CY C L E Quantity Symbol V alue Unit wheel mass 𝑚 10 kg point mass 𝑚 0 5 kg wheel radius 𝑅 0.3 m gravitational acceleration 𝑔 9.81 m/s 2 ( 12 ) as surface in the ( 𝜗 ∗ , ¤ 𝜓 ∗ , ¤ 𝜑 ∗ ) space. Panel (b) sho ws the lev el sets of 𝜗 ∗ in the ( ¤ 𝜓 ∗ , ¤ 𝜑 ∗ ) plane. Similarly , panels (c) and (d) sho w the lev el sets of ¤ 𝜓 ∗ and ¤ 𝜑 ∗ , respectiv ely . Green ( ■ ■ ) and red ( ■ ■ ) areas represent stable and unstable steady state solutions, respectiv ely , while the solid blue line ( — ) is the stability boundary ( 20 ). The dashed blue line ( – – ) in panel (c) sho ws where ¤ 𝜑 ∗ → ±∞ . The dotted blue curves ( · · · ) in panels (b) and (c) indicate the folding of the surface and these bound the hour glass-shaped area in panel (d). In the dark green area ( ■ ■ ) in panel (d), two stable steady state states exist with different yaw rates ¤ 𝜓 ∗ , while in the yellow area ( ■ ■ ) a stable and an unstable steady state exist. In the white area, no steady state exists regardless of the yaw rate due to the folding of the surface. I I I . M O D E L I N G T H E U N I C Y C L E The unicycle is shown in Figure 4 , where the motion is controlled by moving a mass point along the axle of the wheel. In this section, we present the gov erning equations and analyze the steady states of the open loop system. A. Governing equations The model of the unicycle is built upon the rolling wheel model presented in Section II with the added point mass 𝑚 0 moving along the axle. The mass’ relative position is giv en by r GA =  0 𝑟 0  T F 2 , (30) where 𝑟 represents the position of the point mass along the axle while other two coordinates are constrained at zero by Fig. 4. Model of the unicycle for maneuvering the axle, resulting in two geometric constraints. Between the disc and the point mass, there is an internal pair of forces; the control force acting on the point mass is F =  0 𝑢 0  T F 2 , (31) while − F acts on the wheel with opposite sense. The control input 𝑢 is used for balancing and steering the unic ycle. The equations of motion of the unicycle are deriv ed in Appendix C in a similar way as for the rolling wheel; the differences between the two models are highlighted below . The rolling constraints of the wheel still hold, yielding two kinematic constraints ( 2 ) and one geometric constraint ( 3 ). Thus, the unicycle is a nonholonomic mechanical system with 𝑛 g = 3 geometric constraints and 𝑛 k = 2 kinematic constraints. Accordingly , 𝑛 𝑞 = 9 − 𝑛 g = 6 generalized coordinates hav e to be chosen to describe the system unambiguously; let these be  𝑥 G , 𝑦 G , 𝜗, 𝜓 , 𝜑, 𝑟  . (32) Moreov er , 𝑛 𝜎 = 𝑛 𝑞 − 𝑛 k = 4 pseudovelocities ha ve to be cho- sen; let these be the components of the angular velocity 𝝎 of the wheel (see ( 73 ) in Appendix C ) and the axle directional velocity component of the point mass: 𝜔 1 ( 𝑡 ) : = ¤ 𝜗, 𝜔 2 ( 𝑡 ) : = ¤ 𝜓 sin 𝜗 + ¤ 𝜑, 𝜔 3 ( 𝑡 ) : = ¤ 𝜓 cos 𝜗 , 𝜎 ( 𝑡 ) : = ¤ 𝑟 − ¤ 𝜗 𝑅 . (33) The Appellian approach yields the equations of motion                                                      ¤ 𝜔 1 = 1 5 𝑚 𝑅 2 + 4 𝑚 0 𝑟 2  − 4 𝜔 2 1 𝑚 0 𝑅𝑟 − 𝜔 2 3 ( 𝑚 𝑅 2 + 4 𝑚 0 𝑟 2 ) tan 𝜗 − 8 𝜔 1 𝜎 𝑚 0 𝑟 + 𝜔 2 𝜔 3 ( 6 𝑚 𝑅 2 + 4 𝑚 0 𝑅𝑟 tan 𝜗 ) − 4 𝑚 0 𝑔𝑟 cos 𝜗 + 4 𝑚 𝑔 𝑅 sin 𝜗 + 4 𝑅 𝑢  , ¤ 𝜔 2 = 2 3 𝑚 𝑅 2 + 2 𝑚 0 𝑅 2 + 12 𝑚 0 𝑟 2  − 2 𝜔 1 𝜔 2 𝑚 0 𝑅𝑟 + 𝜔 1 𝜔 3 ( 𝑚 0 𝑅 2 − 𝑚 𝑅 2 − 4 𝑚 0 𝑟 2 ) + 2 𝜔 3 𝜎 𝑚 0 𝑅  , ¤ 𝜔 3 = 1 3 𝑚 𝑅 2 + 2 𝑚 0 𝑅 2 + 12 𝑚 0 𝑟 2  − 𝜔 1 𝜔 2 ( 6 𝑚 𝑅 2 + 4 𝑚 0 𝑅 2 ) + 𝜔 1 𝜔 3  ( 3 𝑚 𝑅 2 + 2 𝑚 0 𝑅 2 + 12 𝑚 0 𝑟 2 ) tan 𝜗  − 24 𝜔 3 𝜎 𝑚 0 𝑟 − 20 𝑚 0 𝑅𝑟  , ¤ 𝜗 = 𝜔 1 , ¤ 𝜎 = 𝜔 2 1 𝑟 + 𝜔 2 3 𝑟 − 𝜔 2 𝜔 3 𝑅 − 𝑔 sin 𝜗 + 1 𝑚 0 𝑢 , ¤ 𝑟 = 𝜎 + 𝜔 1 𝑅 ,                ¤ 𝜓 = 𝜔 3 1 cos 𝜗 , ¤ 𝜑 = 𝜔 2 − 𝜔 3 tan 𝜗 , ¤ 𝑥 G = 𝜔 1 𝑅 sin 𝜓 cos 𝜗 + 𝜔 2 𝑅 cos 𝜓 , ¤ 𝑦 G = − 𝜔 1 𝑅 cos 𝜓 cos 𝜗 + 𝜔 2 𝑅 sin 𝜓 , (34) which is a 𝑛 = 9 − 𝑛 g − 𝑛 k / 2 = 5 DoF nonholonomic mechan- ical system. The equations of motion ( 34 ) of the unicycle can be divided into essential dynamics (first six equations) and hidden dynamics (remaining four equations) [ 22 ]. Observ e that the hidden dynamics of the unicycle ( 34 ) are identical to that of the rolling wheel ( 6 ). Moreover , considering 𝑚 0 = 0 and 𝑢 = 0 the essential dynamics in ( 34 ) simplify to that in ( 6 ). The dynamical model ( 34 ) of the unicycle is formulated as a control affine system ¤ x = 𝑓 ( x ) + 𝑔 ( x ) 𝑢 with the state v ector x =  𝜔 1 𝜔 2 𝜔 3 𝜗 𝜎 𝑟 𝜓 𝜑 𝑥 G 𝑦 G  T , (35) and control input 𝑢 ; so this model is ready for control design. 7 B. Steady states Considering 𝑢 ≡ 0 (i.e., the point mass freely moves along the axle), the unicycle model ( 34 ) possesses the steady state motion with essential dynamics 𝜔 1 ( 𝑡 ) ≡ 𝜔 1 ∗ , 𝜔 2 ( 𝑡 ) ≡ 𝜔 2 ∗ , 𝜔 3 ( 𝑡 ) ≡ 𝜔 3 ∗ , 𝜗 ( 𝑡 ) ≡ 𝜗 ∗ , 𝜎 ( 𝑡 ) ≡ 𝜎 ∗ , 𝑟 ( 𝑡 ) ≡ 𝑟 ∗ . (36) Substituting this into the first six equations of ( 34 ) yields 𝜔 2 ∗ 𝜔 3 ∗  6 𝑚 𝑅 2 + 4 𝑚 0 𝑅𝑟 ∗ tan 𝜗 ∗  − 4 𝑚 0 𝑔𝑟 ∗ cos 𝜗 ∗ − 𝜔 2 3 ∗  𝑚 𝑅 2 + 4 𝑚 0 𝑟 2 ∗  tan 𝜗 ∗ + 4 𝑚 𝑔 𝑅 sin 𝜗 ∗ = 0 , 𝜔 2 3 ∗ 𝑟 ∗ − 𝜔 2 ∗ 𝜔 3 ∗ 𝑅 − 𝑔 sin 𝜗 ∗ = 0 , (37) while the hidden motion of the unicycle is the same as the hidden motion of the rolling wheel described in ( 10 ) and ( 11 ). Using ( 33 ), the relations ( 37 ) can be reformulated using the generalized velocities ( ¤ 𝜓 ∗ , ¤ 𝜗 ∗ , ¤ 𝜑 ∗ , ¤ 𝑟 ∗ ): ¤ 𝜓 2 ∗  4 𝑚 0 𝑅𝑟 ∗ sin 2 𝜗 ∗ + ( 5 𝑚 𝑅 2 − 4 𝑚 0 𝑟 2 ∗ ) sin 𝜗 ∗ cos 𝜗 ∗  + ¤ 𝜓 ∗ ¤ 𝜑 ∗  6 𝑚 𝑅 2 cos 𝜗 ∗ + 4 𝑚 0 𝑅𝑟 ∗ sin 𝜗 ∗  − 4 𝑚 0 𝑔𝑟 ∗ cos 𝜗 ∗ + 𝑚 𝑔 𝑅 sin 𝜗 ∗ = 0 , ¤ 𝜓 2 ∗  𝑟 ∗ cos 2 𝜗 ∗ − 𝑅 sin 𝜗 ∗ cos 𝜗 ∗  − ¤ 𝜓 ∗ ¤ 𝜑 ∗ 𝑅 cos 𝜗 ∗ − 𝑔 sin 𝜗 ∗ = 0 , (38) Then one may express the pitch rate ¤ 𝜑 ∗ and point mass position 𝑟 ∗ as a function of the tilt angle 𝜗 ∗ and the yaw rate ¤ 𝜓 ∗ : ¤ 𝜑 ∗ = tan 𝜗 ∗ ¤ 𝜓 ∗ 𝑅 ( 6 ¤ 𝜓 2 ∗ 𝑚 𝑅 cos 3 𝜗 ∗ − 4 𝑚 0 𝑔 )  − 5 ¤ 𝜓 4 ∗ 𝑚 𝑅 2 cos 4 𝜗 ∗ + 4 ¤ 𝜓 2 ∗ 𝑔 𝑅 ( 𝑚 0 cos 𝜗 ∗ − 𝑚 cos 3 𝜗 ∗ ) + 4 𝑚 0 𝑔 2  , 𝑟 ∗ = 𝑚 𝑅 ( ¤ 𝜓 2 ∗ 𝑅 cos 𝜗 ∗ + 2 𝑔 ) sin 𝜗 ∗ cos 𝜗 ∗ 6 ¤ 𝜓 2 ∗ 𝑚 𝑅 cos 3 𝜗 ∗ − 4 𝑚 0 𝑔 , (39) if ¤ 𝜓 ∗ ≠ 0 , ¤ 𝜓 ∗ ≠ ±  2 𝑚 0 𝑔 3 𝑚 𝑅 cos 3 𝜗 ∗ , (40) holds. For the non-generic cases, when ( 40 ) does not hold, the steady state motions are explained belo w . When the yaw rate is zero ¤ 𝜓 ∗ = 0 , straight rolling occurs and ( 38 ) simplify to 4 𝑚 𝑅 sin 𝜗 ∗ − 4 𝑚 0 𝑟 ∗ cos 𝜗 ∗ = 0 , 𝑔 sin 𝜗 ∗ = 0 . (41) These are only satisfied for zero tilt angle 𝜗 ∗ = 0 and centered point mass 𝑟 ∗ = 0 , independent of the pitch rate ¤ 𝜑 ∗ . In case of ¤ 𝜓 ∗ = ±  2 𝑚 0 𝑔 / ( 3 𝑚 𝑅 cos 3 𝜗 ∗ ) , ( 38 ) still results in zero tilt angle 𝜗 ∗ = 0 and ¤ 𝜓 ∗ 1 , 2 = ±  2 𝑚 0 𝑔 3 𝑚 𝑅 , ¤ 𝜑 ∗ 1 , 2 = ±  2 𝑚 0 𝑔 3 𝑚 𝑅 3 𝑟 ∗ , (42) which describe a turning-rolling steady state motion with a non-tilted wheel. This is a new behavior of the unicycle which does not exist for rolling wheel, as the latter must be tilted in order to have a turning-rolling steady state motion, cf. ( 12 ). Note that, this steady state can only occur for the specific yaw rates ¤ 𝜓 ∗ 1 , 2 depending only on physical parameters, while the corresponding pitch rates ¤ 𝜑 ∗ 1 , 2 also depend only on the point mass position 𝑟 ∗ . The non-tilted turning steady state becomes a spinning steady state ( ¤ 𝜑 ∗ 1 , 2 = 0 ) if the point mass is at the center 𝑟 ∗ = 0 . The other spinning steady state can be obtained by substi- tuting zero pitch rate ¤ 𝜑 ∗ = 0 into ( 38 ), which yields ¤ 𝜓 2 ∗  ( 5 𝑚 𝑅 2 − 4 𝑚 0 𝑟 2 ∗ ) sin 𝜗 ∗ cos 𝜗 ∗ + 4 𝑚 0 𝑅𝑟 ∗ sin 2 𝜗 ∗  + 4 𝑚 𝑔 𝑅 sin 𝜗 ∗ − 4 𝑚 0 𝑔𝑟 ∗ cos 𝜗 ∗ = 0 , ¤ 𝜓 2 ∗  𝑟 ∗ cos 2 𝜗 ∗ − 𝑅 sin 𝜗 ∗ cos 𝜗 ∗  − 𝑔 sin 𝜗 ∗ = 0 . (43) Similar to the spinning wheel, such steady state motion exists for zero tilt angle ( 𝜗 ∗ = 0 ) when the point mass is at the center 𝑟 ∗ = 0 . Ho wever , in case of the unicycle, spinning type steady state motion also exists with non-zero tilt angle ( 𝜗 ∗ ≠ 0 ) when ¤ 𝜓 ∗ = ±  4 𝑚 0 𝑔𝑟 ∗ cos 𝜗 ∗ − 4 𝑚 𝑔 𝑅 sin 𝜗 ∗ ( 5 𝑚 𝑅 2 − 4 𝑚 0 𝑟 2 ∗ ) sin 𝜗 ∗ cos 𝜗 ∗ + 4 𝑚 0 𝑅𝑟 ∗ sin 2 𝜗 ∗ , 𝑟 ∗ = 𝑅 tan 𝜗 ∗ 2 𝑚 0  𝑚 cos 2 𝜗 ∗ + 𝑚 0 (44) +  𝑚 2 cos 4 𝜗 ∗ + 3 𝑚 𝑚 0 cos 2 𝜗 ∗ + 𝑚 2 0  . W e refer to this as tilted spinning, which is also a ne w behavior of the unicycle compared to the rolling wheel. According to ( 44 ), the point mass is always above the wheel center since sgn 𝑟 ∗ = sgn 𝜗 ∗ or 𝑟 ∗ 𝜗 ∗ > 0 ; cf. Figure 4 . In vestigating special cases like 𝜗 ∗ = 0 or 𝑟 ∗ = 0 leads to steady states already discussed above. One limiting factor of our unicycle model is that the point mass may slide below the ground which is physically not feasible. T o exclude these cases, one can calculate the 𝑧 coordinate of the position vector r A = r G + r GA in frame F 0 and require it to be positi ve. This yields 𝑅 cos 𝜗 ∗ + 𝑟 ∗ sin 𝜗 ∗ > 0 . All steady state motions of the unicycle are summarized in Figure 5 in the plane of the steady state tilt angle 𝜗 ∗ and the yaw rate ¤ 𝜓 ∗ . The straight rolling, i.e., 𝜗 ∗ = 0 and ¤ 𝜓 ∗ = 0 , is marked by a black square ■ ; in this case the steady state pitch rate ¤ 𝜑 ∗ may be arbitrary and the point mass is at wheel center 𝑟 ∗ = 0 . The non-tilted turning, i.e., 𝜗 ∗ = 0 and ¤ 𝜓 ∗ = ¤ 𝜓 ∗ 1 , 2 , is marked by black dots ; in this case the pitch rate ¤ 𝜑 ∗ and the point mass position 𝑟 ∗ linked as in ( 42 ). The regular (non-tilted) spinning steady state motions, i.e., 𝜗 ∗ = 0 and ¤ 𝜑 ∗ = 0 , are sho wn by the solid black line ( — ); in this case the yaw rate ¤ 𝜓 ∗ may be arbitrary and the point mass is at the wheel center 𝑟 ∗ = 0 . The tilted spinning steady states, i.e., ¤ 𝜑 ∗ = 0 and 𝜗 ∗ ≠ 0 , are shown by dashed black lines ( – – ); the corresponding yaw rate ¤ 𝜓 ∗ and point mass position 𝑟 ∗ are expressed by ( 44 ). The general turning-rolling steady states are di vided into separate areas. The light blue ■ ■ and light purple ■ ■ areas show the physically feasible steady states when the point mass is placed below 0 < 𝑧 A < 𝑧 G or above 𝑧 A > 𝑧 G the wheel center point G , respectively . The white area < = represents the physically unfeasible steady states when the point mass is at or below the ground lev el, that is, 𝑧 A ≤ 0 . Note that Figure 5 was constructed using the numerical parameters in T able I . Using dif ferent parameters, like dif ferent mass ratios 𝑚 0 / 𝑚 , the main structure of the steady states 8 Fig. 5. Summary of the steady state motions of the unicycle remain qualitativ ely similar , but the shape and size of the parameter regions related to the non-feasible steady states may be different. The steady states hav e been categorized with considering zero input force ( 𝑢 ≡ 0 ). Assuming a nonzero constant input ( 𝑢 ≡ 𝑢 ∗ ) may result in further steady states. For example, a straight rolling case may occur with tilted wheel when the mass is held ‘above’ the wheel center . Studying such steady states are out of the scope of the present study . C. Stability of steady states Linear stability analysis is performed to determine the stability of the previously discussed steady state motions of the unicycle. Linearizing ( 34 ) around a steady state x ∗ leads to the form ¤ ˜ x = A ˜ x where ˜ x = x − x ∗ denotes the perturbed states and state matrix A = 𝑓 ′ ( x ∗ ) can be written as A =                     0 𝐴 12 𝐴 13 𝐴 14 0 𝐴 16 0 0 0 0 𝐴 21 0 0 0 𝐴 25 0 0 0 0 0 𝐴 31 0 0 0 𝐴 35 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 𝐴 52 𝐴 53 𝐴 54 0 𝐴 56 0 0 0 0 𝐴 61 0 0 0 1 0 0 0 0 0 0 0 𝐴 73 𝐴 74 0 0 0 0 0 0 0 1 𝐴 83 𝐴 84 0 0 0 0 0 0 𝐴 91 𝐴 92 0 0 0 0 𝐴 97 0 0 0 𝐴 101 𝐴 102 0 0 0 0 𝐴 107 0 0 0                     , (45) whose elements are giv en in Appendix D . Note that the steady state tilt angle 𝜗 ∗ , ya w rate ¤ 𝜓 ∗ , pitch rate ¤ 𝜑 ∗ and point mass position 𝑟 ∗ in the elements of A are not independent of each other , they must satisfy ( 38 ). Therefore, each of the dif ferent steady state motions (turning-rolling, straight rolling, non- tilted turning) must be analyzed separately . In this paper we focus on the straight rolling and turning-rolling motions while the non-tilted turning is left for future research. Also note that, the spinning steady state is a special case of the more general turning-rolling steady state; since it has minor practical relev ance its further analysis is omitted. The characteristic equation det ( 𝜆 I − ˆ A ) = 0 , where ˆ A is the top left 6 by 6 part of A in ( 45 ), takes the form  𝜆 4 − 𝜆 2 ( 𝐴 12 𝐴 21 + 𝐴 13 𝐴 31 + 𝐴 14 + 𝐴 16 𝐴 61 + 𝐴 25 𝐴 52 + 𝐴 35 𝐴 53 + 𝐴 56 ) + 𝐴 12 𝐴 21 𝐴 35 𝐴 53 − 𝐴 12 𝐴 25 𝐴 31 𝐴 53 − 𝐴 12 𝐴 25 𝐴 56 𝐴 61 − 𝐴 13 𝐴 21 𝐴 35 𝐴 52 + 𝐴 13 𝐴 25 𝐴 31 𝐴 52 − 𝐴 13 𝐴 35 𝐴 56 𝐴 61 + 𝐴 16 𝐴 25 𝐴 52 𝐴 61 + 𝐴 16 𝐴 35 𝐴 53 𝐴 61 + 𝐴 12 𝐴 21 𝐴 56 − 𝐴 12 𝐴 25 𝐴 54 + 𝐴 13 𝐴 31 𝐴 56 − 𝐴 13 𝐴 35 𝐴 54 + 𝐴 14 𝐴 25 𝐴 52 + 𝐴 14 𝐴 35 𝐴 53 − 𝐴 16 𝐴 21 𝐴 52 − 𝐴 16 𝐴 31 𝐴 53 + 𝐴 14 𝐴 56 − 𝐴 16 𝐴 54  𝜆 2 = 0 . (46) 1) T urning-r olling steady state: When calculating the char- acteristic roots of ( 46 ) the steady state pitch rate ¤ 𝜑 ∗ and point mass position 𝑟 ∗ must be substituted according to ( 39 ). Due to the highly complicated expressions, numerical ev al- uation is used to obtain the characteristic roots and deter- mine the stability of the turning-rolling steady states. The results are shown in Figure 6 for se veral mass ratios with 𝑚 0 = 0 , 1 , 2 , 5 , 10 , 15 , 20 , 30 kg, while the other physical pa- rameters are presented in T able I . The same coloring scheme is used as in Fig. 3 , namely , green ( ■ ■ ) and red ( ■ ■ ) areas rep- resent stable and unstable steady state solutions, respectiv ely , while the white area ( < = ) represents the ph ysically unfeasible steady states with point mass at or below the ground lev el, cf. Fig. 5 . The case with 𝑚 0 = 0 is identical to the rolling wheel in Figure 3 (c). By the addition the point mass 𝑚 0 > 0 , the stability properties change drastically since the number of degrees of freedom is increased, additional types of equilibria are born (the tilted spinning and the vertical turning), the point mass may be abov e or below the wheel, and it ev en may hit the ground, as w as discussed in Section III-B and Figure 5 . W ith small point masses (such as 𝑚 0 = 1 , 2 kg) the vertical spinning ( 𝜗 ∗ = 0 ) becomes unstable as the point mass is pushed away from the spinning wheel by the centrifugal force. W ith large point masses (such as 𝑚 0 = 15 , 20 , 30 kg) a stable yaw rate ¤ 𝜓 ∗ band appears, but the spinning steady states remains unstable for lar ge enough yaw rates. 2) Straight r olling: The steady state is given as ( 23 ) and all the other states become zero, considering that the initial values 𝑥 0 , 𝑦 0 and ¤ 𝜓 0 are zeros. The state matrix ( 45 ) simplifies to A =                     0 0 6 5 ¤ 𝜑 ∗ 4 𝑔 5 𝑅 0 − 4 𝑚 0 𝑔 5 𝑚 𝑅 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 − 2 ¤ 𝜑 ∗ 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 − 𝑅 ¤ 𝜑 ∗ − 𝑔 0 0 0 0 0 0 𝑅 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 𝑅 0 0 0 0 0 0 0 0 − 𝑅 0 0 0 0 0 𝑅 ¤ 𝜑 ∗ 0 0 0                     . (47) When comparing this with ( 24 ) obtained for the rolling wheel, there are two extra rows and columns, 5 and 6, related to the new states 𝑟 and 𝜎 . The characteristic polynomial ( 46 ) reduces to  𝑎𝜆 4 + 𝑏𝜆 2 + 𝑐  𝜆 2 = 0 , (48) 9 Fig. 6. Stability of steady state motions of the unicycle for various 𝑚 0 values: stable and unstable steady states are denoted by green and red, respectively , while white denotes the unfeasible steady states when the point mass hits the ground. 10 where 𝑎 = 5 𝑚 𝑅 2 , 𝑏 = 12 ¤ 𝜑 2 ∗ 𝑚 𝑅 2 + 4 ( 𝑚 0 − 𝑚 ) 𝑔 𝑅, 𝑐 = 4 𝑚 0 𝑔 ( 2 ¤ 𝜑 2 ∗ 𝑅 − 𝑔 ) . (49) Applying the same ar guments as in Appendix B , the straight rolling is said to be stable if the coef ficients 𝑎 , 𝑏 and 𝑐 hav e the same signs. While 𝑎 > 0 always holds, 𝑏 > 0 and 𝑐 > 0 lead to | ¤ 𝜑 ∗ | >  ( 𝑚 − 𝑚 0 ) 𝑔 3 𝑚 𝑅 and | ¤ 𝜑 ∗ | >  𝑔 2 𝑅 , (50) respectiv ely . The second condition is always stronger , so the critical pitch rate of the unicycle is ¤ 𝜑 crit =  𝑔 2 𝑅 . (51) That is, the straight rolling of the uncontrolled unicycle is stable if | ¤ 𝜑 ∗ | > ¤ 𝜑 crit ; it is unstable if | ¤ 𝜑 ∗ | < ¤ 𝜑 crit . Note that, the critical pitch rate ( 51 ) is independent of the masses 𝑚 and 𝑚 0 . Also, it is larger than the critical pitch rate of the rolling disc ( 26 ). The dependency of the characteristic roots on the steady state pitch rate ¤ 𝜑 ∗ are summarized in Figure 7 for both the rolling wheel ( 𝑚 0 = 0 ) and the unicycle − 5 0 5 Real part of 𝜆 𝑖 𝜆 1 , 2 ( 𝑚 0 = 0 . 1 kg ) 𝜆 3 , 4 ( 𝑚 0 = 0 . 1 kg ) 𝜆 1 , 2 ( 𝑚 0 = 5 kg ) 𝜆 3 , 4 ( 𝑚 0 = 5 kg ) ¤ 𝜑 unicycle crit = p 𝑔 / 2 𝑅 ¤ 𝜑 wheel crit = p 𝑔 / 3 𝑅 𝜆 1 , 2 ( 𝑚 0 = 0 kg ) 0 1 2 3 4 5 6 Steady state pitch rate ¤ 𝜑 ∗ [rad/s] − 10 − 5 0 5 10 Imaginary part of 𝜆 𝑖 Fig. 7. Characteristic roots as a funtion the steady state pitch rate ¤ 𝜑 ∗ for straight rolling comparing the cases of rolling wheel and unicycle ( 𝑚 0 > 0 ). Observe that in case of the rolling wheel the charac- teristic roots 𝜆 1 , 2 (thick gray curves) determine the criticality . Contrarily , in case of the unicycle (e ven for small v alues of 𝑚 0 !), the criticality is determined by the characteristic roots 𝜆 3 , 4 (orange curves) rather than by the characteristic roots 𝜆 1 , 2 (blue curv es). This resolves the seemingly contradictory issue that the critical pitch rate ‘jumps’ from ( 26 ) to ( 51 ) when 𝑚 0 becomes positiv e. On the other hand, this is not surprising in view of the different de grees of freedom in the two systems. 3) Spinning on the spot: The stability analysis of this steady state is covered by the analysis of the turning-rolling motion. Since the motion has minor practical relev ance for unicycles further discussions are omitted here. I V . S T E E R I N G C O N T R O L O F T H E U N I C Y C LE The straight rolling steady state ( 23 ) is used as the basis for designing the steering controller of the unicycle. For simplicity , it is assumed that the unicycle has an initial pitch rate ¤ 𝜑 0 = ¤ 𝜑 ∗ and initially it trav els along the 𝑥 0 axis. Abov e we showed that the straight rolling steady state might be unstable or stable without control. Here we demonstrate that this motion can be made asymptotically stable with appropriate feedback control. The linear state space model of the unicycle assumes the form ¤ ˜ x = A ˜ x + B 𝑢 where ˜ x = x − x ∗ . The state matrix A = 𝑓 ′ ( x ∗ ) is giv en in ( 47 ) and it depends on the steady state pitch rate ¤ 𝜑 ∗ . The input matrix B = 𝑔 ( x ∗ ) is given by B = h 4 5 𝑚 𝑅 0 0 0 1 𝑚 0 0 0 0 0 0 i T . (52) A. Contr ollability The controllability matrix of the unicycle is obtained as M c =  B AB A 2 B . . . A 9 B  . (53) Substituting ( 47 ) and ( 52 ) we obtain M c =                     𝑀 11 0 𝑀 13 0 𝑀 15 0 𝑀 17 0 𝑀 19 0 0 0 0 0 0 0 0 0 0 0 0 𝑐 𝑀 42 0 𝑐 𝑀 44 0 𝑐 𝑀 46 0 𝑐 𝑀 48 0 𝑐 𝑀 410 0 𝑀 42 0 𝑀 44 0 𝑀 46 0 𝑀 48 0 𝑀 410 𝑀 51 0 𝑀 53 0 𝑀 55 0 𝑀 57 0 𝑀 59 0 0 𝑀 62 0 𝑀 64 0 𝑀 66 0 𝑀 68 0 𝑀 610 0 0 𝑀 73 0 𝑀 75 0 𝑀 77 0 𝑀 79 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 𝑀 102 0 𝑀 104 0 𝑀 106 0 𝑀 108 0 𝑀 1010                     , (54) with 𝑐 = − 2 ¤ 𝜑 ∗ . Notice that ro ws 2, 8 and 9 are full of zeros while rows 3 and 4 are linearly dependent. Thus, rank M c = 6 and the controllability matrix is not full row rank. That is, at the linear le vel, the unicycle is not fully controllable with the single input 𝑢 . When spelling out the 3 rd and 4 th equations in the linear system ¤ ˜ x = A ˜ x + B 𝑢 one obtains ¤ 𝜔 3 = − 2 ¤ 𝜑 ∗ 𝜔 1 , ¤ 𝜗 = 𝜔 1 (55) 11 which lead to 𝜔 3 = − 2 ¤ 𝜑 ∗ 𝜗 . (56) This means that the yaw rate ¤ 𝜓 is linearly proportional to the tilt angle 𝜗 , since 𝜔 3 ≈ ¤ 𝜓 for small tilt angles; cf. ( 33 ). That is, to steer the unicycle it is necessary to tilt it accordingly . Note that one may steer the unicycle with a specific yaw rate ¤ 𝜓 ∗ in ( 42 ) ev en with zero tilt angle, but the model must be linearized around the non-tilted turning steady state of the unicycle. This is left for future research. B. Maneuvering In this study , two maneuvers are considered for the steering control of the unicycle: a lane change and a 90 ◦ right turn, while assuming that the unic ycle initially trav els along the 𝑥 0 axis. Considering the steady state velocity 𝑣 G ∗ = 𝑅 ¤ 𝜑 ∗ of the center of gravity leads to the critical speed 𝑣 G , crit =  𝑅 𝑔 / 2 for the critical pitch rate ¤ 𝜑 crit in ( 51 ). For the parameters in T able I , we obtain 𝑣 G , crit ≈ 1 . 21 m/s. Thus, the lane change and right turn maneuvers are in vestigated for two speeds: 𝑣 G1 = 1 m / s which is below the critical speed, and 𝑣 G2 = 5 m / s which is abov e the critical speed. When designing the controllers to execute these maneuvers, the non-reachable states 𝜔 2 , 𝜔 3 , 𝜑 , 𝑥 G are omitted; cf. ( 35 ) and ( 54 ). Then, according to the state matrix ( 47 ) and input matrix ( 52 ) the reduced system becomes ¤ 𝜔 1 =  4 𝑔 5 𝑅 − 12 5 ¤ 𝜑 2 ∗  𝜗 − 4 𝑚 0 𝑔 5 𝑚 𝑅 2 𝑟 + 4 5 𝑚 𝑅 𝑢 , ¤ 𝜗 = 𝜔 1 , ¤ 𝜎 =  2 ¤ 𝜑 2 ∗ 𝑅 − 𝑔  𝜗 + 1 𝑚 0 𝑢 , ¤ 𝑟 = 𝜎 + 𝑅 𝜔 1 , ¤ 𝜓 = − 2 ¤ 𝜑 ∗ 𝜗, ¤ 𝑦 G = − 𝑅 𝜔 1 + ¤ 𝜑 ∗ 𝑅 𝜓 , (57) where the tildes abov e the states (representing perturbations) are omitted since the equilibrium values are zeros for all these states. W e formally define the states in this reduced system as outputs of the full linear system ¤ ˜ x = A ˜ x + B 𝑢 when designing the controllers for the lane change and turning maneuv ers. 1) Lane change: For the lane change maneuver the state 𝑦 G can be utilized to plan the motion since the unicycle should mov e parallel to the 𝑥 0 axis at the beginning and at the end of the maneuver . According to this, we consider the output y = Cx : =  𝜔 1 𝜗 𝜎 𝑟 𝜓 𝑦 G  T , (58) (cf. ( 57 )) with C =             1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1             . (59) The output controllability matrix can be constructed as M oc =  CB C AB CA 2 B . . . CA 9 B  , (60) T ABLE II C O NT RO L G AI N S F O R T H E L A N E C H AN G E M A N EU V E R Below critical speed 𝑣 G1 < 𝑣 G , crit 𝐷 𝜗 − 1407 . 64 Ns 𝑃 𝜗 − 7637 . 29 N 𝐷 𝑟 2116 . 86 Ns/m 𝑃 𝑟 11942 . 04 N/m 𝑃 𝜓 3382 . 02 N 𝑃 𝑦 4509 . 36 N/m Above critical speed 𝑣 G2 > 𝑣 G , crit 𝐷 𝜗 105 . 36 Ns 𝑃 𝜗 777 . 28 N 𝐷 𝑟 99 . 52 Ns/m 𝑃 𝑟 405 . 60 N/m 𝑃 𝜓 676 . 4 N 𝑃 𝑦 180 . 37 N/m T ABLE III C O NT RO L G AI N S F O R T H E T U R NI N G M A N E UV E R Below critical speed 𝑣 G1 < 𝑣 G , crit 𝐷 𝜗 − 511 . 89 Ns 𝑃 𝜗 − 2776 . 88 N 𝐷 𝑟 882 . 53 Ns/m 𝑃 𝑟 4881 . 84 N/m 𝑃 𝜓 536 . 67 N Above critical speed 𝑣 G2 > 𝑣 G , crit 𝐷 𝜗 118 . 88 Ns 𝑃 𝜗 − 128 . 32 N 𝐷 𝑟 41 . 49 Ns/m 𝑃 𝑟 73 . 69 N/m 𝑃 𝜓 112 . 73 N which has a maximal row rank, rank M oc = 6 . That is, the unicycle is output controllable. Then, we apply the linear output feedback law 𝑢 : = − K ( y − y des ) , (61) with control gains K = [ 𝐷 𝜗 𝑃 𝜗 𝐷 𝑟 𝑃 𝑟 𝑃 𝜓 𝑃 𝑦 ] and y des = [ 0 0 0 0 0 𝑦 des ( 𝑡 ) ] T , (62) where 𝑦 des ( 𝑡 ) encodes the desired trajectory in terms of 𝑦 G . Here we consider 𝑦 des ( 𝑡 ) =          0 , 0 ≤ 𝑡 < 2 , 𝐿 2  cos  𝜋 5 ( 𝑡 − 2 )  − 1  , 2 ≤ 𝑡 < 7 , − 𝐿 , 7 ≤ 𝑡 < 10 , (63) and we use 𝐿 = 2 . 5 m for speed 𝑣 G1 = 1 m / s and 𝐿 = 10 m for speed 𝑣 G2 = 2 m / s . The desired trajectories are depicted as a black dashed-dotted curves in the ( 𝑥 0 , 𝑦 0 ) -plane in the bottom left panel in Fig. 8 . For simplicity , we consider the desired yaw angle 𝜓 des ( 𝑡 ) ≡ 0 in ( 62 ). More sophisticated reference trajectories might be constructed utilizing ( 38 ). Note that the proposed linear controller does not ensure the global stability of the nonlinear system. Nonlinear inv es- tigations are outside of the scope of the present study . Also, maneuvers with substantial non-zero ¤ 𝜓 may make the unicycle deviate significantly from straight rolling steady state which may cause large deviation in the observed behavior . The characteristic equation of the closed-loop system is det  𝜆 I − ( A − BK C )  = 0 . (64) By selecting the control gains in K in ( 61 ) appropriately one may place the non-zero characteristic roots to the left half complex plane and guarantee stability for the closed-loop system. In particular , using the control gains in T able II results in 𝜆 𝑘 = − 8 s − 1 , 𝑘 = 1 , . . . , 6 . The performance of the closed loop system is demonstrated by numerical simulations in the left column in Fig. 8 . These were carried out using the full nonlinear equations ( 34 ) ap- plying the control law ( 61 ) with reference trajectory ( 63 ) and control gains in T able II . The time ev olution of tilt angle 𝜗 , 12 − 5 ◦ − 2 . 5 ◦ 0 ◦ 2 . 5 ◦ 5 ◦ 𝜗 Lane c hange − 15 − 7 . 5 0 7 . 5 15 𝑟 [mm] − 90 ◦ − 67 . 5 ◦ − 45 ◦ − 22 . 5 ◦ 0 ◦ 𝜓 0 1 2 3 4 5 6 7 8 9 10 Time 𝑡 [s] − 3 − 1 . 5 0 1 . 5 3 𝑢 [N] 0 5 10 15 20 25 30 35 40 45 50 𝑥 G [m] − 30 − 25 − 20 − 15 − 10 − 5 0 𝑦 G [m] 0 2 4 6 8 10 − 2 0 − 2 . 5 ◦ − 1 . 25 ◦ 0 ◦ 1 . 25 ◦ 2 . 5 ◦ − 300 − 150 0 150 300 − 15 − 7 . 5 0 7 . 5 15 − 5 ◦ − 2 . 5 ◦ 0 ◦ 2 . 5 ◦ 5 ◦ 𝜗 90 ◦ right turn − 15 − 7 . 5 0 7 . 5 15 𝑟 [mm] − 90 ◦ − 67 . 5 ◦ − 45 ◦ − 22 . 5 ◦ 0 ◦ 𝜓 0 1 2 3 4 5 6 7 8 9 10 Time 𝑡 [s] − 3 − 1 . 5 0 1 . 5 3 𝑢 [N] 0 5 10 15 20 25 30 35 40 45 50 𝑥 G [m] − 30 − 25 − 20 − 15 − 10 − 5 0 𝑦 G [m] 0 2 4 6 − 6 − 4 − 2 0 − 2 . 5 ◦ − 1 . 25 ◦ 0 ◦ 1 . 25 ◦ 2 . 5 ◦ − 300 − 150 0 150 300 − 15 − 7 . 5 0 7 . 5 15 desired tra jector y (v G1 = 1 m/s) ( v G2 = 5 m/s) Fig. 8. Maneuvers: lane change (left column) and 90 ◦ right turn (right column) below the critical speed (solid blue) and above (dashed orange) point mass position 𝑟 , yaw angle 𝜓 , and the control input 𝑢 are depicted in the top four panels, while the bottom panel shows the movement of the wheel center G above the ( 𝑥 0 , 𝑦 0 ) - plane. The desired lane change maneuver is followed by the unicycle, and the controller successfully stabilizes the desired motion for both speeds. The gains 𝑃 𝜗 and 𝐷 𝜗 are highlighted in T able II as they hav e opposite signs abov e and below the critical speed. Belo w the critical speed ( 𝑣 G1 < 𝑣 G , crit ), the open loop dynamics is unstable, here the 𝑃 𝜗 and 𝐷 𝜗 gains correspond to positiv e stiffness and damping in the mechanical sense. In contrast, abo ve the critical speed ( 𝑣 G2 > 𝑣 G , crit ) the open loop dynamics is stable and the 𝑃 𝜗 and 𝐷 𝜗 gains correspond to ne gati ve stiffness and damping. This illustrates that executing a lane change maneuver above the critical speed takes more ef fort, since the straight-rolling motion of the open- loop systems is more stable. This can also be observed when comparing the required control input for the dif ferent speeds. Still the required control input remains small ( | 𝑢 | < 10 N) for both speeds. 2) Right turn: For the 90 ◦ right turn maneuver , the unicycle ev entually moves parallel to the 𝑦 0 axis and the state 𝑦 G becomes unreachable by the control input. Thus, the state 𝑦 G may not be used to plan motion in this case. Instead we use the yaw angle 𝜓 to construct the reference trajectory . Considering this, the output v ector y = Cx : =  𝜔 1 𝜗 𝜎 𝑟 𝜓  T , (65) with C =           1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0           (66) is defined; cf. ( 58 ) and ( 59 ). The output controllability matrix can be constructed as in ( 60 ); since it has maximal row rank, rank M oc = 5 , the unic ycle is output controllable. Here the linear output feedback control law of the form ( 61 ) is applied with control gains K = [ 𝐷 𝜗 𝑃 𝜗 𝐷 𝑟 𝑃 𝑟 𝑃 𝜓 ] and y des = [ 0 0 0 0 𝜓 des ( 𝑡 ) ] T , (67) where the yaw angle 𝜓 des ( 𝑡 ) =          0 , 0 ≤ 𝑡 < 2 , 𝜋 4  1 − cos  𝜋 5 ( 𝑡 − 2 )   , 2 ≤ 𝑡 < 7 , 𝜋 2 , 7 ≤ 𝑡 < 10 , (68) 13 encodes the desired trajectory . Such trajectory is depicted as a black dashed-dotted curve in Fig. 8 , (right column, third panel). Choosing the control gains as in T able III for the speeds 𝑣 G1 and 𝑣 G2 the non-zero characteristic roots of the closed-loop system are placed at 𝜆 𝑘 = − 8 s − 1 , 𝑘 = 1 , . . . , 5 . Since these are in the left half complex plane, the resulting closed-loop system is stable. The simulation results are shown in the right side of Figure 8 . These are carried out using the full nonlinear equations ( 34 ) applying the feedback law ( 61 ) with reference trajectory ( 67 ) and control gains in T able III . The right turns are successfully executed for both speeds. Comparing the control gains belo w and above the critical speed for the turning maneuver , only the deriv ative gain 𝐷 𝜗 changes sign as highlighted in T able III . When observing the time ev olution the states and the control input in Figure 8 , one may notice some high-frequency content which is more pronounced for abov e the critical speed (orange curves). These are due to the Coriolis force since the yaw rate is not zero while the point mass moves along the axle. These effects are compensated by the controller while e xecuting the maneuvers. V . C O N C L U S I O N In this study , modeling, analyses and control of an au- tonomous unicycle were considered. The most compact form of the equations of motion were deriv ed by means of the Appellian approach. The steady state motions of the rolling wheel were classified as cases of straight rolling, turning- rolling and spinning on the spot. The stability of these motions were determined by linear analysis and the critical angular velocities (above which the steady state motions are stable) were gi ven in closed form. It was also shown that a turning- rolling steady state is always stable above a critical tilt angle which is independent of the system parameters. Building upon the knowledge gained from the dynamics of the rolling wheel, the simplest possible control strategy is proposed for autonomous unic ycles by means of actuating the position of a point mass normal to the wheel plane, in order to accomplish steering maneuv ers such as lane changes and turns. The nonlinear equations of motions of the unicycle were gi ven in closed form. Apart from finding steady states akin to those of the rolling wheel, additional steady states were identified such as non-tilted turning and tilted spinning. For the open- loop unicycle, stability results were obtained for the straight rolling steady state and the stability of the turning-rolling steady states were determined semi-analytically , numerically . All these stability results rely on linearization, the analysis of the conserved quantities is reserved for future research. Utilizing the equations linearized about the straight rolling steady state, two steering controllers of partial state feedback were proposed for the unicycle to carry out a lane change and a 90 ◦ turn. The resulting linear systems were shown to be output controllable, and when applying the designed controllers to the nonlinear system the unicycle successfully completed the desired maneuvers as confirmed by numerical simulations. The dev eloped modeling, analysis and control framew ork may also allow more sophisticated motion planning as well as nonlinear control design which can take full advantage of the agility of unicycles. Such techniques, which are left for future research, may allow one to achiev e high lev el of maneuverability for autonomous unicycles and to provide steering assist for human-ridden unicycles. The introduced unicycle model incorporates engineering assumptions, which do not necessarily hold in all circum- stances. F or exampl e, sliding may occur violating the rolling assumption once the friction forces are not large enough. This will be analyzed in a separate study . Also, experiments are planned with the human-ridden and the autonomous unic ycles to verify the theoretical findings. A P P E N D I X A M O D E L D E R I V A T I O N F O R T H E R O L L I N G W H E E L V I A T H E A P P E L L I A N A P P RO AC H A vector resolved in frame F 2 can be transformed to the frame F 0 as □ F 0 = T 02 □ F 2 , (69) where the rotation matrix is T 02 =       cos 𝜓 − sin 𝜓 cos 𝜗 sin 𝜓 sin 𝜗 sin 𝜓 cos 𝜓 cos 𝜗 − cos 𝜓 sin 𝜗 0 sin 𝜗 cos 𝜗       . (70) The velocity of the wheel-ground contact point P can be expressed as v P = v G + 𝝎 × r GP , (71) where velocity of the center of gra vity G can resolved in the ground-fixed frame F 0 as v G = ¤ r G =  ¤ 𝑥 G ¤ 𝑦 G ¤ 𝑧 G  T F 0 , (72) while the angular velocity of the wheel can be expressed in the moving frame F 2 as 𝝎 =  ¤ 𝜗 ¤ 𝜓 sin 𝜗 + ¤ 𝜑 ¤ 𝜓 cos 𝜗  T F 2 , (73) using the tilt rate ¤ 𝜗 , the yaw rate ¤ 𝜓 , the pitch rate ¤ 𝜑 and the tilt angle 𝜗 . One may also express the position vector as r GP =  0 0 − 𝑅  T F 2 . (74) Then transforming the cross product in ( 71 ) to the ground- fixed frame F 0 with the help of ( 70 ), the rolling condition ( 1 ) results in the kinematic constraints ( 2 ) and the geometric constraint ( 3 ). The definitions of pseudov elocities are appropriate only if the generalized velocities ( ¤ 𝑥 G , ¤ 𝑦 G , ¤ 𝜓 , ¤ 𝜗 , ¤ 𝜑 ) can unambiguously be expressed by means of the pseudov elocities ( 𝜔 1 , 𝜔 2 , 𝜔 3 ) and generalized coordinates ( 𝑥 G , 𝑦 G , 𝜓 , 𝜗 , 𝜑 ). T o check this, we combine the kinematic constraints ( 2 ) and the definitions of pseudovelocities ( 5 ) into one linear system of equations:           0 0 𝜔 1 𝜔 2 𝜔 3           =           − 1 0 𝑅 sin 𝜗 cos 𝜓 𝑅 sin 𝜓 cos 𝜗 𝑅 cos 𝜓 0 − 1 𝑅 sin 𝜓 sin 𝜗 − 𝑅 cos 𝜓 cos 𝜗 𝑅 sin 𝜓 0 0 0 1 0 0 0 sin 𝜗 0 1 0 0 cos 𝜗 0 0           | {z } C           ¤ 𝑥 G ¤ 𝑦 G ¤ 𝜓 ¤ 𝜗 ¤ 𝜑           . (75) 14 This linear system can be solved if C is in vertible, that is, if its determinant is nonzero. Since det C = cos 𝜗 , the matrix C is singular when the wheel is horizontal ( 𝜗 = ± 𝜋 / 2 ), which is excluded from this analysis. Thus, the generalized velocities ( ¤ 𝑥 G , ¤ 𝑦 G , ¤ 𝜓 , ¤ 𝜗 , ¤ 𝜑 ) can be e xpressed with the pseudovelocities ( 𝜔 1 , 𝜔 2 , 𝜔 3 ) as ¤ 𝑥 G = 𝜔 1 𝑅 sin 𝜓 cos 𝜗 + 𝜔 2 𝑅 cos 𝜓 , ¤ 𝑦 G = − 𝜔 1 𝑅 cos 𝜓 cos 𝜗 + 𝜔 2 𝑅 sin 𝜓 , ¤ 𝜓 = 𝜔 3 1 cos 𝜗 , ¤ 𝜗 = 𝜔 1 , ¤ 𝜑 = 𝜔 2 − 𝜔 3 tan 𝜗 . (76) Also, the velocity of the center of gravity G ( 72 ) and the angular velocity ( 73 ) can be e xpressed in the terms of pseu- dov elocities ( 𝜔 1 , 𝜔 2 , 𝜔 3 ) using ( 76 ); these vectors hav e the most compact form when expressed in frame F 2 : v G =  ¤ 𝜓 𝑅 sin 𝜗 + ¤ 𝜑 𝑅 − ¤ 𝜗 𝑅 0  T F 2 ≡  𝜔 2 𝑅 − 𝜔 1 𝑅 0  T F 2 , 𝝎 =  ¤ 𝜗 ¤ 𝜓 sin 𝜗 + ¤ 𝜑 ¤ 𝜓 cos 𝜗  T F 2 ≡  𝜔 1 𝜔 2 𝜔 3  T F 2 . (77) The acceleration energy of a rigid body is defined as: 𝑆 = 1 2 𝑚 a 2 G + 1 2 𝜶 · J G 𝜶 + 𝜶 · ( 𝝎 × J G 𝝎 ) , (78) where a G = ¤ v G is the acceleration of the center of gravity , 𝜶 = ¤ 𝝎 is the angular acceleration and J G is the mass moment of inertia with respect to the center of gra vity . All these v ectors must be e xpressed based on the generalized coordinates ( 𝑥 G , 𝑦 G , 𝜓 , 𝜗 , 𝜑 ), pseudovelocities ( 𝜔 1 , 𝜔 2 , 𝜔 3 ) and pseudoaccel- erations ( ¤ 𝜔 1 , ¤ 𝜔 2 , ¤ 𝜔 3 ). Also, the vectors hav e the most compact form when expressed in frame F 2 : a G =       𝑅 ( ¤ 𝜔 2 + 𝜔 1 𝜔 3 ) − 𝑅 ( ¤ 𝜔 1 − 𝜔 2 𝜔 3 ) − 𝑅 ( 𝜔 2 1 + 𝜔 2 𝜔 3 tan 𝜗 )       F 2 , J G = 𝑚 𝑅 2 4       1 0 0 0 2 0 0 0 1       F 2 , 𝜶 =       ¤ 𝜔 1 − 𝜔 2 𝜔 3 + 𝜔 2 3 tan 𝜗 ¤ 𝜔 2 ¤ 𝜔 3 + 𝜔 1 𝜔 2 − 𝜔 1 𝜔 3 tan 𝜗       F 2 . (79) The calculation of the accelerations is not trivial in the moving frame F 2 ; still, this giv es the simplest possible algebraic form. These lead to the acceleration energy of the rolling disc: 𝑆 = 𝑚 𝑅 2 8  5 ¤ 𝜔 2 1 + 6 ¤ 𝜔 2 2 + ¤ 𝜔 2 3 + ( 2 𝜔 2 3 tan 𝜗 − 12 𝜔 2 𝜔 3 ) ¤ 𝜔 1 + 8 𝜔 1 𝜔 3 ¤ 𝜔 2 + ( 4 𝜔 1 𝜔 2 − 2 𝜔 1 𝜔 3 tan 𝜗 ) ¤ 𝜔 3  + . . . , (80) where the dots ( . . . ) represent further terms that are indepen- dent of the pseudoaccelerations ( ¤ 𝜔 1 , ¤ 𝜔 2 , ¤ 𝜔 3 ) so they can be neglected. Appell’ s equations are formulated as 𝜕 𝑆 𝜕 ¤ 𝜔 𝑗 = Π 𝑗 , 𝑗 = 1 , . . . , 3 , (81) where Π 𝑗 is the pseudoforce corresponding to the pseudoac- celeration ¤ 𝜔 𝑗 . The pseudoforces can be calculated from the virtual power of the activ e forces: 𝛿 𝑃 = G · 𝛿 v G = 3  𝑗 = 1 Π 𝑗 𝛿 𝜔 𝑗 . (82) Here G = [ 0 0 − 𝑚 𝑔 ] T F 0 represents the gravitational force, the only activ e force in our model, while 𝛿 v G = [ 𝛿 ¤ 𝑥 G 𝛿 ¤ 𝑦 G 𝛿 ¤ 𝑧 G ] T F 0 = [ · · − 𝑅 sin 𝜗 𝛿𝜔 1 ] T F 0 represents the virtual velocity; cf. ( 72 ). These yield the pseudoforces: Π 1 = 𝑚 𝑔 𝑅 sin 𝜗 , Π 2 = 0 , Π 3 = 0 . (83) Based on the acceleration energy ( 80 ), Appell’ s formula ( 81 ) and the pseudoforces ( 83 ), the pseudoaccelerations can be expressed as ¤ 𝜔 1 = 6 5 𝜔 2 𝜔 3 − 1 5 𝜔 2 3 tan 𝜗 + 4 𝑔 5 𝑅 sin 𝜗 , ¤ 𝜔 2 = − 2 3 𝜔 1 𝜔 3 , ¤ 𝜔 3 = − 2 𝜔 1 𝜔 2 + 𝜔 1 𝜔 3 tan 𝜗 , (84) These equations together with the generalized velocities ( 76 ) form the system of eight first order ordinary differential equations, which constitute the equations of the motion of the rolling wheel; cf. ( 6 ). A P P E N D I X B S TAB I L I T Y O F T H E R O L L I N G W H E E L T o prov e the stability of the steady states of the rolling wheel, we follow [ 31 ]. First, time is eliminated from the second and third equations of the essential motion ( 6 ) by expressing 𝜔 ′ 2 = − 2 𝜔 3 / 3 , 𝜔 ′ 3 = − 2 𝜔 2 + 𝜔 3 tan 𝜗 , (85) where 𝜔 ′ 2 = d 𝜔 2 d 𝜗 , 𝜔 ′ 3 = d 𝜔 3 d 𝜗 and ¤ 𝜗 = 𝜔 1 . Since the wheel is a conservati ve mechanical system, the total energy 𝐸 = 𝑇 + 𝑈 is constant, 𝐸 ( 𝑡 ) ≡ 𝐸 ∗ , where 𝑇 = 𝑚 𝑅 2  5 𝜔 2 1 + 6 𝜔 2 2 + 𝜔 2 3  / 8 and 𝑈 = 𝑚 𝑔 𝑅 cos 𝜗 are the kinetic and potential energies of the wheel, respectiv ely . Rearranging the total energy gi ves 𝜔 2 1 ≡ ¤ 𝜗 2 = 1 5 𝑚 𝑅 2  8 𝐸 ∗ − 𝑚 𝑅 ( 6 𝜔 2 2 𝑅 − 𝜔 2 3 𝑅 − 8 𝑔 cos 𝜗 )  | {z } = 𝑓 ( 𝜗 ; 𝜗 ic , 𝜔 2 , ic , 𝜔 3 , ic , 𝐸 ∗ ) . (86) The condition for stability is that 𝑓 has a strict maximum at 𝜗 = 𝜗 ∗ . Expressing d 𝑓 / d 𝜗 = 0 giv es back formula ( 12 ) representing the steady state surface in Figure 3 . The necessary and sufficient condition for steady states to be stable can be obtained with d 2 𝑓 / d 𝜗 2 < 0 , where ( 85 ) is substituted. This leads to ( 20 ) which is obtained from analyzing the linearized systems. A P P E N D I X C M O D E L D E R I V A T I O N F O R T H E U N I C Y C L E V I A T H E A P P E L L I A N A P P RO AC H W ith the added point mass the definition ( 33 ) of pseudove- locities is appropriate because the generalized velocities ( ¤ 𝑥 G , ¤ 𝑦 G , ¤ 𝜓 , ¤ 𝜗 , ¤ 𝜑 , ¤ 𝑟 ) can unambiguously be expressed by means of the pseudovelocities ( 𝜔 1 , 𝜔 2 , 𝜔 3 , 𝜎 ) and the generalized coordinates ( 𝑥 G , 𝑦 G , 𝜓 , 𝜗 , 𝜑 , 𝑟 ) when ( 76 ) is extended with ¤ 𝑟 = 𝜎 + 𝜔 1 𝑅 . The acceleration ener gy 𝑆 of the unicycle consisting of the rigid wheel and the point mass is defined as: 𝑆 = 1 2 𝑚 a 2 G + 1 2 𝜶 · J G 𝜶 + 𝜶 · ( 𝝎 × J G 𝝎 ) + 1 2 𝑚 0 a 2 A , (87) 15 where the first three terms are identical to those in ( 78 ) and ( 79 ). The velocity v A of the point mass 𝑚 0 can be expressed as v A =       ¤ 𝜓 ( 𝑅 sin 𝜗 − 𝑟 cos 𝜗 ) + ¤ 𝜑 𝑅 ¤ 𝑟 − ¤ 𝜗 𝑅 ¤ 𝜗𝑟       F 2 ≡       𝜔 2 𝑅 − 𝜔 3 𝑟 𝜎 𝜔 1 𝑟       F 2 , (88) while the acceleration a A becomes a A =       ¤ 𝜔 2 𝑅 − ¤ 𝜔 3 𝑟 + 𝜔 1 𝜔 3 ( 𝑟 tan 𝜗 − 𝑅 ) − 2 𝜔 3 𝜎 ¤ 𝜎 − 𝜔 2 1 𝑟 − 𝜔 2 3 𝑟 + 𝜔 2 𝜔 3 𝑅 ¤ 𝜔 1 𝑟 + 𝜔 2 1 𝑅 + 2 𝜔 1 𝜎 + 𝜔 2 3 𝑟 tan 𝜗 − 𝜔 2 𝜔 3 𝑅 tan 𝜗       F 2 . (89) These lead to 𝑆 = + 1 8  5 𝑚 𝑅 2 + 4 𝑚 0 𝑟 2  ¤ 𝜔 2 1 + 1 4  3 𝑚 + 2 𝑚 0  𝑅 2 ¤ 𝜔 2 2 + 1 8  𝑚 𝑅 2 + 4 𝑚 0 𝑟 2  ¤ 𝜔 2 3 + 𝑚 0 2 ¤ 𝜎 2 − 𝑚 0 𝑅𝑟 ¤ 𝜔 2 ¤ 𝜔 3 + 1 4  4 𝜔 2 1 𝑚 0 𝑅𝑟 + 8 𝜔 1 𝜎 𝑚 0 𝑟 + 𝜔 2 3 ( 𝑚 𝑅 2 + 4 𝑚 0 𝑟 2 ) tan 𝜗 − 𝜔 2 𝜔 3 ( 6 𝑚 𝑅 2 + 4 𝑚 0 𝑅𝑟 tan 𝜗 )  ¤ 𝜔 1 +  𝜔 1 𝜔 3 ( 𝑚 𝑅 2 − 𝑚 0 𝑅 2 + 𝑚 0 𝑅𝑟 tan 𝜗 ) − 2 𝜔 3 𝜎 𝑚 0 𝑅  ¤ 𝜔 2 + 1 4  2 𝜔 1 𝜔 2 𝑚 𝑅 2 + 𝜔 1 𝜔 3 ( 4 𝑚 0 𝑅𝑟 − 𝑚 𝑅 2 tan 𝜗 − 4 𝑚 0 𝑟 2 tan 𝜗 ) + 8 𝜔 3 𝜎 𝑚 0 𝑟  ¤ 𝜔 3 +  𝜔 2 𝜔 3 𝑚 0 𝑅 − 𝜔 2 1 𝑚 0 𝑟 − 𝜔 2 3 𝑚 0 𝑟  ¤ 𝜎 + . . . . (90) The pseudoforces Π 𝑗 are expressed based on the virtual po wer: 𝛿 𝑃 = G · 𝛿 v G + G A · 𝛿 v A + F · 𝛿 v A − F · 𝛿 v G = 4  𝑗 = 1 Π 𝑗 𝛿 𝜔 𝑗 , (91) where G = [ 0 0 − 𝑚 𝑔 ] T F 0 and G A = [ 0 0 − 𝑚 0 𝑔 ] T F 0 are the gravitational forces acting on the wheel and the point mass, re- spectiv ely , 𝛿 v A = [ 𝑅 𝛿𝜔 2 − 𝑟 𝛿 𝜔 3 𝛿 𝜎 𝑟 𝛿 𝜔 1 ] T F 2 , F is the con- trol force ( 31 ) and we used the notation 𝜔 4 = 𝜎 . This yields the following pseudoforces: Π 1 = ( 𝑚 + 𝑚 0 ) 𝑔 𝑅 sin 𝜗 − 𝑚 0 𝑔𝑟 cos 𝜗 , Π 2 = 0 , Π 3 = 0 , Π 4 = 𝑢 − 𝑚 0 𝑔 sin 𝜗 . (92) In case of 𝑚 0 > 0 , the pseudoaccelerations ( ¤ 𝜔 1 , ¤ 𝜔 2 , ¤ 𝜔 3 , ¤ 𝜎 ) can be obtained from the Appellian equations of the form ( 81 ), while the generalized velocities ( ¤ 𝑥 G , ¤ 𝑦 G , ¤ 𝜓 , ¤ 𝜗 , ¤ 𝜑 , ¤ 𝑟 ) can also be e xpressed from the definition of the pseudovelocities ( 33 ) and the kinematic constraints ( 2 ) for 𝜗 ≠ ± 𝜋 / 2 . These result in the equations of motion of the unicycle ( 34 ). A P P E N D I X D E L E M E N T S O F M A T R I X A I N ( 45 ) Note that in the follo wing formulas the steady state tilt angle 𝜗 ∗ , ya w rate ¤ 𝜓 ∗ , pitch rate ¤ 𝜑 ∗ and point mass position 𝑟 ∗ are not independent of each other but they must satisfy the relations in ( 38 ): 𝐴 12 = ¤ 𝜓 ∗  6 𝑚 𝑅 2 cos 𝜗 ∗ + 4 𝑚 0 𝑅𝑟 ∗ sin 𝜗 ∗  5 𝑚 𝑅 2 + 4 𝑚 0 𝑟 2 ∗ 𝐴 13 = 1 5 𝑚 𝑅 2 + 4 𝑚 0 𝑟 2 ∗  ¤ 𝜑 ∗ 𝑅 ( 6 𝑚 𝑅 + 4 𝑚 0 𝑟 ∗ tan 𝜗 ∗ ) + ¤ 𝜓 ∗ ( 4 𝑚 0 𝑅𝑟 ∗ tan 𝜗 ∗ + ( 4 𝑚 𝑅 2 − 8 𝑚 0 𝑟 2 ∗ ) ) sin 𝜗 ∗  𝐴 14 = 1  5 𝑚 𝑅 2 + 4 𝑚 0 𝑟 2 ∗  cos 𝜗 ∗  4 ¤ 𝜓 ∗ ¤ 𝜑 ∗ 𝑚 0 𝑅𝑟 ∗ + ¤ 𝜓 2 ∗ ( 4 𝑚 0 𝑅𝑟 ∗ tan 𝜗 ∗ − ( 𝑚 𝑅 2 + 4 𝑚 0 𝑟 2 ∗ ) ) cos 𝜗 ∗ + 4 𝑚 0 𝑔𝑟 ∗ sin 𝜗 ∗ cos 𝜗 ∗ + 4 𝑚 𝑔 𝑅 cos 2 𝜗 ∗  𝐴 16 = 𝑚 0  5 𝑚 𝑅 2 + 4 𝑚 0 𝑟 2 ∗  2  ( 16 𝑚 0 𝑔𝑟 2 ∗ − 20 𝑚 𝑔 𝑅 2 ) cos 𝜗 ∗ + ¤ 𝜓 ∗ ¤ 𝜑 ∗ ( ( 20 𝑚 𝑅 3 − 16 𝑚 0 𝑅𝑟 2 ∗ ) sin 𝜗 ∗ − 48 𝑚 𝑅 2 𝑟 ∗ cos 𝜗 ∗ ) − 32 𝑚 𝑔 𝑅𝑟 ∗ sin 𝜗 ∗ + 4 ¤ 𝜓 2 ∗ ( ( 5 𝑚 𝑅 3 − 4 𝑚 0 𝑅𝑟 2 ∗ ) sin 2 𝜗 ∗ − 20 𝑚 𝑅 2 𝑟 ∗ sin 𝜗 ∗ cos 𝜗 ∗ )  𝐴 21 = 1 3 𝑚 𝑅 2 + 2 𝑚 0 𝑅 2 + 12 𝑚 0 𝑟 2 ∗  − 4 𝑅 ¤ 𝜑 ∗ 𝑚 0 𝑟 ∗ + 2 ¤ 𝜓 ∗ ( ( 𝑚 0 𝑅 2 − 𝑚 𝑅 2 − 4 𝑚 0 𝑟 2 ∗ ) cos 𝜗 ∗ − 2 𝑚 0 𝑅𝑟 ∗ sin 𝜗 ∗ )  𝐴 25 = 4 𝑅 ¤ 𝜓 ∗ 𝑚 0 cos 𝜗 ∗ 3 𝑚 𝑅 2 + 2 𝑚 0 𝑅 2 + 12 𝑚 0 𝑟 2 ∗ 𝐴 31 = 1 3 𝑚 𝑅 2 + 2 𝑚 0 𝑅 2 + 12 𝑚 0 𝑟 2 ∗  − ¤ 𝜑 ∗ ( 6 𝑚 𝑅 2 + 4 𝑚 0 𝑅 2 ) + ¤ 𝜓 ∗ ( ( 12 𝑚 0 𝑟 2 ∗ − 3 𝑚 𝑅 2 − 2 𝑚 0 𝑅 2 ) sin 𝜗 ∗ − 20 𝑚 0 𝑅𝑟 ∗ cos 𝜗 ∗ )  𝐴 35 = − 24 ¤ 𝜓 ∗ 𝑚 0 𝑟 ∗ cos 𝜗 ∗ 3 𝑚 𝑅 2 + 2 𝑚 0 𝑅 2 + 12 𝑚 0 𝑟 2 ∗ 𝐴 52 = − ¤ 𝜓 ∗ 𝑅 cos 𝜗 ∗ , 𝐴 53 = ¤ 𝜓 ∗ ( 2 𝑟 ∗ cos 𝜗 ∗ − 𝑅 sin 𝜗 ∗ ) − ¤ 𝜑 ∗ 𝑅 𝐴 54 = − 𝑔 cos 𝜗 ∗ , 𝐴 56 = ¤ 𝜓 2 ∗ cos 2 𝜗 ∗ , 𝐴 61 = 𝑅 𝐴 73 = 1 cos 𝜗 ∗ , 𝐴 74 = ¤ 𝜓 ∗ tan 𝜗 ∗ , 𝐴 83 = − tan 𝜗 ∗ , 𝐴 84 = − ¤ 𝜓 ∗ cos 𝜗 ∗ , 𝐴 91 = 𝑅 sin 𝜓 ∗ cos 𝜗 ∗ , 𝐴 92 = 𝑅 cos 𝜓 ∗ , 𝐴 97 = − 𝑅  ¤ 𝜓 ∗ sin 𝜗 ∗ + ¤ 𝜑 ∗  sin 𝜓 ∗ , 𝐴 101 = − 𝑅 cos 𝜓 ∗ cos 𝜗 ∗ , 𝐴 102 = 𝑅 sin 𝜓 ∗ , 𝐴 107 = 𝑅  ¤ 𝜓 ∗ sin 𝜗 ∗ + ¤ 𝜑 ∗  cos 𝜓 ∗ . R E F E R E N C E S [1] M. Dozza, T . Li, L. Billstein, C. Svernl ¨ ov , and A. Rasch. How do different micro-mobility vehicles affect longitudinal control? results from a field experiment. J ournal of Safety Research , 84:24–32, 2023. [2] C. Ozaka, H. Kano, and M. Masubuchi. Stability of a monocycle- type in verted pendulum. In Thir d V ehicle Automation Symp. of J apan Automation Contr ol Society , pages 63–66, 1980. (in Japanese). [3] A. Schoonwinkel. Design and T est of a Computer Stabilized Unicycle . PhD thesis, Stanford University , 1987. [4] D. W . V os and A. H. von Floto w. Dynamics and nonlinear adaptiv e control of an autonomous unicycle: theory and experiment. In 29th IEEE Confer ence on Decision and Control , volume 1, pages 182–187, 1990. [5] Y . Naveh, P . Z. Bar-Y oseph, and Y . Hale vi. Nonlinear modeling and control of a unicycle. Dynamics and Control , 9:279–296, 1999. [6] D. Zenkov , A. Bloch, and J. Marsden. The Lyapunov-Malkin theorem and stabilization of the unicycle with rider . Systems & Control Letters , 45:293–300, 2002. [7] Y . Isomi and S. Majima. T racking control method for an underactuated unicycle robot using an equilibrium state. In IEEE International Confer ence on Contr ol and Automation , pages 1844–1849, 2009. [8] X. Ruan, J. Hu, and Q. W ang. Modeling with Euler-Lagrang equation and cybernetical analysis for a unic ycle robot. In 2nd International Confer ence on Intelligent Computation T echnology and Automation , volume 2, pages 108–111. IEEE, 2009. [9] S. I. Han and J. M. Lee. Balancing and velocity control of a unicycle robot based on the dynamic model. IEEE T ransactions on Industrial Electr onics , 62(1):405–413, 2014. 16 [10] L. Zhao, X. Zhang, Q. Xu, and J. Ji. Dynamics modeling and postural stability control of a unicycle robot. In International Conference on Fluid P ower and Mechatr onics (FPM) , pages 1123–1127. IEEE, 2015. [11] X. Cao, D. C. Bui, D. T ak ´ acs, and G. Orosz. Autonomous unicycle: mod- eling, dynamics, and control. Multibody System Dynamics , 61(1):43–76, 2023. [12] H. B. Brown and Y . Xu. A single-wheel, gyroscopically stabilized robot. In Pr oceedings of IEEE International Confer ence on Robotics and Automation , volume 4, pages 3658–3663. IEEE, 1996. [13] M. Q. Dao and K. Z. Liu. Gain-scheduled stabilization control of a unicycle robot. JSME International Journal, Series C, Mechanical Systems, Machine Elements and Manufacturing , 48(4):649–656, 2005. [14] K. Ghaff ari T . and J. K ¨ ovecses. Improving stability and performance of digitally controlled systems: The concept of modified holds. In IEEE International Confer ence on Robotics and Automation , pages 5173– 5180. IEEE, 2010. [15] H. Jin, T . W ang, F . Y u, Y . Zhu, J. Zhao, and J. Lee. Unicycle robot sta- bilized by the effect of gyroscopic precession and its control realization based on centrifugal force compensation. IEEE/ASME Tr ansactions on Mechatr onics , 21(6):2737–2745, 2016. [16] Z. Sheng and K. Y amafuji. Study on the stability and motion control of a unicycle: Part I: Dynamics of a human riding a unicycle and its modeling by link mechanisms. JSME International J ournal, Series C, Mechanical Systems, Machine Elements and Manufacturing , 38(2):249– 259, 1995. [17] H. Suzuki, S. Moromugi, and T . Okura. Dev elopment of robotic unicycles. Journal of Robotics and Mec hatr onics , 26(5):540–549, 2014. [18] J. P . Meijaard, J. M. Papadopoulos, A. Ruina, and A. L. Schwab. Linearized dynamics equations for the balance and steer of a bicy- cle: a benchmark and review . Pr oceedings of the Royal Society A , 463(2084):1955–1982, 2007. [19] H. Danko wicz. Multibody Mechanics and V isualization . Springer , 2005. [20] J. Xiong, N. W ang, and C. Liu. Stability analysis for the Whipple bicycle dynamics. Multibody System Dynamics , 48:311–335, 2020. [21] W . B. Qin, Y . Zhang, D. T ak ´ acs, G. St ´ ep ´ an, and G. Orosz. Nonholonomic dynamics and control of road vehicles: moving toward automation. Nonlinear Dynamics , 110(3):1959–2004, 2022. [22] E. J. Routh. The Advanced P art of A T reatise on the Dynamics of a System of Rigid Bodies . MacMillan, 1884. [23] A. V oss. Ueber die Dif ferentialgleichungen der Mechanik (About the differential equations of mechanics). Mathematische Annalen (Mathe- matical Annals) , 25:258–286, 1885. [24] A. M. Bloch, P . S. Krishnaprasad, J. E. Marsden, and R. M. Murray . Nonholonomic mechanical systems with symmetry . Ar chive for Rational Mechanics and Analysis , 136:21–99, 1996. [25] P . Appell. Sur une forme g ´ en ´ erale des ´ equations de la dynamique (On a general form of the equations of dynamics). Journal f ¨ ur die r eine und angewandte Mathematik (Journal for Pur e and Applied Mathematics) , 121:310–319, 1900. [26] J. W . Gibbs. On the fundamental formulae of dynamics. American Journal of Mathematics , 2(1):49–64, 1879. [27] P . V . V oronets. Ob uravnenih dvieni dl neg olonom- nyh sistem (On the equations of motion of nonholonomic systems). Matematiqeski ˘ ı Sbornik (Mathematical Collection) , 22(4):659– 686, 1901. [28] G. Hamel. Nichtholonome Systeme h ¨ oherer Art (Nonholonomic systems of a higher kind). Sitzungsberic hte der Berliner Mathematischen Gesellschaft (Meeting Reports of the Berlin Mathematical Society) , 37:41–52, 1938. [29] T . R. Kane. Dynamics of nonholonomic systems. ASME Journal on Applied Mechanics , 28:574–578, 1961. [30] F . Gantmacher . Lectur es in Analytical Mechanics . MIR Publishers, Moscow , 1970. [31] Ju. I. Neimark and N. A. Fufaev . Dynamics of Nonholonomic Systems , volume 33 of T ranslations of Mathematical Monographs . American Mathematical Society , 1972. [32] W . S. K oon and J. E. Marsden. The Hamiltonian and Lagrangian approaches to the dynamics of nonholonomic systems. Reports on Mathematical Physics , 40(1):21–62, 1997. [33] S. Ostrovskaya and J. Angeles. Nonholonomic systems revisited within the framework of analytical mechanics. Applied Mechanics Reviews , 57(7):415–433, 1998. [34] A. M. Bloch. Nonholonomic Mechanics and Contr ol . Springer , 2003. M ´ at ´ e B. V izi received the BSc degree in Mecha- tronic Engineering and the MSc degree in Me- chanical Engineering Modelling from the Budapest Univ ersity of T echnology and Economics, Hungary , in 2016 and 2018, respectively . He received the PhD degree in the same institution in 2024. Cur- rently he has a postdoctoral position at the University of Michigan. His research interests include nonlinear dynamics, control and time delay systems. G ´ abor Orosz received the MSc degree in Engi- neering Physics from the Budapest Univ ersity of T echnology , Hungary , in 2002 and the PhD degree in Engineering Mathematics from the University of Bristol, UK, in 2006. He held postdoctoral positions at the Uni versity of Ex eter , UK, and at the University of California, Santa Barbara. In 2010, he joined the Univ ersity of Michigan, Ann Arbor where he is currently a Professor in Mechanical Engineering and in Ci vil and En vironmental Engineering. From 2017 to 2018 he was a V isiting Professor in Control and Dynamical Systems at the California Institute of T echnology . In 2022 he was a Distinguished Guest Researcher in Applied Mechanics at the Budapest Univ ersity of T echnology and from 2023 to 2024 he was a Fulbright Scholar at the same institution. His research interests include nonlinear dynamics and control, time delay systems, machine learning, and data-driven systems with applications to connected and automated vehicles, traf fic flow , and biological networks. D ´ enes T ak ´ acs receiv ed his MSc and PhD in Me- chanical Engineering from the Budapest University of T echnology and Economics in 2005 and 2011, respectiv ely . Between 2011 and 2018, he worked in the MT A-BME Research Group on Dynamics of Machines and V ehicles in Budapest, Hungary . Since 2018, he has been an Associate Professor at Budapest Uni versity of T echnology and Economics, Budapest, Hungary . His research interests include tire and vehicle dynamics, nonlinear dynamics and time delay systems. G ´ abor St ´ ep ´ an recei ved the MSc and PhD degrees in mechanical engineering from Budapest University of T echnolgy and Economics, Hungary , in 1978 and 1982, respectively , and the DSc degree from the Hungarian Academy of Sciences, Budapest, Hun- gary , in 1994. He was a V isiting Researcher in the Mechanical Engineering Department of the Univer - sity of Newcastle upon T yne, UK, during 1988– 1989, the Laboratory of Applied Mathematics and Physics of the T echnical University of Denmark in 1991, and the Faculty of Mechanical Engineering of the Delft Univ ersity of T echnology during 1992–1993. He w as a Fulbright V isiting Professor at the Mechanical Engineering Department of the California Institute of T echnology during 1994–1995, and a V isiting Professor at the Department of Engineering Mathematics of Bristol Uni versity in 1996. He is currently a Professor of Applied Mechanics at the Budapest University of T echnology and Economics. He is a fellow of CIRP and SIAM, received the the Delay Systems Lifetime Achievements A ward of IF AC, the Caughey Dynamics A ward and the L yapunov A ward of ASME. He is a member of the Hungarian Academy of Sciences and the Academy of Europe. His research interests include nonlinear vibrations in delayed dynamical systems, and applications in mechanical engineering and biomechanics such as wheel dynamics (rolling, braking, shimmy), robotic force control, machine tool vibrations, human balancing, and traffic dynamics.