CC BY
14Russian Journal of Nonlinear Dynamics, 2021, vol. 17, no. 4, pp. 391-411. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd210403
NONLINEAR PHYSICS AND MECHANICS
MSC 2010: 37J60, 70H14
Symmetry and Relative Equilibria of a Bicycle System
J.Xiong, Y.-B.Jia, C.Liu
In this paper, we study the symmetry of a bicycle moving on a flat, level ground. Applying the Gibbs-Appell equations to the bicycle dynamics, we previously observed that the coefficients of these equations appeared to depend on the lean and steer angles only, and in one such equation, a term quadratic in the rear wheel's angular velocity and a pseudoforce term would always vanish. These properties indeed arise from the symmetry of the bicycle system. From the point of view of the geometric mechanics, the bicycle's configuration space is a trivial principal fiber bundle whose structure group plays the role of a symmetry group to keep the Lagrangian and constraint distribution invariant. We analyze the dimension relationship between the space of admissible velocities and the tangent space to the group orbit, and then employ the reduced nonholonomic Lagrange-d'Alembert equations to directly prove the previously observed properties of the bicycle dynamics. We then point out that the Gibbs-Appell equations give the local representative of the reduced dynamic system on the reduced constraint space, whose relative equilibria are related to the bicycle's uniform upright straight or circular motion. Under the full rank condition of a Jacobian matrix, these relative equilibria are not isolated, but form several families of one-parameter solutions. Finally, we prove that these relative equilibria are Lyapunov (but not asymptotically) stable under certain conditions. However, an isolated asymptotically stable equilibrium may be achieved by restricting the system to an invariant manifold, which is the level set of the reduced constrained energy.
Keywords: bicycle, nonholonomic system, symmetry, reduced system, relative equilibria, Lyapunov stability
Received May 06, 2021 Accepted October 17, 2021
This work was supported by the grant of the National Natural Science Foundation of China (NSFC:11932001).
Jiaming Xiong [email protected] Caishan Liu [email protected]
Peking University
No. 5 Yiheyuan Road Haidian District, Beijing, P.R.China 100871
Yan-Bin Jia [email protected]
Iowa State University Ames, IA 50010, USA
1. Introduction
Studies on bicycle dynamics started at the end of the 19th century with the seminal work by Carvallo [1] and Whipple [2]. Since then this topic has drawn attention from J.Boussi-nesq [3], F.Klein and A. Sommerfeld [4], S. P. Timoshenko and D.H.Young [5], and many recent researchers from different areas including classical mechanics [6-14], applied mathematics [1517], system control [18-20], etc. Meijaard et al. [11] did a comprehensive review of the history of bicycle dynamics, while Kooijman and Schwab [21, 22] presented the state of the art of the literature.
Most of the studies focused on a benchmark model proposed by Whipple [2] and Carvallo [1], in which they simplified the bicycle as a rigid four-body system with knife-edge wheels. When the bicycle moves on a flat, level ground, the non-slipping contact between its two wheels and the ground induces two normal geometric constraints and four tangent velocity constraints. This leads to a nonholonomic system [23, 24], i. e., one whose motion needs to be described using more coordinates than the dimensions of its velocity space. Compared to conventional nonholonomic systems such as the Chaplygin sleigh [25] and a rolling disk [26], 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].
Many methods were developed to establish these constraint equations. 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 [27, 28] to obtain the constraint equations of a bicycle with toroidal tires. Wang et al. [29] 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. Recent derivations by Xiong et al. [30, 31] employed an ordered process proposed by Zhao et al. [26].
An outstanding phenomenon to researchers is that an uncontrolled bicycle can balance itself at the right speed [4, 11, 13, 32, 33]. A proper understanding has to be sought from within the bicycle's dynamic model. From classical mechanics, the bicycle's motion can be characterized using the Newton-Euler approach [10-12], the Euler - Lagrange equations [6, 12, 29], or analytical mechanics such as the Kane method [28] or the Gibbs-Appell equations [30]. In the early years, researchers focused on a linear model obtained by slightly disturbing the bicycle from an upright straight motion [10, 11]. This simplified model offered a proper explanation for the bicycle's self-stability, as verified through experiments [19, 34]. The computer code for multibody dynamics was later used to generate nonlinear models [14, 35-38] based on which the bicycle's circular motion could also be analyzed [7, 39-41]. Basu-Mandal et al. [12] presented for the first time a complete explanation of the circular motion of the Whipple bicycle. They used the Euler - Lagrange method and applied both symbolic manipulation and numerical simulation to obtain a solution to the hands-free circular motion. Although four different one-parameter families exist in the circular motion of this type of bicycle, two of them could be merged into one, leaving three in total. This result was verified by Wang et al. [29] and Xiong et al. [30].
Both the upright straight and circular motions can be regarded as relative equilibria [42, 43] of the bicycle dynamics. As mentioned earlier, these relative equilibria are not isolated, but form several one-parameter solution families. There are many nonholonomic systems that have similar properties. Zenkov et al. [44-46] presented a collection of examples that have nonisolated relative equilibria, such as the Routh problem, a rolling disk, a roller racer and a rattleback. For
the bicycle, this property arises from the special structure of its dynamics. It can be stated as follows in terms of the Gibbs-Appell equations derived in our previous work [30]:
(1) the coefficients of these equations depend on the lean and steer angles only,
(2) in one of the equations, both the term quadratic in the rear wheel's angular velocity and the pseudoforce term vanish.
The above structure has been found via symbolic manipulations [29, 30]. However, no formal proof has been given. Neither the Euler-Lagrange equations nor the Gibbs-Appell equations from classical mechanics seem able to reveal this structure directly, especially its second property, since these equations do not use the symmetry of the bicycle system. The only symmetry that might be considered is the cyclic coordinate [11], which is only a special case where the symmetry group is R1 or S1, i.e., an abelian group [43].
In this paper, we systematically analyze the symmetry of the bicycle system on a flat, level ground utilizing basic concepts of geometric mechanics such as Lie group and Lie algebras [42, 43, 46]. The configuration space of the bicycle is chosen as a (trivial) principal fiber bundle [46-48], which is common in robotic locomotion concerning a two-wheeled planar mobile robot [47, 4951], a snakeboard [49, 50, 52, 53] or a kinematic snake robot composed of three rigid links [54, 55]. Boyer et al. [17, 53] derived the reduced dynamic equations of a bicycle using a general reduction-based approach in the principal fiber bundle, and presented several numerical examples. Nevertheless, these reduced equations did not reveal the special structure of the bicycle system directly.
To prove the existence of the special structure, we consider the left action of the structure group on the bicycle's configuration space, which keeps the Lagrangian and constraint distribution invariant. We analyze the dimension relationship between the tangent space at an arbitrary point and its subspaces, i.e., between the space of admissible velocities and the tangent space to the group orbit. Under this condition, the so-called reduced nonholonomic Lagrange-d'Alembert equations [46, 48] can be employed, and the aforementioned properties can be directly proved. Furthermore, we point out that the Gibbs-Appell equations give the local representative of the reduced dynamic system on the reduced constraint space, i. e., the quotient manifold of the constraint distribution under the action of the structure group. Thus, the bicycle's uniform upright straight and circular motions are related to the equilibrium points of the reduced system. These points are termed as relative equilibria. Finally, we analyze the stability of these relative equilibria with the help of the invariant manifolds of the reduced system under the conserved reduced constrained energy.
The rest of this paper is organized as follows: In Section 2, we present a kinematic analysis of the Whipple bicycle under geometric and velocity constraints. In Section 3, we analyze the symmetry of the bicycle system in terms of the group action. We concretely illustrate the special structure of the bicycle dynamics in Section 4. In Section 5, we prove this structure based on the reduced nonholonomic Lagrange-d'Alembert equations. In Section 6, we study the relative equilibria of the reduced system. We analyze the Lyapunov stability of these relative equilibria in Section 7. Conclusions are offered in Section 8.
2. Whipple bicycle and its constraint equations
In this section, we will present a kinematic analysis for the Whipple bicycle, which is shown in Fig. 1 adapted from [11] as moving on a flat, level ground. The bicycle 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.
Whipple bicycle
Head structure
Fig. 1. Whipple bicycle (in upright configuration)
We choose the bicycle's upright configuration as the reference state. In the configuration, the bicycle rests on the horizontal plane nh such that all of its four bodies are symmetric with respect to a common vertical plane n. Below we define five coordinate frames, one inertial and four body-fixed, each with its axes 1 and 3 inside n, and axis 2 orthogonal to it. The inertial frame Fi = {O; i1; i2, i3} has its origin in nh (which is thus spanned by its axes 1 and 2). The rear and saddle frames Fr = {Or; er 1, er 2, er 3} and Fs = {Os; es 1, es 2, es 3} are oriented the same as Fi and located, respectively, at the center Or of the rear wheel and the center of mass Os of the saddle structure. The head frame Fh = {Oh; eh 1, eh 2, eh 3} is located at the center of mass Oh 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 Ff = {Oh; ef 1, ef 2, ef 3} is located at the center Of of the front wheel and has the same orientation of the head frame Fh. All five frames are illustrated in Fig. 1. The plane n^ of symmetry always coincides with the planes er 1-er 3 and es 1-es 3 in the frames Fr and Fs, respectively.
Without any constraints, the accessible configuration space of the bicycle in a free motion is a nine-dimensional manifold,
Qfree = R3 x SO(3) x 51 x 51 x 51. This manifold has five components described in detail below:
• The first component R3, with three coordinates (x, y, z), describes the motion of a reference point D, which is located at the intersection of the steering axis and the coordinate axis Of ef 1 of the front frame Ff.
• The special orthogonal group SO(3) describes the orientation of the saddle frame Fs relative to the inertial frame Fi. Given a basis {t1, t2, t3} of the Lie algebra so(3) of SO(3), where t1 = (1, 0, 0), t2 = (0, 1, 0), t3 = (0, 0, 1), and 1 R3 ^ so(3) is the hat map, we define the following map y: so(3) ^ SO(3),
Y (&Ti + + #3) = exp(0f3) exp(0?i) exp(yr2). (2.1)
It is easy to verify that T0y = idso(3), which means that y is a local diffeomorphism from a neighborhood of zero in so(3) onto a neighborhood of the identity element in SO(3). Therefore, we can choose (0, d, y) as a set of local coordinates of SO(3), which corresponds to the set of three Euler angles in 312-sequence.
• The last three components S1 x S1 x S1 of Qfree, with coordinates (5, , 0f), represents the rotation of the head frame Fh about the steering axis out of the plane n, and the rotations of the rear and front frames Fr and Ff about the axes 2 of the saddle frame Fs and the head frame Fh, respectively.
In summary, we denote by qfree = (x, y, z, 0, d, y, 5, , ) the set of local coordinates for the bicycle in the free motion.
A movement of the bicycle on the plane nh is subject to two geometric constraints and four velocity constraints to be introduced below. Suppose that the rear and front wheels are in contact with the ground at points Pr and Pf, respectively. On the edges of the two wheels, the two points instantaneously coinciding with Pr and Pf, respectively, are located by parameters ur and Uf of the wheel curves. These curve parameters can be defined as the rotation angles along the directions of the spin axes er>2 and ef 2 of the two wheels. The position vectors of Pr and Pf in the inertial frame Fi, denoted by rr and f, can then be written as
ri = ri (qfree,Ui), i = r,f. (2.2)
The contact points should satisfy the following two geometry conditions [26, 51]: (i) they lie on nh; and (ii) the tangents to the edges of the two wheels at the contact points are perpendicular to the normal vector of nh. These two conditions can be mathematically expressed as follows:
ri • 13 = ^
dvl . i = r,f. (2.3)
The curve parameters ur and uf can be solved from the second equation in (2.3) as the functions of the local coordinates,
ur = ur (qfree), uf = uf(qfree). (2.4)
Substituting (2.4) into the first equation of (2.3), we obtain two geometric constraint equations which can be locally and implicitly written as
z = z(q), y = <p(q), (2.5)
where q = (d, 5, , x, y, 0, ) has seven independent coordinates.
No slip of the two wheels occurs. By calculating the tangential velocities of the contact points Pr and Pf and letting them be zero, we obtain four velocity constraint equations in q,
A(q) ■ q = 0. (2.6)
Under constraints (2.5) and (2.6), the degrees of freedom of the bicycle moving with no slip equals 9 — 2 — 4 = 3.
3. Symmetry of the bicycle system
In this section, we will analyze the symmetry of the bicycle system based on the kinematic analysis in the preceding section. When the two geometric constraints (2.5) are only considered, the configuration space of the bicycle is a seven-dimensional manifold,
Q = 51 x 51 x 51 x SE(2) x 5
with the coordinates (d, S, , (x, y, ) introduced in Section 21.
Denoting by R = 51 x 51 the base space, and by G = 51 x SE(2) x 51 the structure group, we know that the configuration space Q = R x G is a trivial principal fiber bundle [46-48], with the projection n : Q ^ R given by
n(q) = n(r, g) = r,
(3-1)
where q = (r, g) with r G R and g G G.
The bicycle system is also subject to four velocity constraints (2.6). At any configuration q G Q, the set of all possible velocities satisfying these constraints is the subspace of the tangent space to the configuration space at q,
Vq = {v G TqQ | A(q) • v = 0}.
(3.2)
Denote by q = (a, s) another decomposition of the coordinates, where a = (a1, a2, a3) = = (9, 5, ) e R x G1 and s = (s1, s2, s3, s4) = (x, y, ^, 4>f) e G2 x G3. We solve (2.6) for S in terms of a and rewrite the solution as
+ Bj(q)(F =0, i = 1, ..., 4,
(3.3)
where the coefficients Bj, i = 1, ..., 4, j = 1, 2, 3 are the functions of q, and Einstein's summation convention is used for the index j.2 With these notations, is given by
d ■ d
ds%
j = 1, 2, 3 .
(3.4)
1 Strictly speaking, R is a submanifold of S1 x S1 for the reason that the bicycle should be away from the configuration in which it falls to the ground.
2According to this notation, a single term containing repeated indices (no more than twice) implies
summation of that term over all the values of the index.
1
Thus we can obtain a constraint distribution of Q, which is defined as the collection of Dq,
D = U Dq. (3.5)
q
qeQ
The velocity constraints can now be equivalently expressed as q(t) G Dq(t) for all t.
We introduce the left action of the Lie group G on the configuration space Q as follows:
$h(q):= hq = (r,hg), Vh G G, Vq = (r, g) G Q. (3.6)
More specifically, for h = (01, a, b, a, 02) and q = (9, 5, 0r, x, y, 0, 0f), we have
hq = (9, 5, 0r + 01, x cos a — y sin a + a, x sin a + y cos a + b, 0 + a, 0f + 02). (3.7)
It is clear that the quotient space Q/G, whose points are the group orbits, is diffeomorphic to R, i.e., Q/G = R, with the projection n: Q ^ Q/G equivalently given by (3.1).
Note that G is the direct product of its three Lie subgroups, i.e., G — G1 x G2 x G3. where G1 = S1 (rear wheel), G2 = SE(2) and G3 = S1 (front wheel). The action of G2 on the configuration space Q is equivalent to a corresponding translation or rotation of the inertial frame Fi in the plane nh, while the spatial position of the bicycle remains unchanged. The Lagrangian and constraint distribution are invariant under the action of G2 since their expressions are independent of the choice of Fi in the plane nh. These two expressions are also invariant under the actions of G1 and G3 for no appearance of the rotation angles of the two wheels, which are circular (and thus symmetric). We claim that under the action of the Lie group G = G1 x G2 x G3 the bicycle system has the following two invariances:
(1) The Lagrangian is invariant, i.e., L(vq) = L(g*vq), Vq G Q, Vg G G, Vvq G TqQ, where g*: TQ ^ TQ represents the tangent map induced by the action of g on Q.
(2) The constraint distribution is invariant, i.e., g*Dq = Dgq, Vq G Q, Vg G G.
The above claim reveals the symmetry of the bicycle system on a flat, level ground, where the symmetry group G is nonabelian. These invariance properties can also be directly verified via symbolic manipulation, which requires less computation than deriving the whole dynamic equations. Our claim is more rigorous than simply attributing the generalized coordinates 0r, x, y, 0, 0f to the cyclic or ignorable coordinates [11]. When the bicycle moves on a flat, level ground, there is no doubt that 0r, x, y, 0f must be cyclic coordinates that do not appear in the expressions of the Lagrangian and velocity constraint equations. However, since the yaw angle 0 represents the bicycle's orientation in the plane nh, it may be related to the tangent velocities of the two contact points Pr and Pf, and thus may appear in the expressions of the velocity constraint equations [29]. Therefore, we need to treat it more carefully.
4. Properties of bicycle dynamics
In this section, we will reveal the structure of the bicycle dynamic system based on the Gibbs-Appell equations derived in our previous work [30]. There we selected three pseudove-locities a = (<r1, &2, &3) = (9, 5, 0r) and obtained, for an uncontrolled bicycle, a set of three second-order ordinary differential equations (ODEs), called the Gibbs-Appell equations,
mij+ Ci jk&j&k = Pi, i = 1, 2, 3, (4.1) _RUSSIAN JOURNAL OF NONLINEAR DYNAMICS, 2021, 17(4), 391-411_If)
where mijs are the components of the mass matrix, ci jks are the coefficients of the terms which are quadratic in the pseudovelocities, and Pis are the pseudoforces.
Generally speaking, mij, ci -k and Pi may be the functions of all the seven coordinates of the configuration space Q. In the bicycle system, however, they depend on only two of them. And further simplification of the Gibbs-Appell equations is possible, as stated below.
Theorem 1. The coefficients and pseudoforces of the Gibbs-Appell equations (4.1) have the following properties:
(1) They are the functions of the lean angle 9 and steer angle 5, and are independent of the other coordinates. Namely, we have mij = mij (9, 5), ci jk = ci jk(9, 5) and Pi = Pi(9, 5).
(2) In the third equation of (4.1) (i = 3), the term (<j)r)2 (j = k = 3) and the pseudoforce term always vanish, that is, c3,33(9, 5) = 0 and P3(9, 5) = 0.
We initially observed Theorem 1 via symbolic manipulation using Maple [30]. Wang et al. [29] employed the Euler - Lagrange equations to derive the nonlinear equations of a bicycle, eliminated all the dependent variables and all the Lagrange multipliers also using a symbolic manipulation tool, and then obtained the reduced equations that are equivalent to the Gibbs -Appell equations (4.1). They discovered the first property stated in Theorem 1, and gave the following explanation. The yaw angle ^ must appear in the form of the linear combinations of sin(^) and cos(^), and after multiplication and summation over the matrix, the final Euler -Lagrange equations are independent of Indeed, it may be very difficult, if not impossible, to prove the second property using the classical forms of the Gibbs-Appell equations or the Euler - Lagrange equations. We believe that the first and second properties stated in Theorem 1 both follow from the symmetry of the bicycle system. In Section 5, we will prove them based on the reduced equations on the quotient manifold D/G.
The coordinates q = (r, g) of Q, as well as the pseudovelocities a, form a set of local coordinates of the constraint distribution D. Therefore, the Gibbs-Appell equations (4.1), together with the velocity constraints (3.3), determine the dynamic system on D,
n = X (n), (4.2)
where n = (r, g, a) e D. Clearly, this dynamic system agrees with that derived from the Lagrange-d'Alembert principle. Due to the symmetry of the bicycle system under the action of G, the vector field X must be G-invariant [44, 45], i.e., (h*)*X (n) = X (h* rj), Vn eD, Vh e G. In other words, if n(t) is a solution of (4.2) with the initial condition n(0) = no, then h*n(t) is also a solution of (4.2) with the initial condition h*n0.
5. Proof of Theorem 1
In this section, we will prove Theorem 1 using the invariance properties of the bicycle system, i.e., the Lagrangian invariance and constraint distribution invariance introduced in Section 3. The main idea is to employ the reduced nonholonomic Lagrange-d'Alembert equations based on the nonholonomic connection on the principal fiber bundle n: Q Q/G = R [46, 48]. These reduced equations, coming from the symmetry of the nonholonomic mechanics system, will present a more clear structure than the classical Gibbs-Appell equations.
The kernel of the tangent map Tqn is called the vertical space of the principal fiber bundle at a point q G Q, and is denoted by Vq [46]. This space consists of infinitesimal generators of the group action G at the point q, so it is also the tangent space to the group orbit Orb(q). We have
Vq = Ker Tqn = {^(q) | £ G g} = Tq Orb(q), (5.1)
where g = TeG is the Lie algebra of G, and e is the identity element of G.
We select a fixed basis of g: e1 = 3/3§r, e2 = d/dx, e3 = d/dy, e4 = 3/30, e5 = 3/30 f. Then Tq Orb(q) is spanned by {(ea)Q(q) | a = 1, ..., 5}, where
3 3 3
(ei )*(?) = teM?) = ^ (e3)g(?) = ^
3 3 3 d (5.2)
{e^{q) = -VTx+Xd-y + W {e")Q{q) = Wf
According to (3.4), the subspaces Dq and Tq Orb(q) of TqQ have the following relationships:
i Dq + Tq Orb(q) = TqQ,
[Vq n Tq Orb(q) = Sq,
where
(5.3)
Sq = span <j — B\(q)^7 }> (5.4)
d u d
is a one-dimensional subspace of TqQ at q.
Since Sq C Tq Orb(q), a vector subspace gq is defined to be the set of Lie algebra elements in g whose infinitesimal generators, when evaluated at q, lie in Sq [46, 48], i.e.,
gq = {£ G g | Qq) G Sq}. (5.5)
Recall that q = (r, g) with r G R and g G G. Due to the invariance properties under the action of G, it can be proved that
g(r'fl) = Adfl g(r e). (5.6)
According to (5.4), it is clear that g(r'e) = span{e1(r)}, where
d d
Using the kinetic energy metric of the bicycle, e1 (r) can be augmented to form a moving orthogonal basis {ea(r) | a = 1, ..., 5} of the Lie algebra g, i.e.,
«(0,eb(r)), (0,ec(r))»re) = 0, if b = C, (5.8)
where ((•, -))q is the inner product on TqQ induced by the kinetic energy of the bicycle.
With the dimension condition in (5.3) and the invariance properties of the Lagrangian and constraint distribution under the action of G, we can define a principal connection A: TQ ^ g on the principal fiber bundle n: Q ^ Q/G = R. This principal connection, named the nonholo-nomic connection, is determined such that its horizontal space Hq is given by the orthogonal
complement to the space Sq within the space Dq [46, 48]. In the bundle coordinates q = (r, g), the nonholonomic connection takes the form of
A(vq) = Adfl ((g-1 )*g + Aloc(r)r), (5.9)
where vq = (r, g) e TqQ and Aloc(r)r = Aa(r)raea(r) is determined by the orthogonality condition. In the above, Adfl: g — g is the adjoint action of G on g and is invertible for each g e G.
Therefore, at the configuration q = (r,g), the admissible velocity q = (r, g) satisfying the nonholonomic constraints can be uniquely decomposed into two parts:
q = Qq (q)+ rh, (5.10)
where Q e gq, and rh = (r, -g*Aloc(r)r) is the horizontal lift of the tangent vector r on the base space R. According to (5.6), there is a unique vector Qloc e g(r' e) such that Q = Adfl Qloc. Substituting this relationship into (5.10), we obtain the following identity [46, 48]:
Qloc = Aloc(r)r + Z, (5.11)
where Z = (g-1 )*g. Let Qloc = Qaea(r). The nonholonomic constraint equations are equivalent to the following equations:
Qa = 0, a = 2,..., 5. (5.12)
In addition, by comparing (5.7) and (5.11), we find the following relationship between Q1 and the generalized velocities:
Q1 = ^ + ca(r)ra. (5.13)
Due to the invariance property of the Lagrangian, a reduced Lagrangian l can be defined in the reduced velocity phase space TQ/G = TR x g such that
l(r, r, Z) = L(r,e,r,Z). (5.14)
And a constrained reduced Lagrangian lc: D/G — R is defined as
lc(r, r, Qloc) = l(r, r, Qloc - Aloc(r)r), (5.15)
where the quotient manifold D/G is the reduced constraint space.
We define the component of the momentum map in the body representation,
Pi = ei(r))= m = In{r)Ql = /ii(r)^+<{r)fa' (5-16)
where I11(r) is the component of the local form of the locked inertia tensor [46, 48] and da(r) = = I11(r)ca(r). With these notations, we can employ the so-called reduced nonholonomic Lagrange-d'Alembert equations to obtain the reduced equations of the bicycle on D/G.3 One can refer to Theorem 5.7.3 in [46] and Theorem 7.5 in [48] to find the general forms of the
3Since the vector field X on D is G-invariant, it induces a well-defined vector field Y on D/G such that X and Y are -related, i. e., (nvo X = Y o nv, where : D — D/G is the projection.
reduced nonholonomic Lagrange-d'Alembert equations. For the bicycle, given lc as a function of (r, r, p1), these equations are specified as
dt\ dra J dra '
dt -lV \Pl) -r ^laVl' -r ^aßl
dPl = Cl,!11^)2 + V\aPl ra + Val31faf'3, (5.18)
where 111 (r) = (I11 (r)) 1 is the component of the inverse locked inertia tensor. And the other coefficients K^1, , Ka/3l, ^n, D{a, Da^1 come from the combinations of Aa(r), 111(r), the coefficients of l(r, r, Z), the derivatives of Aa(r), 111 (r), ea(r) with respect to r, and the structure constants C0bc(r) of the Lie algebra g defined by [ea(r), ec(r)] = C^bc(r)eb(r) (a, b, c = 1, ..., 5). Thus, these coefficients depend on the base variable r only. Their expressions can be found in [46, 48].
In particular, due to the antisymmetry of the structure constants C^bc(r) with regard to the two subscripts a and c, we must have Cj1(r) = 0. Therefore, the first term of the right hand side of (5.18) vanishes; namely, we have
^± = VloPlfa + Val3lfafl3. (5.19)
Substituting (5.16) into (5.17) and (5.19), we can obtain three second-order ODEs which take the same forms of the Gibbs-Appell equations (4.1). Both sets of equations (together with the velocity constraints (3.3)) determine the vector field X on the constraint distribution D. They must be equivalent. Since the coefficients I11 (r) and da(r) of (5.16) and those of (5.17) and (5.19) are dependent on r only, we know that the coefficients of the Gibbs-Appell equations (4.1) are the functions of d and S, and are independent of the other coordinates , x, y, . The first property stated in Theorem 1 thus follows.
In addition, with (5.16) substituted, (5.19) becomes
dar° + Ink + " 2>Un) ^ + (H " V{adl3 - rV3 = 0. (5.20)
We see that the term (0r)2 and the pseudoforce term vanish in this equation. This proves the second property stated in Theorem 1.
6. Relative equilibria
According to (5.16), Z = (r, a) forms another set of local coordinates of the reduced constraint space D/G. Thus, by transforming the Gibbs-Appell equations (4.1) into the first-order forms and removing the cyclic coordinate , we obtain the local representative of the reduced dynamic system on D/G,
Z = Y (Z), (6.1)
where Y = (Y1, Y2, Y3, Y4, Y5). The equilibrium points of the reduced system (6.1), called the relative equilibria, can be found by setting Y(Z0) = 0. Thus, they must take the form Zo = = (d0, S0, 0, 0, w0), where d0, S0, w0 are three constants determined by the following three equations:
U(e0, S0, W0) = 33(^0,S0)(^0)2 - m,S0)=0, i = 1, 2, 3. (6.2)
By Theorem 1, the third equation in (6.2) (i = 3) is always satisfied. So the three constants d0, , w0 need only satisfy the first two equations, whose independence is determined by the rank of the Jacobian matrix of f1, f2 with respect to d0, 50, w0,
d(h, /2)
'dfi df, df, dd0 c><50 dui0 df2 df2 df2
(6.3)
Theorem 2. Suppose that (0 = (90, S0, 0, 0, w0) is a relative equilibrium of the reduced system (6.1) and the corresponding Jacobian matrix d(f1, f2)/d(90 , 50, w0) has full rank. Then there locally exists a unique one-parameter solution family of relative equilibria that passes through Z0, i. e, (0 is not an isolated relative equilibrium.
Proof. Without loss of generality, we assume that the submatrix d(f1, f2)/d(d0, 50) has full rank. Then, according to the implicit function theorem, the solutions (00, ¿0, W0) of fi = = 0, i = 1, 2 can be written as single-valued functions of W0 in a neighborhood of Z0, i.e., we have d0 = d0(w0), ¿0 = ¿0(w0), where w0 G U(w0) and U(w0) is a neighborhood of w0, such that d0(w0) = d0, ¿0(w0) = ¿0 and fi(d0(W0), ¿0(w0), W0) = 0, i = 1, 2. This means that ((W0) = = (d0(W0), ¿0(W0), 0, 0, W0), w0 G U(w0), is a one-parameter solution family of relative equilibria. Theorem 2 thus follows. □
Note that the full rank condition of the Jacobian matrix d(f1, f2)/d(d0, 50, w0) is satisfied for almost all of triples (d0, 50, w0). Therefore, the relative equilibria of the system (6.1) form a collection of one-parameter solution families in D/G. These relative equilibria can be further divided into two types: one consisting of all the points (0, 0, 0, 0, w0), where w0 is an arbitrary constant, and the other consisting of the relative equilibria whose components d0 and ¿0 satisfy (d0)2 + (¿0)2 = 0. The first type of relative equilibria represents the uniform upright straight motion of the bicycle, while the second type represents the bicycle's uniform circular motion. To illustrate this point, we need to find the solution n(t) = (r(t), g(t), a(t)) of (4.2) whose projection on D/G is the relative equilibrium Z0, i.e., o n(t) = Z0. Clearly, we have r(t) = r0 = (d0, 50) and <r(t) = a0 = (0, 0, w0). The remaining component g(t) is determined by the following left invariant vector field on G:
g = g^ (6.4)
where £0 = w0e1(r0) (see (5.7)) is a constant element of the Lie algebra g.
When r0 = (0, 0), we have £0 = + + J> where £0 = Co = Rr^o, Co =
= Rrw0/Rf = Wf and Wf is the angular velocity of the front wheel. Let us identify each g = = (0r, (x, y, 0), ^f) G G with the following 7 x 7 block diagonal matrix:
g = diag(gi, g2, g3),
(6.5)
where
g1 =
cos — sin sin cos
^cos 0 — sin 0 x^
g2 =
sin 0 cos 0 y 0 0 1
9s =
cos — sin
/
sin
cos
(6.6)
As a result, the solution of (6.4) with the initial condition g(0) = (0r>0, (X0, y0, 00), 0) is
g(t) = g(0) exp(^01) = diag(gi(t), g2(t), g3(t)), (6.7)
d00 do0 a^Q
where
9 (t) = (C°s(Slt + o) — sin(Slt + o)N 1 Vsin(^ot + o) COs(Coi + o) ,
^cos 0o — sin 0o xo + (£o cos 0o)t^
sin 0o Cos 0o yo + (Co sin 0o)t 0 0 1
92 (t) =
/
93 (t) =
+ / o) — sin(Cot + / o)N vsin(Cot + / o) cos(Coi + o) ,
(6.8)
Therefore, n(t) = (r0, g(t), &0) represents the upright straight motion of the bicycle along the direction 0 = 0o and with a constant forward speed v0 = Çfi.
Similarly, when r0 / (0, 0), we know that £0 = ^ + + + ^ +
where £q = wo, £0 = and = 0. The solution of (6.4) with the initial condition g(0) = = (0r,o> (xo, Vo, A), f 0) is given by g(t) = diagg(t), £2(t), gs(t)), where
92 (t) =
/cos(eo4i + 0O) - Sin(e04i + + go(c°s(go^o)-c°s^)+g5(sin(g0t+^0)-sin^0)X
Sin(e04i + 0O) COS($i + 0O) yQ + go3(sin(g^+0n)-sin0n)-^(cos(gn^+0n)-cos0n)
v
0
0
1
(6.9)
/
and the expressions of g^t), g3(t) are the same as those in (6.8). This means that the trajectories of the coordinates x and y are
x(t) = xo —
Co cos 00 + Co sin 00 , V(Co2)2 + (Co3)2
A4 so
+
s 4 so
■cos(s 4t + 01),
y(t) = yo +
Co2 COS 00 ^O3 sin 00 , V(Co2)2 + (Co3)2
ù4 s0
(6.10)
+
s4 s0
■ sin(s41 + 0i),
where 01 is the initial phase angle determined by the following equations:
cos(01 — 0o) =
s3
vW+W
sin(01 — 0o) = —
s2
V(Co2)2 + (e03)2'
(6.11)
Clearly, n(t) = (r0 , g(t), a0) represents the circular motion of the bicycle with a radius sj(Co)2 + (£o3)7I£o4I and a constant angular velocity 0 = ¿¡q-
As an example, we adopt the benchmark parameters of the Whipple bicycle [11, 12]. Table 1 lists the values of these parameters, including the wheel base w, the trail c, the tilt angle A, the positions of the center of mass of the two frames xk, zk (k = s, h), the radii of the two wheels Rk (k = r, f), the masses of the four rigid bodies mk (k = r, s, f, h), and the components of the inertia tensors of the four rigid bodies IkxXX, Ikyy (k = r, s, f, h) and Ik>zz, IkxXZ (k = = s, h). By limiting u0 £ [0, +rc>), ¿0 £ [0, n] and solving (6.2), we obtain the relative equilibria of (6.1). Following Basu-Mandal et al. [12], we describe them by plotting d0 and ¿0, respectively, against a scaled front wheel speed arctan(RfUf/4) in Fig. 2. We see that there a total of five one-parameter solution families of the relative equilibria. Among them, the one family
/
Table 1. Values of the parameters of the benchmark bicycle: w = 1.02 m, c = 0.08 m, A =
ck zk Rk m
k 1k,xx Ik,yy 1k,zz
I
k,xz / / 2.8 -2.4
(a)
cS
1.5 1
0.5 0
-0.5 -1 -1.5
Body
Rear wheel (k = r) / / 0.3 2 0.0603 0.12
Saddle frame (k = s) 0.3 0.9 / 85 9.2 11 Front wheel (k = f) / / 0.35 3 0.1405 0.28 / /
Head frame (k = h) 0.9 0.7 / 4 0.0584 0.06 0.0076 —0.0091 Units: m for length, kg for mass, and kg • m2 for moment of inertia.
.-,-,-,-,-,-n (b)
-first family
second family ■■ third family - - — fourth family -fifth family
G
A, C
F ■K
-0—
E
B
H
3
2.5
2
■s
1.5
1
0.5
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 arctan(J2jW^/4)
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
arctan(_R^ctiy/4)
Fig. 2. Curves of (a) 00 and (b) S0 parameterized with arctan(R^Wf /4)
shown as the curve AB belongs to the first type, while the other four families shown as the curves CD, EFG, III and JKL belong to the second type. Since CD intersects with AB at the point D, we know that there are two solution families passing through D. According to Theorem 2, the Jacobian matrix d(f1, f2)/d(d0, ¿0, w0) at D must be singular. In addition, the two points F and K, in the third and fifth families, respectively, are cusp points, so their Jacobian matrices d(f1, f2)/d(d0, 50, w0) should also be singular.
7. Stability of relative equilibria
In this section, we will analyze the Lyapunov stability of the relative equilibria of (6.1). For a given relative equilibrium Z0 = (d0, S0, 0, 0, w0), we assume full rank of d(f1, f2)/d(d0, 50, w0) such that there is a one-parameter solution family passing through it. In addition, it is easy to prove that Z0 is a regular point of this family curve4, so we can define an arc length s (using the standard metric of R5) to parameterize the curve as Z(s) = ($0(s), ¿0(s), 0, 0, W0(s)) with Z(s0) = Z0. Differentiating Y(Z(s)) = 0 with respect to s at s0, we find that
= (7.1)
dZ \ ds J
4Noting the proof of Theorem 2, we have \\dZ(W0)/dw0|| = 0, so Z0 = Z(w0) is a regular point of the
family curve Z (W0).
This means that the Jacobian matrix J(Z0) = dY(Z0)/dZ of (6.1) at Z0 must have a zero eigenvalue whose eigenvector is parallel to the unit tangent vector dZ(s0)/ds of Z(s) at Z0. In other words, Z0 is not hyperbolic. In this case, the stability of Z0 may be analyzed with the help of the center manifold theorem [30]. Here, an alternative method considering the invariant manifold of (6.1) will be employed to present a more clear structure of the phase flow near Z0.
Due to the Lagrangian invariance and constraint distribution invariance, the total energy of the bicycle, which equals the sum of its kinetic energy and potential energy, induces a well-defined reduced constrained energy E on D/G [44-46, 48]. In the coordinates Z = (r, &), E can be written as follows:
m = ^9ij(r)&i&j + V(r). (7.2)
Since the total energy of the bicycle is conserved, we know that E is a first integral of (6.1), that is, QEi = {Z G D/G | E(Z) = E1} is an invariant manifold for each E1 g R.
Theorem 3. If the relative equilibrium Z0 = (90, §0, 0, 0, w0) satisfies
(1) W = 0,
(2) d(f, f2)/d(90, §0, W0) has full rank,
(3) the other four eigenvalues of J(Z0) have negative real parts,
then the system (6.1), restricted to the invariant manifold Qe0, where E0 = E(Z0), has an asymptotically stable equilibrium at Z0.
Proof. We divide the coordinates of Z into two parts, Z = (Z1, Z2), where Z1 = (r, r) and Z2 = <Pr. Thus, Z0 = (Z1 ,0, Z20), where Z1 ,0 = (90, §0, 0, 0) and Z2,0 = Similarly, the local representative of the vector field Y is divided into two parts, Y = (Y1, Y2), where Y1 = = (Y1, Y2, Y3, Y4) and Y2 = Y5. Let E(Z1, Z2, E1) = E(Z1, Z2)— E1 satisfying E(Z1 , 0, Z2 , 0, E0) = = 0. Since the matrix [gij] is positive definite, we have
dE(Z1,0, Z2,0, E0) fa x \ / n /T o\
-Qjr-= 9-M, ¿oH / 0. (7.3)
Therefore, according to the implicit function theorem, Z( can be written as a single-valued differentiate function of (Z1, E1) in a neighborhood of (Z1 ,0, Z2,0, E0):
Z( = F (Z1 ,E1), (7.4)
such that Z( , 0 = F (Z1 , 0, E0) and
E(Z1, F(Z1, E1), E1) = E(Z1, F(Z1, E1)) — E1 = 0. (7.5)
Differentiating (7.5) with respect to Z1 and E1, respectively, we have
dE dE OF _ dE OF
dc/ds;-1-
(7.6)
The second equation in (7.6) means that dF/dE1 = 0, so the neighborhood of Z0 in D/G can be represented as a union of invariant manifolds = {((1, Z2) £ D/G | Z2 = F(Z1, E1)}. In addition, the invariance of QE^ implies the following identity:
>2(ii, ^(Ci, E,)) = • n(Ci, ^(Ci, ^i)). (7.7)
Differentiating (7.7) with respect to Z1 and E1 at (Z1 0, E0), respectively, we obtain the following equations:
{ dY2(Co) + 2(Co) . ¿^(Ci.o, Eq) = 9F(C1|0, E0) /dY^p) + dY^Q) 9F(Ci,o, E0)
(7.10)
d<1 0^2 d(1 d(1 V dZ1 d(2 d(1
dY2((0) dF((10, E0) = dF((10, E0) (Cq) ^(Ci.o, E0)
9(-2 ' dE1 9Ci ' d(2 ' dE1 *
(7.8)
Here we have used the equation y1(C0) = 0. Since dF(Z1 0, E0)/dE1 = 0, the second equation in (7.8) becomes
dY2((0) dF((10, E0) (Cq) d<2 'X, d(2 • 1 •j
Using (7.9), the first equation in (7.8) can be simplified to
9Y2(Cq) = dF((h0, E0) (CQ) 9Ci 'X, ' 'A, '
In coordinates Z1, the restriction of the system (6.1) to QE^ can be written as follows:
Z1 = n(C1 ,F(C1, E0)). (7.11)
At the equilibrium Z1 = Z1 ,0, the Jacobian matrix of (7.11) is
Tr,, , ^i(Co) I ^i(Co) E0)
Let us define the following invertible matrix:
T = • (7-13)
J(Z0) can be transformed into the following form through the similarity transformation:
T-^(Co)T=(^^ ^ JT=( 0 *o } (7-14)
Since J(Z0) has four eigenvalues with negative real parts and one zero eigenvalue (conditions (2) and (3)), we know that the all four eigenvalues of J(Z\_ 0) must have negative real parts. As a result, Z1 0 is an asymptotically stable equilibrium of (7.11). □
Remark 1. Condition (1) in Theorem 3 is necessary. Otherwise, according to (6.2), we have dE(Zo)/dra = dV(v0)/dva = -Pa(r0) = 0, dE(Zo)/dái = gij(r0)&j = 0, i.e., Zo is a critical point of E. So in this case, QBg is degenerate and it cannot be written in the form of (7.4).
Remark 2. Note that the derivation of (7.14) needs neither the nonisolation property of Z0 (condition (2)) nor the prior assumption of the eigenvalues (condition (3)). As a corollary, for a general dynamic system, we suppose that it has a first integral as well as an equilibrium point. If the equilibrium point is not a critical point of the first integral, then the Jacobian matrix of the dynamic system at this equilibrium point must have a zero eigenvalue.
Theorem 4. Under conditions (1)-(3) of Theorem 3, the system (6.1), restricted to the invariant manifold Qe0, where E0 = E(Z0), has an isolated equilibrium at Z0. In addition, in a neighborhood of Z0, there is a one-to-one correspondence between the set of relative equilibria Z(s) and the set of invariant manifolds Qe .
Proof. According to Theorem 3, the restriction of (6.1) to QE^ has an asymptotically stable equilibrium at Z0 • This means that Z0 must be an isolated equilibrium on QE^.
Noting the similarity transformation (7.14), we know that the eigenvector corresponding to zero eigenvalue of J(Z0) (i.e., the kernel of J(Z0)) is parallel to the following vector:
ß = T[~J = ÔF(C /) ÖCVf(M |. (7.15)
Thus, we have (dZ(s0)/ds)' = k0/, where k0 = 0. The reduced constrained energy E along the family curve Z(s) is a function of s, i.e., E(s) = E(Z(s)). Using the first equation in (7.6), the derivative of E(s) with respect to s at s0 is
dE(sp) (dE(Z0) dE(Z0)\ ( -J-1«!,o) ' ^ \ _ dE(Z0)
(7.16)
According to the implicit function theorem, we confirm that s can be expressed as the single-valued function of E in a neighborhood of E0, i.e., s = s(E). Therefore, in this neighborhood,
there is a one-to-one correspondence between Z(s) and QE under the mapping s ^ E1 = E(s).
□
Based on Theorems 3 and 4, we obtain the following theorem.
Theorem 5. Under conditions (1)-(3) of Theorem 3, Z0 is a Lyapunov stable but not asymptotically stable equilibrium of the system (6.1). Furthermore, the solution of (6.1) under the initial condition Z(0) = Zo which, departs from. Z0 to a sufficiently small extent will converge to a nearby equilibrium ((s0), where s0 = s(E((0)).
Proof. Note that the neighborhood of Zo in D/G can be represented as a union of invariant manifolds QE^. On the one hand, the restriction of (6.1) to QE^ has an asymptotically stable equilibrium at Z0. On the other hand, by continuity the reduced system on each nearby invariant manifold QE^ has an asymptotically stable equilibrium at Z(s(E1)), which is located in a small neighborhood of Z0. This proves that, as an equilibrium of (6.1), Z0 is stable but not asymptotically stable. In addition, the solution ((t) of (6.1) under the initial condition ((0) = Zo satisfies Z(t) £ Q„,t ,, Vi ^ 0. Thus, if Zo is in a small neighborhood of Zo, we know that ((t)
will converge to ((s(E((0))) as t tends to +oo. □
As an example, adopting the benchmark parameters of the Whipple bicycle, we find that the relative equilibrium belonging to the first type (uniform upright straight motion) is Lyapunov stable if w0 G (14.307941787803680, 20.080873384627863) rad/s. As for the second type (uniform circular motion), only small sections of relative equilibria neighboring the cusp points F and K are Lyapunov stable [30].
8. Conclusions
In this paper, we analyze the symmetry of the bicycle system on a flat, level ground based on geometric mechanics. We find that the configuration space of the bicycle is a trivial principal fiber bundle whose structure group G plays the role of a symmetry group to keep the Lagrangian and constraint distribution invariant. This connects to the special structure of the bicycle dynamic system based on the Gibbs-Appell equations obtained in our previous work, that is, the coefficients and pseudoforce terms of these equations are dependent on the lean angle d and steer angle ô only, and the term quadratic in the rear wheel's angular velocity and the pseudoforce term always vanish in one of the equations.
We analyze the dimension relationship between the tangent space at any point and its subspaces, i.e., the admission velocity space and the tangent space to the group orbit. With this dimension condition and the invariance properties, a nonholonomic connection can be defined on the principal fiber bundle. As a result, we can employ the reduced nonholonomic Lagrange -d'Alembert equations to obtain the reduced equations of the bicycle. The aforementioned properties of the bicycle dynamic system can then be directly proved based on the structure of these reduced equations.
Furthermore, we point out that the Gibbs-Appell equations give the local representative of the reduced dynamic system on the reduced constraint space, whose relative equilibria, under the full rank condition of the Jacobian matrix d(f1, f2)/d(d0, ô0, w0), are not isolated but rather form several families of one-parameter solutions. These relative equilibria can be divided into two types, representing the bicycle's uniform upright straight and circular motions, respectively. Finally, we analyze the Lyapunov stability of the relative equilibria. Under certain conditions about the eigenvalues of the Jacobian matrix of the reduced system, a relative equilibrium can be stable but not asymptotically stable. However, the restriction of the dynamic system to an invariant manifold, which is the level set of the reduced constrained energy, may have an isolated asymptotically stable equilibrium.
In summary, we provide a comprehensive analysis of the symmetry and stability of the bicycle moving on a flat, level ground. Our results help to better understand bicycle dynamics and other nonholonomic systems that have invariance properties under the action of Lie groups.
References
[1] Carvallo, E., Théorie du mouvement du monocycle et de la bicyclette: Part 2. Théorie de la bicyclette, J. École polytechnique, Sér. 2, 1901, vol. 6, pp. 1-118.
[2] Whipple, F. J. W., The Stability of the Motion of a Bicycle, Q. J. Pure Appl. Math., 1899, vol. 30, no. 120, pp. 312-348.
[3] Boussinesq, J., Aperçu sur la théorie de la bicyclette, J. Math. Pures Appl., 1899, vol. 5, pp. 117-136.
[4] Klein, F. and Sommerfeld, A., Uber die Theorie des Kreisels, Leipzig: Teubner, 1898.
[5] Timoshenko, S. P. and Young, D.H., Advanced Dynamics, New York: McGraw-Hill, 1948.
[6] Sharp, R. S., The Stability and Control of Motorcycles, J. Mech. Eng. Sci, 1971, vol. 13, no. 5, pp. 316-329.
[7] Psiaki, M. L., Bicycle Stability: A Mathematical and Numerical Analysis, Thesis for Bachelor Degree, Princeton, N.J., Princeton Univ., 1979, 89 p.
[8] Neimark, Ju. I. and Fufaev, N.A., Dynamics of Nonholonomic Systems, Trans. Math. Monogr., vol. 33, Providence, R.I.: AMS, 1972.
[9] Hand, R. S., Comparisons and Stability Analysis of Linearized Equations of Motion for a Basic Bicycle Model, Master's Thesis, Ithaca, N.Y., Cornell Univ., 1988, 200 p.
[10] Papadopoulos, J.M., Bicycle Steering Dynamics and Self-Stability: A Summary Report on Work in Progress: Technical Report, Cornell Bicycle Research Project, Cornell University, Ithaca, NY, 1987.
[11] Meijaard, J. P., Papadopoulos, J. M., Ruina, A., and Schwab, A.L., Linearized Dynamics Equations for the Balance and Steer of a Bicycle: A Benchmark and Review, Proc. R. Soc. A Math. Phys. Eng. Sci., 2007, vol. 463, no. 2084, pp. 1955-1982.
[12] Basu-Mandal, P., Chatterjee, A., and Papadopoulos, J. M., Hands-Free Circular Motions of a Benchmark Bicycle, Proc. R. Soc. A Math. Phys. Eng. Sci., 2007, vol. 463, no. 2084, pp. 1983-2003.
[13] Kooijman, J.D.G., Meijaard, J. P., Papadopoulos, J.M., Ruina, A., and Schwab, A. L., A Bicycle Can Be Self-Stable without Gyroscopic or Caster Effects, Science, 2011, vol. 332, no. 6027, pp. 339342.
[14] Peterson, D. L., Gede, G., and Hubbard, M., Symbolic Linearization of Equations of Motion of Constrained Multibody Systems, Multibody Syst. Dyn, 2015, vol. 33, no. 2, pp. 143-161.
[15] Getz, N.H. and Marsden, J.E., Control for an Autonomous Bicycle, in Proc. of the IEEE Internat. Conf. on Robotics and Automation (Nagoya, May 1995), pp. 1397-1402.
[16] Koon, W. S. and Marsden, J.E., The Hamiltonian and Lagrangian Approaches to the Dynamics of Nonholonomic Systems, Rep. Math. Phys., 1997, vol. 40, no. 1, pp. 21-62.
[17] Boyer, F., Porez, M., and Mauny, J., Reduced Dynamics of the Non-Holonomic Whipple Bicycle, J. Nonlinear Sci., 2018, vol. 28, no. 3, pp. 943-983.
[18] Âstrom, K.J. and Murray, R. M., Feedback Systems: An Introduction for Scientists and Engineers, Princeton: Princeton Univ. Press, 2010.
[19] Baquero-Suârez, M., Cortés-Romero, J., Arcos-Legarda, J., and Coral-Enriquez, H., A Robust Two-Stage Active Disturbance Rejection Control for the Stabilization of a Riderless Bicycle, Multibody Syst. Dyn., 2019, vol. 45, no. 1, pp. 7-35.
[20] Moore, J. K., Human Control of a Bicycle, PhD Thesis, Los Angeles, Calif., Univ. of California, 2012, 316 p.
[21] Kooijman, J. D. G. and Schwab, A. L., A Review on Bicycle and Motorcycle Rider Control with a Perspective on Handling Qualities, Veh. Syst. Dyn., 2013, vol. 51, no. 11, pp. 1722-1764.
[22] Schwab, A. L. and Meijaard, J. P., A Review on Bicycle Dynamics and Rider Control, Veh. Syst. Dyn., 2013, vol. 51, no. 7, pp. 1059-1090.
[23] Sumbatov, A. S., Nonholonomic Systems, Regul. Chaotic Dyn., 2002, vol. 7, no. 2, pp. 221-238.
[24] Borisov, A. V. and Mamaev, I. S., Symmetries and Reduction in Nonholonomic Mechanics, Regul. Chaotic Dyn., 2015, vol. 20, no. 5, pp. 553-604.
[25] Chaplygin, S., On the Theory of Motion of Nonholonomic Systems. The Reducing-Multiplier Theorem, Regul. Chaotic Dyn., 2008, vol. 13, no. 4, pp. 369-376.
[26] Zhao, Z. and Liu, C., Contact Constraints and Dynamical Equations in Lagrangian Systems, Multi-body Syst. Dyn., 2016, vol. 38, no. 1, pp. 77-99.
[27] Peterson, D., Hubbard, M., and Estivalet, M., Analysis of the Holonomic Constraint in the Whipple Bicycle Model (P267), in The Engineering of Sport 7, M. Estivalet, P. Brisson (Eds.), Paris: Springer, 2008, pp. 623-631.
[28] Peterson, D. L., Bicycle Dynamics: Modelling and Experimental Validation, PhD Thesis, Los Angeles, Calif., Univ. of California, 2013.
[29] Wang, E. X., Zou, J., Xue, G., Liu, Y., Li, Y., and Fan, Q., Development of Efficient Nonlinear Benchmark Bicycle Dynamics for Control Applications, IEEE Trans. Intell. Transp. Syst., 2015, vol. 16, no. 4, pp. 2236-2246.
[30] Xiong, J., Wang, N., and Liu, C., Stability Analysis for the Whipple Bicycle Dynamics, Multibody Syst. Dyn, 2020, vol. 48, no. 3, pp. 311-335.
[31] Xiong, J., Wang, N., and Liu, C., Bicycle Dynamics and Its Circular Solution on a Revolution Surface, Acta Mech. Sin, 2020, vol. 36, no. 1, pp. 220-233.
[32] Jones, D.E., The Stability of the Bicycle, Phys. Today, 1970, vol. 23, no. 4, pp. 34-40.
[33] Kirillov, O. N., Locating the Sets of Exceptional Points in Dissipative Systems and the Self-stability of Bicycles, Entropy, 2018, vol. 20, no. 7, pp. 502-517.
[34] Kooijman, J. D. G., Schwab, A. L., and Meijaard, J. P., Experimental Validation of a Model of an Uncontrolled Bicycle, Multibody Syst. Dyn., 2008, vol. 19, nos. 1-2, pp. 115-132.
[35] Schwab, A. L., Meijaard, J. P., and Papadopoulos, J. M., Benchmark Results on the Linearized Equations of Motion of an Uncontrolled Bicycle, J. Mech. Sci. Technol., 2005, vol. 19, no. 1, pp. 292304.
[36] Limebeer, D. J. and Sharp, R. S., Bicycles, Motorcycles, and Models, IEEE Control Syst. Mag., 2006, vol. 26, no. 5, pp. 34-61.
[37] Sharp, R. S. and Limebeer, D. J., A Motorcycle Model for Stability and Control Analysis, Multibody Syst. Dyn., 2001, vol. 6, no. 2, pp. 123-142.
[38] Escalona, J. L. and Recuero, A. M., A Bicycle Model for Education in Multibody Dynamics and RealTime Interactive Simulation, Multibody Syst. Dyn., 2012, vol. 27, no. 3, pp. 383-402.
[39] Franke, G., Suhr, W., and Rieß, F., An Advanced Model of Bicycle Dynamics, Eur. J. Phys., 1990, vol. 11, no. 2, pp. 116-121.
[40] Lennartsson, A., Efficient Multibody Dynamics, PhD Thesis, Stockholm, KTH Royal Institute of Technology, 1999, 156 p.
[41] Äström, K.J., Klein, R.E., and Lennartsson, A., Bicycle Dynamics and Control: Adapted Bicycles for Education and Research, IEEE Control Syst. Mag., 2005, vol. 25, no. 4, pp. 26-47.
[42] Marsden, J.E., Lectures on Mechanics, Cambridge: Cambridge Univ. Press, 1992.
[43] Marsden, J. E., Ratiu, T. S., Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems, 2nd ed., Texts Appl. Math., vol. 17, New York: Springer, 2013.
[44] Zenkov, D. V., Integrability and Stability of Nonholonomic Systems, PhD Thesis, Columbus, Ohio, The Ohio State University, 1998.
[45] Zenkov, D. V., Bloch, A.M., and Marsden, J.E., The Energy-Momentum Method for the Stability of Non-holonomic Systems, Dynam. Stabil. Syst., 1998, vol. 13, no. 2, pp. 123-165.
[46] Bloch, A.M., Nonholonomic Mechanics and Control, J. Interdiscip. Math., vol. 24, New York: Springer, 2003.
[47] Kelly, S.D. and Murray, R.D., Geometric Phases and Robotic Locomotion, J. Robotic Systems, 1995, vol. 12, no. 6, pp. 417-431.
[48] Bloch, A. M., Krishnapasad, P. S., Marsden, J. E., and Murray, R. M., Nonholonomic Mechanical Systems with Symmetry, Arch. Ration. Mech. Anal., 1996, vol. 136, pp. 21-99.
[49] Ostrowski, J., Burdick, J., Lewis, A.D., and Murray, R.M., The Mechanics of Undulatory Locomotion: The Mixed Kinematic and Dynamic Case, in Proc. of the 1995 IEEE Internat. Conf. on Robotics and Automation (Nagoya, Japan, 21-27 May 1995): Vol. 2, pp. 1945-1951.
[50] Ostrowski, J. and Burdick, J., The Geometric Mechanics of Undulatory Robotic Locomotion, Int. J. Robot. Res., 1998, vol. 17, no. 7, pp. 683-701.
[51] Kang, H., Liu, C., and Jia, Y. B., Inverse Dynamics and Energy Optimal Trajectories for a Wheeled Mobile Robot, Int. J. Mech. Sci., 2017, vol. 134, pp. 576-588.
[52] Ostrowski, J., Lewis, A., Murray, R., and Burdick, J., Nonholonomic Mechanics and Locomotion: The Snakeboard Example, in Proc. of the IEEE Internat. Conf. on Robotics and Automation (San Diego, Calif, May 1994), pp. 2391-2397.
[53] Boyer, F. and Belkhiri, A., Reduced Locomotion Dynamics with Passive Internal DoFs: Application to Nonholonomic and Soft Robotics, IEEE Trans. Robot., 2014, vol. 30, no. 3, pp. 578-592.
[54] Shammas, E. A., Choset, H., and Rizzi, A. A., Geometric Motion Planning Analysis for Two Classes of Underactuated Mechanical Systems, Int. J. Robot. Res., 2007, vol. 26, no. 10, pp. 1043-1073.
[55] Gutman, E. and Or, Y., Symmetries and Gaits for Purcell's Three-Link Microswimmer Model, IEEE Trans. Robot, 2015, vol. 32, no. 1, pp. 53-69.
