Towards Global Optimal Control via Koopman Lifts

This paper introduces a framework for solving time-autonomous nonlinear infinite horizon optimal control problems, under the assumption that all minimizers satisfy Pontryagin's necessary optimality conditions. In detail, we use methods from the field of symplectic geometry to analyze the eigenvalues of a Koopman operator that lifts Pontryagin's differential equation into a suitably defined infinite dimensional symplectic space. This has the advantage that methods from the field of spectral analysis can be used to characterize globally optimal control laws. A numerical method for constructing optimal feedback laws for nonlinear systems is then obtained by computing the eigenvalues and eigenvectors of a matrix that is obtained by projecting the Pontryagin-Koopman operator onto a finite dimensional space. We illustrate the effectiveness of this approach by computing accurate approximations of the optimal nonlinear feedback law for a Van der Pol control system, which cannot be stabilized by a linear control law.


Introduction
During the last decades algorithms and software for computing local solutions of nonlinear optimal control problems have reached a high level of maturity [4,9,35], enabling its deployment in industrial applications [29]. By now, such local solutions can be computed within the milli-and microsecond range. Moreover, auto-generated implementations of real-time local optimal control solvers run on embedded hardware systems [24,18]. In contrast to these developments in local optimal control, algorithms for locating global minimizers of non-convex optimal control problems, can hardly be claimed to be ready for widespread industrial application. There are at least two reasons for this. Firstly, generic global optimization methods can often only be applied to problems of modest dimensions. And secondly, their run-time usually exceeds the run-time of local solvers by orders of magnitude. Nevertheless, in the last years there have been a number of promising developments in the field of global optimal control, which are reviewed next.
Dynamic Programming (DP) [3], which proceeds by approximately solving the Hamilton-Jacobi-Bellman (HJB) equation [13], is a historically important method able to find globally optimal feedback laws. In [23] and [15] tailored discretization grids for DP implementations have been developed that can successfully solve nonlinear optimal control problems, as long as the dimension of the state of the system is small. For higher dimensional state-spaces these methods are, however, not applicable as the complexity of storing and processing the value functions during the DP recursion grows exponentially with the number of states. Other deterministic methods for global optimization problems involving nonlinear differential equations are based on Branch-and-Bound (BB) [7,11,22,32] and its variants, including α-BB methods [27,8]. These BB algorithms have in common that they can effectively solve problems with a small number of decision variables. Unfortunately, implementations of branch-and-bound search easily run out of memory in higher dimensional spaces due to exponentially growing search trees. In particular, as the discretization of control inputs typically leads to a large number of optimization variables, BB-methods are usually not suited for solving optimal control problems.
An alternative to Branch-and-Bound are the so-called Branch-and-Lift (BL) methods [16]. In contrast to BB, these methods never discretize the control inputs directly. Instead, these methods branch over orthogonal projections of the control function in lower dimensional subspaces of increasing dimension, until an -suboptimal global control input is found. A rigorous analysis of the mathematical properties of such BL methods can be found in [17]. However, in the context of developing practical implementations of Branch-and-Lift one faces numerical challenges. In particular, these methods require the computation of accurate enclosures of moment-constrained reachable sets of nonlinear differential equations. At the current state of research, the lack of generic methods for computing such moment-constrained reachable set enclosures limits the applicability of Branch-and-Lift-currently, some problems have been solved successfully with BL by using enclosure propagation methods that are tailored to particular applications [12,16].
In order to mitigate the above-mentioned limitations of existing global optimal control methods, this paper proposes a completely new framework for constructing methods for solving nonlinear infinite horizon optimal control problems. Here, the main idea is to analyze the spectrum of a Koopman operator that is associated with Pontryagin's differential equation-under the assumption that Pontryagin's optimality condition [28] is satisfied at the optimal solution. In order to understand the contributions of this paper, it is important to first be aware of the fact that Koopman operators have originally been introduced almost 90 years ago-in fact, by Koopman himself [19]. This theory has later been extended for more general nonlinear differential equations in [26], where the concept of Koopman mode analysis has been introduced. A more complete analysis of these Koopman mode-based methods from the field of ergodic theory has, however, only appeared much later [1]. Notice that a mathematical analysis of finite dimensional approximations of the Koopman operator can be found in [14] and a variety of related applications of approximate Koopman mode analysis methods can be found in [6]. Moreover, in [20] Koopman mode estimation heuristics for data driven control have been introduced.

