Milan Korda
email: milan.korda@epfl.ch

Didier Henrion
email: henrion@laas.fr

Colin N Jones
email: colin.jones@epfl.ch

Controller design and value function approximation for nonlinear dynamical systems

Keywords: Optimal control, nonlinear control, sum-of-squares, semidefinite programming, occupation measures, value function approximation

come

Controller design and value function approximation for nonlinear dynamical systems

Introduction

This paper considers the infinite-time discounted optimal control problem for continuoustime input-affine polynomial dynamical systems subject to polynomial state constraints and box input constraints. This problem has a long history in both control and economics literature. Various methods to tackle this problem have been developed, often based on the analysis of the associated Hamilton-Jacobi-Bellman equation.

In this work we take a different approach: We first lift the problem into an infinite-dimensional space of measures with continuous densities where this problem becomes convex; in fact a linear program (LP). This lifting is a tightening, i.e., its optimal value is greater than or equal to the optimal value of the original problem, and under suitable technical conditions the two optimal values coincide. This infinite-dimensional LP is then further tightened by restricting the class of functions to polynomials of a prescribed degree and replacing nonnegativity constraints by sufficient sum-of-squares (SOS) constraints. This leads to a hierarchy of semidefinite programming (SDP) tightenings of the original problem indexed by the degree of the polynomials. The solutions to the SDPs yield immediately a sequence of rational controllers, and we prove that, under suitable technical assumptions, the value functions associated to these controllers converge from above to the value function of the original problem.

We also describe how to obtain a sequence of polynomial approximations converging from above and from below to the value function associated to each rational controller. Combined with existing techniques to obtain polynomial under approximations of the value function of the original problem (adapted to our setting), this method can be viewed as a design tool providing a sequence of rational controllers asymptotically optimal in the original problem with explicit estimates of suboptimality in each step.

The idea of lifting a nonlinear problem to an infinite-dimensional space dates back at least to the work of L. C. Young [START_REF] Young | Lectures on the calculus of variations and optimal control theory[END_REF] and subsequent works of Warga [START_REF] Warga | Optimal control of differential and functional equations[END_REF], Vinter and Lewis [START_REF] Vinter | The equivalence of strong and weak formulations for certain problems in optimal control[END_REF], Rubio [START_REF] Rubio | Control and optimization: the linear treatment of nonlinear problems[END_REF] and many others, both in deterministic and stochastic settings. These works typically lift the original problem into the space of measures and this lifting is a relaxation (i.e., its optimal value is less than or equal to the optimal value of the original problem) and under suitable conditions the two values coincide.

More recently, this infinite-dimensional lifting was utilized numerically by relaxing the infinitedimensional LP into a finite-dimensional SDP [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF] or finite-dimensional LP [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF]. Whereas the LP relaxations are obtained by classical state-and control-space gridding, the SDP relaxations are obtained by optimizing over truncated moment sequences (i.e., involving only finitely many moments) of the measures and imposing conditions necessary for these truncated moment sequences to be feasible in the infinite-dimensional lifted problem. These finite-dimensional relaxations provide lower bounds on the value function of the optimal control problem and seem to be difficult to use for control design with strong convergence guarantees; a controller extraction from the relaxations is possible although no convergence (e.g., [START_REF] Henrion | Nonlinear optimal control synthesis via occupation measures[END_REF][START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF]) or only very weak convergence can be established (e.g., [START_REF] Korda | Controller design and region of attraction estimation for nonlinear dynamical systems[END_REF][START_REF] Majumdar | Convex optimization of nonlinear feedback controllers via occupation measures[END_REF] in the related context of region of attraction approximation).

Contrary to these works, in this approach we tighten the infinite-dimensional LP by optimizing over polynomial densities of the measures and imposing conditions sufficient for these densities to be feasible in the infinite-dimensional lifted problem, thereby obtaining upper bounds as opposed to lower bounds. Crucially, to ensure that polynomial densities of arbitrarily low degrees exist for our problem (and therefore the resulting SDP tightenings are feasible), we work with free initial and final measures and set up the cost function and constraints such that this additional freedom does not affect optimality. Importantly, we do not assume that the state constraint set is control invariant, a requirement that is often imposed in the existing literature (e.g., [START_REF] Rantzer | Duality between cost and density in optimal control[END_REF]) but rarely met in practice.

The presented approach bears some similarity with the density approach of [START_REF] Prajna | Nonlinear control synthesis by convex optimization[END_REF] for global stabilization later extended to optimal control (in a purely theoretical setting) in [START_REF] Rantzer | Duality between cost and density in optimal control[END_REF] and recently generalized to optimal stabilization of a given invariant set in [START_REF] Raghunathan | Optimal stabilization using Lyapunov measures[END_REF] (providing both theoretical results and a practical computation method). However, contrary to [START_REF] Prajna | Nonlinear control synthesis by convex optimization[END_REF] we consider the problem of optimal control, not stabilization and moreover we work under state constraints. Contrary to [START_REF] Raghunathan | Optimal stabilization using Lyapunov measures[END_REF] we work in continuous time, consider a more general problem (optimal control, not optimal stabilization of a given set) and our approach of finite-dimensional approximation is completely different in the sense that it is based purely on convex optimization and it does not rely on state-space discretization. Moreover, and importantly, our approach comes with convergence guarantees.

Finally, let us mention that this work is inspired by [START_REF] Lasserre | A New Look at Nonnegativity on Closed Sets and Polynomial Optimization[END_REF], where a converging sequence of upper bounds on static polynomial optimization problems was proposed, as opposed to a converging sequence of lower bounds as originally developed in [START_REF] Lasserre | Global Optimization with Polynomials and the Problem of Moments[END_REF].

Preliminaries

Notation

We use L(X; Y ) to denote the space of all Lebesgue measurable functions defined on a set X ⊂ R n and taking values in the set Y ⊂ R m . If the space Y is not specified it is understood to be R. The spaces of integrable functions and essentially bounded functions are denoted by L 1 (X; Y ) and L ∞ (X; Y ), respectively. The spaces of continuous respectively k-times continuously differentiable functions are denoted by C(X; Y ) respectively C k (X; Y ). By a (Borel) measure we understand a countably additive mapping from (Borel) sets to nonnegative real numbers. Integration of a continuous function v with respect to a measure µ on a set X is denoted by X v(x) dµ(x) or also v dµ when the variable and domain of integration are clear from the context. A probability measure is a measure with unit mass (i.e., 1dµ = 1). The support of a measure µ, defined as the smallest closed set whose complement has zero measure, is denoted by spt µ. The ring of all multivariate polynomials in a variable x is denoted by R[x], the vector space of all polynomials of degree no more than d is denoted by R[x] d , and the vector space of m-dimensional polynomial vectors is denoted by R[x] m . The boundary of a set X is denoted by ∂X, the interior by X • and the closure by X. The Euclidean distance of a point x from a set X is denoted by dist X (x). For a possibly matrix-valued function f ∈ C(X; R n×m ) we define f C 0 (X) := sup x∈X max i,j |f i,j (x)| and for a vector-valued function g

∈ C 1 (X; R n ) we define g C 1 (X) := g C 0 (X) + ∂g ∂x C 0 (X)
, where ∂g ∂x denotes the Jacobian of g. If clear from the context we write • C 0 for • C 0 (X) and similarly for the C 1 norm.

SOS programming

Crucial to the material presented in the paper is the ability to decide whether a polynomial p ∈ R[x] is nonnegative on a set

X = {x ∈ R n | g i (x) ≥ 0, i = 1, . . . , n g }, with g i ∈ R[x].
A sufficient condition for p to be nonnegative on X is that it belongs to the truncated quadratic module of degree d associated to X,

