Symmetry and relative equilibria of a bicycle system movingon a revolution surface

Finding the relative equilibria and analyzing their stabil ities are of great significance to revealing the intrinsic properties of mechanical systems and develop ing effective controller. In this paper, we study the symmetry and relative equilibria of a bicycle system moving o a revolution surface. We note that the symmetry group of the bicycle is a three-dimensional Abelian Lie grou p, and the rolling condition of the two wheels produces four time-invariant first-order linear constrain ts to the bicycle system. Therefore, we can classify the bicycle dynamics as a general Voronets system whose Lagr angian and constraint distribution are kept invariant under the action of the symmetry group. Applying t he Voronets equations to the bicycle dynamics, we obtain a seven-dimensional reduced dynamic system on the re duced constraint space. This system takes timereversal and lateral symmetries, and has two kinds of relati v equilibria: the static equilibria and the dynamic equilibria. Further theoretical analysis shows that both k inds of relative equilibria form one-parameter solution families, and their Jacobian matrices take some specific pro perties. We then show that a static equilibrium cannot be stable unless all the eigenvalues of the Jacobian m atrix are located at the imaginary axis of the complex plane. The stability of the dynamic equilibria is st udied by limiting the reduced dynamic system to an invariant manifold, which is established based on the con servation of energy of the system. We prove in a strict mathematical sense that the dynamic equilibria may b e Lyapunov stable, but cannot be asymptotically stable. Finally, we employ symbolic computation to carry ou t n merical simulations in conjunction with the benchmark parameters of a Whipple bicycle. How the revoluti n surface a ffects the relative equilibria and their stabilities is then investigated through our numerical sim ulations.

moves on a ground, the non-slipping 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 take the first-order linear forms. This leads to a nonholonomic Voronets system [23]. At each point of its seven-dimensional configuration space Q, the set of all possible velocities satisfying these constraints is a three-dimensional subspace of the tangent space to Q [24,25]. The collection of these subspaces form a constraint distribution D on Q, which is a ten-dimensional submanifold of the tangent bundle T Q. Thus, the velocity phase trajectories of the bicycle are always restricted on D.
Compared to conventional nonholonomic systems such as a rattleback system [26,27] and a rolling disk [28], 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 [29,30] to obtain the holonomic constraint equations of a bicycle with toroidal tires. Wang et al. [31] parameterized the two wheel profiles and used an extremum condition that the wheel-road contact should be located at the lowest point on the wheel edge. Xiong et al. [32,33] employed an ordered process presented in [28] to obtain the constraint equations of a bicycle moving on a horizontal plane and an arbitrary surface, respectively.
For the bicycle dynamics, its dynamical equations can be established using the Newton-Euler approach in classical mechanics [10][11][12]. The bicycle dynamics can also be analyzed by way of analytical mechanics such as the Euler-Lagrange equations [6,12,31,33], the Kane's method [30] and the Gibbs-Appell equations [32]. In order to analyze the mechanical properties of a Whipple bicycle moving on the horizontal ground, Boyer et al. [17,34] derived the reduced dynamic equations using a general reduction-based approach in geometric mechanics. In this method, the bicycle's configuration space is modelled 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 disadvantage of the method is that the kinematic constraints are not invariant under the left action of the structural group SE(3), thus some reconstruction equations are needed in the reduced system.
It is worth noting that different methods produce different forms of equations of motion, such as the ordinary differential equations (ODEs) [11,32] and the differential algebraic equations (DAEs) [12,33]. This will lead to different difficulties when solving the relative equilibria and analyzing their stabilities. Pappalardo et al. [35] proposed a method for analyzing the stability characteristics of multibody mechanical systems with holonomic and nonholonomic constraints based on the direct linearization of the index-three form of the DAEs. García-Agúndez et al. [36] developed three novel general numerical procedures to obtain the linearized equations of any multibody system: the linearization of the index-one form of the DAEs, the reduction of the linearized equations from DAEs to ODEs form, and a further reduction of the linearized equations by eliminating the holonomic constraints. Using these numerical approaches, they analyzed the linear stability of the bicycle's uniform straight forward motion. Their results coincided with those obtained by Meijaard et al. [11] using a pair of second-order linearized ODEs.
In this paper, we first study the motion of the Whipple bicycle on a surface of revolution. In this case, the bicycle is a system with a three-dimensional Abelian Lie group H whose dimension is less than the number of the velocity constraints. Because all the velocity constraints of the bicycle take the first-order linear forms, the bicycle can be atrributed to a Voronets system, whose dynamics are governed by Voronets equations (VEs) [23,37]. Theoretically, VEs are equivalent to the Gibbs-Appell equations [32] and Kane's equations [30] when they are used to deal with nonholonomic systems. The advantage of VEs is that they can reveal the internal characteristics of the system dynamics by way of clear geometric structures, such as Ehresmann connection and its curvature [25]. Based on VEs, we obtain a seven-dimensional reduced dynamic system on the quotient manifold of D under the action of H [24,25,38], and find that this system takes time-reversal and lateral symmetries [11]. More importantly, we show that the system has two kinds of relative equilibria: one is the static equilibria related to the bicycle's static equilibrium configuration, and the other is dynamic equilibria representing the uniform circular motion of the bicycle. According to the internal characteristics of the bicycle dynamics exposed by VEs, we also prove that these relative equilibria form one-parameter solution families under certain regular conditions. Furthermore, we discuss the Lyapunov stability of these relative equilibria through analyzing the properties of the associated Jacobian matrix. Noting that all the relative equilibria are not hyperbolic, we conclude that a stable equilibrium of the bicycle system cannot be asymptotically stable [39]. Further analysis shows that the disturbed motion of the bicycle system around a stable relative equilibrium will converge to a nearby relative equilibrium.
The rest of this paper is organized as follows: In Sect. 2, we analyze the symmetry of the Whipple bicycle system on the revolution surface 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 will discuss the symmetry of an uncontrolled bicycle system on the revolution surface and derive the reduced dynamic system on the reduced constraint space using the Voronets equations. Figure 1 adapted from [11,31] 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.