Contribution.
The key idea of this paper is to introduce a Koopman operator that is associated with the Pontryagin differential equation of a nonlinear infinite-horizon optimal control problem-in the following we called a Pontryagin-Koopman operator -which can be used to construct globally optimal feedback laws. Section 2 introduces infinite horizon optimal control problems and briefly reviews Pontryagin's necessary condition for optimality. This review section is followed by a presentation the main contributions of this article, which can be outlined as follows.
(1) Section 3 uses ideas from the field of symplectic geometry [2] to characterize the structural properties of Pontryagin-Koopman operators. Section 4 leverages on these structural properties in order to perform a symplectic spectral analysis that finally leads to a complete characterization of globally optimal feedback control laws (see Section 5, Theorem 3). Because this is the first time that the symplecticity of Pontryagin-Koopman operators is analyzed and used to derive practical characterizations of optimal control laws, Theorem 3 can be regarded as the main theoretical contribution of this paper. (2) In order to avoid misunderstandings, it is mentioned clearly that this paper does not claim to solve all the numerical issues regarding the discretization of Pontryagin-Koopman operators that would have to be solved to develop a generic software for global optimal control. However, Section 6 illustrates the practical applicability of the proposed theoretical framework by designing an optimal regulator for a controlled Van der Pol oscillator. Although this implementation is only based on a naive Galerkin projection of the Pontryagin-Koopman operator onto a finite dimensional Legendre basis, this section successfully applies the proposed framework to construct accurate approximations of nonlinear globally optimal feedback laws-in this example, for a system that cannot be stabilized by a linear control law.
The above contributions are relevant for the future of optimal control algorithm development, as they pave the way towards the development of practical procedures for approximating globally optimal control laws of a very general class of nonlinear systems. Therefore, Section 7 does not only conclude the paper, but also outlines the potential and relevance of the proposed framework for future research in control.
Notation. The distance of a point x ∈ R nx to a trajectory φ : R → R n is denoted by We denote with L n p the set of (potentially complex-valued) Lebesgue measurable functions ϕ : R n → C whose p-th power of the absolute value is integrable on R n . Moreover, we use the notation W n k,p to denote the Sobolev space of integrable functions on R n , whose weak derivatives up to order k are all in L n p . The symbol M denotes the Hermitian transpose of the matrix M ∈ C n×m . The symbol ∇ denotes the gradient operator. The associated second order derivative operator is denoted by ∇ 2 = ∇∇ . Last but not least, the support of a function ϕ : R n → C is denoted by supp(ϕ) = {x ∈ R n | ϕ(x) = 0} .

Infinite horizon optimal control
This section introduces infinite horizon optimal control problems and briefly summarizes existing methods for analyzing the local stability properties of optimal periodic limit orbits (see Section 2.2). Moreover, Sections 2.3 and 2.4 review Pontryagin's necessary optimal condition thereby introducing the notation that is used throughout the paper.

Problem formulation
This paper considers infinite horizon optimal control problems of the form Here, x : R → R nx denotes the state trajectory and u : R → R nu denotes the control input. The initial value x 0 ∈ R nx is assumed to be given. Throughout this paper the following assumptions are imposed.
Assumption 1 The function f : R nx × R nu → R nx is twice continuously differentiable and globally Lipschitz continuous in x.
The constant l ∈ R in (1) denotes the optimal average cost, x(0) = x 0 , assuming that this limit exists. Notice that l does not need to be known explicitly in order to solve (1), as adding constant offsets to the objective does not affect the solution of an optimization problem. In this paper, the constant l is merely introduced for mathematical reasons, such that V (x 0 ) remains finite and well-defined on infinite horizons under mild regularity assumptions that shall be introduced further below.

