Cédric M Campos
email: cedricmc@uva@icmat.es

AND Sina Ober-Blöbaum
email: sinaob@math.upb.de

Emmanuel Trélat
email: emmanuel.trelat@upmc.fr

HIGH ORDER VARIATIONAL INTEGRATORS IN THE OPTIMAL CONTROL OF MECHANICAL SYSTEMS

Keywords: 1991 Mathematics Subject Classification. Primary 65P10, Secondary 65L06, 65K10, 49Mxx optimal control, mechanical systems, geometric integration, variational integrator, high order, Runge-Kutta, direct methods, commutation property

In recent years, much effort in designing numerical methods for the simulation and optimization of mechanical systems has been put into schemes which are structure preserving. One particular class are variational integrators which are momentum preserving and symplectic. In this article, we develop two high order variational integrators which distinguish themselves in the dimension of the underling space of approximation and we investigate their application to finite-dimensional optimal control problems posed with mechanical systems. The convergence of state and control variables of the approximated problem is shown. Furthermore, by analyzing the adjoint systems of the optimal control problem and its discretized counterpart, we prove that, for these particular integrators, dualization and discretization commute.

Introduction

In practice, solving an optimal control problem requires the a priori choice of a numerical method. Many approaches do exist, that are either based on a direct discretization of the optimal control problem, resulting into a nonlinear programming problem, or based on the application of the Pontryagin Maximum Principle, reducing into a boundary value problem. The first class of approaches are called direct, and the second ones, based on the preliminary use of the Pontryagin maximum principle, are called indirect. It can be noted that indirect methods, although very precise, suffer from an important drawback: Unless one already has a good knowledge of the optimal solution, they are very difficult to initialize since they are extremely sensitive to the initial guess. Although many solutions exist in order to carefully initialize a shooting method (see [START_REF]Contrôle optimal[END_REF][START_REF]Optimal control and applications to aerospace: some results and challenges[END_REF]), in most of engineering applications direct approaches are preferred due to their simplicity and robustness. Roughly speaking, direct methods consist of (1) discretizing first the cost and the differential system, in order to reduce the optimal control problem to a usual nonlinear minimization problem with constraints, also called nonlinear programming problem (NLP), with dimension inversely proportional to the smoothness of the discretization; (2) and then dualizing, by applying for instance a Lagrange-Newton method to the NLP, deriving the Karush-Kuhn-Tucker equations (KKT), also called discrete adjoint system, and applying a Newton method to solve the resulting optimality system.

Many variants exist, e.g. [START_REF] Betts | Survey of numerical methods for trajectory optimization[END_REF]). In contrast, indirect methods consist of (1) first dualizing the optimal control problem to derive the adjoint system, by applying the Pontryagin Maximum Principle (PMP) (or, equivalently, the Lagrange multipliers necessary condition for optimality in infinite dimension), (2) and then discretizing, by applying a shooting method (that is, a Newton method composed with a numerical integration method).

In shorter words, direct methods consist of 1) discretize, 2) dualize, and indirect methods consist of the converse: 1) dualize, 2) discretize. It is natural to wonder whether this diagram is commutative or not, under usual approximation assumptions.

Since the pioneering works of [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF][START_REF] Ross | A roadmap for optimal control: the right way to commute[END_REF], it is well known by now that, in spite of usual assumptions ensuring consistency and stability, direct methods may diverge. In other words discretization and dualization do not commute in general. This is due to a complex interplay with the mesh, to the appearance of spurious highfrequencies ultimately causing the divergence (see [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF] for very simple finite-dimensional linear quadratic problems and [START_REF] Zuazua | Propagation, observation, and control of waves approximated by finite difference methods[END_REF] for infinite-dimensional wave equations).

Several remedies and solutions exist, from which [START_REF] Bonnans | Computation of order conditions for symplectic partitioned Runge-Kutta schemes with application to optimal control[END_REF][START_REF] Gong | Connections between the covector mapping theorem and convergence of pseudospectral methods for optimal control[END_REF][START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF][START_REF] Ross | Legendre pseudospectral approximations of optimal control problems[END_REF] are a representative sample. For instance, the results of [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF] assert the convergence of optimal control problems under specific smoothness and coercivity assumptions provided that the underlying discretization method be based on a Runge-Kutta method. However, the convergence order of the optimal control solution, which depends on the convergence order of the state and the resulting adjoint scheme, is reduced compared to the order of the Runge-Kutta method applied to the state system. Indeed, the discrete state and adjoint system constitutes a symplectic partitioned Runge-Kutta method for which order conditions on the Runge-Kutta coefficients are derived to preserve the convergence rates. Whereas in [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF] the symplectic partitioned Runge-Kutta scheme for state and adjoint is explicitly derived, in a recent overview article [START_REF] Sanz-Serna | Symplectic Runge-Kutta schemes for adjoint equations, automatic differentiation, optimal control and more[END_REF] a proof is given based on a relation between quadratic invariants (that are also naturally preserved by symplectic partitioned Runge-Kutta integrators) and the fulfillment of the KKT equations. The preservation of convergence rates is referred to as the Covector Mapping Principle (CMP) (see e.g. [START_REF] Gong | Connections between the covector mapping theorem and convergence of pseudospectral methods for optimal control[END_REF]). More precisely, the CMP is satisfied if there exists an order-preserving map between the adjoint variables corresponding to the dualized discrete problem (KKT) and the discretized dual problem (discretized PMP). For the class of Legendre pseudospectral methods the CMP is proven if additional closure conditions are satisfied (see [START_REF] Gong | Connections between the covector mapping theorem and convergence of pseudospectral methods for optimal control[END_REF][START_REF] Ross | Legendre pseudospectral approximations of optimal control problems[END_REF]), whereas for Runge-Kutta methods the CMP holds if the order conditions on the Runge-Kutta coefficients derived in [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF] are satisfied. For a detailed discussion on the commutation properties we refer to [START_REF] Ross | A roadmap for optimal control: the right way to commute[END_REF].

While for general dynamical systems, many studies of discretizations of optimal control problems are based on Runge-Kutta methods (see e.g. [START_REF] Dontchev | The Euler approximation in state constrained optimal control[END_REF][START_REF] Dontchev | Second-order Runge-Kutta approximations in control constrained optimal control[END_REF][START_REF]Numerical analysis in optimal control[END_REF][START_REF] Herty | Implicit-explicit Runge-Kutta schemes for numerical discretization of optimal control problems[END_REF][START_REF] Walther | Automatic differentiation of explicit Runge-Kutta methods for optimal control[END_REF]), particularly for mechanical systems, much effort in designing numerical methods for the simulation and optimization of such systems has been put into schemes which are structure preserving in the sense that important qualitative features of the original dynamics are preserved in its discretization (for an overview on structure preserving integration methods see e.g. [START_REF] Hairer | Geometric numerical integration[END_REF]). One special class of structure preserving integrators is the class of variational integrators, introduced in [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF][START_REF] Suris | Hamiltonian methods of Runge-Kutta type and their variational interpretation[END_REF], which are symplectic and momentum-preserving and have an excellent long-time energy behavior. Variational integrators are based on a discrete variational formulation of the underlying system, e.g. based on a discrete version of Hamilton's principle or of the Lagrange-d'Alembert principle for conservative [START_REF]An overview of variational integrators[END_REF][START_REF]Variational time integrators[END_REF] or dissipative [START_REF] Kane | Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems[END_REF] mechanical systems, respectively. They have been further extended to different kind of systems and applications, e.g. towards constrained [START_REF] Cortés | Non-holonomic integrators[END_REF][START_REF] Kobilarov | Geometric discretization of nonholonomic systems with symmetries[END_REF][START_REF] Leyendecker | Variational integrators for constrained dynamical systems[END_REF], non smooth [START_REF] Fetecau | Nonsmooth Lagrangian mechanics and variational collision integrators[END_REF], multirate and multiscale [START_REF] Leyendecker | A variational approach to multirate integration for constrained systems[END_REF][START_REF] Stern | Implicit-explicit variational integration of highly oscillatory problems[END_REF][START_REF] Tao | Nonintrusive and structure preserving multiscale integration of stiff ODEs, SDEs, and Hamiltonian systems with hidden slow dynamics via flow averaging[END_REF], Lagrangian PDE systems [START_REF] Lew | Asynchronous variational integrators[END_REF][START_REF] Marsden | Multisymplectic geometry, variational integrators, and nonlinear PDEs[END_REF] and electric circuits [START_REF] Ober-Blöbaum | Variational integrators for electric circuits[END_REF]. In the cited works, typically quadrature rules of first or second order are used in order to approximate the action functional of the system. To design high order variational integrators, higher order quadrature rules based on polynomial collocation can be employed. Such so called Galerkin variational integrators have been introduced in [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF] and further studied in [START_REF] Campos | High order variational integrators: a polynomial approach[END_REF][START_REF] Hall | Spectral variational integrators[END_REF][START_REF] Leok | General techniques for constructing variational integrators[END_REF][START_REF] Ober-Blöbaum | Galerkin variational integrators and modified symplectic Runge-Kutta methods[END_REF][START_REF] Ober-Blöbaum | Construction and analysis of higher order galerkin variational integrators[END_REF].

For numerically solving optimal control problems by means of a direct method, the use of variational integrators for the discretization of the problem has been proposed in [START_REF] Campos | Higher order variational time discretization of optimal control problems[END_REF][START_REF] Leyendecker | Discrete mechanics and optimal control for constrained systems[END_REF][START_REF] Ober-Blöbaum | Discrete mechanics and optimal control: an analysis[END_REF]. This approach, denoted by DMOC (Discrete Mechanics and Optimal Control), yields a finite-difference type discretization of the dynamical constraints of the problem which by construction preserves important structural properties of the system, like the evolution of momentum maps associated to symmetries of the Lagrangian or the energy behavior. For one class of Galerkin variational integrators, that is equivalent to the class of symplectic partitioned Runge-Kutta methods, the adjoint system and its convergence rates are analyzed in [START_REF] Ober-Blöbaum | Discrete mechanics and optimal control: an analysis[END_REF]. It is shown that, in contrast to a discretization based on standard Runge-Kutta methods in [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF], the convergence order of the discrete adjoint system is the same as for the state system due to the symplecticity of the discretization scheme. In particular, we obtain the same symplectic-momentum scheme for both state and adjoint system, that means that discretization and dualization commute for this class of symplectic schemes and the CMP is satisfied. For general classes of variational integrators, the commutation property is still an open question. The contribution of this work is twofold:

