GEOMETRIC DISCRETIZATION OF NONHOLONOMIC SYSTEMS WITH SYMMETRIES

. The paper develops discretization schemes for mechanical systems for integration and optimization purposes through a discrete geometric approach. We focus on systems with symmetries, controllable shape (internal variables), and nonholonomic constraints. Motivated by the abundance of important models from science and engineering with such properties, we propose numerical methods speciﬁcally designed to account for their special geometric structure. At the core of the formulation lies a discrete variational principle that respects the structure of the state space and provides a framework for constructing accurate and numerically stable integrators. The dynamics of the systems we study is derived by vertical and horizontal splitting of the variational principle with respect to a nonholonomic connection that encodes the kinematic constraints and symmetries. We formulate a discrete analog of this principle by evaluating the Lagrangian and the connection at selected points along a discretized trajectory and derive discrete momentum equation and discrete reduced Euler-Lagrange equations resulting from the splitting of this principle. A family of nonholonomic integrators that are general, yet simple and easy to implement, are then obtained and applied to two examples-the steered robotic car and the snakeboard. Their numerical advantages are con-ﬁrmed through comparisons with standard methods.

1. Introduction.The goal of this paper is to develop integrators for mechanical systems subject to nonintegrable constraints on the velocities, i.e., nonholonomic constraints.We study systems that evolve on a configuration manifold Q = M × G constructed from a Lie group G whose action leaves the kinetic energy invariant (and so G is a group of symmetries) and a vector space M that describes the system internal shape.This general configuration space applies to systems from several domains, e.g., locomotion systems found in nature [23,28,22], vehicles used in robotics and aerospace [34,35,2,5], systems in molecular dynamics [20,37].
Their dynamics is derived by explicitly factoring out the group invariance through reduction by symmetry, and consequently splitting the equations of motion into vertical-corresponding to symmetries aligned with the constraints and defining the evolution of a momentum, and horizontal-defining the dynamics of the shape space.This has proven not only computationally beneficial, from reducing the dimension and avoiding numerical ill-conditioning, but also crucial in studying the stability, controllability, and motion generation of such systems.In this paper we focus on their proper discretization and propose geometric integrators that respect the statespace structure of the symmetries and constraints, preserve any invariants exhibited by the continuous system, and result in stable and accurate numerical schemes.
We follow the approach of discrete mechanics [30] and derive discrete equations of motion of the system through the discretization of the underlying variational principles governing the dynamics.In particular, we employ a Lagrange-d'Alembert-Pontryagin (LDAP) variational principle [39] that differs from a standard variational principle, such as Lagrange-d'Alembert's, by the presence of a new additional velocity variable v ∈ T q Q at each point q ∈ Q which by definition does not correspond to the rate of change of the configuration but this dependence is indirectly enforced using a kinematic constraint of the form q − v = 0 and a multiplier p ∈ T * q Q corresponding to the momentum.We formulate a discrete version of this principle by varying trajectories with discrete states of the form (q k , v k , p k ) ∈ T Q ⊕ T * Q and by evaluating the continuous Lagrangian and accounting for the constraint distribution along such a discrete path.While conceptually equivalent to using a discrete Lagrangian L d : Q × Q → R (which is the standard way to approximate the action integral in discrete mechanics, e.g. as formulated by Marsden and West [30]) and a discrete nonholonomic distribution D d → Q × Q (introduced in [10]) the Pontryagin formulation has two key practical advantages.The first property that motivated us to employ the approach is that it leads to a straightforward formulation of a reduced principle for nonholonomic systems that involve a nontrivial intersection of the tangent space of symmetries and the constraint distribution.This intersection contains a velocity component whose variation, in the classical variational formulation, must be restricted using the curvature of a nonholonomic connection (e.g.see [9]).Since the notion of discrete curvature of a connection in the discrete setting is still not well understood we believe it is more appropriate to relax such higher order constraints on variations and instead indirectly enforce them through a Pontryagin-type approach.As we will show such a formulation leads to a simple derivation of the discrete mechanics retaining the preservation properties of the dynamics.A second important advantage not explored in this paper lies in the ability to treat degenerate Lagrangian systems, as proposed in the continuous setting by [39], and their discrete reduction by symmetries.Further details about the principle are given in Sec.2.3, and Sec.3.3.
The results presented here build upon previous work on the variational discretization of systems with symmetries, as well as systems with nonholonomic constraints.Bobenko and Suris [3] and Marsden et al. [31] first studied the discrete Euler-Poincaré equations for systems on Lie groups; Bou-Rabee and Marsden [4] extended those ideas in the framework of the Hamilton-Pontryagin principle in order to design more versatile integrators; Jalnapurkar et al. [21] considered the discretization of the more general principle bundle case with abelian symmetry and its reduction using Routh's method.Nonholonomic constraints from a discrete variational viewpoint were first studied by Cortés and Martínez [10] who also considered the invariance of such systems with respect to Lie group actions and derived a momentum equation with properties consistent with the continuous case.M. de Leon et al. [12,13] considered an alternative discretization of nonholonomic systems in terms of generating functions and constructed constraint-preserving integrators.Fedorov and Zenkov [14] extended the reduced discrete approach to systems on a Lie group to include nonholonomic constraints on the group and derived the so called Euler-Poincaré-Suslov equations.McLachlan and Perlmutter [32] studied the general case of systems on vector spaces as well as on a group with nonholonomic constraints focusing on the time-reversibility and the importance of the preservation of invariants.
The framework of Lie groupoids offers a general viewpoint for studying the dynamics of nonholonomic systems and their geometric discretization.Iglesias et al. [18] developed nonholonomic integrators based on this methodology and examined their properties in terms of reversibility and momentum evolution.This approach is also suitable for the type of systems we consider, more specifically by considering an Atiyah Lie groupoid for discrete reduced systems such as the snakeboard.A family of geometric integrators on Lie groups that are not derived from a discrete variational principle were proposed by Ferraro et al. [15] that introduce an extra condition from an elastic impact onto the constraint distribution that, under certain conditions, has energy-preserving properties and can also be folded into a nonholonomic momentum evolution equation to produce an explicit integrator.The construction of these particular integrators is related to the idea of projecting the unconstrained discrete Euler-Lagrange equations onto the constraints to yield a nonholonomic integrator as discussed in [13].