Limit behavior and periodic orbits
In practical instances, one is often interested in whether optimal solutions for the state trajectory of (1) converge to an optimal steady-state or, in more generality, to an optimal periodic limit orbit. These optimal steady states or more general periodic orbits are defined as follows.

Definition 1
The function (x p , u p ) ∈ W nx 1,1 × L nx 1 is called an optimal periodic orbit, if there exists a period T > 0 such that for all t ∈ [0, ∞) For the special case that x p and u p are constant functions, this pair is called an optimal steady-state.
In order to prepare the following analysis, it is helpful to introduce the shorthands assuming that (x p , u p ) is an optimal periodic orbit. Here, f x and f u denote the partial derivatives of f with respect to x and u and an analogous notation is then also used for the mixed second order derivatives of l. The above matrix-valued functions can be used to construct sufficient conditions under which an optimal periodic orbit is locally stabilizable. Here, one relies on the theory of periodic Riccati differential equations [5] of the form It is well-known [5] that if Assumptions 1 and 2 are satisfied and if R has full rank, then the existence of a periodic and positive definite function P satisfying (3) is sufficient to ensure that the solution x of (1) converges to the optimal periodic orbit x p as long as dist (x 0 , x p ) is sufficiently small-that is, if x 0 is in a small neighborhood of the optimal orbit x p . However, this statement is of a rather local nature. This means that, if we wish to understand the global behavior of system (1), an analysis of the periodic Riccati equation is not sufficient. The following sections focus on analyzing global solutions of (1), under the assumption that dist (x 0 , x p ) is not necessarily small.

Pontryagin's differential equations
Pontryagin's maximum principle [28] can be used to derive necessary conditions for the minimizers of (1). The first order variational optimality condition is summarized as follows. Let H : R nx × R nu × R nx → R be the Hamiltonian function of (1), If Assumptions 1 and 2 hold, H is, by construction, twice continuously differentiable in all variables. Next, let denote the associated parametric minimizer of H. At this point, we introduce the following regularity assumption.

Assumption 3
The parametric minimizer u in (5) satisfies the second order sufficient condition, H uu 0.
Notice that for the practically relevant special case that f is affine in u and l strongly convex in u, Assumption (3) always holds whenever Assumptions 1 and 2 are satisfied.
Next, if Assumptions 1, 2, and 3 are satisfied, it is well-known [28] that any minimizer (x, u) of (1) necessarily satisfies Pontryagin's differential equatioṅ for a co-state function λ ∈ W nx 1 . In the following, we summarize this differential equation in the forṁ where y = [x , λ ] denotes the stacked state and F a stacked version of the right-hand side functions of (6) and (7).

Necessary boundary and limit conditions
Besides Pontryagin's differential equation (8), optimal solutions of (1) satisfy necessary boundary conditions. First, of course, the initial condition x(0) = x 0 must hold. Moreover, the co-state λ satisfies a necessary limit condition, which can be summarized as follows.
Proposition 1 Let Assumptions 1, 2, and 3 hold. Moreover, let (x , u ) be a primal optimal solution of (1) converging to an optimal periodic orbit (x p , u p ) ∈ W nx 1,1 × L nx 1 at which the periodic Riccati differential equation (3) admits a positive definite solution. Then, there exists a periodic function λ p ∈ W nx 1,1 such that the associated co-state λ necessarily satisfies lim Notice that Proposition 1 is-at least in very similar versions-well-known in the literature [28,21]. However, as this result is important for understanding the developments in this paper, the following proof briefly recalls the main argument, why this proposition holds.
Proof. Since (3) admits a positive definite solution, the value function V in (1) is well-defined and differentiable in a neighborhood of the optimal orbit x . Thus, (1) can be replaced by an equivalent finite-horizon optimal control problem as long as V is used as a terminal cost. Now, Pontryagin's principle for optimal control problems with Mayer terms [28] yields the boundary condition λ (t) = ∇V (x (t)) for any horizon length t > 0. Clearly, since x converges to the optimal periodic orbit x p , λ must converge to the associated dual limit orbit λ p (t) = ∇V (x p (t)). 2 In the next sections we will see that Proposition 1 implies that the states of Pontryagin's differential equation for stabilizable systems must evolve along certain stable manifolds. As it shall demonstrated, these manifolds can be characterized using a Koopman mode analysis [26].