Q d (X) := s 0 + na i=1 g i (x)s i (x) | s 0 ∈ Σ 2 d 2 , s i ∈ Σ 2 (d-deg g i ) 2
, where Σ 2k is the set of all polynomial sum-of-squares (SOS) of degree at most 2k. Note in particular that

Q d+1 (X) ⊃ Q d (X). If p ∈ Q d (X)
for some d ≥ 0 then clearly p is nonnegative on X, and the following fundamental result shows that a certain converse result holds.

Proposition 1 (Putinar [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF])

Let N -x 2 ∈ Q d (X)
for some d > 0 and N ≥ 0 and let p ∈ R[x] be strictly positive on X.Then p ∈ Q d (X) for some d ≥ 0.

Combining with the Stone-Weierstrass Theorem, as an immediate corollary we get:

Corollary 1 Let f ∈ C(X) be nonnegative on X and let N -x 2 ∈ Q d (X)
for some d > 0 and N ≥ 0. Then for every ≥ 0 there exists d ≥ 0 and

p d ∈ Q d (X) such that f -p d C 0 < .
Corollary 1 says that polynomials in Q d (X) are dense (with respect to the C 0 norm) in the space of continuous functions nonnegative on X when we let d tend to infinity.

In the rest of the text we use standard algebraic operations on sets. For instance if we write that p ∈ gQ d (X) + hR[x] d , then it means that p = gq + hr with q ∈ Q d (X) and r ∈ R[x] d .

The inclusion of p ∈ Q d (X) for a fixed d is equivalent to the existence of a positive semidefinite matrix W such that p

(x) = b(x) W b(x), where b(x) is a basis of R[x] d/2
, the vector space of polynomials of degree at most d/2. Comparing coefficients leads to a set of affine constraints on the coefficients of p and the entries of W . Deciding whether p ∈ Q d (X) therefore translates to the feasibility of a semidefinite programming problem with the coefficients of p entering affinely. As a result, optimization of a linear function of the coefficients of p subject to the constraint p ∈ Q d (X) translates to a semidefinite programming problem (SDP) and hence to a well-understood and widely studied class of convex optimization problems for which powerful algorithms and off-the-shelf software are available. See, e.g., [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] and the references therein for more details.

Problem statement

We consider the continuous-time input-affine 1 controlled dynamical system

ẋ(t) = f (x(t)) + m i=1 f u i (x(t))u i (t), (1)
where x ∈ R n is the state, u ∈ R m is the control input, and the data are polynomial:

f ∈ R[x] n , f u i ∈ R[x] n , i = 1, . . . , m.
The system is subject to semi-algebraic state and box 2 input constraints

x(t) ∈ X := {x ∈ R n | g i (x) ≥ 0, i = 1, . . . , n g }, (2a) u(t) ∈ U := [0, ū] m , (2b)
where g ∈ R[x] ng and ū ≥ 0. The set X is assumed compact and the polynomials defining X are assumed to be such that

ḡ(x) := ng i=1 g i (x) > 0 ∀x ∈ X • . ( 3
)
Since X is assumed compact, we also assume, without loss of generality, that the inequalities defining the sets X contain the inequality N -x 2 ≥ 0 for some N ≥ 0.

The goal of the paper is to (approximately) solve the following optimal control problem (OCP):

V (x 0 ) := inf u(•),τ (•) τ (x 0 ) 0 e -βt [l x (x(t)) + m i=1 l u i (x(t))u i (t)] dt + e -βτ M s.t. x(t) = x 0 + t 0 f (x(s)) + m i=1 f u i (x(s))u i (s) ds, (x(t), u(t)) ∈ X × U ∀t ∈ [0, τ (x 0 )] u ∈ L ∞ ([0, τ (x 0 )]; U ), τ ∈ L(X; [0, ∞]) (4)
where β > 0 is a given discount factor and M is a constant chosen such that

M > β -1 sup x∈X,u∈U {l(x, u)}, (5)
where the joint stage cost

l(x, u) := l x (x) + m i=1 l u i (x)u i (6)
is, without loss of generality, assumed to be nonnegative on X × U . The state and input stage cost functions l x and l u i , i = 1, . . . , m, are assumed to be polynomial. The function τ 1 Any dynamical system ẋ = f (x, u) depending nonlinearly on u can be transformed to the input-affine form by using state inflation ẋ u = f (x, u) v , where u is now a part of the state and v a new control input; constraints on v then correspond to rate constraints on u. Similarly, cost functions depending non-linearly on u in problem (4) can be handled using state inflation in exactly the same fashion.

2 Any box can be affinely transformed to [0, ū].

in OCP ( 4) is referred to as a stopping function; the optimization is therefore both over the control input u and over the final time τ (x 0 ), which can be finite or infinite and can depend on the initial condition x 0 .

The function x → V (x) in ( 4) is called the value function. The reason for choosing the slightly non-standard objective function in ( 4) is because with this objective function the value function V is bounded (by M ) on X and it coincides with the standard 3 discounted infinite-horizon value function for all initial conditions x 0 ∈ X for which the trajectories can be kept within the state constraint set X forever using admissible controls, i.e., for all x 0 in the maximum control invariant set associated to the dynamics (1) and the constraints (2).

To see the first claim, set τ (x 0 ) = 0 for all x 0 ∈ X. To see the second claim notice that with M chosen as in [START_REF] Henrion | Nonlinear optimal control synthesis via occupation measures[END_REF], it is always beneficial to continue the time evolution whenever possible and therefore τ (x 0 ) = +∞ for all x 0 in the maximum controlled invariant set associated to (1) and ( 2).

Remark 1 A constant M satisfying ( 5) can be found either by analytically evaluating the supremum in [START_REF] Henrion | Nonlinear optimal control synthesis via occupation measures[END_REF] or by using the techniques of [START_REF] Lasserre | Global Optimization with Polynomials and the Problem of Moments[END_REF] to find an upper bound.

Given a Lipschitz continuous feedback controller u ∈ C(X; U ) and a stopping function τ ∈ L(X; [0, ∞]), the ODE (1) has a unique solution and we let V u,τ ∈ L(X; [0, ∞]) denote the value function attained by (u, τ ) in ( 4), i.e., setting u(t) = u(x(t)). By V u we denote the value function V u,τ u , where τ u ∈ L(X; [0, ∞]) is the optimal stopping function associated to u. Note that, by the choice of M in ( 5), the optimal stopping function τ u is equal to the first hitting time of the complement of the constraint set X, i.e.,

τ u (x 0 ) = inf{t ≥ 0 | x(t | x 0 ) / ∈ X},
where x(t | x 0 ) is the trajectory of (1) with u(t) = u(x(t)) starting from x 0 . Notice also that V u,τ (x) ≥ V (x) for all x ∈ X and that for any pair (u, τ ) feasible in (4) we have V u,τ (x) ≤ M for all x ∈ X.

Throughout the paper, we make the following technical assumption:

Assumption 1 There exists a sequence of Lipschitz continuous feedback controllers

{u k ∈ C(X; U )} ∞ k=1 and stopping functions {τ k ∈ L(X; [0, ∞])} ∞ k=1 feasible in (4) such that lim k→∞ X (V u k ,τ k (x) -V (x))dx = 0 (7)
and such that for every k ≥ 0 there exist a function

ρ k ∈ C 1 (X) and a scalar γ k > 0 such that ρ k (x) = 0 if dist ∂X (x) < γ k and X τ k (x 0 ) 0 e -βt v(x k (t | x 0 )) dt dx 0 = X v(x)ρ k (x) dx ∀v ∈ C(X), (8)
where x k (• | x 0 ) denotes the solution to (1) controlled by u k .

