CONTROLLED LAGRANGIANS AND STABILIZATION OF DISCRETE MECHANICAL SYSTEMS

Controlled Lagrangian and matching techniques are developed for the stabilization of relative equilibria of discrete mechanical systems with symmetry and equilibria of discrete mechanical systems with broken symmetry. Unexpected phenomena arise in the controlled Lagrangian approach in the discrete context that are not present in the continuous theory. In particular, to make the discrete theory effective, one can make an appropriate selection of momentum levels or, alternatively, introduce a new parameter into the controlled Lagrangian to complete the kinetic shaping procedure. New terms in the controlled shape equation that are necessary for potential shaping in the discrete setting are introduced. The theory is illustrated with the problem of stabilization of the cart-pendulum system on an incline, and the application of the theory to the construction of digital feedback controllers is also discussed.

In the controlled Lagrangian approach, one considers a mechanical system with an uncontrolled (free) Lagrangian equal to kinetic energy minus potential energy. To start with, one considers the case in which the Lagrangian is invariant with respect to the action of a Lie group G on the configuration space. To stabilize a relative equilibrium of interest, the kinetic energy is modified to produce a controlled Lagrangian which describes the dynamics of the controlled closed-loop system. The equations corresponding to this controlled Lagrangian are the closed-loop equations and the new terms appearing in those equations corresponding to the directly controlled variables correspond to control inputs. The modifications to the Lagrangian are chosen so that no new terms appear in the equations corresponding to the variables that are not directly controlled. This process of obtaining controlled Euler-Lagrange equations by modifying the original Lagrangian is referred to as kinetic matching.
One advantage of this approach is that once the form of the control law is derived using the controlled Lagrangian, the stability of a relative equilibrium of the closed-loop system can be determined by energy methods, using any available freedom in the choice of the parameters of the controlled Lagrangian. To obtain asymptotic stabilization, dissipation-emulating terms are added to the control input.
The method is extended in [9] to the class of Lagrangian mechanical systems with potential energy that may break symmetry, i.e., there is still a symmetry group G for the kinetic energy of the system but one may now have a potential energy that need not be G-invariant. Further, in order to define the controlled Lagrangian, a modification to the potential energy is introduced that also breaks symmetry in the group variables. After adding the dissipation-emulating terms to the control input, this procedure allows one to achieve complete statespace asymptotic stabilization of an equilibrium of interest.
The main objective of this paper is to develop the method of controlled Lagrangians for discrete mechanical systems. The discretization is done in the spirit of discrete variational mechanics, as in [20]. In particular, as the closed loop dynamics of a controlled Lagrangian system is itself Lagrangian, it is natural to adopt a variational discretization that exhibits good long-time numerical stability. This study is also motivated by the recent development of structure-preserving algorithms for the numerical simulation of discrete controlled systems, such as recent work on discrete optimization, such as in [11], [15], [16].
The matching procedure is carried out explicitly for discrete systems with one shape and one group degree of freedom to avoid technical issues and to concentrate on the new phenomena that emerge in the discrete setting that have not been observed in the continuous-time theory. In particular, it leads one to either carefully select the momentum levels or introduce a new term in the controlled Lagrangian to perform the discrete kinetic matching. Further, when the potential shaping is carried out, it is necessary to introduce non-conservative forcing in the shape equation associated with the controlled Lagrangian.
It is also shown that once energetically stabilized, the (relative) equilibria of interest can be asymptotically stabilized by adding dissipation emulating terms. The separation of controlled dissipation from physical dissipation remains an interesting topic for future research; even in the continuous theory there are interesting questions remaining, as discussed in [24].
The theoretical analysis is validated by simulating the discrete cart-pendulum system on an incline. When dissipation is added, the inverted pendulum configuration is seen to be asymptotically stabilized, as predicted.
The discrete controlled dynamics is used to construct a realtime model predictive controller with piecewise constant control inputs. This serves to illustrate how discrete mechanics can be naturally applied to yield digital controllers for mechanical systems.
The paper is organized as follows: In Sections II and III we review discrete mechanics and the method of controlled Lagrangians for stabilization of equilibria of mechanical systems. The discrete version of the potential shaping procedure and related stability analysis are discussed in Section IV. The theory is illustrated with the discrete cart-pendulum system in Section V. Simulations and the construction of the digital controller are presented in Sections VI and VII.
In a future publication we intend to treat discrete systems with nonabelian symmetries as well as systems with nonholonomic constraints.

