Discrete Clebsch Optimal Control - Semantic Scholar
Discrete Clebsch Optimal Control Nikolaj Nordkvist, Peter E. Crouch, and Anthony M. Bloch
Abstract— In this paper we analyze a class of discrete optimal control problems. These systems are discretizations of a class of optimal control problems defined on invariant submanifolds which we denote embedded optimal control problems. We analyze a particular subset of these called discrete Clebsch optimal control problems where the invariant manifolds are group orbits. The generating Hamiltonian equations for such systems are analyzed. The analysis provides a large class of geometric integrators for mechanical systems. We apply the theory to two example systems: mechanical systems on matrix Lie groups and mechanical systems on the n-sphere.
I. I NTRODUCTION In this paper we analyze a class of discrete optimal control problems which are discretizations of a class of optimal control problems defined on invariant submanifolds. The latter class of problems we denote embedded optimal control problems. We analyze a particular subset of these called discrete Clebsch optimal control problems. For such systems the invariant manifolds are group orbits. The generating Hamiltonian equations for such systems are analyzed. We apply the theory to two example systems: mechanical systems on matrix Lie groups and mechanical systems on the n-sphere. This geometric approach provides globally defined integrators for a large class of mechanical systems on manifolds. Our theory provides an elegant and unified derivation of a large class of integrators for mechanical systems. This class includes many well known examples such as the MoserVeselov approach (see [1]) to integrating the motion of ndimensional rigid bodies and the discrete variational approach to deriving integrators surveyed in Marsden-West [2]. A globally defined discretization of the n-dimensional rigid body problem was derived and analyzed in [1] from a discrete variational analysis; the discrete equations of motion were analyzed further in [3]. A discrete analog of the EulerPoincaré equations was derived in [4], and closely related globally defined Lie group variational integrators for a rigid body in a conservative force field were studied in [5], [6], [7]. In [8] a similar globally defined discrete variational integrator for mechanical systems on S 2 × S 2 × · · · × S 2 was developed. A globally defined integrator does not depend on local coordinate systems for its implementation, and thus switching between coordinate systems is avoided. The Nikolaj Nordkvist, Mathematics and Sciences Division, Leeward Community College, Pearl City, HI 96782 [email protected] Peter E. Crouch, Electrical Engineering, University of Hawai‘i at M¯anoa, Honolulu, HI 96822 [email protected] Anthony M. Bloch, Department of Mathematics, University of Michigan, Ann Arbor, MI 48109 [email protected]. Supported in part by NSF grants DMS-0806756, DMS-0907949, and DMS-1207693
Crouch-Grossman methods, introduced in [9], provide a global framework for discrete analysis of rigid body equations but these methods fail to be symplectic. In [10] continuous Clebsch optimal control problems were defined and analyzed while [11] and [12] considered the generalization of embedded optimal control problems, the latter also analyzing their discrete definition. Examples of continuous embedded optimal control problems were examined in [13]. In [14] Clebsch integrators were proposed based on a discretization of the Clebsch action principle—as opposed to a discretization of the Clebsch optimal control problem which is the approach of this paper. The theory of Hamilton-Pontryagin integrators, presented in [15], provides a discretization based on a discrete Hamilton-Pontryagin variational principle. In particular in this paper we derive a simple algorithm for integrating geodesic flows on the n-sphere. Simulation are performed for the special case of S 2 which illustrates the geometric properties of the algorithm. II. D ISCRETE C LEBSCH OPTIMAL CONTROL PROBLEMS In this paper we will consider a special case of discrete dynamics of the form qk+1 = f (qk , uk ) where qk are the configuration variables and uk are the control variables. In [12] we analyze discrete optimal control problems with this dynamics. In particular we define discrete embedded optimal control problems. As discussed in [12] embedded optimal control problems are a generalization of optimal control problems on manifolds which have a natural interpretation on invariant submanifolds where the systems are locally controllable. A particular class of embedded optimal control problems are Clebsch optimal control problems (see [10] for the continuous definition, [12] for the discrete definition) where the invariant manifolds are group orbits. Consider the discrete dynamics given by qk+1 = ΦFk (qk ), where Φ : G × Q → Q, (F, q) 7→ ΦF (q), is the action of the (finite dimensional) Lie group G on a smooth manifold Q. We denote by J : T ∗ Q → g∗ the momentum map of Φ. We have that a group orbit orb(q0 ), q0 ∈ Q, by construction is an invariant set under the discrete dynamics. We denote ¯ : G × orb(q0 ) → orb(q0 ) the induced action on by Φ ¯ : T ∗ orb(q0 ) → g∗ its momentum map. orb(q0 ) and by J In order to be able to apply the results obtained in [12] for discrete embedded/associated optimal control problems we furthermore need to assume that G and orb(q0 ) are embedded submanifolds of some vector spaces and that these embeddings are obtained as the preimages at 0 of some smooth submersions (mapping between vector spaces).
We consider the following discrete version of the Clebsch optimal control problem: Problem A (Discrete Clebsch Optimal control problem): Minimize N −1 X
ℓd (qk , Fk )
k=0
over all sequences {qk }N k=0 , qk ∈ Q, satisfying qk+1 = ΦFk (qk ) with Fk ∈ G and fixed endpoints q0 , qN ∈ Q. For the discrete Clebsch optimal control problem to be well defined we need qN ∈ orb(q0 ). The associated optimal control problem is then the following:
Proof: We will only show this for (2). The proof is the same for (4). From Proposition 6.3.2 in [17] we have that ∗ ∗ T ∗ ΦF −1 = T ∗ Φ−1 Fk : T Q → T Q is a symplectic map. The k d map (qk , pk ) 7→ (qk , pk + λ0 δℓ δq (qk , uk )) is also symplectic since this is the time one flow of the Hamiltonian vector field on T ∗ Q with Hamiltonian H(q, p) = −λ0 ℓd (q, uk ). Since (2) is the composition of these two symplectic maps it itself must be symplectic. We define Mk as δℓd Mk := J(αk ) + J (qk , Fk ) , δq
for αk ∈ T ∗ Q. By construction we have Mk ∈ g⊥ qk = {ξ ∈ Problem B (Discrete Clebsch assoc. opt. control problem): g∗ | hξ, ηi = 0 ∀ η ∈ gq } since J(Tq∗ Q) ⊂ g⊥ . q If αk is a Minimize solution to the discrete symmetric representation (1)-(2) then N −1 we get for a left action that Mk satisfies the equation X ℓd (qk , Fk ) δℓd (q , F ) M =J(q , p ) + J k=0 k+1 k+1 k+1 k+1 k+1 δq N over all sequences {qk }k=0 , qk ∈ orb(q0 ), satisfying qk+1 = δℓd ∗ −1 p + =J Φ (q ), T Φ (q , F ) ¯ k Fk k k k ΦFk (qk ) Fk ΦFk (qk ) with Fk ∈ G and fixed endpoints q0 ∈ Q, qN ∈ δq orb(q0 ). δℓd (qk+1 , Fk+1 ) +J From [12] we have the following results regarding the δq extremals of the discrete Clebsch optimal control problems: δℓd ∗ =AdF −1 J qk , pk + (qk , Fk ) k δq Theorem 1: Assume that Φ is a left action. A normal δℓd extremal for the discrete Clebsch optimal control problem (qk+1 , Fk+1 ) . +J as given by applying the discrete maximum principle (see δq [16]) is a solution of These same calculations could be carried out for (3). We can δℓ δℓd therefore express Theorem 1 in terms of qk , Fk , and Mk as d ∗ L Fk (qk , Fk ) = Tid (qk , Fk ) , (1) J(αk ) + J follows: δq δF δℓd Theorem 2: Assume that Φ is a left action. A normal . (2) αk+1 = T ∗ ΦF −1 αk + verαk extremal for the discrete Clebsch optimal control problem k δq as given by applying the discrete maximum principle is a for αk ∈ T ∗ Q. The pullback bundle i∗orb(q0 ) (T ∗ Q) is solution of invariant under the flow of (2). δℓd ∗ An extremal for the discrete Clebsch associated optimal Mk = Tid L Fk (qk , Fk ) , (5) δF control problem is a solution of ¯ qk+1 = ΦFk (qk ), (6) δℓd δℓ ¯ k) + J ¯ δℓd (qk , Fk ) = T ∗ LF J(β (q , F ) , (3) d k k ∗ id k ¯ (qk+1 , Fk+1 ) , (7) Mk+1 = AdF −1 (Mk ) + J δF δq k δq ¯ δℓ d ¯ −1 βk + ver ¯ βk ¯ . (4) for qk ∈ Q and Mk ∈ g⊥ . βk+1 = T ∗ Φ Fk qk δq An extremal (must be normal) for the discrete Clebsch for βk ∈ T ∗ orb(q0 ). The discrete Clebsch associated optimal associated discrete optimal control problem is a solution of control problem does not admit abnormal extremals. δℓd ∗ If αk ∈ i∗orb(q0 ) (T ∗ Q) is a solution to (1)-(2) then βk = Mk = Tid LFk (qk , Fk ) , (8) δF T ∗ iorb(q0 ) αk ∈ T ∗ orb(q0 ) is a solution to (3)-(4). (9) qk+1 = ΦFk (qk ), δℓ For a right action we get very similar results, the only ¯ d ¯ ∗ Mk+1 = Ad∗F −1 (Mk ) + J (10) difference being that Tid LFk in (1) and (3) is replaced with ¯ (qk+1 , Fk+1 ) , k δq ∗ Tid RFk . We implicitly assume that (1) can be inverted to ⊥ give Fk as a function of αk . We refer to (1)-(2) as the for qk ∈ orb(q 0 ) and Mk ∈ gqk . discrete symmetric representation of the discrete extremal If qk , Mk is a solution to (8)-(10) then qk , Mk is a generating equations. solution to (5)-(10). Proposition 1: The discrete maps (2) and (4) are symplecFor a right action we would get very similar equations, the ∗ tic. only difference being that in (5) and (8) Tid LFk is replaced
∗ with Tid RFk and in (7) and (10) Ad∗F −1 (Mk ) is replaced k with Ad∗Fk . We refer to (5)-(7) as the discrete EulerPoincaré representation of the discrete extremal generating equations. In this theorem it is implicitly assumed that (5) can be solved to give Fk as a function of qk and Mk .
III. E XAMPLES We now consider two basic applications which illustrate the general theory. A. Discrete Clebsch optimal control problems on matrix Lie groups Consider Q = gl(n) with the right action Φ of a matrix Lie group G as defined by right matrix multiplication, that is Φg (q) = qg for q ∈ gl(n) and g ∈ G. The orbit that arises when q0 ∈ G is the matrix Lie group G itself. Given a self-adjoint positive definite (with respect to the trace inner product) operator Σ : g → g we let Σd : gl(n) → gl(n) denote a self-adjoint operator that satisfies P ◦ Σd |g = Σ for some positive definite operator Σ : g → g∗ ≃ g; this is equivalent to hu, Σd (u)i = hu, Σ(u)i for all u ∈ g. Such an extension can always be defined as the choice Σd = Σ ◦ P demonstrates. Let us consider the following discrete cost function 11 h2
ℓd (qk , Fk ) =
h(Fk − I), Σd (Fk − I)i − hV (qk ), (11)
for qk ∈ gl(n) and Fk ∈ G. In this problem Fk can be interpreted as an approximation to exp(huk ) where uk ≈ q −1 (tk )q(t ˙ k ) ∈ g is the discrete control at time tk and h = tk+1 − tk is a (fixed) step size satisfying T = N h. This discrete cost function therefore approximates Z tk+1 ℓ q(t), u(t) dt, (12) ℓd (qk , Fk ) ≈ tk
where ℓ(q, u) = 21 hu, Σ(u)i−V (q). When modeling the discrete dynamics of a mechanical system the term 21 hu, Σ(u)i will correspond to the kinetic energy and V (q) will correspond to the potential energy. From straightforward calculations, which can be found in detail in [12], we get that T ∗ ΦF −1 (q, p) = (qg, pg −T ), k
∗ Tid R Fk
J(q, p) = P(q T p), δℓd (qk , Fk ) = h1 P Σd (Fk − I)FkT , δF
where P is the orthogonal projection P : gl(n) → g∗ ≡ g with respect to the trace inner product. Thus from Theorem 1 we get that the discrete symmetric representation for this system becomes T 2 T ∂V (qk ) = P Σd (Fk − I)FkT , (13) hP(qk pk )−h P qk ∂q qk+1 = qk Fk , (14) ∂V pk+1 = qk − h (qk ) Fk−T . (15) ∂q
From Theorem 2 we get that the discrete Euler-Poincaré representation of this system is hMk = P Σd (Fk − I)FkT , qk+1 = qk Fk ,
Mk+1 =
Ad∗Fk (Mk )
− hP
∂V T qk+1
∂q
(qk+1 ) .
For F ∈ G, M ∈ g∗ = g, and η ∈ g we get
hAd∗F (M ), ηi = hM, AdF (η)i = M, F ηF −1
= F T M F −T , η
= P(F T M F −T ), η ,
and thus we have the identity
Ad∗F (M ) = P(F T M F −T ). Using this the discrete Euler-Poincaré representation for this system can be written as (16) hMk = P Σd (Fk − I)FkT , qk+1 = qk Fk , (17) ∂V T (qk+1 ) . (18) Mk+1 = P(FkT Mk Fk−T ) − hP qk+1 ∂q In [12] we show that the Euler-Poincaré representation (16)(18) is at least a first order approximation of the flow of the continuous Euler-Poincaré representation with cost function ℓ(q, u). Let us consider the case where V = 0 and Φ is the action ˜ is the of Gl(n) on Q given by right multiplication and Φ action of Gl(n) on Q given by left multiplication. Then we ˜ g˜ and since ℓd when V = 0 have that Φg commutes with Φ is independent of qk we get that, as in [12], ˜ k , pk ) = pk q T , J(q k is conserved along the flow of the discrete symmetric rep¯˜ , p ) = P(p q T ) is resentation (13)-(15). Therefore J(q k k k k conserved along the flow of (3)-(4). Since Mk = J(qk , pk ) when V = 0 we get that P(pk qkT ) = P(qk−T qkT pk qkT ) = P qk−T P(qkT pk )qkT = P((qk−T Mk qkT ) = Ad∗q−1 (Mk ) and k therefore Ad∗q−1 (Mk ) is conserved along the flow of the disk crete Euler-Poincaré representation (16)-(18). These results regarding conserved quantities were derived directly, using the discrete extremal generating equations only, in [18]. We now consider the special case of the orthgonal group which of course is of particular interest in applications. Example 1: Discrete generalized rigid body equations. Consider the Lie group G = SO(n) with Lie algebra g = so(n) and take Rk ∈ SO(n) and Ω ∈ so(n) to be group and algebra elements. First let us consider the free generalized rigid body, i.e., the case when V = 0. Let Σ : so(n) → so(n) be given by Σ(Ω) = Jd Ω + ΩJd , where Jd is a symmetric n×n matrix which makes Σ positive definite. The orthogonal projection Pso(n) : gl(n) → so(n) is given by Pso(n) (A) = 1 T 2 (A − A ). We choose Σd (A) = 2Jd A,
which gives Pso(n) ◦ Σd (Ω) = 12 (2Jd Ω − (2Jd Ω)T ) = Σ(Ω), showing that this choice of Σd indeed satisfies the condition Pso(n) ◦ Σd |g = Σ. Since Pso(n) (2J(Fk − I)FkT ) = Pso(n) (2J − 2JFkT ) = Fk J − JFkT , the discrete Euler-Poincaré representation (16)-(18) for this system then becomes hMk = Fk Jd − Jd FkT ,
B. Discrete Clebsch optimal control problems on the nsphere Consider Q = Rn+1 with the group G = SO(n+1) acting on Q by left matrix multiplication. Then S n is a group orbit. We define the discrete cost function as ℓd (qk , Fk ) =
T ∗ Φg−1 (q, p) = (gq, g T p),
Mk+1 = FkT Mk Fk .
(Jω)× = ω × J − J ω × , where J = 12 tr(J)I − J is a modified inertia matrix, we have Jd = 2J . We define the discrete angular velocity (in the body fixed frame) as for the continuous model Mk =: Σ(ωk× ) = 12 (Jωk )× . Since (F η)× = F η × F T , for F ∈ SO(3) and η ∈ R3 , we get ∂V (Jωk+1 )× = 2Mk+1 = 2FkT Mk Fk − 2hP (Rk+1 ) ∂R T × × = (Fk ωk ) − τ (Rk+1 ), and therefore the Euler-Poincaré representation (16)-(18) for the rigid body translates to h(Jωk )× = Fk J − J FkT , Rk+1 = Rk Fk , Jωk+1 = FkT Jωk − τ (Rk+1 ). These are the discrete rigid body equations analyzed in [5] in their Hamiltonian form. The discrete rigid body equations in SE(3), as studied in [7], reduce to these equations when there is no translational velocity. In [7] it was shown that this discrete map is a useful integrator as it provides (at least) a first order approximation to the flow of the continuous rigid body equations. The first of these equations is an implicit equation in the unknown Fk ∈ SO(3). The numerical solution of this implicit equation has been discussed in [5]. In [3] explicit solutions to this equation were derived in terms of solutions to an algebraic Riccati equation. The remaining equations are explicit and thus are straightforward to implement numerically.
hFk − I, Fk − Ii − h2V (qk ).
In [12] we show that this discrete cost function approximates the cost of a continuous optimal control problem with solutions being flows of the mechanical system corresponding to a point mass on S n with potential energy V . From [12] we have that
Rk+1 = Rk Fk ,
These are the discrete rigid body equations studied in [1]. Now let us consider a true rigid body, which is the special case when G = SO(3) and V is the potential energy which may be nonzero. The inertia matrix J is a 3×3 matrix which is positive definite and symmetric. Let (·)× : R3 → so(3) be defined by a× b = a × b for a, b ∈ R3 with × being the cross product. As in [7] we define Σ(ω × ) = 12 (Jω)× , which is positive definite by construction, and since we have the identity
1 h
∗ Tid LFk
J(q, p) = Pso(n+1) (pq T ), δℓd (qk , Fk ) = h1 Pso(n+1) (I − FkT ) δF = h1 (Fk − FkT ).
Therefore we get from Theorem 1 that the discrete symmetric representation for this system is T ∂V T 1 p − h2 (F − F ) = hP (q ) qk , (19) k k k so(n+1) k 2 ∂q qk+1 = Fk qk , (20) ∂V pk+1 = Fk pk − h2Fk (qk ), (21) ∂q for qk , pk ∈ Rn+1 . Using Fk = exp(fk ), fk ∈ so(n + 1), we have T 1 1 2 (Fk − Fk ) = 2 exp(fk ) − exp(−fk ) = sinh(fk ).
Thus the discrete symmetric representation for this system formally becomes T ∂V −1 hPso(n+1) pk − h2 Fk = exp sinh , (qk ) qk ∂q qk+1 = Fk qk , ∂V pk+1 = Fk pk − h2Fk (qk ), ∂q Since Ad∗F −1 (Mk ) = FkT Mk Fk = k
=
1 2h (Fk
1 2h Fk (Fk
− FkT )Fk
− FkT ) = Mk ,
we get from Theorem 2 that the discrete Euler-Poincaré representation for this system is hMk = 12 (Fk − FkT ), qk+1 = Fk qk , Mk+1 = Mk − h2Pso(n+1)
(22) (23)
∂V T (qk+1 )qk+1 , ∂q
where M0 ∈ g⊥ q0 . This can formally be rewritten as Fk = exp sinh−1 (hMk ) , qk+1 = Fk qk ,
Mk+1 = Mk − h2Pso(n+1)
∂V T (qk+1 )qk+1 . ∂q
(24)
In [12] we show that the discrete Euler-Poincaré representation (22)-(24) provide a first order integrator approximating the flow of the continuous Euler-Poincaré representation of the equations motion for a point mass on the n-sphere with potential V . If there is symmetry in the w ∈ so(n + 1) direction, i.e., if V (exp(sw)q) = V (q),
0 −0.2
z
−0.4 −0.6 −0.8 1
for all q ∈ S n , then the “angular momentum in the w direction” I = hMk , wi = Pso(n+1) (pk qkT ), w is conserved by the discrete symmetric representation (19)-(21) and the discrete Euler-Poincaré representation (22)-(24). This can also be shown directly as ∂V T ,w hMk+1 , wi = Mk − h2Pso(n+1) (qk+1 )qk+1 ∂q = hMk , wi .
0.5 0.5
0 0 −0.5
−0.5
y
x
−15
10
x 10
8
2
Example 2: The discrete equations of motion for S ⊂ R3 . If we let f ∈ R3 then Rodrigues’ formula gives N ·q
exp(f × ) = I +
6
sin kf k × 1 − cos kf k × 2 f + (f ) . kf k kf k2
4
2
Thus the hyperbolic sine sinh : so(3) → so(3) is given by sinh(f × ) = 21 (exp(f × ) − exp(−f × )) =
sin kf k × f . kf k
We notice that sinh is onto the set {a× ∈ so(3) | kak ≤ 1}. For sinh(f × ) = a× , for a ∈ R3 , we thus get sin kf k f = a, kf k
and thus the inverse is given by sin−1 kak × a . kak
Using Rodrigues’ formula we therefore get 1 − cos(sin−1 kak) × 2 (a ) , kak2 p 1 − 1 − kak2 × 2 = I +a× + (a ) . (25) kak2
exp(sinh−1 (a× )) = I + a× +
If we let (·)∨ : so(3) → R3 be the inverse of (·)× : R3 → so(3) the discrete Euler-Poincaré representation (22)-(24) for this system becomes with Nk× = Mk Fk = exp sinh−1 hNk× , (26) qk+1 = Fk qk ,
(27)
Nk+1 = Nk − h2 Pso(3) = Nk − hqk+1 ×
−4
0
10
20
30
40
50
t
60
70
80
90
100
kf k = sin−1 kak,
⇒
sinh−1 (a× ) =
−2
Fig. 1. A 3D solution (upper figure) of the discrete Euler-Poincaré representation (26)-(28). The trajectory stays on S 2 by construction. Up to machine precision the condition Nk · qk = 0 is satisfied along the flow (lower figure).
which shows that sin kf k = kak
0
∂V T (qk+1 )qk+1 ∂q
∂V (qk+1 ), ∂q
∨ (28)
⊥ 3 where qk ∈ S 2 and Nk× ∈ g⊥ qk . Since gqk = {v ∈ R | v·q0 = × 0} we have that Nk is perpendicular to qk . We notice that using (25) the discrete dynamics is explicitly defined, and thus fast computation is guaranteed. In [8] similar discrete equations of motion were derived for systems on S 2 × S 2 × . . . × S2. We have implemented the discrete Euler-Poincaré representation (22)-(24) for S 2 ⊂ R3 with a potential function V given by
V (x, y, z) = gz, which corresponds to a spherical pendulum in uniform gravity, with g being the acceleration due to gravity. Using the initial condition q0 = (0, 1, 0)T ,
N0 = (1, 0, 0.5)T ,
which satisfies N0 · q0 = 0, we have generated the discrete flow from time 0 to time 100 using as step size h = 0.01. In Figure 1 we have plotted the 3D flow and the quantity N · q. As expected the condition N · q = 0 is satisfied (up to machine precision) along the flow. In Figure 2 we
have plotted the time evolution of the z component of the angular momentum (upper figure) and the Hamiltonian (lower figure). As expected from the previous sections the angular momentum is seen to be conserved. The energy shows a bounded behavior as a symplectic method would, see [19].
1 0.9 0.8 0.7 0.6
I
0.5 0.4 0.3 0.2 0.1 0
0
10
20
30
40
0
10
20
30
40
50
60
70
80
90
100
50
60
70
80
90
100
t
1.8 1.6 1.4 1.2
H
1 0.8 0.6 0.4 0.2 0
t
E D is seen Fig. 2. The z component of the angular momentum I = M, e× z to be conserved along the discrete flow of the discrete generalized EulerPoincaré equations (22)-(24) (upper figure). The discrete time evolution of the Hamiltonian H = 21 hM, M i + 2gz shows a bounded behavior (lower figure).
IV. C ONCLUSION In this paper we analyzed discrete Clebsch optimal control problems. We discussed their solution generating equations, in terms of the discrete symmetric representation and the discrete Euler-Poincaré representation. We applied the theory to two examples, mechanical systems on matrix Lie groups and mechanical systems on the n-sphere. The method provides integrators that are globally defined symplectic maps which are fairly simple in structure (explicit or almost explicit) and thus easy to implement with guaranteed fast computation and good conservation properties.
R EFERENCES [1] J. Moser and A. P. Veselov, “Discrete versions of some classical integrable systems and factorization of matrix polynomials,” Communications in Mathematical Physics, vol. 139, pp. 217–243, 1991. [2] J. Marsden and M. West, “Discrete mechanics and variational integrators,” Acta Numerica, vol. 10, pp. 357–514, 2001. [3] J. R. Cardoso and F. Silva Leite, “The Moser-Veselov equation,” Linear Algebra and its Applications, vol. 360, pp. 237–248, 2003. [4] J. E. Marsden, S. Pekarsky, and S. Shkoller, “Discrete Euler-Poincare and Lie-Poisson equations,” Nonlinearity, vol. 12, pp. 1647–1662, 1999. [5] T. Lee, M. Leok, and N. H. McClamroch, “A Lie group variational integrator for the attitude dynamics of a rigid body with applications to the 3D pendulum,” in IEEE Conf. on Control Applications, 2005, pp. 962–967. [6] T. Lee, M. Leok, and N. McClamroch, “Lie group variational integrators for the full body problem,” Comput. Methods Appl. Mech. Engrg., vol. 196, no. 29-30, pp. 2907–2924, 2007. [7] N. Nordkvist and A. K. Sanyal, “A Lie group variational integrator for rigid body motion in SE(3) with applications to underwater vehicles,” in IEEE Conf. on Decision and Control, Atlanta, GA, Dec. 2010. [8] T. Lee, M. Leok, and N. McClamroch, “Lagrangian mechanics and variational integrators on two-spheres,” International Journal for Numerical Methods in Engineering, vol. 79, no. 9, pp. 1147–1174, 2009. [9] P. E. Crouch and R. Grossman, “Numerical integration of ordinary differential equations on manifolds,” Journal of Nonlinear Science, vol. 3, pp. 1–33, 1993. [10] F. Gay-Balmaz and T. S. Ratiu, “Clebsch optimal control formulation in mechanics,” Journal of Geometric Mechanics, vol. 3, no. 1, pp. 41–79, 2011. [11] A. M. Bloch, P. E. Crouch, N. Nordkvist, and A. K. Sanyal, “Embedded geodesic problems and optimal control for matrix Lie groups,” Journal of Geometric Mechanics, vol. 3, no. 2, pp. 197–223, 2011. [12] N. Nordkvist, P. E. Crouch, and A. M. Bloch, “Continuous and discrete embedded optimal control problems and their application to the analysis of Clebsch optimal control problems and mechanical systems,” 2012, preprint. [13] N. Nordkvist, P. E. Crouch, A. M. Bloch, and A. K. Sanyal, “Embedded optimal control problems,” in IEEE Conf. on Decision and Control, Orlando, FL, Dec. 2011. [14] C. J. Cotter and D. D. Holm, “Continuous and discrete Clebsch variational principles,” Foundations of Computational Mathematics, vol. 9, no. 2, pp. 221–242, 2009. [15] N. Bou-Rabee and J. E. Marsden, “Hamilton-Pontryagin integrators on Lie groups part I: Introduction and structure-preserving properties,” Foundations of Computational Mathematics, vol. 9, pp. 197–219, 2008. [16] V. G. Boltyanskii, Optimal Control of Discrete Systems. New York: John Wiley, 1978. [17] J. E. Marsden and T. S. Ratiu, Introduction to Mechanics and Symmetry, 2nd ed. New York: Springer Verlag, 1999. [18] P. E. Crouch, N. Nordkvist, and A. K. Sanyal, “Optimal control and geodesics on matrix Lie groups,” in Proc. 9th Portuguese Conference on Automatic Control - CONTROLO’2010, Coimbra, Portugal, Sept. 2010. [19] E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration. New York: Springer Verlag, 2002.
Abstract— In this paper we analyze a class of discrete optimal control problems. These systems are discretizations of a class of optimal control problems defined on invariant submanifolds which we denote embedded optimal control problems. We analyze a particular subset of these called discrete Clebsch optimal control problems where the invariant manifolds are group orbits. The generating Hamiltonian equations for such systems are analyzed. The analysis provides a large class of geometric integrators for mechanical systems. We apply the theory to two example systems: mechanical systems on matrix Lie groups and mechanical systems on the n-sphere.
I. I NTRODUCTION In this paper we analyze a class of discrete optimal control problems which are discretizations of a class of optimal control problems defined on invariant submanifolds. The latter class of problems we denote embedded optimal control problems. We analyze a particular subset of these called discrete Clebsch optimal control problems. For such systems the invariant manifolds are group orbits. The generating Hamiltonian equations for such systems are analyzed. We apply the theory to two example systems: mechanical systems on matrix Lie groups and mechanical systems on the n-sphere. This geometric approach provides globally defined integrators for a large class of mechanical systems on manifolds. Our theory provides an elegant and unified derivation of a large class of integrators for mechanical systems. This class includes many well known examples such as the MoserVeselov approach (see [1]) to integrating the motion of ndimensional rigid bodies and the discrete variational approach to deriving integrators surveyed in Marsden-West [2]. A globally defined discretization of the n-dimensional rigid body problem was derived and analyzed in [1] from a discrete variational analysis; the discrete equations of motion were analyzed further in [3]. A discrete analog of the EulerPoincaré equations was derived in [4], and closely related globally defined Lie group variational integrators for a rigid body in a conservative force field were studied in [5], [6], [7]. In [8] a similar globally defined discrete variational integrator for mechanical systems on S 2 × S 2 × · · · × S 2 was developed. A globally defined integrator does not depend on local coordinate systems for its implementation, and thus switching between coordinate systems is avoided. The Nikolaj Nordkvist, Mathematics and Sciences Division, Leeward Community College, Pearl City, HI 96782 [email protected] Peter E. Crouch, Electrical Engineering, University of Hawai‘i at M¯anoa, Honolulu, HI 96822 [email protected] Anthony M. Bloch, Department of Mathematics, University of Michigan, Ann Arbor, MI 48109 [email protected]. Supported in part by NSF grants DMS-0806756, DMS-0907949, and DMS-1207693
Crouch-Grossman methods, introduced in [9], provide a global framework for discrete analysis of rigid body equations but these methods fail to be symplectic. In [10] continuous Clebsch optimal control problems were defined and analyzed while [11] and [12] considered the generalization of embedded optimal control problems, the latter also analyzing their discrete definition. Examples of continuous embedded optimal control problems were examined in [13]. In [14] Clebsch integrators were proposed based on a discretization of the Clebsch action principle—as opposed to a discretization of the Clebsch optimal control problem which is the approach of this paper. The theory of Hamilton-Pontryagin integrators, presented in [15], provides a discretization based on a discrete Hamilton-Pontryagin variational principle. In particular in this paper we derive a simple algorithm for integrating geodesic flows on the n-sphere. Simulation are performed for the special case of S 2 which illustrates the geometric properties of the algorithm. II. D ISCRETE C LEBSCH OPTIMAL CONTROL PROBLEMS In this paper we will consider a special case of discrete dynamics of the form qk+1 = f (qk , uk ) where qk are the configuration variables and uk are the control variables. In [12] we analyze discrete optimal control problems with this dynamics. In particular we define discrete embedded optimal control problems. As discussed in [12] embedded optimal control problems are a generalization of optimal control problems on manifolds which have a natural interpretation on invariant submanifolds where the systems are locally controllable. A particular class of embedded optimal control problems are Clebsch optimal control problems (see [10] for the continuous definition, [12] for the discrete definition) where the invariant manifolds are group orbits. Consider the discrete dynamics given by qk+1 = ΦFk (qk ), where Φ : G × Q → Q, (F, q) 7→ ΦF (q), is the action of the (finite dimensional) Lie group G on a smooth manifold Q. We denote by J : T ∗ Q → g∗ the momentum map of Φ. We have that a group orbit orb(q0 ), q0 ∈ Q, by construction is an invariant set under the discrete dynamics. We denote ¯ : G × orb(q0 ) → orb(q0 ) the induced action on by Φ ¯ : T ∗ orb(q0 ) → g∗ its momentum map. orb(q0 ) and by J In order to be able to apply the results obtained in [12] for discrete embedded/associated optimal control problems we furthermore need to assume that G and orb(q0 ) are embedded submanifolds of some vector spaces and that these embeddings are obtained as the preimages at 0 of some smooth submersions (mapping between vector spaces).
We consider the following discrete version of the Clebsch optimal control problem: Problem A (Discrete Clebsch Optimal control problem): Minimize N −1 X
ℓd (qk , Fk )
k=0
over all sequences {qk }N k=0 , qk ∈ Q, satisfying qk+1 = ΦFk (qk ) with Fk ∈ G and fixed endpoints q0 , qN ∈ Q. For the discrete Clebsch optimal control problem to be well defined we need qN ∈ orb(q0 ). The associated optimal control problem is then the following:
Proof: We will only show this for (2). The proof is the same for (4). From Proposition 6.3.2 in [17] we have that ∗ ∗ T ∗ ΦF −1 = T ∗ Φ−1 Fk : T Q → T Q is a symplectic map. The k d map (qk , pk ) 7→ (qk , pk + λ0 δℓ δq (qk , uk )) is also symplectic since this is the time one flow of the Hamiltonian vector field on T ∗ Q with Hamiltonian H(q, p) = −λ0 ℓd (q, uk ). Since (2) is the composition of these two symplectic maps it itself must be symplectic. We define Mk as δℓd Mk := J(αk ) + J (qk , Fk ) , δq
for αk ∈ T ∗ Q. By construction we have Mk ∈ g⊥ qk = {ξ ∈ Problem B (Discrete Clebsch assoc. opt. control problem): g∗ | hξ, ηi = 0 ∀ η ∈ gq } since J(Tq∗ Q) ⊂ g⊥ . q If αk is a Minimize solution to the discrete symmetric representation (1)-(2) then N −1 we get for a left action that Mk satisfies the equation X ℓd (qk , Fk ) δℓd (q , F ) M =J(q , p ) + J k=0 k+1 k+1 k+1 k+1 k+1 δq N over all sequences {qk }k=0 , qk ∈ orb(q0 ), satisfying qk+1 = δℓd ∗ −1 p + =J Φ (q ), T Φ (q , F ) ¯ k Fk k k k ΦFk (qk ) Fk ΦFk (qk ) with Fk ∈ G and fixed endpoints q0 ∈ Q, qN ∈ δq orb(q0 ). δℓd (qk+1 , Fk+1 ) +J From [12] we have the following results regarding the δq extremals of the discrete Clebsch optimal control problems: δℓd ∗ =AdF −1 J qk , pk + (qk , Fk ) k δq Theorem 1: Assume that Φ is a left action. A normal δℓd extremal for the discrete Clebsch optimal control problem (qk+1 , Fk+1 ) . +J as given by applying the discrete maximum principle (see δq [16]) is a solution of These same calculations could be carried out for (3). We can δℓ δℓd therefore express Theorem 1 in terms of qk , Fk , and Mk as d ∗ L Fk (qk , Fk ) = Tid (qk , Fk ) , (1) J(αk ) + J follows: δq δF δℓd Theorem 2: Assume that Φ is a left action. A normal . (2) αk+1 = T ∗ ΦF −1 αk + verαk extremal for the discrete Clebsch optimal control problem k δq as given by applying the discrete maximum principle is a for αk ∈ T ∗ Q. The pullback bundle i∗orb(q0 ) (T ∗ Q) is solution of invariant under the flow of (2). δℓd ∗ An extremal for the discrete Clebsch associated optimal Mk = Tid L Fk (qk , Fk ) , (5) δF control problem is a solution of ¯ qk+1 = ΦFk (qk ), (6) δℓd δℓ ¯ k) + J ¯ δℓd (qk , Fk ) = T ∗ LF J(β (q , F ) , (3) d k k ∗ id k ¯ (qk+1 , Fk+1 ) , (7) Mk+1 = AdF −1 (Mk ) + J δF δq k δq ¯ δℓ d ¯ −1 βk + ver ¯ βk ¯ . (4) for qk ∈ Q and Mk ∈ g⊥ . βk+1 = T ∗ Φ Fk qk δq An extremal (must be normal) for the discrete Clebsch for βk ∈ T ∗ orb(q0 ). The discrete Clebsch associated optimal associated discrete optimal control problem is a solution of control problem does not admit abnormal extremals. δℓd ∗ If αk ∈ i∗orb(q0 ) (T ∗ Q) is a solution to (1)-(2) then βk = Mk = Tid LFk (qk , Fk ) , (8) δF T ∗ iorb(q0 ) αk ∈ T ∗ orb(q0 ) is a solution to (3)-(4). (9) qk+1 = ΦFk (qk ), δℓ For a right action we get very similar results, the only ¯ d ¯ ∗ Mk+1 = Ad∗F −1 (Mk ) + J (10) difference being that Tid LFk in (1) and (3) is replaced with ¯ (qk+1 , Fk+1 ) , k δq ∗ Tid RFk . We implicitly assume that (1) can be inverted to ⊥ give Fk as a function of αk . We refer to (1)-(2) as the for qk ∈ orb(q 0 ) and Mk ∈ gqk . discrete symmetric representation of the discrete extremal If qk , Mk is a solution to (8)-(10) then qk , Mk is a generating equations. solution to (5)-(10). Proposition 1: The discrete maps (2) and (4) are symplecFor a right action we would get very similar equations, the ∗ tic. only difference being that in (5) and (8) Tid LFk is replaced
∗ with Tid RFk and in (7) and (10) Ad∗F −1 (Mk ) is replaced k with Ad∗Fk . We refer to (5)-(7) as the discrete EulerPoincaré representation of the discrete extremal generating equations. In this theorem it is implicitly assumed that (5) can be solved to give Fk as a function of qk and Mk .
III. E XAMPLES We now consider two basic applications which illustrate the general theory. A. Discrete Clebsch optimal control problems on matrix Lie groups Consider Q = gl(n) with the right action Φ of a matrix Lie group G as defined by right matrix multiplication, that is Φg (q) = qg for q ∈ gl(n) and g ∈ G. The orbit that arises when q0 ∈ G is the matrix Lie group G itself. Given a self-adjoint positive definite (with respect to the trace inner product) operator Σ : g → g we let Σd : gl(n) → gl(n) denote a self-adjoint operator that satisfies P ◦ Σd |g = Σ for some positive definite operator Σ : g → g∗ ≃ g; this is equivalent to hu, Σd (u)i = hu, Σ(u)i for all u ∈ g. Such an extension can always be defined as the choice Σd = Σ ◦ P demonstrates. Let us consider the following discrete cost function 11 h2
ℓd (qk , Fk ) =
h(Fk − I), Σd (Fk − I)i − hV (qk ), (11)
for qk ∈ gl(n) and Fk ∈ G. In this problem Fk can be interpreted as an approximation to exp(huk ) where uk ≈ q −1 (tk )q(t ˙ k ) ∈ g is the discrete control at time tk and h = tk+1 − tk is a (fixed) step size satisfying T = N h. This discrete cost function therefore approximates Z tk+1 ℓ q(t), u(t) dt, (12) ℓd (qk , Fk ) ≈ tk
where ℓ(q, u) = 21 hu, Σ(u)i−V (q). When modeling the discrete dynamics of a mechanical system the term 21 hu, Σ(u)i will correspond to the kinetic energy and V (q) will correspond to the potential energy. From straightforward calculations, which can be found in detail in [12], we get that T ∗ ΦF −1 (q, p) = (qg, pg −T ), k
∗ Tid R Fk
J(q, p) = P(q T p), δℓd (qk , Fk ) = h1 P Σd (Fk − I)FkT , δF
where P is the orthogonal projection P : gl(n) → g∗ ≡ g with respect to the trace inner product. Thus from Theorem 1 we get that the discrete symmetric representation for this system becomes T 2 T ∂V (qk ) = P Σd (Fk − I)FkT , (13) hP(qk pk )−h P qk ∂q qk+1 = qk Fk , (14) ∂V pk+1 = qk − h (qk ) Fk−T . (15) ∂q
From Theorem 2 we get that the discrete Euler-Poincaré representation of this system is hMk = P Σd (Fk − I)FkT , qk+1 = qk Fk ,
Mk+1 =
Ad∗Fk (Mk )
− hP
∂V T qk+1
∂q
(qk+1 ) .
For F ∈ G, M ∈ g∗ = g, and η ∈ g we get
hAd∗F (M ), ηi = hM, AdF (η)i = M, F ηF −1
= F T M F −T , η
= P(F T M F −T ), η ,
and thus we have the identity
Ad∗F (M ) = P(F T M F −T ). Using this the discrete Euler-Poincaré representation for this system can be written as (16) hMk = P Σd (Fk − I)FkT , qk+1 = qk Fk , (17) ∂V T (qk+1 ) . (18) Mk+1 = P(FkT Mk Fk−T ) − hP qk+1 ∂q In [12] we show that the Euler-Poincaré representation (16)(18) is at least a first order approximation of the flow of the continuous Euler-Poincaré representation with cost function ℓ(q, u). Let us consider the case where V = 0 and Φ is the action ˜ is the of Gl(n) on Q given by right multiplication and Φ action of Gl(n) on Q given by left multiplication. Then we ˜ g˜ and since ℓd when V = 0 have that Φg commutes with Φ is independent of qk we get that, as in [12], ˜ k , pk ) = pk q T , J(q k is conserved along the flow of the discrete symmetric rep¯˜ , p ) = P(p q T ) is resentation (13)-(15). Therefore J(q k k k k conserved along the flow of (3)-(4). Since Mk = J(qk , pk ) when V = 0 we get that P(pk qkT ) = P(qk−T qkT pk qkT ) = P qk−T P(qkT pk )qkT = P((qk−T Mk qkT ) = Ad∗q−1 (Mk ) and k therefore Ad∗q−1 (Mk ) is conserved along the flow of the disk crete Euler-Poincaré representation (16)-(18). These results regarding conserved quantities were derived directly, using the discrete extremal generating equations only, in [18]. We now consider the special case of the orthgonal group which of course is of particular interest in applications. Example 1: Discrete generalized rigid body equations. Consider the Lie group G = SO(n) with Lie algebra g = so(n) and take Rk ∈ SO(n) and Ω ∈ so(n) to be group and algebra elements. First let us consider the free generalized rigid body, i.e., the case when V = 0. Let Σ : so(n) → so(n) be given by Σ(Ω) = Jd Ω + ΩJd , where Jd is a symmetric n×n matrix which makes Σ positive definite. The orthogonal projection Pso(n) : gl(n) → so(n) is given by Pso(n) (A) = 1 T 2 (A − A ). We choose Σd (A) = 2Jd A,
which gives Pso(n) ◦ Σd (Ω) = 12 (2Jd Ω − (2Jd Ω)T ) = Σ(Ω), showing that this choice of Σd indeed satisfies the condition Pso(n) ◦ Σd |g = Σ. Since Pso(n) (2J(Fk − I)FkT ) = Pso(n) (2J − 2JFkT ) = Fk J − JFkT , the discrete Euler-Poincaré representation (16)-(18) for this system then becomes hMk = Fk Jd − Jd FkT ,
B. Discrete Clebsch optimal control problems on the nsphere Consider Q = Rn+1 with the group G = SO(n+1) acting on Q by left matrix multiplication. Then S n is a group orbit. We define the discrete cost function as ℓd (qk , Fk ) =
T ∗ Φg−1 (q, p) = (gq, g T p),
Mk+1 = FkT Mk Fk .
(Jω)× = ω × J − J ω × , where J = 12 tr(J)I − J is a modified inertia matrix, we have Jd = 2J . We define the discrete angular velocity (in the body fixed frame) as for the continuous model Mk =: Σ(ωk× ) = 12 (Jωk )× . Since (F η)× = F η × F T , for F ∈ SO(3) and η ∈ R3 , we get ∂V (Jωk+1 )× = 2Mk+1 = 2FkT Mk Fk − 2hP (Rk+1 ) ∂R T × × = (Fk ωk ) − τ (Rk+1 ), and therefore the Euler-Poincaré representation (16)-(18) for the rigid body translates to h(Jωk )× = Fk J − J FkT , Rk+1 = Rk Fk , Jωk+1 = FkT Jωk − τ (Rk+1 ). These are the discrete rigid body equations analyzed in [5] in their Hamiltonian form. The discrete rigid body equations in SE(3), as studied in [7], reduce to these equations when there is no translational velocity. In [7] it was shown that this discrete map is a useful integrator as it provides (at least) a first order approximation to the flow of the continuous rigid body equations. The first of these equations is an implicit equation in the unknown Fk ∈ SO(3). The numerical solution of this implicit equation has been discussed in [5]. In [3] explicit solutions to this equation were derived in terms of solutions to an algebraic Riccati equation. The remaining equations are explicit and thus are straightforward to implement numerically.
hFk − I, Fk − Ii − h2V (qk ).
In [12] we show that this discrete cost function approximates the cost of a continuous optimal control problem with solutions being flows of the mechanical system corresponding to a point mass on S n with potential energy V . From [12] we have that
Rk+1 = Rk Fk ,
These are the discrete rigid body equations studied in [1]. Now let us consider a true rigid body, which is the special case when G = SO(3) and V is the potential energy which may be nonzero. The inertia matrix J is a 3×3 matrix which is positive definite and symmetric. Let (·)× : R3 → so(3) be defined by a× b = a × b for a, b ∈ R3 with × being the cross product. As in [7] we define Σ(ω × ) = 12 (Jω)× , which is positive definite by construction, and since we have the identity
1 h
∗ Tid LFk
J(q, p) = Pso(n+1) (pq T ), δℓd (qk , Fk ) = h1 Pso(n+1) (I − FkT ) δF = h1 (Fk − FkT ).
Therefore we get from Theorem 1 that the discrete symmetric representation for this system is T ∂V T 1 p − h2 (F − F ) = hP (q ) qk , (19) k k k so(n+1) k 2 ∂q qk+1 = Fk qk , (20) ∂V pk+1 = Fk pk − h2Fk (qk ), (21) ∂q for qk , pk ∈ Rn+1 . Using Fk = exp(fk ), fk ∈ so(n + 1), we have T 1 1 2 (Fk − Fk ) = 2 exp(fk ) − exp(−fk ) = sinh(fk ).
Thus the discrete symmetric representation for this system formally becomes T ∂V −1 hPso(n+1) pk − h2 Fk = exp sinh , (qk ) qk ∂q qk+1 = Fk qk , ∂V pk+1 = Fk pk − h2Fk (qk ), ∂q Since Ad∗F −1 (Mk ) = FkT Mk Fk = k
=
1 2h (Fk
1 2h Fk (Fk
− FkT )Fk
− FkT ) = Mk ,
we get from Theorem 2 that the discrete Euler-Poincaré representation for this system is hMk = 12 (Fk − FkT ), qk+1 = Fk qk , Mk+1 = Mk − h2Pso(n+1)
(22) (23)
∂V T (qk+1 )qk+1 , ∂q
where M0 ∈ g⊥ q0 . This can formally be rewritten as Fk = exp sinh−1 (hMk ) , qk+1 = Fk qk ,
Mk+1 = Mk − h2Pso(n+1)
∂V T (qk+1 )qk+1 . ∂q
(24)
In [12] we show that the discrete Euler-Poincaré representation (22)-(24) provide a first order integrator approximating the flow of the continuous Euler-Poincaré representation of the equations motion for a point mass on the n-sphere with potential V . If there is symmetry in the w ∈ so(n + 1) direction, i.e., if V (exp(sw)q) = V (q),
0 −0.2
z
−0.4 −0.6 −0.8 1
for all q ∈ S n , then the “angular momentum in the w direction” I = hMk , wi = Pso(n+1) (pk qkT ), w is conserved by the discrete symmetric representation (19)-(21) and the discrete Euler-Poincaré representation (22)-(24). This can also be shown directly as ∂V T ,w hMk+1 , wi = Mk − h2Pso(n+1) (qk+1 )qk+1 ∂q = hMk , wi .
0.5 0.5
0 0 −0.5
−0.5
y
x
−15
10
x 10
8
2
Example 2: The discrete equations of motion for S ⊂ R3 . If we let f ∈ R3 then Rodrigues’ formula gives N ·q
exp(f × ) = I +
6
sin kf k × 1 − cos kf k × 2 f + (f ) . kf k kf k2
4
2
Thus the hyperbolic sine sinh : so(3) → so(3) is given by sinh(f × ) = 21 (exp(f × ) − exp(−f × )) =
sin kf k × f . kf k
We notice that sinh is onto the set {a× ∈ so(3) | kak ≤ 1}. For sinh(f × ) = a× , for a ∈ R3 , we thus get sin kf k f = a, kf k
and thus the inverse is given by sin−1 kak × a . kak
Using Rodrigues’ formula we therefore get 1 − cos(sin−1 kak) × 2 (a ) , kak2 p 1 − 1 − kak2 × 2 = I +a× + (a ) . (25) kak2
exp(sinh−1 (a× )) = I + a× +
If we let (·)∨ : so(3) → R3 be the inverse of (·)× : R3 → so(3) the discrete Euler-Poincaré representation (22)-(24) for this system becomes with Nk× = Mk Fk = exp sinh−1 hNk× , (26) qk+1 = Fk qk ,
(27)
Nk+1 = Nk − h2 Pso(3) = Nk − hqk+1 ×
−4
0
10
20
30
40
50
t
60
70
80
90
100
kf k = sin−1 kak,
⇒
sinh−1 (a× ) =
−2
Fig. 1. A 3D solution (upper figure) of the discrete Euler-Poincaré representation (26)-(28). The trajectory stays on S 2 by construction. Up to machine precision the condition Nk · qk = 0 is satisfied along the flow (lower figure).
which shows that sin kf k = kak
0
∂V T (qk+1 )qk+1 ∂q
∂V (qk+1 ), ∂q
∨ (28)
⊥ 3 where qk ∈ S 2 and Nk× ∈ g⊥ qk . Since gqk = {v ∈ R | v·q0 = × 0} we have that Nk is perpendicular to qk . We notice that using (25) the discrete dynamics is explicitly defined, and thus fast computation is guaranteed. In [8] similar discrete equations of motion were derived for systems on S 2 × S 2 × . . . × S2. We have implemented the discrete Euler-Poincaré representation (22)-(24) for S 2 ⊂ R3 with a potential function V given by
V (x, y, z) = gz, which corresponds to a spherical pendulum in uniform gravity, with g being the acceleration due to gravity. Using the initial condition q0 = (0, 1, 0)T ,
N0 = (1, 0, 0.5)T ,
which satisfies N0 · q0 = 0, we have generated the discrete flow from time 0 to time 100 using as step size h = 0.01. In Figure 1 we have plotted the 3D flow and the quantity N · q. As expected the condition N · q = 0 is satisfied (up to machine precision) along the flow. In Figure 2 we
have plotted the time evolution of the z component of the angular momentum (upper figure) and the Hamiltonian (lower figure). As expected from the previous sections the angular momentum is seen to be conserved. The energy shows a bounded behavior as a symplectic method would, see [19].
1 0.9 0.8 0.7 0.6
I
0.5 0.4 0.3 0.2 0.1 0
0
10
20
30
40
0
10
20
30
40
50
60
70
80
90
100
50
60
70
80
90
100
t
1.8 1.6 1.4 1.2
H
1 0.8 0.6 0.4 0.2 0
t
E D is seen Fig. 2. The z component of the angular momentum I = M, e× z to be conserved along the discrete flow of the discrete generalized EulerPoincaré equations (22)-(24) (upper figure). The discrete time evolution of the Hamiltonian H = 21 hM, M i + 2gz shows a bounded behavior (lower figure).
IV. C ONCLUSION In this paper we analyzed discrete Clebsch optimal control problems. We discussed their solution generating equations, in terms of the discrete symmetric representation and the discrete Euler-Poincaré representation. We applied the theory to two examples, mechanical systems on matrix Lie groups and mechanical systems on the n-sphere. The method provides integrators that are globally defined symplectic maps which are fairly simple in structure (explicit or almost explicit) and thus easy to implement with guaranteed fast computation and good conservation properties.
R EFERENCES [1] J. Moser and A. P. Veselov, “Discrete versions of some classical integrable systems and factorization of matrix polynomials,” Communications in Mathematical Physics, vol. 139, pp. 217–243, 1991. [2] J. Marsden and M. West, “Discrete mechanics and variational integrators,” Acta Numerica, vol. 10, pp. 357–514, 2001. [3] J. R. Cardoso and F. Silva Leite, “The Moser-Veselov equation,” Linear Algebra and its Applications, vol. 360, pp. 237–248, 2003. [4] J. E. Marsden, S. Pekarsky, and S. Shkoller, “Discrete Euler-Poincare and Lie-Poisson equations,” Nonlinearity, vol. 12, pp. 1647–1662, 1999. [5] T. Lee, M. Leok, and N. H. McClamroch, “A Lie group variational integrator for the attitude dynamics of a rigid body with applications to the 3D pendulum,” in IEEE Conf. on Control Applications, 2005, pp. 962–967. [6] T. Lee, M. Leok, and N. McClamroch, “Lie group variational integrators for the full body problem,” Comput. Methods Appl. Mech. Engrg., vol. 196, no. 29-30, pp. 2907–2924, 2007. [7] N. Nordkvist and A. K. Sanyal, “A Lie group variational integrator for rigid body motion in SE(3) with applications to underwater vehicles,” in IEEE Conf. on Decision and Control, Atlanta, GA, Dec. 2010. [8] T. Lee, M. Leok, and N. McClamroch, “Lagrangian mechanics and variational integrators on two-spheres,” International Journal for Numerical Methods in Engineering, vol. 79, no. 9, pp. 1147–1174, 2009. [9] P. E. Crouch and R. Grossman, “Numerical integration of ordinary differential equations on manifolds,” Journal of Nonlinear Science, vol. 3, pp. 1–33, 1993. [10] F. Gay-Balmaz and T. S. Ratiu, “Clebsch optimal control formulation in mechanics,” Journal of Geometric Mechanics, vol. 3, no. 1, pp. 41–79, 2011. [11] A. M. Bloch, P. E. Crouch, N. Nordkvist, and A. K. Sanyal, “Embedded geodesic problems and optimal control for matrix Lie groups,” Journal of Geometric Mechanics, vol. 3, no. 2, pp. 197–223, 2011. [12] N. Nordkvist, P. E. Crouch, and A. M. Bloch, “Continuous and discrete embedded optimal control problems and their application to the analysis of Clebsch optimal control problems and mechanical systems,” 2012, preprint. [13] N. Nordkvist, P. E. Crouch, A. M. Bloch, and A. K. Sanyal, “Embedded optimal control problems,” in IEEE Conf. on Decision and Control, Orlando, FL, Dec. 2011. [14] C. J. Cotter and D. D. Holm, “Continuous and discrete Clebsch variational principles,” Foundations of Computational Mathematics, vol. 9, no. 2, pp. 221–242, 2009. [15] N. Bou-Rabee and J. E. Marsden, “Hamilton-Pontryagin integrators on Lie groups part I: Introduction and structure-preserving properties,” Foundations of Computational Mathematics, vol. 9, pp. 197–219, 2008. [16] V. G. Boltyanskii, Optimal Control of Discrete Systems. New York: John Wiley, 1978. [17] J. E. Marsden and T. S. Ratiu, Introduction to Mechanics and Symmetry, 2nd ed. New York: Springer Verlag, 1999. [18] P. E. Crouch, N. Nordkvist, and A. K. Sanyal, “Optimal control and geodesics on matrix Lie groups,” in Proc. 9th Portuguese Conference on Automatic Control - CONTROLO’2010, Coimbra, Portugal, Sept. 2010. [19] E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration. New York: Springer Verlag, 2002.