(1) We derive two different kinds of high order variational integrators based on different dimensions of the underling polynomial approximation (Section 3). While the first well-known integrator is equivalent to a symplectic partitioned Runge-Kutta method, the second integrator, denoted as symplectic Galerkin integrator, yields a "new" method which in general, cannot be written as a standard symplectic Runge-Kutta scheme. (2) For the application of high order variational integrators to finite-dimensional optimal control problems posed with mechanical systems, we show the convergence of state and control variables and prove the commutation of discretization and dualization (Sections 4 and 5). The paper is structured as follows: In Section 2 the optimal control problem for a mechanical system is introduced. Its discrete version is formulated in Section 3 based on the derivation of two different kinds of high order variational integrators. The first main result is stated in Section 4: Under specific assumptions on the problem setting we prove the convergence of the primal variables for an appropriate choice of the discrete controls (Theorem 4.1). Along the lines of [START_REF]Numerical analysis in optimal control[END_REF], we demonstrate the influence of the control discretization on the convergence behavior by means of several examples. In Section 5 the discrete adjoint system for the symplectic Galerkin method is derived. It turns out that the reformulation of the variational scheme in Section 3 simplifies the analysis. Whereas commutation of discretization and dualization for symplectic partitioned Runge-Kutta methods has already been shown in [START_REF] Ober-Blöbaum | Discrete mechanics and optimal control: an analysis[END_REF], we prove the same commutation property for the symplectic Galerkin discretization (Theorem 5.2), which is the second main result of this work. In contrast to the discretization with Legendre pseudospectral methods or classical Runge-Kutta methods, no additional closure conditions (see [START_REF] Gong | Connections between the covector mapping theorem and convergence of pseudospectral methods for optimal control[END_REF]) or conditions on the Runge-Kutta coefficients (see [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF]) are required to satisfy the CMP, respectively. Furthermore, using the high order variational integrators presented here, not only the order but also the discretization scheme itself is preserved, i.e. one yields the same schemes for the state and the adjoint system. We conclude with a summary of the results and an outlook for future work in Section 6.

2.

Optimal control for mechanical systems 2.1. Lagrangian dynamics. One of the main subjects of Geometric Mechanics is the study of dynamical systems governed by a Lagrangian. Typically, one considers a mechanical system with configuration manifold Q ⊆ R n together with a Lagrangian function L : T Q → R, where the associated state space T Q describes the position q and velocity q of a particle moving in the system. Usually, the Lagrangian takes the form of kinetic minus potential energy, L(q, q) = K(q, q) -V (q) = 1 2 qT • M (q) • q -V (q), for some (positive definite) mass matrix M (q). A consequence of the principle of least action, also known as Hamilton's principle, establishes that the natural motions q : [0, T ] → Q of the system are characterized by stationary solutions of the action, thus, motions satisfying (1) δ T 0 L(q(t), q(t)) dt = 0 for zero initial and final variations δq(0) = δq(T ) = 0. The resulting equations of motion are the Euler-Lagrange equations (refer to [START_REF] Abraham | Foundations of mechanics[END_REF]),

(2)

d dt ∂L ∂ q - ∂L ∂q = 0 .
When the Lagrangian is regular, that is when the velocity Hessian matrix ∂ 2 L/∂ q2 is nondegenerate, the Lagrangian induces a well defined map, the Lagrangian flow, F t L : T Q → T Q by F t L (q 0 , q0 ) := (q(t), q(t)), where q ∈ C 2 ([0, T ], Q) is the unique solution of the Euler-Lagrange equation (2) with initial condition (q 0 , q0 ) ∈ T Q. By means of the Legendre transform leg L : (q, q) ∈ T Q → (q, p = ∂L ∂ q | (q, q) ) ∈ T * Q, where T * Q is the phase space of positions q and momenta p, one may transform the Lagrangian flow into the Hamiltonian flow F t H : T * Q → T * Q defined by F t H (q 0 , p 0 ) := leg L (q(t), q(t)). Moreover, different preservation laws are present in these systems. For instance the Hamiltonian flow preserves the natural symplectic structure of T * Q and the total energy of the system, typically H(q, p) = K(q, p) + V (q) = 1 2 p T • M (q) -1 • p -V (q) (here K still denotes the kinetic energy, but depending on p rather than on q). Also, if the Lagrangian possess Lie group symmetries, then Noether's theorem asserts that the associated momentum maps are conserved, like for instance the linear momentum and/or the angular momentum.

If external (non conservative) forces F : (q, q) ∈ T Q → (q, F (q, q)) ∈ T * Q are present in the system, Hamilton's principle (1) is replaced by the Lagrange-d'Alembert principle seeking for curves that satisfy the relation

(3) δ T 0 L(q, q) dt + T 0 F (q, q) • δq dt = 0
for zero boundary variations δq(0) = δq(T ) = 0, where the second term is denoted as virtual work. This principle leads to the forced Euler-Lagrange equations

(4) d dt ∂L ∂ q - ∂L ∂q = F (q, q) .
The forced version of Noether's theorem (see e.g. [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF]) states that if the force acts orthogonal to the symmetry action, then momentum maps are still preserved by the flow. Otherwise, the change in momentum maps and energy is determined by the amount of forcing in the system.

2.2. Optimal control problem. Since we are interested in optimally controlling Lagrangian systems, we assume that the mechanical system may be driven by means of some time dependent control parameter u : [0, T ] → U with U ⊂ R m being the control set. Typically, the control appears as an extra variable in the external force such that in the following we consider forces of the form F : T Q × U → T * Q and we replace the right-hand side of (4) by the control dependent force term F (q, q, u).

An optimal control problem for a mechanical system reads (also denoted as Lagranigan optimal control problem in [START_REF] Ober-Blöbaum | Discrete mechanics and optimal control: an analysis[END_REF]) Problem 2.1 (Lagrangian optimal control problem (LOCP)).

(5a) min q, q,u J(q, q, u) = T 0 C(q(t), q(t), u(t)) dt + Φ(q(T ), q(T ))

subject to δ T 0 L(q(t), q(t)) dt + T 0 F (q(t), q(t), u(t)) • δq(t) dt = 0, (5b) (q(0), q(0)) = (q 0 , q0 ), (5c)

with minimization over q ∈ C 1,1 ([0, T ], Q) = W 2,∞ ([0, T ], Q), q ∈ W 1,∞ ([0, T ], T q Q) and u ∈ L ∞ ([0, T ], U ).
The interval length T may either be fixed, or appear as degree of freedom in the optimization problem. Since any optimal control problem with free final time can easily be transformed into a problem with fixed final time (see e.g. [START_REF] Gerdts | Optimal control of ODEs and DAEs[END_REF]), we assume the time T to be fixed from now on. The control set U ⊂ R m is assumed to be closed and convex, and the density cost function C : T Q × U → R and the final cost function Φ : T Q → R n are continuously differentiable, being Φ moreover bounded from below.

Henceforth we should assume that the Lagrangian is regular, i.e. there is a (local) one-to-one correspondence between the velocity q and the momentum p via the Legendre transform and its inverse p = ∂L ∂ q (q, q) and q = ∂L ∂ q -1 (q, p).

Thus, the forced Euler-Lagrange equations ( 4) can be transformed into the partitioned system (6) q(t) = f (q(t), p(t)) , ṗ(t) = g(q(t), p(t), u(t)) with ( 7)

f (q, p) = ∂L ∂ q -1
(q, p) and g(q, p, u) = ∂L ∂q (q, f (q, p)) + F (q, f (q, p), u) .

With some abuse of notation we denote the force and the cost functions defined on T * Q × U and T * Q, respectively, by F (q, p, u) := F (q, f (q, p), u), C(q, p, u) := C(q, f (q, p), u) and Φ(q, p) := Φ(q, f (q, p)) such that Problem 2.1 can be formulated as an optimal control problem for the partitioned system (6).

Problem 2.2 (Optimal control problem (OCP)).

(8a) min q,p,u J(q, p, u) = T 0 C(q(t), p(t), u(t)) dt + Φ(q(T ), p(T ))

subject to q(t) = f (q(t), p(t)) , q(0) = q 0 , (8b) ṗ(t) = g(q(t), p(t), u(t)) , p(0) = p 0 , (8c)

with minimization over q ∈ W 2,∞ ([0, T ], Q), p ∈ W 1,∞ ([0, T ], T * q Q) and u ∈ L ∞ ([0, T ], U ) and the functions f : T * Q → R n , g : T * Q × U → R n are assumed to be Lipschitz continuous.
The first order necessary optimality conditions can be derived by means of the Hamiltonian for the optimal control problem given by ( 9)

H(q, p, u, λ, ψ, ρ 0 ) = ρ 0 C(q, p, u) + λ • f (q, p) + ψ • g(q, p, u)
with ρ 0 ∈ R and λ and ψ are covectors in R n .

Theorem 2.3 (Minimum Principle, e.g. [START_REF] Gerdts | Optimal control of ODEs and DAEs[END_REF]).

Let (q * , p * , u * ) ∈ W 2,∞ ([0, T ], Q)×W 1,∞ ([0, T ], T * q * Q)× L ∞ ([0, T ], U ) be an optimal solution to Problem 2.2. Then there exist functions λ ∈ W 1,∞ ([0, T ], R n ) and ψ ∈ W 1,∞ ([0, T ], R n ) and a constant ρ 0 ≥ 0 satisfying (ρ 0 , λ, ψ) = (0, 0, 0) for all t ∈ [0, T ] such that (10a)
H(q * (t), p * (t), u * (t), λ(t), ψ(t), ρ 0 ) = min u∈U H(q(t), p(t), u, λ(t), ψ(t), ρ 0 ) , for t ∈ [0, T ], and (ρ 0 , λ, ψ) solves the following initial value problem:

λ = -∇ q H(q * , p * , u * , λ, ψ, ρ 0 ), λ(T ) = ρ 0 ∇ q Φ(q * (T ), p * (T )), (10b)
ψ = -∇ p H(q * , p * , u * , λ, ψ, ρ 0 ), ψ(T ) = ρ 0 ∇ p Φ(q * (T ), p * (T )). (10c)
The vectors λ(t) and ψ(t) are the costate or the adjoint variables of the Hamiltonian equations of optimal control. The scalar ρ 0 is called the abnormal multiplier. In the abnormal case, it holds ρ 0 = 0, and otherwise the multiplier can be normalized to ρ 0 = 1. Since no final constraint on the state is present in the optimal control problem, the above principle holds true with ρ 0 = 1 (as proved for instance in [START_REF]Contrôle optimal[END_REF]).

