Symmetry and relative equilibria of a bicycle system moving on a surface of revolution

Finding out the relative equilibria and analyzing their stability is of great significance for revealing the intrinsic characteristics of mechanical system and developing effective controllers to improve system performance. In this paper, we study the symmetry and relative equilibria of a bicycle system moving on a surface of revolution. We notice that the symmetry group existing in the bicycle configuration description is a three-dimensional Abelian Lie group, and the rolling condition of the two wheels produces four time-invariant first-order linear constraints on the bicycle system. Therefore, the bicycle dynamics can be classified as a general Voronets system in which the Lagrangian and constraint distribution remain unchanged under the action of the symmetry group. Applying the Voronets equations to bicycle dynamics modeling, we obtain a seven-dimensional reduced dynamic system on the reduced constraint space. Theoretical analysis for the reduced dynamic system shows that it has the properties of time-reversal and lateral symmetries. In addition, two types of relative equilibrium points, the static equilibria and the dynamic equilibria, exist. Further theoretical analysis shows that the two kinds of relative equilibria both form one-parameter solution families, and the Jacobian matrix at an equilibrium point has some specific properties that support the relevant stability analysis. The necessary condition responsible for a stable static equilibrium point is that all the eigenvalues of the Jacobian matrix at the equilibrium point must lie on the imaginary axis of the complex plane. Due to the existence of zero eigenvalues of the Jacobian matrix, the stability of the dynamic equilibria is studied by limiting the reduced dynamic system to an invariant manifold based on the conservation of system energy. We prove in a strict mathematical sense that the dynamic equilibria may be Lyapunov stable, but cannot be asymptotically stable. Finally, symbolic computation combined with Whipple bicycle benchmark parameters was used for numerical simulations. We then use our numerical simulation to study how the parameter of the surface of revolution affects the relative equilibrium solution and its stability.

Most of the studies focused on a Whipple bicycle model proposed by Whipple [2] and Carvallo [1], who reduced the bicycle to a rigid four-body system with knife-edge wheels. When the bicycle moves on a ground, the nonslipping contact between its two wheels and the ground induces two normal geometric constraints and four tangent velocity constraints. All the constraints are time-invariant, and the velocity constraints adopt first-order linear forms. This leads to a nonholonomic Voronets system confined to a sevendimensional configuration space Q [23 -25]. The set of all possible velocities satisfying the four tangent velocity constraints is a three-dimensional subspace of the tangent space to Q [26,27]. The collection of these subspaces forms a constraint distribution D on Q, corresponding to a ten-dimensional submanifold of the tangent bundle T Q.
Compared to conventional nonholonomic systems such as a rattleback system [28,29] and a rolling disk [30], the bicycle's constraint equations assume more complicated forms due to the dependencies among individual contact points and the kinematic coupling between the saddle structure attitude and the handlebar steering [11]. Historically, many methods were developed to establish these constraint equations, especially for the case of horizontal ground. Psiaki [7] applied analytic geometry to derive by hand the geometric constraints of a bicycle moving on a horizontal plane. He stated for the first time that the algebraic equations of these geometric constraints are implicit and highly nonlinear. Basu-Mandal et al. [12] derived the constraint equations based on a condition for finding the lowest points of the two wheels. A similar method was used by Peterson and Hubbard [31,32] to obtain the holonomic constraint equations of a bicycle with toroidal tires. Wang et al. [33] parameterized the two wheel profiles and derived the constraint equations by an extremum condition that the wheel-road contact should be located at the lowest point on the wheel edge. Xiong et al. [24,25] employed an ordered process presented in [30] to obtain the constraint equations of a bicycle moving on a horizontal plane and an arbitrary surface, respectively.
Once the constraint equations are specified, the bicycle dynamical equations can be established by many different methods, such as the Newton-Euler approach in classical mechanics [10][11][12], and the methods of analytical mechanics including the Euler-Lagrange equations [6,12,25,33], the Kane's method [32] and the Gibbs-Appell equations [24]. In addition, geometric mechanics theory can also be used to analyze bicycle dynamics. Boyer et al. [17,34] used a general reduction-based approach to derive the reduced dynamic equations for the Whipple bicycle moving on the level ground. In this method, the bicycle's configuration space is modeled as a principal fiber bundle SE(3) × (S 1 ) 3 , through which the geometric nonlinearities caused by shifting the roll and pitch coordinates to the base space can be eliminated. The advantage of this method is that it can clearly expose some of the mathematical characteristics of bicycle dynamics. The disadvantage is that the kinematic constraints are not invariant under the left action of the structural group SE(3), so some reconstruction equations are needed in the reduced system. It is worth noting that different methods will produce different forms of motion equations, such as the ordinary differential equations (ODEs) [11,24] and the differential algebraic equations (DAEs) with different indices [12,25]. Therefore, applying these different forms of motion equations to solve the relative equilibria and analyzing their stability requires different mathematical methods and skills. For the nonlinear ODEs systems, it is well known that the motion stability can be analyzed by linearization around an equilibrium point [35]. For DAEs systems, Pappalardo et al. [36] proposed a direct linearization analysis method based on index-three form of DAEs to analyze the stability characteristics of multibody systems with holonomic and nonholonomic constraints. García-Agúndez et al. [37] obtained the linearized equations of any multibody system through the following three general numerical steps: (a) the linearization of the index-one form of the DAEs; (b) the reduction of the linearized equations from DAEs to ODEs form; and (c) a further reduction of the linearized equations by eliminating the holonomic constraints. Meijaard et al. [11] used the stability theory of linear systems, based on a linearized model with two degrees of freedom obtained by linear approximation to the constraint equations, to study the stability of the Whipple bicycle in an upright straight motion. Based on the linearization of the bicycle's nonlinear ODEs model obtained from the Gibbs-Appell equations, Xiong et al. [24] studied the stability of the bicycle's straight and circular motions on the flat ground. Basu-Mandal et al. [12] developed a DAEs model for the Whipple bicycle moving on the flat ground, then numerically obtained the solutions of the bicycle's circular motion and used a method similar to step (a) in [37] to study the stability around the circular motion. Xiong et al. [25] applied the Lagrangian equations of the first type to establish the nonlinear DAEs for the bicycle moving on a surface of revolution and studied the stability of the bicycle's circular motion by using the similar method in [37], with consideration for the effect of cyclic coordinates.
In this paper, we re-examine the dynamics of the Whipple bicycle moving on a surface of revolution. The bicycle in this case is a system with a three-dimensional Abelian symmetry group H . The Lagrangian and constraint distribution of the bicycle system, under the action of H , remain invariant. This property allows us to reduce the dynamic equations of the bicycle system through the Voronets equation (VE), which is a set of unified equations in geometric mechanics for dealing with dynamic system subject to time-invariant first-order linear constraints [23,38]. Compared to the Lagrangian equations [25], the VE method is free of Lagrangian multipliers and produces a set of ODEs with the minimum dimension. Compared to the Gibbs-Appell equations [24], the set of ODEs obtained from VE method has clear geometric structures, such as Ehresmann connection and its curvature [27], to reveal the internal characteristics of the system dynamics. The purpose of this paper is to reveal the intrinsic structures and mathematical properties of the relative equilibria based on the bicycle ODEs system obtained from the VE method.
By using the VE method to deal with the bicycle system, we obtain a seven-dimensional reduced dynamic system on the quotient manifold of D under the action of H [26,27,39]. More importantly, from the reduced dynamic system we can prove that the bicycle system has some specific mathematical properties. For example, the bicycle system takes time-reversal and lateral symmetries as discussed by Meijaard et al. [11]; two kinds of relative equilibria, the static equilibria related to the bicycle's static equilibrium configuration and the dynamic equilibria representing the uniform circular motion of the bicycle, exist in the system; the two kinds of equilibria under certain regular conditions both form one-parameter solution families. Furthermore, we discuss the Lyapunov stability of these relative equilibria through analyzing the properties of the associated Jacobian matrix. The necessary condition responsible for a stable static equilibrium point is that all the eigenvalues of the Jacobian matrix at the equilibrium point must lie on the imaginary axis of the complex plane. We studied the stability of the dynamic equilibria by limiting the reduced dynamic system to an invariant manifold based on the conservation of system energy. We prove that all the relative equilibria are not hyperbolic, and thus conclude that a stable equilibrium of the bicycle system cannot be asymptotically stable [35]. Further analysis shows that the disturbed motion of the bicycle system around a stable relative equilibrium will converge to a nearby relative equilibrium. Finally, we perform numerical simulations via symbolic computation combined with Whipple bicycle benchmark parameters to show how the parameter of the surface of revolution affects the relative equilibrium solution and its stability.
The rest of this paper is organized as follows: In Sect. 2, we analyze the symmetry of the Whipple bicycle system on the surface of revolution and derive the reduced dynamic system using the Voronets equations. In Sect. 3, we study the relative equilibria of the reduced dynamic system. We analyze the properties of the Jacobian matrix at a relative equilibrium in Sect. 4 and then discuss its Lyapunov stability in Sect. 5. In Sect. 6, we apply symbolic computation to carry out numerical simulations for the Whipple bicycle in conjunction with the benchmark parameters. Conclusions are offered in Sect. 7.