Symplectic Koopman operators
This section analyzes the symplectic structure of the flow associated with Pontryagin's differential equation. These symplectic structures are needed to understand the properties of the Pontryagin Koopman operators. In fact, Section 3.2 introduces a symplectic test space-that is, a space of observables-in which the Pontryagin-Koopman operator inherits certain symplecticity properties of its underlying incompressible flow field. Moreover, Section 3.3 leverages on ideas from the field of symplectic geometry [2] in order to work out the symplectic dual of the Pontryagin-Koopman operator. Notice that these developments are the basis for the developments in Section 4, which uses a symplectic spectral analysis in order to characterize globally optimal control laws.

Symplectic Flows
Let Γ t : R 2nx → R 2nx denote the flow of Pontryagin's differential equations such that d dt Γ t (z) = F (Γ t (z)) and Γ 0 (z) = z for all z ∈ R 2nx . If Assumptions 1, 2, and 3 hold, then Γ t is a well-defined, continuously differentiable function. In the mathematical literature [2] it is well-known that Γ t is a so-called symplectic flow. In order to reveal this symplectic structure it is helpful to introduce the block matrix Now, the structural properties of Γ t can be summarized as follows.
Lemma 1 Let Assumptions 1, 2, and 3 be satisfied. Then, the function ∂ ∂z Γ t satisfies the equation that is, Γ t is symplectic function.
Proof. Let us first recall that the right-hand side F of Pontryagin's differential equation is given by (6) and (7), which can also be summarized as This implies that the derivative of F with respect y can be written in the form Here, we have used Assumptions 1 and 2 to ensure that the second derivatives of H with respect to x and u exist. Moreover, we have used Assumption 3, which implies that H uu is invertible and, as a consequence, that the implicit function theorem holds. In particular, we have In this form, it becomes clear that the matrix ΩF y is symmetric and we arrive at the intermediate result In order to proceed, we write the first order variational differential equation for Γ t in the form d dt Now, the main idea of this proof is to show that the function vanishes, ∆(t) = 0. Here, we first have ∆(0) = 0 by construction, since ∂ ∂z Γ 0 = I. Moreover, the derivative of ∆ with respect to time is given by∆ = 0 Thus, in summary, we must have ∆(t) = 0 for all t ∈ R, which yields the statement of this lemma. 2 The following corollary summarizes two immediate consequences of Lemma 1 which are both equivalent to stating that Γ t is an incompressible flow.
Corollary 1 Let Assumptions 1, 2, and 3 be satisfied. Then, Γ t is an incompressible flow; that is, we have Moreover, the divergence of the associated vector field F vanishes, div(F ) = ∇ F = 0.
Proof. By taking the determinant on both sides of (9) and using that det(Ω) = 1, we find that Since ∂ ∂z Γ 0 = I and since ∂ ∂z Γ t is continuous, this is only possible if (15) holds. Next, by taking the logarithm on both sides of (15) and differentiating with respect to time, we find 1 This completes the proof of the corollary. 2