Remark 2.4. If g is affine w.r.t. u, then the topologies can be taken as L 2 for the controls, H 2 = W 2,2 on q and H 1 = W 1,2 on p, and the PMP would still be valid for these classes. Besides, optimal control problems where the optimal control is in L 2 but not in L ∞ are very seldom. For instance, if one is able to express u in function of (q, p, λ, ψ), as for the assumptions of Theorem 5.2, then u is clearly in L ∞ .

Discretization

Since we are interested in solving optimal control problems by some kind of direct method, a discrete approximation of Problem 2.2 is required. To this end, we first introduce two different variational integrators that we employ for the approximation of the control system given in (8b)-(8c). Based on these discrete schemes, we derive the discrete approximations of the optimal control problem that can be solved by standard numerical optimization methods. The controls play no role in the derivations of the variational integrators, therefore we will omit temporarily the dependence of the external force F on u, which will lighten the notation. The discrete schemes including the approximation of the controls are given in Section 3.3.

Discrete Mechanics and Variational

Integrators. Discrete Mechanics is, roughly speaking, a discretization of Geometric Mechanics theory. As a result, one obtains a set of discrete equations corresponding to the Euler-Lagrange equation (4) above but, instead of a direct discretization of the ODE, the latter are derived from a discretization of the base objects of the theory, the state space T Q, the Lagrangian L, etc. In fact, one seeks for a sequence {(t 0 , q 0 ), (t 1 , q 1 ), . . . , (t N , q N )} that approximates the actual trajectory q(t) of the system (q k ≈ q(t k )), for a constant time-step h = t k+1t k > 0.

A variational integrator is an iterative rule that outputs this sequence and it is derived in an analogous manner to the continuous framework. Given a discrete Lagrangian L

d : Q × Q → R and discrete forces F ± d : Q × Q → T * Q,
which are in principle thought to approximate the continuous Lagrangian action and the virtual work, respectively, over a short time

L d (q k , q k+1 ) ≈ t k+1 t k L(q(t), q(t))dt , (11a) F - d (q k , q k+1 ) • δq k + F + d (q k , q k+1 ) • δq k+1 ≈ t k+1 t k F (q(t), q(t)) • δq(t)dt , (11b)
one applies a variational principle to derive the well-known forced discrete Euler-Lagrange (DEL) equation, ( 12)

D 1 L d (q k , q k+1 ) + D 2 L d (q k-1 , q k ) + F - d (q k , q k+1 ) + F + d (q k-1 , q k ) = 0 , for k = 1, . . . , N -1
, where D i stands for the partial derivative with respect to the i-th component. The equation defines an integration rule of the type (q k-1 , q k ) → (q k , q k+1 ), however if we define the pre-and post-momenta (also denoted as discrete Legendre transforms)

p - k := -D 1 L d (q k , q k+1 ) -F - d (q k , q k+1 ), k = 0, . . . , N -1, and (13a)
p + k := D 2 L d (q k-1 , q k ) + F + d (q k-1 , q k ), k = 1, . . . , N, (13b)
the discrete Euler-Lagrange equation ( 12) is read as the momentum matching p - k = p + k =: p k and defines an integration rule of the type (q k , p k ) → (q k+1 , p k+1 ).

The nice part of the story is that the integrators derived in this way naturally preserve (or nearly preserve) the quantities that are preserved in the continuous framework, the symplectic form, the total energy (for conservative systems) and, in presence of symmetries, the linear and/or angular momentum (for more details, see [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF]). Furthermore, other aspects of the continuous theory can be "easily" adapted, symmetry reduction [START_REF] Campos | Discrete Lagrange-d'Alembert-Poincaré equations for Euler's disk[END_REF][START_REF] Colombo | Discrete second-order Euler-Poincaré equations. Applications to optimal control[END_REF][START_REF] Iglesias | Discrete nonholonomic Lagrangian systems on Lie groupoids[END_REF], constraints [START_REF] Johnson | Dangers of two-point holonomic constraints for variational integrators[END_REF][START_REF] Kobilarov | Geometric discretization of nonholonomic systems with symmetries[END_REF], control forces [START_REF] Campos | Higher order variational time discretization of optimal control problems[END_REF][START_REF] Ober-Blöbaum | Discrete mechanics and optimal control: an analysis[END_REF], etc.

High Order Variational Integrators. High order variational integrators for time dependent or independent systems (HOVI[t]

) are a class of integrators that, by using a multi-stage approach, aim at a high order accuracy on the computation of the natural trajectories of a mechanical system while preserving some intrinsic properties of such systems. In particular, symplectic-partitioned Runge-Kutta methods (spRK) and, what we call here, symplectic Galerkin methods (sG) are s-stage variational integrators of order up to 2s.

The derivation of these methods follows the general scheme that comes next, the specifics of each particular case are detailed in the following subsections. For a fixed time step h, one considers a series of points q k , refereed as macro-nodes. Between each couple of macro-nodes (q k , q k+1 ), one also considers a set of micro-data, the s stages: For the particular cases of sG and spRK methods, we consider micro-nodes Q 1 , . . . , Q s and micro-velocities Q1 , . . . , Qs , respectively. Both macro-nodes and micro-data (micro-nodes or micro-velocities) are required to satisfy a variational principle, giving rise to a set of equations, which properly combined, define the final integrator.

Here and after, we will use the following notation: Let 0 ≤ c 1 < . . . < c s ≤ 1 denote a set of collocation points and consider the associated Lagrange polynomials and nodal weights, that is,

l j (t) := i =j t -c i c j -c i and b j := 1 0 l j (t)dt ,
respectively. Note that the pair of (c i , b i )'s define a quadrature rule and that, for appropriate c i 's, this rule may be a Gaussian-like quadrature, for instance, Gauss-Legendre, Gauss-Lobatto, Radau or Chebyshev. Now, for the sake of simplicity and independently of the method, we will use the same notation for the nodal coefficients. We define for spRK and sG, respectively, ( 14)

a ij := ci 0 l j (t)dt and a ij := dl j dt ci .
Moreover, for spRK, we will also use the nodal weights and coefficients ( bj , āij ) given by Equation ( 22) and, for sG, the source and target coefficients α j := l j (0) and β j := l j (1) .

Finally, if L denotes a Lagrangian from R n × R n to R coupled with an external force F : (q, q) ∈ R n × R n → (q, F (q, q)) ∈ R n × R n , then we define

P i := ∂L ∂ q i = ∂L ∂ q (Qi, Qi) and Ṗi := ∂L ∂q i + F i = ∂L ∂q (Qi, Qi) + F (Q i , Qi ) ,
where (Q i , Qi ) are couples of micro-nodes and micro-velocities given by each method. Besides, D i will stand for the partial derivative with respect to the i-th component.

3.2.1. Symplectic-Partitioned Runge-Kutta Methods. Although the variational derivation of spRK methods in the framework of Geometric Mechanics is an already known fact (see [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF] for an "intrinsic" derivation, as the current, or [START_REF] Hairer | Geometric numerical integration[END_REF] for a "constrained" one), both based on the original works of [START_REF] Suris | Hamiltonian methods of Runge-Kutta type and their variational interpretation[END_REF][START_REF] Sanz-Serna | Numerical Hamiltonian problems[END_REF], we present it here again in order to ease the understanding of and the comparison with sG methods below. Given a point q 0 ∈ R n and vectors { Qi } i=1,...,s ⊂ R n , we define the polynomial curves

Q(t) := s j=1 l j (t/h) Qj and Q(t) := q 0 + h s j=1 t/h 0 l j (τ )dτ Qj .
We have

(15) Qi = Q(h • c i ) and Q i := Q(h • c i ) = q 0 + h s j=1 a ij Qj .
Note that the polynomial curve Q is uniquely determined by q 0 and { Qi } i=1,...,s . In fact, it is the unique polynomial curve Q of degree s such that Q(0) = q 0 and Q(h • c i ) = Qi . However, if we define the configuration point ( 16)

q 1 := Q(h • 1) = q 0 + h s j=1 b j Qj
and consider it fixed, then Q is uniquely determined by q 0 , q 1 and the Qi 's but one. Namely, take any 1 ≤ i 0 ≤ s such that b i0 = 0 and fix it, then

Qi0 =   q 1 -q 0 h - j =i0 b j Qj   /b i0 .
We now define the multi-vector discrete Lagrangian

(17) L d ( Q1 , . . . , Qs ) := h s i=1 b i L(Q i , Qi )
and the multi-vector discrete force

F d ( Q1 , . . . , Qs ) • (δQ 1 , . . . , δQ s ) := h s i=1 b i F (Q i , Qi ) • δQ i .
Although not explicitly stated, they both depend also on q 0 . If we write the micro-node variations δQ i in terms of the micro-velocity variations δ Qi (by definition ( 15)), we have that the multi-vector discrete force is

F d ( Q1 , . . . , Qs ) • (δQ 1 , . . . , δQ s ) = h 2 s j=1 s i=1 b i a ij F (Q i , Qi ) • δ Qj .
The two-point discrete Lagrangian is then (18) L d (q 0 , q 1 ) := ext

P s ([0,h],R n ,q0,q1) L d ( Q1 , . . . , Qs )
where P s ([0, h], R n , q 0 , q 1 ) is the space of polynomials Q of order s from [0, h] to R n such that Q(0) = q 0 and Q(h) = q 1 and the vectors Qi 's determine such polynomials as discussed above. The so called "extremal" is realized by a polynomial Q ∈ P s ([0, h], R n , q 0 , q 1 ) such that for any variations (δ Q1 , . . . , δ Qs ), taking into account that δq 0 = δq 1 = 0 and that δ Qi0 = j =i0 ∂ Qi0 /∂ Qj δ Qj . For convenience, the previous equation ( 19) is developed afterwards. The two-point discrete forward and backward forces are then (20)

F ± d (q 0 , q 1 ) • δ(q 0 , q 1 ) := h s i=1 b i F (Q i , Qi ) • ∂Q i ∂q ± δq ± ,
where q -= q 0 and q + = q 1 . Using the previous relations, we may write

F - d = h s i=1 b i (1 -a ii0 /b i0 )F i and F + d = h s i=1 b i a ii0 /b i0 F i .
By the momenta-matching rule ( 13), we have that

-p 0 = -D i0 L d ( Q1 , . . . , Qs )/(hb i0 ) + D q0 L d ( Q1 , . . . , Qs ) + F - d , p 1 = D i0 L d ( Q1 , . . . , Qs )/(hb i0 ) + F + d .
where D q0 stands for the partial derivative with respect to q 0 . Combining both equations, we obtain that

D i0 L d + h 2 s i=1 b i a ii0 F i = hb i0 p 1 and p 1 = p 0 + D q0 L d + h s i=1 b i F i .
Coming back to Equation ( 19), we have that