Contributions.
Our work provides a systematic and practical approach to the design of structure-respecting integrators for nonholonomic mechanical systems.While related to the work of several authors noted above, our discrete approach to nonholonomic systems with symmetries captures the geometry of general systems (defined in terms of principle bundle and nonholonomic connections) with arbitrary group structure, constraints, and shape dynamics, and is not restricted to a configuration space that is either solely a group or has a Chaplygin-type symmetry.As a result our formulation contains a discrete momentum equation and a set of discrete reduced Euler-Lagrange equations analogous to the continuous case (e.g. as described in [1,9]) that explicitly account for and respect the interaction between symmetries and constraints in the vehicle dynamics.
The constructed integrators are applied to two examples: the steered car with simple dynamics and the snakeboard.Our method is compared to standard Runge Kutta methods in terms of its numerical accuracy, stability, and run-time efficiency.In addition, for the snakeboard we include comparisons with the discrete Lagrange-d'Alembert (DLA) integrator [10] and with the geometric integrator proposed by [15].Interestingly, all three nonholonomic integrators studied have similar accuracy for short term integrations but start to differ in stability as the time step and integration duration increase.

2.
Overview of mechanical integrators.The integrators employed in this paper are based on the discretization of geometric, variational principles.We start with a brief review of these variational integrators, as well as their extensions that handle group structure and symmetries.
2.1.Variational integrators.Variational integrators [30] are based on the idea that the update rules for a discrete mechanical system should come from a global "least action" principle such as Hamilton's principle.Variational integrators first approximate the time integral of the continuous Lagrangian L(q, q) by a function of two consecutive states q k and q k+1 : Equipped with such a discrete Lagrangian, one can now formulate a discrete version of Hamilton's principle according to where variations are taken with respect to each position q k along the path.Thus, if we use D i to denote the partial derivative w.r.t the i th variable, one must have: for every three consecutive positions q k−1 , q k , q k+1 of the mechanical system: this equation thus defines an integration scheme which solves for q k+1 knowing the two previous positions q k and q k−1 .
Simple example.Consider a continuous Lagrangian of the form L(q, q) = 1 2 qT M q − V (q) (here, V is a potential function) and define the discrete Lagrangian L d (q k , q k−1 ) = hL q k+ 1 2 , q k+1 − q k h , using the notation q k+ 1 2 := (q k + q k+1 )/2.The resulting equations are , which is a discrete analog of Newton's law M q = −∇V (q).For controlled (i.e., non conservative) systems, forces can be added using a discrete version of Lagranged'Alembert principle in a similar manner [30].
2.2.Lie group integrators.Classical integrators (including the variational ones we just reviewed) are formulated to compute a displacement in a vector space (e.g.R n ) added to the current configuration in order to advance the numerical solution in time.If the configuration space has special structure such as a Lie group or such that arise from other holonomic constraints then this solution must be projected onto the constraint manifold.This might be computationally expensive and causes energy dissipation which leads to inaccuracies that multiply over time.A typical example is the integration of rigid body dynamics either using rotation matrices whose projection requires costly orthogonal decomposition [17], or using quaternions that require increased resolution in order to perform stably even in a short-time integration (see [25] for a numerical comparison for systems on SE(3)).
In order to avoid such problems Lie group integrators have been proposed in the mechanical literature to automatically enforce that the updated poses remain within the proper group.These special integrators often express the updated configuration in terms of a retraction map τ , i.e., a map that expresses changes in the group in terms of elements in its Lie algebra.The exponential map is such a map proposed for integration purposes in [36].Retaining the Lie group structure and motion invariants under discretization has, since then, proven to be not only a nice mathematical property, but also key to improved numerics, as they capture the right dynamics (even in long-time integration) and exhibit increased accuracy [19,4].
More abstractly, Lie group integrators preserve symmetry and group structure for systems with motion invariants.Throughout this paper we will use a configuration manifold Q = M ×G where G is a Lie group (with Lie algebra g) whose action leaves the system invariant.In our case of vehicle dynamics, G = SE(3) is typically the group of rigid body motions of an articulated body while M is a space of internal variables of the vehicles.The idea of Lie group integrators is to transform the equations of motion from the original state space T Q into equations on the reduced space T M × g-elements of T G are translated to the origin and expressed in the algebra g.This reduced space being a linear space, standard integration methods can then be used.The inverse of this transformation is then used to map curves in the algebra back onto the group.Two standards retraction maps τ have been commonly used to achieve this transformation for any Lie group G: • Exponential map exp : g → G, defined by exp(ξ) = γ(1), with γ : R → G is the integral curve through the identity of the left invariant vector field associated with ξ ∈ g (hence, with γ(0) = ξ); • Canonical coordinates of the second kind ccsk : where {e i } is the Lie algebra basis.A third choice for τ , valid only for certain quadratic matrix groups [6] (which include the rigid motion groups SO(3), SE(2), and SE(3)), is the Cayley map cay : g → G, cay(ξ) = (e − ξ/2) −1 (e + ξ/2).Although this last map provides only an approximation to the integral curve defined by exp, we include it as one of our choices since it is very easy to compute and thus results in a more efficient implementation.