3 By standard we mean a discounted optimal control problem with cost

∞ 0 e -βt [l x (x(t)) + m i=1 l ui (x(t))u i (t)
] dt and no stopping function.

Remark 2 Note that V u k ,τ k ≥ V on X by construction and therefore ( 7) is equivalent to the

L 1 convergence of V u k ,τ k to V .
Assumption 1 says that the optimal control inputs and stopping functions for OCP (4) can be well approximated by Lipschitz continuous feedback controllers and measurable stopping functions such that the resulting densities of the discounted occupation measures are continuously differentiable and vanish near the boundary of X. Note that the existence of an optimal feedback controller, as well as whether it can be well approximated by Lipschitz controllers, are subtle issues. Similarly it is a subtle issue whether asymptotically optimal stopping functions can be found such that the associated densities ρ k in ( 8) are continuously differentiable and vanish near the boundary of X (note, however, that the left hand side of ( 8) can always be represented as X v(x)dµ k (x) for some nonnegative measure µ k ). This problem is of rather technical nature and has been studied in the literature (e.g., [3, Section 1.4] or [START_REF] Rantzer | An converse theorem for density functions[END_REF]), where affirmative results have been established in related settings. We do not undertake a study of this problem here and rely on Assumption 1, which is, for ease of reading, not stated in its most general form. For example, the functions ρ k do not need to be C 1 but only weakly differentiable and the integration on the left-hand side of ( 8) can be weighted by a nonnegative function

ρ k 0 ∈ L 1 (X) satisfying ρ k 0 ≥ 1 on X and ρ k 0 → 1 in L 1 (X).
In addition, we conjecture that it is enough to require ρ k = 0 on ∂X and not necessarily on some neighborhood of ∂X; this is in particular the case when X is a box or a ball but we expect all the results of the paper to hold with a general semialgebraic set for which the defining functions satisfy (3).

The main result of this paper is a hierarchy of sum-of-squares (SOS) problems providing an explicit sequence of rational feedback controllers u k ∈ C ∞ (X; U ) such that, under Assumption 1, (7) holds with τ k = τ u k , i.e., a sequence of asymptotically optimal rational controllers in the sense of the L 1 convergence of the associated value functions (see Remark 2).

Converging hierarchy of solutions

In this section we present an infinite-dimensional linear program (LP) in the space of continuous functions whose sum-of-squares (SOS) approximations provide a sequence of rational controllers u k satisfying [START_REF] Korda | Convex computation of the maximum controlled invariant set for polynomial control systems[END_REF]. This infinite-dimensional LP is closely related (and in a weak sense equivalent) to OCP (4); the rationale behind the derivation of the LP and its relation to OCP (4) is detailed in Section 5.

The infinite-dimensional LP reads inf

ρ, ρ 0 , ρ T , σ X l x (x)ρ(x) dx + m i=1 X l u i (x)σ i (x) dx + M X ρ T (x) dx s.t. ρ T -ρ 0 + βρ + div(ρf ) + m i=1 div(σ i f u i ) = 0 ρ ≤ 0 on ∂X ρ 0 ≥ 1 on X ūρ ≥ σ i on X, i = 1, . . . , m. ρ T ≥ 0 on X σ i ≥ 0 on X, i = 1, . . . , m. (9)
The optimization in ( 9) is over functions (ρ,

ρ 0 , ρ T , σ) ∈ C 1 (X) × C(X) × C(X) × C 1 (X) m with σ = (σ 1 , . . . , σ m ).
The optimal value of (9) will be denoted by p . The value attained in ( 9) by any tuple of densities (ρ, ρ 0 , ρ T , σ) feasible in (9) will be denoted by p(ρ, ρ 0 , ρ T , σ).

Remark 3 (Non-uniform weighting) Note that we could have imposed ρ 0 ≥ ρ0 for any polynomial ρ0 nonnegative on X. Choosing a different ρ0 has no impact on the asymptotic convergence of the value functions established in the rest of the paper as long as ρ0 is strictly positive on X. It may, however, influence the speed of convergence in different subsets of X.

In general we expect faster convergence where ρ0 is large and slower convergence where it is small. Choosing a non-constant ρ0 therefore allows to assign a different importance to different subsets of X.

The infinite-dimensional LP ( 9) is then approximated by a hierarchy of sum-of-squares (SOS) problems, which immediately translate to finite-dimensional semidefinite programs (SDPs).

The SOS approximation of degree d of (9) reads inf

(ρ,ρ 0 ,ρ T ,σ)∈R[x] 3+m d X l x (x)ρ(x) dx + m i=1 X l u i (x)σ i (x) dx + M X ρ T (x) dx s.t. ρ T -ρ 0 + βρ + div(ρf ) + m i=1 div(σ i f u i ) = 0 -ρ ∈ Q d (X) + g i R[x] d-deg g i + ḡR[x] d-deg ḡ i = 1, . . . , n g ρ 0 -1 ∈ Q d (X) ūρ -σ i ∈ Q d (X) + ḡQ d-deg ḡ(X ) i = 1, . . . , m ρ T ∈ Q d (X) σ i ∈ Q d (X) + ḡQ d-deg ḡ(X ), i = 1, . . . , m. (10)
Once a basis for R[x] d is fixed (e.g., the standard monomial basis), the objective becomes linear in the coefficients of polynomials ρ, σ and ρ T , and the equality constraint is imposed by equating the coefficients. The inclusions in the quadratic modules translate to semidefinite constraints and affine equality constraints; see Section 2.2. Optimization problem [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] therefore immediately translates to an SDP.

Remark 4 (Feasibility) Trivially, any feasible solution to [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] is feasible in [START_REF] Korda | Controller design and region of attraction estimation for nonlinear dynamical systems[END_REF]. Also, problem [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] is feasible for any d ≥ 0. Indeed (ρ, ρ 0 , ρ T , σ) = (0, 1, 1, 0) is always feasible in [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]. See also Remark 6 below.

If non-uniform weighting of initial conditions (see Remark 3) was required, the constraint ρ 0 -1 ∈ Q d (X) would be replaced by ρ 0 -ρ0 ∈ Q d (X) for a polynomial weighting function ρ0 nonnegative on X.

Given an optimal solution (ρ d , ρ d 0 , ρ d T , σ d ) to [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF], we define a rational control law u d by

u d i (x) := σ d i (x) ρ d (x) ∀x ∈ X, i = 1, . . . , m. (11)
The main result of the paper is the following theorem stating that the controllers u d are asymptotically optimal:

Theorem 1 For all d ≥ 0 we have u d (x) ∈ U for all x ∈ X and if Assumption 1 holds, then

lim d→∞ X (V u d (x) -V (x)) dx = 0, (12)
that is,

V u d → V in L 1 (X) (note that V u d ≥ V on X).
5 Rationale behind the LP formulation ( 9) and proof of the main Theorem 1

This section explains the rationale behind the LP problem [START_REF] Korda | Controller design and region of attraction estimation for nonlinear dynamical systems[END_REF] and its relation to the OCP (4) and gives the proof of Theorem 1. First, we lift the original problem (4) into the space of measures with nonnegative densities in C(X); this lifting is problem [START_REF] Korda | Controller design and region of attraction estimation for nonlinear dynamical systems[END_REF]. Next we tighten the problem by considering only polynomials of prescribed degree and with nonnegativity constraints enforced via SOS conditions; this is problem [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]. Importantly, the lifting ( 9) is a tightening of the original problem (4) as show in Theorem 3 below. This is in contrast with [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF] where the original problem was lifted into the space of measures and this lifting was a relaxation.