0 = δL d ( Q1 , . . . , Qs ) • (δ Q1 , . . . , δ Qs ) + F d ( Q1 , . . . , Qs ) • (δQ 1 , . . . , δQ s ) = j =i0 D j L d + h 2 s i=1 b i a ij F i + ∂ Qi0 ∂ Qj D i0 L d + h 2 s i=1 b i a ii0 F i δ Qj .
Therefore, for j = i 0 , we have that

D j L d + h 2 s i=1 b i a ij F i = b j /b i0 • D i0 L d + h 2 s i=1 b i a ii0 F i .
Thus, the integrator is defined by

D j L d ( Q1 , . . . , Qs ) + h 2 s i=1 b i a ij F i = hb j p 1 , (21a)
q 1 = q 0 + h s j=1 b j Qj , (21b)
p 1 = p 0 + D q0 L d ( Q1 , . . . , Qs ) + h s i=1 b i F i . (21c)
Besides, using the definition of the discrete Lagrangian, we have

D j L d ( Q1 , . . . , Qs ) + h 2 s i=1 b i a ij F i = h 2 s i=1 b i a ij Ṗi + hb j P j , D q0 L d ( Q1 , . . . , Qs ) + h s i=1 b i F i = h s i=1 b i Ṗi .
Therefore, we may write

P j = p 0 + h s i=1 b i (1 -a ij /b j ) Ṗi = p 0 + h s i=1 āji Ṗi , p 1 = p 0 + h s i=1 b i Ṗi = p 0 + h s i=1 bi Ṗi ,
were āij and bi are given by [START_REF] Iglesias | Discrete nonholonomic Lagrangian systems on Lie groupoids[END_REF] b

i āij + bj a ji = b i bj , b i = bi .
In summary, the equations that define the spRK integrator (with forces), are together with ( 22)

q 1 =q 0 + h s j=1 b j Qj , p 1 =p 0 + h s j=1 bj Ṗj , (23a)
Q i =q 0 + h s j=1 a ij Qj , P i =p 0 + h s j=1 āij Ṗj , (23b)
P i = ∂L ∂ q (Q i , Qi ) , Ṗi = ∂L ∂q (Q i , Qi ) + F (Q i , Qi ) . (23c)
3.2.2. Symplectic Galerkin Methods. Galerkin methods are a class of methods to transform a problem given by a continuous operator (such as a differential operator) to a discrete problem. As such, spRK methods fall into the scope of this technique and could be also classified as "symplectic Galerkin" methods. However, we want to stress on the difference between what is called spRK in the literature and what we refer here as sG. The wording should not be confused by the one used in [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF].

Given points {Q i } i=1,...,s ⊂ R n , we define the polynomial curves

Q(t) := s j=1 l j (t/h)Q j and Q(t) := 1 h s j=1 lj (t/h)Q j .
We have

Q i = Q(h • c i ) and Qi := Q(h • c i ) = 1 h s j=1 a ij Q j .
Note that the polynomial curve Q is uniquely determined by the points {Q i } i=1,...,s . In fact, it is the unique polynomial curve Q of degree s -

1 such that Q(h • c i ) = Q i .
However, if we define the configuration points (24)

q 0 := Q(h • 0) = s j=1 α j Q j and q 1 := Q(h • 1) = s j=1 β j Q j
and consider them fixed, then Q is uniquely determined by q 0 , q 1 and the Q i 's but a couple. For instance, we may consider Q 1 and Q s as functions of the others, since the relations (24) define a system of linear equations where the coefficient matrix has determinant γ := α 1 β sα s β 1 = 0 (if and only if c 1 = c s ). More precisely,

Q 1 Q s = 1 γ β s -α s -β 1 α 1 q 0 - s-1 j=2 α j Q j q 1 - s-1 j=2 β j Q j .
We now define the multi-point discrete Lagrangian (25)

L d (Q 1 , . . . , Q s ) := h s i=1 b i L(Q i , Qi )
and the multi-vector discrete force

F d (Q 1 , . . . , Q s ) • (δQ 1 , . . . , δQ s ) := h s i=1 b i F (Q i , Qi ) • δQ i .
The two-point discrete Lagrangian is then (26) L d (q 0 , q 1 ) := ext

P s ([0,h],R n ,q0,q1) L d (Q 1 , . . . , Q s )
where P s ([0, h], R n , q 0 , q 1 ) is the space of polynomials Q of order s from [0, h] to R n such that the points Q i 's determine such polynomials as discussed above. The so called "extremal" is realized by a polynomial Q ∈ P s ([0, h], R n , q 0 , q 1 ) such that

(27) δL d (Q 1 , . . . , Q s ) • (δQ 1 , . . . , δQ s ) + F d (Q 1 , . . . , Q s ) • (δQ 1 , .
. . , δQ s ) = 0 for any variations (δQ 1 , . . . , δQ s ), taking into account that δq 0 = δq 1 = 0 and that δQ i = s-1 j=2 ∂Q i /∂Q j δQ j , i = 1, s. For convenience, the previous equation ( 27) is developed afterwards.

The two-point discrete forward and backward forces are then formally defined by Equation [START_REF] Hall | Spectral variational integrators[END_REF]. Using the previous relations, we may write

F - d = h(β s b 1 F 1 -β 1 b s F s )/γ and F + d = h(α 1 b s F s -α s b 1 F 1 )/γ
. By the momenta-matching rule (13), we have that

-p 0 = β s /γ • (D 1 L d + hb 1 F 1 ) -β 1 /γ • (D s L d + hb s F s ) and p 1 = -α s /γ • (D 1 L d + hb 1 F 1 ) + α 1 /γ • (D s L d + hb s F s ) .
By a linear transformation of both equations, we obtain

D 1 L d (Q 1 , . . . , Q s ) + hb 1 F 1 = -α 1 p 0 + β 1 p 1 and D s L d (Q 1 , . . . , Q s ) + hb s F s = -α s p 0 + β s p 1 .
Coming back to Equation ( 27), we have that

0 = (δL d (Q 1 , . . . , Q s ) + F d (Q 1 , . . . , Q s )) • (δQ 1 , . . . , δQ s ) = s-1 j=2 (D 1 L d + hb 1 F 1 ) ∂Q 1 ∂Q j + (D j L d + hb j F j ) + (D s L d + hb s F s ) ∂Q s ∂Q j δQ j .
Therefore, for j = 2, . . . , s -1, we obtain

γ(D j L d + hb j F j ) = (α j β s -α s β j )(D 1 L d + hb 1 F 1 ) + (α 1 β j -α j β 1 )(D s L d + hb s F s ) = (α 1 β s -α s β 1 )(β j p 1 -α j p 0 ) .
Thus, the integrator is defined by

D j L d (Q 1 , . . . , Q s ) + hb j F j = -α j p 0 + β j p 1 , j = 1, . . . , s; (28a)
q 0 = s j=1 α j Q j and q 1 = s j=1 β j Q j (28b)
Besides, using the definition of the discrete Lagrangian, we have

D j L d (Q 1 , . . . , Q s ) = h s i=1 b i ∂L ∂q i ∂Q i ∂ Qj + ∂L ∂ q i ∂ Qi ∂ Qj D j L d (Q 1 , . . . , Q s ) + hb j F j = hb j Ṗj + s j=1 b i a ij P i .
Therefore, we may simply write

hb j Ṗj + s j=1 b i a ij P i = -α j p 0 + β j p 1 .
In summary and for a proper comparison, we write the equations that define the sG integrator (with forces) in a pRK way, that is

q 0 = s j=1 α j Q j , q 1 = s j=1 β j Q j , (29a) Qi = 1 h s j=1 a ij Q j , Ṗi = β i p 1 -α i p 0 h bi + 1 h s j=1 āij P j , (29b)
P i = ∂L ∂ q (Q i , Qi ) , Ṗi = ∂L ∂q (Q i , Qi ) + F (Q i , Qi ) , (29c)
where b i a ij + bj āji = 0 and b i = bj .

We remark that Equation (28a) generalizes the ones obtained in [START_REF] Campos | Higher order variational time discretization of optimal control problems[END_REF][START_REF] Leok | Foundations of computational geometric mechanics[END_REF], where the collocation points are chosen such that c 1 = 0 and c s = 1, which is a rather particular case.

3.2.3. Similarities and differences between spRK and sG. As already mentioned, both methods can be considered of Galerkin type. In this sense, spRK and sG could be refereed as a symplectic Galerkin integrators of 1st and 0th kind, respectively, since spRK is derived from the 1st derivative of an extremal polynomial and sG from the polynomial itself. At this point, a very natural question could arise: Are spRK and sG actually two different integrator schemes? Even though the derivations of both methods are quite similar, they are in general different (although they could coincide for particular choices of the Lagrangian, the collocation points and the integral quadrature). A weak but still fair argument to support this is that, at each step, spRK relies on the determination of the micro-velocities Qi , while sG does so on the micro-nodes Q i . All the other "unknowns" are then computed from the determined micro-data.

In the simplest of the cases, that is, the case where one considers a Lagrangian of the form kinetic minus potential energy, L(q, q) = 1 2 qT M q -U (q), with M a constant mass matrix; s = 2 micro-nodes (inner-stages); and Lobatto's quadrature, c 1 = 0, c 2 = 1; one may show that both schemes, spRK [START_REF] Johnson | Dangers of two-point holonomic constraints for variational integrators[END_REF] and sG [START_REF]An overview of variational integrators[END_REF], reduce to the well-known leap-frog or Verlet method. They will differ when the previous main assumptions are violated, for instance if M is not constant or the quadrature is other than Lobatto's.

Example 3.1. We consider a Lagrangian with a scalar mass matrix dependent on the configuration, that is, a Lagrangian of the form L(q, q) = 1 2 λ(q) q 2 -V (q), with λ : Q → R. Under this assumption and noting λ 1/2 := λ0+λ1 2 , (∇)λ i := (∇)λ(q i ), (∇)V i := (∇)V (q i ), i = 0, 1, the spRK scheme [START_REF] Johnson | Dangers of two-point holonomic constraints for variational integrators[END_REF] as well as the the sG scheme (29) reduce to