Pontryagin-Koopman operators
The developments from the previous section can be used to analyze the structural properties of the Pontryagin-Koopman operator, which are defined to be the Koopman operator that is associated with the flow Γ t of Pontryagin's differential equation, or, more formally: Definition 2 (Pontryagin-Koopman Operator) We use the notation U t : W 2nx 1,2 → W 2nx 1,2 to denote the Pontryagin Koopman operator of Γ t , which is defined such that with • denoting the composition operator.
It is well known [19,26] that U t is a linear operator satisfying U t1+t2 = U t1 U t2 for all t 1 , t 2 ∈ R. Moreover, U t satisfies U −1 t = U −t , which follows by substituting t = t 1 = −t 2 in the previous equation and using that Γ 0 = id.
In order to introduce a notion of symplecticity in the space of observables of the Pontryagin-Koopman operator, we introduce the bilinear form for all ϕ, Φ ∈ W 2nx 1,2 . Because we have Ω = −Ω, the bilinear form ω is skew-symmetric, and W 2nx 1,2 , ω is a symplectic space [2].
Theorem 1 Let Assumptions 1, 2, and 3 be satisfied. Now, U t is for all t ∈ R a symplectic operator in the space W 2nx 1,2 , ω . This means that we have Proof. Assumptions 1, 2, and 3 imply that the function Γ t and its inverse Γ −1 t = Γ −t are both continuously differentiable. This implies in particular that the equivalence holds. Next, the definition of the operator U t and the chain rule for differentiation imply that for all Φ ∈ W 2nx 1,2 . Furthermore, because Assumptions 1, 2, and 3 hold, it follows from Lemma 1 that Let us multiply this equation with ∂ ∂z Γ t Ω from the left and substitute Ω 2 = −I. This yields Next, we multiply the latter equation with ∂ ∂z Γ t −1 Ω from the right and use the relation Ω 2 = −I once more, in order to arrive at the equation Thus, by substituting the previous equations, we get ∇ϕ(Γ t (z)) Ω∇Φ(Γ t (z)) dz .
In order to further transform this integral, we need to introduce the change of variables z = Γ t (z). Because Corollary 1 ensures that this change of variables is volume preserving. Thus, this implies that for all ϕ, Φ ∈ W 2nx 1,2 . Thus, U t is a symplectic operator. 2

Symplectic duality
In order to prepare the analysis of the properties of the spectrum of symplectic Koopman operators, it is helpful to introduce a notion of duality. However, instead of defining duality in a Hilbert space, we propose to introduce the following notion of symplectic duality. for all ϕ, Φ ∈ W 2nx 1,2 , then we call A the symplectic adjoint operator of A.
In the following, we consider a linear differential operator L : W 2nx 1,2 → W 2nx 1,2 that is given by for all Φ ∈ W 2nx 1,2 . Notice that this operator is found by differentiating the definition of the Pontryagin-Koopman operator U t with respect to time, such that L =U t .
Lemma 2 Let Assumptions 1, 2, and 3 hold. The operator L admits a symplectic adjoint L in the symplectic space (W 2nx 2,2 , ω). This adjoint is given by Proof. For any function Φ ∈ W 2nx 2,2 , we have which follows by differentiating (18) on both sides. Next, we know from Theorem 1 that U t is a symplectic operator. Thus, we can differentiate the equation ω(U t ϕ, U t Φ) = ω(ϕ, Φ) on both sides with respect to t, which yields 2,2 , we can expand the derivative on the right as and, after changing variables, z ← Γ t (z), recalling that det ∂ ∂z Γ t = 1, and resorting terms, Moreover, we recall that the equation also holds (see Equation (12)). Thus, in summary, we have Because the latter equation holds for all ϕ, Φ ∈ W 2nx 2,2 , we find that L = −L is the symplectic adjoint of L, as claimed by the statement of this lemma. 2 Remark 1 In many articles on Koopman operators, a so-called transport differential equation of the form is considered. If one introduces suitable boundary conditions, this advection PDE can be interpreted as the derivative of a Perron-Frobenius operator that is dual to the Koopman operator U t with respect to the standard L 2 -scalar product [10]. Notice that in our case, the symplectic adjoint operator happens to coincide with the right-hand operator of the above advection PDE. Thus, physically, one could interpret the operator L as an advection operator of the incompressible flow field F recalling that div(F ) = 0.

Spectral Analysis
In this section, we show that the symplectic structures of the Koopman operator, as analyzed in the previous section, have important consequences on its set of eigenvalues. Moreover, Section 5 uses these spectral properties of the Koopman operator to characterize global optimal control laws that are associated with the infinite horizon optimal control problem (1).