To be more concrete, observe that any initial measure µ 0 , stopping function τ ∈ L(X; [0, ∞]), and family of trajectories {x(• | x 0 )} x 0 ∈X of (1) generated by a Lipschitz controller u ∈ C(X; U ) give rise to a triplet of measures defined by

X v(x)dµ(x) = X τ (x 0 ) 0 e -βt v(x(t | x 0 )) dt dµ 0 (x 0 ), ( 13a
) X v(x)dµ T (x) = X e -βτ (x 0 ) v(x(τ (x 0 ) | x 0 )) dµ 0 (x 0 ), ( 13b
) X v(x)dν i (x) = X τ (x 0 ) 0 e -βt v(x(t | x 0 ))u i (x(t | x 0 )) dt dµ 0 (x 0 ). ( 13c
)
The measure µ is called discounted occupation measure, the measure µ T terminal measure and the measures ν i , i = 1, . . . , m, control measures. These measures satisfy the discounted Liouville equation

X v dµ T (x) = X v dµ 0 (x) + X (∇v • f -βv) dµ(x) + m i=1 X ∇v • f u i dν i (x) (14)
for all v ∈ C 1 (X). This follows by direct computation; see, e.g., [START_REF] Korda | Controller design and region of attraction estimation for nonlinear dynamical systems[END_REF]. Notice also that dν i (x) = u i (x)dµ(x), i.e., ν i is absolutely continuous with respect to µ with Radon-Nikodým derivative equal to u i .

Crucially, the converse statement is also true, although we have to go from stopping functions to stopping measures:

Theorem 2 (Superposition) If measures µ, µ 0 , µ T and ν i , i = 1, . . . , m, satisfy ( 14) with spt µ 0 ⊂ X, spt µ ⊂ X and spt µ T ⊂ X and dν i = u i dµ for some Lipschitz u ∈ C(X, U ), then there exists an ensemble of probability measures (i.e., measures with unit mass) {τ x 0 } x 0 ∈X and an ensemble of trajectories {x(• | x 0 )} x 0 ∈X of the system (1) controlled with u(t) = u(x(t)) such that x(t | x 0 ) ∈ X for all t ∈ spt τ x 0 and

X v(x) dµ 0 (x) = X v(x(0 | x 0 )) dµ 0 (x 0 ), ( 15a
) X v(x) dµ(x) = X ∞ 0 τ 0 e -βt v(x(t | x 0 )) dt dτ x 0 (τ ) dµ 0 (x 0 ), ( 15b
) X v(x) dµ T (x) = X ∞ 0 e -βτ v(τ (x 0 )) dτ x 0 (τ ) dµ 0 (x 0 ), ( 15c
) X v(x) dν i (x) = X ∞ 0 τ 0 e -βt v(x(t | x 0 ))u i (x(t | x 0 )) dt dτ x 0 (τ ) dµ 0 (x 0 ) (15d)
for all v ∈ C 1 (X).

Proof: See Appendix B.

Remark 5 (Interpretation of Theorem 2) Theorem 2 says that any measures satisfying ( 14) are generated by a superposition of the trajectories of the dynamical system ẋ =

f (x) + m i=1 f u i (x)u i (x)
, where the superposition is over the final time of the trajectories. Note that there is a unique trajectory corresponding to each initial condition (since the vector field f (x)+ m i=1 f u i (x)u i (x) is Lipschitz) but this unique trajectory can be stopped at multiple times (in fact at a whole continuum of times) allowing for superposition; this superposition is captured by the stopping measures {τ x 0 } x 0 ∈X . For example, if the τ x 0 is a Dirac measure at a given time, then there is no superposition; if τ x 0 has a discrete distribution, then there is a superposition of finitely or countably many overlapping trajectories starting at x 0 stopped at different time instances; if τ x 0 has a continuous distribution then there is a superposition of a continuum of overlapping trajectories starting from x 0 stopped at different time instances.

If in addition the measures µ 0 , µ, µ T satisfying the discounted Liouville equation ( 14) are absolutely continuous with respect to the Lebesgue measure with densities ρ 0 ∈ C(X), ρ ∈ C 1 (X), ρ T ∈ C(X) such that ρ = 0 on ∂X, then these densities satisfy

ρ T -ρ 0 + βρ + div(f ρ) + m i=1 div(f u i σ i ) = 0 ( 16
)
with σ i = u i ρ, i = 1, . . . , m. This follows directly by substituting dµ 0 = ρ 0 dx, dµ = ρdx, dµ T = ρ T dx and dν i = u i dµ = u i ρdx = σ i dx in ( 14) and using integration by parts. Equation ( 16) holds almost everywhere in X with a Lipschitz controller u, since Lipschitz functions are differentiable almost everywhere and the integration by parts formula applies to them, and everywhere with u ∈ C 1 (X; U ).

Remark 6 (Role of the terminal measure) An important feature of the SOS tightenings [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] is that they are feasible for arbitrarily low degrees (see Remark 4), which is crucial from a practical point of view and is not satisfied with other, more obvious, formulations (e.g., those not involving a stopping function in (4)); the reason for this is that, in the absence of a terminal measure (i.e., ρ T = 0), the discounted Liouville equation ( 16) may not have a solution with a polynomial ρ even though ρ 0 and the dynamics are polynomial. Indeed, for example with f = -x, f u i = 0, β = 1, ρ 0 = 1 on X = [-1, 1] and zero elsewhere, the only solution to ( 16) with ρ T = 0 is ρ(x) = -ln(|x|).

Theorem 2 immediately enables us to prove a representation of the cost of problem ( 9) in terms of trajectories of (1).

Lemma 1 If (ρ, ρ 0 , ρ T , σ) is feasible in ( 9) and u = σ/ρ, then

p(ρ, ρ 0 , ρ T , σ) = X ∞ 0 τ 0 e -βt l x (x(t | x 0 ))dt dτ x 0 (τ )ρ 0 (x 0 ) dx 0 + m i=1 X ∞ 0 τ 0 e -βt l u i (x(t | x 0 ))u i (x(t | x 0 )) dt dτ x 0 (τ )ρ 0 (x 0 )dx 0 + M X ∞ 0 e -βτ dτ x 0 (τ )ρ 0 (x 0 )dx 0 , (17)
where x(• | x 0 ) are trajectories of (1) controlled by u(t) = u(x(t)) and τ x 0 are stopping probability measures with support spt

τ x 0 included in [0, ∞]. Moreover the state-control trajectories x(• | x 0 ) and u(x(• | x 0 )) are feasible in (4) in the sense that x(t | x 0 ) ∈ X and u(t | x 0 ) ∈ U for all t ∈ spt τ x 0 .
Proof: Let (ρ, ρ 0 , ρ T , σ) be feasible in [START_REF] Korda | Controller design and region of attraction estimation for nonlinear dynamical systems[END_REF] and let p(ρ, ρ 0 , ρ T , σ) denote the value attained by (ρ, ρ 0 , ρ T , σ) in [START_REF] Korda | Controller design and region of attraction estimation for nonlinear dynamical systems[END_REF]. The equality constraint of ( 9) is exactly [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF]. Since the constraint of ( 9) implies ρ = 0 on ∂X, equation ( 14) holds with dµ 0 = ρ 0 dx, dµ = ρdx, dµ T = ρ T dx and dν i = u i dµ = u i ρdx = σ i dx, where u i = σ i ρ ∈ C 1 (X; U ), i = 1, . . . , m. By Theorem 2 (setting v(x) = l x (x) in (15b), v(x) = 1 in (15c) and v(x) = l u i (x) in (15d)) we obtain the result (noticing that the constraints of (9) imply that u(x) ∈ U for all x ∈ X).

