# Optimal Control of State Constrained Systems via Measure Relaxations and Polynomial Optimization

Saroj Prasad Chhatoi, Aneel Tanwani, Didier Henrion

► To cite this version:

Saroj Prasad Chhatoi, Aneel Tanwani, Didier Henrion. Optimal Control of State Constrained Systems via Measure Relaxations and Polynomial Optimization. 2025. hal-05345286

HAL Id: hal-05345286

https://hal.science/hal-05345286v1

Preprint submitted on 8 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimal Control of State Constrained Systems via Measure Relaxations and Polynomial Optimization

Saroj Prasad Chhatoi Aneel Tanwani Didier Henrion

***Abstract***—We address the optimal control problem for a class of dynamical systems with constrained state trajectories. These systems are modeled by a differential inclusion with a drift term and a normal cone mapping associated with the constraint set. The optimal control problem is considered in continuous-time and discrete-time, where the latter provides a computational advantage over the former. In both cases, the nonlinear problem is reformulated as an infinite-dimensional linear program over occupation measures. We show that this does not introduce any relaxation gap, that is, the optimal value remains the same for the reformulated linear program. Using appropriate tools from functional analysis and optimal transport, we also show the convergence of the optimal value of the discrete problem to the optimal value of the continuous problem. We propose finite-dimensional convex optimization algorithms based on the moment-sum-of-squares hierarchy to provide numerical approximations of the proposed infinite-dimensional linear programs.

***Index Terms***—Nonsmooth systems; Optimal control; Measure relaxations; Young measure; Occupation measure; Moment sum-of-squares approximation.

## I. INTRODUCTION

Optimal control of nonlinear systems has remained a problem of interest for the control community over several decades. Developments on the theoretical side as well as numerical methods for efficiently solving such problems continue to attract attention of researchers from different community. In this article, we consider the optimal control problem for a class of nonsmooth systems, and propose numerically tractable algorithms for our approach. In particular, the systems that we consider are described by a differential inclusion where the right-hand side comprises a drift term and a normal cone mapping associated to a time-varying constraint set; see [1] for an overview of such models and their relevance in various applications that include robotics, mechanics, electrical circuits, market equilibrium, and power systems with switching.

Let us recall that the optimal control of nonlinear systems with well-posedness properties is a well-studied problem [2], [3]. In contrast, optimal control for nonsmooth dynamical systems with set-valued right-hand side is more complicated and has been an area of active research over the last two decades because of its theoretical and practical relevance [4].

Among the existing work on optimal control problems for the class of nonsmooth systems (which are closely related to the ones studied here), we note that the authors in [5]

consider an evolution variational inequality with a closed convex constraint set. By deriving the first-order optimality conditions for the regularized problem, and taking the limit as the regularization parameter tends to zero, the corresponding optimality conditions for the original nonsmooth system are obtained. Pontryagin Maximum Principle based optimality conditions for discrete approximations of such systems are derived in [6] and sophisticated tools from variational analysis are used to derive the limit equation. In [7], the authors use exact penalization-based techniques to extend the optimality conditions to more general problems with weaker assumptions on the system data. In a more recent article, [8] uses distance to the constraint set as the penalization which allows the authors to obtain uniform, w.r.t. control estimates and they show strong convergence of the approximating trajectories in $W^{1,2}$ metric.

Direct approach based on discretizing the problem and solving the resulting static optimization have several fundamental issues. In [9], the authors pointed out that the error in gradients of the simulation results are independent of the step sizes, which means that optimization problems based on such discretization schemes often fail. The authors propose an optimal control strategy that approximates the discontinuous system with a smooth right-hand side and use the Euler integrator, ensuring that the gradients approach the true values under sufficiently small step sizes relative to the smoothing parameter. Finite-element-based discretization along with time-freezing methods have been recently proposed in software packages [10]. Further in [11], the authors extend the finite element with switch detection method to address optimal control problems with nonsmooth differential equations by transforming them into dynamic complementarity systems. For a class of linear complementarity systems (closely related to the formalism adopted in our present work), but where the right-hand side of the differential equation is a Lipschitz continuous map, the paper [12] presents first-order optimality conditions and some numerical results.

On the other hand, global methods for nonlinear control based on occupation measures have gained popularity in the last decade because of their powerful modeling capabilities and the availability of efficient algorithms and semidefinite programming (SDP) solvers [13]. These methods reformulate a finite-dimensional nonlinear optimal control problem into a primal/dual pair of infinite-dimensional linear problems. The primal problem is expressed in the cone of non-negative Borel measures and the dual problem is expressed in the cone of non-negative continuous functions [14], [15]. The infinite-dimensional linear problems admit a sequence of finite dimensional relaxations using the moment-sums-of-

<footnote>S. P. Chhatoi, A. Tanwani, and D. Henrion are with LAAS-CNRS, University of Toulouse, France. D. Henrion is also with the Faculty of Electrical Engineering of the Czech Technical University in Prague, Czechia. spchhatoi@laas.fr, tanwani@laas.fr, henrion@laas.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1: Outline of the approach developed in this paper for the optimal control of a class of nonsmooth dynamical systems. For the continuous-time case, we first relax the problem using Young measures and then the occupation measures. It is shown that this relaxation does not introduce any relaxation gap. For the discrete-time version, the problem is lifted to the cone of non-negative Borel measures. We further show the convergence in the limit as $\tau \to 0$ (no gap).</figcaption></figure>

square(SOS) hierarchy [13]. Each relaxation is formulated as an SDP and can be solved using off-the-shelf SDP solvers. Under mild assumptions, the sequence yields a monotonically non-decreasing sequence of lower bounds converging to the true optima. Occupation measures were also studied in the context of control of stochastic systems [16], [17]. One more advantage of occupation measure-based methods is the ability to easily handle uncertainty in the initial distribution.

Some tools in our approach come from the theory of optimal transport – a research topic dealing with the optimal way to transport an initial distribution of mass to some final distribution. In [18], the authors solve an optimal transport problem for a Lagrangian-based cost obtained from a dynamical system with nonholonomic constraints. The problem of uniqueness and the existence of solutions for Lagrangian cost derived from linear quadratic regulator were studied in [19]. Similar ideas have been pursued in [20]. In these works, the problem of mass transport for Lagrangian systems is seen as a decoupled problem where the lower level computes the optimal control for the finite-dimensional particle system and then this cost is used to find an optimal map to transport the initial measures to the final configuration.

In the present work, we address the problem of optimal control of nonsmooth dynamical systems. We use occupation measures to relax both the discrete-time problem and the continuous-time problem into an infinite-dimensional linear problem defined in the cone of non-negative Borel measures. Under some mild assumptions, we show that this does not produce any relaxation gap. This means that the value of the linear problem is equal to the value of the original optimal control problem. Further, we use tools from optimal transport

theory to show the convergence of the discrete-time problem to the continuous-time problem when the time step size goes to zero. The overall layout is summarized in Figure 1. Apart from the independent theoretical interest, the convergence of the discrete-time problem to the continuous-time problem is also of practical importance, as solving the discrete-time problem provides an approximation to the continuous-time problem. We also propose numerical techniques which provide finite dimensional approximation of the relaxed discrete-time problem.

### A. Problem Setting

We consider a controlled nonsmooth dynamical system modeled as an evolution variational inequality:

$$ \dot{x}(t) \in f(x(t), u(t)) - N_{S(t)}(x(t)), \quad x(0) = x_0 \quad (1) $$

where $f: \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}^n$ is a given map (called vector field, or drift), and $N_{S(t)}(x(t))$ denotes the outward normal cone to the compact convex set $S(t)$ at the point $x(t) \in \mathbb{R}^n$ at each time $t \in [0, T]$. The control $u(t) \in \mathbb{R}^m$ at each time $t$ belongs to a convex compact set $U \subset \mathbb{R}^m$. We make the following assumption for the existence and uniqueness of the solutions.

**Assumption 1.** There exist a compact set $\Omega$ and constant $L_f > 0$ such that, $S(t) \subset \Omega$ for each $t \in [0, T]$, and for every $x, x_1, x_2 \in \Omega$, we have

$$ \begin{gathered} \max_{u \in U} |f(x, u)| \le L_f(1 + |x|) \\ |f(x_1, u) - f(x_2, u)| \le L_f |x_1 - x_2| \quad \forall u \in U. \end{gathered} $$

**Assumption 2.** The set-valued mapping $S: [0, T] \to \mathbb{R}^n$ satisfies the following:
• For every $t \in [0, T]$, $S(t)$ is nonempty, closed and convex.

• $S(\cdot)$ is Lipschitz continuous w.r.t. the Hausdorff distance: there exists $L_s \ge 0$ such that $d_H(S(t), S(s)) \le L_s|t-s|$ for $t, s \in [0, T]$ where $d_H(A, B) = \max\{\sup_x \text{dist}(x, A), \sup_x \text{dist}(x, B)\}$ is the Hausdorff distance between sets $A, B \subset \mathbb{R}^n$.

An optimal control problem for such a system, over a finite interval $[0, T]$, can be formulated as

$$ J^*(x_0) = \inf_{u \in L^\infty([0,T]; \mathbb{R}^m)} \int_0^T \mathcal{L}(x(t), u(t))dt \quad (2) \\ \text{s.t.} \quad \dot{x}(t) \in -\mathcal{N}_{S(t)}(x(t)) + f(x(t), u(t)), \\ x(0) = x_0, \quad u(t) \in U, \quad \forall t \in [0, T] $$

where $\mathcal{L}(x, u)$ is a continuous cost function in $x, u$. Further assumptions on the cost function and the drift term, that are necessary for our purposes, are formulated as follows:

**Assumption 3.** For each $x \in \mathbb{R}^n$, the image set $f(x, U)$ is convex. The cost $\mathcal{L}(x, u)$ is convex in $u$ for each $x$, and satisfies the growth bound:

$$ \mathcal{L}(x, u) \le l(x)(1 + |u|^2) \quad \forall x \in \mathbb{R}^n \text{ and } u \in \mathbb{R}^m \quad (3) $$

where $l$ is locally bounded and Borel-measurable.

In order to study the discrete-time version of (2), consider a partition of $[0, T] = \{0 = t_0, t_1, \dots, t_N = T\}$ with $t_k - t_{k-1} = \tau$ being the time step between two samples. We use a time-stepping-based algorithm to define the evolution of states [1]. Let $x_k^\tau$ resp. $x_{k+1}^\tau$ be the state at time $k\tau$ resp. $(k+1)\tau$, and let $u_k$ be the control at time $k\tau$. They are related by

$$ x_{k+1}^\tau = P_{S_k}(\tau f_k(x_k^\tau, u_k) + x_k^\tau) \quad (4) $$

where $P_{S_k}$ is the (unique) projection onto the compact convex set $S_k := S(k\tau)$.

The discrete-time optimal control problem can then be written as

$$ J_{\tau}^{*}(x_{0}) = \inf_{x_{k}^{\tau} \in \mathbb{R}^{n}} \sum_{k=0}^{N} \tau \mathcal{L}(x_{k}^{\tau}, u_{k}) \quad (5) \\ \text{s.t.} \quad x_{k+1}^{\tau} = P_{S_{k}}(\tau f_{k}(x_{k}^{\tau}, u_{k}) + x_{k}) \\ x(0) = x_{0}, \quad u_{k} \in U, \quad \forall k = 0, 1, \dots, N. $$

## B. Contribution

The primary problem studied in this paper concerns the continuous-time problem (2) and discrete-time problem (5) for the nonsmooth system (1). In particular, we provide an alternative formulation for these problems and study the relation between them as the sampling time tends to zero. The main contributions can be summarized as follows:

• We study the optimal control problem in both continuous-time (2) and discrete-time (5) by formulating or relaxing them to linear optimization problems in occupation measures.

• We prove that this relaxation introduces no gap, meaning the value of the linear problem on measures matches the value of the original optimal control problem.

• We analyze the convergence of the relaxed discrete-time problem to the continuous-time problem using interpolation schemes based on optimal transport.

• We propose numerical techniques for solving approximately the infinite-dimensional linear programs using the moment-SOS hierarchy, a family of finite-dimensional semidefinite programs of increasing size.

A preliminary version of this work was presented in the conference paper [21]. In this article, the theoretical framework is developed in a more rigorous manner with complete proofs and the results have been further generalized to accommodate time-varying constraints. Additional contributions include the analysis of the convergence of the value function, as well as new sections on numerical methods and an illustrative example. Due to space constraints, proofs of certain technical statements have been included in the extended version [22].

## C. Outline

The remainder of the paper is organized as follows: In Section II, we discuss the relaxation of the continuous-time optimal control problem, first using Young measures and then using occupation measures. We also examine key aspects of this relaxation process, such as the absence of a relaxation gap. In Section III, we formulate the relaxation of the discrete-time problem in the space of measures, ensuring that the relaxation does not introduce any gap between the optimal values of the problems. In Section IV, we analyze the convergence of the relaxed discrete-time optimal control problem to its continuous-time counterpart. Finally, in Section V, we introduce relaxations based on semidefinite programming for the relaxed discrete-time problem and illustrate the proposed approach with an academic example.

## D. Notation

Here we introduce the key mathematical symbols and definitions used throughout the paper. The Borel $\sigma$-algebra $\mathcal{B}(X)$ of a topological space $X$ is the smallest collection of sets that can be formed by countable unions, intersections, and complements of open sets in $X$. We denote by $\mathcal{M}(X)$ the space of signed Borel measures on $X$, which are mappings from $\mathcal{B}(X)$ to $\mathbb{R}$, satisfying $\mu(\emptyset) = 0$ and $\mu(\bigcup_{i=1}^\infty A_i) = \sum_{i=1}^\infty \mu(A_i)$ for every countable pairwise disjoint collection $\{A_i\}$. The space of continuous functions from $X$ to $\mathbb{R}$, resp. $\mathbb{R}^n$, is denoted by $\mathcal{C}(X)$, resp. $\mathcal{C}(X; \mathbb{R}^n)$. If $X$ is compact, the dual space to $\mathcal{C}(X)$ consists of all bounded linear functionals on $\mathcal{C}(X)$ and it can be identified with the space of signed Borel measures $\mathcal{M}(X)$. We use $\mathcal{C}_+(X)$ and $\mathcal{M}_+(X)$ to denote the sets of nonnegative continuous functions and nonnegative Borel measures on $X$, respectively. We further denote by $\mathcal{C}^k(X)$ the space of $k$-times continuously differentiable functions on $X$, by $\mathcal{C}_0(X)$ the space of continuous functions that vanish at infinity, and by $\mathcal{C}_b(X)$ the space of continuous bounded functions.

The duality pairing between a measure $\mu \in \mathcal{M}(X)$ and a continuous function $f \in \mathcal{C}(X)$ corresponds to integration: $\langle \mu, f \rangle = \int_X f d\mu$.
For $1 \le p < \infty$, the space of $p$-integrable functions with respect to a measure $\mu$ on $X$ is denoted by

$$L^p(X, \mu) = \left\{ f : X \to \mathbb{R} \mid \int_X |f|^p d\mu < \infty \right\}.$$

When $\mu$ is the Lebesgue measure, we simply write $L^p(X)$. For $p = \infty$, $L^\infty(X, \mu)$ denotes the space of essentially bounded functions.

We say that a sequence $\{f_j\}_{j \in \mathbb{N}}$ of functions in $L^p(X)$ converges weakly to a function $f \in L^p(X)$, denoted $f_j \rightharpoonup f$, whenever

$$\lim_{j \to \infty} \int_X f_j g \, d\mu = \int_X fg \, d\mu \quad \text{for all } g \in L^q(X),$$