II. AN OVERVIEW OF DISCRETE MECHANICS
A discrete analogue of Lagrangian mechanics can be obtained by considering a discretization of Hamilton's principle; this approach underlies the construction of variational integrators. See Marsden and West [20], and references therein, for a more detailed discussion of discrete mechanics.
Consider a Lagrangian mechanical system with configuration manifold Q and Lagrangian L : T Q → R. A key notion is that of a discrete Lagrangian, which is a map L d : Q×Q → R that approximates the action integral along an exact solution of the Euler-Lagrange equations joining the configurations q k , q k+1 ∈ Q, where C([0, h], Q) is the space of curves q : [0, h] → Q with q(0) = q k , q(h) = q k+1 , and ext denotes extremum.
In the discrete setting, the action integral of Lagrangian mechanics is replaced by an action sum where q k ∈ Q, k = 0, 1, . . . , N , is a finite sequence of points in the configuration space. The equations are obtained by the discrete Hamilton principle, which extremizes the discrete action given fixed endpoints q 0 and q N . Taking the extremum over q 1 , . . . , q N −1 gives the discrete Euler-Lagrange equations for k = 1, . . . , N − 1. This implicitly defines the update map Φ : Q × Q → Q × Q, where Φ(q k−1 , q k ) = (q k , q k+1 ) and Q × Q replaces the phase space T Q of Lagrangian mechanics.
Since we are concerned with control, we need to consider the effect of external forces on Lagrangian systems. In the context of discrete mechanics, this is addressed by introducing the discrete Lagrange-d'Alembert principle (see Kane, Marsden, Ortiz, and West [17]), which states that for all variations δq of q that vanish at the endpoints. Here, q denotes the vector of positions (q 0 , q 1 , . . . , q N ), and δq = (δq 0 , δq 1 , . . . , δq N ), where δq k ∈ T q k Q. The discrete oneform F d on Q × Q approximates the impulse integral between the points q k and q k+1 , just as the discrete Lagrangian L d approximates the action integral. We define the maps F d 1 , F d 2 : The discrete Lagrange-d'Alembert principle may then be rewritten as for all variations δq of q that vanish at the endpoints. This is equivalent to the forced discrete Euler-Lagrange equations

A. Controlled Euler-Lagrange Equations
This paper focuses on systems with one shape and one group degree of freedom. It is further assumed that the configuration space Q is the direct product of a one-dimensional shape space S and a one-dimensional Lie group G.
The configuration variables are written as q = (φ, s), with φ ∈ S, and s ∈ G. The velocity phase space, T Q, has coordinates (φ, s,φ,ṡ). The Lagrangian is the kinetic minus potential energy with G-invariant kinetic energy. The corresponding controlled Euler-Lagrange dynamics is where u is the control input.

B. Continuous-Time Kinetic Shaping
Assume that the potential energy is G-invariant, i.e., V (q) = V (φ), and that the relative equilibria φ = φ e ,ṡ = const are unstable and given by non-degenerate critical points of V (φ). To stabilize the relative equilibria φ = φ e ,ṡ = const with respect to φ, kinetic shaping is used. The controlled Lagrangian in this case is defined by where τ (φ) = κβ(φ). This velocity shift corresponds to a new choice of the horizontal space (see [8] for details). The dynamics is just the Euler-Lagrange dynamics for controlled Lagrangian (5), Lagrangian (5) satisfies the simplified matching conditions of [9] when the kinetic energy metric coefficient γ in (2) is constant. Setting u = −d γτ (φ)φ /dt defines the control input, makes equations (4) and (7) identical, and results in controlled momentum conservation by dynamics (3) and (4). Setting σ = −1/γκ makes equations (3) and (6) reduced on the controlled momentum level identical.
A very interesting feature of systems (3), (4) and (6), (7) is that the reduced dynamics are the same on all momentum levels, which follows from the independence of equations (3) and (6) of the group velocityṡ. We will see in Section IV that this property does not hold in the discrete setting, and one has to carefully select the momentum levels when performing discrete kinetic shaping.

C. Continuous-Time Potential Shaping
Now, consider the case when the kinetic energy is group invariant, but the potential is not. Consider the special case when the potential energy is having a local non-degenerate maximum at φ e , and the goal is to stabilize the equilibrium φ = φ e , s = s e . As it becomes necessary to shape the potential energy as well, the controlled Lagrangian is defined by the formula where V ε (y) is an arbitrary negative-definite function, and (φ e , s e ) is the equilibrium of interest. Below we assume that φ e = 0, which can be always accomplished by an appropriate choice of local coordinates for each (relative) equilibrium.