p 1/2 = p 0 + h 2 ∇λ 0 2λ 2 a p 1/2 2 -∇V 0 , q 1 = q 0 + h 2 1 λ a + 1 λ b p 1/2 , p 1 = p 1/2 + h 2 ∇λ 1 2λ 2 b p 1/2 2 -∇V 1 ,
with a slight difference in the subindexes appearing in the framed factors. While in the spRK scheme, a = 0 and b = 1; in the sG scheme a = b = 1/2. It is important to note that, even though the difference is small, it makes both schemes certainly different. Besides notice that the first two equations define p 1/2 and q 1 implicitly and that the whole set reduces to the Verlet method for a constant λ. Indeed, it is shown in [START_REF] Campos | High order variational integrators: a polynomial approach[END_REF][START_REF] Ober-Blöbaum | Galerkin variational integrators and modified symplectic Runge-Kutta methods[END_REF] that for a Lagrangian with constant mass matrix and Lobatto quadrature rule, the sG and the spRK method coincide.

3.2.4. Order of the schemes. With respect to the accuracy of the schemes, for any Gaussian-like quadrature (Gauss-Legendre, Gauss-Lobatto, Radau and Chebyshev) and any method (spRK and sG), the schemes have a convergence order up to 2s (which is only attained by the combination Gauss-Lobatto together with spRK) but no lower than 2s-2, being s the number of internal stages, see Table 1. We emphasize that these orders have been determined numerically experimenting with several "toy examples" for which exact solutions are known, e.g. the harmonic oscillator and the 2 body problem (see Section 5), however they coincide with the known analytic results when available, that is spRK together with the Gauss-Legendre or Gauss-Lobatto quadratures (see e.g.

[19]).

3.3.

Discrete optimal control problem. For the discretization of the optimal control problem 2.1, we employ the class of high order variational integrators. By choosing an appropriate approximation J d of the cost functional J, the general discrete Lagrangian optimal control problem as discretization of the Langrangian optimal control problem 2.1 reads (see also [START_REF] Ober-Blöbaum | Discrete mechanics and optimal control: an analysis[END_REF]) Problem 3.2 (Discrete Lagrangian optimal control problem).

(30a) min subject to

{q k ,u k } N k=0 J d ({q k , u k } N k=0 ) spRK sG micro-data Qi Q i polynomial degree s s -
q 0 = q 0 , (30b) D 2 L(q 0 , q0 ) + D 1 L d (q 0 , q 1 ) + F - 0 = 0, (30c) D 2 L d (q k-1 , q k ) + D 1 L d (q k , q k+1 ) + F + d (q k-1 , q k ) + F - d (q k , q k+1 ) = 0, (30d)
for k = 1, . . . , N -1 and where Equation (30d) is the forced discrete Euler-Lagrange equation defined in [START_REF] Dontchev | The Euler approximation in state constrained optimal control[END_REF] and Equations (30b)-(30c) correspond to the initial condition (5c) expressed by means of the discrete Legendre transform [START_REF] Dontchev | Second-order Runge-Kutta approximations in control constrained optimal control[END_REF]. Here, the control trajectory u is approximated by the discrete values u k , k = 0, . . . , N , such that u k ≈ u(t k ). Note that for the controlled case, the F ± d are dependent on u k . To specify the discrete optimal control problem, in particular, the approximation of the control parameter u and the cost functional J, we focus on the high order variational integrators derived in Section 3.2, namely the spRK and the sG method, and find a discrete version of the more general optimal control problem 2.2.

As for the approximation of q(t) and q(t), we also use a polynomial for the approximation of the control function u(t) on [0, h]. For a given set of collocation points 0 ≤ c 1 < . . . < c s ≤ 1 and given control points {U i } i=1,...,s ⊂ U we define the polynomial of degree s -1

U(t) := s j=1 l j (t/h)U j such that U i = U(h • c i ), i = 1, . . . , s.
Note that the control polynomial U(t) has the same degree as the polynomial Q(t) for the sG integrator, whereas for the spRK integrator it coincides with the polynomial degree of Q(t). To take in consideration the control dependent force in the previous derivation of the spRK and the sG schemes into account, we replace in the definitions for F ± d in Equation [START_REF] Hall | Spectral variational integrators[END_REF] the external force

F i = F (Q i , Qi ) by F i = F (Q i , Qi , U i ).
Furthermore, for a regular Lagrangian and by using the definitions for f and g in ( 7), we can write Equations (23c) and (29c) as

f (Q i , P i ) = ∂L ∂ q -1 (Q i , P i ) , g(Q i , P i , U i ) = ∂L ∂q (Q i , Qi ) + F (Q i , Qi , U i )
such that the spRK scheme can be written as

q 1 =q 0 + h s j=1 b j f (Q j , P j ) , p 1 =p 0 + h s j=1 bj g(Q j , P j , U j ) , (31a)
Q i =q 0 + h s j=1 a ij f (Q j , P j ) , P i =p 0 + h s j=1 āij g(Q j , P j , U j ) , (31b)
with b i āij + bj a ji = b i bj and b i = bj and the sG scheme reduces to

q 0 = s j=1 α j Q j , q 1 = s j=1 β j Q j , (32a) f (Q i , P i ) = 1 h s j=1 a ij Q j , g(Q i , P i , U i ) = β i p 1 -α i p 0 h bi + 1 h s j=1 āij P j , (32b)
where b i a ij + bj āji = 0 and b i = bj . Remember that the coefficients a ij are different for the two schemes [START_REF] Leyendecker | Variational integrators for constrained dynamical systems[END_REF] and [START_REF] Leyendecker | A variational approach to multirate integration for constrained systems[END_REF] (see ( 14)). To approximate the cost functional h 0 C(q(t), p(t), u(t)) dt in (8a) we employ the same quadrature rule that we use to approximate the action on [0, h] (cf. ( 17) and ( 25)) such that the discrete density cost function C d is defined by

C d ({Q i , P i , U i } s i=1 ) := h s i=1 b i C(Q i , P i , U i ) ≈ h 0 C(q(t), p(t), u(t)) dt.
So as to prevent a proliferation of symbols and alleviate the notation, along a time step interval [t k , t k+1 ], we write q k h , p k h and u k h instead of {q k , {Q k i } s i=1 , q k+1 }, {p k , {P k i } s i=1 , p k+1 } and {U k i } s i=1 , respectively. We drop the superscript k if we consider an arbitrary time step interval [0, h]. With some abuse, along the whole interval [0, T ], we equally write q h , p h and u h instead of {{q k , Q k i } k=0,...,N -1 i=1,...,s , q N }, {{p k , P k i } k=0,...,N -1 i=1,...,s , p N } and {U k i } k=0,...,N -1 i=1,...,s

, respectively. With this notation we define the discrete cost function J d as

J d (q h , p h , u h ) := N -1 k=0 C d (q k h , p k h , u k h ) + Φ(q N , p N )
and introduce the following two discrete optimal control problems, where the discretization in Problem 3.3 is based on the spRK integrator and in Problem 3.4 on the sG integrator.

Problem 3.3 (Discrete optimal control problem: the spRK case). (33a) min

q h ,p h ,u h J d (q h , p h , u h ) subject to q k+1 =q k + h s j=1 b j f (Q k j , P k j ) , p k+1 =p k + h s j=1 bj g(Q k j , P k j , U k j ) , (33b)
Q k i =q k + h s j=1 a ij f (Q k j , P k j ) , P k i =p k + h s j=1 āij g(Q k j , P k j , U k j ) (33c) k = 0, . . . , N -1, i = 1, . . . , s, with b i āij + bj a ji = b i bj and b i = bj , (q 0 , p 0 ) = (q 0 , p 0 ), U k i ∈ U. (33d)
Problem 3.4 (Discrete optimal control problem: the sG case). (34a) min

q h ,p h ,u h J d (q h , p h , u h ) subject to q k = s j=1 α j Q k j , q k+1 = s j=1 β j Q k j , (34b) f (Q k i , P k i ) = 1 h s j=1 a ij Q k j , g(Q k i , P k i , U k i ) = β i p k+1 -α i p k h bi + 1 h s j=1 āij P k j , (34c)
k = 0, . . . , N -1, i = 1, . . . , s, with b i a ij + bj āji = 0 and b i = bj , (q 0 , p 0 ) = (q 0 , p 0 ), U k i ∈ U. (34d) Since Problem 3.3 has been extensively studied in [START_REF] Ober-Blöbaum | Discrete mechanics and optimal control: an analysis[END_REF] (as discussed in Section 3.4), in this work we focus on Problem 3.4.

Comparison of different solution methods.

In Figure 1 (see also [START_REF] Ober-Blöbaum | Discrete mechanics and optimal control: an analysis[END_REF]) we present schematically different discretization strategies for optimal control problems. Starting with the Lagrangian optimal control problem 2.1, we obtain via variation (for the derivation of the Euler-Lagrange equations) the optimal control problem 2.2. For its solution, direct or indirect methods can be employed (the differences of direct and indirect methods are already discussed in Section 1).

In the DMOC approach, rather than discretizing the differential equations arising from the Lagrange-d'Alembert principle, we discretize in the earliest stage, namely already at the level of the variational principle. Then, we perform the variation only on the discrete level which results in a nonlinear programming problem (in particular we obtain the discrete Lagrangian optimal control problem 3.2). Its necessary optimality conditions are derived by a dualization step as for a standard direct method. This approach that uses the concept of discrete mechanics leads to a special discretization of the system equations based on variational integrators. Thus, the discrete optimal control problem inherits special properties exhibited by variational integrators as extensively discussed in [START_REF] Ober-Blöbaum | Discrete mechanics and optimal control: an analysis[END_REF]. In this work we are interested in the question under which conditions the discretized necessary optimality conditions (discrete PMP) and the KKT conditions resulting from the discrete optimal control problems (3.3) and (3.4) ( KKT) are identical. To this end, we summarize by now already known equivalence relations (A), (B) and (C) as indicated in Figure 1. Equivalence (A). This equivalence corresponds to the commutation of variation and discretization. For particular variational integrators their equivalence to other well-known integration methods has been shown, e.g. the equivalence to the Störmer-Verlet method ( [START_REF] Mackay | Some aspects of the dynamics and numerics of Hamiltonian systems[END_REF]), to the Newmark algorithm ( [START_REF] Kane | Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems[END_REF]) or, more generally, to spRK methods ( [START_REF] Suris | Hamiltonian methods of Runge-Kutta type and their variational interpretation[END_REF][START_REF] Marsden | Discrete mechanics and variational integrators[END_REF][START_REF] Hairer | Geometric numerical integration[END_REF]) applied to the corresponding Hamiltonian system. Whereas the variational derivation of spRK methods is an already known fact and was presented in Section 3.2.1 from a slightly different point of view, we found in Section 3.2.2 a new class of integrators, the sG methods, that applied to the Hamiltonian systems are equivalent to variational integrators based on the discrete Lagrangian [START_REF] Leok | Foundations of computational geometric mechanics[END_REF]. Obviously, if the variation of the discrete Lagrange-d'Alembert principle (the variational integrator) and the discretization of the Euler-Lagrange or, equivalently, the Hamiltonian equations is the same, then, using the same discretization for the cost functional provides the same NLP in the middle and the right branch in Figure 1.