Symmetry and reduced dynamic system of a bicycle
In this section, we first describe the configuration of the Whipple bicycle system moving on a surface of revolution. Based on the configuration description, we identify a three-dimensional Abelian Lie group under the action of which the Lagrangian and constraint distribution of the bicycle system remain invariant. Since the bicycle system is subject to four time-invariant firstorder linear constraints, we can use the Voronets equations to reduce the bicycle dynamics into a minimum ODEs system. The structure of the ODEs system can be clearly exposed based on the symmetries of the bicycle system. Similar to the case where the bicycle system moving on a horizontal plane [11], the dynamics of the bicycle moving a surface of revolution has the timereversal and lateral symmetries. Figure 1 adapted from [11,33] shows a Whipple bicycle that consists of four rigid bodies including a rear wheel, a saddle structure, a head structure and a front wheel. The two wheels are circularly symmetric and make ideal knife-edge rolling contact with the ground. The saddle and head structures have lateral symmetries in their shape and mass distributions. The effects of structural compliance, joint friction and rolling friction are assumed to be negligible. We choose the bicycle's upright configuration as the reference state. In the configuration, the bicycle is assumed to rest on the horizontal plane h (which is not necessarily be the ground) such that all of its four bodies are symmetric with respect to a common vertical plane v . Below we define five coordinate frames, one inertial and four body-fixed, each with its axes 1 and 3 inside v , and axis 2 orthogonal to it. The inertial frame F i = {O; i 1 , i 2 , i 3 } has its origin in h (which is thus spanned by its axes 1 and 2). The rear and saddle frames F r = {O r ; e r,1 , e r,2 , e r,3 } and F s = {O s ; e s,1 , e s,2 , e s,3 } are oriented the same as F i and located, respectively, at the center O r of the rear wheel and the center of mass O s of the saddle structure. The head frame F h = {O h ; e h,1 , e h,2 , e h,3 } is located at the center of mass O h of the head structure; the frame's axis 3 is aligned with the steering axis (i.e., the straight axis of the head tube) and points upward.