A. Discrete Controlled Dynamics
In discretizing the method of controlled Lagrangians, we combine formulae (1), (2), and (5). In the rest of this paper, we will adopt the notations This allows us to construct a second-order accurate discrete Lagrangian Thus, for a system with one shape and one group degree of freedom the discrete Lagrangian is given by the formula The discrete dynamics is governed by the equations where u k is the control input.

B. Kinetic Shaping
At first, it will be assumed that the potential energy is Ginvariant, i.e., V (q) = V (φ), and that relative equilibria φ k = 0, ∆s k = const of (11) and (12) in the absence of control input are unstable. We will see that one needs to either appropriately select the momentum levels or introduce a new parameter into the controlled Lagrangian to complete the matching procedure.
Motivated by the continuous-time matching procedure (see Section III), we define the discrete controlled Lagrangian by the formula where L τ,σ (q,q) is the continuous-time controlled Lagrangian. The dynamics associated with (13) is Equation (15) is equivalent to the discrete controlled momentum conservation: where Setting makes equations (12) and (15) identical and allows one to represent the discrete momentum equation (12) as the discrete momentum conservation law Theorem 1: The dynamics determined by equations (11) and (12) restricted to the momentum level p k = p is equivalent to the dynamics of equations (14) and (15) restricted to the momentum level p k = µ if and only if the matching conditions hold. Proof: Solve equations (16) and (19) for ∆s k and substitute the solutions in equations (11) and (14), respectively. This process is a simple version of discrete reduction [14]. A computation shows that the equations obtained this way are equivalent if and only if Since β(φ) = 0 and ∆φ k = 0 generically, equations (11) and (12) are equivalent if and only if which is equivalent to (20) Note that the momentum levels p and µ are not the same.
Remark. As h → 0, formulae (18) and (21) become respectively. That is, as h → 0, one recovers the continuoustime control input and the continuous-time matching condition, σ = −1/γκ. Condition µ = p/(1 + γκ) becomes redundant after taking the limit, i.e., the reduced dynamics can be matched on arbitrary momentum levels in the continuous-time case, which agrees with observations made in Section III. We now discuss an alternative matching procedure. Define the discrete controlled Lagrangian Λ d τ,σ,λ (q k , q k+1 ) by the formula The discrete controlled momentum is given by formula and equation (23) is equivalent to the discrete momentum conservation (19). Theorem 2: The dynamics (11) and (12) restricted to the momentum level p k = p is equivalent to the dynamics (22) and (23) hold.
Proof: Similar to the proof of Theorem 1, solve equation (19) for ∆s k and substitute the solution in equations (11) and (22), respectively. A computation shows that the equations obtained this way are equivalent if and only if which implies (25). Note that in this case we add an extra term to the controlled Lagrangian which eliminates the need for adjusting the momentum level.
Remark. The ratio Λ d τ,σ,λ /h becomes L d τ,σ + λκβ(ϕ)φ as h → 0. That is, as we let the time step go to 0, we obtain the continuous-time controlled Lagrangian modified by a term which is a derivative of the function λκ β(φ) dφ with respect to time. It is well-known that adding such a derivative term to a Lagrangian does not change the dynamics associated with this Lagrangian.
The stability properties of the relative equilibria φ k = 0, s k = const of equations (11) and (12) are now investigated.
Theorem 3: The relative equilibria φ k = 0, ∆s k = const of equations (11) and (12), with u k defined by (18), are spectrally stable if Proof: Let V ′′ (0) = −C, where C > 0 (see Section III). The linearization of the reduced dynamics (11) and (12) Observe that the value of p does not affect the linearized dynamics.
The linearized dynamics preserves the quadratic approximation of the discrete energy The equilibrium φ k = 0 of (28) is stable if and only if the function (29) is negative-definite at φ k = φ k+1 = 0. The latter requirement is equivalent to condition (27).
Remark. The stability condition (27) is identical to the stability condition of the continuous-time cart-pendulum system, and it can be rewritten as The spectrum of the linear map (φ k−1 , φ k ) → (φ k , φ k+1 ) defined by (28) belongs to the unit circle. Spectral stability in this situation is not sufficient to conclude nonlinear stability.
We now modify the control input (18) by adding the kinetic discrete dissipation-emulating term D(∆φ k−1 + ∆φ k ) 2h in order to achieve the asymptotic stabilization of the upward position of the pendulum. In the above, D is a positive constant. The discrete momentum conservation law becomes Straightforward calculation shows that the spectrum of the matrix of the linear map (φ k−1 , φ k ) → (φ k , φ k+1 ) defined by the reduced discrete dynamics belongs to the open unit disc. This implies that the equilibrium φ = 0 is asymptotically stable.