Eigenfunctions and Eigenvalues
As already mentioned in the introduction, the spectrum of general Koopman operators has been analyzed by many authors; for example in [1,6,26]. In the context of this paper, we call a weakly differentiable function Ψ an eigenfunction of the differential Pontryagin-Koopman operator L, if the Lebesgue measure of the set supp(Ψ) is either unbounded, or, if it exists, is not equal to 0, and for a potentially complex eigenvalue κ ∈ C. The fact that the symplectic adjoint of the operator L is given by L = −L has important consequences on its spectrum, which can be summarized as follows.

Remark 2
The statement of Theorem 2 is formally not directly applicable for eigenfunctions Ψ 1 , Ψ 2 of L that are locally weakly twice differentiable but whose derivatives are not square integrable, such that Ψ 1 , Ψ 2 are not elements of the Sobolev space W 2nx 2,2 . However, in practice, one is usually interested in constructing eigenfunctions of L on a compact domain C ⊂ R 2nx , in the following called the region of interest, such that This region of interest C can, for example, model a-priori bounds on the optimal primal and dual trajectories of (1).
Consequently, because we are simply not interested in how the flow Γ t and the associated eigenfunctions Ψ are defined outside of the set C, these functions can simply be redefined arbitrarily for x / ∈ C, for example, such that the desired eigenfunctions of L satisfy Ψ ∈ W 2nx 2,2 by construction. Thus, for the purpose of this paper, it is not restrictive at all to assume that the derivatives of the eigenfunctions of L are square integrable. Notice that this technique is also illustrated by our tutorial case study in Section 6, where we explain how to choose C and how to discretize (24) on C.
Let σ(L) ⊆ C denote the spectrum of L; that is, the set of eigenvalues of the linear operator L. In the following, we use the notation Ψ κ ∈ W 2nx 2,2 to denote an eigenfunction that is associated with an eigenvalue κ ∈ σ(L). Moreover, for a function q : R 2nx → R 2nx , we use the shorthand notation to denote the matrix that is obtained by evaluating the skew symmetric bilinear form ω for all possible combinations of the components of q. We call q a skew-orthogonal function in the symplectic space (W 2nx 1,2 , ω), if it satisfies ω(q, q) = I. Notice that the construction of such skew orthogonal functions is straightforward by using the standard skew-symmetric variant of the Gram-Schmidt algorithm [2].

Definition 4
The operator L is said to admit a spectral decomposition with respect to a given function q, if there exists a generalized function a : σ(L) → C 2nx such that Notice that the above definition uses the "control engineering notation" for generalized functions, which means that we use a as if it was a standard function, although this notation suppresses the distributional nature of a. Thus, in mathematical terms, this notation has to be translated as "a represents a distribution of order 0; that is, a linear operator on σ(L) that is Lipschitz continuous with respect to the L ∞ -norm". 2 The following statement is a consequence of Theorem 2.
Corollary 2 Let Assumptions 1, 2, and 3 hold. Let the operator L admit a spectral decomposition with respect to a skew-orthogonal function q. Then, there exist at least 2n x eigenvalues such that κ + i = −κ − i , where κ + i has a non-negative real part for all i ∈ {1, 2, . . . , n x }.
Proof. By substituting (25) in the equation ω(q, q) = I, we find that the equation holds. Let us have a closer look at the terms in (26). Clearly, the unit matrix on the left has full rank. But, on the other side, we have an integral over the rank 1 matrices a(κ )a(κ) . This integral term can only have full rank, if there are at least 2n x pairs of eigenvalues (κ , κ) ∈ σ(L) × σ(L) for which But now Theorem 2 implies that all these pairs must be such that κ = −κ and, after sorting all eigenvalues with respect to their real-part, we find that there must be at least n x eigenvalues with non-negative real part and n x associated mirrored eigenvalues with non-positive real-part, as claimed by the statement of this corollary. 2 Remark 3 Notice that Corollary 2 makes a statement about the spectrum of L under the assumption that this operator admits a spectral decomposition with respect to at least one skew orthogonal function. Thus, one should further ask the question under which assumptions the existence of such a spectral decomposition can be guaranteed for at least one such skew orthogonal function. In full generality, this question is difficult to answer, but sufficient conditions for the existence of (much more general and, in our context, sufficient) spectral decompositions can be found in [1,6,25,26], which use ideas from the field of ergodic theory [33] as well as Yoshida's theorem [34]. For example, if the monodromy matrix, ∂ ∂z Γ T (x p (0)), of a periodic optimal orbit with period length T > 0 is diagonalizable, one can ensure that a spectral decomposition is possible. This sufficient condition follows by applying Proposition 3 in [25] to the Pontryagin differential equation. A more complete review of such results from the field of functional analysis would, however, go beyond the scope of the present paper.