Corollary 2 If (ρ, ρ 0 , ρ T , σ) is feasible in [START_REF] Korda | Controller design and region of attraction estimation for nonlinear dynamical systems[END_REF] and u = σ/ρ, then

p(ρ, ρ 0 , ρ T , σ) ≥ X V u (x 0 )ρ 0 (x 0 )dx 0 . ( 18
)
If in addition the stopping measures {τ x 0 } x 0 ∈X in ( 17) are equal to the Dirac measures

{δ τ (x 0 ) } x 0 ∈X for some stopping function τ ∈ L(X; [0, ∞]), then p(ρ, ρ 0 , ρ T , σ) = X V u,τ (x 0 )ρ 0 (x 0 )dx 0 . ( 19
)
Proof: Let (ρ, ρ 0 , ρ T , σ) be feasible in [START_REF] Korda | Controller design and region of attraction estimation for nonlinear dynamical systems[END_REF]. Using Lemma 1, p(ρ, ρ 0 , ρ T , σ) has representation [START_REF] Seiler | SOSOPT: A toolbox for polynomial optimization[END_REF], where the state-control trajectories in [START_REF] Seiler | SOSOPT: A toolbox for polynomial optimization[END_REF] are feasible in [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF]. Since the measures τ x 0 in [START_REF] Seiler | SOSOPT: A toolbox for polynomial optimization[END_REF] have unit mass for all x 0 ∈ X, we obtain [START_REF] Rantzer | An converse theorem for density functions[END_REF]. If τ x 0 = δ τ (x 0 ) for some stopping function τ ∈ L(X; [0, ∞]), then the integrals with respect to τ x 0 in ( 17) become evaluations at τ (x 0 ) and hence [START_REF] Rantzer | Duality between cost and density in optimal control[END_REF] holds.

Corollary 2 immediately implies that the problem (9) (and hence problem [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]) is a tightening of the original problem (4):

Proof: Since ḡ > 0 on X • , we can factor ρ = ḡh with h ∈ C 1 (X) given by

h(x) := ρ(x)/ḡ(x) if dist ∂X (x) ≥ ζ 0 otherwise.
Since polynomials are dense in C 1 there exists for every δ > 0 a polynomial ĥ > 0 such that ĥ

-h C 1 (X) < δ. ( 22
)
Applying Proposition 1 to ĥ we see that there exists p d ∈ Q d(X ) for some d ≥ 0 such that ĥ

-p d C 1 (X) < δ. ( 23
ρ -p d C 0 = hḡ -p d ḡ C 0 ≤ ḡ C 0 h -p d C 0 < 2δ ḡ C 0 and ∇ρ -∇p d C 0 = ḡ ∇h + h ∇ḡ -ḡ ∇p d + p d ∇ḡ C 0 ≤ ∇ḡ C 0 h -p d C 0 + ḡ C 0 ∇h -∇p d C 0 ≤ 2δ ∇ḡ C 0 + ḡ C 0 .
Therefore choosing δ such that 2δ ∇ḡ C 0 + 2 ḡ C 0 < gives the desired result. Now we are ready to prove our main result, Theorem 1.

Proof (of Theorem 1): Consider the sequences

{u k ∈ C 1 (X; U )} ∞ k=1 , {τ k ∈ L(X; [0, ∞])} ∞ k=1
from Assumption 1. By the first part of Lemma 2 the sequence of associated densities (ρ k , ρ k 0 , ρ k T , σ k ) generated by (u k , τ k ) is feasible in [START_REF] Korda | Controller design and region of attraction estimation for nonlinear dynamical systems[END_REF] and satisfies [START_REF] Prajna | Nonlinear control synthesis by convex optimization[END_REF]. By Assumption 1,

ρ k = 0 and σ k = 0 on {x ∈ X : dist ∂X (x) < γ k } (since σ k = u k ρ k ) with γ k > 0.
Hence by Lemma 3 there exist polynomial densities ρ k,pol ∈ ḡQ

d k -deg ḡ(X ), σ k,pol ∈ ḡQ d k -deg ḡ(X ) m for some degrees d k ≥ 0 such that ρ k -ρ k,pol C 1 (X) < 1/k ( 24
)
σ k i -σ k,pol i C 1 (X) < 1/k (25) ūρ k,pol -σ k,pol i ∈ ḡQ d k -deg ḡ(X ) for all i = 1, . . . , m (since u k (x) ∈ U = [0, ū] m for all x ∈ X and hence ūρ k ≥ σ k i on X). Notice also that since ρ k,pol ∈ ḡQ d k -deg ḡ(X ), we have -ρ k,pol ∈ ḡR d k -deg ḡ. Next, since ρ k
0 ≥ 1 and ρ k T ≥ 0, we can find, by Corollary 1, polynomial densities ρk,pol

0 ∈ 1 + Q d k (X) and ρk,pol T ∈ Q d k (X) such that ρ k 0 -ρk,pol 0 C 0 (X) < 1/k, ( 26
)
ρ k T -ρk,pol T C 0 (X) < 1/k. ( 27
)
Since (ρ k , ρ k 0 , ρ k T , σ k ) satisfy the equality constraint of ( 9) we have ρk,pol

T + βρ k,pol -ρk,pol 0 + div(ρ k,pol f ) + m i=1 div(σ k,pol i f u i ) = ω k
where

ω k := ρk,pol T -ρ k T + β(ρ k,pol -ρ k ) -(ρ k,pol 0 -ρ k 0 ) + div[(ρ k,pol -ρ k )f ] + m i=1 div[(σ k,pol i -σ k i )f u i ]
is a polynomial such that ω k C 0 → 0 as k → ∞ in view of ( 24)-( 27). Defining the constants k = 1/k + ω k C 0 and setting ρ k,pol T := ρk,pol

T + k ρ k,pol 0 := ρk,pol 0 + k + ω k
we see that

ρ k,pol T + βρ k,pol -ρ k,pol 0 + div(ρ k,pol f ) + m i=1 div(σ k,pol i f u i ) = 0,
and ρ k,pol 0 -1 and ρ k,pol T are strictly positive on X and hence belong to Q d k (X). The densities (ρ k,pol , ρ k,pol 0 , ρ k,pol T , σ k,pol ) are therefore feasible in [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] for some d k ≥ 0. In addition, by construction, ρ k,pol

0 -ρ k 0 C 0 → 0 and ρ k,pol T -ρ k T C 0 → 0 as k → ∞.
Therefore we have obtained a sequence of polynomial densities (ρ k,pol , ρ k,pol 0 , ρ k,pol T , σ k,pol ) that are feasible in [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] and such that

ρ k,pol 0 -ρ k 0 C 0 → 0, ρ k,pol T -ρ k T C 0 → 0, ρ k,pol -ρ k C 1 → 0, σ k,pol -σ k C 1 → 0 as k → ∞. This implies that |p(ρ k,pol , ρ k,pol 0 , ρ k,pol T , σ k,pol ) -p(ρ k , ρ k 0 , ρ k T , σ k )| → 0
and hence (ρ k,pol , ρ k,pol 0 , ρ k,pol T , σ k,pol ) satisfies [START_REF] Prajna | Nonlinear control synthesis by convex optimization[END_REF] and so p(ρ k,pol , ρ k,pol 0 , ρ k,pol T , σ k,pol ) → p by Theorem 3. Therefore [START_REF] Korda | Convex computation of the maximum controlled invariant set for polynomial control systems[END_REF] holds with the rational controllers u k := σ k,pol /ρ k,pol by the second part of Lemma 2. This finishes the proof.

Value function approximations

In this section we propose a converging hierarchy of approximations from below and from above to the value function V u associated to a rational controller u = σ/ρ with σ ∈ R[x] m and ρ ∈ R[x] satisfying 0 ≤ σ i ≤ ūρ on X. In addition we describe a hierarchy of approximations from below to the optimal value function V . This is useful as a post-processing step, once a rational control law has been computed as described in Section 4, providing an explicit bound on the suboptimality of the controller.

Remark 7

The question whether V converges from below to V as degree d in (31) tends to infinity is open (although likely to hold). A proof would require an extension of the superposition Theorem 7 to non-Lipschitz vector fields (in the spirit of the finite-time superposition result of [START_REF] Ambrosio | Existence, uniqueness, stability and differentiability properties of the flow associated to weakly differentiable vector fields[END_REF]Theorem 4.4]) or an extension of the argument of [START_REF] Gaitsgory | Linear programming approach to deterministic infinite horizon optimal control problems with discounting[END_REF] to the case of µ T = 0, either of which is beyond the scope of this paper.