C. Potential Shaping
Recall that the the discrete dynamics associated with discrete Lagrangian (10) is governed by equations (11) and (12), where u k is the control input. The goal of the procedure developed in this section is to stabilize the equilibrium (φ, s) = (0, 0) of (11) and (12).
The dynamics associated with (30) is amended by the term w k in the discrete shape equation: This term w k is important for matching systems (11), (12) and (31), (32). The presence of the terms w k represents an interesting (but manageable) departure from the continuous theory. Let and where y k is obtained by substituting φ k and s k in formula (9). Remark. Equations (11), (12) define closed-loop dynamics when u k is given by formula (33). The terms w k vanish when β(φ) = const as they become proportional to the left-hand side of equation (32).
As in the case of kinetic shaping, the stability analysis is done by means of an analysis of the spectrum of the linearized discrete equations. We assume that the equilibrium to be stabilized is (φ k , s k ) = (0, 0).
The linearized dynamics preserves the quadratic approximation of the discrete energy E k,k+1 defined by where Since V ′′ 1 (0) is negative, the equilibrium (φ k , s k ) = (0, 0) of equations (35) and (36) is stable if the quadratic approximation of the discrete controlled energy (37) is negative-definite. The latter requirement is equivalent to conditions (34). The spectrum of the linearized discrete dynamics in this case belongs to the unit circle. Remarks. Spectral stability in this situation is not sufficient to conclude nonlinear stability. The stability conditions (34) are identical to the stability conditions of the corresponding continuous-time system.
Following [9], we now modify the control input (33) by adding the potential discrete dissipation-emulating term in order to achieve the asymptotic stabilization of the equilibrium (φ k , s k ) = (0, 0). In the above, D is a constant. The linearized discrete dynamics becomes where x k is obtained by substituting φ k and s k in formula (38).
where E k,k+1 is the quadratic approximation of the discrete energy (29). Recall that E k,k+1 is negative-definite (see the proof of Theorem 5). It is possible to show that, in some neighborhood of (φ k , s k ) = (0, 0), the quantity ∆x k−1 + ∆x k ≡ 0 along a solution of equations (40) and (41) unless this solution is the equilibrium (φ k , s k ) = (0, 0). Therefore, E k,k+1 increases along non-equilibrium solutions of (40) and (41). Since equations (40) and (41) are linear, this is only possible if the spectrum of (40) and (41) is inside the open unit disk, which implies asymptotic stability of the equilibrium of both linear system (40) and (41) and nonlinear system (11) and (12) with potenital discrete dissipation-emulating term (39) added to u k .

V. STABILIZATION OF THE DISCRETE PENDULUM ON THE CART
A basic example treated in earlier papers in the smooth setting is the pendulum on a cart. Let s denote the position of the cart on the s-axis, φ denote the angle of the pendulum with the upright vertical, and ψ denote the elevation angle of the incline, as in Figure 1. The configuration space for this system is Q = S × G = S 1 × R, with the first factor being the pendulum angle φ and the second factor being the cart position s. The symmetry group G of the kinetic energy of the pendulum-cart system is that of translation in the s variable, so G = R.
The length of the pendulum is l, the mass of the pendulum is m and that of the cart is M .
Since the Lagrangian for the cart-pendulum system is of the form (2), the discrete control laws (18) and (33) stabilize the upward vertical equilibrium of the pendulum. As in the continuous-time setting, the cart is stabilized by symmetrybreaking controller (33) and is not stabilized by symmetrypreserving controller (18).
Simulations of the discrete cart-pendulum system are shown in the next section.