Optimal Feedback Control Laws
The theoretical results from the previous sections can be used to derive optimality conditions, which, in turn, can be used to develop practical numerical algorithms for solving (1). Let us introduce the set σ + (L) = {κ ∈ σ(L) | Re(κ) > 0} of unstable eigenvalues and its associated invariant manifold Because we assume that the eigenfunctions Ψ κ are weakly differentiable, we may assume without loss of generality that the functions Ψ κ are also continuous-otherwise there exists a continuous functionΨ κ (x) = Ψ κ (x) for almost every x ∈ R 2nx and we can useΨ κ instead of Ψ κ recalling that we work with Sobolev spaces, in which such arguments are indeed possible.
In the following, we say that µ is a regular optimal control law of (1), if the closed-loop trajectories, are minimizers of (1) at which Pontryagin's necessary optimality conditions are satisfied such that x converges to an optimal periodic orbit at which the conditions of Proposition 1 are satisfied.
Proof. Since L is a time-autonomous infinitesimal generator of U t , the eigenfunctions of L also satisfy [26] U t Ψ κ = e κt Ψ κ (27) for all κ ∈ σ(L). The remainder of the proof is divided into two parts. In the first part, we show that (27) implies that the optimal periodic limit orbit is a subset of the manifold M + and, in the second part, we use this property of the limit orbit to construct a globally optimal feedback law.
Part I: Let (x p , λ p ) denote the optimal periodic limit orbit of (1) with period time T . If we would have (x p (t), λ p (t)) / ∈ M + for a time t ∈ [0, T ], we must have Ψ κ (x p (t), λ p (t)) = 0 for at least one κ ∈ σ + (L). Because the optimal periodic orbit satisfies Pontryagin's differential equation, this implies in particular that as κ has a strictly positive real part. But this is a contradiction, since (x p , λ p ) is periodic. Thus, in summary, we must have Part II: Let (x, λ) denote an optimal solution of (1) with x(0) = x 0 . Now, if we would have (x(0), λ(0)) / ∈ M + , then there would have to exist at least one κ ∈ σ + (L) for which Ψ κ (x(0), λ(0)) = 0 .
But if this would be the case, then we would also have as Ψ κ is strictly unstable. But this limit statement is in conflict with (28), since Ψ κ is continuous and we assume that (x(t), λ(t)) converges to the optimal periodic limit orbit. Thus, in summary, there exists for every x 0 a λ 0 = λ(0) with (x 0 , λ 0 ) ∈ M + and the corresponding map from x 0 to λ 0 can be denoted by Λ. The associated optimal control input is given by This is already sufficient to establish the statement of this theorem, as the optimal feedback law must be timeautonomous. 2 Notice that Theorem 3 can be used to systematically search for globally optimal solutions of (1). Here, one of the key observations is that if, in addition to the assumptions of Theorem 3, the assumptions of Corollary 2 are also satisfied and if none of the non-negative eigenvalues κ + 1 , κ + 2 , . . . , κ + nx happens to be on the imaginary axis, then the parametric equation system consists of at least n x independent equations while the number of variables, dim(λ) = n x , is equal to n x , too. Thus, this equation system can, in many practical instances, be expected to admit a finite number of parametric solutions λ 0 = Λ(x 0 ) only. This means that, if all the above regularity assumptions are satisfied, Theorem 3 singles out a finite number of candidate control laws µ (x) = u (x, Λ(x))-at least one of which must be globally optimal.
This section illustrates how Theorem 3 can be used to construct accurate approximations of globally optimal control laws of (1) for a Van der Pol oscillator system with and l(x, u) = 1 2 The associated system,ẋ(t) = f (x(t), u(t)), has a linearly uncontrollable equilibrium at x = (0, 0) . Notice that for this particular example an explicit solution of the Hamilton-Jacobi-Bellman equation is known [30]: it turns out that the optimal value function V and the globally optimal feedback law µ are given, respectively, by In the following, this explicit solution is, however, only used to assess the accuracy of the proposed method; that is, the numerical procedure below neither knows the above expression for V nor for µ .