where $1 \le p < \infty$, $\frac{1}{p} + \frac{1}{q} = 1$ and $L^q(X)$ is the dual space to $L^p(X)$. In dual spaces, weak-$\star$ convergence is used, which tests convergence against functions in the pre-dual space. For example, we say that a sequence $\{f_j\}$ of functions in $L^\infty(X)$ converges weak-$\star$ to a function $f$, denoted $f_j \stackrel{\star}{\rightharpoonup} f$, whenever

$$\lim_{j \to \infty} \int_X f_j g = \int_X fg \quad \text{for all } g \in L^1(X),$$

where we note that $L^1(X)$ is the pre-dual of $L^\infty(X)$. Similarly, we say that a sequence $\{\mu_j\}$ of measures in $M(X)$ converges weak-$\star$ to a measure $\mu \in M(X)$, denoted $\mu_j \stackrel{\star}{\rightharpoonup} \mu$, whenever $\lim_{j \to \infty} \int_X f d\mu_j = \int_X f d\mu$ for all $f \in C(X)$. Given a measurable function $f: X \to Y$, the push-forward of a measure $\mu$ is denoted $f_\#\mu$ and is defined by $(f_\#\mu)(B) = \mu(f^{-1}(B))$ for all $B \in \mathcal{B}(Y)$. For comprehensive treatment of related notions from functional analysis and measure theory, see [23], [24].

## II. RELAXED CONTINUOUS-TIME OPTIMAL CONTROL

In this section, our aim is to describe the relaxation of the continuous time optimal control problem (2) in the space of non-negative Borel measures. As a first step in this direction, in Theorem 1, we study weak convergence of the sequence of inputs $(u_j)$ using Young measures, and show that the corresponding sequence of state trajectories $(x_j)$ is strongly convergent. We then introduce the relaxation of problem (2) using Young measures and demonstrate in Theorem 2 that this does not alter the optimal value, i.e. there is no relaxation gap. Then, we introduce occupation measures, which allow us to further relax the optimal control problem so that the resulting problem becomes linear with respect to the measures. Finally, in Theorem 3, we prove that this last reformulation does not introduce any relaxation gap.

### A. Minimizing sequences and Young measures

Let us consider a minimizing sequence $(x_j(\cdot), u_j(\cdot))_{j \in \mathbb{N}}$ of control and state trajectories. By that, we mean that each state trajectory $x_j(\cdot)$ is a solution to (1) for the given control $u_j(\cdot)$, and the value of the problem (2) decreases when $j$ increases. Next we consider the control sequence of this minimizing sequence to be weakly convergent and we study the convergence properties of the sequence of associated state trajectories.

Note however that the nonlinearity of the dynamics in (1) does not preserve weak convergence. For a simple illustrating example, consider a sequence of functions $(g_j)_{j \in \mathbb{N}}$ defined as follows:

$$g_j(t) = \sum_{i=0}^{j-1} \mathbf{1}_{[\frac{2i}{2j}, \frac{2i+1}{2j}]}(t) - \mathbf{1}_{[\frac{2i+1}{2j}, \frac{2i+2}{2j}]}(t) \quad (6)$$

where $\mathbf{1}_I(t)$ is equal to one if $t \in I$ and zero otherwise. Given $p \in [1, \infty)$, the sequence $(g_j)_{j \in \mathbb{N}}$ converges weakly to zero in $L^p([0, 1])$, i.e. $\lim_{j \to \infty} \int \phi(t)g_j(t)dt = 0$ for all $\phi \in L^q([0, 1])$, where $q \in [1, \infty)$ satisfies $\frac{1}{p} + \frac{1}{q} = 1$. But if we consider a nonlinearity $f(u) = |u|$, then $f(g_j(t)) = 1$ for all $j \in \mathbb{N}$. So, we observe that $\lim_{j \to \infty} f(g_j) \neq f(\lim_{j \to \infty} g_j) = f(0) = 0$. Rather, the oscillating behavior of the weakly converging sequence can be captured by Young measures [25].

**Definition 1 (Young measure).** Let $I = [0, T]$ and $X = \mathbb{R}^n$. A Young measure on $I$ with values in $X$ is a mapping $t \mapsto \nu_t \in \mathcal{P}(X)$ such that for every $\phi \in C_b(X)$ the function

$$t \mapsto \int_X \phi(x) \, d\nu_t(x)$$

is Borel measurable on $I$.

In the previous example, the sequence generates a family of homogeneous (i.e. independent of $t$) Young measures $\omega_t = \frac{\delta_{-1} + \delta_1}{2}$ which capture correctly the weak-$\star$ convergence

$$\lim_{j \to \infty} \int_0^T \phi(t) f(g_j(t)) dt = \int_0^T \int_{\mathbb{R}} \phi(t) f(u) d\omega_t(u) dt \\ = \int_0^T \int_{\mathbb{R}} \phi(t) |u| \frac{\delta_{-1} + \delta_1}{2} dt = \int_0^T \phi(t) dt$$

where $\phi \in L^1([0, T])$. Thus, we obtain the correct weak-$\star$ limit of the compositions, $f(g_j(\cdot)) \stackrel{\star}{\rightharpoonup} 1$. The next result formalizes this mode of convergence: Young measures encode the limits of compositions $\psi(t, x_j(t), u_j(t))$ when $(x_j, u_j)$ converge weakly and $\psi$ is a Carathéodory integrand. A function $\psi: [0, T] \times \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}$ is called *Carathéodory* if it is measurable in $t$ for each fixed $(x, u)$, and continuous in $(x, u)$ for almost every $t \in [0, T]$.

**Proposition 1** ([25, Theorem 6.2]). Let $(x_j, u_j)_{j \in \mathbb{N}}$ be a sequence of measurable and equibounded functions, and let $\psi: [0, T] \times \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}$ be a Carathéodory function. For every weakly convergent sequence $(\psi(t, x_j(t), u_j(t)))_{j \in \mathbb{N}}$ in $L^1([0, T])$, there exists a subsequence (without relabelling), and a family of probability measures $(\nu_t)_{t \in [0, T]}$, such that

$$\lim_{j \to \infty} \int_0^T \psi(t, x_j(t), u_j(t)) dt = \int_0^T \int_{\mathbb{R}^n \times \mathbb{R}^m} \psi(t, x, u) d\nu_t(x, u) dt.$$

The above result can be understood as a statement about the convergence properties of sequences of compositions of non-linear functions, with Young measures encoding the limiting behaviour. In this context:
• The sequence $(x_j, u_j)_{j \in \mathbb{N}}$ can be seen as generating a sequence of Dirac measures $(\delta_{(x_j, u_j)})_{j \in \mathbb{N}}$ in the space $L^\infty([0, T]; \mathcal{M}(\mathbb{R}^n \times \mathbb{R}^m))$. <fnref n="1" />

• The weak-⋆ convergence of $\delta_{(x_j, u_j)}$ in this measure space yields a Young measure $(\omega_t)_{t \in [0, T]}$ which captures the limiting behavior of $(x_j, u_j)$ when tested against functions in $L^1([0, T]; C_0(\mathbb{R}^n \times \mathbb{R}^m))$.

Young measures also help to establish strong convergence of state trajectories for certain class of dynamical systems even under weak convergence of control sequences. For example, consider the following dynamical system:

$$ \dot{x}(t) = |u(t)| \quad (7) $$

with the sequence of inputs $(u_j)_{j \in \mathbb{N}} = (g_j)_{j \in \mathbb{N}}$ where $g_j$ is defined in (6) and $x(0) = 0$. The sequence converges weakly to zero in $L^p([0, T])$ for each $p \in [1, \infty)$, and in weak-⋆ sense for $p = \infty$. Thus, we might expect that, in the limit, $\dot{x}(t) = |u(t)| = 0$ and hence $x(t) = 0$ for $t \in [0, T]$. But the sequence of solutions to (7), corresponding to $u_j$, is given by $x_j(t) = t$ which uniformly converges to $x(t) = t$. This sort of convergence can be explained by the use of Young measures, where we describe the limiting control by a measure and system (7) with such inputs is interpreted as follows:

$$ \dot{x}(t) = \int |u| d\omega_t(u) = \int |u| \frac{\delta_{-1} + \delta_1}{2} = 1. \quad (8) $$

The last equation clearly implies that $x(t) = t$. Note that system (8) is a relaxation of system (7) in the sense that system (8) corresponds to the particular choice of Young measure $d\omega_t(u) = \delta_{u(t)}$ in system (8). We therefore expect system (8) to have a larger set of solutions than system (7).

We next show that, given a weakly convergent control sequence $(u_j(\cdot))_{j \in \mathbb{N}}$, the corresponding sequence of state trajectories $(x_j(\cdot))_{j \in \mathbb{N}}$ is strongly convergent when interpreting limiting controls in the sense of Young measure. Moreover, each pair $(x_j(\cdot), u_j(\cdot))$ satisfies (1) and the limit trajectory $x(\cdot)$ satisfies

$$ \dot{x}(t) \in -\mathcal{N}_{S(t)}(x(t)) + \int_U f(x(t), u) d\omega(u|t). \quad (9) $$

**Theorem 1.** Let $(u_j)_{j \in \mathbb{N}}$ be a weakly converging sequence, in $L^\infty([0, T]; \mathbb{R}^m)$, which generates the limiting Young measure $d\omega_t$. Then, up to a subsequence, the sequence $(x_j)_{j \in \mathbb{N}}$ of associated solution trajectories (1) converges uniformly to $x$ satisfying (9).

*Proof.* Let us consider the weakly converging sequence $(u_j)_{j \in \mathbb{N}}$ with the Young measure $\omega_t$ as the limit (in the sense of Proposition 1) and associate with it a sequence of absolutely continuous functions $(x_j)_{j \in \mathbb{N}}$ where each $x_j: [0, T] \to \mathbb{R}^n$ is obtained as a solution to (1) by applying the control input $u_j$, $j \in \mathbb{N}$. We aim to prove that the sequence $(x_j)_{j \in \mathbb{N}}$ is uniformly

$$ \langle \nu, \phi \rangle = \int_0^T \int_{\mathbb{R}^n} \phi(t, x) \nu_t(dx) dt, \\ \text{for } \nu \in L^\infty([0, T]; \mathcal{M}(\mathbb{R}^n)) \text{ and } \phi \in L^1([0, T]; C_0(\mathbb{R}^n)). $$

<figure><img src="image_3.png" /><figcaption>Fig. 2: Schematic diagram of the proof of Theorem 1.</figcaption></figure>

convergent to $x$, where $x$ is described by (9). The proof consists of three steps: (1) approximation by regularization, (2) convergence to Young measure, (3) uniform convergence of the relaxed dynamics.

1) Approximation by regularization: We consider the Moreau-Yosida regularization of the dynamics in (1)

$$ \dot{x}_j^\lambda(t) = f(x_j^\lambda(t), u_j(t)) - \frac{1}{\lambda}(x_j^\lambda(t) - P_{S(t)}(x_j^\lambda(t))). \quad (10) $$

For $\lambda$ small enough, $x_j^\lambda(t) \in \Omega^2$ and we borrow the following inequalities from [26, Theorem 3.1]:

$$ |x_j^\lambda(t)| \le \exp(2L_fT)|x_j^\lambda(0)| + (\exp(2L_fT) - 1) \frac{2L_f + L_s}{2L_f}, \quad (11) $$

$$ |\dot{x}_j^\lambda(t)| \le 2L_f + L_f |x_j^\lambda(t)| + L_f \max_{0 \le s \le t} |x_j^\lambda(s)| + L_s \quad (12) $$

which establish the bounds on $x_j^\lambda(t)$ and $\dot{x}_j^\lambda(t)$ are uniform in $j$ and $\lambda$. So, considering a sequence of $x_j^\lambda(\cdot)$ for $\lambda \to 0$, we can use the Arzelà-Ascoli theorem to conclude that the sequence converges uniformly to $x_j(\cdot)$ satisfying (1) with $u = u_j(t)$.

2) Convergence to Young measure: Now we address the convergence of sequence $x_j^\lambda(\cdot) \to x^\lambda(\cdot)$ as $j \to \infty$. Using the bounds in (11), we conclude that the sequences $x_j^\lambda(\cdot), \dot{x}_j^\lambda(\cdot)$ are uniformly bounded in $j$. So, using the Arzelà-Ascoli theorem, we conclude that $x_j^\lambda$ converges uniformly (up to a subsequence) to some $x^\lambda$. To characterize the limit $x^\lambda(t)$, we use Proposition 1 to infer that the sequence $\{x_j^\lambda(\cdot), u_j(\cdot)\}$ generates a Young measure $\nu^\lambda(\cdot, \cdot|t) \in \mathcal{P}(\mathbb{R}^n \times U)$ such that

$$ \dot{x}^{\lambda}(t) = \int_{\mathbb{R}^n \times U} \left[ -\frac{1}{\lambda}(x - P_{S(t)}(x)) + f(x, u) \right] d\nu^{\lambda}(x, u|t). \quad (13) $$

<footnote n="1">¹The spaces $L^1([0, T]; C_0(\mathbb{R}^n))$ and $L^\infty([0, T]; \mathcal{M}(\mathbb{R}^n))$ form a dual pair, with their duality relation given by</footnote>

<footnote n="2">²In Assumption 1, the set $\Omega$ can be chosen such that $x_j^\lambda(t) \in \Omega$ holds for small enough $\lambda$ and all $t \in [0, T]$.</footnote>
In order to show that the trajectories $x_j^\lambda(t)$ converge pointwise to $x^\lambda(t)$ at each time $t$, we consider

$$|x^\lambda(t) - x_j^\lambda(t)| = |x_0 + \int_0^t \int_{\mathbb{R}^n \times U} [f(x, u) - \frac{1}{\lambda}(x - P_{S(t)}(x))] d\nu^\lambda(x, u|s) ds - x_0 - \left[ \int_0^t f(x_j^\lambda(t), u_j(t)) - \frac{1}{\lambda}(x_j^\lambda(t) - P_{S(t)}(x_j^\lambda(t))) \right]|$$

where we have used the fact that $x_j^\lambda(0) = x^\lambda(0) = x_0$. The pointwise convergence $\dot{x}_j^\lambda(t)$ to $\dot{x}^\lambda(t)$, implied by Proposition 1, and it being integrable allow us to use the dominated convergence theorem to conclude that $|x_j^\lambda(t) - x^\lambda(t)| \to 0$ pointwise in $t$ as $j \to \infty$. Moreover, since we have pointwise convergence of the whole sequence and uniform convergence up to a subsequence, we can conclude that each subsequence converges uniformly to the same limit $x^\lambda(t)$. The Young measure $\nu_t^\lambda$ can thus be written as $\nu^\lambda(dx, du|t) = \delta_{x^\lambda(t)}(dx)\omega(du|t)$ [25, Theorem 6.2]. Consequently, the dynamics in (13) reduces to

$$\dot{x}^{\lambda}(t) = -\frac{1}{\lambda}(x^{\lambda}(t) - P_{S(t)}(x^{\lambda}(t))) + \int_{U} f(x^{\lambda}(t), u) d\omega(u|t). \quad (14)$$

3) Uniform convergence of the relaxed dynamics: Next, we consider the convergence of the sequence $x^\lambda(t) \to x(t)$ as $\lambda \to 0$. We note that the term $\int_U f(x^\lambda(t), u) d\omega(u|t)$ in (14) is Lipschitz and it satisfies the linear growth condition. Thus, the dynamics in (14) satisfies bounds similar to (11) which ensure uniform convergence of $x^\lambda(t) \to x(t)$ using the Arzelà-Ascoli theorem. Using arguments similar to those in [26, Lemma 3.4], we conclude that the limit trajectory $x(\cdot)$ satisfies (9).

Thus we have the following inequality