Remark 8 Besides closed-loop cost function with respect to the OCP (4), one can assess other aspects of the closed-loop behavior of the dynamical system (1) controlled by the rational controller u = σ/ρ. In particular, regions of attraction or maximum controlled invariant sets can be estimated by methods of [START_REF] Henrion | Convex computation of the region of attraction of polynomial control systems[END_REF][START_REF] Korda | Controller design and region of attraction estimation for nonlinear dynamical systems[END_REF][START_REF] Korda | Convex computation of the maximum controlled invariant set for polynomial control systems[END_REF], which extend readily to the case of rational systems.

Numerical examples

This section demonstrates the approach on numerical examples. To improve the numerical conditioning of the SDPs solved, we use the Chebyshev basis to parametrize all polynomials. More specifically, we use tensor products of univariate Chebyshev polynomials of the first kind to obtain a multivariate Chebyshev basis. We note, however, that similar results, albeit slightly less accurate could be obtained with the standard multivariate monomial basis (in which case the SDPs can be readily formulated using high level modelling tools such as Yalmip [START_REF] Löfberg | YALMIP : A toolbox for modeling and optimization in MATLAB[END_REF] or SOSOPT [START_REF] Seiler | SOSOPT: A toolbox for polynomial optimization[END_REF]). The resulting SDPs were solved using MOSEK.

Nonlinear double integrator

As our first example we consider the nonlinear double integrator ẋ1 = x 2 + 0.1x 3 1 ẋ2 = 0.3u subject to the constraints u ∈ [-1, 1] and x ∈ X := {x : x 2 < 1} and stage costs l x (x) = x x and l u (x) = 0. The discount factor β was set to 1; the constant M to 1.01 > sup x∈X {x x}/β = 1. First we obtain a rational controller of degree six by solving ( 9) with d = 6. The graph of the controller is shown in Figure 1. Next we obtain a polynomial upper bound V u of degree 14 on the value function associated to u by solving (28) with d = 14. To assess suboptimality of the controller u we compare it with a lower bound V on the optimal value function of the problem (4) obtained by solving (31) with d = 14. The graphs of the two value functions are plotted in Figure 2. We see that the gap between the upper bound on V u and lower bound on V is relatively small, verifying a good performance of the extracted controller. Quantitatively, the average performance gap defined as 100 X (V u -V )dx/ X V dx is equal to 19.5%.

x 1 x 2 u(x)

Controlled Lotka-Volterra

In our second example we apply the proposed method to a population model governed by n-dimensional controlled Lotka-Volterra equations

ẋ = r • x • (e -Ax) + u + -u -,
where e ∈ R n is the vector of ones and • denotes the componentwise (Hadamard) product. Each component x i of the state x ∈ R n represents the size of the population of species i. The vector r ∈ R n contains the intrinsic growth rates of each species and the matrix A ∈ R n×n captures the interaction between the species. If A i,j > 0, then species j is harmful to species i (e.g., competes for resources) and if A i,j < 0, then species j is helpful to species i (e.g., species i feeds on species j); the diagonal components A i,i are normalized to one. The control inputs u + ∈ [0, 1] n and u -∈ [0, 1] n represent, respectively, the inflow and outflow of new species from the outside. For our numerical example we select n = 4 and model parameters

r =     1 0.6 0.4 0.2     , A =     1 0.3 0.4 0.2 -0.2 1 0.4 -0.1 -0.1 -0.2 1 0.3 -0.1 -0.2 -0.3 1     ,
which results in a system with four states and eight control inputs. The economic objective is to harvest species number one while ensuring that no species goes extinct. More specifically the cost function is l u (x) = (-1.0, 0.5, 0.6, 0.8, 1.1, 2, 4, 6) and l x (x) = 1, where the vector l u (x) is associated with the control input vector u = (u -, u + ). Therefore there is a reward for harvesting species number one and cost associated with both introduction and hunting of all other species, the cost of hunting being lower than the cost of introduction. The reason for choosing l x (x) = 1 is in order to make the joint stage cost l(x, u) (6) nonnegative; this choice does not affect optimality since l x (x(t)) = 1 irrespective of the control input applied. The non-extinction constraint is expressed as g

(x) = 1 -(Q -1 x -q) (Q -1 x - q) ≥ 0 with Q = diag(0.475 • e
) and q = 0.525 • e. We choose β = 1 and M = 16.16 > sup x∈X,u∈U {l(u, x)}/β = 16. We apply the coordinate transformation x = Qx + q and solve solve obtain a rational controller of degree eight by solving [START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF]. Figure 3 we shows plots for two different initial conditions, one with low population size of the first species and one with high. Finally, we evaluate the suboptimality of the extracted controller using the polynomial lower bound on the optimal value function of degree 11 obtained from (31). Using Monte Carlo simulation with 1000 samples of initial conditions drawn from a uniform distribution over the constraint set we obtain average cost of the extracted controller to be 0.89 whereas the lower bounds predicts average cost of 0.72; hence the extracted controller is no more than 23.6 % suboptimal (modulo the statistical estimation error). Note that we could also obtain a deterministic suboptimality estimate using the upper bound on the value function of the extracted controller obtained from (28). In this case, however, the upper bound (28) is not informative. Nevertheless, the Monte Carlo simulation along with the lower bound (31) is a viable alternative in this case, since the extracted controller is simple and hence trajectories of the controlled system can be simulated rapidly.

Conclusion

This paper presented a method to obtain a sequence of rational controllers asymptotically optimal (under suitable technical assumptions) in a discounted optimal control problem and a method to explicitly estimate suboptimality of each controller. The rational controller of a given degree is obtained by solving a single sum-of-squares problem with no extraction step.

The SOS problem solved is feasible for any degree and therefore this method allows to trade off complexity of the controller against performance.

The approach is based on lifting the nonconvex optimal control problem into an infinite dimensional space of measures with continuous densities, where this problem becomes linear. Crucially, this problem is a tightening of the original problem, which follows immediately from the representation result for solutions of the discounted Liouville's equation with a terminal measure (Theorem 7). Asymptotic optimality of the extracted controllers then follows by approximating the asymptotically optimal continuous densities (guaranteed to exist by Assumption 1) with polynomial densities in such a way that these densities correspond to the densities of the dynamical system (this is the essence of the proof of Theorem 1).

Appendix A

This Appendix contains the proof of Theorem 4; we use the same notation as in Section 6.

The inequalities V u d ≥ V u ≥ V u d follow from Gronwall's Lemma by noticing that the constraints of (28) and (29) imply that

∇V u d • (f + m i=1 f u i u i ) ≤ βV u d -(l x + m i=1 l u i u i ), (32)
∇V u d • (f + m i=1 f u i u i ) ≥ βV u d -(l x + m i=1 l u i u i ) (33)
on X and V u ≥ M , V u ≤ M on ∂X. We detail the argument for the inequality V u d ≥ V u , the inequality V u ≥ V u d being similar. Given x 0 ∈ X the inequality (32) implies that