Configuration description
The front frame F f = {O h ; e f,1 , e f,2 , e f,3 } is located at the center O f of the front wheel and has the same orientation of the head frame F h . All five frames are illustrated in Fig. 1.
In this paper, the ground is assumed to be a surface of revolution. In the inertial frame F i , it is a submanifold of the three-dimensional Euclidean space R 3 : where ρ = (x 1 ) 2 + (x 2 ) 2 . Clearly, a horizontal plane is a special surface of revolution where f takes the form f (ρ) ≡ 0, ∀ρ ∈ [0, +∞).
A movement of the bicycle on is subject to two geometric constraints and four velocity constraints. The two types of constraints are both time-invariant, and the velocity constraints take the first-order linear forms. Therefore, the bicycle system is a nonholonomic Voronets system [23]. A detailed description of these constraints was given in our previous work [25]. As a result, the configuration space of the bicycle is a sevendimensional manifold: with the coordinates q = (θ, δ, φ r , (x, y, ψ), φ f ), where θ is the lean angle of the rear frame, δ is the steer angle, φ r is the rotation angle of the rear wheel, x, y are the first two coordinates of a reference point D (which is located at the intersection of the steering axis and the coordinate axis O f e f,1 of the front frame F f ) in the inertial frame F i , ψ is the yaw angle of the rear frame, and φ f is the rotation angle of the front wheel. At any configuration q ∈ Q, the set D q of all possible velocities satisfying the four first-order linear constraints is a three-dimensional subspace of the tangent space T q Q. The collection of D q forms a constraint distribution D on Q: Clearly, D is a ten-dimensional submanifold of the tangent bundle T Q.

Symmetry of Lagrangian and constraint distribution
We introduce the symmetry group for the bicycle system on as follows: where the three components of product H 1 , H 2 , H 3 represent the rotation of the rear wheel φ r , rotation of the rear frame in the horizontal plane ψ (i.e., yaw motion) and rotation of the front wheel φ f , respectively. Symmetry group H is a three-dimensional Abelian Lie group as will be stated below.
The left action of H on Q is defined as a smooth mapping : For each h ∈ H , Eq. (5) induces a diffeomorphism h from Q onto itself: and we denote by hq = h (q) for simplicity. To simplify the expression of the action, we use another set of coordinates of Q, i.e.,q = (θ, δ, φ r ,x,ỹ, ψ, φ f ), where x = x cos(ψ) + y sin(ψ), y = −x sin(ψ) + y cos(ψ).
Clearly,q consists of two parts: the coordinates (θ, δ,x, y) of the quotient manifold Q/H , and the coordinates (φ r , ψ, φ f ) of the Lie group H . Therefore, the left action of H on Q can be expressed in these coordinates as follows: Denote by L : T Q → R the Lagrangian of the bicycle system on . Using the set of coordinatesq, it can be written as follows: where K i j s are the components of the mass matrix, V is the gravitational potential energy, and Einstein's summation convention is used for the indices i, j = 1, · · · , 7 1 . In addition, denote byq = (σ ,s), wherẽ σ = (σ 1 ,σ 2 ,σ 3 ) = (θ, δ, φ r ) ands = (s 1 , · · · ,s 4 ) = (x,ỹ, ψ, φ f ). The four first-order linear constraint equations can then be written as follows: whereσ = (θ,δ,φ r ) are three independent velocities. Thus, we have which means that η = (q,σ ) forms a set of coordinates of D, whereq represents the base coordinates andσ represents the fiber coordinates. Note that the surface of revolution is rotational symmetric about the axis i 3 of F i and the two wheels of the bicycle are rotational symmetric about their spin axes. We claim that under the action of H , the bicycle system has the following two invariances [27]: According to Eq. (8), the Lagrangian and constraint distribution invariances mean that φ r , ψ, φ f do not appear explicitly in the expressions of Lagrangian L and velocity constraint equations (10), i.e., we have Therefore, φ r , ψ, φ f are cyclic coordinates, consistent with the fact that H is an Abelian Lie group [40]. Due to the constraint distribution invariance, we can define the left action of H on D as follows: ∀η = (q,σ ) ∈ D, ∀h ∈ H , we have The corresponding quotient manifold D/H , termed as reduced constraint space in geometric mechanics [26,27,39], is seven-dimensional, whose coordinates For the sake of a brief description, we defineξ = (ξ 1 , · · · ,ξ 4 ) = (x,ỹ, θ, δ) as the configuration variables included in v, and thus v = (ξ ,σ ).

Reduced dynamic system
We employ the Voronets equations [27] to obtain the motion equations of the uncontrolled Whipple bicycle on . Denote by L c = L|D : D → R the constrained Lagrangian. According to Eq. (12), L c takes the following form: where m αβ , α, β = 1, 2, 3 are functions of K i j and A a α . The above expression shows that L c is also invariant under the action defined in Eq. (13), which is a necessary consequence of the Lagrangian and constraint distribution invariances.
The Voronets equations take the following forms where with b = 1, · · · , 4, α, β = 1, 2, 3, are anti-symmetric about the two indices α and β. From the point of geometric mechanics view [27], A a α can be considered as the coefficients of an Ehresmann connection, and constraint distribution D is given by the horizontal subbundle associated with this connection. Then, B b αβ represent the coefficients of the curvature of this connection. By combining similar terms, Eq. (15) reads where c αβγ and P α are functions ofξ , i.e., c αβγ = c αβγ (ξ ), α, β, γ = 1, 2, 3, Noting the Lagrangian L and the constrained Lagrangian L c in Eq. (9) and Eq. (14), respectively, and applying the Voronets equations in Eq. (15), some specific coefficients in Eq. (18), which are important for revealing the internal characteristics of the bicycle dynamics, can be expressed as follows: By transforming Eq. (17) into first-order forms, and combining with the first two velocity constraint equations in Eq. (10) (i.e., a = 1, 2), we obtain the reduced dynamic system on D/H as follows: where where M = [m αβ ] ∈ R 3×3 , and Now, let us discuss the symmetries of reduced dynamic system (20) about time t and the lateral direction of the bicycle. Defineṽ = (x,ỹ, θ, δ, −θ, −δ, −φ r ) to represent the system state as time t is reversal. Clearly, we have Meijaard et al. [11] claimed that the bicycle system moving on a horizontal plane takes lateral symmetry. This symmetry still remains in the bicycle moving on a surface of revolution. Letv = (x, −ỹ, −θ, −δ, −θ, −δ,φ r ). Then, if v(t) is a solution of Eq. (20), so isv(t). This symmetry is equivalent to the invariance of system (20) under the left action of a discrete group Z 2 on D/H and is attributed to the assumption that both the four bodies of the Whipple bicycle and the ground have reflection symmetries about the vertical plane v in the reference configuration (see Fig. 1). As a result, the tangent vector field Y has the following properties:

Relative equilibria
In geometric mechanics, we have the following definition for relative equilibrium [27].
Definition 1 A relative equilibrium is an equilibrium that is defined on the reduced dynamic system.
Based on this definition, we can find the relative equilibria of reduced dynamic system (20) by where there are five constants to be determined. From Eq. (21), we have the following five algebraic equations: Noting that reduced dynamic system (20) takes the time-reversal and lateral symmetries, we can claim that, if v 0 is a relative equilibrium, so areṽ 0 = We can classify the relative equilibria into two types, depending on whether ω 0 is zero or not. We give the following definition: Clearly, the static equilibrium corresponds to a scenario where the bicycle stands statically on surface . Under the condition ω 0 = 0, we can find that the first two equations in Eq. (25) are always satisfied, and the static equilibrium configurationξ s 0 of the bicycle can be solved through the last three equations, namely, It is clear that the number of the independent algebraic equations in Eq. (26) is less than the number of the elements inξ s 0 . Thus in general, the solution is not unique. Let us define J s as the Jacobian matrix of Eq. (26) around a static equilibrium v s 0 = (ξ s 0 , 0, 0, 0). Namely, Under a regular condition that J s is of full rank, we suppose that, without loss of generality, the submatrix consisting of the last three columns of J s is nonsingular. According to the implicit function theorem, in a neighborhood U (ξ s 0 ) ofξ s 0 we can use equations P α (ξ ) = 0, to determine single-valued functions ofỹ, θ, δ with respect tox, namely, as a dynamic equilibrium, whereξ c 0 corresponds to the dynamic equilibrium configuration. In principle, we can uniquely determine the five values of (ξ c 0 , ω 0 ) since there are five equations in Eq. (25). However, further analysis shows that the last equation in Eq. (25) is an identical equation. In fact, if ω 0 = 0, the first two equa- Thus the last equation in Eq. (25) (α = 3) is an identical equation, and Eq. (25) can then be reduced as follows: where there are five variables involved within four equations.
Let us define J c as the Jacobian matrix of Eq. (29) around v c 0 . Namely, Similar to the case of static equilibria, we can also conclude that, if J c is of full rank, there exists a neighborhood of v c 0 , in which all the dynamic equilibria form a one-parameter solution family. Therefore, we have the following theorem for the solution family of relative equilibria.
represents the relative equilibrium of reduced dynamic system (20). Whatever v 0 is at static or dynamic equilibrium, under a regular condition that the associated Jacobian matrix J s or J c has full rank, there always exists a neighborhood U (v 0 ) of v 0 , in which all the relative equilibria of the same kind form a one-parameter solution family that passes through v 0 .
We claim that the dynamic equilibrium v c 0 = (ξ c 0 , 0, 0, ω 0 ) is related to a uniform circular motion of the bicycle on surface . This point is illustrated as follows. According to constraint equations (10), we easily find that v c 0 should be associated with a constant precession velocityψ 0 = −A 3 3 (ξ c 0 )ω 0 . Furthermore, we can apply Eq. (7) to obtain the motion equations of x(t) and y(t) under the initial condition ψ(0) = 0: x(t) = ρ 0 cos ψ 0 t + 0 , where ρ 0 is the radius of the trajectory of point D, 0 is the phase difference, and these two constants are given by the following equations:

Mathematical properties of the Jacobian matrices
In this section, we discuss the essential properties of the Jacobian matrix J r (v 0 ) of Eq. (20) at a given relative equilibrium v 0 . It can be proved that J r (v 0 ) must have a zero eigenvalue. For a pair of relative equilib-ria (v 0 ,ṽ 0 ) satisfying the time-reversal symmetry, the corresponding eigenvalues of J r (v 0 ) and J r (ṽ 0 ) must be negative of each other. For a pair of relative equilibria (v 0 ,v 0 ) satisfying the lateral symmetry, the corresponding J r (v 0 ) and J r (v 0 ) must take the same eigenvalues. Now, we will prove each of these properties in turn. Generally, the Jacobian matrix of Eq. (20) at a relative equilibrium v 0 = (ξ 0 , 0, 0, ω 0 ) can be written as Noting the expressions of Y i listed in Eq. (21) and the structure of v 0 , we can present the concrete expression of J r (v 0 ) as follows: where Property 1 For each relative equilibrium v 0 , the Jacobian matrix J r (v 0 ) must have at least one zero eigenvalue.
Proof If J r (v 0 ) is rank-deficient, the property is naturally true. Based on Eq. (34), the rank of matrix J r (v 0 ) can be equally computed by Under a regular condition that M −1 0 is nonsingular, the following identity is true: [ J] 3×4 [L] 3×1 5×5 .
When v 0 is a static equilibrium, we have [ In this case, J r (v 0 ) has at least two zero eigenvalues.

Noting that [ J]
j 3 = ∂ P 3 ∂v j − ∂c 333 ∂v j ω 2 0 , we can find that the last row of J (v 0 ) is the linear combination of its first two rows. Therefore, using the elementary transformation of matrix, we have where J c takes the same form as shown in Eq. (30). This implies the maximum value of rank J (v 0 ) is 4, thus rank J r (v 0 ) ≤ 6. Clearly, whatever v 0 is at static or dynamic equilibrium, the Jacobian matrix J r (v 0 ) must be rankdeficient; thus, it should have at least one zero eigenvalue.
Proof Based on Eq. (23), if v 0 andṽ 0 correspond to a pair of relative equilibria satisfying the time-reversal symmetry, we can generally represent the Jacobian matrix J r (v 0 ) and J r (ṽ 0 ) as the following partitioned matrices, Let us denote by μ an eigenvalue of J r (v 0 ), which satisfies the following characteristic equation: After performing four row operations and three column operations with (−1) on the above equation, we get Noting the expression of J r (ṽ 0 ), we should have namely, −μ is an eigenvalue of J r (ṽ 0 ). Conversely, ifμ is an eigenvalue of J r (ṽ 0 ), we can prove that −μ is an eigenvalue of J r (v 0 ). Therefore, we can conclude that the eigenvalues of matrices J r (v 0 ) and J r (ṽ 0 ) must be negative of each other.
are two relative equilibria satisfying the lateral symmetry. This symmetry makes J r (v 0 ) and J r (v 0 ) with identical eigenvalues.
Proof According to Eq. (24), the entries of the Jacobian matrices around v 0 andv 0 satisfy the following relationships, j, k = 2, · · · , 6. This means that there exists a transform matrix Namely, J r (v 0 ) and J r (v 0 ) are similarity matrices. The property thus follows.

Lyapunov stability of relative equilibria
In this section, we will analyze the Lyapunov stability of the relative equilibria of reduced dynamic system (20) in a strict sense of mathematics. We first discuss the stability of static equilibria based on the timereversal symmetry of Eq. (20). For the stability analysis of a dynamic equilibrium, the zero eigenvalue of matrix J r (v 0 ) causes complexity. For this reason and noting that the reduced constrained energy is conservative, we perform the stability analysis based on an invariant manifold of Eq. (20), which corresponds to a level set of the reduced constrained energy.

Stability of static equilibria
Because bicycle dynamics is similar to an inverted pendulum system, intuitively any static equilibrium of the bicycle cannot be stable on the surface. Here, based on the second property of the Jacobian matrix discovered in the previous section, we have the following proposition:

Proposition 1 The necessary condition responsible for a stable static equilibrium is that all the eigenvalues of its Jacobian matrix must be located at the imaginary axis of the complex plane.
Proof Suppose that v 0 is a static equilibrium point of the bicycle system. It is clear that the corresponding static equilibria satisfying the time-reversal symmetry should be given by v 0 =ṽ 0 . Accordingly, we have J r (v 0 ) = J r (ṽ 0 ); thus, both of them should take the same eigenvalues. In addition, the second mathematical property of the Jacobian matrix indicates that time-reversal symmetry makes the eigenvalues of the matrices J r (v 0 ) and J r (ṽ 0 ) negative of each other. This implies that, if J r (v 0 ) takes an eigenvalue with a negative real part, there is another eigenvalue that must take a positive real part. Therefore, v 0 must be an unstable static equilibrium if J r (v 0 ) has at least one eigenvalue that is not located at the imaginary axis of the complex plane.

Stability of dynamic equilibria
In this section, we study the stability of a dynamic We know from our previous analysis that the Jacobian matrix J r (v 0 ) must have a zero eigenvalue. This implies that v 0 is not a hyperbolic equilibrium point, and its Lyapunov stability cannot be directly determined via the linearized system around it [35]. In this case, we can claim v 0 is unstable if its Jacobian matrix has at least one eigenvalue with positive real part. However, we cannot say anything about v 0 even though all nonzero eigenvalues of J r (v 0 ) have negative real parts. We notice that the bicycle system takes a specific property as follows. Under the rolling condition of the two wheels, the total energy of the bicycle system, which equals the sum of kinetic and potential energies, is conserved. Due to the Lagrangian and constraint distribution invariances, it induces a well-defined reduced constrained energy E : D/H → R on D/H : Clearly, E(v) is also conserved, and it thus can be considered as a first integral of Eq. (20). Equivalently, each of its level set, is an invariant set. Therefore, for a given dynamic equilibrium v 0 , we can always define an invariant set Q E 0 , where E 0 = E(v 0 ), and the stability of v 0 can then be analyzed by restricting the motion of system (20) on Q E 0 . We divide the coordinates of D/H into two parts: v  = (u,φ r ), where u = (x,ỹ, θ, δ,θ,δ), and divide the components of the tangent vector field Y into two parts: Y = (X, Y 7 ), where X = (Y 1 , · · · , Y 6 ). Let us define a function:Ẽ(u,φ r , Clearly we haveẼ(u 0 , ω 0 , E 0 ) = 0, and ∂Ẽ(u 0 , ω 0 , E 0 )/∂φ r = m 33 (ξ 0 )ω 0 = 0, where u 0 = (ξ 0 , 0, 0). According to the implicit function theorem, the equa-tionẼ(u,φ r , E 1 ) = 0 determines a single-valued function ofφ r with respect to (u, E 1 ) in a neighborhood of (u 0 , ω 0 , where U (u 0 ) and U (E 0 ) are neighborhoods of u 0 and E 0 , respectively, and ω 0 = F(u 0 , E 0 ). Therefore, for each E 1 ∈ U (E 0 ), the invariant set Q E 1 can be locally expressed as follows: which means that in local, Q E 1 is an invariant manifold of Eq. (20) with the coordinates u.
According to the theory of nonlinear dynamics [35,41], if J r (v 0 ) does not have a positive real part eigenvalue, then the stability of v 0 is unclear. Next, our stability analysis focuses only on the case where J r (v 0 ) has six negative real part eigenvalues.
On the other hand, if J r (v 0 ) has six nonzero eigenvalues, we must have rank J r (v 0 ) = 6. According to Eq. (35), J c has full rank, and thus according to Theorem 1, v 0 must belong to a one-parameter solution family of dynamic equilibria. We have the following corollary for the curves of the solution families. Proof According to the second property of the Jacobian matrix of a given equilibrium point v 0 , the distribution of eigenvalues of J r (v 0 ) must be symmetrical about the imaginary axis of the complex plane if ω 0 = 0. As a result, v 0 must be a dynamic equilibrium when J r (v 0 ) has six eigenvalues with negative real parts.
Denote byv(s) the one-parameter solution family of the relative equilibria that passes through v 0 , where s is an arc length parameter. Without loss of generality, we assume thatv(0) = v 0 . Clearly, we have Y (v(s)) ≡ 0. Differentiating the left-and right-hand sides of this equation with respect to s at s = 0, we have Therefore, the eigenvector of J r (v 0 ) in associated with its zero eigenvalue should be parallel to the tangent vector of the parameterized curvev(s) of the solution family at point v 0 .
Furthermore, we have the following corollary for a six-dimensional reduced system of system (20) when its motion is limited on the invariant manifold Q E 0 .
Corollary 2 Suppose that v 0 is a relative equilibrium of system (20), and all eigenvalues of its Jacobian matrix J r (v 0 ), except for the inherent zero eigenvalue, have negative real parts. Then, we can define a sixdimensional reduced system in coordinates u by limiting original system (20) on the invariant manifold Q E 0 , where E 0 = E(v 0 ). The reduced system has an isolated asymptotically stable equilibrium at u 0 .
Proof Under the restriction of Eq. (20) on Q E 0 , we establish a six-dimensional reduced system with coordinates u, which can be expressed as follows: Clearly, u 0 is an equilibrium of Eq. (41), at which the Jacobian matrix of Eq. (41) is given by Let us consider a perturbed system in coordinates u and with variable energy E 1 , restricted in a small region u ∈ U (u 0 ) and E 1 ∈ U (E 0 ). The neighborhood of v 0 in D/H can then be represented as a union of invariant manifolds Q E 1 . In addition, the invariance of Q E 1 implies that ∀(u, E 1 ) ∈ U (u 0 ) × U (E 0 ), the following identity holds: Differentiating both sides of Eq. (43) with respect to u and E 1 at point ζ 0 = (u 0 , E 0 ), respectively, we obtain the following equations: Here we have used the equation X(v 0 ) = 0. Noting that we always have ∂ F(ζ 0 )/∂ E 1 = 0 according to the differential of Eq. (38) with respect to variable E 1 , Eq. (44) can be simplified as It is worth noting that the Jacobian matrix J r (v 0 ) of original system (20) can be written by, Let us define an invertible matrix then, J r (v 0 ) undergoes by the similarity transformation of T as follows: Since J r (v 0 ) has six eigenvalues with negative real parts as well as one zero eigenvalue, we know that the all six eigenvalues of J (u 0 ) must have negative real parts. As a result, u 0 is an asymptotically stable equilibrium of reduced system (41), and thus it is isolated.
Based on the above analysis, we have the following theorem to illustrate the stability of a dynamic equilibrium v 0 .
Theorem 2 Suppose that v 0 is a relative equilibrium, and the other six eigenvalues of the Jacobian matrix J r (v 0 ) have negative real parts, then v 0 is a Lyapunov stable but not asymptotically stable equilibrium of system (20). Furthermore, the solution v(t) of Eq. (20) under the initial condition v(0) which departs from v 0 to a sufficiently small extent will converge to a nearby equilibrium v 0 =v(s 0 ) in the one-parameter solution familyv(s), where s 0 is uniquely determined by the equation E(v(s 0 )) = E(v(0)).
Proof In order to prove this theorem, let us first illustrate the relationship between the set of relative equilibriav(s) and the set of invariant manifolds Q E 1 in the neighborhood of v 0 .
On the one hand, according to corollary 1, the zero eigenvector of J r (v 0 ) is parallel to dv(0)/ds. On the other hand, we know from similarity transformation (46) that this zero eigenvector should be parallel to the following vector: Thus, we have (dv(0)/ds) = k 0β , where k 0 = 0. Note that the reduced constrained energy E 1 along the family curvev(s) is a function of s, i.e., E 1 = E(s) = E(v(s)). The derivative ofÊ(s) with respect to s at s = 0 is given by Differentiating Eq. (38) with respect to u at point ζ 0 , we have Substituting Eqs. (47) and (49) into Eq. (48), we get which means thatÊ(s) will take different values for different s, when s is limited in a neighborhood U (0) of s = 0. Accordingly, a neighbor-hoodÛ (E 0 ) ⊂ U (E 0 ) of E 0 exists, such that the restriction ofÊ on U (0) gives a diffeomorphism from U (0) ontoÛ (E 0 ). We denote byf :Û (E 0 ) → U (0) the inverse map ofÊ| U (0) . We can then claim that there is a neighborhood U (v 0 ) of v 0 , in which a one-to-one correspondence between the set of relative equilibriav(s) and the set of invariant manifolds Q E 1 exists under the mapping s → E 1 =Ê(s).
Note that the neighborhood of v 0 in D/H can be represented as a union of invariant manifolds Q E 1 . On the one hand, the restriction of Eq. (20) on Q E 0 has an asymptotically stable equilibrium at v 0 . On the other hand, by continuity the reduced system on each nearby invariant manifold Q E 1 has an asymptotically stable equilibrium atv(f (E 1 )), which is located in a small neighborhood of v 0 . Therefore, v 0 is stable but not asymptotically stable. In addition, the solution v(t) of Eq. (20) under the initial condition v(0) will stay on the invariant manifold Q E(v(0)) for all time. Thus, if v(0) is in a sufficiently small neighborhood of v 0 , we know that v(t) will converge tov(s 0 ), where s 0 =f (E(v(0))), as t tends to +∞.

Numerical results
In this section, we apply the theory developed in this paper to numerically investigate the relative equilibria of a Whipple bicycle with the benchmark parameters shown in [11,12]. Table 1 lists the values of these parameters, including the wheel base w, the trail c, the tilt angle λ, the positions of the center of mass of the two frames x k , z k (k = s, h), the radii of the two wheels R k (k = r, f ), the masses of the four rigid bodies m k (k = r, s, f, h) and the components of the inertia tensors of the four rigid bodies I k,x x , I k,yy (k = r, s, f, h) and I k,zz , I k,xz (k = s, h). The acceleration of gravity is g = 9.81m/s 2 .
We suppose that the bicycle is set on the surface of a paraboloid of revolution with the following function where κ > 0 is a coefficient, and ρ is the distance of a point on the surface to the rotation axis (see Eq. (1)).
In terms of the parameters involved in the bicycle and the paraboloid of revolution, we can apply symbolic computation to establish the expressions of the system Lagrangian function (see Eq. (9)), the nonholonomic constraint equations (see Eq. (10)) as well as the Voronets equations (see Eq. (17)). The coefficients involved in Eq. (25) can then be determined, and  Eqs. (26) and (29) are employed to find the solutions of the static and dynamic equilibria, respectively. Based on Theorem 1, we know that these relative equilibria are not isolated, but belong to a family of oneparameter solutions. Therefore, these relative equilibria can be expressed as functions with respect to an independent variable selected from the state variables of the relative equilibria. Due to the time-reversal and lateral symmetries of reduced dynamic system (20), we compute the solution families by limiting ω 0 ∈ [0, +∞) and δ 0 ∈ [0, π] without loss of generality.
For the static equilibria, we find two solution families. Among them, one takes the special expression v s (ỹ) = (x 0 (κ),ỹ, 0, 0, 0, 0, 0), wherex 0 (κ) is a function of κ, andỹ ∈ R is an independent variable. Clearly, the lean and steer angles in this family are zero. Whenỹ = 0, the static equilibrium v s (0) = (x 0 (κ), 0, 0, 0, 0, 0, 0) represents the configuration where the bicycle is positioned in a radial direction. Figure 3 shows the curve ofx 0 with respect to κ in the range of κ ∈ (0, 1]m −1 . We can see thatx 0 increases first and then decreases as κ increases. Another solution family corresponds to a curve passing through an accumulation point. In order to present a clear description for the equilibrium configuration of the bicycle, we use the variables ρ 0 and 0 defined in Eq. (32) to replace the variablesx 0 andỹ 0 . Clearly, (ρ 0 , 0 ) are polar coordinates of (x 0 ,ỹ 0 ). Figure 4 shows this kind of solution family in the case of κ = 1/200m −1 , where the curves of 0 , θ 0 , δ 0 are plotted with respect to the independent variable ρ 0 . In addition, the variable ρ 0 is scaled by arctan (ρ 0 ), such that its values are limited in a scope [0, π/2]. For each of these relative equilibria, we can see that θ 0 is close to zero, but steer angle δ 0 is always large. This means that the bicycle always takes an almost upright configuration with a large steer angle. At the same value of ρ 0 , multiple solutions for 0 , θ 0 , δ 0 are allowed in the solution families of the static equilibria. The accumulation point, induced by the intersection of the static and dynamic solution families, is designated by point A in the curves. Figure 5 shows the second type of the solution family of the static equilibria in the cases of κ = 1/20m −1 , 1m −1 , respectively. These solution curves show that, when ρ 0 is large enough, 0 , θ 0 , δ 0 will take a unique value, respectively. Another feature of these curves is that in the case of a large κ, the static equilibrium configuration seems to be more sensitive to ρ 0 .
Basu-Mandal et al. [12] investigated the dynamic equilibria of a bicycle moving on a horizontal plane and found that there are four solution families relating to the dynamic equilibria. These solution families were plotted as the curves of θ 0 and δ 0 with respect to arctan (R f ω f /4), where ω f is the angular velocity of the front wheel. Noting that the relationship between ω 0 and ω f is given by ω f = −A 4  Figure 6 shows the four families of the dynamic equilibria of the bicycle in the case of κ = 1/200m −1 , where the configuration values (ρ 0 , 0 , θ 0 , δ 0 ) are plotted with respect to arctan (R f ω f /4). We find that ω f should be always negative for the second and fourth solution families. Just for graphical purposes, the curves of these two families are plotted with respect to the absolute value of ω f . The first solution family is related to the curve AB, where δ 0 < π/2. Point A is the accumulation point with an equilibrium configurationξ a 0 = (x a 0 ,ỹ a 0 , θ a 0 , δ a 0 ) = (0.9896, −0.2720, 0, 1.3448) (see Fig. 4), which can be uniquely obtained from (29) under ω 0 = 0. The second one is represented by curve C DE, where point D is a so-called cusp point [12]. The third and fourth solution families are related to the curves FG and H I J , respectively. We can see that points G and H take the same ρ 0 , 0 as well as the same lean angle θ 0 and steer angle δ 0 (but the rotation directions of the front wheel are reverse). Thus these two solution families can be merged into a combined family if we do not distinguish the sign of ω f . Point I in the fourth solution family is also a cusp point.
The structures of solution families of the dynamic equilibria vary as κ increases. Figure 7 shows the four solution families in the case of κ = 1/20m −1 . Compared to the case of κ = 1/200m −1 , there is no cusp point in the solution curves. This means that the appearance of cusp point strictly depends on κ. Our numerical investigation shows that the second and fourth solution families approach each other accompanying the increase of κ. When κ is large enough, the two families will separate from each other in another direction.
We also investigate the solution families of the dynamic equilibria when κ is sufficiently large. Figure 8 shows the solution curves for the case of κ = 1m −1 . Generally, the curve profiles change a lot. We find that the third solution family is still merged with the fourth solution family at point G(H ). Additionally, point B in the first family and point F in the third family take the same values for θ 0 and δ 0 , but different values for ρ 0 and 0 . We still use symbol A to represent the intersection point between the solution families of the static and dynamic equilibria. This point is always located at start of the first family of the dynamic equilibria. Now, let us numerically investigate the Lyapunov stability of the relative equilibria. Based on Proposi- tion 1, we know that a static equilibrium v 0 is always unstable unless all the eigenvalues of the Jacobian matrix J r (v 0 ) are located at the imaginary axis of the complex plane. We confirm this point by checking the eigenvalues of J r (v 0 ) for each static equilibrium in the cases of κ = 1/200m −1 , 1/20m −1 , 1m −1 , and find that it indeed always has an eigenvalue with positive real part. According to Theorem 2, a sufficient condition for a stable dynamic equilibrium v 0 is that the associated Jacobian matrix J r (v 0 ) has six eigenvalues with negative real parts. Based on our numerical investigations, we find that this condition can be satisfied in a small section neighboring the cusp point D in the second solution family when κ = 1/200m −1 . To demonstrate the stable behaviors of the relative equilibria, we select a dynamic equilibrium: v 1 0 = (0.9894, 2.1544, 0.7374, 2.8370, 0, 0, 18), whose Jacobian matrix J r (v 1 0 ) takes the eigenvalues as follows: 0, −12.2845, −0.4277, −0.1407±6.7094i, −0.0111 ± 2.3502i. It is clear that a zero eigenvalue appears in J r (v 1 0 ), and the other six eigenvalues take negative real parts. We can confirm that v 1 0 should be a stable equilibrium point. Meanwhile, we can also conclude that, according to Property 3 induced by the lateral symmetry of Eq.  Figure 9 shows the evolutions ofx,ỹ, θ, δ,φ r with respect to time t. We can see that as t tends to +∞, all of these variables converge to constants, corresponding to a nearby relative equilibrium v 0 = (0.9896, 2.1490, 0.7351, 2.8357, 0, 0, 17.9363) in the same solution family. These results clearly show that v 1 0 is indeed stable, but not asymptotically stable. As an example for the unstable relative equilibrium, we select a dynamic equilibrium in the first solution family when κ = 1/200m −1 : v 2 0 = (0.9894, −6.7698, −0.3887, 0.1431, 0, 0, 18). The eigenvalues of the Jacobian matrix J r (v 2 0 ) are computed as follows: 0, −13.3958, 0.3845, −1.7317 ± 5.7338i, −0.0057±0.9769i. Clearly, there is an eigenvalue with positive real part, and thus v 2 0 is unstable. To see the unstable behavior of system (20) around v 2 0 , we set an initial condition v(0) = (0.9894, −6.7698, −0.3887, 0.1431, −0.1, 0, 18) close to v 2 0 . Figure 10 shows the evolutions ofx,ỹ, θ, δ,φ r with respect to  Fig. 10 Evolutions ofx,ỹ, θ, δ,φ r with respect to time t under an initial condition close to an unstable relative equilibrium in the case of κ = 1/200m −1 time t in the range of t ∈ [0, 13.5]s. We can see that the absolute value of the lean angle θ increases quickly, meaning that the bicycle tends to a fall state. In addition, although the value ofx changes a little, the absolute value ofỹ decreases quickly. This implies that the altitude of the bicycle on surface decreases, and thus its gravitational potential energy also decreases. Due to the conservation of the total energy, the bicycle's kinetic energy must increase. This is consistent with the quick increase ofφ r as shown in Figure 10.

Conclusions
In this paper, we study the symmetry and relative equilibria of the uncontrolled Whipple bicycle on a surface of revolution. This system is a nonholonomic Voronets system, and its symmetry group H is a three-dimensional Abelian Lie group. By applying the Voronets equations to the bicycle dynamics, we characterize the bicycle dynamics as a seven-dimensional reduced dynamic system with time-reversal and lateral symmetries on the reduced constraint space D/H . The relative equilibria of the reduced dynamic system can be divided into two types, corresponding to the bicycle's static equilibrium state and uniform circular motion, respectively. These two kinds of equilibria, under a regular condition that their associated Jacobian matrices have full ranks, should form one-parameter solution families. Additionally, there may be an accumulation point induced by the solution family of the static equilibria intersecting with that of the dynamic equilibria. In order to study the Lyapunov stability of a given relative equilibrium, we first discuss the basic properties of the Jacobian matrix at the equilibrium point. One of the most important points is that the Jacobian matrix must have a zero eigenvalue, which means that the relative equilibrium is not hyperbolic. Noting that the reduced constrained energy of the system is conservative, which can be considered as a first integral of the reduced dynamic system, we can then analyze the stability of the relative equilibrium by restricting the motion of the system on the invariant manifold. For a relative equilibrium whose Jacobian matrix takes six eigenvalues with negative real parts, we prove that it is Lyapunov stable, but cannot be asymptotically stable. And the solution of the reduced dynamic system under the initial condition close to this relative equilibrium will converge to a nearby relative equilibrium.
Based on symbolic computation, we carry out numerical simulations using the reduced dynamic equations in conjunction with the benchmark parameters of the Whipple bicycle. How the parameter of the surface of revolution affects the relative equilibria and their stabilities of the bicycle system is then numerically investigated in this paper.
In summary, we perform a comprehensive analysis about the symmetry and relative equilibria of the uncontrolled Whipple bicycle system on the surface of revolution. These results provide theoretical bases for a better understanding of the intrinsic properties of the bicycle dynamics and other nonholonomic systems. The theoretical results may also help to develop controllers for these systems.