Equivalence (B). Of course, if NLP and NLP are identical (equivalence (A)), then also the KKT and the KKT conditions are identical. Equivalence (C). This equivalence corresponds to the commutation of discretization and dualization. As already summarized in Section 1, for the class of Legendre pseudospectral methods the commutation property has been proved if additional closure conditions are satisfied (see [START_REF] Gong | Connections between the covector mapping theorem and convergence of pseudospectral methods for optimal control[END_REF][START_REF] Ross | Legendre pseudospectral approximations of optimal control problems[END_REF]). In [START_REF] Ober-Blöbaum | Discrete mechanics and optimal control: an analysis[END_REF] the commutation of discretization and dualization is proved for spRK methods, that means, the dualization of Problem 3.3 is the same as the discretization of the necessary optimality conditions using a spRK method. Due to equivalence (A), this commutation property also holds for high order variational integrators based on the discrete Lagrangian [START_REF]Numerical analysis in optimal control[END_REF]. However, for general classes of variational integrators (in particular there exist high order variational integrators that are not equivalent to RK methods, see e.g. [START_REF] Ober-Blöbaum | Galerkin variational integrators and modified symplectic Runge-Kutta methods[END_REF]), the commutation property is still an open question. As new contribution in this paper, we show in Section 5 that discretization and dualization also commute for the sG integrator (see Theorem 5.2) (i.e. the dualization of Problem 3.4 is the same as the discretization of the necessary optimality conditions using an sG method). We therefore find (besides the spRK methods) another class of variational integrators that fulfills the commutation property.

Analysis of convergence of the primal variables

In this section, we present one of the main results of the paper, that is the convergence of the primal variables. Before that, a couple of comments are necessary to clarify the assumptions and notation that appear in the statement. Examples will enlighten the result.

We say that a function f : H → R is coercive, where H is a Hilbert space with norm • (in our case of study H is either R m or L 2 ([0, T ], R m )), if there exists a scalar factor α > 0 for which

lim inf x →∞ f (x) x 2 ≥ α .
If f depends on a further variable, f = f (x, y), we say that f is uniformly coercive in x (with respect to y) if the coercivity factor α does not depend on y.

In the next result, a discrete trajectory q h , either over a time step interval [0, h] or over the whole time interval [0, T ], should be understood as a continuous trajectory along [0, h] or [0, T ], respectively. To that, on [0, h], q h = {q 0 , {Q i } s i=1 , q 1 } can be also viewed as its own linear interpolation, that is as the piecewise-linear continuous function q h : [0, h] → R n determined by q h (0) = q 0 , q h (i • h/(s + 1)) = Q i , i = 1, . . . , s , q h (h) = q 1 and linear in between. One proceeds similarly on [0, T ] and as well for p h and u h . (H2) F is affine on the controls, i.e. F (q, q, u) = F 0 (q, q) + u • F 1 (q, q); (H3) C is uniformly coercive in u and smooth in (q, p); (H4) (OCP ), the continuous Problem 2.2, has a unique solution (q, p, ū); (H5) b i > 0 for i = 1, . . . , s; and (H6) the associated spRK or sG scheme is convergent (for L, F and any fixed u).

Then (q h , ph , ūh ) converges (up to subsequence) to (q, p, ū) as h → 0 (N → ∞), strongly in (q, p) and weakly in u, where (q h , ph , ūh ) is the solution to (OCP ) h , the corresponding discrete Problem 3.3 or 3.4.

Proof.

The assumption of coercivity H3 on the density cost function C (recall Φ is assumed to be bounded from below) implies the coercivity of the total cost functional J and, together with the assumption of positiveness H5 on the weight coefficients b i , ensures also the coercivity of the discrete density cost C d and of the discrete total cost J d : Indeed, from the definition and for u h big enough,

C d (q h , p h , u h ; h) u h 2 ≥ h min i b i α and J d (q h , p h , u h ; h) u h 2 ≥ T min i b i α ,
where (q h , p h , u h ) represents a one step trajectory on the left while a full trajectory on the right. The coercivity of J follows similarly by integration. We stress the fact that J d is uniformly coercive with respect to the time step h. Besides we have from H1 that, for h small enough, the discrete Lagrangian L d is regular. Therefore for each control set u h , a unique solution (q h , p h ) to the discrete Euler-Lagrange equations exists. The minimization in (OCP ) h is then nothing but a finite dimensional minimization problem whose solution existence is provided by the coercivity of C d . We denote (q h , ph , ūh ) such solution.

Since the discrete cost J d is uniformly coercive with respect to h, the sequence (ū h ), thought in L 2 ([0, T ], U ), is bounded. Hence, up to subsequence, it converges to some control ũ ∈ L 2 ([0, T ], U ), for the weak topology of L 2 . It only remains to show that (q, p, ũ), where (q, p) is the unique solution of the mechanical system corresponding to the control ũ (hypothesis H1), is in fact the same point as (q, p, ū).

Firstly, by H2, the controls enter the continuous and discrete dynamical equations (5b), (23a), (23b), (29a) and (29b) linearly. It follows that (q h , ph ) converges to (q, p) (more technical details on this standard reasoning may be found in [START_REF] Trélat | Some properties of the value function and its level sets for affine control systems with quadratic cost[END_REF]). Secondly, we observe that the control ũ is optimal: Indeed, since the discrete cost J d is an approximation to the continuous one J, we have that for some exponent r ≥ 1

J(q h , ph , ūh ) = J d (q h , ph , ūh ; h) + O(h r ) ≤ J d ((q, p, ū) h ; h) + O(h r ) = J(q, p, ū) + O(h r )
, where (q, p, ū) h , for each h, represents simply the evaluation of (q, p, ū) at the collocation points of each time interval. Passing to the limit, from the lower semi-continuity of J with respect to u (given by the integral expression and hypothesis H3), it follows that J(q, p, ũ) ≤ J(q, p, ū) . Hence (q, p, ũ) is optimal and by uniqueness coincides with (q, p, ū). Remark 4.2. Even though the rather "classical" definition of coercivity used here is perhaps less restrictive than the one given in [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF][START_REF]Numerical analysis in optimal control[END_REF], there is no direct relation between them in the sense that none of them implies or is implied by the other. They have non-empty intersection. In there, the coercivity is a classical sufficient second-order assumption ensuring the absence of conjugate points. In here, the coercivity implies the existence of a (global) optimal control. The framework that 4.1 proposes has some advantages: It permits to prove without difficulty the existence of solutions for the discretized problem and it ensures the convergence from a simple topological argument. Moreover, the proof itself has the potential of being more general, for instance to consider final constraints (by means of finer arguments, the concept of end-point mapping, the general conjugate point theory).

Remark 4.3. It is worth to note that the previous proof withstands some easy generalizations. If we now take more general dynamics (nonlinear force in u) and costs, then the above reasoning works as well provided the cost is coercive in some L p and the dynamics satisfy, for instance,

lim sup u →∞ Ψ(q, p, u) u r = 0 , p > r ,
where Ψ = 0 stands for the dynamical constraints (5b).

Remark 4.4. A convergence proof (including consistency and stability) for general variational integrators (as assumed in H6) is topic of ongoing research. For particular classes, the convergence is proven by showing that the variational integrator is equivalent to another well-known convergent method, as for example for symplectic partitioned Runge-Kutta methods. For a recent convergence analysis for Galerkin variational integrators by means of variational error analysis we refer to [START_REF] Hall | Spectral variational integrators[END_REF]. The assumption H4, the uniqueness of the solution of (OCP ), is a classical one. It can be weakened by stating the result in terms of closure points as follows. We assume that C 0 ([0, T ], T * Q) is endowed with the uniform convergence topology and that L 2 ([0, T ], U ) is endowed with the weak topology. Then every closure point of the family of solutions (q h , ph , ūh )

of (OCP ) h in C 0 ([0, T ], T * Q) × L 2 ([0, T ], U
) is a solution of (OCP ).

Remark 4.5. In the formulation of the previous discrete optimal control problems 3.3 and 3.4, we have chosen to discretize the control parameter and the cost functional in accordance to the dynamical discretization. Other possibilities are available which must be pondered. Let's assume temporarily that, besides the original set of collocation points 0 ≤ c 1 < . . . < c s ≤ 1, we have a couple of extra sets of them: 0 ≤ d 1 < . . . < d r ≤ 1 and 0 ≤ e 1 < . . . < e t ≤ 1, for which U is determined by the former and J is discretized by the quadrature rule associated to the latter. That is, U : [0, h] → R is a polynomial of order r -1 determined by r points Ūi = U(d i h), i = 1, . . . , r, and for which U j := U(c j h), j = 1, . . . , s, and Ûk := U(e k h), k = 1, . . . , t, are mere evaluations. And the cost function J is discretized by the formula h

N -1 k=0 t i=1 bi C( Qk i , P k i , Û k i )
, where with a similar notation the "hat" stands for weights and evaluations related to the e's. Now, different cases arise:

• If r > t, one does an underevaluation of the controls within the discrete cost functional which will allow these to diverge (for instance a control could not appear explicitly in the discrete cost and therefore take any arbitrary value). • If on the contrary r < t, one does an overevaluation of the controls which, in the case of a coercive discrete cost function, will only increase the computational cost with, a priori, no convergence benefits. • Therefore, the case r = t seems to be the best fit, which is the case where there is a minimal number of evaluations of the controls (each control is evaluated just once in the discrete cost) so to have convergence under the assumption of coercivity. Assuming the last situation and continuing with the discussion, further cases arise:

• On the one hand, if r > s, the convergence rate of the controls will be limited by the accuracy of the discrete dynamics. • On the other hand, if r < s, the convergence rate of the controls will suffer from a computational payload by the high accuracy of the dynamics. • Therefore, the case r = s seems again to be the best fit, which is the case that equates accuracy of the dynamics with convergence rate of the controls. Finally, under the assumption r = s = t, choosing a unique set of collocation points 0 ≤ c 1 < . . . < c s ≤ 1 (we drop the "hats", "bars", e's and d's), one minimizes the number of polynomial evaluations and therefore the total computational cost (besides of simplifying the problem).

In the following example, we solve a simple optimal control problem with a linear dynamical constraint and a quadratic cost function. The numerical experiments show, in the spirit of [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF] and the previous discussion before it, how a good choice and proper combination of the discretization gives a "fast" convergence of the scheme, while other combinations show "slow" convergence or even divergence of the controls, all of it exemplifying Theorem 4.1.