Configuration description
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. 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 revolution surface where f takes the form f (ρ) ≡ 0, ∀ρ ∈ [0, +∞).
A movement of the bicycle on Σ is subject to two geometric constraints and four velocity constraints. Both types of constraints are 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 were given in our previous work [33]. As a result, the configuration space of the bicycle is a seven-dimensional 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, rotation of the rear frame in the horizontal plane (i.e., yaw motion) and rotation of the front wheel, respectively.
The left action of H on Q is defined as a smooth mapping Φ : For each h ∈ H, (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(ψ), Clearly,q consists of two parts: the coordinates (θ, δ,x,ỹ) 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 byσ = (θ,δ,φ r ) the set of three independent velocities. The four first-order linear constraint equations can then be written as follows: whereṡ = (ẋ,ẏ,ψ,φ f ). 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 Σ 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 [25]: (1) The Lagrangian is invariant, i.e., L(vq) = L(h * vq), ∀vq = (q,q) ∈ T Q, ∀h ∈ H, where h * : T Q → T Q represents the tangent map induced by the action of h on Q. (2) The constraint distribution is invariant, i.e., h * Dq = D hq , ∀q ∈ Q, ∀h ∈ H.
According to (8), the Lagrangian and constraint distribution invariances mean that φ r , ψ, φ f do not appear explicitly in the expressions of the 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 hη = hq,σ .

Reduced dynamic system
We employ the Voronets equations [25] to obtain the motion equations of the uncontrolled Whipple bicycle on Σ. Denote by L c = L|D : D → R the constrained Lagrangian. According to (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 (13), which is a necessary consequence of the Lagrangian and constraint distribution invariances.
The Voronets equations take the following forms [23]: where with b = 1, · · · , 4, α, β = 1, 2, 3, are anti-symmetric about the two indices α and β. From the point of geometric mechanics view [25], 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, (15) can be written as follows: where c αβγ and P α are functions ofξ, i.e., Noting the Lagrangian L and the constrained Lagrangian L c in (9) and (14), respectively, and applying the Voronets equations (15), some specific coefficients in (18), which are important for revealing the internal characteristics of the bicycle dynamics, can be expressed as follows: By transforming (17) into first-order forms, and combining the first two velocity constraint equations in (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 the 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 This means that (20) has time-reversal symmetry, i.e., if v(t) is a solution of (20), so isṽ(−t).
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 revolution surface. (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 [25].
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 the reduced dynamic system (20) by setting Y(v 0 ) = 0. According to (21), they must take the form v 0 = (x 0 ,ỹ 0 , θ 0 , δ 0 , 0, 0, ω 0 ) = (ξ 0 , 0, 0, ω 0 ), where there are five constants to be determined. From (21), we have the following five algebraic equations: Noting that the reduced dynamic system (20) takes the time-reversal and lateral symmetries, we can claim that, if v 0 is a relative equilibrium, so We can classify the relative equilibria into two types, depending on weather ω 0 is zero or not. We give the following definition: (20) is called a static equilibrium if ω 0 = 0, and is called a dynamic equilibrium if ω 0 0. 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 (25) are always satisfied, and the static equilibrium configurationξ s 0 of the bicycle can be solved through the last three equations in (25), namely, It is clear that the number of the independent algebraic equations in (26) is less than the number of the elements inξ s 0 . Thus in general, the solution of (26) is not unique.
Let us define J s as the Jacobian matrix of (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,ỹ c 0 is the dynamic equilibrium configuration. In principle, we can uniquely determine the five values of (ξ c 0 , ω 0 ) since there are five equations in (25). However, further analysis shows that the last equation in (25) is an identical equation. This point can be illustrated below.
If ω 0 0, the first two equations of (25) can be equivalent to Thus the last equation in (25) (α = 3) is an identical equation, and the equilibrium equations (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 (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.
Theorem 1 Suppose that v 0 = (ξ 0 , 0, 0, ω 0 ) represents the relative equilibrium of the 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 below.
According to the 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 (7) to obtain the motion equations of x(t) and y(t) under the initial condition ψ(0) = 0: where ρ 0 is the radius of the trajectory of point D, Θ 0 is the phase difference between x(t), y(t) and ψ(t) =ψ 0 t, and these two constants are given by the following equations: Obviously, this corresponds to a situation where the bicycle with a constant angular velocityψ 0 and a constant configuration of lean angle θ 0 and steer angle δ 0 moves along a circular trajectory on surface Σ. It is obvious that equation (29), under the condition ω 0 = 0, can give a unique solution for the dynamic equilibrium configuration. We denote this configuration byξ a 0 , and we can easily verify thatξ a 0 also satisfies (26), the conditions relating to a static equilibrium configuration (A 1 3 = A 2 3 = 0 corresponds to P 3 = 0). Therefore,ξ a 0 is a common configuration that belongs to bothξ as ω 0 approaches zero if J c has full rank at v a 0 . Figure 2 shows the intersection structure of the two kinds of relative equilibria. The blue line represents the one-parameter solution family of the static equilibria, wherex 0 is considered as a single parameter. The red line represents the one-parameter solution family of the dynamic equilibria, where ω 0 is set as the single parameter. These two solution families intersect at point v a 0 on the hyperplane ω 0 = 0.

Mathematical properties of the Jacobian matrices
In this section, we discuss the essential properties of the Jacobian matrix J r (v 0 ) of (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 equilibria (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 (20) at a relative equilibrium v 0 = (ξ 0 , 0, 0, ω 0 ) can be written as Noting the expressions of Y i listed in (21) and the structure of v 0 , we can present the concrete expression of J r (v 0 ) as follows: where Proof. If J r (v 0 ) is rank-deficient, the property is naturally true. Based on (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: When v 0 is a static equilibrium, we have [ where J s is given by (27). Thus, meaning that J r (v 0 ) in this case is rank-deficient. In particular, if v 0 is an accumulation point, we have [A 3 ] 2×1 = [0] 2×1 , thus rank J (v 0 ) ≤ 3. In this case, J r (v 0 ) has at least two zero eigenvalues.
If ω 0 0, we know from (28) and (29) (19), the partial derivatives of c 333 and P 3 take the following forms, 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 (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 rank-deficient, thus it should have at least one zero eigenvalue.
Proof. Based on (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.
Proof. According to (24), the entries of the Jacobian matrices around v 0 andv 0 satisfy the following relationships, This means that there exists a transform matrix T(v 0 ,v 0 ), . 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 the reduced dynamic system (20) in a strict sense of mathematics. We first discuss the stability of static equilibria based on the time-reversal symmetry of (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 (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 For a static equilibrium whose Jacobian matrix has at least one eigenvalue that is not located at the imaginary axis of the complex plane, it must be unstable.
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 equilibrium v 0 = (ξ 0 , 0, 0, ω 0 ), where ω 0 0. 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 [39]. 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 (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 the 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 , E 1 ) = E(u,φ r )− E 1 . 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 equationẼ(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 , E 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 (20) with the coordinates u.
According to the theory of nonlinear dynamics [39,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 (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.

Corollary 1
For a relative equilibrium v 0 , suppose that J r (v 0 ) has six negative real part eigenvalues, then, (i) v 0 must be a dynamic equilibrium, i.e., ω 0 0; (ii) the zero eigenvector of J r (v 0 ) is parallel to the tangent vector of the family curve at v 0 .
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 the 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 (20) on Q E 0 , we establish a six-dimensional reduced system with coordinates u, which can be expressed as follows: u   = X(u, F(u, E 0 )). (41) Clearly, u 0 is an equilibrium of (41), at which the Jacobian matrix of (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 (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 equation (38) with respect to variable E 1 , (44) can be simplified as It is worth noting that the Jacobian matrix J r (v 0 ) of the 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 the 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 the system (20). Furthermore, the solution v(t) of (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 the 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 =Ê(s) = E(v(s)). The derivative ofÊ(s) with respect to s at s = 0 is given by Differentiating (38) with respect u at point ζ 0 , we have Substituting equations (47) and (49) into equation (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 neighborhoodÛ(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 (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 (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  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 (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 (9), the nonholonomic constraint equations (10) as well as the Voronets equations (17). The coefficients involved in (25) can then be determined, and equations (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 one-parameter 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. In order to present a clear description for the equilibrium configuration of the bicycle, we use the variables ρ 0 and Θ 0 defined in (32) to replace the variablesx 0 andỹ 0 . Clearly, (ρ 0 , Θ 0 ) are polar coordinates of (x 0 ,ỹ 0 ).
In the following, we select ρ 0 and ω 0 as the independent parameter for the solution families of static and dynamic equilibria, respectively. Due to the time-reversal and lateral symmetries of the reduced dynamic system (20), we compute the solution families by limiting ω 0 ∈ [0, +∞) and δ 0 ∈ [0, π] without loss of generality. Figure 3 shows the solution family of the static equilibria in the case of κ = 1/200, where the curves of Θ 0 , θ 0 , δ 0 are plotted with respect to ρ 0 . 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  Figure 4 shows the solution families of the static equilibria in the cases of κ = 1/20, 1. Compared to the case of a small κ, these solution curves show that, when ρ 0 is large enough, Θ 0 , θ 0 , δ 0 will take a unique value, respectively. Another feature 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 3 (x 0 ,ỹ 0 , θ 0 , δ 0 )ω 0 , we use the same way as in [12] to plot the curves of solution families of the dynamic equilibria. Figure 5 shows the four families of the dynamic equilibria of the bicycle in the case of κ = 1/200, 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 values of ω f . The first solution family is related to the curve AB, where δ 0 < π/2. The second one is represented by curve CDE, where point D is a so-called cusp point [12]. The third and fourth solution families are related to the curves FG and HIJ, 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.
Point A is an accumulation point with an equilibrium configurationξ a 0 = (x a 0 ,ỹ a 0 , θ a 0 , δ a 0 ) = (0.9896, −0.2720, 0, 1.3448), which can be uniquely obtained from (29) under ω 0 = 0. Noting that the point also satisfies the condition of the static equilibrium, it corresponds to the intersection of the curves of the static equilibria and the first solution family of the dynamic equilibria (see Fig. 3). The structures of solution families of the dynamic equilibria vary as κ increases. Figure 6 shows the four solution families in the case of κ = 1/20. Compared to the case of κ = 1/200, there is no cusp points 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 7 shows the solution curves for the case of κ = 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 Proposition 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/200, 1/20, 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/200. 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 (20), pointv 1 0 = (0.9894, −2.1544, −0.7374, −2.8370, 0, 0, 18), should also be a stable dynamic equilibrium.  In order to demonstrate the dynamical behavior of the bicycle around a stable relative equilibrium, we simulate (20) using the ode45 function of MATLAB with a fixed time step △t = 0.001s by setting an initial condition v(0) = (0.9894, 2.1544, 0.7374, 2.8370, −0.1, 0, 18) close to v 1 0 . Figure 8 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/200:  Figure 9 shows the evolutions ofx,ỹ, θ, δ,φ r with respect to 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.  Fig. 9 Evolutions ofx,ỹ, θ, δ,φ r with respect to time t under an initial condition close to an unstable relative equilibrium in the case of κ = 1/200.
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 9.

Conclusions
In this paper, we study the symmetry and relative equilibria of the uncontrolled Whipple bicycle on the revolution surface. 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 carried out numerical simulations using the reduced dynamic equations in conjunction with the benchmark parameters of the Whipple bicycle. How the revolution surface affects the relative equilibria and their stabilities of a 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 controls of these systems.