d dt V u d (x(t | x 0 )) ≤ βV u d (x(t | x 0 )) -l x (x(t | x 0 )) + m i=1 l u i (x(t | x 0 ))u i (x(t | x 0 )) ,
and therefore by Gronwall's Lemma

V u d (x(t | x 0 )) ≤ e βt V u d (x 0 ) - t 0 e β(t-s) l x (x(s | x 0 )) + m i=1 l u i (x(s | x 0 ))u i (x(s | x 0 )) ds computation then gives: ∇v • f (x(t | x 0 )) = d dt v(x(t | x 0 )) = - d dt ∞ 0 e -βs w(x(s | x(t | x 0 ))) ds = - d dt ∞ 0 e -βs w(x(t + s | x 0 )) ds = - ∞ 0 e -βs ∇w(x(t + s | x 0 )) • f (x(t + s | x 0 )) ds = - ∞ 0 e -βs ∇w(t + s | x 0 )) • f (x(t + s | x 0 )) ds = - ∞ 0 e -βs d ds w(x(t + s | x 0 )) ds = -β ∞ 0 e -βs w(x(t + s | x 0 )) ds -[e -βs w(x(t + s | x 0 ))] ∞ 0 = -β ∞ 0 e -βs w(x(s | x(t | x 0 ))) ds + w(x(t | x 0 )) = βv(x(t | x 0 )) + w(x(t | x 0 )).
Setting t = 0, we arrive at (37).

Proof (of Theorem 7) We will proceed in several steps.

Two Diracs. We start with the simplest case of µ 0 = δ x 0 and µ T = aδ x T , a > 0, x T ∈ R n , and some µ ≥ 0. First we will show that if (µ T , µ 0 , µ) solves (35) then there exists a time τ ≥ 0 such that x(τ | x 0 ) = x T . Consider now any w ∈ C 1 c (R n ), w ≥ 0 and the associated v ∈ C 1 (R n ) solving (37). Then we have

av(x T ) -v(x 0 ) = R n (∇v • f -βv) dµ = R n w dµ ≥ 0. Therefore, by Lemma 4, av(x T ) ≥ v(x 0 ) = -
Now pick S ≥ 0 (to be specified later) and consider the traces

X 0 = {x(t | x 0 ) | 0 ≤ t ≤ S}. X T = {x(t | x T ) | 0 ≤ t ≤ S}.
Assuming there is no τ ≥ 0 such that x(τ | x 0 ) = x T we have X 0 ∩ X T = ∅ and since X 0 and X T are compact there exist (by Uryshon's Lemma with mollification) a function w ∈ C 1 c (R d ; [0, 1]) such that w = 0 on X 0 and w = 1 on X T . Then the left hand side of (39) is greater than or equal to a(1 -e -βS )/β whereas the right hand side is less than or equal to e -βS /β. Since a > 0 and β > 0 we arrive at a contradiction with (39) by picking a sufficiently large S. Therefore there exists a τ ≥ 0 such that x(τ | x 0 ) = x T (i.e., x T and x 0 are on the same trace of the flow associated to ẋ = f (x)). Now we prove that a ≤ e -βτ . Since x T = x(τ ) and x 0 are on the same trace we have (40)

Consider the set

X τ = {x(t | x 0 ) | 0 ≤ t ≤ τ }.
Since x 0 and x T are on the same trace (and x T follows x 0 ) there exists w ∈ C 1 c (X), w ≥ 0, such that w = 0 on X τ and w > 0 elsewhere (e.g., let w(x) = min(dist(x, X τ ), 1) with appropriate mollification). With this choice of w the equation ( 40 i.e., µ is indeed generated by trajectories of ẋ = f (x) (in this case by two trajectories, both starting at x 0 , one stopping at τ , the other one continuing to infinity with weights given by the ratio of masses of µ 0 and µ T ). That is the measure τ x 0 is given by τ x 0 = ae βτ δ τ + (1 -ae βτ )δ ∞ as expected.

Dirac at x 0 , sum of Diracs for µ T . Next we treat the case where µ T = ∞ i=1 a i δ x i for some a i ≥ 0 and x i ∈ R n . Using the same argument as in the previous case we can show that

x i ∈ X 0 = {x(t | x 0 ) | t ≥ 0}
for all i and that the condition

∞ i=1
a i e βτ i ≤ 1, holds with τ i being the times to reach x i from x 0 . Then we have

τ x 0 = ∞ i=1 a i e βτ i δ τ i + (1 - ∞ i=1
a i e βτ i )δ ∞ .

Dirac at x 0 arbitrary µ T . In the same way as before we can show that the support of µ T must be on the trace X 0 . Then we can define the measure τx 0 by τx 0 (A) := µ T (x(A | x 0 )), A ⊂ [0, ∞)

and show that it has to satisfy the condition ∞ 0 e βt dτ x 0 (t) ≤ 1. Next, using the fact that the mapping t → x(t | x 0 ) is invertible, we obtain

R n v dµ T = ∞ 0 v(x(t | x 0 )) dτ x 0 (t).
The conclusion of the theorem then holds with τ x 0 defined by Arbitrary µ 0 , arbitrary µ T . The general case follows by approximating µ 0 by a sum of Dirac measures, using the fact that any measure is the weak limit of a sequence of Dirac measures.

τ x 0 (A) =

)

  Defining p d := p d ḡ we see that p d ∈ ḡQ d-deg ḡ(X ) with d = d + deg(ḡ) and that p d = 0 on ∂X. Finally,

Figure 1 :

 1 Figure 1: Nonlinear double integrator -rational controller of degree six.

x 1 x 2 Figure 2 :

 22 Figure 2: Nonlinear double integrator -upper bound on the value function V u associated to the extracted controller (red); lower bound on the optimal value function V (blue).

Figure 3 :

 3 Figure 3: Controlled Lotka-Volterra -(blue) trajectory starting from a high initial population of the first species and low initial population of the other species; (red) trajectory starting from low initial population of the first species and high initial population of the other species.

∞ 0 e

 0 -βt w(x(t | x 0 )) dt. Using (38) again on v(x T ) we get -a ∞ 0 e -βt w(x(t | x T )) dt ≥ -∞ 0 e -βt w(x(t | x 0 )) dt, or a ∞ 0 e -βt w(x(t | x T )) dt ≤ ∞ 0 e -βt w(x(t | x 0 )) dt.

  v(x 0 ) = e -βτ v(x T ) v(x(τ )) -τ 0 w(x(t | x 0 )) dt. Using again av(x T ) ≥ v(x 0 ) if w ≥ 0 we get av(x T ) ≥ e -βτ v(x T ) -τ 0 w(x(t | x 0 )) dt, or (e -βτ -a) ∞ 0 e -βt w(x(t | x T ))dt ≥ -τ 0 w(x(t | x 0 )) dt.

0 e 0 e 0 e= ae βτ τ 0 e 0 e

 00000 ) gives (e -βτ -a)∞ 0 e -βt w(x(t | x T ))dt ≥ 0 and therefore a ≤ e -βτ since the integral is strictly positive. This proves the first two claims.To finish we observe that any solution to (37) satisfiese -βτ v(x T ) = v(x 0 ) + τ 0 e -βt w(x(t | x 0 )) dt. Therefore av(x T ) = v(x 0 )ae βτ + ae βτ τ -βt w(x(t | x 0 )) dt.Using (38) we getav(x T ) = v(x 0 ) + ae βτ ≥0 τ -βt w(x(t | x 0 )) dt + (1 -ae βτ ) ≥0 ∞ -βt w(x(t | x 0 )) dt. Since av(x T ) -v(x 0 ) = -βt w(x(t | x 0 )) dt + (1 -ae βτ )∞ -βt w(x(t | x 0 )) dt,