Example 4.6. Consider the problem min q, q,u T 0 ( q2 + u 2 ) dt (35a) s.t. q = 1 + u , (q(0), q(0)) = (0, 0) (35b) for which the functions

q(t) = cosh(t) -1 cosh(T ) and u(t) = cosh(t) cosh(T ) - 1
are the unique solution. We identify from the forced Euler-Lagrange equation (35b) the Lagrangian function L(q, q) = 1 2 q2 + q and the control force F (q, q, u) = u. The density cost function is obviously C(q, q, u) = q2 + u 2 .

We discretize the mechanical system by using a symplectic Galerkin approach together with a Lobatto quadrature for s = 3 points. We initially assume that the controls are also discretized by r = 3 nodes. Then the right-hand side equations of (29b) or (34c) are

Ṗ1 = -p0 h/6 + 1 h (4 Q1 + 2 Q2 ) = 1 + U 1 , (36a) Ṗ2 = 1 h (-Q1 + Q3 ) = 1 + U 2 , (36b) Ṗ3 = p1 h/6 + 1 h (2 Q2 -4 Q3 ) = 1 + U 3 , (36c)
where the micro-veloticies Qi are given by the left equations of (29b) or (34c), which are in this particular case

  Q1 Q2 Q3   = 1 h   -3 4 -1 -1 0 1 1 -4 3     Q 1 Q 2 Q 3  
For the cost function, we consider four different discretizations with Lobatto's quadrature rule for t = 1, 2, 3, 4 quadrature points. These are respectively

C d (q h , p h , u h ) = h Q2 2 + U 2 2 , (37a) C d (q h , p h , u h ) = h 2 Q2 1 + Q2 3 + U 2 1 + U 2 3 , (37b) C d (q h , p h , u h ) = h 6 Q2 1 + 4 Q2 2 + Q2 3 + U 2 1 + 4U 2 2 + U 2 3 , (37c) C d (q h , p h , u h ) = h 12 (3 - √ 5) Q2 1 + (3 - √ 5) Q2 2 + Q2 3 + (37d) +(1 -1/ √ 5) ( Q1 + Q2 ) 2 + ( Q2 + Q3 ) 2 + h 30 2(U 1 + U 2 ) 2 + (U 1 -U 3 ) 2 + 2(U 2 + U 3 ) 2 +U 2 1 + 12U 2 2 + U 2 3 .
The first two discretizations, Equations (37a) and (37b), are clearly not coercive with respect to the controls (U 1 and U 3 are missing in (37a) and U 2 is missing in (37b)), which will be allowed to diverge (see Figures 2a and2b). Nonetheless the last two discretizations, Equations (37c) and (37d), are indeed coercive, still one outperforms the other in terms of convergence (see Figures 2c and2d). The discrete cost function (37d), besides of having a higher computational cost, shows a slower convergence rate. The discrete cost function (37c) corresponds to the method presented in Problem [START_REF] Mackay | Some aspects of the dynamics and numerics of Hamiltonian systems[END_REF] and Theorem 4.1.

We continue by assuming that the approximated control U is determined only by two points, that is

U(t) = U 1 + t( U 3 -U 1 )
(we note U 3 instead of U 2 to make the notation more appealing). The previous set of Equations [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF] and [START_REF] Maurer | Second order sufficient conditions for time-optimal bang-bang control[END_REF] are then updated by merely substituting the controls by

U 1 = U 1 , U 2 = 1 2 ( U 1 + U 3 )
, and U 3 = U 3 , which leads to

C d (q h , p h , u h ) = h Q2 2 + ( U 1 + U 3 ) 2 /4 , (38a) C d (q h , p h , u h ) = h 2 Q2 1 + Q2 3 + U 2 1 + U 2 3 , (38b) C d (q h , p h , u h ) = h 6 Q2 1 + 4 Q2 2 + Q2 3 + U 2 1 + ( U 1 + U 3 ) 2 + U 2 3 , (38c) C d (q h , p h , u h ) = h 12 (3 - √ 5) Q2 1 + (3 - √ 5) Q2 2 + Q2 3 + (38d) +(1 -1/ √ 5) ( Q1 + Q2 ) 2 + ( Q2 + Q3 ) 2 + h 6 U 2 1 + ( U 1 + U 3 ) 2 + U 2 3 .
In this occasion, only the first discretization, Equation (38a), defines a non-coercive discrete cost function (see Figure 3a). From the rest (see Figures 3b-3d), Equations (38c) and (38d) show the fastest convergence rate, although slow on the controls and with a computational payload for (38d). Equation (38c) corresponds to a discretization of the cost with three quadrature points.

Commutation of discretization and dualization

In this section we investigate the equivalence (C) in Figure 1 for the special choice of sG discretization of the optimal control problem 2.2. To this end, we analyze and compare the adjoint systems for the continuous and the discrete optimal control problems.

Throughout the section we assume that all controls under consideration do not saturate the constraints. In other words, we assume that the optimal control is in the interior of the set of constraints on controls. This is obviously the case if U = R m . This assumption allows us to avoid the typical situation of bang-bang controls, and under slight extra conditions, to derive (from the Pontryagin Maximum Principle) extremal controls that are smooth functions of the state and costate. The necessary optimality conditions [START_REF] Colombo | Discrete second-order Euler-Poincaré equations. Applications to optimal control[END_REF] are

λ = -∇ q C -λ • ∇ q f -ψ • ∇ q g , λ(T ) = ∇ q Φ(q(T ), p(T )) , (39a) ψ = -∇ p C -λ • ∇ p f -ψ • ∇ p g , ψ(T ) = ∇ q Φ(q(T ), p(T )) , (39b) 0 = ∇ u C + ψ • ∇ u g . (39c)
If we use the sG integrator for the discretization of Problem 2.2, we obtain the discretized optimal control problem 3.4. To derive the necessary optimality conditions for the discretized optimal control problem, we introduce the discrete adjoint vectors (covectors in R n ) λ 0 , . . . , λ N , µ 0 , . . . , µ N -1 ,

ψ 0 , Λ 0 i , . . . , Λ N -1 i , Ψ 0 i , . . . , Ψ N -1 i
, i = 1, . . . , s, and define the discrete optimal control Lagrangian as

L d = N -1 k=0 s i=1 hb i C k i + Φ(q N , p N ) -λ 0 • (q 0 -q 0 ) -ψ 0 • (p 0 -p 0 ) + N -1 k=0   µ k •   q k - s j=1 α j Q k j   -λ k+1 •   q k+1 - s j=1 β j Q k j   + s i=1 Λ k i •   hf k i - s j=1 a ij Q k j   + Ψ k i •   hg k i - β i p k+1 -α i p k bi - s j=1 āij P k j     , (40)
where C k i is a short notation for C(Q k i , P k i , U k i ) (analogously for f k i and g k i ). The necessary optimality conditions (KKT equations) are derived by differentiation w.r.t. the discrete variables q k , p k , k = 0, . . . , N and Q k i , P k i , U k i , k = 0, . . . , N -1, i = 1, . . . , s, which leads to for k = 0, . . . , N -1 :

µ k -λ k = 0, (41a) ∇ q Φ(q N , p N ) -λ N = 0, (41b) -ψ 0 + s i=1 α i bi Ψ 0 i = 0, (41c) for k = 1, . . . , N -1 : s i=1 α i bi Ψ k i - s i=1 β i bi Ψ k-1 i = 0, (41d) ∇ p Φ(q N , p N ) - s i=1 β i bi Ψ N -1 i = 0, (41e) 1 h -α i µ k + β i λ k+1 - s j=1 a ji Λ k j + b i ∇ q C k i + Λ k i • ∇ q f k i + Ψ k i • ∇ q g k i = 0, (41f) - 1 h s j=1 āji Ψ k j + b i ∇ p C k i + Λ k i • ∇ p f k i + Ψ k i • ∇ p g k i = 0, (41g) b i ∇ u C k i + Ψ k i • ∇ u g k i = 0, (41h)
with k = 0, . . . , N -1, i = 1, . . . , s, for the last three equations. We transform the necessary optimality conditions by defining

Γ k i := Λ k i / bi and χ k i := Ψ k i / bi for k = 0, . . . , N -1, i = 1, . . . , s, (42a) ψ - k := s i=1 α i χ k i , k = 0, . . . , N -1, and ψ + k := s i=1 β i χ k-1 i , k = 1, . . . , N, (42b)
such that Equation (41d) reduces to ψ - k = ψ + k := ψ k . By eliminating the variables µ 0 , . . . , µ N -1 with Equation (41a) and by exploiting the conditions on the coefficients b i a ij + bj āji = 0 and b i = bi , we obtain the following discrete adjoint system

ψ k = s j=1 α j χ k j , (43a)
ψ k+1 = s j=1 β j χ k j , (43b) -∇ q C k i -Γ k i • ∇ q f k i -χ k i • ∇ q g k i = β i λ k+1 -α i λ k h bi + 1 h s j=1 āij Γ k j , (43c) -∇ p C k i -Γ k i • ∇ p f -χ k i • ∇ p g k i = 1 h s j=1 a ij χ k j , (43d)
∇ u C k i + χ k i • ∇ u g k i =0 , (43e)
for k = 0, . . . , N -1, i = 1, . . . , s, and with final conditions

(43f) λ N = ∇ q Φ(q N , p N ) and ψ N = ∇ p Φ(q N , p N ) , where (43g)
b i a ij + bj āji = 0 and b i = bj .

Note that the adjoint scheme (43a)-(43d) together with the final constraints (43f) and the conditions on the coefficients (43g) is exactly the symplectic Galerkin integrator applied to the adjoint system (39a)-(39b). To ensure that the discrete adjoint system (43) is indeed equivalent to the necessary optimality conditions defined in [START_REF] Ober-Blöbaum | Variational integrators for electric circuits[END_REF] we show the following proposition.

Proposition 5.1. If bi > 0 for each i, then the necessary optimality conditions (41) and the dis-

crete adjoint system (43) are equivalent. That is, if (µ 0 , . . . , µ N -1 , Λ 0 i , . . . , Λ N -1 i , Ψ 0 i , . . . , Ψ N -1 i ), i = 1, . . . , s, satisfy (41), then (43) hold for (ψ k , Γ k i , χ k i ) defined in (42a) and (42b). Conversely, if (ψ 1 , . . . , ψ N , Γ 0 i , . . . , Γ N -1 i , χ 0 i , . . . , χ N -1 i ), i = 1, . . . , s, satisfy (43)
, then (41) hold for (µ k , Λ k i , Ψ k i ) defined in (41a) and (42a).