Unified view.
The integration algorithms proposed in this paper are based on a discrete version of the Lagrange-d'Alembert-Pontryagin (LDAP) principle [39,16].The LDAP viewpoint unifies the Lagrangian and Hamiltonian descriptions of mechanics [4] and extends to systems with symmetries and constraints.
The Lagrange-d'Alembert-Pontryagin principle.We briefly recall the general formulation of the continuous LDAP principle for a system with Lagrangian L : for variations that vanish at the endpoints.The curve v(t) describes the velocity determined from the dynamics of the system.In view of the formulation, v does not necessarily correspond to the rate of change of the configuration q.The additional variable p, though, indirectly enforces this dependence and corresponds to both Lagrange multipliers and the momenta of the system.Thus (1) generalizes the Lagrange-d'Alembert principle and is linked to the Pontryagin maximum principle of optimal control.
The LDAP principle is conceptually equivalent to the Lagrange-d'Alembert principle.For systems evolving on Lie groups our formulation could be alternatively derived along the lines of the discrete Euler-Poincaré (DEP) approach [31].For more general systems with principal bundle structure and regular Lagrangians it is possible to derive the discrete mechanics by restricting variations along sections consistent with the constraints and symmetries, as described in the groupoid framework [18].We choose to employ a discrete Pontryagin approach (Sec.4.2) approach since it is self-contained, does not require restrictions on the variations, and leads to simple construction of integrators that can be easily implemented.This is linked to the fact that the additional, seemingly redundant, multiplier variables in fact have concrete physical meaning -i.e. they denote the momenta -and fit naturally into the resulting integration algorithms.
3. Nonholonomic systems with symmetry.This section considers nonholonomic systems with symmetries in the continuous setting.We start by recalling standard concepts used to define the state space geometry.Then we formulate the nonholonomic LDAP principle and obtain the continuous nonholonomic equations of motion.
Assume that π : Q → Q/G is a principle bundle on the manifold Q with group G.The system has Lagrangian L : T Q → R and is subject to nonholonomic constraints defined by the regular distribution D consisting of subspaces D q ⊂ T q Q defined at each q ∈ Q.A group orbit is a submanifold denoted by Orb(q) := {gq | g ∈ G}.
If g is the Lie algebra of G, then T q Orb(q) = {ξ Q (q) | ξ ∈ g}, where ξ Q is the infinitesimal generator corresponding to the Lie algebra element ξ defined by Define the subspaces V q , S q , H q according to V q = T q Orb(q), S q = D q ∩ V q , D q = S q ⊕ H q .
These definitions have the following physical meaning (see [2] for a more detailed description): • V q -space of tangent vectors parallel to symmetry directions, i.e. the vertical space; • S q -space of symmetry directions that satisfy the constraints; • H q -space of tangent vectors that satisfy the constraints but are not aligned with any directions of symmetry, i.e. the horizontal space.We make the following additional assumptions that are standard in the literature (see [2,9]) • Dimension Assumption: For each q ∈ Q, we have Define the vector space s q ⊂ T q Q/G to be the set of Lie algebra element whose infinitesimal generators lie in S q , i.e. the space of symmetry directions that satisfy the constraints, by s q = {ξ(q) | ξ Q (q) ∈ S q }.The bundle with fibers s q at all q ∈ Q is denoted s.
Since our main interest is in a configuration space that is by construction of the form Q = M × G we will restrict any further derivations to the trivial bundle case.While the more general case (introduced so far) can be treated in an analogous manner with slight modification we stick to the trivial case for clarity without loosing the general applicability of our results.Using (global) trivial bundle coordinates (r, g) ∈ M × G we have ξ Q (r, g) = (0, ξ(r, g)g) ∈ S q .If we denote the basis for s q by {e b (r, g)}, for b = 1, ..., dim(S) then since D is G-invariant g can be factored out from this basis, i.e. e b (r, g) = Ad g e b (r), where {e b (r)} is the body-fixed basis.We denote s r the space spanned by {e b (r)}.
Lastly, the system is subject to control force f : [0, T ] → T * M restricted to the shape space.
3.1.Nonholonomic connection.With these definitions we can define a principle connection A : T Q → g with horizontal distribution that coincides with H q at point q.This connection is called the nonholonomic connection and is constructed according to A = A kin + A sym , where A kin is the kinematic connection enforcing the nonholonomic constraints and A sym is the mechanical connection corresponding to symmetries in the constrained directions (i.e. the group orbit directions satisfying the constraints).These maps satisfy where Ω ∈ s r is called the locked angular velocity, i.e. the velocity resulting from instantaneously locking the joints described by the variables r.Intuitively, when the joints stop moving the system continues its motion uniformly along a curve (with tangent vectors in S) with body-fixed velocity Ω and a corresponding spatial momentum that is conserved.By definition the principle connection can be expressed as where A(r) is the local form and the two components in (2) can be added to get The vector ver r q = (0, Ω) ∈ (T M × s) r is the vertical component relative the combined connection A and hor r q = ( ṙ, −A(r) ṙ) ∈ (T M × g) r is the horizontal component.Velocity vectors on T Q/G ∼ T M × g are split according to ( ṙ, g −1 ġ) r = ver r q + hor r q = (0, Ω) + ( ṙ, −A(r) ṙ).
3.2.Vertical and horizontal variations.We now consider variations of the configuration variables in the vertical and horizontal directions.Following [8,2] define the following Definition 3.1.Vertical variations δq = (δr, δg) are such that δr = 0 and δgg Since vertical and horizontal variations can be taken independently we can consider variational principles based on vertical and horizontal variations separately.The next section presents these principles and the derivation of the resulting nonholonomic equations of motion.
Lie algebra basis.In practice, the constrained symmetry space S (r,g) often must be generated from Lie algebra basis that depends on the shape r ∈ M .Therefore, we assume that the basis {e a (r) | a = 1, ..., dim(G)} spans g r , in such a way that {e b (r) | b = 1, ..., dim(s)} is an orthogonal basis for s r at each r.Then Ω ∈ s r in this basis is Ω = Ω b e b (r).
While both reduced Lagrangians capture the group invariance of the system, using the constrained reduced Lagrangian l c has several advantages.One is that, unlike ξ, the locked angular velocity Ω diagonalizes the kinetic energy which has important implications in studying the stability of the system [29].Another is that in the resulting equations of motion the rate of change of the generalized momentum decouples from that of the shape variables which is key in exploiting the holonomy of the system for locomotion and motion planning purposes.
Next, we begin from the general formulation of the LDAP principle (1) and extend it to the principle bundle setting introduced earlier.Then we formulate two equivalent reduced principles, first in terms of the reduced Lagrangian ℓ and then in terms of the constrained reduced Lagrangian l c .A related cotangent bundle reduction approach in the more general setting of Dirac structures (without explicitly focusing on nonhonolonmic constraints) is developed in [38].
The LDAP principle (1) can now be written in a reduced form by substituting the constrained reduced Lagrangian, by enforcing the nonholonomic connection, and by expressing the momentum components at point (r, g) , where µ ∈ g * is the body fixed momentum.Next, we formulate the reduced LDAP principle directly and state the resulting equations of motion.See [24] for details of the derivation.(3) The principle (3) now contains all information necessary to derive the equations of motion that explicitly account for the symmetries, nonholonomic constraints, and their interaction.The resulting equations of motion are for b = 1, ..., dim(s).These equations are standard in the literature on nonholonomic systems with symmetries (e.g.[1,26,9,2]) and were obtained here directly from a reduced variational principle by restricting the variations on the configuration variables only.This is equivalent to the approach of separately studying the evolution of a momentum map (e.g. as taken in [1]) or by additionally restricting the allowable variations on the velocity variables ξ or Ω (explored in [9,7]).Yet, our main motivation for this alternative intrinsic formulation is that such a self-contained principle can be easily cast in a discrete framework and we expect that the resulting discrete equations of motion would most closely preserve the variational, geometric structure of the original system.We develop the discrete framework in the next section.
4. Geometric discretization of nonholonomic systems.In this section we formulate a discrete variational principle and derive a family of simple nonholonomic integrators that account for the group structure and constraint distribution, respect the work-energy balance, and have discrete momentum equation and a corresponding momentum map with properties analogous to the continuous case.
4.1.Discrete approximation.As we noted in Sec. 2 the discrete mechanics approach is based on varying discrete trajectories in order to find critical values of an action integral approximated through quadrature.The approximation scheme can be simple, i.e. by joining the discrete points along the path with simple local interpolation and few quadrature evaluations along the segments; or they can be higher-order by further discretizing each segment and performing multiple quadrature computations.Here, for clarity we focus on a simple type of discretization using one quadrature point per segment (as depicted in Fig. 4.1); higher order discretization can also be achieved following the Lie group discretization proposed by [27] or [4].
Discrete trajectory.The discrete LDAP framework is defined in terms of a discrete trajectory whose states are elements of the tangent spaces in the reduced bundle of velocities and momenta.The trajectory is formally defined as follows (see also Fig  where ξ k = Ω k − A(r k+α )u k , with r k+α := (1 − α)r k + αr k+1 for a chosen α ∈ [0, 1] and the map τ : g → G represents the difference between two configurations in the group by an element in its algebra (see Sec. 2.2).
The discrete control force is f d : {t k } N k=0 → T * M approximating a force controlling the shape.
Constraint consistency.It is important to note that the discretization constraints between configurations and velocities from Dfn. 4.1 are invariant to left translations of the discrete trajectory.Left-translating a pair of configurations (g k , g k+1 ) used to define velocity ξ k is equivalent to applying the lifted left action to ξ k itself, i.e.
Therefore, the approximation (9) remains valid in a left trivialization e, g −1 (t k+α ) ġ(t k+α ) ≈ e, τ −1 (g −1 k g k+1 )/h , and the left-invariant discrete body-fixed velocity ξ k can be used for discrete reduction analogous to the continuous case.
Connection equivariance.Similarly to (9) the nonholonomic connection is approximated as for all k = 0, ..., N − 1 and α ∈ [0, 1].Note that the discretization constraint is also consistent with the required equivariance of the connection: In the above formulation variations δu k , δΩ k , δp k , δµ k are free.Allowing the Lagrangian and the connection to be evaluated at r k+α gives design freedom.Standard values of α are 0, 0.5, 1. Setting α = 0.5 provides more accurate approximation of the base dynamics while α = 0, 1 results in less accurate, but simpler integrators.As noted earlier, there are more general ways to define the discrete principle that allows arbitrary high approximation order but here the exposition is limited to lower order integrators for clarity.
The discrete momentum map.Next, we define a discrete momentum map, examine its properties, and compare it to its continuous analog.It is well-known that for the nonholonomic systems that we consider the momentum, even in the direction of constrained symmetries, is not conserved in general.Instead, the momentum equation defines how the momentum components evolve in time.In the discrete setting the vertical equation (or the discrete momentum equation) is its analog.Similar to the continuous setting the discrete vertical equation can be viewed as defining the evolution of a discrete momentum map that we define next.
and the spatial discrete momentum map J : With these definitions we can compute the evolution of the discrete momentum map along symmetry directions that are allowed by the constraints, i.e. along the elements of the basis {e b (r, g)} at point (r, g) ∈ Q, for b = 1, ..., dim(S).Note that this basis is constructed from a body-fixed basis {e b (r)} according to e b (r, g) = Ad g e b (r).For all such e b : Q → s we define the momentum map components Proposition 1. Discrete Momentum Map Change.The momentum components J nh b evolve along discrete LDAP solution trajectories according to Proof.Rewriting the momentum equation ( 21), derived form the discrete LDAP principle, in terms of the momentum map we obtain ), e b (r k ) = 0, which is a momentum map balance equation depicted in Fig. 4.2.1.In spatial frame it reads which yields the component difference.17) into (11) we get for k = 1, ..., N − 1.While in the continuous case the shape acceleration can be written independently from the momentum change evolution (see (4)) this is not generally the case in the geometric discretization resulting from (22).The reason is that, generally, momentum preservation results into an implicit numerical condition.For numerical purposes it is then easier to work with the horizontal equations expressed in terms of the reduced Lagrangian ℓ instead of the the constrained reduced Lagrangian l c .The equations are then The case of a linear connection.Next we study the special case when the connection A(r) is linear in the base point r.This case is useful in comparing the resulting integrator to the continuous case in order to gain insight into the effect of discretization.
Assume that A(r) is linear.The following expressions then trivially hold and after substituting them into (22) and using τ = exp the horizontal equations become where the curvature covariant derivative dA is defined by dA(u, δr) = DA(u, δr) − DA(δr, u), and B i are the Bernoulli numbers with the first few given by B 1 = −1/2, B 2 = 1/6, B 3 = 0, ...There are several special cases that lead to further simplification of the horizontal equations.The Lie bracket in (24) vanishes, for instance, when G is abelian; when g is one-dimensional; or whenever Ω and A(r)u lie in the same one-dimensional vector space for all Ω, r, and u (as in the snakeboard example from Sec. 5.2).Proposition 2. If the connection A(r) is linear and the ad operator in (24) maps to 0 along the path, then the non-conservative forces on the right-hand side of the reduced discrete Euler-Lagrange equations (22) match the continuous case exactly.
Proof.If the ad operator maps to 0 then the curvature equals the covariant derivative, i.e.B = dA.Then if we denote the continuous gyroscopic force by F A (r, u, µ) β = µ, dA(u, δr β ) , the discrete forces on the right-hand side of (24) become αF A (r k−1+α , u k−1 , µ k−1 ) + (1 − α)F A (r k+α , u k , µ k ) exactly representing the continuous force acting on the left and the right (depending on the value of α) of the fiber at r k .This claim is analogous to the result obtained by Cortés [11] for Chaplygintype symmetries and, as noted by the same author, if the gyroscopic forces vanish then the horizontal equations become a decoupled variational integrator on their own.Prop. 2 asserts that under similar conditions (linearity of the connection and vanishing of the bracket) this is also true for the systems considered in this paper.Numerical formulation.For numerical purposes it is convenient to write the discrete dynamics in terms of vector-matrix notation, by treating the Lie algebra variables ξ and Ω as vectors of coordinates with respect to a chosen canonical basis (see App. B for an example).The equations ( 23) and ( 21) are expressed as where ℓ k := ℓ(r k+α , u k , ξ k ) and ξ k = Ω k − A(r k+α )u k .These equations along with the reconstruction equations constitute the complete discrete evolution.5.1.Car with simple dynamics.We study the kinematic car model defined in [23] with added simple dynamics (Fig. 3).The configuration space is Q = S 1 × S 1 × SE(2) with coordinates q = (ψ, σ, θ, x, y), where (θ, x, y) are the orientation and position of the car, ψ is the rolling angle of the rear wheels, and σ is defined by σ = tan(φ) where φ is the steering angle.The car has mass m, rear wheel inertia I, rotational inertia K, and we assume that the steering inertia is negligible.The car is controlled by rear wheels torque f ψ and steering velocity u σ .The Lagrangian is then expressed as: and the constraints (see [23]) are cos θdx + sin θdy = ρdψ, − sin θdx + cos θdy = 0, where l is the distance between front and rear wheel axles, and ρ is the radius of the wheels.These constraints simply encode the fact that the car must turn in the direction in which the front wheels are pointing, that the car cannot slide sideways, and that the change in orientation is proportional to the steering angle and turning rate of the wheels.Note now that for any element g = (α, a, b) of SE(2), the action Φ g (q) = (φ ρ , φ L , θ + α, a + cos(α)x − sin(α)y, b + sin(α)x + cos(α)y) leaves the Lagrangian and constraints invariant.As the shape coordinates are r = (ψ, σ), the reduced Lagrangian thus becomes where ξ is used as a vector in R 3 of coordinates with respect to the standard Lie algebra basis (see App.B).
The matrix representation of the connection A dependent on r becomes: This model is an example of the principle kinematic case in which the constraint distribution complements the space tangent to the group orbits.This is easily seen noting that Thus, S = ∅ and there is no momentum equation.
Continuous equations of motion.The resulting continuous equations of motion are ẋ = ρ cos θ ψ, ẏ = ρ sin θ ψ, θ = ρσ l ψ, Car integrator.The discrete equations of motion will be derived by substituting the Lagrangian and the connection of the steered car into (25).Define u = (u ψ , u σ ) and ξ = −A(r) • u and pick τ = exp (defined in App.B).The discrete dynamics involves the term dexp −1 ξ defined in (28).Observing that in the case of the car ad * A(r)•u µ, A(r) • δ = 0 for any µ ∈ h * and u, δ ∈ T M and therefore (dexp −1 ξ ) * µ, A(r) • δ = µ, A(r) • δ which leads to the simplified equations of motion The exact equations of motion can now be derived by substituting µ = ∂ℓ ∂ξ = (Kξ 1 , mξ 2 , 0), ∂ℓ ∂r = (0, 0), ∂ℓ ∂u = (Iu ψ , 0), and become where v k = ρu ψ k , ω k = (ρ/l)σ k+α u ψ k .Thus, the integrator is easily computed as it is fully explicit for any choice of quadrature point α.We would like to refer the reader to the numerical comparisons with standard methods developed in [25] which verify that the proposed integrator is at least as accurate as the second order Runge-Kutta (RK2) at a fraction of its runtime.5.2.The snakeboard.The snakeboard (Fig. 3) is a wheeled board closely resembling the popular skateboard with the main difference that both the front and the rear wheels can be steered independently.This feature causes an interesting interplay between momentum conservation and the nonholonomic constraints: the rider is able build up velocity without pushing off the ground by transferring the momentum generated by a twist of the torso into motion of the board coupled with steering of the wheels through pivoting of the feet.When the steering wheels stop turning the systems moves along a circular arc and the momentum around the center of this rotation is conserved.A robotic version of the snakeboard also exists, equipped with a momentum-generating rotor and steering servos [33].
The shape space variables of the snakeboard are r = (ψ, φ) ∈ S × S denoting the rotor angle and the steering wheels angle, while its configuration is defined by (θ, x, y) denoting orientation and position of the board (see Figure 3).This corresponds to a configuration space Q = S ×S ×SE(2) with shape space M = S ×S and group G = SE(2).Additional parameters are its mass m, distance l from its center to the wheels, and moments of inertia I and J of the board and the steering.The kinematic constraints of the snakeboard are: − l cos φdθ − sin(θ + φ)dx + cos(θ + φ)dy = 0, l cos φdθ − sin(θ − φ)dx + cos(θ − φ)dy = 0, enforcing the fact that the system must move in the direction in which the wheels are pointing and spinning.The constraint distribution is spanned by three covectors: where a = −2l cos θ cos 2 φ, b = −2l sin θ cos 2 φ, c = sin 2φ.The group directions defining the vertical space are: and therefore the constrained symmetry space becomes: Since D q = S q ⊕ H q , we have H q = span ∂ ∂ψ , ∂ ∂φ .Finally, the Lagrangian of the system is L(q, q) = 1 2 qT M q where .
The reduced Lagrangian is, therefore: ℓ(r, u, ξ) = (u, ξ) T M (u, ξ).There is only one direction along which snakeboard motions lead to momentum conservation: it is defined by the basis element and, hence, there is only one momentum variable µ 1 = ∂ℓ ∂ξ , e 1 (r) .Using this variable we can derive the connection according to [33,1] as Continuous equations of motion.The dynamics of the system can be derived either in terms of Ω or in terms of µ as unknown variables.Here, we provide the resulting equations of motion based on µ since this has been the choice in previous work and will be easier to compare against.The continuous equations of motion can be derived as A nonholonomic integrator.The discrete equations of motion will be derived by substituting the Lagrangian and the connection of the snakeboard into (25) and choosing the map τ = exp.Observing that in the case of the snakeboard ad * ξ µ, η = 0 for any µ ∈ h * and ξ, η ∈ s and therefore for reduced d'Alembert-Pontryagin) and GNI perform as accurately as RK2, at a fraction of the computational time-due to their explicit update scheme.The accuracy of the DLA also almost coincides with RDP but the method is more expensive due to its implicit nature that requires iterative root-finding.Note that all methods are implemented using simplified and optimized C-code for maximum efficiency.In the second case of longer trajectories (Fig. 5) as the number of time steps used is reduced, the RK methods start to become unstable and accumulate large errors.The RDP integrator at this resolution now becomes more accurate than RK4.The implicit DLA scheme is not included in this graph since at such low resolutions it frequently failed to converge (using Euler-type initialization) either because of Jacobian ill-conditioning or after converging to a point that is not a root.This is due to the fact that when not initialized close to the real root the iterative Newton-type method can quickly become unstable.We plan to investigate these limitations of DLA in more depth in our future work.RDP and GNI still require only a fraction of the RK4 runtime with GNI being slightly more efficient.
6. Conclusion.This paper was concerned with the discretization of nonholonomic mechanical systems through a discrete variational approach.Our main contribution was the derivation of reduced integrators for systems with Lie group structure and internal controllable shape dynamics.The geometric nature of the integrators enabled us to identify a discrete nonholonomic momentum map and discrete constraint forces which in certain cases have properties similar to their continuous analogs.Numerical comparisons with standard integration methods as well as other nonholonomic integrators revealed that the variational and reduced nature of the proposed algorithms contributes to stable and accurate integration even at larger time steps.It would be useful to investigate the nature of these numerical results further through backward error analysis and to also establish a notion of optimality of the chosen discretization.Further insight is still necessary to precisely define the notions of a discrete curvature of a connection as well as discrete constraint forces.It is interesting to determine how the proposed integrators fit in the more general framework of Lie groupoids [18] and whether some of the raised issues can be explained through this more general viewpoint.