∞ 0 I 0 e

 00 A (t)e βt dτ x 0 (t) + 1 -∞ βt dτ x 0 (t) I A (∞), A ⊂ [0, ∞].

Such modification is always possible. For instance let f (x) = min y∈X {f (y) + m i=1 f ui (y)u i (y) + L xy }, where L is the Lipschitz constant of f + m i=1 f ui u i on X.

 Theorem 3The optimal value of (9) of p satisfies

Proof: Follows from Corollary 2 since ρ 0 ≥ 1 and V u ≥ V ≥ 0. Now we are in a position to prove the following crucial lemma linking problems (4) and [START_REF] Korda | Controller design and region of attraction estimation for nonlinear dynamical systems[END_REF].

Lemma 2 If {u k ∈ C(X; U )} ∞ k=1 and {τ k ∈ L(X; [0, ∞])} ∞ k=1 are respectively sequences of controllers and stopping functions satisfying the conditions of Assumption 1, then the corresponding densities {ρ k , ρ k 0 , ρ k T , σ k } ∞ k=1 with ρ k 0 = 1 are feasible in [START_REF] Korda | Controller design and region of attraction estimation for nonlinear dynamical systems[END_REF] and satisfy

Conversely, if {ρ k , ρ k 0 , ρ k T , σ k } ∞ k=1 is a sequence such that lim k→∞ p(ρ k , ρ k 0 , ρ k T , σ k ) = p and if Assumption (1) holds, then equation ( 7) holds with u k = σ k /ρ k .

Proof: To prove the first part of the statement consider the controllers u k , stopping functions τ k and densities ρ k from Assumption [START_REF] Ambrosio | Existence, uniqueness, stability and differentiability properties of the flow associated to weakly differentiable vector fields[END_REF]. Setting ρ k 0 = 1 and defining [START_REF] Korda | Controller design and region of attraction estimation for nonlinear dynamical systems[END_REF]. In addition, in view of ( 8), the representation (17) holds with τ x 0 = δ τ k (x 0 ) . Therefore by Lemma 1

and hence [START_REF] Prajna | Nonlinear control synthesis by convex optimization[END_REF] holds since {V u k ,τ k } ∞ k=1 satisfies [START_REF] Korda | Convex computation of the maximum controlled invariant set for polynomial control systems[END_REF] and ρ k 0 = 1 for all k ≥ 0. To prove the second part of the statement, let {ρ k , ρ k 0 , ρ k T , σ k } ∞ k=1 be any sequence such that lim k→∞ p(ρ k , ρ k 0 , ρ k T , σ k ) = p . Then this sequence satisfies [START_REF] Prajna | Nonlinear control synthesis by convex optimization[END_REF] by Theorem 3 and by the first part of Lemma 2 just proven. Therefore [START_REF] Korda | Convex computation of the maximum controlled invariant set for polynomial control systems[END_REF] holds with

We will also need the following result showing that nonnegative C 1 functions vanishing on a neighborhood of ∂X can be approximated by polynomials in Q d (K) vanishing on ∂X.

Then for any > 0 there exists d ≥ 0 and a polynomial

and p d = 0 on ∂X.

Note that, trivially, approximations from above to V u provide approximations from above to

, the degree d polynomial upper and lower bounds are given by min

and max

respectively. Fixing a basis of R[x] d , the objective functions of ( 28) and (29) become linear in the coefficients of V u respectively V u in this basis. Problems (28) and ( 29) are therefore convex SOS problems and immediately translate to SDPs (see Section 2.2).

Theorem 4 Let V u d and V u d denote solutions to (28) and (29) of degree d.

Proof: See Appendix A.

As a simple corollary we obtain a converging sequence of polynomial over-approximations to V , the optimal value function of (4):

u d 1 denote the degree d 2 polynomial approximation from above to the value function associated to the rational controller u d 1 obtained from (10) using [START_REF] Lasserre | Global Optimization with Polynomials and the Problem of Moments[END_REF]. Then V

Now we describe a hierarchy of lower bounds on V : max

Proof: Follows by similar arguments based on Gronwall's Lemma as in the proof of Theorem 4. and hence

is the first exit time of X. Next we observe that V u , the value function associated to u, is equal to

In view of (34), we conclude that

is polynomial and hence bounded on X (and hence e -βt V u d (x(t | x 0 )) → 0); and we conclude that

Convergence of the upper and lower bounds (30) follows from Theorem 7 using infinitedimensional LP duality and standard results on the convergence of moment relaxations.

The proof is similar to the proof of Theorem 5 in [START_REF] Korda | Convex computation of the maximum controlled invariant set for polynomial control systems[END_REF] or Theorem 3.6 in [START_REF] Lasserre | Nonlinear optimal control via occupation measures and LMI relaxations[END_REF] and therefore we only outline it. The hierarchy of SOS programming problems ( 28) and ( 29) is dual to the hierarchy of moment relaxations of an infinite-dimensional LP in the cone of nonnegative measures whose dual is an infinite-dimensional LP in C 1 (X) and feasible solutions of this dual provide upper or lower bounds on V u . Crucial to applying infinite-dimensional duality results (e.g., [2, Theorem 3.10]) is the boundedness of measures satisfying the discounted Liouville equation ( 14) with ν i ≤ ūµ and µ 0 = λ X , where λ X is the restriction of the Lebesgue measure to X. Plugging v = 1 in [START_REF] Löfberg | YALMIP : A toolbox for modeling and optimization in MATLAB[END_REF] we have µ T (X) + βµ(X) = µ 0 (X). Since µ 0 (X) = λ X (X) = vol X < ∞ and β > 0 we conclude that µ T and µ are indeed bounded, which implies that ν i is also bounded for i = 1, . . . , m. Equally important is the absence of duality gap between the finite-dimensional moment relaxations and SOS tightenings (which are both SDP problems); this follows immediately from the presence of the constraint g i = N -x 2 2 among the constraints describing X, which implies the boundedness of the truncated moment sequences feasible in the moment relaxations. The absence of duality gap then follows from [START_REF] Trnovská | Strong duality conditions in semidefinite programming[END_REF]Lemma 2].

Appendix B

This appendix presents a proof of Theorem 2. We will prove a slightly more general version of the result from which Theorem 2 immediately follows:

Theorem 7 Let f : R n → R n be globally Lipschitz and let the nonnegative measures µ, µ 0 , µ T on R n satisfy

for all v ∈ C 1 (R n ). Then there exists an ensemble of probability measures {τ x 0 } x 0 ∈X with spt τ x 0 ⊂ [0, ∞] and an ensemble of trajectories {x(• | x 0 )} x 0 ∈X of the ODE ẋ = f (x) and

Theorem 2 follows from Theorem 7 by setting f = f + m i=1 f u i u i and modifying f and f u i outside the compact set X such that f is globally Lipschitz 4 . The conclusion that x(t | x 0 ) ∈ X for all t ∈ spt τ x 0 follows by taking v(x) = e -x 2 I R n \X (x) in (36), where I A is the indicator function of a set A, i.e. I A (x) = 1 if x ∈ A and I A (x) = 0 otherwise. Suppose therefore that (36) holds. First we will prove a simple result. In the rest of this Appendix we will use the notation C k c for the space of all compactly supported k-times continuously differentiable functions. Proof: Since f is globally Lipschitz the solution x(t | x 0 ) is defined for all x 0 ∈ R n and all t ≥ 0. Therefore (38) is well defined (notice that w is bounded and β > 0). Direct