Proof. We already derived the adjoint system (43) starting from the necessary optimality conditions [START_REF] Ober-Blöbaum | Variational integrators for electric circuits[END_REF]. We now suppose that (ψ 1 , . . . , ψ N , Γ 0

i , . . . , Γ N -1 i , χ 0 i , . . . , χ N -1 i
), i = 1, . . . , s, satisfy the adjoint system [START_REF] Ross | Legendre pseudospectral approximations of optimal control problems[END_REF]. Equation (41a) holds by assumption. The condition for λ N in (43f) and (41b) are identical. The condition for ψ N in (43f) together with Equation (43b) for k = N -1 and the definition (42a) yields (41e) whereas Equation (43a) for k = 0 together with definition (42a) yields (41c). By subtracting Equations (43a) and (43b) for the same index k and using definition (42a) we obtain (41d). Finally, by taking the condition (43g) on the coefficients into account, (43c)-(43e) and definition (42a) yield (41f)-(41h), respectively.

With the classical Legendre assumption, i.e. (∂ 2 H/∂u 2 )(q * , p * , u * , λ, ψ, 1) is a positive definite symmetric matrix, with Equation (39c) u can be expressed as function of the states and the adjoints, u = u(q, p, λ, ψ). We denote by ν and η the functions defined by ν(q, p, λ, ψ) = (-∇ q C(q, p, u)λ • ∇ q f (q, p)ψ • ∇ q g(q, p, u)) | u=u(q,p,λ,ψ) , η(q, p, λ, ψ) = (-∇ p C(q, p, u)λ • ∇ p f (q, p)ψ • ∇ p g(q, p, u)) | u=u(q,p,λ,ψ) .

With some abuse of notation, let g(q, p, λ, ψ) denote the function g(q, p, u(q, p, λ, ψ)). In the case where the control has the form

U k i = u(Q k i , P k i , Γ k i , χ k i )
, the state and adjoint scheme based on the symplectic Galerkin integrator can be expressed as

q k = s j=1 α j Q k j , q k+1 = s j=1 β j Q k j , ψ k = s j=1 α j χ k j , ψ k+1 = s j=1 β j χ k j , (44a)
f k i = 1 h s j=1 a ij Q k j , g k i = β i p k+1 -α i p k h bi + 1 h s j=1 āij P k j , (44b)
η k i = 1 h s j=1 a ij χ k j , ν k i = β i λ k+1 -α i λ k h bi + 1 h s j=1 āij Γ k j , (44c) k = 0, . . . , N -1, i = 1, . . . , s, (44d) q 0 = q 0 , p 0 = p 0 , λ N = ∇ q Φ(q N , p N ), ψ N = ∇ p Φ(q N , p N ),
where b i a ij + bj āji = 0 and b i = bj and where f k i and g k i are short notations for f (Q k i , P k i ) and g(Q k i , P k i , Γ k i , χ k i ) (analogously for η k i and ν k i ). Scheme (44) can be viewed as symplectic Galerkin discretization of the two-point boundary value problem q = f (q, p) , q(0) = q 0 , (45a) ṗ = g(q, p, λ, ψ) , p(0) = p 0 , (45b) λ = ν(q, p, λ, ψ) , λ(T ) = ∇ q Φ(q(T ), p(T )) , (45c) ψ = η(q, p, λ, ψ) , ψ(T ) = ∇ q Φ(q(T ), p(T )) ,

where the variables (q, ψ) and (p, λ) are treated in the same way, respectively. Since the same discrete scheme is used for state and adjoint system, the orders of approximation coincide. This leads to the following statement. Theorem 5.2 (Commutation property). Given the (OCP ) 2.2, besides of H1, H3, H4 from Theorem 4.1, we assume that (∂ 2 H/∂u 2 )(q * , p * , u * , λ, ψ, 1) is a positive definite symmetric matrix. If a convergent symplectic Galerkin method with b i > 0, i = 1, . . . , s, is used for the discretization of the state system, dualization and discretization commute, i.e. the dualization of Problem 3.4 coincides with the sG discretization of the boundary value problem [START_REF] Sanz-Serna | Numerical Hamiltonian problems[END_REF].

Remark 5.3. In spirit of the Covector Mapping Principle (see [START_REF] Gong | Connections between the covector mapping theorem and convergence of pseudospectral methods for optimal control[END_REF]), the order-preserving map between the adjoint variables corresponding to the dualized discrete problem (KKT) and the discretized dual problem (discrete PMP) is given by Equation (42a).

Remark 5.4. If the controls do not saturate the constraints, i.e. control constraints are active, the optimal solution is typically only Lipschitz continuous. Then we expect analogously to [START_REF] Dontchev | Second-order Runge-Kutta approximations in control constrained optimal control[END_REF] that convergence rates are limited to order two even for higher order approximation schemes.

Conclusions

In this work, we investigate the application of high order variational integrators to the numerical solution of optimal control problems of mechanical systems. We derive two different schemes of high order variational integrators, the spRK and the sG method, which are both used for the discretization of a Lagrangian optimal control problem. The convergence of the primal variables of the resulting discrete optimal control problems is proven. Furthermore, the commutation of dualization and discretization for the sG method is shown, which extends the result in [START_REF] Ober-Blöbaum | Discrete mechanics and optimal control: an analysis[END_REF] to another class of variational integrators that fulfills this commutation property and directly implies that the Covector Mapping Principle is satisfied. In particular, due to the commutation, not only the order of the adjoint scheme but also the discretization method itself is preserved and in contrast to Legendre pseudospectral methods or classical Runge-Kutta methods, no additional closure conditions (see [START_REF] Gong | Connections between the covector mapping theorem and convergence of pseudospectral methods for optimal control[END_REF]) or conditions on the Runge Kutta coefficients (see [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF]), respectively, are required.

The fulfillment of the Covector Mapping Principle provides a convenient way to prove the convergence of the dual variables (as done, for example, in [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF]). With Theorem 5.2 the solution of the discretized direct problem coincides with the discrete solution of a shooting method applied to the necessary conditions of optimality. By showing the convergence of the shooting approach, we can conclude directly the convergence of the direct approach. The convergence proof for the adjoint variables is left for future work.

In the present paper we restricted ourselves to optimal control problems without any constraint on the final state. If we consider more general optimal control problems, involving constraints on the final state, then we expect that we will have to use the general conjugate point theory (see [START_REF] Bonnard | Second order optimality conditions in the smooth case and applications in optimal control[END_REF]), in order to provide second-order conditions for optimality, related with the classical sensitivity analysis along a given optimal trajectory. Also, if there are some constraints on the final point then abnormal extremals may occur in the application of the Pontryagin Maximum Principle, which may raise a major problem in the analysis. Fortunately, it is known that abnormal minimizers do not exist under generic assumptions on the system and on the cost (see [START_REF] Chitour | Singular trajectories of control-affine systems[END_REF]), and we expect that the results presented in this paper may hold in such a generic context. Otherwise the possible presence of abnormal minimizers is responsible for a loss of compactness (in particular, adjoint vectors do not stay in a compact anymore, see [START_REF] Trélat | Some properties of the value function and its level sets for affine control systems with quadratic cost[END_REF]), which may imply the failure of our method of proof, and it is not very clear then if we can expect that the Covector Mapping Principle hold true in that case. These general considerations will be investigated in future work. We stress again that, in the present paper, we have restricted ourselves to a more simple and tractable case. Once again, note that we have considered control-affine systems with (quasi)-quadratic costs, with assumptions ensuring the smoothness of optimal controls (as functions of the state and of the costate). As in [START_REF] Hager | Runge-Kutta methods in optimal control and the transformed adjoint system[END_REF], our theory can be applied whenever there exist control constraints, provided the optimal controls under consideration belong to the interior of the set of constraints. Otherwise, in the general case controls may saturate the constraints, typically bang-bang controls do appear and then additional assumptions must be done in order to ensure a nice regularity of controls as functions of the state and of the costate. For instance, it is desirable to avoid chattering phenomena, in which a given optimal bang-bang control can have an infinite number of commutations over a compact interval. Also, we expect that one has to use the corresponding conjugate point theory in the bang-bang case. Such a theory does exist in the purely bang-bang case (see [START_REF] Agrachev | A hamitonian approach to strong minima in optimal control[END_REF][START_REF] Maurer | Second order sufficient conditions for time-optimal bang-bang control[END_REF]) but then requires additional assumptions (ruling out, in particular, chattering). It can be noted that a general conjugate point theory, involving all possible subarcs -free, bang, singular, boundary -is still to be done. In order to get a general Covector Mapping Principle, such a complete theory is certainly required.

Another interesting issue is the application of more general symplectic and structure preserving methods in optimal control. Whereas the commutation property is shown for an already rich class of symplectic integrators (spRK and sG) it is still an open question, if it is satisfied for any variational and thus, any symplectic integrator (note that the classes of variational and symplectic integrators are identical, see e.g. [START_REF] Hairer | Geometric numerical integration[END_REF]).

( 19 )

 19 δL d ( Q1 , . . . , Qs ) • (δ Q1 , . . . , δ Qs ) + F d ( Q1 , . . . , Qs ) • (δQ 1 , . . . , δQ s ) = 0

Figure 1 .

 1 Figure 1. Optimal control for mechanical systems: the order of variation, dualization and discretization for deriving the necessary optimality conditions.

Theorem 4 . 1 .

 41 Given a Lagrangian function L : T Q → R, an external control force F : T Q × U → T * Q, a density cost function C : T * Q×U → R and a set of collocation points 0 ≤ c 1 < . . . < c s ≤ 1 defining a quadrature rule (c i , b i ), let us assume that (H1) L is regular;

  q, slope=3.4634 error in u, slope=0.97128(d) cost function (37d)

Figure 2 .

 2 Figure 2. Convergence behavior of the discrete solution for the discrete cost functions given in (37).

Figure 3 .

 3 Figure 3. Convergence behavior of the discrete solution for the discrete cost functions given in[START_REF] Ober-Blöbaum | Galerkin variational integrators and modified symplectic Runge-Kutta methods[END_REF].

Table 1 .

 1 Comparison of s-stage variational integrators.

	1

Acknowledgments

This work has been partially supported by the European Union under the 7th Framework Programme FP7-PEOPLE-2010-ITN, grant agreement number 264735-SADCO, the Spanish MINECO under the National Reasearch Programme I+D+i 2010-2013, grant agreement number MTM2010-21186-C02-01. One of the authors, C.M.C., thanks the Universidad de Valladolid for its current postdoc position cofunded by the European Social Fund and the Junta de Castilla y León.