$$\|x_j(\cdot) - x(\cdot)\| \le \|x_j(\cdot) - x_j^\lambda(\cdot)\| + \|x_j^\lambda(\cdot) - x^\lambda(\cdot)\| + \|x^\lambda(\cdot) - x(\cdot)\|$$

and we have already established uniform convergence for all the three terms on the right-hand side of the inequality. By applying the $\epsilon/3$ argument, we establish that $x_j(\cdot)$ converges uniformly to $x(\cdot)$, where $x(\cdot)$ is the solution of system (9). □

## B. Relaxation with Young measures

In the previous section, the formalism of Young measures was introduced to describe the selection of optimal control inputs. The dynamical system (1) also contains a multi-valued mapping for which we need to consider a selection rule. Following [27, Section 3], we make a selection from $\mathcal{N}_{S(t)}(x)$ by choosing a probability measure $\eta(\cdot|t, x(t)) \in \mathcal{P}(\mathcal{N}_{S(t)}(x(t)))$. Then, (1) can be expressed as

$$\dot{x}(t) = - \int_{\mathcal{N}_{S(t)}(x(t))} \zeta d\eta(\zeta|t, x(t)) + f(x(t), u(t)) \\ =: g(x(t), u(t)) \quad \text{s.t. } x(0) = x_0. \quad (15)$$

Let $(x_j(\cdot), u_j(\cdot))_{j \in \mathbb{N}}$ be a minimizing sequence of state and controls of the problem defined in (2) such that $\lim_{j \to \infty} \int \mathcal{L}(x_j(t), u_j(t)) dt \to J^*(x_0)$. Under Assumption 3, the function $\mathcal{L}(x, u)$ is upper bounded by $|u|^2$ and the uniform boundedness of $u_j(\cdot)$ implies that $(\mathcal{L}(x_j(\cdot), u_j(\cdot)))_{j \in \mathbb{N}}$ is uniformly integrable. Thus using Proposition 1, we get

$$\lim_{j \to \infty} \int_0^T \mathcal{L}(x_j(t), u_j(t)) dt = \int_0^T \int_{\mathbb{R}^n \times U} \mathcal{L}(x, u) d\nu_t(x, u) dt.$$

Following the results stated in Theorem 1 and [25, Theorem 6.2], $\nu_t$ can be written as $\nu_t = \delta_{x(t)}\omega_t$. Thus,

$$\lim_{j \to \infty} \int_0^T \mathcal{L}(x_j(t), u_j(t)) dt = \int_0^T \int_U \mathcal{L}(x(t), u) d\omega(u|t) dt.$$

Using this relaxation of the cost function and Theorem 1, the optimal control problem defined in (2) can be relaxed as follows

$$J_r^*(x_0) := \inf_{x(\cdot), \omega(\cdot|\cdot)} \int_0^T \int_U \mathcal{L}(x(t), u) \, d\omega(u|t) \, dt \quad (16)$$

subject to,

$$\dot{x}(t) = \int_U f(x(t), u) \, d\omega(u|t) - \int_{\mathcal{N}_{S(t)}(x(t))} \zeta \, d\eta(\zeta|t, x(t)), \quad (17)$$

$$x(0) = x_0, \\ \eta(\cdot|t, x(t)) \in \mathcal{P}(\mathcal{N}_{S(t)}(x(t))), \quad \omega(\cdot|t) \in \mathcal{P}(U), \quad \forall t \in [0, T].$$

For this relaxed optimal control problem, our next result shows that this relaxation does not introduce any gap in the optimal value.

**Theorem 2.** The optimal value $J^*(x_0)$ of optimal control problem (2) and the optimal value $J_r^*(x_0)$ of Young's measure relaxed optimal control problem (16) are equal.

*Proof.* (Proof of $J_r^*(x_0) \le J^*(x_0)$): We observe that if $(x(\cdot), u(\cdot))$ is feasible for (2), then it is feasible for (17). So if we consider an optimal solution $(x^*(\cdot), u^*(\cdot))$ of (2) then $(x^*(\cdot), \delta_{u^*(\cdot)})$ is feasible for (17), so

$$J^*(x_0) = \int_0^T \mathcal{L}(x^*(t), u^*(t)) dt \\ = \int_0^T \mathcal{L}(x^*(t), u) \delta_{u^*(\cdot)} dt \ge J_r^*(x_0)$$

where for the inequality we have used the fact that $(x^*(\cdot), \delta_{u^*(\cdot)})$ is a feasible solution to (16).

(Proof of $J^*(x_0) \le J_r^*(x_0)$): Since we know that $\mathcal{L}(\cdot, u)$ is a convex function in $u$ (based on Assumption 3), we can use Jensen's inequality to obtain

$$J_r^*(x_0) = \int_0^T \int_U \mathcal{L}(x(t), u) d\omega(u|t) dt \ge \int_0^T \mathcal{L}(x(t), u(t)) dt.$$

Now we invoke Assumption 3 which ensures that $g(x, U)$ is convex, so $\int g(x(t), u) d\omega(u|t) \in g(x(t), U)$. There exists a measurable selection $u(t) \in U$ such that $\dot{x}(t) =$
$g(x(t), u(t)) = \int g(x(t), u)d\omega(u|t)$ [28, Theorem 8.1.3]. The pair $(x(\cdot), u(\cdot))$ thus obtained is feasible for (2). So, we get

$$J_r^*(x_0) = \int_0^T \int_U \mathcal{L}(x(t), u)d\omega(u|t)dt \\ \geq \int_0^T \mathcal{L}(x(t), u(t))dt \geq J^*(x_0). \quad \square$$

The problem (16)–(17) above is still nonlinear in $x(\cdot)$, and in the next subsection, we use occupation measures which will lead us to a linear program in the space of non-negative Borel measures. To establish this, we make a passage from the Young measure to occupation measures by embedding the Young measure in the space of linear functionals.

### C. Occupation measures

Each triplet $(x(t), \eta(\cdot|t, x(t)), \omega(\cdot|t))$ for $t \in [0, T]$ which satisfy (17) can be associated with a pair of measures $(\hat{\mu}, \xi) \in (\mathcal{M}(\mathfrak{B} \times U), \mathcal{M}(S(T)))$ where $\mathcal{M}(\mathfrak{B} \times U)$ denotes the space of Borel measures on $\mathfrak{B} \times U$ for

$$\mathfrak{B} := \{(t, x, \zeta) \mid t \in [0, T], x \in \Omega \subset \mathbb{R}^n, \zeta \in \mathcal{N}_{S(t)}(x)\}, \quad (18)$$

and $\Omega \supset \bigcup_{t \in [0, T]} S(t)$. Let the pair $(\hat{\mu}, \xi)$ satisfy the following relation

$$\langle \hat{\mu}, h \rangle + \langle \xi, g \rangle = \int_0^T \int_U h(t, x(t), u)d\omega(u|t)dt + g(T, x(T)) \quad (19)$$

for all $h \in \mathcal{C}([0, T] \times \Omega \times U)$ and $g \in \mathcal{C}(S(T))$. Let us define

$$\mathfrak{R} := \left\{ \begin{array}{l} (\hat{\mu}, \xi) \in \mathcal{M}_+(\mathfrak{B} \times U) \times \mathcal{M}_+(S(T)) \quad (19) \text{ holds for} \\ \text{some } \omega(\cdot|t) \in \mathcal{P}(U) \text{ and } x(\cdot), \\ \text{with } x(\cdot) \text{ satisfying (17) for the same } \omega(\cdot|t) \\ \text{for some } \eta(\cdot|t, x(t)) \in \mathcal{P}(\mathcal{N}_S(x(t))) \quad \forall t \in [0, T]. \end{array} \right\} \quad (20)$$

The set $\mathfrak{R}$ provides an embedding of the triplet $(x(t), \eta(\cdot|t, x(t)), \omega(\cdot|t))$ for $t \in [0, T]$ into the space of measures $\mathcal{M}(\mathfrak{B} \times U) \times \mathcal{M}(S(T))$. A canonical embedding is given by

$$d\hat{\mu}(t, x, \zeta, u) = dt \delta_{x(t)}(dx) \nu(d\zeta, du | t, x(t)), \quad \xi := \delta_{x(T)},$$