VI. SIMULATIONS
Simulating the behavior of the discrete controlled Lagrangian system involves viewing equations (11) and (32) as an implict update map Φ : (q k−2 , q k−1 ) → (q k−1 , q k ). This presupposes that the initial conditions are given in the form (q 0 , q 1 ); however it is generally preferable to specify the initial conditions as (q 0 ,q 0 ). This is achieved by solving the boundary condition ∂L ∂q (q 0 ,q 0 ) + D 1 L d (q 0 , q 1 ) + F d 1 (q 0 , q 1 ) = 0 for q 1 . Once the initial conditions are expressed in the form (q 0 , q 1 ), the discrete evolution can be obtained using the implicit update map Φ.
We first consider the case of kinetic shaping on a level surface (with ψ = 0), when κ is twice the critical value, and without dissipation. Here, h = 0.05 sec, m = 0.14 kg, M = 0.44 kg, and l = 0.215 m. As shown in Figure 2, the φ dynamics is stabilized, but since there is no dissipation, the oscillations are sustained. The s dynamics exhibits both a drift and oscillations, as potential shaping is necessary to stabilize the translational dynamics.
When dissipation is added, the φ dynamics is asymptotically stabilized, as shown in Figure 3. However, even though the oscillations are damped, the s dynamics retains a drift motion, as expected.
We next consider the case of potential shaping on an inclined surface (with ψ = π 9 radians) without dissipation, with the other physical parameters as before. Here, our goal is to regulate the cart at s = 0 and the pendulum at φ = 0. We set V ε = − ε 2 y 2 . The control gains are chosen to be κ = 20, ρ = −0.02, and ε = 0.00001. It is worth noting that the discrete dynamics remain bounded near the desired equilibrium, and this behavior persists even for significantly longer simulation runs involving 10 6 time-steps. To more clearly visualize the dynamics, we only include a 4000 timestep segment of this computation in Figure 4. The exceptional stability of the discrete controlled trajectory can presumably be understood in terms of the bounded energy oscillations characteristic of symplectic and variational integrators.  Discrete controlled dynamics with kinetic shaping and without dissipation. The discrete controlled system stabilizes the φ motion about the equilibrium, but the s dynamics is not stabilized; since there is no dissipation, the oscillations are sustained. stabilizing control law, as illustrated in Figure 5. This is consistent with the stability analysis of Section V.

VII. MODEL PREDICTIVE CONTROLLER
We now explore the use of the forced discrete Euler-Lagrange equations as the model in a real-time model predictive controller, with piecewise constant control forces. Algorithm 1 below describes the details of the procedure.
The digital controller uses the position information it senses for t = −2h, −h to estimate the positions at t = 0, h during the time interval t = [−h, 0]. This allows it to compute a symmetric finite difference approximation to the continuous control force u(φ, s,φ,ṡ) at t = h/2 using the approximation where the overbar indicates that the position variable is being estimated by the numerical model. This control is then applied as a constant control input for the time interval [0, h]. This algorithm can be implemented in real-time if the two forward solves can be computed within the time interval h. On a 2.5 GHz PowerPC G5 running MATLAB, two forward solves take 631.2 µsec, which is sufficiently fast to drive a digital controller with a frequency in excess of 1.5 kHz. In our simulation, the digital controller had a frequency of 20 Hz, which involves a computational load that is easily accommodated by an embedded controller. The initialization of the discrete controller is somewhat involved, since the system is unforced during the time interval [0, 2h] while the controller senses the initial states, and computes the appropriate control forces. Consequently, a combination of the forced and unforced discrete Euler-Lagrange equations are used to predict the initial evolution of the system.
We present the numerical simulation results for the digital controller in both the case of kinetic shaping ( Figure 6) and potential shaping (Figure 7). We see that in the case of kinetic shaping, the system is asymptotically stabilized in only the φ variable, and the s dynamics exhibits a drift, whereas in the case of potential shaping, the system is asymptotically stabilized in both the φ and s variables. Notice that the use of a piecewise constant control introduces dissipation-like effects, which are reduced as the time-step is decreased.

VIII. CONCLUSIONS
In this paper we have introduced potential shaping techniques for discrete systems and have shown that these lead to an effective numerical implementation for stabilization in the case of the discrete cart-pendulum model. The method in this paper is related to other discrete methods in control that have a long history; recent papers that use discrete mechanics in the context of optimal control and celestial navigation are [11], [15], [16], and [23]. The method of discrete controlled Lagrangians for systems with higher-dimensional configuration space and with non-commutative symmetry will be developed in a forthcoming paper.

IX. ACKNOWLEDGMENTS
The research of AMB was supported by NSF grants DMS-0305837, DMS-0604307, and CMS-0408542. The research of ML was partially supported by NSF grant DMS-0504747 and a University of Michigan Rackham faculty grant. The research of JEM was partially supported by AFOSR Contract FA9550-05-1-0343. The research of DVZ was partially supported by NSF grants DMS-0306017 and DMS-0604108.