Galerkin discretization
In order to illustrate the practical applicability of the developments in Section 5, we introduce a simple Galerkin discretization of the operator L. Let ϕ 1 , . . . , ϕ N ∈ W 2nx 1,2 be orthogonal functions with respect to the standard L 2 -scalar product ·, · on W 2nx 1,2 . Next, if we compute the coefficients for all i, j ∈ {1, . . . , N }, the matrix M can be interpreted as a Galerkin approximation of the operator L over the subspace spanned by ϕ 1 , . . . , ϕ N . In our implementation, we set ϕ i to the ith multivariate Legendre polynomial on the 4-dimensional compact interval box C = [− 1 2 , 1 2 ] 4 and we set ϕ i (x, λ) = 0 outside of this domain; that is, for (x, λ) / ∈ C. Consequently, our discretization can only be expected to be accurate inside our region of interest C, but other choices for C and for the basis functions would be possible, too.
Remark 4 Standard Galerkin methods are, in general, numerically unstable when applied to advection operators and, consequently, although this method happens to yield reasonable approximations for our particular example, such naive discretization schemes cannot be recommended in general. More advanced discretization schemes for linear advection operators can be found in the modern PDE literature; see, for example [31]. A more complete discussion of such numerical discretization methods for Pontryagin-Koopman operators, is, however, beyond the scope of this paper.

Approximations of the optimal feedback law
The above Galerkin approximation of the operator L can be used to construct approximate eigenfunctions. In detail, if a ∈ R N is a left eigenvector of M with eigenvalue κ ∈ C, then is an approximation of an eigenfunction of L. The right plot in Figure 1 shows the spectrum of the matrix M for different choices of N (blue circles: N = 4, red squares: N = 15, and black circles with a cross: N = 35). In order to understand the structure of this spectrum it is helpful to recall Corollary 2, which predicts that there exits at least 2 eigenvalues κ 1 , κ 2 of L such that −κ 1 and −κ 2 are also eigenvalues of L. Notice that such symmetric eigenvalue pairs are indeed present in the spectrum of M , although M is only a Galerkin approximation of L. In order to further illustrate how the above spectral analysis of M can be used to construct approximations of the globally optimal control law µ , we can compute an approximation of the manifold M + by using the approximate eigenfunctions instead of the exact ones (see Theorem 3). For example, for N = 4, the matrix M has the eigenvalues κ ± 1,2 = 1 48 (β m ± β p √ −1) with β p = +983 + 96 √ 143 and β m = −983 + 96 √ 143 (34) and the associated Galerkin approximation of the globally optimal control law is given by The left plot in Figure 1 shows µ and compares it to the optimal feedback law µ . In fact, the squared integral error over C is approximately 6 × 10 −5 . Quite remarkably, this approximate feedback law can even be used to control the system for initial values outside of C. The plot in the middle of Figure 1 shows the corresponding trajectories of the closed-loop system that are obtained by using the approximately optimal feedback law µ.

Conclusions
This paper has presented an analysis of infinite horizon nonlinear optimal control problems, whose minimizers satisfy Pontryagin's necessary conditions of optimality. The proposed formalism is based on Pontryagin-Koopman operators, which have been shown to possess a symplectic structure, as revealed by Theorem 1. Moreover, Theorem 2 and Corollary 2 have established conditions under which the spectrum of the differential Pontryagin-Koopman operator contains at least 2n x mirrored eigenvalues. This spectral structure is used in Theorem 3 to characterize optimal control laws.