with $\nu(\cdot|t, x(t)) \in \mathcal{P}(\mathcal{N}_{S(t)}(x(t)) \times U)$, whose conditional marginals satisfy $\pi_{\zeta\#} = \eta$ and $\pi_{u\#} = \omega$. Specifically, using (19), we derive the following bound on the norm of $\hat{\mu}$:

$$\|\hat{\mu}\| = \sup_{\|h\|_{\infty} \le 1} \left| \int_0^T \int_U h(t, x(t), u) d\omega(u|t) dt \right| \\ \le \int_0^T \|\omega\|_{TV} dt \le T$$

where $\|\omega\|_{TV}$ is the total variation norm of $\omega(\cdot|t) \in \mathcal{P}(U)$ (see [23] for the definition of total variation norm). Similarly, $\xi$ satisfies the following bound:

$$\|\xi\| = \sup_{\|g\|_{\infty} \le 1} |\langle \xi, g \rangle| = \sup_{\|g\|_{\infty} \le 1} |g(T, x(T))| \le 1$$

and thus any pair $(\hat{\mu}, \xi) \in \mathfrak{R}$ satisfies:

$$\|\hat{\mu}\| \le T; \quad \|\xi\| \le 1; \quad \hat{\mu}, \xi \ge 0. \quad (21)$$

Because the pair $(\hat{\mu}, \xi)$ satisfy (9) and (19), and both relations depend nonlinearly on the state trajectory $x(\cdot)$, convex combinations of admissible triplets need not be admissible; hence the set $\mathfrak{R}$ in (20) is not convex.

We observe that for every triplet $(x(t), \eta(\cdot|t, x(t)), \omega(\cdot|t))$ satisfying Young measure relaxed dynamics (17), the associated measure $d\mu(t, x, \zeta, u) = dt\delta_{x(t)}(u)d\omega(u|t)d\eta(\zeta|x(t))$ satisfies the following equation for every $\phi \in C^1([0, T] \times \Omega)$:

$$\int_{[0,T] \times S(t)} \int_{\mathcal{N}_{S(t)}(x)} \int_U \left[ \partial_t \phi(t, x) d\mu(t, x, \zeta, u) + \partial_x \phi(t, x) \cdot (-\zeta + f(x, u)) d\mu(t, x, \zeta, u) \right] = \int_{S(T)} \phi(T, x) d\mu_T(x) - \phi(0, x_0). \quad (22)$$

where $\mu \in \mathcal{M}(\mathfrak{B} \times U)$, $\mu_T \in \mathcal{M}(S(T))$. Equation (22) is called continuity equation or Liouville equation and the equation obtained is linear in $(\mu, \mu_T)$.

Next, we characterize all possible solutions to (22). Fix a constant $C_\eta > 0$ such that every admissible triplet $(x(t), \eta(\cdot|t, x(t)), \omega(\cdot|t))$ that satisfies the dynamics (17) also satisfies $\int_0^T \int_{\mathcal{N}_{S(t)}(x(t)) \times U} |\zeta| d\eta(\zeta, u|t, x(t)) dt \le C_\eta$ (the constant can be derived from the fact that $\dot{x}(\cdot)$ has uniform norm bound and the fact that we have $\int_U f(x(t), u)d\omega(u|t)$ is bounded uniformly from Assumption 1). Consider the set of pairs of measures $(\mu, \mu_T)$ such that

$$\mathcal{D} := \left\{ (\mu, \mu_T) \in \mathcal{M}_+(\mathfrak{B} \times U) \times \mathcal{M}_+(S(T)) \mid \begin{array}{l} \mu, \mu_T \ge 0, \quad \|\mu\| \le T, \quad \|\mu_T\| \le 1, \\ \text{(22) holds, } \int_{\mathfrak{B} \times U} |\zeta| d\mu \le C_\eta \end{array} \right\}. \quad (23)$$

Since the measure $d\mu(t, x, \zeta, u) = dt\delta_{x(t)}(d\nu(\zeta, u|t, x(t)))$ with $\pi_{\zeta\#} = \eta$ and $\pi_{u\#} = \omega$, satisfies the continuity equation (22) and the uniform first order moment bound on $\eta$, this establishes the inclusion $\mathfrak{R} \subset \mathcal{D}$. Subsequently, we show that the set $\mathcal{D}$ is convex and weak-\* compact. The convexity of set $\mathcal{D}$ follows from the linearity of the governing equation in $(\mu, \mu_T)$.

**Lemma 1.** The set $\mathcal{D}$ defined in (23) is weak-\* sequentially compact.

The proof can be found in Appendix A.

### D. Relaxation with occupation measures

In the following subsection, we present an occupation measure-based reformulation of the problem (2). We note that any occupation measure satisfying (22) can be disintegrated as $d\mu(t, x, \zeta, u) = dt d\mu(x|t) d\nu(\zeta, u|t, x)$ with $\nu \in \mathcal{P}(\mathcal{N}_{S(t)}(x) \times U)$ has conditional marginals $\pi_{\zeta\#} = \eta$, $\pi_{u\#} = \omega$. Using this occupation measure framework, the problem in (2) can be reformulated as follows:

$$J_o^*(x_0) = \inf_{\mu, \mu_T} \int_{\mathfrak{B} \times U} \mathcal{L}(x, u) d\mu(t, x, \zeta, u) \quad (24a)$$
subject to

$$
\begin{align}
\partial_t \mu(t, x, \zeta, u) + \nabla_x \cdot [(-\zeta + f(x, u))\mu(t, x, \zeta, u)] &= \delta_T \otimes \mu_T - \delta_0 \otimes \delta_{x_0}, \tag{24b} \\
\mu \in \mathcal{M}_+(\mathfrak{B} \times U), \quad \mu_T \in \mathcal{M}_+(S(T)), \nonumber
\end{align}
$$

where $\nabla_x \cdot$ is the divergence operator, and the dynamics
constraint in (24), has to be interpreted in the weak form,
as in (22), i.e., when integrated against test functions $\phi \in$
$C^1([0, T] \times \Omega)$. Note that the search for the minimizer $(\mu, \mu_T)$
in the above program is over the set $\mathcal{D}$, defined in (23),
which is weak-⋆ compact. Further, the cost function is weak-⋆
continuous (which can be checked using the fact that $\mathcal{L}(\cdot, \cdot)$
is continuous), so the existence of a minimizer follows from
the direct method of calculus of variations. Before stating the
main no-relaxation-gap theorem, we record a lemma used in
the proof of Theorem 3.

**Lemma 2.** Consider the sets $\mathcal{R}$ and $\mathcal{D}$ given in (20) and (23), respectively, and let $(\mu, \mu_T) \in \mathcal{D}$ be a feasible occupation measure pair satisfying (22). Then there exists a probability measure $\theta \in \mathcal{P}(\mathcal{R})$ such that

$$
\int_{\mathfrak{B} \times U} \mathcal{L}(x, u) \, d\mu = \int_{\mathfrak{R}} \left[ \int_{\mathfrak{B} \times U} \mathcal{L}(x, u) \, d\hat{\mu} \right] \, d\theta(\hat{\mu}) \\
\text{for every } \mathcal{L} : \mathbb{R}^n \times U \to \mathbb{R} \text{ satisfying Assumption 3.}
$$

The proof of this lemma is given in the Appendix B.

**Theorem 3.** For a fixed initial condition $x_0$, the optimal value $J_o^*(x_0)$ of occupation measure relaxed optimal control problem (24) and the optimal value $J_r^*(x_0)$ of Young's measure relaxed optimal control problem (16) are equal, i.e., there is no relaxation gap when relaxing to occupation measures.

Proof. One direction of the inequality, namely $J_o^*(x_0) \le J_r^*(x_0)$, holds trivially since the admissible set $\mathcal{R}$ for (16) is a subset of $\mathcal{D}$ which is an admissible set of (24).

Next, we show that there exists $\hat{\mu}$ such that the objective function value in both problems are equal. We consider an optimal pair $(\mu, \mu_T) \in \mathcal{D}$ that solves (24). By applying the superposition principle as stated in Lemma 2, we can represent the objective value of (24) as

$$
\int_{\mathfrak{B} \times U} \mathcal{L}(x, u) \, d\mu = \int_{\mathfrak{R}} \left[ \int_{\mathfrak{B} \times U} \mathcal{L}(x, u) \, d\hat{\mu} \right] \, d\theta(\hat{\mu}), \quad (25)
$$

where $(\hat{\mu}, \delta_{x(T)}) \in \mathcal{R}$ and $\theta \in \mathcal{P}(\mathcal{R})$ is a probability measure over admissible relaxed trajectories.

If it were the case that $\int \mathcal{L}(x, u) d\hat{\mu} > \int \mathcal{L}(x, u) d\mu$ for all $\hat{\mu}$ in the support of $\theta$, then the right-hand side of (25) would be strictly greater than the left-hand side, leading to a contradiction. Therefore, there must exist some $\hat{\mu} \in \mathcal{R}$ such that

$$
\int_{\mathfrak{B} \times U} \mathcal{L}(x, u) \, d\hat{\mu} = \int_{\mathfrak{B} \times U} \mathcal{L}(x, u) \, d\mu,
$$

which implies that the optimal values of both problems are
equal.
□

The generalization of problem (24) to account for an initial
measure other than a Dirac distribution concentrated at a point
$x_0$ is straightforward. Specifically, when the initial condition

is described by $\mu_0 \in \mathcal{P}(S(0))$, the objective and the dynamics
in (24) are integrated with respect to $\mu_0$. To generalize the
solution of (24), which depends on a fixed initial condition
$x_0$, we introduce average occupation measures by integrating
over an initial distribution $\mu_0$. Let us define the set

$$
\mathcal{E} := \{(t, \zeta, x, u) \mid t \in [0, T], x \in S(t), \\
\qquad \qquad \zeta \in \mathcal{N}_{S(t)}(x), u \in U\}.
$$

We then define the average occupation measure $\bar{\mu} \in \mathcal{P}(\mathcal{E})$ as

$$
\bar{\mu}(A) = \int_{S(0)} \mu(A \mid x_0) \, d\mu_0(x_0), \quad (26)
$$

for all measurable sets $A \subseteq \mathcal{C}$, where $\mu(\cdot | x_0) \in \mathcal{P}(\mathcal{C})$ is the
occupation measure associated with the initial condition $x_0$.
Similarly, the averaged terminal measure $\bar{\mu}_T$ is defined as

$$
\bar{\mu}_T(B) = \int_{S(0)} \mu_T(B | x_0) \, d\mu_0(x_0), \quad (27)
$$

for all measurable sets $B \subseteq S(T)$.

Since (24) corresponds to the special case of a Dirac initial
condition $\mu_0 = \delta_{x_0}$, we will use the same notation, $\mu$ and
$\mu_T$, to denote the occupation and terminal measures averaged
over the initial distribution $\mu_0$. With this notation, the optimal
control problem of a nonsmooth dynamical system becomes
the following linear program:

$$
\begin{align}
J_o^*[\mu_0] &= \inf_{\mu, \mu_T} \int \mathcal{L}(x, u) \, d\mu(t, x, \zeta, u) \nonumber \\
\text{s.t.} \quad & \partial_t \mu(t, x, \zeta, u) + \nabla_x \cdot [(\zeta + f(x, u))\mu(t, x, \zeta, u)] \nonumber \\
& = \delta_T \otimes \mu_T - \delta_0 \otimes \mu_0 \nonumber \\
& \qquad \mu \in \mathcal{M}_+(\mathfrak{B} \times U), \quad \mu_T \in \mathcal{M}_+(S(T)) \tag{28}
\end{align}
$$

where we recall $\mathfrak{B}$ from (18).

The dual program to (28) is described as follows:

$$
\begin{align}
D^*[\mu_0] &= \sup_v \int_{S(0)} v(0, x) d\mu_0(x) \nonumber \\
\text{s.t.} \quad & \mathcal{L} + \frac{\partial v}{\partial t} + \sum_{i=1}^{n} \frac{\partial v}{\partial x_i} (-\zeta + f)_i \in \mathcal{C}_+(\mathfrak{B} \times U) \nonumber \\
& v(T, \cdot) \in \mathcal{C}_+(S(T)). \tag{29}
\end{align}
$$

Using the arguments in [29, Theorem 7.2, Chapter 4], it can
be shown that strong duality holds between (28) and (29), that
is, $J_o^*[\mu_0] = D^*[\mu_0]$.

III. RELAXED DISCRETE-TIME OPTIMAL CONTROL

In this section, we relax the discrete-time optimal control problem (5) with dynamics (4) to the space of measures. We begin with the case where the constraint set varies with time. For this setting, we formulate an equivalent linear program over measures, given in (38), that captures the evolution of the state distribution across time steps. We then specialize our construction to the static-set case. Under the assumption that S is time-independent, we derive in Lemma 3 a discrete-time counterpart of the continuity equation in (43). Theorem 4 then establishes that for any solution of (43), there exists a stochastic kernel whose propagation through (4) generates the
corresponding sequence of state distributions. Together, these
results enable a reformulation of the optimal control problem
with a static constraint set, see (45), as a linear program in
the space of measures in which the only unknowns are the
occupation measure $\mu^\tau$ and the terminal measure $\mu_T^\tau$ (with
$\tau > 0$ the time step), in contrast to (38), which introduces a
distinct measure variable for each time instant.

We consider the uniform partition of the time interval $[0, T]$

$$
\{0 = t_0, t_1, \dots, t_N = T\}, \quad t_k = k\tau, \quad N = T/\tau \in \mathbb{N}. \tag{30}
$$

We know from (4) that the states at successive time instants
$k\tau$ and $(k+1)\tau$ are related as

$$
x_{k+1}^{\tau} = P_{S_{k+1}} \circ (\tau f(x_k^{\tau}, u_k^{\tau}) + x_k) := G_k^{\tau}(x_k, u_k) \quad (31)
$$

where $P_{S_k}$ is as before the projection mapping on the compact
convex set $S_k := S(t_k)$. To relax the problem in (5) to the
space of measures we will rely on probabilistic arguments
since we are modeling the initial condition and the control as
random variables [16]. For a fixed value of $\tau$, we define the
measures $\{\tilde{\mu}_k^\tau\}_{k\in\mathbb{N}}$, such that $\tilde{\mu}_k^\tau \in \mathcal{P}(S_k)$, at time instants $t_k$
which are a result of recursive evolution of measures through
dynamics (31). We introduce a time-varying stochastic kernel,
also called Markov control policy, $\omega_k^\tau(u|x) \in \mathcal{P}(U)$, which
is the probability measure on controls $u$ at time instant $k$
conditioned over the state $x$. Using this stochastic kernel we
define a transition kernel which is the probability measure of
states at time $k\tau$ given the state was at $x_{k-1}$ at time $(k-1)\tau$:

$$
Q_{\omega}^{\tau}(A|x_{k-1}, k-1) \triangleq \int_{U} \mathbb{I}_{A}(G_{k-1}^{\tau}(x_{k-1}, u))d\omega_{k-1}^{\tau}(u|x_{k-1}) \quad (32)
$$

where, $\mathbb{I}_A$ denotes the indicator function on set $A \subset S_k$.
Kernel $Q_\omega^\tau(\cdot|\cdot, \cdot)$ captures the effect of feedback $\omega_k^\tau(\cdot|\cdot)$ at
each time step. Now, given $\tilde{\mu}_0^\tau$ is the measure at time 0, we can
compute the measure at the next time step, $\tilde{\mu}_1^\tau$, as its successor.
This is given by:

$$
\tilde{\mu}_1^\tau(A) = \int_{S_0} \int_U \mathbb{I}_A(G_0^\tau(x, u)) d\omega_0^\tau(u|x) d\tilde{\mu}_0^\tau(x) \quad (33)
$$

where $A \subset S_1$. We can recursively compute the measures at
time $k$ using the following equation:

$$
\tilde{\mu}_k^\tau(A) = \int_{S_{k-1}} \int_U \mathbb{I}_A(G_{k-1}^\tau(x, u)) d\omega_{k-1}^\tau(u|x) d\tilde{\mu}_{k-1}^\tau(x). \quad (34)
$$

Next, we identify a measure $\mu_k^\tau \in \mathcal{P}(S_k, U)$

$$
\mu_k^\tau(dx, du) := \omega_k^\tau(du|x) \tilde{\mu}_k^\tau(dx) \quad (35)
$$

at each time instant. We also note that $\tilde{\mu}_k^\tau$ is the push-forward
of $\mu_{k-1}^\tau$ through $G_k^\tau$. Therefore, the following relationships
hold:

$$
x_k^\tau \xleftarrow{G_{k-1}^\tau} (x_{k-1}^\tau, u_{k-1}^\tau), \quad (36)
$$

$$
\tilde{\mu}_k^\tau(dx) \xleftarrow{G_{k-1}^\tau \#} \mu_{k-1}^\tau(dx, du). \quad (37)
$$

Let $\pi_x$ be a projection mapping such that $\pi_x : \mathbb{R}^n \times U \rightarrow \mathbb{R}^n$.
Using (37), we have the following linear program in the space
of Borel measures for (5):

$$
J^{*\tau}[\mu_0] = \inf_{\mu_k^\tau} \sum_{k=1}^{N} \tau \int_{S_k \times U} \mathcal{L}(x, u) d\mu_k^\tau(x, u) \quad (38)
$$

$$
\begin{equation} \tag{39}
\begin{split}
\text{s.t. } \pi_{x\#} \mu_k^\tau &= [P_{S_k} \circ (\tau f(\cdot, \cdot) + \pi_x)]_{\#} \mu_{k-1}^\tau, \\
\mu_k^\tau &\in \mathcal{P}(S_k \times U), \quad \forall k = 1, \dots, N.
\end{split}
\end{equation}
$$

The dual problem to the above can be derived as follows:

$$
\begin{equation}
\begin{split}
J_d^{*\tau}[\mu_0] &= \max_{V_i(\cdot)} \int_{S_0} V_0(x) d\mu_0(x) \\
\text{s.t.} \quad &V_i \le \tau \mathcal{L} + V_{i+1}(G_k^\tau(x, u)) \quad \text{for } i = 0, \dots, N-1, \\
&V_N = 0 \quad \text{on } S_N
\end{split}
\tag{40}
\end{equation}
$$

where, $V_i \in C_b(\mathbb{R}^n)$. Following the techniques used in [29, Theorem 7.2, Chapter 4], it can be shown that, for a given $\mu_0 \in \mathcal{P}(S_0)$, there is no duality gap between (38) and (40), i.e., $J_d^{*\tau}[\mu_0] = J^{*\tau}[\mu_0]$.

Next, we introduce a discrete-time continuity equation, which has a similar structure to (22), but applies when the set is static, i.e., $S(t) = S$. In this case, the update rule becomes

$$
G^{\tau}(x_k, u_k) := P_S \circ (\tau f(x_k, u_k) + x_k). \quad (41)
$$

We will later highlight that considering this specific case
separately can be numerically advantageous, especially when
the set defining the normal cone constraint remains static. For
the following discussion, we define discrete-time occupation
measure $\mu^\tau \in \mathcal{M}_+(S \times U)$ satisfying the following relation-
ship

$$
\int_{A \times B} \phi(x, u) d\mu^\tau(x, u) := \sum_{k=0}^{N-1} \int_{A \times B} \phi(x, u) d\mu_k^\tau(x, u) \quad (42)
$$

for all $\phi(x, u)$ bounded measurable function and $A \subset S$, $B \subset U$. This measure captures the time spent by all possible trajectories in some subset of state space and admissible control values.

**Lemma 3.** Let $\omega_k^\tau(\cdot|\cdot) \in \mathcal{P}(U)$ be a stochastic kernel and $\mu_k^\tau$, $\mu^\tau$ be as defined in (35), (42). Let $\tilde{\mu}_0^\tau$, $\tilde{\mu}_T^\tau$ be such that $x_0 \sim \tilde{\mu}_0^\tau$ and $x_N \sim \tilde{\mu}_T^\tau$ respectively. Then, $\mu^\tau$, $\tilde{\mu}_0^\tau$, $\tilde{\mu}_T^\tau$ satisfy the discrete-time continuity equation,

$$
\begin{align}
& \int_{S \times U} v(x) d\mu^\tau(x, u) + \int_{S} v(x) d\tilde{\mu}_T^\tau(x) \nonumber \\
& = \int_{S \times U} v(G^\tau(x, u)) d\mu^\tau(x, u) + \int_{S} v(x) d\tilde{\mu}_0^\tau(x) \tag{43}
\end{align}
$$

for all $v(\cdot) : \mathbb{R}^n \to \mathbb{R}$ bounded, Borel-measurable functions and $G^\tau$ defined in (41).
**Proof.** Let $v(x)$ be some bounded measurable function, then using (42) and (35),

$$ 
\begin{align*}
\int_{S \times U} v(x) d\mu^\tau(x, u) &= \sum_{k=0}^{N-1} \int_{S \times U} v(x) d\mu_k^\tau(x, u) \\
&= \sum_{k=0}^{N-1} \int_{S \times U} v(x) d\omega_k^\tau(u|x) d\tilde{\mu}_k^\tau(x) = \sum_{k=0}^{N-1} \int_S v(x) d\tilde{\mu}_k^\tau(x) \\
&= \sum_{k=0}^{N-1} \int_S v(x) dG_\#^\tau \mu_{k-1}^\tau(x, u)
\end{align*}
$$

where in the third and fourth equality we have used the fact that $\omega_k^\tau(\cdot|x) \in \mathcal{P}(U)$ and (37) respectively. Now using the change of variables formula for push-forward measures we get,

$$ 
\int_{S \times U} v(x) d\mu^\tau(x, u) = \sum_{k=0}^{N-1} \int_{S \times U} v(G^\tau(x, u)) d\mu_k^\tau(x, u) \\ 
+ \int_S v(x) d\tilde{\mu}_0^\tau(x) - \int_S v(x) d\tilde{\mu}_T^\tau(x). \quad \square
$$

Next, we prove that for any solution of (43), we have a stochastic kernel such that if we propagate this stochastic kernel through (31), we obtain a sequence of measures defining the probability of state at each time instant.

**Theorem 4.** Let $\tilde{\mu}_0^\tau, \tilde{\mu}_T^\tau, \mu^\tau$ be the measures which satisfy discrete-time continuity equation (43). Then there exists a stochastic kernel $\omega_k^\tau(\cdot|x) \in \mathcal{P}(U)$ at each time $k\tau$ which defines the evolution of measure $\tilde{\mu}_0^\tau$ through (34) to $\tilde{\mu}_T^\tau \in \mathcal{P}(S)m$ and $\mu^\tau \in \mathcal{M}_+(S \times U)$ is the corresponding discrete-time occupation measure satisfying (42).

**Proof.** Let us consider nonnegative measurable functions $p_k \in L^1(\mu^\tau)$ such that

$$ 
\sum_{k=0}^{N-1} p_k(x, u) = 1 (\mu^\tau \text{-a.s.}), \quad \int_{S \times U} p_k(x, u) d\mu^\tau(x, u) = 1, \\ 
\forall k = 0, \dots, N-1.
$$

Define $\mu_k^\tau := p_k \mu^\tau$. Then, for any measurable $A \subset S$ and $B \subset U$,

$$ 
\int_{S \times U} \mathbb{I}_{A \times B}(x, u) d\mu^\tau(x, u) = \\ 
\int_{S \times U} \mathbb{I}_{A \times B}(x, u) \sum_{k=0}^{N-1} p_k(x, u) d\mu^\tau(x, u) = \\ 
\sum_{k=0}^{N-1} \int_{S \times U} \mathbb{I}_{A \times B}(x, u) d\mu_k^\tau(x, u)
$$

where in the last equality we have used the definition of $\mu_k^\tau$. The set of such $\{p_k\}_k$ is nonempty as $p_k = \frac{1}{N}$ for $k = 0, \dots, N-1$ is a feasible point. Using [30, Corollary

10.4.13], each $\mu_k^\tau$ can be disintegrated as $\mu_k^\tau(dx, du) = \omega_k^\tau(du|x) \tilde{\mu}_k^\tau(dx)$, then (43) becomes

$$ 
\begin{align*}
\int_{S \times U} \mathbb{I}_A(x) d\mu^\tau(x, u) + \tilde{\mu}_T^\tau(A) \\
= \sum_{k=0}^{N-1} \int_{S \times U} \mathbb{I}_A(G^\tau(x, u)) d\omega_k^\tau(u|x) d\tilde{\mu}_k^\tau(x) + \tilde{\mu}_0^\tau(A). \quad (44)
\end{align*}
$$

Using $\omega_k^\tau$ as the time-varying stochastic kernel we can define $Q_\omega^\tau(A|x_k, k)$ using (32). Thus (44) can be equivalently written as,

$$ \mu^\tau(A, U) + \tilde{\mu}_T^\tau(A) = \sum_{k=0}^{N-1} \int_S Q_\omega^\tau(A|x_k, k) d\tilde{\mu}_k^\tau(x_k) + \tilde{\mu}_0^\tau(A). $$

$$ \text{So, } \mu^\tau(A, U) = \sum_{k=0}^{N} \tilde{\mu}_k^\tau(A) - \tilde{\mu}_T^\tau(A) = \sum_{k=0}^{N-1} \tilde{\mu}_k^\tau(A). \quad \square $$

Using Lemma 3 and Theorem 4, we define the following linear program in the cone of non-negative Borel measures,

$$ 
\begin{align}
J^{*\tau}[\mu_0] &= \inf_{\mu^\tau, \tilde{\mu}_T^\tau} \int \mathcal{L}(x, u) d\mu^\tau(x, u) \nonumber \\
\text{s.t.} \quad &\pi_{x\#} \mu^\tau + \tilde{\mu}_T^\tau = [P_S \circ (\tau f(\cdot, \cdot) + \pi_x)]_\# \mu^\tau + \tilde{\mu}_0^\tau, \nonumber \\
&\mu^\tau \in \mathcal{P}(S \times U), \quad \tilde{\mu}_T^\tau \in \mathcal{P}(S). \tag{45}
\end{align}
$$

Here, the dynamics constraint in (45) is interpreted in weak sense, as in (43). As we observe in the above program, we only solve for $\mu^\tau$ and $\tilde{\mu}_T^\tau$, which are the aggregated occupation measure and the terminal measure respectively. The derivation in the proof of Lemma 3 requires the push-forward map $G^\tau$ to be time-invariant. For the case with time-varying set as (5), we cannot derive an equation analogous to (43) rather we directly use the relations defined in (36) to obtain (39).

## IV. CONVERGENCE OF THE DISCRETE-TIME OPTIMAL CONTROL PROBLEM TO THE CONTINUOUS TIME PROBLEM

Having defined the relaxation for the continuous-time problem in (28) and the discrete-time problem in (38), in this section, we study the convergence $J^{*\tau} \xrightarrow{\tau \to 0} J_o^*$, where $J^{*\tau}$ and $J_o^*$ are defined in (38) and (28) respectively, and $\tau$ is the time step. There are three main ingredients of this convergence problem: (1) Construction of suitable interpolations for the admissible solutions to (38), (2) Limiting behavior of the interpolated curves as $\tau \to 0$, and (3) Convergence of the optimal value $J^{*\tau}$ of (38) towards the optimal value $J_o^*$ of (28) as $\tau \to 0$. Before addressing each of these steps, we present some important technical definitions and tools that will be used in the various parts of the proof.

### A. Wasserstein space

We will denote $\mathcal{P}_2(\mathbb{R}^n)$ as the space of measures with finite second order moment. In such measure spaces, Wasserstein metric is commonly used to evaluate the distance between two measures. It is defined as

$$ W_2(\mu, \nu) := \min_{\theta \in \Theta(\mu, \nu)} \left( \int_{\mathbb{R}^n \times \mathbb{R}^n} |x - y|^2 d\theta(x, y) \right)^{1/2} $$
where $\Theta(\mu, \nu)$ is the set of joint probability measures on $\mathbb{R}^n \times \mathbb{R}^n$ with given marginals $\mu$ and $\nu$, or in other words $\theta \in \Theta(\mu, \nu)$ satisfies $\theta(A, \mathbb{R}^n) = \mu(A)$ and $\theta(\mathbb{R}^n, A) = \nu(A)$, for every measurable $A \subset \mathbb{R}^n$. We refer to $(\mathcal{P}_2(\mathbb{R}^n), W_2)$ as the (2-)Wasserstein space; it is a geodesic metric space, so for any $\mu_0, \mu_1 \in \mathcal{P}_2(\mathbb{R}^n)$, there exists a constant-speed geodesic $(\mu_t)_{t \in [0,1]}$ with $W_2(\mu_s, \mu_t) = |t-s|W_2(\mu_0, \mu_1)$. Every such geodesic is absolutely continuous and solves the continuity equation $\partial_t \mu_t + \nabla \cdot (v_t \mu_t) = 0$ with $\|v_t\|_{L^2(\mu_t)} = |\mu|(t) = W_2(\mu_0, \mu_1)$ a.e. $t$, for a (minimal-norm) velocity field $v_t \in L^2(\mu_t; \mathbb{R}^n)$ (see [31, Theorem 8.3.1]).

## B. Interpolated curves

Consider the partition of the time interval $[0, T]$ introduced in (30). Recall from (37) and (34) that $\tilde{\mu}_k^\tau$ is the image measure of $\mu_{k-1}^\tau$ under the push-forward map $P_{S_k} \circ (\tau f(\cdot, \cdot) + \text{id}_x)$. Equivalently, for every Borel set $A \subset \mathbb{R}^n$,

$$ \tilde{\mu}_k^\tau(A) = \int_{A, U} \left[ P_{S_k} \circ \left( \tau f(y, u) + \pi_x \right) \right]_{\#} d\mu_{k-1}^\tau(y, u). $$

Using the disintegration of $\mu_{k-1}^\tau$ from (35), we obtain:

$$ \tilde{\mu}_k^\tau = \left[ P_{S_k} \circ \left( \tau \int_U f(\cdot, u) d\tilde{w}(u \mid k\tau, \cdot) + \pi_x(\cdot) \right) \right]_{\#} \tilde{\mu}_{k-1}^\tau. $$

In the following, we denote:

$$ G_k^\tau(x) := P_{S_k} \left( \tau \int_U f(x, u) d\tilde{w}(u \mid k\tau, x) + x \right). $$

Thus, at each time $k\tau$ for $k \in \{0, 1, \dots, \lceil \frac{T}{\tau} \rceil\}$, we have:

$$ \tilde{\mu}_{k+1}^\tau = G_k^\tau \# \tilde{\mu}_k^\tau, \quad (46) $$

$$ v_{k+1}^{\tau} = \frac{G_{k}^{\tau}(x) - x}{\tau}. \quad (47) $$

Given the discrete sequence $\{\mu_k^\tau\}_{k=0}^{[T/\tau]}$, we associate two interpolations that will be crucial in the convergence analysis.

(1) McCann (constant speed) interpolation: Let us define a constant speed geodesic interpolated curve $\tilde{\mu}_t^\tau$ as

$$ \tilde{\mu}_t^\tau := (tx + (1-t)y)_{\#} \theta(dx, dy) \quad \text{for } t \in (k\tau, (k+1)\tau] \quad (48) $$

such that $\theta$ is the optimal transport plan between $\tilde{\mu}_{k+1}^\tau$ and $\tilde{\mu}_k^\tau$. The velocity at time $t$ is defined as

$$ \tilde{v}_t^\tau := v_{k+1}^\tau \circ (G_t^\tau)^{-1}(x) \text{ for } t \in (k\tau, (k+1)\tau]. $$

One important characterization of these curves [32, Chapter 5] is that,

$$ \|\tilde{v}_t^\tau\|_{L^2(\tilde{\mu}_t^\tau)} = \frac{W_2(\tilde{\mu}_k^\tau, \tilde{\mu}_{k+1}^\tau)}{\tau} = |(\tilde{\mu}^\tau)'|(t) \quad \text{for all } t \in (k\tau, (k+1)\tau]. \quad (49) $$

This interpolated curve satisfies the continuity equation

$$ \partial_t \tilde{\mu}_t^\tau + \nabla \cdot (\tilde{E}_t^\tau) = 0 $$

where $(\tilde{E}_t^\tau) := \tilde{v}_t^\tau \tilde{\mu}_t^\tau$.

(2) Piecewise constant interpolation curves such that

$$ \hat{\mu}_t^\tau = \tilde{\mu}_{k+1}^\tau, \quad \text{for } t \in (k\tau, (k+1)\tau] \quad (50) $$

and velocity $\hat{v}_t^\tau = v_{k+1}^\tau$ for $t \in (k\tau, (k+1)\tau]$.

We also define the corresponding momentum vector $\hat{E}_t := \hat{v}_t^\tau \hat{\mu}_t$.

In what follows, we will use the following notation

$$ g^{\tau}(k\tau, x) = \int_{U} f(x, u) \, d\omega^{\tau}(u|k\tau, x). \quad (51) $$

**Remark 1.** Two interpolations are used for distinct but complementary purposes. The convergence of the curve $\mu_t^\tau$ ensures convergence to the limit $\mu_t$, which solves the continuity equation. On the other hand, the convergence of velocities of the piecewise constant interpolated curve $\hat{\mu}_t^\tau$ ensures that the limiting velocity field of the continuity equation satisfies

$$ v_t(x) \in -\mathcal{N}_{S(t)}(x) + \int_U f(x, u) \, d\omega(u \mid t, x). $$

## C. Convergence of dynamics constraint

Next, we will show that for $t \in [0, T]$, the interpolated trajectories $\mu_t^\tau$ and $\hat{\mu}_t^\tau$ converge to $\mu_t$ which is the absolutely continuous solution to the continuity equation

$$ \partial_t \mu_t + \nabla \cdot (v_t \mu_t) = 0 \quad (52) $$

with vector field $v_t: \mathbb{R}^n \to \mathbb{R}^n$ defined in (9).

**Theorem 5.** For system (1) and measure $\mu_0 \in \mathcal{P}(S(0))$, suppose that Assumptions 1 and 3 hold. For $t \in [0, T]$, the curves $\{\tilde{\mu}_t^\tau\}_{t \in [0, T]}$, $\{\hat{\mu}_t^\tau\}_{t \in [0, T]}$ defined in (48) and (50), respectively, converge uniformly, with respect to the $W_2$ metric, to $\{\mu_t\}_{t \in \mathcal{P}(S(t))}$ which satisfies (52). The momentum vectors $\hat{E}_t^\tau = v_t^\tau \hat{\mu}_t^\tau$ and $\hat{E}_t^\tau = \hat{v}_t^\tau \hat{\mu}_t^\tau$ also weak-\* converge to $E_t = v_t \mu_t$ where, $v_t(x) \in -\mathcal{N}_{S(t)}(x) + \int_U f(x, u) d\omega(u|t, x)$.

*Proof.* The proof is divided in four parts: (1) Proof of convergence of $\tilde{\mu}_t^\tau \to \mu_t$ and $\hat{\mu}_t^\tau \to \mu_t$; (2) Proof of convergence of $\hat{E}_t^\tau$ and $\hat{E}_t^\tau$ to $E_t$; (3) Proof of absolute continuity of $E_t$ w.r.t. $\mu_t$ such that $E_t = v_t \mu_t$; (4) Proof of convergence of $\hat{v}_t^\tau$ to $v_t$ such that $v_t \in -\mathcal{N}_{S(t)}(x) + \int_U f(x, u) d\omega(u|t, x)$.

In order to prove the first part, we need the following lemma; its proof can be found in [27, Lemma 2]. Let $C_{\max} := \sup_{\Omega} \|x_k\|$. By [1, Section 5], $C_{\max} < \infty$; it depends on the system data and the initial state $x_0$, and is independent of time instants $k$.

**Lemma 4.** Consider the definitions of $\tilde{\mu}_k^\tau$ and $\tilde{\mu}_{k+1}^\tau$ as given in (46). Then, the following inequality holds:

$$ W_2(\tilde{\mu}_k^\tau, \tilde{\mu}_{k+1}^\tau) \le \tau(L_f C_{\max} + L_s) \quad (53) $$

where $L_f$ is defined in Assumption 1, and $L_s$ is defined in Assumption 2.

(1) *Proof of convergence of $\tilde{\mu}_t^\tau \to \mu_t$ and $\hat{\mu}_t^\tau \to \mu_t$:* Using Lemma 4, we can derive the following bound (see [27] for details),

$$ W_2(\tilde{\mu}_t^\tau, \tilde{\mu}_s^\tau) \le (t-s)^{1/2} (L_f C_{\max} + L_s) T^{1/2} \quad \forall s, t \in [0, T]. \quad (54) $$
The same estimate holds for the curve $\hat{\mu}_t$; specifically; $W_2(\hat{\mu}_t^\tau, \hat{\mu}_s^\tau)$ has the same bound as (54). Thus, the curves $\hat{\mu}_t^\tau$, and $\tilde{\mu}_t^\tau$ are $\frac{1}{2}$-Hölder continuous. Using the Arzelà-Ascoli theorem, we conclude that there exists a subsequence $\{\tau_j\}_{j}$ for which $\tilde{\mu}_t^\tau \to \tilde{\mu}_t$ uniformly in $W_2$ space and the limit curve $\tilde{\mu}_t$ is absolutely continuous. Same conclusion holds for curves $\hat{\mu}_t^\tau$ i.e., $\hat{\mu}_t^{\tau j} \to \hat{\mu}_t$.

Moreover, observe that the curve $\hat{\mu}_t^\tau$ coincides with $\tilde{\mu}_t^\tau$ at times $k\tau$ and is constant on the interval $(k\tau, (k+1)\tau]$. Thus, we conclude from (54) that $W_2(\hat{\mu}_t^\tau, \tilde{\mu}_t^\tau) \le \tau^{1/2}(L_f C_{\max} + L_s)T^{1/2}$. So, both curves converge to the same limit curve as $\tau \to 0$, which we denote by $\mu_t$.

(2) Convergence of $\tilde{E}_t^\tau$ to $E_t$: We define $\tilde{m}^\tau \in \mathcal{M}^n([0, T] \times \mathbb{R}^n)$ as

$$ \tilde{m}^{\tau} = \tilde{v}_{t}^{\tau} \tilde{\mu}_{t}^{\tau} dt \quad (55) $$

which is a vector-valued measure.<fnref n="3" />

**Lemma 5.** An estimate on the norm of $\tilde{m}^\tau$ defined in (55) is as follows:

$$ |\tilde{m}^{\tau}|([0, T] \times \Omega) \le T^{\frac{3}{2}} (L_f C_{\max} + L_s). \quad (56) $$

where $C_{\max}$ represents the uniform bound on $|x_k|$, $L_f$ is defined in Assumption 1, and $L_s$ is defined in Assumption 2.

The proof of the estimate (56) can again be found in [27, Lemma 3]. Using the same proof, one can obtain the same bounds on $|\hat{m}^\tau|$. The sequences $\tilde{m}^\tau$ and $\hat{m}^\tau$ are uniformly bounded and therefore relatively compact in the weak-\* topology. As a result, we have (up to a subsequence) $\tilde{m}^\tau \xrightarrow{\ast} \hat{m}$, $\hat{m}^\tau \rightrightarrows \hat{m}$, $\tilde{E}_t^\tau \rightrightarrows \tilde{E}_t$ and $\tilde{E}_t^\tau \rightrightarrows \tilde{E}_t$. Finally, by applying [32, Lemma 8.9], we conclude that $\tilde{m} = \hat{m} =: m$ and thus $\tilde{E}_t = \hat{E}_t =: E_t$.

(3) Absolute continuity of $E_t$ w.r.t. $\mu_t$: To prove the absolute continuity of $E_t$ we use the properties of Benamou-Brenier functional $\mathcal{B}(\mu dt, E dt)$ (for more details refer [32, Chapter 5]). A key property of $\mathcal{B}(\tilde{\mu}_t^\tau dt, \tilde{E}_t^\tau dt)$ which is of importance here is the following: $\tilde{E}_t^\tau \ll \tilde{\mu}_t^\tau$ such that $\tilde{E}_t^\tau = \tilde{v}_t^\tau \tilde{\mu}_t^\tau$, then $\mathcal{B}(\tilde{\mu}_t^\tau dt, \tilde{E}_t^\tau dt) = \int_0^T \int_\Omega |\tilde{v}_t^\tau|^2 d\tilde{\mu}_t^\tau dt$. Now the uniform bound on $|m^\tau|$ also implies an uniform bound on $\int_0^T \int_\Omega |v_t^\tau|^2 d\tilde{\mu}_t^\tau dt$ from (56). Moreover, nonnegativity, convexity and lower semi-continuity of $\mathcal{B}(\cdot, \cdot)$ [32, Proposition 5.18] implies

$$ \mathcal{B}(\mu_t dt, dm) \le \liminf_{k \to \infty} \mathcal{B}(\tilde{\mu}_t^\tau dt, d\tilde{m}^\tau) < \infty, $$

where $\mu_t$ and $m$ are the weak-\* limit of $\tilde{\mu}_t$ and $\tilde{m}^\tau$ as $\tau \to 0$. The finiteness of $\mathcal{B}(\mu_t dt, m)$ implies that $E_t \ll \mu_t$ and $E_t = v_t \mu_t$ with $v_t \in L^2(\mu_t)$. Applying the same argument to the sequence $(\hat{\mu}_t^\tau, \hat{E}_t^\tau)$ gives $\hat{E}_t = \hat{v}_t \hat{\mu}_t$. From the previous step we already have $\hat{\mu}_t = \mu_t$ and $\hat{E}_t = E_t$; hence $\hat{v}_t = v_t \mu_t$-a.e.

(4) Proof that $v_t(x) \in -\mathcal{N}_{S(t)}(x) + \int_U f(x, u) d\omega(u|t, x)$: To prove this statement, we consider the sequence of curves

$\hat{v}_t^\tau$. Using (47) for $g^\tau(k\tau, x)$ (defined in (51)) and for every $x \in S_k$, we have

$$ \left\langle y - P_{S_{k+1}}(\tau g^\tau(k\tau, x) + x), \frac{(P_{S_{k+1}} \circ (\tau g^\tau(k\tau, \cdot) + \text{id})x - x)}{\tau} \right\rangle \ge 0 $$

which follows from the construction of the discrete-time velocities $\hat{v}_k^\tau$. Then

$$ \langle y - P_{S_{k+1}} \circ (\tau g^\tau(k\tau, \cdot) + \text{id})x, v_{k+1}^\tau(x) - g^\tau(k\tau, x) \rangle \ge 0, $$

where the brackets denote the inner product between vectors in $\mathbb{R}^n$. This is the condition for $v_{k+1}^\tau(x) - g^\tau(k\tau, x)$ to be in the normal cone to set $S_{k+1}$. In the integral form the above condition is expressed as

$$ \int_{\Omega} \left\langle h(t_k, x)(y - P_{S_{k+1}} \circ (\tau g^\tau(k\tau, \cdot) + \text{id})x), (v_{k+1}^\tau(x) - g^\tau(k\tau, x))d\hat{\mu}_{k+1}^\tau(x) \right\rangle \ge 0, \quad \forall y \in S_{k+1} \quad (57) $$

for any positive $h(t, x)$. We define the interpolation of $g^\tau(k\tau, x)$ as

$$ g^\tau(t, x) = g^\tau(k\tau, x), \quad \text{for } t \in (k\tau, (k+1)\tau]. \quad (58) $$

Using the piecewise constant interpolations defined in (50), we have

$$ \begin{aligned} d\hat{\mu}_t^\tau &= d\hat{\mu}_{k+1}^\tau \\ d\hat{E}_t^\tau &= \hat{v}_t^\tau d\hat{\mu}_t^\tau = v_{k+1}^\tau d\hat{\mu}_{k+1}^\tau \end{aligned} $$

for $t \in (k\tau, (k+1)\tau]$. In order to study the convergence of (57), we will need the following result.

**Lemma 6.** For any positive measurable function $h(t, x)$ we have

$$ \lim_{\tau \to 0} \int h(t, x)(y - P_{S(t)}(x))g^\tau(t, x)d\hat{\mu}_t^\tau(x) = \int h(t, x)(y - x)g(t, x)d\mu_t(x) $$

where $g^\tau$ is defined in (58), $\hat{\mu}_t^\tau$ defined in (50) and $t \in (k\tau, (k+1)\tau]$.

The proof can be found in Appendix C.

Using Lemma 6 and the convergence results for $\hat{\mu}_t^\tau \to \mu_t$, $\hat{E}_t^\tau \to E_t$, for any $h(t, x)$ positive function we have,

$$ \int_0^T \int_\Omega \langle h(t, x)(y-x), dE_t(x)dt \rangle - \int_0^T \int_\Omega \langle h(t, x)(y-x), g(t, x)d\mu_t(x)dt \rangle \ge 0, \quad \forall y \in S(t). $$

Since $h(t, x)$ is an arbitrary positive function we get,

$$ \langle y - x, v(t, x) - g(t, x) \rangle \ge 0, \quad \forall y \in S(t). $$

Thus, $v(t, x) - \int_U f(x, u)d\omega(u|t, x) \in -\mathcal{N}_{S(t)}(x)$ which can be re-expressed in terms of some selection $\eta \in \mathcal{P}(\mathcal{N}_{S(t)}(x))$ as,

$$ v(t, x) = - \int_{\mathcal{N}_{S(t)}(x)} \zeta d\eta(\zeta|t, x) + \int_U f(x, u)d\omega(u|t, x). $$

<footnote n="3">³The space of vector-valued measures $\mathcal{M}^n(\Omega)$ is dual to $C_0(\Omega; \mathbb{R}^n)$ and is endowed with weak-\* convergence with respect to this duality.</footnote>
This completes the proof of Theorem 5.

<figure><img src="image_4.png" /></figure>

As a result of Theorem 5, we obtain the continuity equation in the following form:

$$ \int_{S(T)} \varphi(T, x) d\mu_T(x) - \int_{S(0)} \varphi(0, x) d\mu_0(x) = \int_{[0,T] \times S(t)} \int_U \left( \partial_t \varphi(t, x) + \partial_x \varphi(t, x) \cdot \left[ -\int_{N_{S(t)}(x)} \zeta \, d\eta(\zeta|t, x) \right. \right. \\ \left. \left. + \int_U f(x, u) d\omega(u|t, x) \right] \right) d\tilde{\mu}_t(x) dt, \quad (59) $$

for all test functions $\varphi \in C^1([0, T], \Omega)$. Moreover, Lemma 2 shows that any solution $\mu$ to the problem (24) admits a disintegration as $d\mu(t, x, \zeta, u) = d\nu(\zeta, u|t, x)d\tilde{\mu}_t(x)dt$ such that $\nu \in \mathcal{P}(N_{S(t)}(x) \times U)$ has marginals $\pi_{\zeta \#} \nu = \eta$, $\pi_{u \#} \nu = \omega$ and the resulting continuity equation is (74) with the vector being (75) which is equivalent to (59).

### D. Convergence of value function

In the previous subsection, we proved that for any sequence of feasible solutions $\mu^{\tau_j}$, to the family of problems (38), we obtain a feasible solution $\mu$ of (28). We now establish convergence of the optimal values as $\tau_j \to 0$. For point-mass initial data $\mu_0 = \delta_{x_0}$, we write $J^{\tau_j^*}(x_0) := J^{\tau_j^*}[\mu_0]$ and $J_o^*(x_0) := J_o^*[\mu_0]$ for the optimal values of the discrete and continuous problems, respectively.

**Theorem 6.** Let $\mu_0 = \delta_{x_0}$ with $x_0 \in S(0)$ be fixed. The optimal value of discrete-time LP (38) converges to the optimal value of continuous-time LP (28) as $\tau_j \to 0$, i.e., $J^{\tau_j^*}(x_0) \to J_o^*(x_0)$. In other words, $J^{\tau_j^*}$ converges pointwise to $J_o^*$ for fixed $x_0$.

*Proof.* The discrete-time objective function in (38) is

$$ V^{\tau_j}[\omega^{\tau_j}, \mu_0] = \sum_{k=0}^{N} \tau_j \int_{S_k \times U} \mathcal{L}(x, u) d\omega^{\tau_j}(u|k\tau, x) d\tilde{\mu}_k^{\tau_j}(x). $$

We first prove that the objective function value $V^{\tau_j}$ for each admissible solution $\omega^{\tau_j}$ converges to the objective function value of (28) obtained for $\omega = \lim_{\tau_j \to 0} \omega^{\tau_j}$. We study this convergence for $\mu_0 = \delta_{x_0}$. Using the convergence results proved in Theorem 5 we have

$$ \begin{aligned} & \sum_{k=0}^{N} \tau_j \int_{S_k \times U} \mathcal{L}(x, u) d\omega^{\tau_j}(u|k\tau_j, x) d\tilde{\mu}_k^{\tau_j}(x) \\ &= \sum_{k=0}^{N} \int_{S_k \times U} \mathcal{L}(x, u) d\omega^{\tau_j}(u|k\tau_j, x) d\tilde{\mu}_k^{\tau_j}(x) \tau_j \\ & \xrightarrow{\tau_j \to 0} \int_0^T \int_{S(t) \times U} \mathcal{L}(x, u) d\omega(u|t, x) d\tilde{\mu}_t(x) dt \\ &= \int_0^T \int_{S(t) \times U} \int_{N_{S(t)}(x)} \mathcal{L}(x, u) d\mu(t, x, \zeta, u) \end{aligned} $$

where we have used the weak-\* convergence results for $\omega^{\tau_j} \xrightarrow{\ast} \omega$ and $d\tilde{\mu}_t^{\tau_j} \xrightarrow{\ast} d\tilde{\mu}_t dt$. This establishes that for each admissible sequence $\{(\omega^{\tau_j}, \tilde{\mu}_k^{\tau_j})\}_{j \in \mathbb{N}}$, the objective function value converges to objective function value of (24).

To study the convergence of value function, we consider two quantities: $\limsup_{\tau_j \to 0} J^{*\tau_j}$ and $\liminf_{\tau_j \to 0} J^{*\tau_j}$. We aim to prove the following inequality

$$ \liminf_{\tau_j \to 0} J^{*\tau_j} \geq J^* \geq \limsup_{\tau_j \to 0} J^{*\tau_j}. $$

*Proof of (liminf$_{\tau_j \to 0}$ $J^{*\tau_j} \geq J_o^*$):* Let us consider a sequence of problems, as defined in (38), indexed by $\tau_j$ as $\tau_j \to 0$, and let $\omega^{\tau_j}$ be the sequence of optimal stochastic kernels for each of the problems. As the sequence of $\omega^{\tau_j} \in \mathcal{P}(U)$ is bounded, the Banach Alaoglu theorem ensures the existence of a subsequence (without relabelling) which converges to some $\omega \in \mathcal{P}(U)$, i.e., $\omega^{\tau_j} \to \omega$. Furthermore from Theorem 5 we saw that the interpolated trajectories $\tilde{\mu}_t^{\tau_j}$, associated with each $\tau_j$ in the convergent subsequence, converge to a feasible trajectory of (24). Thus,

$$ \begin{aligned} \liminf_{\tau_j \to 0} J^{*\tau_j} &:= \liminf_{\tau_j \to 0} \sum_{k=0}^{N} \tau_j \int_{S(t) \times U} \mathcal{L}(x, u) d\omega^{\tau_j}(u|k\tau_j, x) d\tilde{\mu}_k^{\tau_j}(x) \\ &= \int_0^T \int_{S(t) \times U} \mathcal{L}(x, u) d\omega(u|t, x) d\tilde{\mu}_t(x) dt \geq J_o^* \end{aligned} $$

where for the last inequality we have used the fact that $d\mu = d\eta(\zeta|x, t)d\omega(u|t, x)d\tilde{\mu}_t(x)dt$ is a feasible solution to (24). So, we obtain the desired inequality $\liminf_{\tau_j \to 0} J^{*\tau_j} \geq J_o^*$.

*Proof of (limsup$_{\tau_j \to 0}$ $J^{*\tau_j} \leq J_o^*$):* Let $\mu$ be an optimal solution of (24) which can be decomposed as

$$ d\mu(t, x, \zeta, u) = d\eta(\zeta|x, t)d\omega(u|t, x)d\tilde{\mu}_t(x)dt. \quad (60) $$

Now, let us consider sequence of time steps $\tau_j \to 0$ and consider partitions of $[0, T]$ for these time steps. We then discretize the optimal feedback $\omega$ in (60), with time steps $\tau_j$, to obtain a sequence of stochastic kernels $\{\omega^{\tau_j}\}_{j \in \mathbb{N}}$. We have

$$ \begin{aligned} J_o^* &= \int_0^T \int_{S(t) \times U} \mathcal{L}(x, u) d\omega(u|t, x) d\mu_t(x) dt \\ &= \lim_{\tau_j \to 0} \sum_{k=0}^N \tau_j \int_{S_k \times U} \mathcal{L}(x, u) d\omega^{\tau_j}(u|k\tau_j, x) d\tilde{\mu}_k^{\tau_j}(x). \end{aligned} $$

We note that $\omega^{\tau_j}, \tilde{\mu}^{\tau_j}$ are admissible for (38) with time step $\tau_j$. So,

$$ \sum_{k=0}^{N} \tau_j \int_{S_k \times U} \mathcal{L}(x, u) d\omega^{\tau_j}(u|k\tau, x) d\tilde{\mu}_k^{\tau_j}(x) \geq J^{\tau*}. $$

Thus we have the desired inequality $\limsup_{\tau_j \to 0} J^{*\tau_j} \leq J_o^*$.

So, $\lim_{\tau_j \to 0} J^{*\tau_j}$ exists, $\lim_{\tau_j \to 0} J^{*\tau_j} = J_o^*$ for each initial condition $x_0 \in S(0)$. $\square$

<figure><img src="image_5.png" /></figure>

**Corollary 1.** Let $\mu_0 \in \mathcal{P}(S(0))$. Then, as $\tau_j \to 0$,

$$ \lim_{\tau_j \to 0} \int_{S(0)} |J^{*\tau_j}(x) - J_o^*(x)| d\mu_0(x) = 0, $$

with convergence in the $L^1$-norm w.r.t. $\mu_0$.
V. COMPUTATIONAL ASPECTS

To provide a numerical approximation of the optimal control problem (38) studied in this paper, which is the measure relaxation of the discrete-time optimal control problem (5), we consider a finite-dimensional relaxation of the problem. This relaxation is obtained by reformulating the infinite-dimensional problem over measures as a finite-dimensional problem over their moments, truncated to a certain degree. The resulting truncated moment problem can then be relaxed into a family of semidefinite programs (SDP), which can be solved using some off-the-shelf solvers.

Let $\mathbb{R}[x, u]$ be the ring of real polynomials in the variables $x \in \mathbb{R}^n$, $u \in \mathbb{R}^m$ where $x$ and $u$ are the state and the control vectors, respectively. Let $\mathbb{R}[x, u]_d \subset \mathbb{R}[x, u]$ be the subset of polynomials of degree up to $d$. We make the following assumptions:

**Assumption 4.** Each set $S_k \subset \mathbb{R}^n$ in equation (39) is a basic semialgebraic set, i.e.,

$$S_k := \{x \in \mathbb{R}^n \mid g_i^{(k)}(x) \ge 0, i = 1, \dots, \ell\},$$

where each $g_i^{(k)}$ is a polynomial. For notational convenience, let $g_0^{(k)}(x) = 1$, and that one of the inequalities is of the form $R - \|x\|^2 \ge 0$ for some sufficiently large constant $R > 0$. Similarly, the control input set $U \subset \mathbb{R}^m$ is a basic semialgebraic set defined as

$$U := \{u \in \mathbb{R}^m \mid h_j(u) \ge 0, j = 1, \dots, r\},$$

where each $h_j$ is a polynomial function, $h_0(u) = 1$, and one of the constraints is $R - \|u\|^2 \ge 0$. Furthermore, the cost function $\mathcal{L}(x, u)$ in (38) is a polynomial function.

For integers $n, m \ge 1$ and a fixed degree $d \in \mathbb{N}$, let

$$\mathcal{A}_d := \{\alpha = (\alpha_1, \alpha_2) \in \mathbb{N}^{n+m} : |\alpha_1| + |\alpha_2| \le d\}.$$

Define the vector of monomial basis functions

$$b_d(x, u) := [x^{\alpha_1} u^{\alpha_2}]_{\alpha \in \mathcal{A}_d} \in \mathbb{R}^{n_d}$$

where $n_d := |\mathcal{A}_d| = \binom{n+m+d}{d}$. Given a moment sequence $y = \{y_\alpha\}_{\alpha \in \mathbb{N}^{n+m}}$, define the Riesz functional

$$\ell_y : \mathbb{R}[x, u]_d \to \mathbb{R}, \quad \ell_y \left( \sum_{\alpha \in \mathcal{A}_d} f_\alpha x^{\alpha_1} u^{\alpha_2} \right) = \sum_{\alpha \in \mathcal{A}_d} f_\alpha y_\alpha.$$

The truncated moment matrix of order $d$ associated with $y$ is defined by

$$M_d(y) := \ell_y(b_d(x, u) b_d(x, u)^\top) = [y_{\alpha+\beta}]_{\alpha, \beta \in \mathcal{A}_d}$$

where $\alpha+\beta := (\alpha_1+\beta_1, \alpha_2+\beta_2)$. Given a non-negative Borel measure $\nu \in \mathcal{M}_+(\mathbb{R}^{n+m})$, we can naturally associate its moment sequence $y$ as $y = \int x^{\alpha_1} u^{\alpha_2} d\nu(x, u)$; the corresponding Riesz functional $\ell_y$ represents integration against $\nu$.

Now let us consider the program in (38)-(39) and let $y_k$ be the moment sequence associated to $\mu_k^\tau$ at time step $k$. Then, the objective function in (38) can be written as

$$\sum_{k=1}^{N} \int \mathcal{L}(x, u) d\mu_k^\tau(x, u) = \sum_k c_l^T y_k \quad (61)$$

where $c_l$ is the coefficient vector for the cost function expressed in monomial basis. For the dynamics constraint in program (38) and any test function $v(x)$, we observe that

$$\begin{aligned} \int_{S_k} v(x) d\mu_k^\tau(x, u) &= \int_{S_k} v(x) d\omega^\tau(du|x) d\tilde{\mu}_k^\tau(x) \\ &= \int_{S_{k-1}} \int_U v(x) dG_{k\#}^\tau \mu_{k-1}^\tau(x, u) \\ &= \int_{S_{k-1}} \int_U v(G_k^\tau(x, u)) d\mu_{k-1}^\tau(x, u) \end{aligned} \quad (62)$$

where we have used (35) and (37) in the first and second equality. When $v(x)$ are monomials, we obtain

$$y_k - c_f^T y_{k-1} = 0 \quad (63)$$

where $c_f$ is coefficient of $v(G_k^\tau(x, u))$ in the RHS of (62). For every $h \in \mathbb{R}[x, u]_d$ such that $h$ is the vector of coefficients of $h$ we have

$$\begin{aligned} \ell_{y_k}(h^2) &= \int \mathbf{h} b_d(x, u) b_d(x, u)^\top \mathbf{h}^\top d\mu_k^\tau(x, u) \\ &= \mathbf{h}^\top M_d(y_k) \mathbf{h} = \int h^2(x, u) d\mu_k^\tau(x, u) \ge 0. \end{aligned}$$

Hence, for arbitrary $h(x)$ we have $\mathbf{h}^\top M_d(y_k) \mathbf{h} \ge 0$, and we conclude that

$$M_d(y_k) \ge 0, \quad \text{for each order } d \ge 0, \quad (64)$$

where the notation $A \succeq 0$ denotes that the matrix $A$ is positive semidefinite. The localizing matrix associated with the polynomial $g_i^{(k)}$ and the moment sequence $y_k$ is defined as

$$M_{d-\deg(g_i^{(k)})}(g_i^{(k)} y_k) = \int b(x, u)_{d-\deg(g_i^{(k)})} b(x, u)_{d-\deg(g_i^{(k)})}^\top g_i^{(k)}(x) d\mu_k^\tau(x, u),$$

where $\deg(g_i^{(k)})$ denotes the degree of the polynomial $g_i^{(k)}$ introduced in Assumption 4. To enforce that the support of the measure $\mu_k^\tau$ is contained in $S$, we require the associated localizing matrix to be positive semidefinite:

$$\begin{aligned} M_{d-\deg(g_i^{(k)})}(g_i^{(k)} y_k) \ge 0 \quad &\text{for } d \ge \max_i \deg(g_i^{(k)}) \\ M_{d-\deg(h_j^{(k)})}(h_j^{(k)} y_k) \ge 0 \quad &\text{for } d \ge \max_j \deg(h_j^{(k)}). \end{aligned} \quad (65)$$

Furthermore, the following theorem provides necessary and sufficient condition for the existence of representing measure corresponding to an infinite moment sequence $y_k$.

**Theorem 7.** Under Assumption 4, let $y_k$ be a moment sequence which satisfies (64) and (65), then $y_k$ has a representing measure with support $S_k \times U$.

Note, however, that for any fixed truncation order $d$, a truncated vector $y_k$ that satisfies (64) and the localizing conditions (65) is not guaranteed to admit a representing measure. Consequently, the constraints (64) and (65) are only necessary conditions; imposing them therefore yields a truncated moment relaxation (an outer approximation) of the original problem.
Under Assumption 3, an order d truncation of the program in (38) is given by

$$J_d^* = \min_{\{y_k\}} \sum_{k=1}^{N} c_l^T y_k \quad (66)$$

$$
\begin{align*}
& \text{s.t. } y_k - c_f^T y_{k-1} = 0 \quad \forall k = 1, \dots, N \\
& M_d(y) \ge 0, \quad M_{d-\deg(g_i^{(k)})}(g_i^{(k)}y) \ge 0, \forall i = 1, \dots, \ell \\
& M_{d-\deg(h_j)}(h_jy) \ge 0, \quad \forall j = 1, \dots, r \tag{67}
\end{align*}
$$

where $c_l$ and $c_f$ are defined in (61)-(63), and $y_{\mu_0}$ is the moment of the initial distribution up to degree $d$.

**Theorem 8.** Under Assumption 4, the sequence $J_d^*$ (66) is a monotonically non-decreasing sequence and converges to $J^{\tau^*}(38)$ as $d \to \infty$.

For the proof of Theorem 7, and Theorem 8 and a systematic treatment of such moment–SOS relaxations, see [13], [33]. For rate of convergence of the sequence, see [34].

## VI. EXAMPLE

In this section, we illustrate the SDP relaxations introduced in (66) using an academic example of a two dimensional system. Consider the following example

$$\dot{\mathbf{x}}(t) \in (1 + u_x, u_y) - \mathcal{N}_S(\mathbf{x}(t)) \quad (68)$$

where, $S := \{(x, y) : x^2 + y^2 \le 1\}$ and $(u_x, u_y)$ is the control vector. We have the following discretization of (68):

$$(x_{k+1}, y_{k+1}) = P_S((x_k, y_k) + \tau(1 + u_x, u_y)). \quad (69)$$

To be able to apply the techniques described in the previous subsection and obtain moment relaxations of the form (66) for the problem defined in (5), we need to obtain a polynomial expression for the projection map onto the disc. This can be achieved by first noting that the expression for the projection map onto the set $S$ in (69) is

$$P_S(w) = \begin{cases} w, & \text{if } \|w\| \le 1, \\ \frac{w}{\|w\|}, & \text{if } \|w\| \ge 1. \end{cases}$$

Next, we introduce variable $z = 1/\sqrt{x^2 + y^2}$ to obtain an equivalent definition of projection map:

$$P_S(x, y) = \begin{cases} (x, y), & x^2 + y^2 \le 1, \\ \{(zx, zy) : z^2(x^2+y^2)=1\}, & x^2 + y^2 \ge 1. \end{cases}$$

We decompose the measure $\mu_k^\tau$ in (62) into two measures, $\mu_k^{\tau S}$, supported on the interior of the disc and $\mu_k^{\tau \bar{S}}$, supported on the exterior of the disc. Then,

$$\text{supp}(\mu_k^{\tau S}) = \{(x, y) \mid (x + \tau(u_x + 1))^2 + (y + \tau u_y)^2 \le 1\} \quad (70)$$

$$\text{supp}(\mu_k^{\tau \bar{S}}) = \{(x, y) \mid (x + \tau(u_x + 1))^2 + (y + \tau u_y)^2 \ge 1\}. \quad (71)$$

In this case, the dynamics constraint in (39) becomes

$$\mu_k^{\tau S} + \mu_k^{\tau \bar{S}} = \left[ P_S \circ (\tau f(\cdot, \cdot) + \text{id}(\cdot)) \right]_{\#} \mu_{k-1}^{\tau \bar{S}} + (\tau f(\cdot, \cdot) + \text{id}(\cdot))_{\#} \mu_{k-1}^{\tau S}. \quad (72)$$

We consider cost function as $\mathcal{L}(x, u) = (x - x_T)^2 + |u|^2$, the initial distribution is $\mu_0 = \delta_{(0.7,0.5)}$ and target distribution $\mu_T = \delta_{(0.7,0.1)}$. For this cost and dynamics (72), we obtain a SDP relaxation of (38) as described in (66). The convergence result of the values of the associated relaxations follow from Theorem 8.

In this example, we consider the case where the initial and the desired final distributions are Dirac measures. Consequently, the distribution is expected to evolve as a Dirac measure over time. In this setting, the first-order approximate moments $y_k$, also called pseudo-moments, of the distribution $\mu_k^\tau$ serve as an approximation of the first-order moment of the measure at each discrete time step $k$. These moments provide directly information about the state and control trajectories across different time instants, as $\int x \mu_k^\tau(x, u) = \int x \delta_{(x_k^\tau, u_k^\tau)} = x_k^\tau$, $\int u \mu_k^\tau(x, u) = \int u \delta_{(x_k^\tau, u_k^\tau)} = u_k^\tau$.

For more general initial conditions, the extraction of state and control trajectories boil down to evaluating the support of the measure. This can be achieved using Christoffel-Darboux kernels, as discussed in [35], or [26] for a more specific exposition. Specifically, the method proposed in [35] employs a semialgebraic approximant based on the Christoffel polynomial, which is constructed using the moment information of the measure at each degree $d$. For the uncontrolled case, [27] proposed SDP relaxation to simulate the dynamical system (68) in the continuous time. And similar approach can be followed for the SDP relaxation optimal control problem in the continuous time as well.

The simulation results are displayed in Figure 3, where we plot the sum of first order moments of $\mu_k^{\tau S}$ and $\mu_k^{\tau \bar{S}}$ over the time horizon of the problem. The results are computed for a fixed relaxation order using GloptiPoly and Mosek. We observe that the program finds an optimal solution which prefers to slide along the boundary of the set $S$ before the control kicks the particle back into the set in order to reach the target configuration. Moreover, in Figure 4 we observe that during the sliding motion zero controls are applied.

## VII. CONCLUSION

We addressed the problem of optimal control for nonsmooth dynamical systems by formulating it as a linear problem in the cone of non-negative Borel measures. We showed that relaxing the problem from finite dimensional space to infinite dimensional does not produce any relaxation gap. We further showed the convergence of the discrete-time problem to continuous-time problem in the space of measures using tools from optimal transport theory. We defined interpolated curves and then showed that as the time step goes to zero, the interpolated curves converge to the solutions of the continuity equation where the corresponding velocity vector defines the nonsmooth dynamical system. We further establish the convergence of the value function of the discrete-time optimal control problem to the value function of continuous-time problem as the time step goes to zero. We observe that in comparison to the continuous-time problem, the discrete-time case has fewer variables as there is no extra variable corresponding to the selection of the vector field from the normal cone, which makes it computationally less demanding.
<figure><img src="image_6.png" /><figcaption>Fig. 3: Approximate first order moments solving linear problem (66) for nonsmooth dynamical system (68).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4: Approximate control in x and y directions, where ux/uy interior denote controls from the interior measure $\mu_k^{\tau S}$ (70) and ux/uy exterior from the exterior measure $\mu_k^{\tau S}$ (71).</figcaption></figure>

## ACKNOWLEDGEMENT

S. P. Chhatoi would like to thank Jérôme Bertrand, Elsa Cazelles and Simone DiMarino for discussions related to optimal transport theory.

## REFERENCES

[1] B. Brogliato and A. Tanwani. Dynamical systems coupled with mono-tone set-valued operators: Formalisms, applications, well-posedness, and stability. *SIAM Review*, 62(1):3–129, 2020.

[2] M. Athans and P.L. Falb. *Optimal Control: An Introduction to the Theory and Its Applications*. Dover Books on Engineering. Dover Publications, 2007.

[3] D. Liberzon. *Calculus of Variations and Optimal Control Theory: A Concise Introduction*. 2012.

[4] R. Vinter. *Optimal Control*. Birkhäuser Boston, 2010.

[5] M. Brokate and P. Krejčí. Optimal control of ODE systems involving a rate independent variational inequality. *Discrete and Continuous Dynamical Systems-series B*, 18:331–348, 2012.

[6] G. Colombo, R. Henrion, D.H. Nguyen, and B.S. Mordukhovich. Optimal control of the sweeping process over polyhedral controlled sets. *Journal of Differential Equations*, 260(4):3397–3447, 2016.

[7] C. Hermosilla and M. Palladino. Optimal control of the sweeping process with a nonsmooth moving set. *SIAM Journal on Control and Optimization*, 60(5):2811–2834, 2022.

[8] P. Bettiol, G. Colombo, and P. Gidoni. Penalization and necessary optimality conditions for a class of nonsmooth sweeping processes. *Set-Valued and Variational Analysis*, 33(1):33, 2025.

[9] D. Stewart and M. Anitescu. Optimal control of systems with discontinuous differential equations. *Numerische Mathematik*, 114:653–695, 2012.

[10] A. Nurkanović and M. Diehl. NOSNOC: A software package for numerical optimal control of nonsmooth systems. *IEEE Control Systems Letters*, 6:3110–3115, 2022.

[11] A. Nurkanović, A. Pozharskiy, J. Frey, and M. Diehl. Finite elements with switch detection for numerical optimal control of nonsmooth dynamical systems with set-valued heaviside step functions. *Nonlinear Analysis: Hybrid Systems*, 54:101518, 2024.

[12] A. Vieira, B. Brogliato, and C. Prieur. Quadratic optimal control of linear complementarity systems: First-order necessary conditions and numerical analysis. *IEEE Transactions on Automatic Control*, 65(6):2743–2750, 2020.

[13] D. Henrion, M. Korda, and J.B. Lasserre. *The Moment-SOS Hierarchy*. Series On Optimization And Its Applications. World Scientific Publishing Company, 2020.

[14] R. Vinter. Convex duality and nonlinear optimal control. *SIAM Journal on Control and Optimization*, 31(2):518–538, 1993.

[15] J.-B. Lasserre, D. Henrion, C. Prieur, and E. Trélat. Nonlinear optimal control via occupation measures and LMI-relaxations. *SIAM Journal on Control and Optimization*, 47(4):1643–1666, 2008.

[16] O. Hernández-Lerma and J.-B. Lasserre. *Discrete-Time Markov Control Processes: Basic Optimality Criteria*. Stochastic Modelling and Applied Probability. Springer New York, 2012.

[17] A. G. Bhatt and V. S. Borkai. Occupation measures for controlled Markov processes: Characterization and optimality. *Annals of Probability*, 24:1531–1562, 1996.

[18] A. Agrachev and P. Lee. Optimal transportation under nonholonomic constraints. *Transactions of the American Mathematical Society*, 361:6019–6047, 2007.

[19] A. Hindawi, L. Rifford, and J.B. Pomet. Mass transportation with LQ cost functions. *Acta Applicandae Mathematicae*, 113:215–229, 2011.

[20] Y. Chen, T. T. Georgiou, and M. Pavon. Optimal transport over a linear dynamical system. *IEEE Transactions on Automatic Control*, 62(5):2137–2152, 2017.

[21] S.P. Chhatoi, A. Tanwani, and D. Henrion. Optimal control of nonsmooth dynamical systems using measure relaxations. In *2024 IEEE 63rd Conference on Decision and Control (CDC)*, pages 888–893, 2024.

[22] S.P. Chhatoi, A. Tanwani, and D. Henrion. Optimal control of state constrained systems via measure relaxations and polynomial optimization. Available online: https://hal.science/hal-05345286, https://spykspeigel.github.io/OCNSD/OCNSd.pdf, 2025.

[23] P. Billingsley. *Probability and Measure*. Wiley, New York, 3rd edition, 1995.

[24] H. Brézis. *Functional Analysis, Sobolev Spaces and Partial Differential Equations*. Springer, New York, 2011.

[25] P. Pedregal. *Parametrized Measures and Variational Principles*. Birkhäuser Architectural Guide. Springer Nature, 1997.

[26] M. Souaiby, A. Tanwani, and D. Henrion. Ensemble approximations for constrained dynamical systems using Liouville equation. *Automatica*, 149:110836, 2023.

[27] S.P. Chhatoi, A. Tanwani, and D. Henrion. Evolution of measures in nonsmooth dynamical systems: Formalisms and computation. *Automatica*, 179(112402), 2025.

[28] J.-P. Aubin and H. Frankowska. *Set-Valued Analysis*. Systems & Control: Foundations & Applications. Birkhäuser, Boston, 1990.

[29] A.I. Barvinok. A course in convexity. In *Graduate Studies in Mathematics*, 2002.

[30] V. I. Bogachev. *Measure Theory*. Springer Berlin Heidelberg, 2007.

[31] L. Ambrosio, N. Gigli, and G. Savaré. *Gradient Flows: In Metric Spaces and in the Space of Probability Measures*. Birkhäuser, Basel, 2005.

[32] F. Santambrogio. *Optimal Transport for Applied Mathematicians. Calculus of Variations, PDEs and Modeling*. Birkhäuser Cham, 2015.

[33] J.-B. Lasserre. *Moments, Positive Polynomials and Their Applications*, volume 1 of *Series on Optimization and Its Applications*. Imperial College Press, London, 2009.

[34] C. Schlosser, M. Tacchi-Bénard, and A. Lazarev. Convergence rates for the moment-SoS hierarchy. *Numerical Algebra, Control and Optimization*, 2025.

[35] J.-B. Lasserre, E. Pauwels, and M. Putinar. *The Christoffel-Darboux Kernel for Data Analysis*. Cambridge University Press, 2022.

# APPENDIX A

## PROOF OF LEMMA 1

*Proof.* By Markov's inequality and the coercivity constraint
$$ \int_{\mathcal{B} \times U} |\zeta| d\mu \le C_\eta, $$

$$ \sup_{(\mu, \mu_T) \in \mathcal{D}} \mu(\{|\zeta| > R\}) \le \frac{C_\eta}{R} \xrightarrow{R \to \infty} 0, $$
so {$\mu$} is tight in the $\zeta$-coordinate. Since $[0, T] \times \Omega \times U$ is compact, this implies tightness of {$\mu$} on $\mathfrak{B} \times U$. Likewise, the compactness of $S(T)$ implies tightness of {$\mu_T$}. With the uniform mass bounds $||\mu|| \le T$ and $||\mu_T|| \le 1$, Prokhorov's theorem [23, Thm. 5.1] yields relative weak-$\star$ sequential compactness of $\mathcal{D}$.

Now we show weak-$\star$ closedness. Let $(\mu_k, \mu_{T,k}) \in \mathcal{D}$ and $(\mu_k, \mu_{T,k}) \stackrel{*}{\to} (\mu, \mu_T)$ (along a subsequence given by the relative weak-$\star$ compactness of $\mathcal{D}$). Then: (i) positivity and the mass bounds are preserved under weak-$\star$ limits; (ii) the map $\nu \mapsto \int_{\mathfrak{B} \times U} |\zeta| d\nu$ is weak-$\star$ lower semicontinuous, hence $\int |\zeta| d\mu \le C_\eta$; (iii) the constraint (22) is stable under weak-$\star$ limits. Indeed, for $\varphi \in C^1([0, T] \times \Omega)$ set $g_1(t, x, \zeta, u) := \partial_t \varphi(t, x)$, $g_2(t, x, \zeta, u) := \nabla_x \varphi(t, x) \cdot (-\zeta + f(x, u))$. The terms with $g_1$ and with $\nabla_x \varphi \cdot f$ are bounded continuous on $\mathfrak{B} \times U$ and therefore converge by narrow convergence. For the $-\nabla_x \varphi \cdot \zeta$ term, fix $R > 0$ and write

$$ \int \nabla_x \varphi \cdot (-\zeta) \, d\mu_k = \int \nabla_x \varphi \cdot (-\zeta) \mathbf{1}_{\{|\zeta| \le R\}} \, d\mu_k + \int \nabla_x \varphi \cdot (-\zeta) \mathbf{1}_{\{|\zeta| > R\}} \, d\mu_k. $$

The first part converges to the corresponding expression with $\mu$ because the truncated integrand is bounded and continuous; the second part is uniformly small, since

$$ \sup_k \left| \int \nabla_x \varphi \cdot (-\zeta) \mathbf{1}_{\{|\zeta| > R\}} \, d\mu_k \right| \le \|\nabla_x \varphi\|_{\infty} \sup_k \int_{|\zeta| > R} |\zeta| \, d\mu_k \le \|\nabla_x \varphi\|_{\infty} \frac{C_\eta}{R}. $$

Let $k \to \infty$ and then $R \to \infty$ to conclude. The terminal term $\int_{S(T)} \varphi(T, x) d\mu_{T,k}$ also converges to $\int_{S(T)} \varphi(T, x) d\mu_T$. Hence (22) holds for $(\mu, \mu_T)$.

Therefore $\mathcal{D}$ is weak-$\star$ closed; combined with relative compactness, this proves weak-$\star$ sequential compactness. $\square$

## APPENDIX B
### PROOF OF LEMMA 2

*Proof.* The occupation measure $\mu$ satisfying (22) admits the following disintegration:

$$ d\mu(t, x, \zeta, u) = d\nu(\zeta, u | t, x) d\bar{\mu}_t(x) dt, \quad (73) $$

where $\nu(\cdot | t, x) \in \mathcal{P}(\mathcal{N}_{S(t)}(x) \times U)$ and $\bar{\mu}_t$ is the marginal of $\mu$ at time $t$. Substituting the disintegration of $\mu$ from (73) in (22), we obtain

$$ \int_{[0,T] \times S(t)} \left[ \partial_t \phi(t, x) d\bar{\mu}_t(x) dt + \partial_x \phi(t, x) \cdot b(t, x) d\bar{\mu}_t(x) dt \right] \\ = \int_{S(T)} \phi(T, x) d\mu_T(x) - \phi(0, x_0). \quad (74) $$

where the averaged vector field $b(t, x)$ is

$$ b(t, x) := \int_{\mathcal{N}_{S(t)}(x) \times U} (-\zeta + f(x, u)) d\nu(\zeta, u | t, x). \quad (75) $$

Let $\pi_{\zeta\#\nu} = \eta$, $\pi_{u\#\nu} = \omega$ (i.e. the marginals of $\nu \in \mathcal{P}(\mathcal{N}_{S(t)}(x) \times U)$), then the vector field in (75)

can be written as, $b(t, x) = -\int_{\mathcal{N}_{S(t)} \times U} \zeta d\nu(\zeta, u|t, x) + \int_{\mathcal{N}_{S(t)} \times U} f(x, u) d\nu(\zeta, u|t, x) = -\int_{\mathcal{N}_{S(t)}} \zeta d\eta(\zeta|t, x) + \int_U f(x, u) d\omega(u|t, x)$. Using Assumption 1, we have $\int_{[0,T] \times S(t)} \int_{\mathcal{N}_{S(t)}} \int_U |f(x, u)| d\nu(\zeta, u|t, x) d\bar{\mu}_t(x) dt < \infty$. The uniform moment bound in $\zeta$ direction (which is part of the definition in set (23)) yields $\int_{[0,T] \times S(t)} \int_{\mathcal{N}_{S(t)}} \int_U |\zeta| d\nu(\zeta, u|t, x) d\bar{\mu}_t(x) dt < \infty$. Using triangle's inequality, we have $\int_0^T \int_{S(t)} |b(t, x)| d\bar{\mu}_t(x) dt < \infty$. So, the integrability of $b(t, x)$ allows the use of superposition principle [31, Theorem 8.2.1]; there exists a probability measure $\Pi \in \mathcal{P}(\mathcal{AC}([0, T]; \mathbb{R}^n))$ such that for each $t$, the marginal $\bar{\mu}_t = (e_t)_{\#}\Pi$, and $\Pi$-almost every path $x(\cdot)$ satisfies the ODE $\dot{x}(t) = b(t, x(t))$. We note that $\dot{x}(t) = b(t, x(t))$ is the same as the Young measure relaxed dynamics in (17) but driven by state dependent $\omega$. For each triple $(x(\cdot), \eta(\cdot|t, x), \omega(\cdot|t, x))$ satisfying (17), define the pair $(\hat{\mu}, \xi)$ as

$$ d\hat{\mu}(t, x, \zeta, u) = dt \delta_{x(t)}(dx) \eta(d\zeta|t, x(t)) \omega(du|t, x(t)), \quad \xi := \delta_{x(T)}. $$

Then $(\hat{\mu}, \xi) \in \mathfrak{R}$. Also, define the push-forward measure $\theta := \Phi_{\#}\Pi$, where $\Phi(x(\cdot)) := (\hat{\mu}, \xi)$. For any admissible cost function $\mathcal{L}(x, u)$, we compute:

$$ \begin{aligned} \int \mathcal{L}(x, u) d\mu &= \int \mathcal{L}(x, u) d\bar{\mu}_t(x) d\eta(\zeta|t, x) d\omega(u|t, x) dt \\ &= \int \mathcal{L}(x, u) \delta_{x(t)}(\Pi(x(\cdot)) d\eta(\zeta|x) d\omega(u|t, x) dt \\ &= \int_{\mathfrak{R}} \left[ \int \mathcal{L}(x, u) d\hat{\mu} \right] d\tilde{\theta}(\hat{\mu}), \end{aligned} $$

where $\tilde{\theta}$ is the marginal of $\theta$. $\square$

## APPENDIX C
### PROOF OF LEMMA 6

*Proof.* We already have shown that $\hat{\mu}_t^\tau \to \mu_t$ uniformly in the $W_2$ metric when $\tau \to 0$. The term

$$ (y - P_{S(t)}(x)) g^{\tau}(t, x) \xrightarrow{\tau \to 0} (y - x) g(t, x) $$

converges pointwise, and the integrand $(y - P_{S(t)}(x)) g^{\tau}(t, x)$ is uniformly bounded (because $y - P_{S(t)}(x)$ lies in a compact set, while $g^{\tau}$ is bounded by Assumption 1). Hence, by the dominated convergence theorem we get

$$ \left| \lim_{\tau \to 0} \int h(t, x)(y - P_{S(t)}(x))g^{\tau}(t, x)d\hat{\mu}_{t}^{\tau}(x) - \int h(t, x)(y - x)g(t, x)d\mu_{t}(x) \right| \to 0. \quad \square $$

## APPENDIX D
### PROOF OF COROLLARY 1

*Proof.* We consider an initial condition $\mu_0 \in \mathcal{P}(S(0))$ in (24) which are different from Dirac distribution $\delta_{x_0}$ just considered above. To prove $L^1$ convergence of $J^{*\tau}[\mu_0]$ as $\tau \to 0$ for notational convenience, we write $\tau$ instead of
the sequence $\tau_j$), we first establish a uniform bound on $J^{*\tau} = \tau \sum_{k=0}^{N} \mathcal{L}(x_k^\tau, u_k^\tau)$. Using Assumption 3, we know that $\mathcal{L}(x_k^\tau, u_k^\tau) \le l(x_k^\tau)(1 + |u_k^\tau|^2)$ for $l(\cdot) \in L_{\text{loc}}^\infty(\mathbb{R}^n)$. Next, we recall from [1] that the following bound holds for $x_k^\tau$ in the discrete-time dynamics (31):

$$|x_k^\tau| \le C_1|x_0| + C_2,$$

where $C_1$ and $C_2$ depend on the system data. Consequently, we have

$$ \begin{aligned} J^{*\tau}(x_0) &= \tau \sum_{k=0}^{N} \mathcal{L}(x_k^\tau, u_k^\tau) \le \tau \sum_{k=0}^{N} l(x_k^\tau)(1 + |u_k^\tau|^2) \\ &\le \tau \sum_{k=0}^{N} l_b(C_1|x_0| + C_2)U_b \le Th_b(C_1\max_{x \in S(0)}|x| + C_2)U_b \end{aligned} $$

where, $l_b$ is the bound on $l(\cdot)$ and $U_b$ is the bound on the controls. Thus we have uniform bounds on the $J_o^*[\mu_0]$ for $\mu_0 \in \mathcal{P}(S(0))$. The uniform bounds and the pointwise convergence of $J_o^*(x_0)$, established in the previous part of the proof, allows us to use the dominated convergence theorem to obtain

$$ \lim_{\tau \to 0} \int |J^{*\tau}(x) - J_o^*(x)| d\mu_0(x) = 0. \quad \square $$