Figure 1 .
Figure 1.Discrete approximation (dashed) of continuous trajectories (solid) in the shape space (left) using linear interpolation, and in the group (right) using local geodesics defined by the flow of the map τ .The discrete velocity vectors shown approximate the average velocity along the segment and satisfy the constraint as defined underneath the figures.These velocity vectors are attached at quadrature points determined by the choice of α ∈ [0, 1].

Figure 2 . 1 k
Figure 2. Evolution of the discrete momentum map.At point r k−1 the map is computed by projecting the covector J loc k−1 onto s r k−1 defined by the basis e b (r k−1 ); then in the Lie algebra basis attached at r k the covector J loc k−1 transforms by Ad * g −1 k g k+1 and the change in the momentum map is computed by subtracting it from the next momentum J loc k and projecting onto s r k (the notation J loc k := J loc (r k , u k , ξ k ) was used with covectors drawn pointing towards the vectors that they act on).

Corollary 1 .
Properties of the momentum map 1.The map components J nh b are not conserved in general.2. If e b (r) are independent of r, then the discrete momentum equations are the discrete Euler-Poincaré equations projected onto the constraint symmetry space s. 3.If e b (r) are independent of r and if G is abelian then J nh b are constant along the discrete trajectory.This follows from the equality e(r k , g k ) = e(r k−1 , g k−1 ) since in this special case e b (r k ) = e b (r k+1 ) and Ad g = Id.Horizontal equations.Horizontal variations (Sec.3.2) are constrained according to g −1 k δg k = η k , where η k = −A(r k )δr k for variations δr k in the base.With this definition of η k and after substituting p k from (