# Minimal-time nonlinear control via semi-infinite programming

Antoine Oustry, Matteo Tacchi

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

## To cite this version:

Antoine Oustry, Matteo Tacchi. Minimal-time nonlinear control via semi-infinite programming. 2023. [hal-04145402v2](https://hal.science/hal-04145402v2)

HAL Id: hal-04145402

https://hal.science/hal-04145402v2

Preprint submitted on 6 Jul 2023

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_2.png" /></figure>
# Minimal-time nonlinear control via semi-infinite programming

Antoine Oustry, Matteo Tacchi

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

## To cite this version:

Antoine Oustry, Matteo Tacchi. Minimal-time nonlinear control via semi-infinite programming. 2023. [hal-04145402v3](https://hal.science/hal-04145402v3)

HAL Id: hal-04145402

https://hal.science/hal-04145402v3

Preprint submitted on 2 Dec 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
MINIMAL TIME NONLINEAR CONTROL VIA SEMI-INFINITE PROGRAMMING

ANTOINE OUSTRY<fnref n="1" />,<fnref n="2" />, MATTEO TACCHI<fnref n="3" />,<fnref n="4" />

**Abstract.** We address the problem of computing a control for a time-dependent nonlinear system to reach a target set in a minimal time. To solve this minimal time control problem, we introduce a hierarchy of linear semi-infinite programs, the values of which converge to the value of the control problem. These semi-infinite programs are increasing restrictions of the dual of the nonlinear control problem, which is a maximization problem over the subsolutions of the Hamilton-Jacobi-Bellman (HJB) equation. Our approach is compatible with generic dynamical systems and state constraints. Specifically, we use an oracle that, for a given differentiable function, returns a point at which the function violates the HJB inequality. We solve the semi-infinite programs using a classical convex optimization algorithm with a convergence rate of $O(1/k)$, where $k$ is the number of calls to the oracle. This algorithm yields subsolutions of the HJB equation that approximate the value function and provide a lower bound on the optimal time. We study the closed-loop control built on the obtained approximate value functions, and we give theoretical guarantees on its performance depending on the approximation error for the value function. We show promising numerical results for three non-polynomial systems with up to 6 state variables and 5 control variables.

**Key words.** Nonlinear control, Minimal time control, Weak formulation, Semi-infinite programming.

## 1. Introduction.

### 1.1. Motivation and related works.
This paper deals with the control of a deterministic dynamical system to reach a target set in a minimal time. We consider a general case of a time-dependent nonlinear system under nonlinear state constraints. Several applications in various fields, such as robotics [22], aerospace [38], maritime routing [30] or medicine [42], can be formulated as minimal time control problems. Minimal time control, also known as time optimal control, can be seen as a special case of the general framework of Optimal Control Problems (OCP). Solving an OCP for such generic dynamics and constraints is a difficult challenge, although deep theoretical tools are available such as the Pontryagin Maximum Principle (PMP) [8, 13, 35] and the Hamilton-Jacobi-Bellman (HJB) equation [14, 15]. Those theoretical tools, initially developed in the unconstrained setting, have been extended to the case of state constraints [10, 37]. From a numerical point of view, the *multiple shooting* techniques [34, 41] are based on the PMP, and reduce to the solution of a two-point boundary value problem. The *direct methods* reduce to the solution a nonlinear programming problem after discretizing the time space, or parameterizing the control $u(t)$ in a finite dimensional subspace [38, 41]. The celebrated Model Predictive Control (MPC) approach belongs to the category of direct methods [9]. Another approach is to compute the value function of the problem as a maximal subsolution of the HJB equation [40]. This approach is related to the weak formulation of the OCP, which is an infinite dimensional linear program (LP) involving occupation measures. The dual problem of this LP is exactly the problem of finding a maximal subsolution of the HJB equation [20, 26]. In [19, 26], the Moment Sum-of-Squares (SoS) hierarchy is used to approximate the solution of the resulting infinite dimensional LPs, in the case where the dynamics and the constraints of the OCP are defined by polynomials. The convergence rate of this numerical scheme is studied in [25] for infinite-time discounted polynomial control problems. Still in the context of polynomial control problems, a work [23] based on the dual LP and the SoS hierarchy also studies the design of a closed-loop controller based on the approximate value function that is computed. In [5], an extension of the SoS hierarchy based on Kernel methods is employed to extend this computation to general nonlinear system. Regarding the methods specifically dedicated to time-optimal control, we find the same categories: direct methods such as MPC [39], indirect methods based on the PMP and the bang-bang property [28, 31] or methods based on convex optimization [27].

### 1.2. Contribution.
In this paper, we focus on the problem of computing a control to reach a target set in a minimal time. We follow the line of works that use convex optimization to solve the dual problem of the nonlinear control problem, over the subsolutions of the HJB equation [19, 20, 25, 26, 40]. In contrast to several works using the Moment-SoS hierarchy [23, 25, 26, 33, 36], the dynamical system and the state constraints considered here are generic and, in particular, are not assumed to be defined by polynomials. Instead of using polynomial optimization theory and the associated positivity certificates, our approach relies on the existence of a *separation oracle* capable of returning, for a given differentiable function $V$, a point $(t, x)$ where the function

<footnote n="1">¹Ecole des Ponts, Marne-la-Vallée, France.</footnote>

<footnote n="2">²Laboratoire d'informatique de l'École polytechnique, Institut Polytechnique de Paris, Palaiseau, France.</footnote>

<footnote n="3">³Univ. Grenoble Alpes, CNRS, Grenoble INP (Institute of Engineering Univ. Grenoble Alpes), GIPSA-lab, 38000 Grenoble, France.</footnote>

<footnote n="4">⁴Corresponding author. matteo.tacchi@gipsa-lab.fr</footnote>
MINIMAL TIME NONLINEAR CONTROL VIA SEMI-INFINITE PROGRAMMING

ANTOINE OUSTRY<fnref n="1" />,<fnref n="2" />,<fnref n="3" />, MATTEO TACCHI<fnref n="4" />,<fnref n="5" />

**Abstract.** We address the problem of computing a control for a time-dependent nonlinear system to reach a target set in a minimal time. To solve this minimal time control problem, we introduce a hierarchy of linear semi-infinite programs, the values of which converge to the value of the control problem. These semi-infinite programs are increasing restrictions of the dual of the nonlinear control problem, which is a maximization problem over the subsolutions of the Hamilton-Jacobi-Bellman (HJB) equation. Our approach is compatible with Lipschitz dynamical systems and state constraints. Specifically, we use an oracle that, for a given differentiable function, returns a point at which the function violates the HJB inequality. We solve the semi-infinite programs using a classical convex optimization algorithm with a convergence rate of $O(1/k)$, where $k$ is the number of calls to the oracle. This algorithm yields subsolutions of the HJB equation that approximate the value function and provide a lower bound on the optimal time. We study the closed-loop control built on the obtained approximate value functions, and we give theoretical guarantees on its performance depending on the approximation error for the value function. We show promising numerical results for three non-polynomial systems with up to 6 state variables and 5 control variables.

**Key words.** Nonlinear control, Minimal time control, Weak formulation, Semi-infinite programming.

# 1. Introduction.

## 1.1. Motivation and related works.
This paper deals with the control of a dynamical system to reach a target set in a minimal time. We consider a general case of a time-dependent nonlinear Lipschitz system under nonlinear state constraints. Several applications in various fields, such as robotics [30], aerospace [48], maritime routing [39] or medicine [52], can be formulated as minimal time control problems. Minimal time control, also known as time optimal control, can be seen as a special case of the general framework of Optimal Control Problems (OCP). Solving an OCP for such a general case is a difficult challenge, although deep theoretical tools are available such as the Pontryagin Maximum Principle (PMP) [9, 14, 45] and the Hamilton-Jacobi-Bellman (HJB) equation [15, 17]. Those theoretical tools, initially developed in the unconstrained setting, have been extended to the case of state constraints [11, 47]. From a numerical point of view, the *multiple shooting* techniques [44, 51] are based on the PMP, and reduce to the solution of a two-point boundary value problem. The *direct methods* reduce to the solution of a nonlinear programming problem after discretizing the time space, or parameterizing the control $u(t)$ in a finite dimensional subspace [48, 51]. The celebrated Model Predictive Control (MPC) approach belongs to the category of direct methods [10]. Another approach is to compute the value function of the problem as a maximal subsolution of the HJB equation [50]. This approach is related to the weak formulation of the OCP, which is an infinite dimensional linear program (LP) involving occupation measures [34]. The dual problem of this LP is exactly the problem of finding a maximal subsolution of the HJB equation [27, 34]. In [26, 34], the Moment Sum-of-Squares (SoS) hierarchy is used to approximate the solution of the resulting infinite dimensional LPs, in the case where the dynamics and the constraints of the OCP are defined by polynomials. The convergence rate of this numerical scheme is studied in [33] for infinite-time discounted polynomial control problems. Still in the context of polynomial control problems, a work [31] based on the dual LP and the SoS hierarchy also studies the design of a closed-loop controller based on the approximate value function that is computed. In [5], an extension of the SoS hierarchy based on Kernel methods is employed to extend this computation to general nonlinear system. Regarding the methods specifically dedicated to time-optimal control, we find the same categories: direct methods such as MPC [49], indirect methods based on the PMP and the bang-bang property [36, 40] or methods based on convex optimization [35].

## 1.2. Contribution and paper organization.
In this paper, we focus on the problem of computing a control to reach a target set in a minimal time. We follow the line of works that use convex optimization to solve the dual problem of the nonlinear control problem, over the subsolutions of the HJB equation [26, 27, 33, 34, 50]. In contrast to several works using the Moment-SoS hierarchy [31, 33, 34, 43, 46], the dynamical system and the state constraints considered here are only assumed to be Lipschitz and, in particular, are not assumed to be defined by polynomials. Instead of using polynomial optimization theory and the associated positivity certificates, our approach relies on the existence of a *separation oracle* capable of returning, for a given differentiable function V, a point $(t, x)$ where the function V does not satisfy the HJB inequality. Such

<footnote n="1">¹Ecole des Ponts, Marne-la-Vallée, France.</footnote>

<footnote n="2">²Laboratoire d'informatique de l'École polytechnique, Institut Polytechnique de Paris, Palaiseau, France.</footnote>

<footnote n="3">³Réseau de Transport d'Électricité, La Défense, France.</footnote>

<footnote n="4">⁴Univ. Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab, Grenoble, France.</footnote>

<footnote n="5">⁵Corresponding author. matteo.tacchi@gipsa-lab.fr</footnote>
V does not satisfy the HJB inequality. Such an oracle can be provided by a global optimization solver or by a sampling scheme in a black-box optimization approach. In particular, our approach is compatible with the sampled-data control paradigm [5, 8, 24]. Our contribution is manifold

* We introduce a hierarchy of linear semi-infinite programs, the values of which converge to the value of the control problem. After regularization, we solve these semi-infinite programs using a classical algorithm with a convergence rate in $O(1/k)$, where $k$ is the number of calls to the oracle. This yields subsolutions of the HJB equation that lower-approximate the value function and provide a certified lower bound on the minimum time.

* It is known that one can leverage any function $V(t, x)$ approximating the value function, to design a closed-loop, i.e., feedback controller [18, 23]. In this paper, we study the existence of trajectories generated by such a controller.

* We study the performance of such a closed-loop controller, depending on how well $V(t, x)$ approximates the value function, in a way distinct from the analysis in [23]. In particular, this novel analysis enables us to give a sufficient condition for the closed-loop controller to effectively generate a trajectory reaching the target set within the considered time horizon.

* We perform numerical experiments on three non-polynomial controlled systems and compute lower and upper bounds on the minimum time.

**1.3. Mathematical notation.** For any $p \in \mathbb{N}^*$, and $k \in \mathbb{N} \cup \{\infty\}$, we denote by $C^k(\mathbb{R}^p) = C^k(\mathbb{R}^p, \mathbb{R})$ the vector space of real-valued functions with $k$ continuous derivatives over $\mathbb{R}^p$. For a given set $A \subset \mathbb{R}^p$, for any function $f \in C^k(\mathbb{R}^p)$, we denote by $f|_A$ the restriction of $f$ on $A$; moreover, we define the vector space $C^k(\mathbb{R}^p|A) = \{f|_A: f \in C^k(\mathbb{R}^p)\}$ of restrictions on $A$ of $C^k$ functions. For any locally Lipschitz function $f$, we denote by $\partial^c f$ its Clarke subdifferential [12], to be distinguished from $\partial_{x_i} g$, the partial derivative of a differentiable function $g$ with respect to $x_i$.

For any two Lebesgue integrable functions $f, g \in L^1(\mathbb{R}^p)$, we define the convolution product $f \star g = g \star f$ as $f \star g(x) = \int_{\mathbb{R}^p} f(x)g(x-h)dh$. We emphasize that this convolution product is also well defined if $f$ is supported on a compact set, and $g$ is locally integrable. We denote by $\mathbb{R}[x_1, \dots, x_p]$ the vector space of real multivariate polynomials with variables $x_1, \dots, x_p$, and $\mathbb{R}_d[x_1, \dots, x_p]$ the vector space of such real multivariate polynomials with degree at most $d$.

For any set $A \subset \mathbb{R}^p$, we write $\text{conv}(A)$ for the convex hull of the set $A$. For any nonempty set $A$, and any $x \in \mathbb{R}^p$, we denote $d(x, A) = \inf_{a \in A} \|x - a\|_2$ the distance between the set $A$ and the point $x$. We also define the contingent cone to $A$ at $x \in A$, denoted $T_A(x)$ as the set of directions $d \in \mathbb{R}^p$, such that there exist a sequence $(t_k) \in \mathbb{R}_+^n$, and a sequence $(d_k) \in (\mathbb{R}^p)^n$, satisfying $t_k \to 0$, $d_k \to d$, and $x + t_k d_k \in A$, for all $k \in \mathbb{N}$. Finally, we say that a property $P$ holds “almost everywhere” (a.e.) on $A$, or equivalently “for almost all $x \in A$”, to denote that there exists a set $N$ of Lebesgue measure zero such that the property $P$ holds for all $x \in A \setminus N$.

## 2. Problem statement and Linear Programming formulations.

**2.1. Definition of the minimal time control problem.** Let $n$ and $m$ be nonzero integers. We consider on $\mathbb{R}^n$ the control system

$$ (2.1) \qquad \dot{x}(t) = f(t, x(t), u(t)), $$

where $f: \mathbb{R} \times \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}^n$ is Lipschitz continuous, and where the controls are bounded measurable functions, defined on intervals $[t_0, t_1] \subset [0, T]$, and taking their values in a compact set $U$ of $\mathbb{R}^m$. Let $X$ and $K \subset X$ be compact sets of $\mathbb{R}^n$ and $x_0 \in \mathbb{R}$. For $t_0, t_1 \ge 0$, a control $u$ is said admissible on $[t_0, t_1]$ whenever the solution $x(\cdot)$ of (2.1), such that $x(t_0) = x_0$, is well defined on $[t_0, t_1]$ and satisfies the constraints

$$ (2.2) \qquad (x(t), u(t)) \in X \times U, \quad \text{a. e. on } [t_0, t_1], $$

and satisfies the terminal state constraint

$$ (2.3) \qquad x(t_1) \in K. $$

We denote by $\mathcal{U}(t_0, t_1, x_0)$ the set of admissible controls on $[t_0, t_1]$. We consider the question of the minimal time problem from $x_0$ to $K$,

$$ (2.4) \qquad V^*(t_0, x_0) = \inf_{\substack{t_1 \in [t_0, T] \\ u(\cdot) \in \mathcal{U}(t_0, t_1, x_0)}} t_1 - t_0. $$
an oracle can be provided by a global optimization solver or by a sampling scheme in a black-box optimization approach. In particular, our approach is compatible with the sampled-data control paradigm [5, 9, 32]. Our contribution is manifold

* We introduce a hierarchy of linear semi-infinite programs, the values of which converge to the value of the control problem. After regularization, we solve these semi-infinite programs using a classical algorithm with a convergence rate in $O(1/k)$, where $k$ is the number of calls to the oracle. This yields subsolutions of the HJB equation that lower-approximate the value function and provide a certified lower bound on the minimum time.

* It is known that one can leverage any function $V(t, x)$ approximating the value function, to design a closed-loop, i.e., feedback controller [25, 31]. In this paper, we study the existence of trajectories generated by such a controller.

* We study the performance of such a closed-loop controller, depending on how well $V(t, x)$ approximates the value function, in a way distinct from the analysis in [31]. In particular, this novel analysis enables us to give a sufficient condition for the closed-loop controller to effectively generate a trajectory reaching the target set within the considered time horizon.

* We perform numerical experiments on three non-polynomial controlled systems and compute lower and upper bounds on the minimum time.

The paper is organized as follows. Section 2 is devoted to the statement of the minimal time control problem and the presentation of an infinite dimensional Linear Programming formulations with the same value. Section 3 presents a hierarchy of semi-infinite linear programming problems that approach the infinite dimensional problems and allow to obtain approximate value functions of the control problem. Section 4 presents the feedback controller obtained from approximate value functions and studies the performance of such a controller. Section 5 is devoted to numerical experiments.

### 1.3. Mathematical notation.
$\mathbb{N}$ denotes the set of natural integers and $\mathbb{N}^* = \mathbb{N} \setminus \{0\}$, $\mathbb{R}$ stands for real numbers, $\mathbb{R}_+$ is the set of nonnegative real numbers and $\mathbb{R}_{++} = \mathbb{R}_+ \setminus \{0\}$. For any $p \in \mathbb{N}^*$, and $k \in \mathbb{N} \cup \{\infty\}$, we denote by $C^k(\mathbb{R}^p) = C^k(\mathbb{R}^p, \mathbb{R})$ the vector space of real-valued functions with $k$ continuous derivatives over $\mathbb{R}^p$. For a given set $A \subset \mathbb{R}^p$, for any function $f \in C^k(\mathbb{R}^p)$, we denote by $f|_A$ the restriction of $f$ on $A$; moreover, we define the vector space $C^k(\mathbb{R}^p|A) = \{f|_A : f \in C^k(\mathbb{R}^p)\}$ of restrictions on $A$ of $C^k$ functions. For any locally Lipschitz function $f$, we denote by $\partial^c f$ its Clarke subdifferential [13, Definition 1.1], to be distinguished from $\partial_{x_i} g$ and $\nabla_x V$, the partial derivative of a differentiable function $g$ with respect to $x_i$ and the gradient with respect to the space variable $x$ of a differentiable function $(t, x) \mapsto V(t, x)$, respectively. For any two measurable functions $f, g$, we denote the convolution product by $f \star g = g \star f$ as $f \star g(x) = \int_{\mathbb{R}^p} f(h)g(x - h)dh$, whenever this integral is finite. We denote by $\mathbb{R}[x_1, \dots, x_p]$ the vector space of real multivariate polynomials with variables $x_1, \dots, x_p$, and $\mathbb{R}_d[x_1, \dots, x_p]$ the vector space of such real multivariate polynomials with degree at most $d$. For any polynomial $w \in \mathbb{R}[x_1, \dots, x_p]$, we denote by $\deg(w)$ its degree. For any set $A \subset \mathbb{R}^p$, we write $\text{conv}(A)$ for the convex hull of the set $A$. For any nonempty set $A$, and any $x \in \mathbb{R}^p$, we denote $d(x, A) = \inf_{a \in A} \|x - a\|_2$ the distance between the set $A$ and the point $x$. We also define the contingent cone to $A$ at $x \in A$, denoted $T_A(x)$ as the set of directions $d \in \mathbb{R}^p$, such that there exist a sequence $(t_k) \in \mathbb{R}_{++}^N$, and a sequence $(d_k) \in (\mathbb{R}^p)^N$, satisfying $t_k \to 0$, $d_k \to d$, and $x + t_k d_k \in A$, for all $k \in \mathbb{N}$. Moreover, we say that a property $P$ holds “almost everywhere” (a.e.) on $A$, or equivalently “for almost all $x \in A$”, to denote that there exists a set $N$ of Lebesgue measure zero such that the property $P$ holds for all $x \in A \setminus N$. Finally, for any optimization problem $(P)$ we denote by $\text{val}(P)$ its optimal value.

## 2. Problem statement and Linear Programming formulations.

### 2.1. Definition of the minimal time control problem.
Let $n$ and $m$ be nonzero integers. We consider on $\mathbb{R}^n$ the control system

$$ (2.1) \qquad \dot{x}(t) = f(t, x(t), u(t)), $$

where $f: \mathbb{R} \times \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}^n$ is locally Lipschitz continuous, and where the controls are bounded measurable functions, defined on intervals $[t_0, t_1] \subset [0, T]$, and taking their values in a compact set $U$ of $\mathbb{R}^m$. Let $X$ and $K \subset X$ be compact sets of $\mathbb{R}^n$ and $x_0 \in X$. For $t_0, t_1 \ge 0$, a control $u$ is said admissible on $[t_0, t_1]$ whenever the solution $x(.)$ of (2.1), such that $x(t_0) = x_0$, is well defined on $[t_0, t_1]$ (in the sense that it is unique and absolutely continuous) and satisfies the constraints

$$ (2.2) \qquad (x(t), u(t)) \in X \times U, \quad \text{a.e. on } [t_0, t_1], $$
This is a particular case of the OCP with free final time [26], associated with the cost $\int_{t_0}^{t_1} \ell(t, x(t), u(t))dt$ for $\ell(t, x(t), u(t)) = 1$. The function $V^*$ is called the value function of this minimal time control problem: this describes the smallest time to reach the target set $K$, starting from $x_0$ at time $t_0$.

**ASSUMPTION 1.** For any $(t, x) \in [0, T] \times X$, the set $f(t, x, U)$ is convex.

We underline that we do not have any convexity assumption on the constraint set $X$ and on the target set $K$.

**REMARK 1.** Even if the dynamical system of interest does not satisfy Assumption 1, we can apply the present analysis to the convexified inclusion $\dot{x}(t) \in \text{conv } f(t, x(t), U)$. According to the Filippov-Ważewski relaxation Theorem [2, Th. 10.4.4], the trajectories of the original control problem are dense in the set of trajectories of the convexified inclusion. The trajectories of the convexified inclusion may be seen as the limit of chattering trajectories, i.e., when the control oscillates infinitely fast and where the constraint set is infinitesimally dilated.

**THEOREM 2.1.** Under Assumption 1, the minimal time control problem (2.1)-(2.4) associated with a starting point $(t_0, x_0) \in [0, T] \times X$ is either infeasible or admits an optimal trajectory.

*Proof.* We consider the case where a feasible trajectory exists. This is a direct application of [40, Th. 2.1], which, among others, characterizes the existence of an optimal trajectory for a control problem over a differential inclusion. To emphasize the correspondence with the notation of [40], we highlight that we apply the theorem with: the running cost function $\ell(t, x, p) = 1$, the terminal cost function $g(t, x) = 0$, the set-valued map $F(t, x) = f(t, x, U)$, the constraint set $A = [0, T] \times X$ and the target set $C = [0, T] \times K$. We underline that the assumptions (H1)-(H5) in [40] are satisfied here; more precisely, we highlight that our Assumption 1 enforces (H2) and the hypothesis that a feasible trajectory exists enforces (H4). □

**2.2. Hamilton-Jacobi-Bellman equation and subsolutions.** In optimal control theory, a well-known sufficient condition for a function $V$ to be the value function $V^*$ is to satisfy the Hamilton-Jacobi-Bellman (HJB) Partial Differential Equation (PDE). This PDE may be seen as a continuous time generalization of Bellman's dynamic programming optimality principle in discrete time [4]. In our minimal time control setting, the HJB PDE reads

$$ (2.5) \qquad \partial_t V(t, x) + \min_{u \in U} \{1 + \nabla_x V(t, x)^\top f(t, x, u)\} = 0, \quad \forall (t, x) \in [0, T] \times X $$

$$ (2.6) \qquad V(t, x) = 0, \quad \forall (t, x) \in [0, T] \times K. $$

In general, differentiable solutions of this PDE may not exist, so the concept of viscosity solutions is typically used [14]. Another approach to get around the lack of a differentiable solution to the HJB PDE consists in leveraging the concept of subsolutions [40], i.e., functions $V \in C^1(\mathbb{R}^{n+1})$ satisfying the following inequalities:

$$ (2.7) \qquad \partial_t V(t, x) + \min_{u \in U} \{1 + \nabla_x V(t, x)^\top f(t, x, u)\} \geq 0, \quad \forall (t, x) \in [0, T] \times X $$

$$ (2.8) \qquad V(t, x) \leq 0, \quad \forall (t, x) \in [0, T] \times K. $$

The following lemma states that any subsolution of the HJB PDE is an under-approximation of the value function.

**LEMMA 2.2.** For any $V \in C^1(\mathbb{R}^{n+1})$ satisfying Eqs. (2.7)-(2.8), the following holds:

$$ V(t, x) \leq V^*(t, x), \quad \forall (t, x) \in [0, T] \times X. $$

*Proof.* We take any $(t, x) \in [0, T] \times X$ and we consider that $V^*(t, x) < \infty$, since the case $V^*(t, x) = \infty$ is trivial. Hence, according to Th. 2.1, there exists an admissible control $\bar{u}(\cdot) \in \mathcal{U}(t, t_1, x)$ for $t_1 \in [t, T]$ such that $V^*(t, x) = t_1 - t$, and an associate trajectory $\bar{x}(t)$ such that $\bar{x}(t) = x$ and $\bar{x}(t_1) \in K$. We observe that $\frac{d}{dt}[V(t, \bar{x}(t))] = \partial_t V(t, \bar{x}(t)) + \nabla_x V(t, \bar{x}(t))^\top f(t, \bar{x}(t), \bar{u}(t)) \geq -1$ a.e. on $[t, t_1]$, the inequality holding since $V$ satisfies Eq. (2.7). By integration, we observe that $V(t_1, \bar{x}(t_1)) - V(t, x) \geq t - t_1 = -V^*(t, x)$, i.e., $V(t_1, \bar{x}(t_1)) + V^*(t, x) \geq V(t, x)$. As $V$ satisfies Eq. (2.8) and as $\bar{x}(t_1) \in K$, we observe that $0 \geq V(t_1, \bar{x}(t_1))$, and therefore $V^*(t, x) \geq V(t, x)$. □

**2.3. Infinite dimensional Linear Programming formulations.** In the rest of the paper, we consider a given point $x_0 \in X$, and we raise the issue of computing the minimal time from $x_0$ to $K$ and the associated control. We make the following assumption:
and satisfies the terminal state constraint

$$ (2.3) \qquad x(t_1) \in K. $$

We denote by $\mathcal{U}(t_0, t_1, x_0)$ the set of admissible controls on $[t_0, t_1]$. We consider the question of the minimal time problem from $x_0$ to $K$,

$$ (2.4) \qquad V^*(t_0, x_0) = \inf_{\substack{t_1 \in [t_0, T] \\ u(\cdot) \in \mathcal{U}(t_0, t_1, x_0)}} t_1 - t_0. $$

This is a particular case of the OCP with free final time [34], associated with the cost $\int_{t_0}^{t_1} \ell(t, x(t), u(t))dt$ for $\ell(t, x(t), u(t)) = 1$. The function $V^*$ is called the value function of this minimal time control problem: this describes the smallest time to reach the target set $K$, starting from $x_0$ at time $t_0$.

**ASSUMPTION 1.** For any $(t, x) \in [0, T] \times X$, the set $f(t, x, U)$ is convex.

We underline that we do not have any convexity assumption on the constraint set $X$ and on the target set $K$.

**REMARK 1.** Even if the dynamical system of interest does not satisfy Assumption 1, we can apply the present analysis to the convexified inclusion $\dot{x}(t) \in \text{conv} f(t, x(t), U)$. According to the Filippov-Ważewski relaxation Theorem [2, Th. 10.4.4], the trajectories of the original control problem are dense in the set of trajectories of the convexified inclusion [18]. Hence, as recalled in [24, Section 3.2], the trajectories of the convexified inclusion may be seen as the limit of chattering trajectories, i.e., when the control oscillates infinitely fast and where the constraint set is infinitesimally dilated.

**THEOREM 2.1.** Under Assumption 1, the minimal time control problem (2.1)-(2.4) associated with a starting point $(t_0, x_0) \in [0, T] \times X$ is either infeasible or admits an optimal trajectory.

*Proof.* We consider the case where a feasible trajectory exists. This is a direct application of [50, Th. 2.1], which, among others, characterizes the existence of an optimal trajectory for a control problem over a differential inclusion. To emphasize the correspondence with the notation of [50], we highlight that we apply the theorem with: the running cost function $\ell(t, x, p) = 1$, the terminal cost function $g(t, x) = 0$, the set-valued map $F(t, x) = f(t, x, U)$, the constraint set $A = [0, T] \times X$ and the target set $C = [0, T] \times K$. We underline that the assumptions (H1)-(H5) in [50] are satisfied here; more precisely, we highlight that our Assumption 1 enforces (H2) and the hypothesis that a feasible trajectory exists enforces (H4). □

**2.2. Hamilton-Jacobi-Bellman equation and subsolutions.** In optimal control theory, a well-known sufficient condition for a function $V$ to be the value function $V^*$ is to satisfy the Hamilton-Jacobi-Bellman (HJB) Partial Differential Equation (PDE). This PDE may be seen as a continuous time generalization of Bellman's dynamic programming optimality principle in discrete time [4]. In our minimal time control setting, the HJB PDE reads

$$ (2.5) \qquad \partial_t V(t, x) + \min_{u \in U} \{1 + \nabla_x V(t, x)^\top f(t, x, u)\} = 0, \quad \forall (t, x) \in [0, T] \times X $$

$$ (2.6) \qquad V(t, x) = 0, \quad \forall (t, x) \in [0, T] \times K. $$

In general, differentiable solutions of this PDE may not exist, so the concept of viscosity solutions is typically used, guaranteeing existence and uniqueness of a solution under Lemma 2.2 below [15]. Another approach to get around the lack of a differentiable solution to the HJB PDE consists in leveraging the concept of subsolutions [50], i.e., functions $V \in C^1(\mathbb{R}^{n+1})$ satisfying the following inequalities:

$$ (2.7) \qquad \partial_t V(t, x) + \min_{u \in U} \{1 + \nabla_x V(t, x)^\top f(t, x, u)\} \geq 0, \quad \forall (t, x) \in [0, T] \times X $$

$$ (2.8) \qquad V(t, x) \leq 0, \quad \forall (t, x) \in [0, T] \times K. $$

The following lemma states that any subsolution of the HJB PDE is an under-approximation of the value function.

**LEMMA 2.2.** For any $V \in C^1(\mathbb{R}^{n+1})$ satisfying Eqs. (2.7)-(2.8), the following holds:

$$ V(t, x) \leq V^*(t, x), \quad \forall (t, x) \in [0, T] \times X. $$
ASSUMPTION 2. There exists an admissible control $u \in U(0, t_1, x_0)$ associated with $t_1 \in [0, T]$. In other words, $V^*(0, x_0) \le t_1 < \infty$.

We consider the optimization problem of finding the subsolution of the HJB PDE that maximizes the evaluation in $(0, x_0)$. This problem may be cast as an infinite dimensional linear program:

$$ (D_F) \qquad \begin{array}{lll} \sup_{V \in \mathcal{F}} & V(0, x_0) & \\ \text{s.t.} & \partial_t V(t, x) + 1 + \nabla_x V(t, x)^\top f(t, x, u) \ge 0 & \forall (t, x, u) \in [0, T] \times X \times U \\ & V(t, x) \le 0 & \forall (t, x) \in [0, T] \times K, \end{array} $$

with $\mathcal{F} \in \{C^1(\mathbb{R}^{n+1}), C^\infty(\mathbb{R}^{n+1}), \mathbb{R}[t, x_1, \dots, x_n]\}$. For a given $V \in \mathcal{F}$, the feasibility in $(D_F)$ is clearly equivalent to the satisfaction of Eqs. (2.7)-(2.8). We also note that this infinite dimensional LP formulation corresponds to the dual LP formulation in [26]; in fact, this is the dual problem of an infinite dimensional LP formulation of the control problem based on occupation measures. According to the next theorem, the problem $(D_F)$ on $C^1$ functions has the same value as the minimal time control problem.

THEOREM 2.3. Under Assumption 1 and Assumption 2, and for $\mathcal{F} = C^1(\mathbb{R}^{n+1})$, the value of the LP formulation $(D_F)$ equals $V^*(0, x_0)$.

*Proof.* As for the proof of Th. 2.1, this is a direct application of [40, Th. 2.1], which also states the absence of duality gap between a control problem over a differential inclusion and a maximization problem over subsolutions of the HJB equation. We underline that the assumptions (H1)-(H5) in [40] are satisfied here; more precisely, our Assumption 1 enforces (H2) and our Assumption 2 enforces (H4). □

Theorem 2.5 extends this result by stating that we can require the subsolutions of the HJB equation to be in $C^\infty(\mathbb{R}^{n+1})$, while preserving the value of $(D_F)$. Before stating this theorem, we introduce an auxiliary lemma.

LEMMA 2.4. For any $V \in C^1(\mathbb{R}^{n+1})$ satisfying Eqs. (2.7)-(2.8) with feasibility error less or equal than $\eta \ge 0$, $V(t, x) + \eta(t - 1 - T)$ satisfies Eqs. (2.7)-(2.8).

*Proof.* We introduce $\tilde{V}(t, x) = V(t, x) + \eta(t - 1 - T)$. By assumption on $V(t, x)$, we have $\partial_t V(t, x) + 1 + \nabla_x V(t, x)^\top f(t, x, u) \ge -\eta$, for all $(t, x, u) \in [0, T] \times X \times U$. By linearity, and since $\partial_t (t - 1 - T) = 1$ and $\nabla_x (t - 1 - T) = 0$, $\partial_t \tilde{V}(t, x) + 1 + \nabla_x \tilde{V}(t, x)^\top f(t, x, u) \ge 0$. By assumption on $V(t, x)$, we have $V(t, x) \le \eta$, for all $(t, x) \in [0, T] \times K$. Hence, $\tilde{V}(t, x) \le \eta + \eta(t - 1 - T) \le \eta + \eta(T - 1 - T) = 0$ for all $(t, x) \in [0, T] \times K$. □

THEOREM 2.5. Under Assumption 1 and Assumption 2, and for $\mathcal{F} = C^\infty(\mathbb{R}^{n+1})$, the value of the LP formulation $(D_F)$ equals $V^*(0, x_0)$.

*Proof.* We consider $\mathcal{F} = C^\infty(\mathbb{R}^{n+1})$ and we use the notation $Y$ to denote the compact set $[0, T] \times X$. We fix $\epsilon > 0$, and we will prove that there exists $V \in C^\infty(\mathbb{R}^{n+1})$ that is feasible in $(D_F)$ and such that $V(0, x_0) \ge V^*(0, x_0) - \epsilon$. According to Th. 2.3, there exists $V_1 \in C^1(\mathbb{R}^{n+1})$ that is feasible in $(D_F)$ and such that $V_1(0, x_0) \ge V^*(0, x_0) - \frac{\epsilon}{2}$. For any $\sigma \in (0, 1]$, we introduce the mollified function $V_{1\sigma} = V_1 * \phi_\sigma \in C^\infty(\mathbb{R}^{n+1})$, where $\omega_\sigma$ is the standard mollifier defined as $\omega_\sigma(y) = \frac{1}{\sigma^{n+1}}\omega(y/\sigma)$, where $\omega(y) = \begin{cases} \xi e^{-\frac{1}{1-\|y\|^2}} & \text{if } \|y\| < 1 \\ 0 & \text{if } \|y\| \ge 1 \end{cases}$ for a given constant $\xi > 0$ such that $\int_{\mathbb{R}^{n+1}} \omega(y)dy = 1$. Hence, a simple change of variable shows that $\int_{\mathbb{R}^{n+1}} \omega_\sigma(y)dy = 1$. We also underline that $\omega_\sigma$ is non-negative and supported on the ball $B(0, \sigma)$. For any $y = (t, x) \in Y$ and any $\sigma \in (0, 1]$, we have that $|V_1(y) - V_{1\sigma}(y)| = |V_1(y) - \int_{B(0,\sigma)} V_1(y-h)\omega_\sigma(h)dh| = |\int_{B(0,\sigma)} (V_1(y) - V_1(y-h))\omega_\sigma(h)dh|$ as $\int_{B(0,\sigma)} \omega_\sigma(h)dh = 1$. We denote by $L_V$ an upper bound for the continuous function $\|\nabla V_1(y)\|_2$ over the compact set $\hat{Y} = \{y \in \mathbb{R}^{n+1}: d(y, Y) \le 1\}$, which is a Lipschitz constant for the function $V_1$. We deduce, by
**Proof.** We take any $(t, x) \in [0, T] \times X$ and we consider that $V^*(t, x) < \infty$, since the case $V^*(t, x) = \infty$ is trivial. Hence, according to Th. 2.1, there exists an admissible control $\bar{u}(\cdot) \in \mathcal{U}(t, t_1, x)$ for $t_1 \in [t, T]$ such that $V^*(t, x) = t_1 - t$, and an associate trajectory $\bar{x}(t)$ such that $\bar{x}(t) = x$ and $\bar{x}(t_1) \in K$. We observe that $\frac{d}{dt}[V(t, \bar{x}(t))] = \partial_t V(t, \bar{x}(t)) + \nabla_x V(t, \bar{x}(t))^{\top} f(t, \bar{x}(t), \bar{u}(t)) \ge -1$ a.e. on $[t, t_1]$, the inequality holding since $V$ satisfies Eq. (2.7). By integration, we observe that $V(t_1, \bar{x}(t_1)) - V(t, x) \ge t - t_1 = -V^*(t, x)$, i.e., $V(t_1, \bar{x}(t_1)) + V^*(t, x) \ge V(t, x)$. As $V$ satisfies Eq. (2.8) and as $\bar{x}(t_1) \in K$, we observe that $0 \ge V(t_1, \bar{x}(t_1))$, and therefore $V^*(t, x) \ge V(t, x)$. $\square$

**2.3. Infinite dimensional Linear Programming formulations.** In the rest of the paper, we consider a given point $x_0 \in X$, and we raise the issue of computing the minimal time from $x_0$ to $K$ and the associated control. We make the following assumption:

**ASSUMPTION 2.** There exists an admissible control $u \in \mathcal{U}(0, t_1, x_0)$ associated with $t_1 \in [0, T]$. In other words, $V^*(0, x_0) \le t_1 < \infty$.

We consider the optimization problem of finding the subsolution of the HJB PDE that maximizes the evaluation in $(0, x_0)$. This problem may be cast as an infinite dimensional linear program:

$$ (D_{\mathcal{F}}) \qquad \begin{array}{ll} \displaystyle \sup_{V \in \mathcal{F}} & V(0, x_0) \\ \text{s.t.} & \partial_t V(t, x) + 1 + \nabla_x V(t, x)^{\top} f(t, x, u) \ge 0 \quad \forall (t, x, u) \in [0, T] \times X \times U \\ & V(t, x) \le 0 \quad \forall (t, x) \in [0, T] \times K, \end{array} $$

with $\mathcal{F} \in \{C^1(\mathbb{R}^{n+1}), C^\infty(\mathbb{R}^{n+1}), \mathbb{R}[t, x_1, \dots, x_n]\}$. For a given $V \in \mathcal{F}$, the feasibility in $(D_{\mathcal{F}})$ is clearly equivalent to the satisfaction of Eqs. (2.7)-(2.8). We also note that this infinite dimensional LP formulation corresponds to the dual LP formulation in [34]; in fact, this is the dual problem of an infinite dimensional LP formulation of the control problem based on occupation measures. According to the next theorem, the problem $(D_{\mathcal{F}})$ on $C^1$ functions has the same value as the minimal time control problem.

**THEOREM 2.3.** Under Assumptions 1 and 2, and for $\mathcal{F} = C^1(\mathbb{R}^{n+1})$, the value of the LP formulation $(D_{\mathcal{F}})$ equals $V^*(0, x_0)$.

**Proof.** As for the proof of Th. 2.1, this is a direct application of [50, Th. 2.1], which also states the absence of duality gap between a control problem over a differential inclusion and a maximization problem over subsolutions of the HJB equation. We underline that the assumptions (H1)-(H5) in [50] are satisfied here; more precisely, our Assumption 1 enforces (H2) and our Assumption 2 enforces (H4). $\square$

Theorem 2.5 extends this result by stating that we can require the subsolutions of the HJB equation to be in $C^\infty(\mathbb{R}^{n+1})$, while preserving the value of $(D_{\mathcal{F}})$. Before stating this theorem, we introduce an auxiliary lemma.

**LEMMA 2.4.** For any $V \in C^1(\mathbb{R}^{n+1})$ and $\eta \ge 0$, if $\partial_t V(t, x) + 1 + \nabla_x V(t, x)^{\top} f(t, x, u) \ge -\eta$ for all $(t, x, u) \in [0, T] \times X \times U$ and $V(t, x) \le \eta$ for all $(t, x) \in [0, T] \times K$, then $V(t, x) + (t - 1 - T)\eta$ satisfies Eqs. (2.7)-(2.8).

**Proof.** We introduce $\tilde{V}(t, x) = V(t, x) + (t - 1 - T)\eta$. We have $\partial_t V(t, x) + 1 + \nabla_x V(t, x)^{\top} f(t, x, u) \ge -\eta$, for all $(t, x, u) \in [0, T] \times X \times U$. By linearity, and since $\partial_t (t - 1 - T) = 1$ and $\nabla_x (t - 1 - T) = 0$, it holds $\partial_t \tilde{V}(t, x) + 1 + \nabla_x \tilde{V}(t, x)^{\top} f(t, x, u) \ge 0$. Since we also have $V(t, x) \le \eta$, for all $(t, x) \in [0, T] \times K$ then $\tilde{V}(t, x) \le \eta + (t - 1 - T)\eta \le \eta + (t - 1 - T)\eta \le 0$ for all $(t, x) \in [0, T] \times K$. $\square$

**THEOREM 2.5.** Under Assumptions 1 and 2, and for $\mathcal{F} = C^\infty(\mathbb{R}^{n+1})$, the value of the LP formulation $(D_{\mathcal{F}})$ equals $V^*(0, x_0)$.

**Proof.** We consider $\mathcal{F} = C^\infty(\mathbb{R}^{n+1})$ and we use the notation $Y$ to denote the compact set $[0, T] \times X$. We fix $\epsilon > 0$, and we will prove that there exists $V \in C^\infty(\mathbb{R}^{n+1})$ that is feasible in $(D_{\mathcal{F}})$ and such that $V(0, x_0) \ge V^*(0, x_0) - \epsilon$. According to Th. 2.3, there exists $V_1 \in C^1(\mathbb{R}^{n+1})$ that is feasible in $(D_{\mathcal{F}})$ and such that $V_1(0, x_0) \ge V^*(0, x_0) - \frac{\epsilon}{2}$. For any $\sigma \in (0, 1]$, we introduce the mollified function [19] $V_{1\sigma} = V_1 * \omega_\sigma \in C^\infty(\mathbb{R}^{n+1})$, where $\omega_\sigma$ is the standard mollifier [21, p. 11] defined as $\omega_\sigma(y) = \frac{1}{\sigma^{n+1}} \omega(y/\sigma)$, where $\omega(y) = \begin{cases} \xi e^{-\frac{1}{1-\|y\|^2}} & \text{if } \|y\| < 1 \\ 0 & \text{if } \|y\| \ge 1 \end{cases}$ for a given constant $\xi > 0$ such that $\int_{\mathbb{R}^{n+1}} \omega(y)dy = 1$. Hence, a simple change of variable shows that $\int_{\mathbb{R}^{n+1}} \omega_\sigma(y)dy = 1$. We also underline that $\omega_\sigma$ is non-negative and supported on the ball $B(0, \sigma)$. For any $y = (t, x) \in Y$ and any $\sigma \in (0, 1]$, we have that $|V_1(y) - V_{1\sigma}(y)| = |V_1(y) - \int_{B(0,\sigma)} V_1(y-h)\omega_\sigma(h)dh| = |\int_{B(0,\sigma)} (V_1(y) - V_1(y-h))\omega_\sigma(h)dh|$ as $\int_{B(0,\sigma)} \omega_\sigma(h)dh = 1$. We denote by $L_V$ an upper bound for the continuous function $\|\nabla V_1(y)\|_2$ over the
triangular inequality and non-negativity of $\omega_\sigma$ that for any $y \in Y$,

$$ (2.9) \qquad |V_1(y) - V_{1\sigma}(y)| \le \int_{B(0,\sigma)} |V_1(y) - V_1(y-h)|\omega_\sigma(h)dh $$

$$ (2.10) \qquad \le \int_{B(0,\sigma)} L_V \|h\| \omega_\sigma(h) dh $$

$$ (2.11) \qquad \le L_V \sigma \int_{B(0,\sigma)} \|h/\sigma\| \omega(h/\sigma) \frac{1}{\sigma^{n+1}} dh $$

$$ (2.12) \qquad \le L_V \sigma \underbrace{\int_{B(0,1)} \|\tilde{h}\| \omega(\tilde{h}) d\tilde{h}}_{\text{constant, denoted } \mathcal{I}}. $$

By property of the mollifiers [21], we have $\partial_t V_{1\sigma}(y) = \partial_t (V_1 * \omega_\sigma) = (\partial_t V_1 * \omega_\sigma)$ for any $i \in \{t, x_1, \dots, x_n\}$. Therefore, $\partial_t V_{1\sigma}(y) = \int_{B(0,\sigma)} \partial_t V_1(y-h)\omega_\sigma(h)dh$ and $\nabla_x V_{1\sigma}(y) = \int_{B(0,\sigma)} \nabla_x V_1(y-h)\omega_\sigma(h)dh$. Using the equality $\int_{B(0,\sigma)} \omega_\sigma(h)dh = 1$, we deduce that for any $y \in Y$,

$$ (2.13) \qquad \partial_t V_{1\sigma}(y) + 1 + (\nabla_x V_{1\sigma}(y))^{\top} f(y, u) = \int_{B(0,\sigma)} (\partial_t V_1(y-h) + 1 + (\nabla_x V_1(y-h))^{\top} f(y, u))\omega_\sigma(h)dh $$

$$ (2.14) \qquad = \int_{B(0,\sigma)} (\partial_t V_1(y-h) + 1 + (\nabla_x V_1(y-h))^{\top} f(y-h, u))\omega_\sigma(h)dh $$

$$ (2.15) \qquad + \int_{B(0,\sigma)} \nabla_x V_1(y-h)^{\top} (f(y, u) - f(y-h, u))\omega_\sigma(h)dh. $$

We compute lower bounds for the two terms of the sum. We start with the second term: using Cauchy-Schwarz inequality, we notice that $\int_{B(0,\sigma)} \nabla_x V_1(y-h)^{\top} (f(y, u) - f(y-h, u))\omega_\sigma(h)dh \ge -\int_{B(0,\sigma)} \|V_1(y-h)\| \|f(y, u) - f(y-h, u)\|\omega_\sigma(h)dh$. Noticing that $\|\nabla_x V_1(y-h)\| \le L_V$, since $y-h \in \hat{Y}$ for any $h \in B(0,\sigma) \subset B(0,1)$, and introducing the Lipschitz constant $L_f$ for $f$, we have

$$ (2.16) \qquad \int_{B(0,\sigma)} \nabla_x V_1(y-h)^{\top} (f(y, u) - f(y-h, u))\omega_\sigma(h)dh \ge -L_V L_f \int_{B(0,\sigma)} \|h\| \omega_\sigma(h)dh = -L_V L_f \sigma \mathcal{I}. $$

We define $\eta = \frac{\epsilon}{2(T+2)}$. We introduce the compact set $Z = [0, T] \times X \times U$ and the family of compact sets $Z_\delta = \{z \in \mathbb{R}^N : d(z, Z) \le \delta\}$ for $\delta \in (0,1]$. For any $z = (y, u) \in Z_1$, we introduce $\psi(z) = \partial_t V_1(y) + 1 + (\nabla_x V_1(y))^{\top} f(y, u)$. The function $\psi(z)$ is continuous and according to Lemma A.1, there exists $\sigma_1 > 0$ such that $\min_{z \in Z_\sigma} \psi(z) \ge \min_{z \in Z} \psi(z) - \frac{\eta}{2}$ for any $\sigma \in (0, \sigma_1]$. By feasibility of $V_1$ in $(D_\mathcal{I})$, we know that $\min_{z \in Z} \psi(z) \ge 0$, which yields that $\psi(z) \ge -\frac{\eta}{2}$ for any $z \in Z_{\sigma_1}$. We deduce that

$$ (2.17) \qquad \int_{B(0,\sigma)} (\partial_t V_1(y-h) + 1 + (\nabla_x V_1(y-h))^{\top} f(y-h, u))\omega_\sigma(h)dh \ge -\int_{B(0,\sigma)} \frac{\eta}{2}\omega_\sigma(h)dh = -\frac{\eta}{2}, $$

since $(y-h, u) \in Z_\sigma$ for any $h \in B(0,\sigma)$. Combining the decomposition of Eqs. (2.13)-(2.15), with the lower bounds of Eq. (2.16) and (2.17), we deduce that

$$ (2.18) \qquad \partial_t V_{1\sigma}(y) + 1 + (\nabla_x V_{1\sigma}(y))^{\top} f(y, u) \ge -(L_V L_f \sigma \mathcal{I} + \frac{\eta}{2}), $$

for any $(y, u) = (t, x, u) \in [0, T] \times X \times U$ and $\sigma \in (0, \sigma_1]$. We define $\tilde{\sigma} = \min\{\sigma_1, \frac{\eta}{2L_V L_f \mathcal{I}}, \frac{\eta}{L_V \mathcal{I}}\}$. From Eq. (2.12) and Eq. (2.18), we deduce that

$$ (2.19) \qquad V_{1\tilde{\sigma}}(0, x_0) \ge V_1(0, x_0) - \eta \ge V^*(0, x_0) - \frac{\epsilon}{2} - \eta $$

$$ (2.20) \qquad V_{1\tilde{\sigma}}(t, x) \le V_1(t, x) + \eta \le \eta, \quad \forall (t, x) \in [0, T] \times K $$

$$ (2.21) \qquad \partial_t V_{1\tilde{\sigma}}(t, x) + 1 + \nabla_x V_{1\tilde{\sigma}}(t, x)^{\top} f(t, x, u) \ge -\eta, \quad \forall (t, x, u) \in [0, T] \times X \times U. $$

From Lemma 2.4, we deduce that $V(t, x) = V_{1\tilde{\sigma}}(t, x) + \eta(t-1-T) \in C^\infty(\mathbb{R}^{n+1})$ is feasible in $(D_\mathcal{I})$. From Eq. (2.19), we deduce that $V(0, x_0) \ge V^*(0, x_0) - \frac{\epsilon}{2} - \eta - (1+T)\eta$, and by definition of $\eta$, $V(0, x_0) \ge V^*(0, x_0) - \epsilon$. □
compact set $\hat{Y} = \{y \in \mathbb{R}^{n+1}: d(y, Y) \le 1\}$, which is a Lipschitz constant for the function $V_1$. We deduce, by triangular inequality and non-negativity of $\omega_\sigma$ that for any $y \in Y$,

$$ (2.9) \qquad \begin{aligned} |V_1(y) - V_{1\sigma}(y)| &\le \int_{B(0,\sigma)} |V_1(y) - V_1(y-h)|\omega_\sigma(h)dh \\ &\le \int_{B(0,\sigma)} L_V \|h\|\omega_\sigma(h)dh \\ &\le L_V \sigma \int_{B(0,\sigma)} \|h/\sigma\|\omega(h/\sigma) \frac{1}{\sigma^{n+1}}dh \\ &\le L_V \sigma \underbrace{\int_{B(0,1)} \|\tilde{h}\|\omega(\tilde{h})d\tilde{h}}_{\text{constant, denoted } \mathcal{I}}. \end{aligned} $$

By property of the mollifiers [29], we have $\partial_i V_{1\sigma}(y) = \partial_i(V_1 \star \omega_\sigma) = (\partial_i V_1 \star \omega_\sigma)$ for any $i \in \{t, x_1, \dots, x_n\}$. Therefore, $\partial_t V_{1\sigma}(y) = \int_{B(0,\sigma)} \partial_t V_1(y-h)\omega_\sigma(h)dh$ and $\nabla_x V_{1\sigma}(y) = \int_{B(0,\sigma)} \nabla_x V_1(y-h)\omega_\sigma(h)dh$. Using the equality $\int_{B(0,\sigma)} \omega_\sigma(h)dh = 1$, we deduce that for any $y \in Y$,

$$ (2.10) \qquad \begin{aligned} \partial_t V_{1\sigma}(y) + 1 + (\nabla_x V_{1\sigma}(y))^{\top} f(y, u) &= \int_{B(0,\sigma)} (\partial_t V_1(y-h) + 1 + (\nabla_x V_1(y-h))^{\top} f(y, u))\omega_\sigma(h)dh \\ &= \int_{B(0,\sigma)} (\partial_t V_1(y-h) + 1 + (\nabla_x V_1(y-h))^{\top} f(y-h, u))\omega_\sigma(h)dh \\ &\quad + \int_{B(0,\sigma)} \nabla_x V_1(y-h)^{\top} (f(y, u) - f(y-h, u))\omega_\sigma(h)dh. \end{aligned} $$

We compute lower bounds for the two terms of the sum. We start with the second term: using Cauchy-Schwarz inequality, we notice that $\int_{B(0,\sigma)} \nabla_x V_1(y-h)^{\top} (f(y, u) - f(y-h, u))\omega_\sigma(h)dh \ge -\int_{B(0,\sigma)} \|V_1(y-h)\|\|f(y, u) - f(y-h, u)\|\omega_\sigma(h)dh$. Noticing that $\|\nabla_x V_1(y-h)\| \le L_V$, since $y-h \in \hat{Y}$ for any $h \in B(0,\sigma) \subset B(0,1)$, and introducing the Lipschitz constant $L_f$ for $f$, we have

$$ (2.11) \qquad \int_{B(0,\sigma)} \nabla_x V_1(y-h)^{\top} (f(y, u) - f(y-h, u))\omega_\sigma(h)dh \ge -L_V L_f \int_{B(0,\sigma)} \|h\|\omega_\sigma(h)dh = -L_V L_f \sigma \mathcal{I}. $$

We define $\eta = \frac{\epsilon}{2(T+2)}$. We introduce the compact set $Z = [0, T] \times X \times U$ and the family of compact sets $Z_\delta = \{z \in \mathbb{R}^N: d(z, Z) \le \delta\}$ for $\delta \in (0,1]$. For any $z = (y, u) \in Z_1$, we introduce $\psi(z) = \partial_t V_1(y) + 1 + (\nabla_x V_1(y))^{\top} f(y, u)$. The function $\psi(z)$ is continuous and according to Lemma A.1, there exists $\sigma_1 > 0$ such that $\min_{z \in Z_\sigma} \psi(z) \ge \min_{z \in Z} \psi(z) - \frac{\eta}{2}$ for any $\sigma \in (0, \sigma_1]$. By feasibility of $V_1$ in $(D_\mathcal{I})$, we know that $\min_{z \in Z} \psi(z) \ge 0$, which yields that $\psi(z) \ge -\frac{\eta}{2}$ for any $z \in Z_{\sigma_1}$. We deduce that

$$ (2.12) \qquad \int_{B(0,\sigma)} (\partial_t V_1(y-h) + 1 + (\nabla_x V_1(y-h))^{\top} f(y-h, u))\omega_\sigma(h)dh \ge -\int_{B(0,\sigma)} \frac{\eta}{2}\omega_\sigma(h)dh = -\frac{\eta}{2}, $$

since $(y-h, u) \in Z_\sigma$ for any $h \in B(0,\sigma)$. Combining the decomposition of Eqs. (2.10), with the lower bounds of Eq. (2.11) and (2.12), we deduce that

$$ (2.13) \qquad \partial_t V_{1\sigma}(y) + 1 + (\nabla_x V_{1\sigma}(y))^{\top} f(y, u) \ge -(L_V L_f \sigma \mathcal{I} + \frac{\eta}{2}), $$

for any $(y, u) = (t, x, u) \in [0, T] \times X \times U$ and $\sigma \in (0, \sigma_1]$. We define $\tilde{\sigma} = \min\{\sigma_1, \frac{\eta}{2L_V L_f \sigma}, \frac{\eta}{L_V \mathcal{I}}\}$. From Eq. (2.9) and Eq. (2.13), we deduce that

$$ (2.14) \qquad \begin{aligned} V_{1\tilde{\sigma}}(0, x_0) &\ge V_1(0, x_0) - \eta \ge V^*(0, x_0) - \frac{\epsilon}{2} - \eta \\ V_{1\tilde{\sigma}}(t, x) &\le V_1(t, x) + \eta \le \eta, \quad \forall (t, x) \in [0, T] \times K \\ \partial_t V_{1\tilde{\sigma}}(t, x) + 1 + \nabla_x V_{1\tilde{\sigma}}(t, x)^{\top} f(t, x, u) &\ge -\eta, \quad \forall (t, x, u) \in [0, T] \times X \times U. \end{aligned} $$

From Lemma 2.4, we deduce that $V(t, x) = V_{1\tilde{\sigma}}(t, x) + (t-1-T)\eta \in C^\infty(\mathbb{R}^{n+1})$ is feasible in $(D_\mathcal{I})$. From Eq. (2.14), we deduce that $V(0, x_0) \ge V^*(0, x_0) - \frac{\epsilon}{2} - \eta - (1+T)\eta$, and by definition of $\eta$, $V(0, x_0) \ge V^*(0, x_0) - \epsilon$. □
The next theorem underlies the convergence proof of the hierarchy of semi-infinite problems in Sect. 3: if we restrict to polynomials HJB subsolutions, the value of the problem ($D_F$) remains unchanged.

**THEOREM 2.6.** *Under Assumption 1 and Assumption 2, and for $\mathcal{F} = \mathbb{R}[t, x_1, \dots, x_n]$, the value of the LP formulation ($D_F$) equals $V^*(0, x_0)$.*

*Proof.* We consider $\mathcal{F} = \mathbb{R}[t, x_1, \dots, x_n]$. For a given $\epsilon > 0$, and we will prove that there exists $V \in \mathbb{R}[t, x_1, \dots, x_n]$ that is feasible in ($D_F$) and such that $V(0, x_0) \ge V^*(0, x_0) - \epsilon$. According to Th. 2.5, there exists a function $Q \in C^\infty(\mathbb{R}^{n+1})$ which is a subsolution of the HJB equation and such that $Q(0, x_0) \ge V^*(0, x_0) - \frac{\epsilon}{2}$. We notice that $Q$ has a locally Lipschitz gradient. Therefore, we can apply Lemma A.2. This yields, in particular, that for any $\nu > 0$, there exists a polynomial $w \in \mathbb{R}[t, x_1, \dots, x_n]$ such that for all $(t, x) \in [0, T] \times X$, $|w(y) - w(y)| \le \nu$ and $|\partial_i w(t, x) - \partial_i V(t, x)| \le \nu$, $i \in \{t, x_1, \dots, x_N\}$. We deduce that $|\partial_t Q(t, x) + \nabla_x Q(t, x)^\top f(t, x, u) - \partial_t w(t, x) + \nabla_x w(t, x)^\top f(t, x, u)| \le |\partial_t Q(t, x) - \partial_t w(t, x)| + \sum_{i=1}^n |\partial_{x_i} Q(t, x) - \partial_{x_i} w(t, x)| M_i \le \nu(1 + \sum_{i=1}^n M_i)$, where $M_i = \max_{(t,x,u) \in [0,T] \times X \times U} |f_i(t,x,u)|$. Therefore, we observe that for all $(t, x, u) \in [0, T] \times X \times U$,

$$ (2.22) \qquad \partial_t w(t, x) + 1 + \nabla_x w(t, x)^\top f(t, x, u) \ge \partial_t Q(t, x) + 1 + \nabla_x Q(t, x)^\top f(t, x, u) - \nu \left(1 + \sum_{i=1}^n M_i\right) $$

$$ (2.23) \qquad \ge -\nu \left(1 + \sum_{i=1}^n M_i\right), $$

as $Q$ is a subsolution of the HJB equation. In summary, for $\nu = \eta(1 + \sum_{i=1}^n M_i)^{-1} \le \eta$,

$$ (2.24) \qquad w(0, x_0) \ge Q(0, x_0) - \eta \ge V^*(0, x_0) - \frac{\epsilon}{2} - \eta $$

$$ (2.25) \qquad w(t, x) \le Q(t, x) + \eta \le \eta, \quad \forall (t, x) \in [0, T] \times K $$

$$ (2.26) \qquad \partial_t w(t, x) + 1 + \nabla_x w(t, x)^\top f(t, x, u) \ge -\eta, \quad \forall (t, x, u) \in [0, T] \times X \times U, $$

the last inequality following from Eq. (2.23). Based on Eqs. (2.25)-(2.26) and Lemma 2.4, we notice that the polynomial $V(t, x) = w(t, x) + \eta(t - T - 1) \in \mathbb{R}[t, x_1, \dots, x_n]$ is feasible in the problem ($D_F$). Having defined $\eta = \frac{\epsilon}{2(T+2)}$, we see, based on Eq. (2.24), that it satisfies $V(0, x_0) \ge V^*(0, x_0) - \frac{\epsilon}{2} - \eta - \eta(T+1) = V^*(0, x_0) - \epsilon$.□

**3. Convex semi-infinite programming to compute near-optimal subsolutions.** For $\mathcal{F}$ being either $C^1(\mathbb{R}^{n+1})$, $C^\infty(\mathbb{R}^{n+1})$, or $\mathbb{R}[t, x_1, \dots, x_n]$, the linear program ($D_F$) is infinite dimensional, and thus, not tractable as it stands. Therefore, we next present a hierarchy of convex SIP problems that are solvable with a dedicated algorithm, to compute subsolutions to the HJB equation that are near optimal in the problem ($D_F$).

**3.1. A hierarchy of linear semi-infinite programs.** Instead of having an optimization space $\mathcal{F}$ that is infinite dimensional, we suggest to restrict to the finite dimensional subspaces $\mathbb{R}_d[t, x_1, \dots, x_n]$ of polynomials of degree bounded by $d$. This restricted dual problem is:

$$ (R_d) \qquad \begin{array}{ll} \displaystyle \sup_{V \in \mathbb{R}_d[t, x_1, \dots, x_n]} & V(0, x_0) \\ \text{s.t.} & \partial_t V(t, x) + 1 + \nabla_x V(t, x)^\top f(t, x, u) \ge 0 \quad \forall (t, x, u) \in [0, T] \times X \times U \\ & V(t, x) \le 0 \quad \forall (t, x) \in [0, T] \times K. \end{array} $$

In the rest of the paper, we will denote by $N$ the dimension of the vector space $\mathbb{R}_d[t, x_1, \dots, x_n]$ and $\Phi(t, x) \in \mathbb{R}^N$ a basis of this space. For both objects, there is indeed a dependence of $d$, that is implicit here for readability reasons. For any $V \in \mathbb{R}_d[t, x_1, \dots, x_n]$, we introduce the vector $\theta$ of the coordinates of $V$ in the basis $\Phi$. Hence, we have the relation

$$ (3.1) \qquad V(t, x) = \theta^\top \Phi(t, x) \in \mathbb{R}_d[t, x_1, \dots, x_n]. $$

Expressing problem ($R_d$) as an optimization problem over the vector of coefficients, it appears clearly that this is a linear semi-infinite program.

**PROPOSITION 3.1.** For $d \in \mathbb{N}^*$, problem ($R_d$) is a linear semi-infinite program, i.e, a linear program with a finite number of variables and an infinite number of constraints. More precisely, there exist a vector $c \in \mathbb{R}^N$, and a compact set $\mathcal{Y} \subset \mathbb{R}^{N+1}$ such that ($R_d$) reads

$$ (SIP) \qquad \begin{array}{ll} \displaystyle \sup_{\theta \in \mathbb{R}^N} & c^\top \theta \\ \text{s.t.} & a^\top \theta + b \le 0 \quad \forall (a, b) \in \mathcal{Y}. \end{array} $$
The next theorem underlies the convergence proof of the hierarchy of semi-infinite problems in Sect. 3: if we restrict to polynomials HJB subsolutions, the value of the problem ($D_F$) remains unchanged.

**THEOREM 2.6.** Under Assumptions 1 and 2, and for $\mathcal{F} = \mathbb{R}[t, x_1, \dots, x_n]$, the value of the LP formulation ($D_F$) equals $V^*(0, x_0)$.

*Proof.* We consider $\mathcal{F} = \mathbb{R}[t, x_1, \dots, x_n]$. For a given $\epsilon > 0$, and we will prove that there exists $V \in \mathbb{R}[t, x_1, \dots, x_n]$ that is feasible in ($D_F$) and such that $V(0, x_0) \ge V^*(0, x_0) - \epsilon$. According to Th. 2.5, there exists a function $Q \in C^\infty(\mathbb{R}^{n+1})$ which is a subsolution of the HJB equation and such that $Q(0, x_0) \ge V^*(0, x_0) - \frac{\epsilon}{2}$. We notice that $Q$ has a locally Lipschitz gradient. Therefore, we can apply Lemma A.2. This yields, in particular, that for any $\nu > 0$, there exists a polynomial $w \in \mathbb{R}[t, x_1, \dots, x_n]$ such that for all $(t, x) \in [0, T] \times X$, $|w(y) - w(y)| \le \nu$ and $|\partial_i w(t, x) - \partial_i V(t, x)| \le \nu$, $i \in \{t, x_1, \dots, x_N\}$. We deduce that $|\partial_t Q(t, x) + \nabla_x Q(t, x)^\top f(t, x, u) - \partial_t w(t, x) + \nabla_x w(t, x)^\top f(t, x, u)| \le |\partial_t Q(t, x) - \partial_t w(t, x)| + \sum_{i=1}^n |\partial_{x_i} Q(t, x) - \partial_{x_i} w(t, x)| M_i \le \nu(1 + \sum_{i=1}^n M_i)$, where $M_i = \max_{(t,x,u) \in [0,T] \times X \times U} |f_i(t,x,u)|$. Therefore, we observe that for all $(t, x, u) \in [0, T] \times X \times U$,

$$ (2.15) \qquad \partial_t w(t, x) + 1 + \nabla_x w(t, x)^\top f(t, x, u) \ge \partial_t Q(t, x) + 1 + \nabla_x Q(t, x)^\top f(t, x, u) - \nu \left(1 + \sum_{i=1}^n M_i\right) $$

$$ (2.16) \qquad \ge -\nu \left(1 + \sum_{i=1}^n M_i\right), $$

as $Q$ is a subsolution of the HJB equation. In summary, for $\nu = \eta(1 + \sum_{i=1}^n M_i)^{-1} \le \eta$,

$$ (2.17) \qquad w(0, x_0) \ge Q(0, x_0) - \eta \ge V^*(0, x_0) - \frac{\epsilon}{2} - \eta $$

$$ (2.18) \qquad w(t, x) \le Q(t, x) + \eta \le \eta, \quad \forall (t, x) \in [0, T] \times K $$

$$ (2.18) \qquad \partial_t w(t, x) + 1 + \nabla_x w(t, x)^\top f(t, x, u) \ge -\eta, \quad \forall (t, x, u) \in [0, T] \times X \times U, $$

the last inequality following from Eq. (2.16). Based on Eq. (2.18) and Lemma 2.4, we notice that the polynomial $V(t, x) = w(t, x) + (t - T - 1)\eta \in \mathbb{R}[t, x_1, \dots, x_n]$ is feasible in the problem ($D_F$). Having defined $\eta = \frac{\epsilon}{2(T+2)}$, we see, based on Eq. (2.17), that it satisfies $V(0, x_0) \ge V^*(0, x_0) - \frac{\epsilon}{2} - \eta - \eta(T + 1) = V^*(0, x_0) - \epsilon$. $\square$

**3. Convex semi-infinite programming to compute near-optimal subsolutions.** For $\mathcal{F}$ being either $C^1(\mathbb{R}^{n+1})$, $C^\infty(\mathbb{R}^{n+1})$, or $\mathbb{R}[t, x_1, \dots, x_n]$, the linear program ($D_F$) is infinite dimensional, and thus, not tractable as it stands. Therefore, we next present a hierarchy of convex semi-infinite programming (SIP) problems that are solvable with a dedicated algorithm, to compute subsolutions to the HJB equation that are near optimal in the problem ($D_F$).

**3.1. A hierarchy of linear semi-infinite programs.** Instead of having an optimization space $\mathcal{F}$ that is infinite dimensional, we suggest to restrict to the finite dimensional subspaces $\mathbb{R}_d[t, x_1, \dots, x_n]$ of polynomials of degree bounded by $d$. This restricted dual problem is:

$$ (D_d) \qquad \begin{array}{ll} \displaystyle \sup_{V \in \mathbb{R}_d[t, x_1, \dots, x_n]} & V(0, x_0) \\ \text{s.t.} & \partial_t V(t, x) + 1 + \nabla_x V(t, x)^\top f(t, x, u) \ge 0 \quad \forall (t, x, u) \in [0, T] \times X \times U \\ & V(t, x) \le 0 \quad \forall (t, x) \in [0, T] \times K. \end{array} $$

In the rest of the paper, we will denote by $N$ the dimension of the vector space $\mathbb{R}_d[t, x_1, \dots, x_n]$ and $\Phi(t, x) \in \mathbb{R}_d[t, x_1, \dots, x_n]^N$ a basis of this space. For both objects, there is indeed a dependence of $d$, that is implicit here for readability reasons. For any $V \in \mathbb{R}_d[t, x_1, \dots, x_n]$, we introduce the vector $\theta$ of the coordinates of $V$ in the basis $\Phi$. Hence, we have the relation

$$ (3.1) \qquad V(t, x) = \theta^\top \Phi(t, x) \in \mathbb{R}_d[t, x_1, \dots, x_n]. $$

Expressing problem ($D_d$) as an optimization problem over the vector of coefficients, it appears clearly that this is a linear semi-infinite program. For a comprehensive description of linear semi-infinite programming, the reader is referred to [22], which provides detailed insights into its mathematical formulation and dedicated algorithms. For a more general point of view on semi-infinite programming, the reader is referred to [16, 28, 38].
Proof. We define the vector $c = \Phi(0, x_0)$, and the compact sets

$$ (3.2) \qquad \mathcal{Y}_1 = \{(-\partial_t \Phi(t, x) - \nabla_x \Phi(t, x)^\top f(t, x, u), -1), (t, x, u) \in [0, T] \times X \times U\} $$

$$ (3.3) \qquad \mathcal{Y}_2 = \{(\Phi(t, x), 0), (t, x) \in [0, T] \times K\} $$

$$ (3.4) \qquad \mathcal{Y} = \mathcal{Y}_1 \cup \mathcal{Y}_2. $$

We see that for any $V_\theta(t, x) = \theta^\top \Phi(t, x) \in \mathbb{R}_d[t, x_1, \dots, x_n]$, $V_\theta(0, x_0) = c^\top \theta$, and $V_\theta(t, x)$ is feasible in $(R_d)$ if and only if $a^\top \theta + b \le 0$, for all $(a, b) \in \mathcal{Y}$. $\square$

We will see in the next section how to efficiently solve those semi-infinite programs. Prior to that, we state the convergence of this hierarchy of semi-infinite programs.

**THEOREM 3.2.** *The sequence $\text{val}(R_d)$ converges to $V^*(0, x_0)$ when $d \to \infty$.*

*Proof.* On the one hand, we introduce the notation $v_d = \text{val}(R_d)$. This sequence is obviously an increasing sequence, bounded above by $V^*(0, x_0)$. Hence, it converges to a value $\ell$, and any subsequence converges to $\ell \le V^*(0, x_0)$. On the other hand, Th. 2.6 guarantees that there exists a sequence of polynomials $w_k \in \mathbb{R}[t, x_1, \dots, x_n]$ that are feasible in $(D_\mathcal{F})$ and such that $w_k(0, x_0) \to_k V^*(0, x_0)$. By definition, we have $v_{d_k} \ge w_k(0, x_0)$, where $d_k = \text{deg}(w_k)$. Up to the extraction of a subsequence of $(w_k)$, we can assume that the sequence $d_k$ increasing, therefore $(v_{d_k})_{k \in \mathbb{N}}$ is a subsequence of $(v_d)_{d \in \mathbb{N}}$. As $v_{d_k} \to \ell$ and $w_k(0, x_0) \to_k V^*(0, x_0)$, we deduce that $\ell \ge V^*(0, x_0)$, which yields the equality $\ell = V^*(0, x_0)$. $\square$

**3.2. Regularization and solution of the semi-infinite programs.** We introduce a quadratic regularization in the semi-infinite program (SIP), yielding the following formulation depending on $\mu \in \mathbb{R}_{++}$:

$$ (SIP_\mu) \qquad \begin{array}{ll} \max_{\theta \in \mathbb{R}^N} & c^\top \theta - \frac{\mu}{2} \|\theta\|^2 \\ \text{s.t.} & a^\top \theta + b \le 0 \quad \forall (a, b) \in \mathcal{Y}. \end{array} $$

**PROPOSITION 3.3.** For any $\mu > 0$, the semi-infinite program $(SIP_\mu)$ has a unique optimal solution with value $\text{val}(SIP_\mu) \le V^*(0, x_0)$. Moreover, $\text{val}(SIP_\mu) \to_{\mu \to 0} \text{val}(SIP)$.

*Proof.* The feasible set of $(SIP_\mu)$ being convex, and the objective function being strongly concave, this optimization problem admits a unique maximum $\theta$. By definition, $\text{val}(SIP_\mu) = c^\top \theta - \frac{\mu}{2} \|\theta\|^2 \le c^\top \theta \le \text{val}(SIP)$, since $\theta$ is also feasible in the maximization problem (SIP). Additionally, $\text{val}(SIP) \le V^*(0, x_0)$, since any function $V$ feasible in $(R_d)$ satisfies $V(0, x_0) \le V^*(0, x_0)$. We also notice that the function $\mu \mapsto \text{val}(SIP_\mu)$ is decreasing, so it admits a limit $\ell$ at $0^+$, due to the aforementioned inequalities, $\ell \le \text{val}(SIP)$. For any $\mu, \epsilon > 0$, if we take $\theta_\epsilon$ an $\epsilon$-optimal solution in the problem (SIP), we see that $\text{val}(SIP) - \epsilon - \frac{\mu}{2} \|\theta_\epsilon\|^2 \le c^\top \theta_\epsilon - \frac{\mu}{2} \|\theta_\epsilon\|^2 \le \text{val}(SIP_\mu)$. For a fixed $\epsilon$, and taking $\mu \to 0^+$, we obtain $\text{val}(SIP) - \epsilon \le \ell$. This being true for any $\epsilon > 0$, we deduce that $\text{val}(SIP) \le \ell$, which proves the equality. $\square$

Setting the regularization parameter $\mu$ in practice implies a trade-off between the computational tractability of the semi-infinite program $(SIP_\mu)$ and the accuracy of the approximation of the original problem (SIP). To solve the formulation $(SIP_\mu)$, we propose to use a standard algorithm for convex semi-infinite programming, called cutting-plane (CP) algorithm [11]. To that extent, we assume to have an *separation oracle* computing, for any $\theta \in \mathbb{R}^N$,

$$ (3.5) \qquad \phi(\theta) = \max_{(a,b) \in \mathcal{Y}} a^\top \theta + b, $$

and an associate argmaximum. Solving the optimization problem in Eq. (3.5) may be computationally intensive, since the compact set $\mathcal{Y}$ may not be convex. Therefore, we only assume to have an oracle with relative optimality gap $\delta \in [0, 1)$ computing $(a, b) \in \mathcal{Y}$, such that $\phi(\theta) - (a^\top \theta + b) \le \delta |\phi(\theta)|$. We treat this oracle as a black box, regardless of its implementation, via global optimization, gridding, interval arithmetics or sampling for instance.
**PROPOSITION 3.1.** For $d \in \mathbb{N}^*$, problem $(D_d)$ is a linear semi-infinite program, i.e, a linear program with a finite number of variables and an infinite number of constraints. More precisely, there exist a vector $c \in \mathbb{R}^N$, and a compact set $\mathcal{Y} \subset \mathbb{R}^{N+1}$ such that $(D_d)$ reads

$$ (SIP) \qquad \begin{array}{ll} \sup_{\theta \in \mathbb{R}^N} & c^\top \theta \\ \text{s.t.} & a^\top \theta + b \le 0 \quad \forall (a, b) \in \mathcal{Y}. \end{array} $$

*Proof.* We define the vector $c = \Phi(0, x_0)$, and the compact sets

$$ (3.2) \qquad \mathcal{Y}_1 = \{(-\partial_t \Phi(t, x) - \nabla_x \Phi(t, x)^\top f(t, x, u), -1), (t, x, u) \in [0, T] \times X \times U\} $$

$$ (3.3) \qquad \mathcal{Y}_2 = \{(\Phi(t, x), 0), (t, x) \in [0, T] \times K\} $$

$$ (3.4) \qquad \mathcal{Y} = \mathcal{Y}_1 \cup \mathcal{Y}_2. $$

We see that for any $V_\theta(t, x) = \theta^\top \Phi(t, x) \in \mathbb{R}_d[t, x_1, \dots, x_n]$, $V_\theta(0, x_0) = c^\top \theta$, and $V_\theta(t, x)$ is feasible in $(D_d)$ if and only if $a^\top \theta + b \le 0$, for all $(a, b) \in \mathcal{Y}$. $\square$

We state the convergence of this hierarchy of semi-infinite programs.

**THEOREM 3.2.** *The sequence $\text{val}(D_d)$ converges to $V^*(0, x_0)$ when $d \to \infty$.*

*Proof.* On the one hand, we introduce the notation $v_d = \text{val}(D_d)$. This sequence is obviously an increasing sequence, bounded above by $V^*(0, x_0)$. Hence, it converges to a value $\ell$, and any subsequence converges to $\ell \le V^*(0, x_0)$. On the other hand, Th. 2.6 guarantees that there exists a sequence of polynomials $w_k \in \mathbb{R}[t, x_1, \dots, x_n]$ that are feasible in $(D_\ell)$ and such that $w_k(0, x_0) \to_k V^*(0, x_0)$. By definition, we have $v_{d_k} \ge w_k(0, x_0)$, where $d_k = \text{deg}(w_k)$; hence for all $d \ge d_1$, $v_d \ge v_{d_1} > -\infty$. Up to the extraction of a subsequence of $(w_k)$, we can assume that the sequence $d_k$ increasing, therefore $(v_{d_k})_{k \in \mathbb{N}}$ is a subsequence of $(v_d)_{d \in \mathbb{N}}$. As $v_{d_k} \to \ell$ and $w_k(0, x_0) \to_k V^*(0, x_0)$, we deduce that $\ell \ge V^*(0, x_0)$, which yields the equality $\ell = V^*(0, x_0)$. $\square$

The asymptotic convergence of $\text{val}(D_d)$ to $V^*(0, x_0)$ motivates the solution of such problems $(D_d)$. However, in practice, such problems can only be solved for relaxation degrees of limited size ($d \le 20$). For the end of Section 3, we set a degree $d \in \mathbb{N}$, and study how to solve the resulting problem $(D_d)$ efficiently; as mentioned above, this problem is denoted by (SIP). The question of the asymptotic complexity of solving problem $(D_d)$ as a function of degree $d$, as well as the rate of convergence of the hierarchy, are theoretical questions that are outside the scope of this article and would require further research. Instead, this article demonstrates the relevance and tractability of the approach through its numerical results in Section 5.

**3.2. Regularization and solution of the semi-infinite programs.** For reasons of computational tractability, we introduce a quadratic regularization in the semi-infinite program (SIP), yielding the following formulation depending on $\mu \in \mathbb{R}_{++}$:

$$ (SIP_\mu) \qquad \begin{array}{ll} \max_{\theta \in \mathbb{R}^N} & c^\top \theta - \frac{\mu}{2} \|\theta\|^2 \\ \text{s.t.} & a^\top \theta + b \le 0 \quad \forall (a, b) \in \mathcal{Y}. \end{array} $$

**PROPOSITION 3.3.** For any $\mu > 0$, the convex semi-infinite program $(SIP_\mu)$ has a unique optimal solution with value $\text{val}(SIP_\mu) \le V^*(0, x_0)$. Moreover, $\text{val}(SIP_\mu) \underset{\mu \to 0}{\to} \text{val}(SIP)$.

*Proof.* The feasible set of $(SIP_\mu)$ being convex, and the objective function being strongly concave, this optimization problem admits a unique maximum $\theta$. By definition, $\text{val}(SIP_\mu) = c^\top \theta - \frac{\mu}{2} \|\theta\|^2 \le c^\top \theta \le \text{val}(SIP)$, since $\theta$ is also feasible in the maximization problem (SIP). Additionally, $\text{val}(SIP) \le V^*(0, x_0)$, since any function $V$ feasible in $(D_d)$ satisfies $V(0, x_0) \le V^*(0, x_0)$. We also notice that the function $\mu \mapsto \text{val}(SIP_\mu)$ is decreasing, so it admits a limit $\ell$ at $0^+$, due to the aforementioned inequalities $\ell \le \text{val}(SIP)$. For any $\mu, \epsilon > 0$, if we take $\theta_\epsilon$ an $\epsilon$-optimal solution in the problem (SIP), we see that $\text{val}(SIP) - \epsilon - \frac{\mu}{2} \|\theta_\epsilon\|^2 \le c^\top \theta_\epsilon - \frac{\mu}{2} \|\theta_\epsilon\|^2 \le \text{val}(SIP_\mu)$. For a fixed $\epsilon$, and taking $\mu \to 0^+$, we obtain $\text{val}(SIP) - \epsilon \le \ell$. This being true for any $\epsilon > 0$, we deduce that $\text{val}(SIP) \le \ell$, which proves the equality. $\square$

This theorem shows that even by choosing a sufficiently small parameter $\mu$, the unique solution of $(SIP_\mu)$ is a quasi-optimal solution of (SIP). In practice, setting the regularization parameter $\mu$ implies a trade-off between the computational tractability of the semi-infinite program $(SIP_\mu)$ and the precision of the approximation of the original problem (SIP). To solve the formulation $(SIP_\mu)$, we propose to use a standard algorithm for convex
**Algorithm 3.1** Cutting-plane algorithm for (SIP<sub>μ</sub>)

Input: An oracle with parameter $\delta \in [0, 1)$, a tolerance $\epsilon \in \mathbb{R}_+$, a finite set $\mathcal{Y}^0 \subset \mathcal{Y}$, $k \leftarrow 0$

0: **while** true **do**

1: Compute $\theta^k$, the solution of the convex Quadratic Programming problem

$$ (3.6) \qquad \begin{array}{ll} \max_{\theta \in \mathbb{R}^N} & c^\top \theta - \frac{\mu}{2} \|\theta\|^2 \\ \text{s.t.} & a^\top \theta + b \le 0 \quad \forall (a, b) \in \mathcal{Y}^k. \end{array} $$

2: Call the oracle to compute $(a^k, b^k)$ an approximate solution of (3.5) with relative optimality gap $\delta$.

3: **if** $(a^k)^\top \theta^k + b^k \le \epsilon$ **then**

4: Return $\theta^k$.

5: **else**

6: $\mathcal{Y}^{k+1} \leftarrow \mathcal{Y}^k \cup \{(a^k, b^k)\}$

7: $k \leftarrow k + 1$

8: **end if**

9: **end while**

Before stating the termination and the convergence of Algorithm 3.1, we introduce the vector $\hat{\theta} \in \mathbb{R}^N$ of coordinates of the polynomial $\hat{v}(t, x) = t - 1 - T$ in the basis $\Phi(t, x)$, and we notice that this element helps obtaining feasible solutions since $\phi(\hat{\theta}) = -1$: Due to Lemma 2.4, we observe that if $\theta$ has a feasibility error less or equal than $\eta \ge 0$ in (SIP) and (SIP<sub>μ</sub>) then, $\theta + \eta \hat{\theta}$ is feasible in (SIP) and (SIP<sub>μ</sub>). For any $\mu > 0$, we define the convex and compact set $\mathcal{X}_\mu = \{\theta \in \mathbb{R}^N: c^\top \theta - \frac{\mu}{2} \|\theta\|^2 \ge c^\top \hat{\theta} - \frac{\mu}{2} \|\hat{\theta}\|^2\}$, and we define $R_\mu = \sup_{\theta \in \mathcal{X}_\mu} \|\theta\|$. Finally, we define the function $r_\mu(e) = e(1 + T + \mu R_\mu^2(1 + \frac{\epsilon}{2}))$. Note that $r_\mu(e) \underset{e \to 0}{\to} 0$.

**THEOREM 3.4.** If $\epsilon > 0$, Algorithm 3.1 stops after a finite number K of iterations, and $\theta^K + \frac{\epsilon}{1-\delta} \hat{\theta}$ is a feasible and $r_\mu(\frac{\epsilon}{1-\delta})$-optimal in (SIP<sub>μ</sub>). If $\epsilon = 0$, the alternative holds: (a) Algorithm 3.1 either stops after a finite number of iterations, and the last iterate is the optimal solution of (SIP<sub>μ</sub>), (b) Or it generates an infinite sequence, and the optimality gap and the feasibility error converge towards zero with an asymptotic rate in $O(\frac{1}{k})$.

**Proof.** First of all, we notice that during the execution of Algorithm 3.1, we necessarily have $\theta^k \in \mathcal{X}_\mu$, since $\hat{\theta}$ is a feasible solution in (3.6) with value $c^\top \hat{\theta} - \frac{\mu}{2} \|\hat{\theta}\|^2$, therefore by optimality of $\theta^k$ in (3.6), $c^\top \theta^k - \frac{\mu}{2} \|\theta^k\|^2 \ge c^\top \hat{\theta} - \frac{\mu}{2} \|\hat{\theta}\|^2$. The finite convergence of Algorithm 3.1 if $\epsilon > 0$, and the convergence rate in the case $\epsilon = 0$ (if no finite convergence) follows from [32, Th.1.1-1.2] (which is an extension of [11]): we apply these theorems to the problem (SIP<sub>μ</sub>) with the additional constraint $\theta \in \mathcal{X}_\mu$. As previously explained, this additional constraint does not change the execution of the algorithm, but it enables us to satisfy the compactness assumption of [32, Th.1.1-1.2]. We also note that the objective function is $\mu$-strongly concave, and that $\hat{\theta} \in \mathcal{X}_\mu$ is a strictly feasible point with respect to the semi-infinite constraints.

We finish the proof by showing that if Algorithm 3.1 stops at iteration K, then $\tilde{\theta} = \theta^K + \frac{\epsilon}{1-\delta} \hat{\theta}$ is feasible and $r_\mu(\frac{\epsilon}{1-\delta})$-optimal in (SIP<sub>μ</sub>). If Algorithm 3.1 stops at iteration K, this means that $(a^K)^\top \theta^K + b^K \le \epsilon$. If $\phi(\theta^K) \le 0$, then $\theta^K$ is feasible in (SIP<sub>μ</sub>), and so is $\tilde{\theta}$ due to Lemma 2.4. If $\phi(\theta^K) > 0$, then by property of the $\delta$-oracle, $(1 - \delta)\phi(\theta^K) \le (a^K)^\top \theta^K + b^K \le \epsilon$, and we deduce that the feasibility error is $\phi(\theta^K) \le \frac{\epsilon}{1-\delta}$. With Lemma 2.4, we deduce that $\tilde{\theta}$ is feasible in (SIP<sub>μ</sub>). We also note that

$$ (3.7) \qquad c^\top \tilde{\theta} - \frac{\mu}{2} \|\tilde{\theta}\|^2 = c^\top \theta^K + \frac{\epsilon}{1-\delta} c^\top \hat{\theta} - \frac{\mu}{2} \|\theta^K + \frac{\epsilon}{1-\delta} \hat{\theta}\|^2 $$

$$ (3.8) \qquad \ge c^\top \theta^K - \frac{\epsilon}{1-\delta}(1+T) - \frac{\mu}{2} \left( \|\theta^K\|^2 + \frac{2\epsilon}{1-\delta} \|\theta^K\| \|\hat{\theta}\| + \frac{\epsilon^2}{(1-\delta)^2} \|\hat{\theta}\|^2 \right), $$

since $c^\top \hat{\theta} = V_{\hat{\theta}}(0, x_0) = -(1 + T)$, and due to the Cauchy-Schwartz inequality. By optimality of $\theta^K$ in (3.6),
semi-infinite programming, called cutting-plane (CP) algorithm, that offers convergence rate guarantees in $O(\frac{1}{k})$ [12, 41, 42]. To that extent, we assume to have an *separation oracle* computing, for any $\theta \in \mathbb{R}^N$,

$$ (3.5) \qquad \phi(\theta) = \max_{(a,b) \in \mathcal{Y}} a^\top \theta + b, $$

and an associate argmaximum. Solving the optimization problem in Eq. (3.5) may be computationally intensive, since the compact set $\mathcal{Y}$ may not be convex. Therefore, we only assume to have an oracle with relative optimality gap $\delta \in [0, 1)$ computing $(a, b) \in \mathcal{Y}$, such that $\phi(\theta) - (a^\top \theta + b) \le \delta|\phi(\theta)|$. We treat this oracle as a black box, regardless of its implementation, via global optimization, gridding, interval arithmetics or sampling for instance.

**Algorithm 3.1** Cutting-plane algorithm for (SIP$_\mu$)

**Input:** An oracle with parameter $\delta \in [0, 1)$, a tolerance $\epsilon \in \mathbb{R}_+$, a finite set $\mathcal{Y}^0 \subset \mathcal{Y}$

0: $k \leftarrow 0$

1: **while** true **do**

2: Compute $\theta^k$, the solution of the convex Quadratic Programming problem

$$ (3.6) \qquad \begin{array}{ll} \max_{\theta \in \mathbb{R}^N} & c^\top \theta - \frac{\mu}{2} \|\theta\|^2 \\ \text{s.t.} & a^\top \theta + b \le 0 \quad \forall (a, b) \in \mathcal{Y}^k. \end{array} $$

3: Call the oracle to compute $(a^k, b^k)$ an approximate solution of (3.5) for $\theta = \theta^k$, with relative optimality gap $\delta$.

4: **if** $(a^k)^\top \theta^k + b^k \le \epsilon$ **then**

5: Return $\theta^k$.

6: **else**

7: $\mathcal{Y}^{k+1} \leftarrow \mathcal{Y}^k \cup \{(a^k, b^k)\}$

8: $k \leftarrow k + 1$

9: **end if**

10: **end while**

Before stating the termination and the convergence of Algorithm 3.1, we introduce the vector $\hat{\theta} \in \mathbb{R}^N$ of coordinates of the polynomial $\hat{v}(t, x) = t - 1 - T$ in the basis $\Phi(t, x)$, and we notice that this element helps obtaining feasible solutions since $\phi(\hat{\theta}) = -1$: Due to Lemma 2.4, we observe that if $\theta$ has a feasibility error less or equal than $\eta \ge 0$ in (SIP) and (SIP$_\mu$) then, $\theta + \eta \hat{\theta}$ is feasible in (SIP) and (SIP$_\mu$). For any $\mu > 0$, we define the convex and compact set $\mathcal{X}_\mu = \{\theta \in \mathbb{R}^N : c^\top \theta - \frac{\mu}{2} \|\theta\|^2 \ge c^\top \hat{\theta} - \frac{\mu}{2} \|\hat{\theta}\|^2\}$, and we define $R_\mu = \sup_{\theta \in \mathcal{X}_\mu} \|\theta\|$. Finally, we define the function $r_\mu(e) = e(1 + T + \mu R_\mu^2(1 + \frac{e}{2}))$. Note that $r_\mu(e) \xrightarrow{e \to 0} 0$.

**THEOREM 3.4.** If $\epsilon > 0$, Algorithm 3.1 stops after a finite number K of iterations, and $\theta^K + \frac{\epsilon}{1-\delta}\hat{\theta}$ is a feasible and $r_\mu(\frac{\epsilon}{1-\delta})$-optimal in (SIP$_\mu$). If $\epsilon = 0$, the alternative holds: (a) Algorithm 3.1 either stops after a finite number of iterations, and the last iterate is the optimal solution of (SIP$_\mu$), (b) Or it generates an infinite sequence, and the optimality gap and the feasibility error converge towards zero with an asymptotic rate in $O(\frac{1}{k})$.

*Proof.* First of all, we notice that during the execution of Algorithm 3.1, we necessarily have $\theta^k \in \mathcal{X}_\mu$, since $\hat{\theta}$ is a feasible solution in (3.6) with value $c^\top \hat{\theta} - \frac{\mu}{2} \|\hat{\theta}\|^2$, therefore by optimality of $\theta^k$ in (3.6), $c^\top \theta^k - \frac{\mu}{2} \|\theta^k\|^2 \ge c^\top \hat{\theta} - \frac{\mu}{2} \|\hat{\theta}\|^2$. The convergence rate of Algorithm 3.1, and the fact that $(a^k)^\top \theta^k + b^k \to 0$ if no finite convergence, follow from [42, Th.5-6] : we apply these theorems to the problem (SIP$_\mu$) with the additional constraint $\theta \in \mathcal{X}_\mu$. As previously explained, this additional constraint does not change the execution of the algorithm, but it enables us to satisfy the compactness assumption of [42, Th.5-6]. We also note that the objective function is $\mu$-strongly concave, and that $\hat{\theta} \in \mathcal{X}_\mu$ is a strictly feasible point with respect to the semi-infinite constraints. From $(a^k)^\top \theta^k + b^k \to 0$ we deduce that, if no finite convergence, then $\epsilon = 0$. By contrapositive, $\epsilon > 0$ implies finite convergence.

We finish the proof by showing that if Algorithm 3.1 stops at iteration K, then $\tilde{\theta} = \theta^K + \frac{\epsilon}{1-\delta}\hat{\theta}$ is feasible and $r_\mu(\frac{\epsilon}{1-\delta})$-optimal in (SIP$_\mu$). If Algorithm 3.1 stops at iteration K, this means that $(a^K)^\top \theta^K + b^K \le \epsilon$. If $\phi(\theta^K) \le 0$, then $\theta^K$ is feasible in (SIP$_\mu$), and so is $\tilde{\theta}$ due to Lemma 2.4. If $\phi(\theta^K) > 0$, then by property of the $\delta$-oracle, $(1 - \delta)\phi(\theta^K) \le (a^K)^\top \theta^K + b^K \le \epsilon$, and we deduce that the feasibility error is $\phi(\theta^K) \le \frac{\epsilon}{1-\delta}$. With
which is a relaxation of $(SIP_\mu)$, we know that $\text{val}(SIP_\mu) \le c^\top \theta^K - \frac{\mu}{2} \|\theta^K\|^2$. Applying this, we deduce that

$$ (3.9) \qquad c^\top \tilde{\theta} - \frac{\mu}{2} \|\tilde{\theta}\|^2 \ge \text{val}(SIP_\mu) - \frac{\epsilon}{1-\delta}(1+T) - \frac{\mu}{2} \left( \frac{2\epsilon}{1-\delta} \|\theta^K\| \|\hat{\theta}\| + \frac{\epsilon^2}{(1-\delta)^2} \|\hat{\theta}\|^2 \right) $$

$$ (3.10) \qquad \ge \text{val}(SIP_\mu) - \frac{\epsilon}{1-\delta}(1+T) - \mu R_\mu^2 \left( \frac{\epsilon}{1-\delta} + \frac{\epsilon^2}{2(1-\delta)^2} \right) $$

$$ (3.11) \qquad \ge \text{val}(SIP_\mu) - r_\mu \left( \frac{\epsilon}{1-\delta} \right), $$

the second inequality following from the fact that $\|\hat{\theta}\| \le R_\mu$ and $\|\theta^K\| \le R_\mu$, as $\hat{\theta}, \theta^K \in \mathcal{X}_\mu$.

□

**4. Feedback control based on approximate value functions.** In the previous section, we have seen how to compute subsolutions of the HJB equation based on convex semi-infinite programming, and how to deduce a lower bound on the minimal travel time. In this section, we focus on how subsolutions of the HJB equation, which approximate the value function $V^*$, enable one to recover a near-optimal control for the minimal time control problem (2.1)-(2.4).

**4.1. Controller design and existence of trajectories.** For a given continuously differentiable function $V \in C^1(\mathbb{R}^{n+1})$, we define the set-valued maps

$$ (4.1) \qquad \mathcal{U}_V(t, x) = \underset{u \in U}{\operatorname{argmin}} \nabla_x V(t, x)^\top f(t, x, u) $$

$$ (4.2) \qquad \mathcal{I}_V(t, x) = \{u \in \mathcal{U}_V(t, x) : f(t, x, u) \in T_X(x)\}, $$

where $T_X(x)$ is the contingent cone to $X$ at point $x$ (see Introduction). In line with previous works designing feedback controllers based on approximate value functions [18, 23], we are interested in the trajectories satisfying the following differential inclusion depending on the function $V \in C^1(\mathbb{R}^{n+1})$:

$$ (CL_V) \qquad \dot{x}_V(t) = f(t, x_V(t), u_V(t)) \text{ with } u_V(t) \in \mathcal{U}_V(t, x_V(t)). $$

Intuitively, such a feedback control pushes the system towards the descent direction of the function $V$. The following proposition confirms that, should the function $V \in C^1(\mathbb{R}^{n+1})$ be optimal in problem $(D_F)$, then any minimal time trajectory satisfies the differential inclusion $(CL_V)$ with respect to $V$.

**PROPOSITION 4.1.** *Under Assumptions 1-2, we consider an optimal trajectory $(x^*(\cdot), u^*(\cdot))$ of the minimal time control problem (2.1)-(2.4) starting from $(0, x_0)$, with hitting time $\tau^* = V^*(0, x_0)$. If the linear program $(D_F)$, for $\mathcal{F} = C^1(\mathbb{R}^{n+1})$, admits an optimal solution $V$, then, for almost every $t \in [0, \tau^*]$, *Proof.**

$$ (4.3) \qquad u^*(t) \in \mathcal{I}_V(t, x^*(t)) \subset \mathcal{U}_V(t, x^*(t)). $$

In particular, the trajectory $(x^*(\cdot), u^*(\cdot))$ satisfies the differential inclusion $(CL_V)$.

*Proof.* We define the function $\alpha(t) = V(t, x^*(t)) + t$, which is differentiable. We have that $\alpha'(t) = \partial_t V(t, x^*(t)) + 1 + \nabla_x V(t, x^*(t))^\top f(t, x^*(t), u^*(t))$, for almost all $t \in [0, \tau^*]$. Since $V$ is feasible in $(D_F)$, therefore satisfies Eq. (2.7), and since $(x^*(t), u^*(t)) \in X \times U$ a. e. on $[0, \tau^*]$, we know that $\alpha'(t) \ge 0$ a. e. on $[0, \tau^*]$. This proves that the differentiable function $\alpha(t)$ is non-decreasing function over $[0, \tau^*]$. By optimality of $V$ in $(D_F)$, and due to Th. 2.3 (Assumptions 1-2 are satisfied), $\alpha(0) = V(0, x_0) = \text{val}(D_F) = \tau^*$. Moreover, $\alpha(\tau^*) = \tau^* + V(\tau^*, x^*(\tau^*)) \le \tau^*$, since $V$ satisfies Eq. (2.8) and $x^*(\tau^*) \in K$. From $\alpha(\tau^*) \le \alpha(0)$, we obtain that $\alpha(t)$ is constant. Hence, $\partial_t V(t, x^*(t)) + 1 + \nabla_x V(t, x^*(t))^\top f(t, x^*(t), u^*(t)) = 0$, meaning

$$ (4.4) \qquad \nabla_x V(t, x^*(t))^\top f(t, x^*(t), u^*(t)) = -(\partial_t V(t, x^*(t)) + 1), \text{ a.e. on } [0, \tau^*]. $$

As $V$ satisfies Eq. (2.7), we have that $\nabla_x V(t, x^*(t))^\top f(t, x^*(t), u) \ge -(\partial_t V(t, x^*(t)) + 1)$ for all $t \in [0, \tau^*]$ and for all $u \in U$. Together with Eq. (4.4), we deduce that $u^*(t) \in \mathcal{U}_V(t, x^*(t))$ for almost all $t \in [0, \tau^*]$. Based on this fact, Lemma A.3 yields that for almost all $t \in [0, \tau^*]$, $f(t, x^*(t), u^*(t)) \in T_X(x^*(t))$. Therefore, for almost all $t \in [0, \tau^*]$, $u^*(t) \in \mathcal{I}_V(t, x^*(t))$. □

We just saw that whenever $V \in C^1(\mathbb{R}^{n+1})$ is optimal in the linear program $(D_F)$, any minimal time trajectory is a solution of the differential inclusion $(CL_V)$ associated with the function $V$. However, we may not be able to compute exactly such an optimal function in practice, especially because it may not exist. The next theorem states the existence of closed-loop trajectories following $(CL_V)$, for any function $V \in C^1(\mathbb{R}^{n+1})$.
Lemma 2.4, we deduce that $\tilde{\theta}$ is feasible in $(SIP_{\mu})$. We also note that

$$ (3.7) \qquad \begin{aligned} c^{\top}\tilde{\theta} - \frac{\mu}{2}\|\tilde{\theta}\|^2 &= c^{\top}\theta^K + \frac{\epsilon}{1-\delta}c^{\top}\hat{\theta} - \frac{\mu}{2}\|\theta^K\| + \frac{\epsilon}{1-\delta}\hat{\theta}\|^2 \\ &\geq c^{\top}\theta^K - \frac{\epsilon}{1-\delta}(1+T) - \frac{\mu}{2}\left(\|\theta^K\|^2 + \frac{2\epsilon}{1-\delta}\|\theta^K\|\|\hat{\theta}\| + \frac{\epsilon^2}{(1-\delta)^2}\|\hat{\theta}\|^2\right), \end{aligned} $$

since $c^{\top}\hat{\theta} = V_{\hat{\theta}}(0, x_0) = -(1+T)$, and due to the Cauchy-Schwartz inequality. By optimality of $\theta^K$ in (3.6), which is a relaxation of $(SIP_{\mu})$, we know that $\text{val}(SIP_{\mu}) \leq c^{\top}\theta^K - \frac{\mu}{2}\|\theta^K\|^2$. Applying this, we deduce that

$$ (3.8) \qquad \begin{aligned} c^{\top}\tilde{\theta} - \frac{\mu}{2}\|\tilde{\theta}\|^2 &\geq \text{val}(SIP_{\mu}) - \frac{\epsilon}{1-\delta}(1+T) - \frac{\mu}{2}\left(\frac{2\epsilon}{1-\delta}\|\theta^K\|\|\hat{\theta}\| + \frac{\epsilon^2}{(1-\delta)^2}\|\hat{\theta}\|^2\right) \\ &\geq \text{val}(SIP_{\mu}) - \frac{\epsilon}{1-\delta}(1+T) - \mu R_{\mu}^2 \left(\frac{\epsilon}{1-\delta} + \frac{\epsilon^2}{2(1-\delta)^2}\right) \\ &\geq \text{val}(SIP_{\mu}) - r_{\mu}\left(\frac{\epsilon}{1-\delta}\right), \end{aligned} $$

the second inequality following from the fact that $\|\hat{\theta}\| \leq R_{\mu}$ and $\|\theta^K\| \leq R_{\mu}$, as $\hat{\theta}, \theta^K \in \mathcal{X}_{\mu}$. $\square$

**4. Feedback control based on approximate value functions.** In the previous section, we have seen how to compute subsolutions of the HJB equation based on convex semi-infinite programming, and how to deduce a lower bound on the minimal travel time. In this section, we focus on how subsolutions of the HJB equation, which approximate the value function $V^*$, enable one to recover a near-optimal control for the minimal time control problem (2.1)-(2.4).

**4.1. Controller design and existence of trajectories.** For a given continuously differentiable function $V \in C^1(\mathbb{R}^{n+1})$, we define the set-valued maps

$$ (4.1) \qquad \mathcal{U}_V(t, x) = \underset{u \in U}{\operatorname{argmin}} \nabla_x V(t, x)^{\top} f(t, x, u) $$

$$ (4.2) \qquad \mathcal{I}_V(t, x) = \{u \in \mathcal{U}_V(t, x) : f(t, x, u) \in T_X(x)\}, $$

where $T_X(x)$ is the contingent cone to $X$ at point $x$ (see Introduction). In line with previous works designing feedback controllers based on approximate value functions [25, 31], we are interested in the trajectories satisfying the following differential inclusion depending on the function $V \in C^1(\mathbb{R}^{n+1})$:

$$ (CL_V) \qquad \dot{x}_V(t) = f(t, x_V(t), u_V(t)) \text{ with } u_V(t) \in \mathcal{U}_V(t, x_V(t)). $$

Intuitively, such a feedback control pushes the system towards the descent direction of the function $V$. The rest of this section is devoted to proving that this feedback control has interesting properties. First, the following proposition confirms that, should the function $V \in C^1(\mathbb{R}^{n+1})$ be optimal in problem $(D_{\mathcal{F}})$, then any minimal time trajectory satisfies the differential inclusion $(CL_V)$ with respect to $V$. Theorem 4.2 proves the existence of closed-loop trajectories following $(CL_V)$, and Theorem 4.3 the performance of the such trajectories as a function of the approximation quality of $V$: these results motivate the choice of such a feedback control.

**PROPOSITION 4.1.** *Under Assumptions 1-2, we consider an optimal trajectory $(x^*(\cdot), u^*(\cdot))$ of the minimal time control problem (2.1)-(2.4) starting from $(0, x_0)$, with hitting time $\tau^* = V^*(0, x_0)$. If the linear program $(D_{\mathcal{F}})$, for $\mathcal{F} = C^1(\mathbb{R}^{n+1})$, admits an optimal solution $V$, then, for almost every $t \in [0, \tau^*]$,*

$$ (4.3) \qquad u^*(t) \in \mathcal{I}_V(t, x^*(t)) \subset \mathcal{U}_V(t, x^*(t)). $$

In particular, the trajectory $(x^*(\cdot), u^*(\cdot))$ satisfies the differential inclusion $(CL_V)$.

*Proof.* We define the function $\alpha(t) = V(t, x^*(t)) + t$, which is differentiable. We have that $\alpha'(t) = \partial_t V(t, x^*(t)) + 1 + \nabla_x V(t, x^*(t))^{\top} f(t, x^*(t), u^*(t))$, for almost all $t \in [0, \tau^*]$. Since $V$ is feasible in $(D_{\mathcal{F}})$, therefore satisfies Eq. (2.7), and since $(x^*(t), u^*(t)) \in X \times U$ a.e. on $[0, \tau^*]$, we know that $\alpha'(t) \geq 0$ a.e. on $[0, \tau^*]$. This proves that the differentiable function $\alpha(t)$ is non-decreasing function over $[0, \tau^*]$. By optimality of $V$ in $(D_{\mathcal{F}})$, and due to Th. 2.3 (Assumptions 1-2 are satisfied), $\alpha(0) = V(0, x_0) = \text{val}(D_{\mathcal{F}}) = \tau^*$. Moreover,
**THEOREM 4.2.** Under Assumptions 1-2, if $V \in C^1(\mathbb{R}^{n+1})$ is such that for any $(t, x) \in \mathbb{R}_+ \times X$, $\mathcal{I}_V(t, x) \neq \emptyset$, then there exists a trajectory $(x_V(\cdot), u_V(\cdot))$ starting at $(0, x_0)$, satisfying the differential inclusion $(CL_V)$ over $[0, \infty)$ and such that $x_V(t) \in X$ for almost all $t \in [0, \infty)$.

*Proof.* We introduce an auxiliary control system to reduce to a time-invariant system with a convex control set, so as to fit in the setting of [1, Th. 6.6.6]. In what follows, we use the notation $y = (t, x)$ again. We introduce two objects: the set-valued map $\hat{U}(y) = \{f(y, u), u \in U\}$ and the function $\hat{f}(y, v) = \binom{1}{v}$ for $y \in \mathbb{R}^{n+1}$ and $v \in \mathbb{R}^n$. According to the terminology introduced in [1, Def. 6.1.3], $(\hat{U}, \hat{f})$ is a Marchaud control system, as (i) $\{(y, v) \in \mathbb{R}^{2n+1}: v \in \hat{U}(y)\}$ is closed, (ii) $\hat{f}$ is continuous, (iii) the velocity set $\{1\} \times f(y, U)$ is convex according to Assumption 1 and (iv) $\hat{f}$ has a linear growth, and so has $\hat{U}$ due to the fact that $f$ is Lipschitz continuous and $\mathcal{U}$ is bounded. We introduce $\mathcal{C} = \mathbb{R}_+ \times X$ and define the regulation map $\text{REG}(y) = \{v \in \hat{U}(y): \hat{f}(y, v) \in T_{\mathcal{C}}(y)\}$. We also introduce the set-valued map $\text{SEL}(y) = \underset{v \in \hat{U}(y)}{\operatorname{argmin}} \nabla_x V(y)^\top v$. We prove now that the graph of $\text{SEL}$ is closed.

For any converging sequence $(y_k, v_k) \to (\bar{y}, \bar{v})$ with $v_k \in \text{SEL}(y_k)$, we see that for all $k \in \mathbb{N}$, $v_k = f(y_k, u_k)$ for a given $u_k \in U$ and $\nabla_x V(y_k)^\top f(y_k, u_k) = h(y_k)$, where $h(y_k) = \min_{u \in U} \nabla_x V(y_k)^\top f(y_k, u)$. Up to extracting a subsequence of $u_k$, we can assume that $u_k \to \bar{u}$, as $U$ is compact. Note that $h$ is continuous, by application of the Maximum Theorem [1, Th. 2.1.6], in so far as (i) $(y, u) \mapsto \nabla_x V(y)^\top f(y, u)$ is continuous, therefore lower and upper semicontinuous, (ii) the set-valued map $M(y) = U$ is compact-valued, and lower and upper semicontinuous since it is constant. By continuity of $h$, $\nabla V$ and $f$, we conclude that $\nabla_x V(\bar{y})^\top \bar{v} = \nabla_x V(\bar{y})^\top f(\bar{y}, \bar{u}) = h(\bar{y}) = \min_{v \in \hat{U}(\bar{y})} \nabla_x V(\bar{y})^\top v$, meaning that $\bar{v} = f(\bar{y}, \bar{u}) \in \text{SEL}(\bar{y})$.

We notice that if $u \in \mathcal{I}_V(y)$, then $v = f(y, u) \in \text{REG}(y) \cap \text{SEL}(y)$. As $\mathcal{I}_V(y) \neq \emptyset$ for all $y \in \mathcal{C}$ (by assumption), $\text{REG}(y) \cap \text{SEL}(y) \neq \emptyset$. Together with the closedness of the graph of $\text{SEL}$, this means $\text{SEL}$ is a selection procedure of $\text{REG}$, according to the terminology of [1, Def. 6.5.2], and has convex values. We underline that $\text{REG}(y) \neq \emptyset$, for all $y \in \mathcal{C}$, i.e., $\mathcal{C}$ is a viability domain for $(\hat{U}, \hat{f})$. As $(0, x_0) \in \mathcal{C}$, [1, Th. 6.6.6] yields the existence of a solution $(y(\cdot), v(\cdot))$ such that $y(t) \in \mathcal{C}$, $v(t) \in \text{REG}(y(t))$ and

$$ (4.5) \qquad v(t) \in \text{SEL}(y(t)) \cap \text{REG}(y(t)), $$

for almost all $t \in [0, \infty)$. We notice first that $y_1(0) = 0$ and $\dot{y}_1(t) = 1$ for almost all $t \ge 0$, thus $y_1(t) = t$. Hence, we can indeed see $y(t)$ as $(t, x(t))$, with $x(0) = x_0$ and $\dot{x}(t) = v(t)$. Moreover, $v(t) = f(t, x(t), u(t))$ for a given $u(t) \in U$, since $v(t) \in \hat{U}(y(t)) = f(t, x(t), U)$ a.e. on $[0, \infty)$. We deduce from $v(t) \in \text{SEL}(y(t))$, which comes from (4.5), that $u(t) \in \mathcal{U}_V(t, x(t))$. Moreover, we deduce from $y(t) \in \mathcal{C}$ that $x(t) \in X$ a.e. on $[0, \infty)$. $\square$

**REMARK 2.** The condition $\mathcal{I}_V(t, x) \neq \emptyset$ in Th. 4.2 may appear restrictive, because it is not evident why a vector $f(t, x, u_V)$ minimizing $\nabla_x V(t, x)^\top f(t, x, u)$ over $u \in U$ would belong to $T_X(x)$. However, we have seen that under the hypotheses of Prop. 4.1, Eq. (4.3) yields $\mathcal{I}_V(t, x) \neq \emptyset$. Moreover, should the condition $\mathcal{I}_V(t, x) \neq \emptyset$ not be satisfied, we could enlarge the definition of $\mathcal{U}_V(t, x)$ in $\mathcal{U}_{V,\epsilon}(t, x) = \underset{u \in U}{\operatorname{argmin}}_\epsilon \nabla_x V(t, x)^\top f(t, x, u)$, so that for $\epsilon > 0$ large enough, $\mathcal{I}_{V,\epsilon}(t, x) = \{u \in \mathcal{U}_{V,\epsilon}(t, x): f(t, x, u) \in T_X(x)\} \neq \emptyset$.

**4.2. Performance of the feedback controller depending on the value function approximation error.** Previously, we introduced closed-loop trajectories satisfying the differential inclusion $(CL_V)$ with respect to a function $V \in C^1(\mathbb{R}^{n+1})$. In this section, we state some performance guarantees on those trajectories, depending on some properties of the function $V$. In the following, we assume that, up to an enlargement of the time horizon, the system can reach the target set starting from any initial condition $(t, x) \in [0, T] \times X$, and that the associated value function is Lipschitz.

**ASSUMPTION 3.** There exists a time $T^\sharp \ge T$ such that the minimal time control problem (2.1)-(2.4) defined over $[0, T^\sharp]$ has a value function $V^\sharp$ which takes finite values over $Y = [0, T] \times X$, and is Lipschitz continuous.

We emphasize that, under Assumption 3, $V^*(t, x) < \infty$ implies $V^*(t, x) = V^\sharp(t, x)$ for any $(t, x) \in Y$. Since $V^\sharp$ is Lipschitz continuous over $Y \subset \mathbb{R}^{n+1}$, it admits a Lipschitz continuous extension over $\mathbb{R}^{n+1}$ [16, Chap. 3, Th. 1]. We assimilate the value function and its extension on $\mathbb{R}^{n+1}$, such that we can speak about the Clarke's generalized derivative $\partial^c V^\sharp(y)$ of $V^\sharp$ at $y \in Y$. For any $V \in C^1(\mathbb{R}^{n+1})$, we introduce the notation

$$ (4.6) \qquad \|\nabla V - \nabla V^\sharp\|_{\infty} = \sup_{y \in Y} \sup_{g \in \partial^c V^\sharp(y)} \|\nabla V(y) - g\|_2. $$

We also define the constant $C_f = \sup_{(t,x,u) \in Y \times U} \|f(t, x, u)\| < \infty$.
$\alpha(\tau^*) = \tau^* + V(\tau^*, x^*(\tau^*)) \le \tau^*$, since $V$ satisfies Eq. (2.8) and $x^*(\tau^*) \in K$. From $\alpha(\tau^*) \le \alpha(0)$, we obtain that $\alpha(t)$ is constant. Hence, $\partial_t V(t, x^*(t)) + 1 + \nabla_x V(t, x^*(t))^{\top} f(t, x^*(t), u^*(t)) = 0$, meaning

$$ (4.4) \qquad \nabla_x V(t, x^*(t))^{\top} f(t, x^*(t), u^*(t)) = -(\partial_t V(t, x^*(t)) + 1), \text{ a.e. on } [0, \tau^*]. $$

As $V$ satisfies Eq. (2.7), we have that $\nabla_x V(t, x^*(t))^{\top} f(t, x^*(t), u) \ge -(\partial_t V(t, x^*(t)) + 1)$ for all $t \in [0, \tau^*]$ and for all $u \in U$. Together with Eq. (4.4), we deduce that $u^*(t) \in \mathcal{U}_V(t, x^*(t))$ for almost all $t \in [0, \tau^*]$. Based on this fact, Lemma A.3 yields that for almost all $t \in [0, \tau^*]$, $f(t, x^*(t), u^*(t)) \in T_X(x^*(t))$. Therefore, for almost all $t \in [0, \tau^*]$, $u^*(t) \in \mathcal{I}_V(t, x^*(t))$. $\square$

We just saw that whenever $V \in C^1(\mathbb{R}^{n+1})$ is optimal in the linear program $(D_F)$, any minimal time trajectory is a solution of the differential inclusion $(CL_V)$ associated with the function $V$. However, we may not be able to compute exactly such an optimal function in practice, especially because it may not exist. The next theorem states the existence of closed-loop trajectories following $(CL_V)$, for any function $V \in C^1(\mathbb{R}^{n+1})$.

**THEOREM 4.2.** Under Assumptions 1-2, if $V \in C^1(\mathbb{R}^{n+1})$ is such that for any $(t, x) \in \mathbb{R}_+ \times X$, $\mathcal{I}_V(t, x) \ne \emptyset$, then there exists a trajectory $(x_V(\cdot), u_V(\cdot))$ starting at $(0, x_0)$, satisfying the differential inclusion $(CL_V)$ over $[0, \infty)$ and such that $x_V(t) \in X$ for almost all $t \in [0, \infty)$.

*Proof.* We introduce an auxiliary control system to reduce to a time-invariant system with a convex control set, so as to fit in the setting of [1, Th. 6.6.6]. In what follows, we use the notation $y = (t, x)$ again. We introduce two objects: the set-valued map $\hat{U}(y) = \{f(y, u), u \in U\}$ and the function $\hat{f}(y, v) = \begin{pmatrix} 1 \\ v \end{pmatrix}$ for $y \in \mathbb{R}^{n+1}$ and $v \in \mathbb{R}^n$. According to the terminology introduced in [1, Def. 6.1.3], $(\hat{U}, \hat{f})$ is a Marchaud control system, as (i) $\{(y, v) \in \mathbb{R}^{2n+1}: v \in \hat{U}(y)\}$ is closed, (ii) $\hat{f}$ is continuous, (iii) the velocity set $\{1\} \times f(y, U)$ is convex according to Assumption 1 and (iv) $\hat{f}$ has a linear growth, and so has $\hat{U}$ due to the fact that $f$ is Lipschitz continuous and $\mathcal{U}$ is bounded. We introduce $\mathcal{C} = \mathbb{R}_+ \times X$ and define the regulation map $\text{REG}(y) = \{v \in \hat{U}(y): \hat{f}(y, v) \in T_C(y)\}$. We also introduce the set-valued map $\text{SEL}(y) = \underset{v \in \hat{U}(y)}{\operatorname{argmin}} \nabla_x V(y)^{\top} v$. We prove now that the graph of $\text{SEL}$ is closed.

For any converging sequence $(y_k, v_k) \to (\bar{y}, \bar{v})$ with $v_k \in \text{SEL}(y_k)$, we see that for all $k \in \mathbb{N}$, $v_k = f(y_k, u_k)$ for a given $u_k \in U$ and $\nabla_x V(y_k)^{\top} f(y_k, u_k) = h(y_k)$, where $h(y_k) = \min_{u \in U} \nabla_x V(y_k)^{\top} f(y_k, u)$. Up to extracting a subsequence of $u_k$, we can assume that $u_k \to \bar{u}$, as $U$ is compact. Note that $h$ is continuous, by application of the Maximum Theorem [1, Th. 2.1.6], in so far as (i) $(y, u) \mapsto \nabla_x V(y)^{\top} f(y, u)$ is continuous, therefore lower and upper semicontinuous, (ii) the set-valued map $M(y) = U$ is compact-valued, and lower and upper semicontinuous since it is constant. By continuity of $h$, $\nabla V$ and $f$, we conclude that $\nabla_x V(\bar{y})^{\top} \bar{v} = \nabla_x V(\bar{y})^{\top} f(\bar{y}, \bar{u}) = h(\bar{y}) = \min_{v \in \hat{U}(\bar{y})} \nabla_x V(\bar{y})^{\top} v$, meaning that $\bar{v} = f(\bar{y}, \bar{u}) \in \text{SEL}(\bar{y})$.

We notice that if $u \in \mathcal{I}_V(y)$, then $v = f(y, u) \in \text{REG}(y) \cap \text{SEL}(y)$. As $\mathcal{I}_V(y) \ne \emptyset$ for all $y \in \mathcal{C}$ (by assumption), $\text{REG}(y) \cap \text{SEL}(y) \ne \emptyset$. Together with the closedness of the graph of $\text{SEL}$, this means $\text{SEL}$ is a selection procedure of $\text{REG}$, according to the terminology of [1, Def. 6.5.2], and has convex values. We underline that $\text{REG}(y) \ne \emptyset$, for all $y \in \mathcal{C}$, i.e., $\mathcal{C}$ is a viability domain for $(\hat{U}, \hat{f})$. As $(0, x_0) \in \mathcal{C}$, [1, Th. 6.6.6] yields the existence of a solution $(y(\cdot), v(\cdot))$ such that $y(t) \in \mathcal{C}$, $v(t) \in \text{REG}(y(t))$ and

$$ (4.5) \qquad v(t) \in \text{SEL}(y(t)) \cap \text{REG}(y(t)), $$

for almost all $t \in [0, \infty)$. We notice first that $y_1(0) = 0$ and $\dot{y}_1(t) = 1$ for almost all $t \ge 0$, thus $y_1(t) = t$. Hence, we can indeed see $y(t)$ as $(t, x(t))$, with $x(0) = x_0$ and $\dot{x}(t) = v(t)$. Moreover, $v(t) = f(t, x(t), u(t))$ for a given $u(t) \in U$, since $v(t) \in \hat{U}(y(t)) = f(t, x(t), U)$ a.e. on $[0, \infty)$. We deduce from $v(t) \in \text{SEL}(y(t))$, which comes from (4.5), that $u(t) \in \mathcal{U}_V(t, x(t))$. Moreover, we deduce from $y(t) \in \mathcal{C}$ that $x(t) \in X$ a.e. on $[0, \infty)$. $\square$

**REMARK 2.** The condition $\mathcal{I}_V(t, x) \ne \emptyset$ in Th. 4.2 may appear restrictive, because it is not evident why a vector $f(t, x, u_V)$ minimizing $\nabla_x V(t, x)^{\top} f(t, x, u)$ over $u \in U$ would belong to $T_X(x)$. However, we have seen that under the hypotheses of Prop. 4.1, Eq. (4.3) yields $\mathcal{I}_V(t, x) \ne \emptyset$. Moreover, should the condition $\mathcal{I}_V(t, x) \ne \emptyset$ not be satisfied, we could enlarge the definition of $\mathcal{U}_V(t, x)$ in

$$ \mathcal{U}_{V,\epsilon}(t, x) = \left\{ u \in U : \nabla_x V(t, x)^{\top} f(t, x, u) \le \epsilon + \min_{u' \in U} \nabla_x V(t, x)^{\top} f(t, x, u') \right\}. $$

Then, for $\epsilon > 0$ large enough, we have $\mathcal{I}_{V,\epsilon}(t, x) = \{u \in \mathcal{U}_{V,\epsilon}(t, x): f(t, x, u) \in T_X(x)\} \ne \emptyset$.
**THEOREM 4.3.** Let $V \in C^1(\mathbb{R}^{n+1})$ be a continuously differentiable function, and let $(x_V(\cdot), u_V(\cdot))$ be a closed-loop trajectory starting at $(0, x_0)$ satisfying the differential inclusion $(CL_V)$ and the state constraints over $[0, T]$. We define $t_V = \sup\{t \in [0, T]: x_V([0, t]) \subset X \setminus K\}$. Then, under Assumptions 1-3,

$$ (4.7) \qquad V^{\sharp}(t, x_V(t)) \leq (\tau^* - t) + t 2(1 + C_f) \|\nabla V - \nabla V^{\sharp}\|_{\infty} \quad \forall t \in [0, t_V], $$

where $\tau^* = V^*(0, x_0) = V^*(0, x_0) \leq t_V$. In particular, we notice that

$$ (4.8) \qquad V^{\sharp}(\tau^*, x_V(\tau^*)) \leq 2\tau^*(1 + C_f) \|\nabla V - \nabla V^{\sharp}\|_{\infty}. $$

In Eq. (4.8), $V^{\sharp}(\tau^*, x_V(\tau^*))$ measures how far the closed-loop trajectory $(x_V(\cdot), u_V(\cdot))$ is from the target set $K$ at the moment when the time-optimal trajectory reaches $K$. As a corollary, we give a condition for the closed-loop trajectory $(x_V(\cdot), u_V(\cdot))$ to effectively reach the target set $K$, with a bounded delay compared to the time-optimal trajectory.

**COROLLARY 4.4.** Under the same hypotheses as Th. 4.3, if $\|\nabla V - \nabla V^{\sharp}\|_{\infty} \leq \frac{1-\tau^*/T}{2(1+C_f)}$, then

$$ (4.9) \qquad x_V(t_V) \in K \text{ with } t_V \in [\tau^*, \frac{1}{1-2(1+C_f)\|\nabla V - \nabla V^{\sharp}\|_{\infty}}\tau^*]. $$

We underline that the hitting time $t_V \geq \tau^*$ converges to the minimal time $\tau^*$, when the approximation error $\|\nabla V - \nabla V^{\sharp}\|_{\infty}$ vanishes.

**Proof of Th. 4.3 and Cor. 4.4.** For any Lipschitz continuous function $F: \mathbb{R}^{n+1} \to \mathbb{R}$, we recall that $\partial^c F(y)$ denote the Clarke's generalized derivative at $y$, and we define $H_F$ as

$$ (4.10) \qquad H_F(y) = 1 + \min_{\substack{u \in U \\ g \in \partial^c F(y)}} \{g^\top \begin{pmatrix} 1 \\ f(y, u) \end{pmatrix}\}. $$

The minimum is attained by continuity of the objective, and by compactness of $U$ and $\partial^c F(y)$ (see [12]). Note also that for any $V \in C^1(\mathbb{R}^{n+1})$, for any $y = (t, x) \in \mathbb{R}^{n+1}$, $H_V(t, x) = 1 + \partial_t V(t, x) + \min_{u \in U} \nabla_x V(t, x)^\top f(t, x, u)$, and the argmin is $\mathcal{U}_V(t, x)$. By application of the Maximum Theorem [1, Th. 2.1.6], we know that $H_F$ is lower semi-continuous, since (i) $\partial^c F(y)$ is a compact-valued and upper semi-continuous set-valued map [12], therefore so is $y \mapsto U \times \partial^c F(y)$, and (ii) $(y, u, g) \mapsto g^\top \begin{pmatrix} 1 \\ f(y, u) \end{pmatrix}$ is continuous.

First, we take any $y_1 = (t_1, x_1) \in [0, T] \times X \setminus K$, and we prove that $H_{V^{\sharp}}(t_1, x_1) \leq 0$. According to Assumption 3, $V^{\sharp}(t_1, x_1) < \infty$, and according to Th. 2.1 applied to the control system (2.1)-(2.4) on the interval $[0, T^{\sharp}]$, there exists an optimal trajectory $(x(\cdot), u(\cdot))$ over $[t_1, t_2]$ (with $t_2 > t_1$ since $x_0 \notin K$) starting from $(t_1, x_1)$. By definition, $V^{\sharp}(t_1, x_1) = t_2 - t_1$. We can also prove that for all $t \in [t_1, t_2]$, $V^{\sharp}(t, x(t)) = t_2 - t$: (i) the trajectory restricted to $[t, t_2]$, yields an admissible trajectory starting from $(t, x(t))$, therefore $V^{\sharp}(t, x(t)) \leq t_2 - t$, and (ii) for an optimal trajectory $(\tilde{x}(\cdot), \tilde{u}(\cdot))$ starting from $(t, x(t))$ over $[t, t_3]$, the trajectory following $(x(\cdot), u(\cdot))$ over $[t_1, t]$ and $(\tilde{x}(\cdot), \tilde{u}(\cdot))$ over $[t, t_3]$ is admissible and starting from $(t_1, x_1)$, therefore, $V^{\sharp}(t, x(t)) + (t - t_1) \geq V^{\sharp}(t_1, x_1) = t_2 - t_1$, giving $V^{\sharp}(t, x(t)) \geq t_2 - t$. As $\alpha(t) = V^{\sharp}(t, x(t)) = t_2 - t$ for all $t \in [t_1, t_2]$, we deduce that

Moreover, since $V^{\sharp}$ is Lipschitz continuous by assumption, and $t \mapsto (t, x(t))$ is Lipschitz continuous as $x(t)$ is differentiable a.e. with a bounded derivative, Lemma A.4 gives: $\alpha'(t) = \frac{d(V^{\sharp}(t, x(t)))}{dt} \geq \min_{g \in \partial^c V^{\sharp}(t, x(t))} g^\top \begin{pmatrix} 1 \\ \dot{x}(t) \end{pmatrix}$ a.e. on $[t_1, t_2]$. Using that $\dot{x}(t) = f(t, x(t), u(t))$ a.e. on $[t_1, t_2]$, and the definition of $H_{V^{\sharp}}$:

$$ (4.12) \qquad \alpha'(t) \geq \min_{g \in \partial^c V^{\sharp}(t, x(t))} g^\top \begin{pmatrix} 1 \\ f(t, x(t), u(t)) \end{pmatrix} \geq H_{V^{\sharp}}(t, x(t)) - 1, $$

a.e. on $[t_1, t_2]$. Combining this with Eq. (4.11), we deduce that for almost all $t \in [t_1, t_2]$, $H_{V^{\sharp}}(t, x(t)) \leq 0$. By lower semi-continuity of $H_{V^{\sharp}}$ (see above), and by continuity of $x(\cdot)$

$$ (4.13) \qquad H_{V^{\sharp}}(t_1, x_1) \leq 0. $$
### 4.2. Performance of the feedback controller depending on the value function approximation error.
Previously, we introduced closed-loop trajectories satisfying the differential inclusion $(CL_V)$ with respect to a function $V \in C^1(\mathbb{R}^{n+1})$. In this section, we state some performance guarantees on those trajectories, depending on some properties of the function $V$. In the following, we assume that, up to an enlargement of the time horizon, the system can reach the target set starting from any initial condition $(t, x) \in [0, T] \times X$, and that the associated value function is Lipschitz.

**ASSUMPTION 3.** There exists a time $T^\# \ge T$ such that the minimal time control problem (2.1)-(2.4) defined over $[0, T^\#]$ has a value function $V^\#$ which takes finite values over $Y = [0, T] \times X$, and is Lipschitz continuous.

We emphasize that, under Assumption 3, $V^*(t, x) < \infty$ implies $V^*(t, x) = V^\#(t, x)$ for any $(t, x) \in Y$. Since $V^\#$ is Lipschitz continuous over $Y \subset \mathbb{R}^{n+1}$, it admits a Lipschitz continuous extension over $\mathbb{R}^{n+1}$ [20, Chap. 3, Th. 1]. We assimilate the value function and its extension on $\mathbb{R}^{n+1}$, such that we can speak about the Clarke's generalized derivative $\partial^c V^\#(y)$ of $V^\#$ at $y \in Y$. For any $V \in C^1(\mathbb{R}^{n+1})$, we introduce the notation

$$ (4.6) \qquad \|\nabla V - \nabla V^\#\|_{\infty} = \sup_{y \in Y} \sup_{g \in \partial^c V^\#(y)} \|\nabla V(y) - g\|_2. $$

We also define the constant $C_f = \sup_{(t,x,u) \in Y \times U} \|f(t, x, u)\| < \infty$.

**THEOREM 4.3.** Let $V \in C^1(\mathbb{R}^{n+1})$ be a continuously differentiable function, and let $(x_V(\cdot), u_V(\cdot))$ be a closed-loop trajectory starting at $(0, x_0)$ satisfying the differential inclusion $(CL_V)$ and the state constraints over $[0, T]$. We define $t_V = \sup\{t \in [0, T]: x_V([0, t]) \subset X \setminus K\}$. Then, under Assumptions 1-3,

$$ (4.7) \qquad V^\#(t, x_V(t)) \le (\tau^* - t) + t(1 + C_f) \|\nabla V - \nabla V^\#\|_{\infty} \quad \forall t \in [0, t_V], $$

where $\tau^* = V^*(0, x_0) = V^\#(0, x_0) \le t_V$. In particular, we notice that

$$ (4.8) \qquad V^\#(\tau^*, x_V(\tau^*)) \le 2\tau^*(1 + C_f) \|\nabla V - \nabla V^\#\|_{\infty}. $$

In Eq. (4.8), $V^\#(\tau^*, x_V(\tau^*))$ measures how far the closed-loop trajectory $(x_V(\cdot), u_V(\cdot))$ is from the target set $K$ at the moment when the time-optimal trajectory reaches $K$. As a corollary, we give a condition for the closed-loop trajectory $(x_V(\cdot), u_V(\cdot))$ to effectively reach the target set $K$, with a bounded delay compared to the time-optimal trajectory.

**COROLLARY 4.4.** Under the same hypotheses as Th. 4.3, if $\|\nabla V - \nabla V^\#\|_{\infty} \le \frac{1-\tau^*/T}{2(1+C_f)}$, then

$$ (4.9) \qquad x_V(t_V) \in K \text{ with } t_V \in [\tau^*, \frac{1}{1-2(1+C_f)\|\nabla V - \nabla V^\#\|_{\infty}}\tau^*]. $$

We underline that the hitting time $t_V \ge \tau^*$ converges to the minimal time $\tau^*$, when the approximation error $\|\nabla V - \nabla V^\#\|_{\infty}$ vanishes.

*Proof of Th 4.3 and Cor. 4.4.* For any Lipschitz continuous function $F: \mathbb{R}^{n+1} \to \mathbb{R}$, we recall that $\partial^c F(y)$ denote the Clarke's generalized derivative at $y$, and we define $H_F$ as

$$ (4.10) \qquad H_F(y) = 1 + \min_{\substack{u \in U \\ g \in \partial^c F(y)}} \left\{ g^\top \begin{pmatrix} 1 \\ f(y, u) \end{pmatrix} \right\}. $$

The minimum is attained by continuity of the objective, and by compactness of $U$ and $\partial^c F(y)$ (see [13]). Note also that for any $V \in C^1(\mathbb{R}^{n+1})$, for any $y = (t, x) \in \mathbb{R}^{n+1}$, $H_V(t, x) = 1 + \partial_t V(t, x) + \min_{u \in U} \nabla_x V(t, x)^\top f(t, x, u)$, and the argmin is $\mathcal{U}_V(t, x)$. By application of the Maximum Theorem [1, Th. 2.1.6], we know that $H_F$ is lower semi-continuous, since (i) $\partial^c F(y)$ is a compact-valued and upper semi-continuous set-valued map [13], therefore so is $y \mapsto U \times \partial^c F(y)$, and (ii) $(y, u, g) \mapsto g^\top \begin{pmatrix} 1 \\ f(y, u) \end{pmatrix}$ is continuous.

First, we take any $y_1 = (t_1, x_1) \in [0, T] \times X \setminus K$, and we prove that $H_{V^\#}(t_1, x_1) \le 0$. According to Assumption 3, $V^\#(t_1, x_1) < \infty$, and according to Th. 2.1 applied to the control system (2.1)-(2.4) on the interval $[0, T^\#]$, there exists an optimal trajectory $(x(\cdot), u(\cdot))$ over $[t_1, t_2]$ (with $t_2 > t_1$ since $x_0 \notin K$) starting from
Second, still for any $(t_1, x_1) \in [0, T] \times X \setminus K$, we observe that there exists $(g, u_1) \in \partial^c V^\sharp(t_1, x_1) \times U$ such that $H_{V^\sharp}(t_1, x_1) = 1 + g^\top \begin{pmatrix} 1 \\ f(t_1, x_1, u_1) \end{pmatrix}$; indeed, we already mentioned that the minimum in (4.10) is attained. Therefore, for any $V \in C^1(\mathbb{R}^{n+1})$

$$ (4.14) \qquad 1 + \partial_t V(t_1, x_1) + \nabla_x V(t_1, x_1)^\top f(t_1, x_1, u_1) = H_{V^\sharp}(t_1, x_1) + (\nabla V(t_1, x_1) - g)^\top \begin{pmatrix} 1 \\ f(t_1, x_1, u_1) \end{pmatrix} \\ (4.15) \qquad \leq H_{V^\sharp}(t_1, x_1) + \|\nabla V - \nabla V^\sharp\|(1 + C_f), $$

the inequality being due to Cauchy-Schwartz inequality, and the definition of $\|\nabla V - \nabla V^\sharp\|$. We know that $H_V(t_1, x_1) \le \partial_t V(t_1, x_1) + 1 + \nabla_x V(t_1, x_1)^\top f(t_1, x_1, u_1)$ by definition of $H_V(t_1, x_1)$ (as $u_1 \in U$), therefore Eq. (4.15) gives $H_V(t_1, x_1) \le H_{V^\sharp}(t_1, x_1) + (1 + C_f)\|\nabla V - \nabla V^\sharp\|$. Using this inequality and Eq. (4.13), we deduce that for all $(t_1, x_1) \in [0, T] \times X \setminus K$,

$$ (4.16) \qquad H_V(t_1, x_1) \le (1 + C_f) \|\nabla V - \nabla V^\sharp\|. $$

Third, according to the hypotheses of the theorem, we take any $V \in C^1(\mathbb{R}^{n+1})$, and any closed-loop trajectory $(x_V(\cdot), u_V(\cdot))$ starting at $(0, x_0)$ satisfying the differential inclusion $(CL_V)$ and the state constraints over $[0, T]$. We, then, study the evolution of $V^\sharp$ over this trajectory. As $x_V(t)$ is Lipschitz continuous, Lemma A.4 yields the existence of $g(t) \in \partial^c V^\sharp(t, x_V(t))$ for almost all $t \in [0, T]$, such that

$$ (4.17) \qquad \frac{d}{dt} (V^\sharp(t, x_V(t))) \leq g(t)^\top \begin{pmatrix} 1 \\ f(t, x_V(t), u_V(t)) \end{pmatrix} \quad \text{a.e. on } [0, T]. $$

As $u_V(t) \in U_V(t, x_V(t))$, we know that $H_V(t, x_V(t)) = 1 + \nabla V(t, x_V(t))^\top \begin{pmatrix} 1 \\ f(t, x_V(t), u_V(t)) \end{pmatrix}$, and therefore,

$$ (4.18) \qquad \frac{d}{dt} (V^\sharp(t, x_V(t))) \le -1 + H_V(t, x_V(t)) + (g(t) - \nabla V(t, x_V(t)))^\top \begin{pmatrix} 1 \\ f(t, x_V(t), u_V(t)) \end{pmatrix} $$

We deduce, using Cauchy-Schwartz inequality and the definition of $\|\nabla V - \nabla V^\sharp\|$,

$$ (4.19) \qquad \frac{d}{dt} (V^\sharp(t, x_V(t))) \le -1 + H_V(t, x_V(t)) + (1 + C_f) \|\nabla V - \nabla V^\sharp\|, $$

for almost all $[0, T]$. Moreover, for all $t \in [0, t_V)$, $x_V(t) \notin K$. Therefore, we can apply Eq. (4.16) to deduce, in combination with Eq. (4.19), that for almost all $[0, t_V]$, $\frac{d}{dt} (V^\sharp(t, x_V(t))) \le -1 + 2(1 + C_f)\|\nabla V - \nabla V^\sharp\|$. By integration, we deduce that for all $t \in [0, t_V]$, $V^\sharp(t, x_V(t)) - \tau^* \le -t + 2t(1 + C_f)\|\nabla V - \nabla V^\sharp\|$, as $V^\sharp(0, x_V(0)) = V^\sharp(0, x_0) = \tau^*$. This proves Eq. (4.7).

Fourth and finally, we prove the corollary. Due to the definition of $t_V$, the following (non-exclusive) alternative holds: either $x_V(t_V) \in K$ or $t_V = T$. Moreover, if $\|\nabla V - \nabla V^\sharp\|_\infty \le \frac{1-\tau^*/T}{2(1+C_f)}$, then $V^\sharp(t, x_V(t)) - \tau^* \le -t + t(1-\tau^*/T)$ and $V^\sharp(t, x_V(t)) \le \tau^*(1-t/T)$ for all $t \in [0, t_V]$. We notice that if $t_V = T$, then $V^\sharp(t_V, x_V(t_V)) \le 0$, i.e., $x_V(t_V) \in K$. Coming to the aforementioned alternative, we deduce that $x_V(t_V) \in K$. Moreover, this fact combined with Eq. (4.7) gives us that $0 \le (\tau^* - t_V) + t_V 2(1 + C_f)\|\nabla V - \nabla V^\sharp\|_\infty$, hence

$$ (4.20) \qquad t_V (1 - 2(1 + C_f) \|\nabla V - \nabla V^\sharp\|_\infty) \le \tau^*. $$

By assumption, $1 - 2(1 + C_f)\|\nabla V - \nabla V^\sharp\|_\infty \ge \tau^*/T > 0$, we can thus divide Eq. (4.20) by this quantity to obtain the result of the corollary: $t_V \le \tau^*/(1 - 2(1 + C_f)\|\nabla V - \nabla V^\sharp\|_\infty)$. $\square$

In the previous theorem and the corollary, we saw that the suboptimality, in terms of hitting time, of a closed-loop trajectory $(x_V(\cdot), u_V(\cdot))$ satisfying the differential inclusion $(CL_V)$ decreases as the approximation error $\|\nabla V - \nabla V^\sharp\|_\infty$ decreases. Furthermore, we see that the closed-loop trajectory comes close to optimality when the approximation error vanishes. We now study a sufficient condition under which the approximation $\|\nabla V_d - \nabla V^\sharp\|_\infty$ can be made arbitrarily small, using a polynomial $V_d(t, x)$ of sufficiently large degree $d \in \mathbb{N}$.
(t₁, x₁). By definition, $V^{\#}(t_1, x_1) = t_2 - t_1$. We can also prove that for all $t \in [t_1, t_2]$, $V^{\#}(t, x(t)) = t_2 - t$: (i) the trajectory restricted to $[t, t_2]$, yields an admissible trajectory starting from $(t, x(t))$, therefore $V^{\#}(t, x(t)) \le t_2 - t$, and (ii) for an optimal trajectory $(\tilde{x}(\cdot), \tilde{u}(\cdot))$ starting from $(t, x(t))$ over $[t, t_3]$, the trajectory following $(x(\cdot), u(\cdot))$ over $[t_1, t]$ and $(\tilde{x}(\cdot), \tilde{u}(\cdot))$ over $[t, t_3]$ is admissible and starting from $(t_1, x_1)$, therefore, $V^{\#}(t, x(t)) + (t - t_1) \ge V^{\#}(t_1, x_1) = t_2 - t_1$, giving $V^{\#}(t, x(t)) \ge t_2 - t$. As $\alpha(t) = V^{\#}(t, x(t)) = t_2 - t$ for all $t \in [t_1, t_2]$, we deduce that

$$ (4.11) \qquad \alpha'(t) = -1 \quad \text{a.e. on } [t_1, t_2]. $$

Moreover, since $V^{\#}$ is Lipschitz continuous by assumption, and $t \mapsto (t, x(t))$ is Lipschitz continuous as $x(t)$ is differentiable a.e. with a bounded derivative, Lemma A.4 gives: $\alpha'(t) = \frac{d(V^{\#}(t, x(t)))}{dt} \ge \min_{g \in \partial^c V^{\#}(t, x(t))} g^{\top} \begin{pmatrix} 1 \\ \dot{x}(t) \end{pmatrix}$ a.e. on $[t_1, t_2]$. Using that $\dot{x}(t) = f(t, x(t), u(t))$ a.e. on $[t_1, t_2]$, and the definition of $H_{V^{\#}}$:

$$ (4.12) \qquad \alpha'(t) \ge \min_{g \in \partial^c V^{\#}(t, x(t))} g^{\top} \begin{pmatrix} 1 \\ f(t, x(t), u(t)) \end{pmatrix} \ge H_{V^{\#}}(t, x(t)) - 1, $$

a.e. on $[t_1, t_2]$. Combining this with Eq. (4.11), we deduce that for almost all $t \in [t_1, t_2]$, $H_{V^{\#}}(t, x(t)) \le 0$. By lower semi-continuity of $H_{V^{\#}}$ (see above), and by continuity of $x(\cdot)$

$$ (4.13) \qquad H_{V^{\#}}(t_1, x_1) \le 0. $$

Second, still for any $(t_1, x_1) \in [0, T] \times X \setminus K$, we observe that there exists $(g, u_1) \in \partial^c V^{\#}(t_1, x_1) \times U$ such that $H_{V^{\#}}(t_1, x_1) = 1 + g^{\top} \begin{pmatrix} 1 \\ f(t_1, x_1, u_1) \end{pmatrix}$; indeed, we already mentioned that the minimum in (4.10) is attained. Therefore, for any $V \in C^1(\mathbb{R}^{n+1})$

$$ (4.14) \qquad 1 + \partial_t V(t_1, x_1) + \nabla_x V(t_1, x_1)^\top f(t_1, x_1, u_1) = H_{V^{\#}}(t_1, x_1) + (\nabla V(t_1, x_1) - g)^\top \begin{pmatrix} 1 \\ f(t_1, x_1, u_1) \end{pmatrix} \\ \le H_{V^{\#}}(t_1, x_1) + \|\nabla V - \nabla V^{\#}\|(1 + C_f), $$

the inequality being due to Cauchy-Schwartz inequality, and the definition of $\|\nabla V - \nabla V^{\#}\|$. We know that $H_V(t_1, x_1) \le \partial_t V(t_1, x_1) + 1 + \nabla_x V(t_1, x_1)^\top f(t_1, x_1, u_1)$ by definition of $H_V(t_1, x_1)$ (as $u_1 \in U$), therefore Eq. (4.14) gives $H_V(t_1, x_1) \le H_{V^{\#}}(t_1, x_1) + (1 + C_f)\|\nabla V - \nabla V^{\#}\|$. Using this inequality and Eq. (4.13), we deduce that for all $(t_1, x_1) \in [0, T] \times X \setminus K$,

$$ (4.15) \qquad H_V(t_1, x_1) \le (1 + C_f) \|\nabla V - \nabla V^{\#}\|. $$

Third, according to the hypotheses of the theorem, we take any $V \in C^1(\mathbb{R}^{n+1})$, and any closed-loop trajectory $(x_V(\cdot), u_V(\cdot))$ starting at $(0, x_0)$ satisfying the differential inclusion $(CL_V)$ and the state constraints over $[0, T]$. We, then, study the evolution of $V^{\#}$ over this trajectory. As $x_V(t)$ is Lipschitz continuous, Lemma A.4 yields the existence of $g(t) \in \partial^c V^{\#}(t, x_V(t))$ for almost all $t \in [0, T]$, such that

$$ (4.16) \qquad \frac{d}{dt} (V^{\#}(t, x_V(t))) \le g(t)^{\top} \begin{pmatrix} 1 \\ f(t, x_V(t), u_V(t)) \end{pmatrix} \quad \text{a.e. on } [0, T]. $$

As $u_V(t) \in U_V(t, x_V(t))$, we know that $H_V(t, x_V(t)) = 1 + \nabla V(t, x_V(t))^{\top} \begin{pmatrix} 1 \\ f(t, x_V(t), u_V(t)) \end{pmatrix}$, and therefore,

$$ (4.17) \qquad \frac{d}{dt} (V^{\#}(t, x_V(t))) \le -1 + H_V(t, x_V(t)) + (g(t) - \nabla V(t, x_V(t)))^{\top} \begin{pmatrix} 1 \\ f(t, x_V(t), u_V(t)) \end{pmatrix}, $$

We deduce, using Cauchy-Schwartz inequality and the definition of $\|\nabla V - \nabla V^{\#}\|$,

$$ (4.18) \qquad \frac{d}{dt} (V^{\#}(t, x_V(t))) \le -1 + H_V(t, x_V(t)) + (1 + C_f) \|\nabla V - \nabla V^{\#}\|, $$

for almost all $t \in [0, T]$. Moreover, for all $t \in [0, t_V)$, $x_V(t) \notin K$. Therefore, we can apply Eq. (4.15) to deduce, in combination with Eq. (4.18), that for almost all $[0, t_V]$, $\frac{d}{dt} (V^{\#}(t, x_V(t))) \le -1 + 2(1 + C_f) \|\nabla V - \nabla V^{\#}\|$. By
integration, we deduce that for all $t \in [0, t_V]$, $V^\#(t, x_V(t)) - \tau^* \le -t + 2t(1 + C_f) \|\nabla V - \nabla V^\#\|$, as $V^\#(0, x_V(0)) = V^\#(0, x_0) = \tau^*$. This proves Eq. (4.7).

Fourth and finally, we prove the corollary. Due to the definition of $t_V$, the following (non-exclusive) alternative holds: either $x_V(t_V) \in K$ or $t_V = T$. Moreover, if $\|\nabla V - \nabla V^\#\|_\infty \le \frac{1-\tau^*/T}{2(1+C_f)}$, then $V^\#(t, x_V(t)) - \tau^* \le -t + t(1-\tau^*/T)$ and $V^\#(t, x_V(t)) \le \tau^*(1-t/T)$ for all $t \in [0, t_V]$. We notice that if $t_V = T$, then $V^\#(t_V, x_V(t_V)) \le 0$, i.e., $x_V(t_V) \in K$. Coming to the aforementioned alternative, we deduce that $x_V(t_V) \in K$. Moreover, this fact combined with Eq. (4.7) gives us that $0 \le (\tau^* - t_V) + t_V(2(1+C_f) \|\nabla V - \nabla V^\#\|_\infty)$, hence

$$ (4.19) \qquad t_V (1 - 2(1 + C_f) \|\nabla V - \nabla V^\#\|_\infty) \le \tau^*. $$

By assumption, $1 - 2(1 + C_f) \|\nabla V - \nabla V^\#\|_\infty \ge \tau^*/T > 0$, we can thus divide Eq. (4.19) by this quantity to obtain the result of the corollary: $t_V \le \tau^*/(1 - 2(1 + C_f) \|\nabla V - \nabla V^\#\|_\infty)$. $\square$

In the previous theorem and the corollary, we saw that the suboptimality, in terms of hitting time, of a closed-loop trajectory $(x_V(\cdot), u_V(\cdot))$ satisfying the differential inclusion $(CL_V)$ decreases as the approximation error $\|\nabla V - \nabla V^\#\|_\infty$ decreases. Furthermore, we see that the closed-loop trajectory comes close to optimality when the approximation error vanishes. We now study a sufficient condition under which the approximation $\|\nabla V_d - \nabla V^\#\|_\infty$ can be made arbitrarily small, using a polynomial $V_d(t, x)$ of sufficiently large degree $d \in \mathbb{N}$.

### 4.3. A sufficient regularity condition for the existence of near-optimal controllers based on polynomials.
In the case where the value function is twice differentiable, there exist polynomials $V_d$ with such a vanishing approximation error $\|\nabla V_d - \nabla V^\#\|_\infty$, and that are near optimal solutions in the hierarchy of semi-infinite programs $(D_d)$.

**THEOREM 4.5.** Under Assumptions 1-3, if the value function $V^\#$ belongs to $C^2(\mathbb{R}^p|Y)$, and is a subsolution to the HJB equation, then there exist a sequence of polynomials $(V_d(t, x))_{d \in \mathbb{N}^*}$, with $V_d(t, x) \in \mathbb{R}_d[t, x_1, \dots, x_n]$, and two constants $c_1, c_2 > 0$, such that for all $d \in \mathbb{N}^*$,

* The polynomial $V_d(t, x)$ is feasible, and $c_1/d$-optimal in the problems $(D_F)$ and $(D_d)$,

* The following inequality holds: $\|\nabla V_d - \nabla V^\#\|_\infty \le c_2/d$.

Under these hypotheses, the polynomials $V_d(t, x)$ are subsolutions to the HJB equation, and form a maximizing sequence of the problem $(D_F)$; we also notice that the hierarchy of semi-infinite programs $(D_d)$ converges in $O(1/d)$ in terms of objective value. Moreover, according to Cor. 4.4, for any sequence of closed-loop trajectories $(x_{V_d}(\cdot), u_{V_d}(\cdot))$, the associated hitting times converge to the minimal time $\tau^*$: this is a minimizing sequence of trajectories for the optimal time control problem (2.1)-(2.4).

*Proof.* By definition of $C^2(\mathbb{R}^{n+1}|Y)$, there exists a function $Q \in C^2(\mathbb{R}^{n+1})$ such that $V^\#(y) = Q(y)$ and $\nabla V^\#(y) = \nabla Q(y)$ for all $y \in Y$. In application of Lemma A.2, as $Q$ has a locally Lipschitz gradient since it is twice differentiable, there exists a constant $A > 0$, and a sequence of polynomials $(w_d(t, x))_{d \in \mathbb{N}^*}$ with $w_d(t, x) \in \mathbb{R}_d[t, x_1, \dots, x_n]$ and such that for all $(t, x) \in Y$, $|w_d(t, x) - Q(t, x)| \le \frac{A}{d}$ and $\|\nabla w_d(t, x) - \nabla Q(t, x)\|_2 \le \frac{A}{d}$. With $\alpha_d = \frac{A(1+C_f)}{d}$, and $\beta_d = \frac{A}{d}(1 + T + TC_f)$, we define the polynomial $V_d(t, x) = w_d(t, x) + \alpha_d t - \beta_d \in \mathbb{R}_d[t, x_1, \dots, x_n]$. First, we notice that $\|\nabla V_d - \nabla V^\#\|_\infty \le \|\nabla w_d - \nabla V^\#\|_\infty + \alpha_d \le \frac{A(2+C_f)}{d}$ for all $d \ge 1$. This proves the second point of the theorem, having defined the constant $c_2 = A(2 + C_f)$, which is independent from $d$. We prove now the first point. For all $d \ge 1$, and $(t, x, u) \in [0, T] \times X \times U$,

$$ (4.20) \qquad \begin{aligned} \partial_t V_d(t, x) + 1 + \nabla_x V_d(t, x)^\top f(t, x, u) &= \alpha_d + \partial_t V^\#(t, x) + 1 + \nabla_x V^\#(t, x)^\top f(t, x, u) \\ & \qquad + (\nabla w_d(t, x) - \nabla V^\#(t, x))^\top \begin{pmatrix} 1 \\ f(t, x, u) \end{pmatrix} \\ & \qquad \ge \alpha_d + (\nabla w_d(t, x) - \nabla V^\#(t, x))^\top \begin{pmatrix} 1 \\ f(t, x, u) \end{pmatrix}, \end{aligned} $$

as $V^\#$ is a subsolution to the HJB equation, hence satisfies Eq. (2.7). Using the Cauchy-Schwartz inequality, we obtain $\partial_t V_d(t, x) + 1 + \nabla_x V_d(t, x)^\top f(t, x, u) \ge \alpha_d - \|\nabla w_d(y) - \nabla V^\#(y)\|_2(1 + C_f) \ge \alpha_d - \frac{A}{d}(1 + C_f) = 0$. This
simplex algorithm of the commercial solver Gurobi 10.0 [17]. At each iteration, we add a maximum of 100 points to the set $\mathcal{Y}^k$. The separation oracle (3.5) is implemented with a random sampling scheme (with 500,000 samples at each iteration to detect violated constraints), and with the global optimization solver SCIP 8 [6], for the certification at the last iterate. This solver is used with a relative tolerance $\delta = 10^{-4}$, and with a time limit of 10,000s. We also precise that we compute a heuristic trajectory based on the particularities of each problem; this heuristic is not optimal, but provides an upper bound $T$ on the minimum time, and therefore, a relevant time horizon $[0, T]$. The trajectory resulting from the heuristic is used to initialize the set $\mathcal{Y}^0$, in the sense that we enforce the HJB inequality for some points of this trajectory. During the iterations of the algorithm, we obtain functions $V_{\theta^k}(t, x)$ and we simulate the associate feedback trajectory defined by the differential inclusion $(CL_V)$; if the obtained trajectory reaches the target set, it gives us an upper bound. Those trajectories are also used to enrich the set $\mathcal{Y}^k$. For all the numerical experiments, the regularization parameter is $\mu = 10^{-5}$, and we use the tolerance $\epsilon = 10^{-3}$.

Table 1, Table 3 and Table 4 present the numerical results for three different applications. The different columns of these tables are the following:

* $d \in \mathbb{N}$ is the degree of the polynomial basis used.

* "Estimated value of $(R_d)$" stands for the value $V_{\theta}(0, x_0)$, where $\theta$ is the output of Algorithm 3.1, using the sampling oracle. This estimated value of $(R_d)$ is not an exact lower bound, since this sampling oracle does not provide the guarantee that $\theta$ is indeed feasible in (SIP).

* "Certified lower bound for $(R_d)$" stands for $V_{\theta}(0, x_0) - \hat{\phi}(\theta)(1+T)$, where $\theta$ is as defined above and $\hat{\phi}(\theta)$ is a guaranteed upper bound on $\phi(\theta)$, the feasibility error of $\theta$ in (SIP), computed by the global optimization solver SCIP 8. As $V_{\theta}(t, x) + \hat{\phi}(\theta)(t-1-T)$ is therefore feasible in $(R_d)$, the value $V_{\theta}(0, x_0) - \hat{\phi}(\theta)(1+T)$ is a guaranteed lower bound on $\mathrm{val}(R_d)$, and, therefore, on $V^*(0, x_0)$.

* "Value feedback control $(CL_V)$" is the hitting time of the best feasible control generated along the iterations: either with the heuristic control at the first iteration, or the closed-loop controlled trajectory defined by $(CL_V)$ associated with $V = V_{\theta^k}$ at iteration $k$ of Algorithm 3.1.

* "Solution time (in $s$)" is the total computational time of the heuristic control, of the iterations of Algorithm 3.1 including the sampling oracle, and of the closed-loop trajectory simulation. Therefore, this is the computational time needed to obtain the estimated value of $(R_d)$ (second column), and the best feasible control (fourth column).

* "Iterations number" is the total number of iterations of Algorithm 3.1.

* "Certification time (in $s$)" is the computational time of the global optimization solver SCIP 8, playing the role of $\delta$-oracle, to compute the aforementioned bound $\hat{\phi}(\theta)$, and deduce the certified lower bound (third column).

**5.1. A time-dependent Zermelo problem.** We consider a time-dependent nonlinear system with $n = 2$ and $m = 2$, defined by

$$ (5.1) \qquad \dot{x}_1(t) = u_1(t) + \frac{1}{2}(1+t) \sin(\pi x_2(t)) $$

$$ (5.2) \qquad \dot{x}_2(t) = u_2(t), $$

with the state constraint set $X = [-1, 1] \times [-1, 0]$, the control set $U = B(0, 1)$. This is the celebrated Zermelo problem, but with a river flow gaining in intensity over time. Fig. 1 gives a representation of this flow. The initial condition is $x(0) = (0, -1)$, and the target set is $K = B(0, r)$, for $r = 0.05$. The travel time associated with the heuristic control, consisting in following a straight trajectory, is 1.261 (see Fig 2). Table 1 presents the numerical results for different values of $d$. We see that the value of the linear semi-infinite program $(R_d)$ quickly converges as $d$ increases: starting from $d = 6$, the 4 first digits of the estimated value (second column) reach a plateau which corresponds to the value (1.100) of the best feasible trajectory we generate with our feedback control. As regards the certified lower bound, the best value (1.092) is obtained for $d = 5$. For greater $d$, we see that increasing $d$ deteriorates the tightness of the best certified bound. This is due to the fact that the separation problem becomes more difficult, with two consequences: (i) the sampling fails to detect unsatisfied constraints, so Algorithm 3.1 stops with a solution that has a real infeasibility $\phi(\theta)$ larger than $\epsilon$ (targeted tolerance), and (ii) the global optimization solver called afterwards does not manage to solve the separation problem to global optimality within the time limit (case $d \in \{7, 8\}$), given only a large upper bound $\hat{\phi}(\theta)$ on the true infeasibility $\phi(\theta)$. We notice that as soon as $d \ge 3$, the feedback control defined by $(CL_V)$ (see Sect. 4) yields a trajectory that is 13% faster than the heuristic trajectory. In summary, we obtain a certified optimization gap of 0.7% for this minimal time control problem.
proves that $V_d$ satisfies Eq. (2.7). It also satisfies Eq. (2.8), because for any $(t, x) \in [0, T] \times K$,

$$ (4.21) \qquad \begin{aligned} V_d(t,x) &= w_d(t,x) + \alpha_d t - \beta_d \\ &\le V^\#(t,x) + \frac{A}{d} + \alpha_d t - \beta_d \\ &\le V^\#(t,x) + \frac{A}{d} + \alpha_d T - \beta_d \\ &\le V^\#(t,x) = 0. \end{aligned} $$

since $\frac{A}{d} + \alpha_d T - \beta_d = 0$ by definition of $\beta_d$, and since $x \in K$. We deduce that $V_d$ is feasible in $(D_d)$. Its objective value is $V_d(0, x_0) \ge w_d(0, x_0) - \beta_d \ge V^\#(0, x_0) - \frac{c_1}{d} - \beta_d = V^\#(0, x_0) - \frac{c_1}{d}$, where $c_1 = A(2 + T + TC_f)$. As $V^\#(0, x_0) = V^*(0, x_0)$ due to Assumption 2, $V_d(0, x_0) \ge V^*(0, x_0) - \frac{c_1}{d} \ge \text{val}(D_\mathcal{F}) - \frac{c_1}{d} \ge \text{val}(D_d) - \frac{c_1}{d}$, and we therefore conclude that $V_d$ is $\frac{c_1}{d}$-optimal in $(D_\mathcal{F})$ and $(D_d)$. $\square$

REMARK 3. Admittedly, the hypothesis in Th. 4.5 that the value function $V^\#$ belongs to $C^2(\mathbb{R}^p|\Sigma)$ is stringent. It is worth noting, however, that there exist systems that satisfy this hypothesis. Here is an example: $\dot{x}(t) = u(t)$, $x(t) \in X = [0, 1]^2$, $\|u(t)\| \le 1$ and $K = \{0\} \times [0, 1]$. The value function associated with the horizon $T = \infty$ is $V^\#(t, x) = x_1$.

5. **Illustrative examples.** We implemented and tested the proposed methodology on three Minimal Time Control Problems: a generalization of the Zermelo problem, a regatta problem and a generalization of the Brockett integrator. The numerical examples in this section were processed with our Julia [7] package `MinTimeControl.jl`<fnref n="1" />. In this implementation of Algorithm 3.1, the master problem (3.6) is solved with the simplex algorithm of the commercial solver `Gurobi 10.0` [23]. At each iteration, we add a maximum of 100 points to the set $\mathcal{Y}^k$. The separation oracle (3.5) is implemented with a random sampling scheme (with 500 000 samples at each iteration to detect violated constraints), and with the global optimization solver `SCIP 8` [6], for the certification at the last iterate. This solver is used with a relative tolerance $\delta = 10^{-4}$, and with a time limit of 10 000 s. We also precise that we compute a heuristic trajectory based on the particularities of each problem; this heuristic is not optimal, but provides an upper bound $T$ on the minimum time, and therefore, a relevant time horizon $[0, T]$. The trajectory resulting from the heuristic is used to initialize the set $\mathcal{Y}^0$, in the sense that we enforce the HJB inequality for some points of this trajectory. During the iterations of the algorithm, we obtain functions $V_{\theta^k}(t, x)$ and we simulate the associate feedback trajectory defined by the differential inclusion $(CL_V)$; if the obtained trajectory reaches the target set, it gives us an upper bound. Those trajectories are also used to enrich the set $\mathcal{Y}^k$. For all the numerical experiments, the regularization parameter is $\mu = 10^{-5}$, and we use the tolerance $\epsilon = 10^{-3}$.

Tables 1, 3 and 4 present the numerical results for three different applications. The different columns of these tables are the following:

* “$d \in \mathbb{N}$” is the degree of the polynomial basis used.

* “Estimated value of $(D_d)$” stands for the value $V_\theta(0, x_0)$, where $\theta$ is the output of Algorithm 3.1, using the sampling oracle. This estimated value of $(D_d)$ is not an exact lower bound, since this sampling oracle does not provide the guarantee that $\theta$ is indeed feasible in (SIP).

* “Certified lower bound for $(D_d)$” stands for $V_\theta(0, x_0) - \hat{\phi}(\theta)(1+T)$, where $\theta$ is as defined above and $\hat{\phi}(\theta)$ is a guaranteed upper bound on $\phi(\theta)$, the feasibility error of $\theta$ in (SIP), computed by the global optimization solver `SCIP 8`. As $V_\theta(t, x) + \hat{\phi}(\theta)(t-1-T)$ is therefore feasible in $(D_d)$, the value $V_\theta(0, x_0) - \hat{\phi}(\theta)(1+T)$ is a guaranteed lower bound on $\text{val}(D_d)$, and, therefore, on $V^*(0, x_0)$.

* “Value feedback control $(CL_V)$” is the hitting time of the best feasible control generated along the iterations: either with the heuristic control at the first iteration, or the closed-loop controlled trajectory defined by $(CL_V)$ associated with $V = V_{\theta^k}$ at iteration $k$ of Algorithm 3.1.

* “Solution time (in s)” is the total computational time of the heuristic control, of the iterations of Algorithm 3.1 including the sampling oracle, and of the closed-loop trajectory simulation. Therefore, this is the computational time needed to obtain the estimated value of $(D_d)$ (second column), and the best feasible control (fourth column).

* “Iterations number” is the total number of iterations of Algorithm 3.1.

* “Certification time (in s)” is the computational time of the global optimization solver `SCIP 8`, playing the role of $\delta$-oracle, to compute the aforementioned bound $\hat{\phi}(\theta)$, and deduce the certified lower bound

<footnote n="1">¹This package is available at github.com/aoustry/MinTimeControl.jl</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1: Representation of the water flow in the Zermelo problem</figcaption></figure>

<table><caption>Table 1: Time-dependent Zermelo problem: lower and upper bounds, and computational times for various degrees of the SIP hierarchy (R<sub>d</sub>)</caption><thead><tr><th>d ∈ N</th><th>Estimated value of (R<sub>d</sub>)</th><th>Certified LB for (R<sub>d</sub>)</th><th>Value feedback control (CL<sub>V</sub>)</th><th>Solution time (in s)</th><th>Iterations number</th><th>Certification time (in s)</th></tr></thead><tbody><tr><td>2</td><td>0.952</td><td>0.945</td><td>1.261</td><td>2</td><td>4</td><td>1</td></tr><tr><td>3</td><td>1.064</td><td>1.044</td><td>1.101</td><td>12</td><td>14</td><td>12</td></tr><tr><td>4</td><td>1.096</td><td>1.084</td><td>1.100</td><td>22</td><td>17</td><td>1530</td></tr><tr><td>5</td><td>1.099</td><td>1.092</td><td>1.100</td><td>54</td><td>22</td><td>4000</td></tr><tr><td>6</td><td>1.100</td><td>1.059</td><td>1.100</td><td>60</td><td>18</td><td>2330</td></tr><tr><td>7</td><td>1.100</td><td>1.051</td><td>1.100</td><td>105</td><td>24</td><td>TL</td></tr><tr><td>8</td><td>1.100</td><td>0.690</td><td>1.100</td><td>215</td><td>31</td><td>TL</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 2: Time-dependent Zermelo problem: heuristic control and feedback control (d = 6)</figcaption></figure>
(third column).

**5.1. A time-dependent Zermelo problem.** We consider a time-dependent nonlinear system with $n = 2$ and $m = 2$, defined by

$$ (5.1) \qquad \dot{x}_1(t) = u_1(t) + \frac{1}{2}(1+t) \sin(\pi x_2(t)) $$

$$ (5.2) \qquad \dot{x}_2(t) = u_2(t), $$

with the state constraint set $X = [-1, 1] \times [-1, 0]$, the control set $U = B(0, 1)$. This is the celebrated Zermelo problem, but with a river flow gaining in intensity over time. Fig. 1 gives a representation of this flow. The initial condition is $x(0) = (0, -1)$, and the target set is $K = B(0, r)$, for $r = 0.05$. The travel time associated with the heuristic control, consisting in following a straight trajectory, is 1.261 (see Fig 2). Table 1 presents the numerical results for different values of $d$. We see that the value of the linear semi-infinite program ($D_d$) quickly converges as $d$ increases: starting from $d = 6$, the 4 first digits of the estimated value (second column) reach a plateau which corresponds to the value (1.100) of the best feasible trajectory we generate with our feedback control. As regards the certified lower bound, the best value (1.092) is obtained for $d = 5$. For greater $d$, we see that increasing $d$ deteriorates the tightness of the best certified bound. This is due to the fact that the separation problem becomes more difficult, with two consequences: (i) the sampling fails to detect unsatisfied constraints, so Algorithm 3.1 stops with a solution that has a real infeasibility $\phi(\theta)$ larger than $\epsilon$ (targeted tolerance), and (ii) the global optimization solver called afterwards does not manage to solve the separation problem to global optimality within the time limit (case $d \in \{7, 8\}$), given only a large upper bound $\dot{\phi}(\theta)$ on the true infeasibility $\phi(\theta)$. We notice that as soon as $d \ge 3$, the feedback control defined by ($CL_V$) (see Sect. 4) yields a trajectory that is 13% faster than the heuristic trajectory. In summary, we obtain a certified optimization gap of 0.7% for this minimal time control problem.

<figure><img src="image_3.png" /><figcaption>Fig. 1: Representation of the water flow in the Zermelo problem</figcaption></figure>

In the special case of this non-polynomial controlled system, a polynomial reformulation exists, at the price of increasing the dimension of the system to $n = 4$:

$$ (5.3) \qquad \dot{x}_1(t) = u_1(t) + \frac{1}{2}(1+t) x_4(t) $$

$$ (5.4) \qquad \dot{x}_2(t) = u_2(t) $$

$$ (5.5) \qquad \dot{x}_3(t) = -\pi x_4(t) u_2(t) $$

$$ (5.6) \qquad \dot{x}_4(t) = \pi x_3(t) u_2(t), $$

with the state constraint set $\hat{X} = [-1, 1] \times [-1, 0] \times [-1, 1] \times [-1, 0]$, the control set $\hat{U} = B(0, 1)$, the terminal set $\hat{K} = B(0, r) \times [-1, 1] \times [-1, 0]$, and the initial condition $x_0 = (0, -1, -1, 0)$. The dynamics maintain the equalities
In the special case of this non-polynomial controlled system, a polynomial reformulation exists, at the price of increasing the dimension of the system to $n = 4$:

$$ (5.3) \qquad \dot{x}_1(t) = u_1(t) + \frac{1}{2}(1+t) \sin(\pi x_2(t)) $$

$$ (5.4) \qquad \dot{x}_2(t) = u_2(t) $$

$$ (5.5) \qquad \dot{x}_3(t) = -\pi x_4(t)u_2(t) $$

$$ (5.6) \qquad \dot{x}_4(t) = \pi x_3(t)u_2(t), $$

with the state constraint set $\hat{X} = [-1, 1] \times [-1, 0] \times [-1, 1] \times [-1, 0]$, the control set $\hat{U} = B(0, 1)$, the terminal set $\hat{K} = B(0, r) \times [-1, 1] \times [-1, 0]$, and the initial condition $x_0 = (0, -1, -1, 0)$. The dynamics maintain the equalities $x_3(t) = \cos(\pi x_2(t))$ and $x_4(t) = \sin(\pi x_2(t))$. We are therefore able to compare our approach with the sum-of-squares (SOS) hierarchy, which consists of replacing SIP inequalities in $(R_d)$ with SOS positivity certificates. For each order $k$ of the hierarchy, i.e., for a maximal degree $d = 2k$ of the polynomial basis, this yields a semi-definite programming problem that we solve with the solver CSDP, used with the package SumOfSquares.jl. We obtain a polynomial $V(t, x_1, x_2, x_3, x_4)$ that is solution of the corresponding relaxation. Based on this polynomial, we can also generate a feedback controlled trajectory solution of the differential inclusion $(CL_V)$. Table 2 compares

<table><caption>Table 2: Time-dependent Zermelo problem: comparing the SIP and the SOS approaches</caption><thead><tr><th rowspan="2">d ∈ N</th><th colspan="3">SIP hierarchy</th><th colspan="3">SOS hierarchy</th></tr><tr><th>Est./Cert. LB</th><th>Val. feedback control (CL<sub>V</sub>)</th><th>Sol./Cert. time (in s)</th><th>Cert. LB</th><th>Val. feedback control (CL<sub>V</sub>)</th><th>Sol. time (in s)</th></tr></thead><tbody><tr><td>2</td><td>0.952/0.945</td><td>1.261</td><td>2/1</td><td>0.533</td><td>1.261</td><td>≤ 1</td></tr><tr><td>4</td><td>1.096/1.084</td><td>1.101</td><td>22/1530</td><td>1.064</td><td>1.105</td><td>1</td></tr><tr><td>6</td><td>1.100/1.059</td><td>1.100</td><td>60/2330</td><td>1.099</td><td>1.100</td><td>12</td></tr><tr><td>8</td><td>1.100/0.690</td><td>1.100</td><td>215/TL</td><td>1.100</td><td>1.100</td><td>190</td></tr></tbody></table>

the performance of the SIP and the SOS approaches. We see that for low-degree polynomials ($d \le 4$), the semi-infinite hierarchy gives better lower bounds than the SOS hierarchy, although at a higher computational time in the case $d = 4$. For $d \in \{6, 8\}$, the lower bound of the SOS hierarchy is tight, while only the estimated lower-bound of the SIP hierarchy is tight: to obtain a certified lower bound, the SOS hierarchy performs better. For these values of the degree $d$, this optional certification step (calling to the global optimization solver) is costly in the proposed approach. For this first example, where the SOS hierarchy is applicable since a polynomial reformulation of the dynamical system exists, the SIP approach is slower than the SOS hierarchy.

**5.2. A regatta toy-model.** We consider a time-dependent nonlinear (and non-polynomial) system with $n = 2$ and $m = 1$, defined by

$$ (5.7) \qquad \dot{x}_1(t) = \text{windspeed}(t) \text{ polar}[u(t)] \cos(u(t) + \text{windangle}(t)) $$

$$ (5.8) \qquad \dot{x}_2(t) = \text{windspeed}(t) \text{ polar}[u(t)] \sin(u(t) + \text{windangle}(t)), $$

where $\text{windspeed}(t) = 2 + t$, $\text{windangle}(t) = \frac{\pi}{2}(1 - 0.4t)$ and $\text{polar}[u] = |\sin(\frac{2u}{3})|$. In this model, the control $u(t)$ represents the relative angle between the heading of the boat and the (origin) direction of the wind. The evolution of the wind direction over time is depicted in Fig. 3. The polar curve of this toy model of a sailing boat is represented in Fig. 4; this figure clearly shows that this model does not satisfy Assumption 1. Although the absence of duality gap between the control problem and the LP problem $(D_F)$ is, therefore, not guaranteed, we see in Table 3 that if this gap exists in this case, it is low (below 1.6%). The state constraint set is $X = [-1, 1]^2$, and the control set $U = [-\pi, \pi]$. The initial condition is $x(0) = (0, -1)$, and the target set is $K = B(0, r)$, for $r = 0.05$. The travel time associated with the heuristic control, consisting in following a straight trajectory, is 1.278 (see Fig. 5).

We see that the highest estimated value of $(R_d)$, for $d = 7$ and $d = 8$, is 0.5% lower than the value (0.913) of the best feasible trajectory obtained with the feedback controller for $d = 6$. This feedback controller yields a trajectory which is 29% faster than the heuristic trajectory. As regards the certified lower bound, $d = 4$ yields the best result (0.896), at a price of a running time of 498s for the exact oracle (SCIP 8). For the same reasons
<table><caption>Table 1: Time-dependent Zermelo problem: lower and upper bounds, and computational times for various degrees of the SIP hierarchy ($D_d$)</caption><thead><tr><th>d &isin; N</th><th>Estimated value of (D<sub>d</sub>)</th><th>Certified LB for (D<sub>d</sub>)</th><th>Value feedback control (CL<sub>V</sub>)</th><th>Solution time (in s)</th><th>Iterations number</th><th>Certification time (in s)</th></tr></thead><tbody><tr><td>2</td><td>0.952</td><td>0.945</td><td>1.261</td><td>2</td><td>4</td><td>1</td></tr><tr><td>3</td><td>1.064</td><td>1.044</td><td>1.101</td><td>12</td><td>14</td><td>12</td></tr><tr><td>4</td><td>1.096</td><td>1.084</td><td>1.100</td><td>22</td><td>17</td><td>1 530</td></tr><tr><td>5</td><td>1.099</td><td>1.092</td><td>1.100</td><td>54</td><td>22</td><td>4 000</td></tr><tr><td>6</td><td>1.100</td><td>1.059</td><td>1.100</td><td>60</td><td>18</td><td>2 330</td></tr><tr><td>7</td><td>1.100</td><td>1.051</td><td>1.100</td><td>105</td><td>24</td><td>TL</td></tr><tr><td>8</td><td>1.100</td><td>0.690</td><td>1.100</td><td>215</td><td>31</td><td>TL</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 2: Time-dependent Zermelo problem: heuristic control and feedback control ($d = 6$)</figcaption></figure>

$x_3(t) = \cos(\pi x_2(t))$ and $x_4(t) = \sin(\pi x_2(t))$. We are therefore able to compare our approach with the sum-of-squares (SOS) hierarchy, which consists of replacing SIP inequalities in ($D_d$) with SOS positivity certificates. For each order $k$ of the hierarchy, i.e., for a maximal degree $d = 2k$ of the polynomial basis, this yields a semi-definite programming problem that we solve with the solver CSDP, used with the package `SumOfSquares.jl`. We obtain a polynomial $V(t, x_1, x_2, x_3, x_4)$ that is solution of the corresponding relaxation. Based on this polynomial, we can also generate a feedback controlled trajectory solution of the differential inclusion ($CL_V$). Table 2 compares

<table><caption>Table 2: Time-dependent Zermelo problem: comparing the SIP and the SOS approaches</caption><thead><tr><th rowspan="2">d &isin; N</th><th colspan="3">SIP hierarchy</th><th colspan="3">SOS hierarchy</th></tr><tr><th>Est./Cert. LB</th><th>Val. feedback control (CL<sub>V</sub>)</th><th>Sol./Cert. time (in s)</th><th>Cert. LB</th><th>Val. feedback control (CL<sub>V</sub>)</th><th>Sol. time (in s)</th></tr></thead><tbody><tr><td>2</td><td>0.952/0.945</td><td>1.261</td><td>2/1</td><td>0.533</td><td>1.261</td><td>&le; 1</td></tr><tr><td>4</td><td>1.096/1.084</td><td>1.101</td><td>22/1530</td><td>1.064</td><td>1.105</td><td>1</td></tr><tr><td>6</td><td>1.100/1.059</td><td>1.100</td><td>60/2330</td><td>1.099</td><td>1.100</td><td>12</td></tr><tr><td>8</td><td>1.100/0.690</td><td>1.100</td><td>215/TL</td><td>1.100</td><td>1.100</td><td>190</td></tr></tbody></table>

the performance of the SIP and the SOS approaches. We see that for low-degree polynomials ($d \le 4$), the semi-infinite hierarchy gives better lower bounds than the SOS hierarchy, although at a higher computational time in the case $d = 4$. For $d \in \{6, 8\}$, the lower bound of the SOS hierarchy is tight, while only the estimated lower-
<figure><img src="image_5.png" /><figcaption>Fig. 3: Regatta problem: wind direction at different times</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4: Regatta problem: the polar curve of the sailing boat</figcaption></figure>

<table><caption>Table 3: Regatta problem: lower and upper bounds, and computational times for various degrees of the SIP hierarchy (R<sub>d</sub>)</caption><thead><tr><th>d ∈ N</th><th>Estimated value of (R<sub>d</sub>)</th><th>Certified LB for (R<sub>d</sub>)</th><th>Value feedback control (CL<sub>V</sub>)</th><th>Solution time (in s)</th><th>Iterations number</th><th>Certification time (in s)</th></tr></thead><tbody><tr><td>2</td><td>0.834</td><td>0.829</td><td>1.278</td><td>6</td><td>6.0</td><td>2</td></tr><tr><td>3</td><td>0.896</td><td>0.880</td><td>0.912</td><td>16</td><td>10.0</td><td>56</td></tr><tr><td>4</td><td>0.904</td><td>0.896</td><td>0.915</td><td>31</td><td>13.0</td><td>498</td></tr><tr><td>5</td><td>0.907</td><td>0.774</td><td>0.912</td><td>52</td><td>16.0</td><td>1020</td></tr><tr><td>6</td><td>0.907</td><td>0.799</td><td>0.912</td><td>100</td><td>22.0</td><td>1930</td></tr><tr><td>7</td><td>0.908</td><td>0.691</td><td>0.912</td><td>190</td><td>29.0</td><td>7600</td></tr><tr><td>8</td><td>0.908</td><td>0.000</td><td>0.911</td><td>312</td><td>33.0</td><td>TL</td></tr></tbody></table>

as in the previous application, a larger *d* does not necessarily mean a better certified lower bound obtained within the time limit. In summary, we obtain a certified optimization gap of 1.6% for this minimal time control problem.

**5.3. A generalized Brockett integrator.** For $n \in \mathbb{N}^*$ and $m = n - 1$ and given a continuous mapping $q: \mathbb{R}^n \to \mathbb{R}^m$, we consider the following generalization of the Brockett integrator [29],

$$ (5.9) \qquad \dot{x}_i(t) = u_i(t) \quad \forall i \in \{1, \dots, n-1\} $$

$$ (5.10) \qquad \dot{x}_n(t) = q(x(t))^{\top} u(t). $$
bound of the SIP hierarchy is tight: to obtain a certified lower bound, the SOS hierarchy performs better. For these values of the degree $d$, this optional certification step (calling to the global optimization solver) is costly in the proposed approach. For this first example, where the SOS hierarchy is applicable since a polynomial reformulation of the dynamical system exists, the SIP approach is slower than the SOS hierarchy.

**5.2. A regatta toy-model.** We consider a time-dependent nonlinear (and non-polynomial) system with $n = 2$ and $m = 1$, defined by

$$ (5.7) \qquad \dot{x}_1(t) = \text{windspeed}(t) \text{ polar}[u(t)] \cos(u(t) + \text{windangle}(t)) $$

$$ (5.8) \qquad \dot{x}_2(t) = \text{windspeed}(t) \text{ polar}[u(t)] \sin(u(t) + \text{windangle}(t)), $$

where $\text{windspeed}(t) = 2 + t$, $\text{windangle}(t) = \frac{\pi}{2}(1 - 0.4t)$ and $\text{polar}[u] = |\sin(\frac{2u}{3})|$. In this model, the control $u(t)$ represents the relative angle between the heading of the boat and the (origin) direction of the wind. The evolution of the wind direction over time is depicted in Fig. 3. The polar curve of this toy model of a sailing boat is represented in Fig 4; this figure clearly shows that this model does not satisfy Assumption 1. Although the absence of duality gap between the control problem and the LP problem ($D_\mathcal{F}$) is, therefore, not guaranteed, we see in Table 3 that if this gap exists in this case, it is low (below 1.6%). The state constraint set is $X = [-1, 1]^2$, and the control set $U = [-\pi, \pi]$. The initial condition is $x(0) = (0, -1)$, and the target set is $K = B(0, r)$, for $r = 0.05$. The travel time associated with the heuristic control, consisting in following a straight trajectory, is 1.278 (see Fig. 5).

<figure><img src="image_5.png" /><figcaption>Fig. 3: Regatta problem: wind direction at different times</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4: Regatta problem: the polar curve of the sailing boat</figcaption></figure>

We see that the highest estimated value of ($D_d$), for $d = 7$ and $d = 8$, is 0.5% lower than the value (0.913) of the best feasible trajectory obtained with the feedback controller for $d = 6$. This feedback controller yields a
<figure><img src="image_7.png" /><figcaption>Fig. 5: Regatta problem: heuristic control and feedback control ($d = 6$)</figcaption></figure>

In particular, we study this system for $n = 6$, and $q(x) = (2/(2 + x_4), -x_1, -\cos(x_1x_3), \exp(x_2), x_1x_2x_6)$. The state constraint set is $X = [-1, 1]^5$, and the control set is $U = B(0, 1)$. The initial condition is $(x_0) = \frac{1}{2}1$, and the target set is $K = B(0, r)$, for $r = 0.05$. The travel time associated with the heuristic control is 1.377.

<table><caption>Table 4: Generalized Brockett integrator: lower and upper bounds, and computational times for various degrees of the SIP hierarchy ($R_d$)</caption><thead><tr><th>d ∈ N</th><th>Estimated value of (R<sub>d</sub>)</th><th>Certified LB for (R<sub>d</sub>)</th><th>Value feedback control (CL<sub>V</sub>)</th><th>Solution time (in s)</th><th>Iterations number</th><th>Certification time (in s)</th></tr></thead><tbody><tr><td>2</td><td>1.071</td><td>0.763</td><td>1.071</td><td>220</td><td>28</td><td>73</td></tr><tr><td>3</td><td>1.072</td><td>0.000</td><td>1.071</td><td>5630</td><td>133</td><td>TL</td></tr><tr><td>4</td><td>1.072</td><td>0.000</td><td>1.070</td><td>147400</td><td>319</td><td>TL</td></tr></tbody></table>

Since this system has a larger dimension than the other two examples, we see that the computation times are longer for the same degree $d$. Already for $d = 2$, we obtain an estimated value of ($R_d$) that is within 0.1% of the value of the feedback control (1.070). This feedback control yields an improvement of 22% over the heuristic trajectory. Note that the estimated values of ($R_d$) computed by Algorithm 3.1 with the (inexact) sampling oracle are slightly larger than the value of the best trajectory we computed: thus, these estimates are not valid lower bounds, but only estimates of the value of the minimum time control problem. Regarding the certification of lower bounds, the global optimization solver `SCIP 8` fails to produce tight upper and lower bounds on $\phi(\theta)$, the infeasibility of the solution $\theta$ returned by Algorithm 3.1. Therefore, the resulting certified lower bounds are not tight either. In summary, we obtain a certified optimization gap of 29% for this minimal time control problem.

**6. Discussion.** We apply the dual approach in minimal time control, that consists in searching for maximal subsolutions of the HJB equation, to generic nonlinear, even non-polynomial, controlled systems. The basis functions used to generate these subsolutions are polynomials, that are subject to semi-infinite constraints. We prove the theoretical convergence of the resulting hierarchy of semi-infinite linear programs, and our numerical tests on three different systems show good convergence properties in practice. These results show that the use of a random sampling oracle allows a good approximation of the value of the control problem. For small systems, it is even possible to obtain tight and certified lower bounds, based on a global optimization solver. Finally,
<table><caption>Table 3: Regatta problem: lower and upper bounds, and computational times for various degrees of the SIP hierarchy ($D_d$)</caption><thead><tr><th>d &isin; N</th><th>Estimated value of (D<sub>d</sub>)</th><th>Certified LB for (D<sub>d</sub>)</th><th>Value feedback control (CL<sub>V</sub>)</th><th>Solution time (in s)</th><th>Iterations number</th><th>Certification time (in s)</th></tr></thead><tbody><tr><td>2</td><td>0.834</td><td>0.829</td><td>1.278</td><td>6</td><td>6.0</td><td>2</td></tr><tr><td>3</td><td>0.896</td><td>0.880</td><td>0.912</td><td>16</td><td>10.0</td><td>56</td></tr><tr><td>4</td><td>0.904</td><td>0.896</td><td>0.915</td><td>31</td><td>13.0</td><td>498</td></tr><tr><td>5</td><td>0.907</td><td>0.774</td><td>0.912</td><td>52</td><td>16.0</td><td>1 020</td></tr><tr><td>6</td><td>0.907</td><td>0.799</td><td>0.912</td><td>100</td><td>22.0</td><td>1 930</td></tr><tr><td>7</td><td>0.908</td><td>0.691</td><td>0.912</td><td>190</td><td>29.0</td><td>7 600</td></tr><tr><td>8</td><td>0.908</td><td>0.000</td><td>0.911</td><td>312</td><td>33.0</td><td>TL</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 5: Regatta problem: heuristic control and feedback control ($d = 6$)</figcaption></figure>

trajectory which is 29% faster than the heuristic trajectory. As regards the certified lower bound, $d = 4$ yields the best result (0.896), at a price of a running time of 498s for the exact oracle (SCIP 8). For the same reasons as in the previous application, a larger $d$ does not necessarily mean a better certified lower bound obtained within the time limit. In summary, we obtain a certified optimization gap of 1.6% for this minimal time control problem.

**5.3. A generalized Brockett integrator.** For $n \in \mathbb{N}^*$ and $m = n - 1$ and given a continuous mapping $q: \mathbb{R}^n \to \mathbb{R}^m$, we consider the following generalization of the Brockett integrator [37],

$$ (5.9) \qquad \dot{x}_i(t) = u_i(t) \quad \forall i \in \{1, \dots, n-1\} $$

$$ (5.10) \qquad \dot{x}_n(t) = q(x(t))^{\top} u(t). $$

In particular, we study this system for $n = 6$, and $q(x) = (2/(2 + x_4), -x_1, -\cos(x_1x_3), \exp(x_2), x_1x_2x_6)$. The state constraint set is $X = [-1, 1]^5$, and the control set is $U = B(0, 1)$. The initial condition is $(x_0) = \frac{1}{2}\mathbf{1}$, and the target set is $K = B(0, r)$, for $r = 0.05$. The travel time associated with the heuristic control is 1.377.

Since this system has a larger dimension than the other two examples, we see that the computation times are longer for the same degree $d$. Already for $d = 2$, we obtain an estimated value of ($D_d$) that is within 0.1% of the value of the feedback control (1.070). This feedback control yields an improvement of 22% over the heuristic
the numerical experiments also show that the computed subsolutions of the HJB equations help to recover
near-optimal controls in a closed-loop form. As illustrated in these numerical experiments, the advantage of
our approach based on semi-infinite programming, compared to the sum-of-squares approach, is the ability to
handle non-polynomial systems. In the numerical example where a polynomial reformulation of the system was
possible, the sum-of-squares approach was, however, faster.

A promising avenue for continuing this work is to investigate the use of a other basis of functions to search
for an approximate value function, resulting in other semi-infinite programming hierarchies with convergence
guarantees. In particular, it would be relevant to use non-differentiable functions in the basis to improve
approximation capabilities for non-differentiable value functions. Another avenue of research is to extend the
approach and theoretical results to a generic optimal control problem.

**Acknowledgments.** The authors would like to thank Maxime Dupuy, Leo Liberti and Claudia D'Ambrosio for fruitful discussions and advice.

REFERENCES

[1] Jean-Pierre Aubin. *Viability Theory*. Systems & control: foundations & applications. Springer Science+Business Media, 1991.

[2] Jean-Pierre Aubin and Hélène Frankowska. *Set-Valued Analysis*. Birkhäuser, 2009.

[3] T. Bagby, L. Bos, and N. Levenberg. Multivariate simultaneous approximation. *Constructive Approximation*, 18(4):569–577, 2002.

[4] Richard Bellman. Dynamic programming. *Science*, 153(3731):34–37, 1966. Publisher: American Association for the Advancement of Science.

[5] Eloïse Berthier, Justin Carpentier, Alessandro Rudi, and Francis Bach. Infinite-dimensional sums-of-squares for optimal control. In *2022 IEEE 61st Conference on Decision and Control (CDC)*, pages 577–582. IEEE, 2022.

[6] Ksenia Bestuzheva, Mathieu Besançon, Wei-Kun Chen, Antonia Chmiela, Tim Donkiewicz, Jasper van Doornmalen, Leon Eifler, Oliver Gaul, Gerald Gamrath, Ambros Gleixner, Leona Gottwald, Christoph Graczyk, Katrin Halbig, Alexander Hoen, Christopher Hojny, Rolf van der Hulst, Thorsten Koch, Marco Lübbecke, Stephen J. Maher, Frederic Matter, Erik Mühmer, Benjamin Müller, Marc E. Pfetsch, Daniel Rehfeldt, Steffan Schlein, Franziska Schlösser, Felipe Serrano, Yuji Shinano, Boro Sofranac, Mark Turner, Stefan Vigerske, Fabian Wegscheider, Philipp Wellner, Dieter Weninger, and Jakob Witzig. The SCIP optimization suite 8.0. Technical report, Optimization Online, 2021.

[7] Jérôme Bolte and Edouard Pauwels. Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning. *Mathematical Programming*, 188(1):19–51, 2021.

[8] Loïc Bourdin and Emmanuel Trélat. Pontryagin maximum principle for optimal sampled-data control problems. *IFAC-PapersOnLine*, 48(25):80–84, 2015.

[9] Eduardo F. Camacho and Carlos Bordons Alba. *Model predictive control*. Springer science & business media, 2013.

[10] I. Capuzzo-Dolcetta and P.-L. Lions. Hamilton-jacobi equations with state constraints. *Transactions of the American Mathematical Society*, 318(2):643–683, 1990.

[11] Martina Cerulli, Antoine Oustry, Claudia D'Ambrosio, and Leo Liberti. Convergent algorithms for a class of convex semi-infinite programs. *SIAM Journal on Optimization*, 32(4):2493–2526, 2022. Publisher: Society for Industrial and Applied Mathematics.

[12] Frank H Clarke. Generalized gradients and applications. *Transactions of the American Mathematical Society*, 205:247–262, 1975.

[13] Frank H. Clarke and Richard B. Vinter. The relationship between the maximum principle and dynamic programming. *SIAM Journal on Control and Optimization*, 25(5):1291–1311, 1987. Publisher: Society for Industrial and Applied Mathematics.

[14] Michael G. Crandall and Pierre-Louis Lions. Viscosity solutions of hamilton-jacobi equations. *Transactions of the American Mathematical Society*, 277(1):1–42, 1983. Publisher: American Mathematical Society.

[15] Halina Frankowska. Optimal trajectories associated with a solution of the contingent hamilton-jacobi equation. *Applied Mathematics and Optimization*, 19(1):291–311, 1989.

[16] Lawrence Craig Evans Gariepy, Ronald F. *Measure Theory and Fine Properties of Functions, Revised Edition*. Chapman and Hall/CRC, 2015.

[17] LLC Gurobi Optimization. Gurobi optimizer reference manual, 2021.

[18] Didier Henrion, Jean B. Lasserre, and Carlo Savorgnan. Nonlinear optimal control synthesis via occupation measures. In *2008 47th IEEE Conference on Decision and Control*, pages 4749–4754, 2008. ISSN: 0191-2216.

[19] Didier Henrion and Edouard Pauwels. Linear conic optimization for nonlinear optimal control, July 2014. arXiv:1407.1650 [math].

[20] Daniel Hernández-Hernández, Onésimo Hernández-Lerma, and Michael Taksar. The linear programming approach to deterministic optimal control problems. *Applicationes Mathematicae*, 24(1):17–33, 1996.

[21] Lars Hörmander. *The Analysis of Linear Partial Differential Operators I*. Classics in Mathematics. Springer, 2003.

[22] Reza Jazar. *Theory of Applied Robotics: Kinematics, Dynamics, and Control (2nd Edition)*. Springer, 2010.

[23] Morgan Jones and Matthew M. Peet. Polynomial approximation of value functions and nonlinear controller design with performance bounds.

[24] Milan Korda. Computing controlled invariant sets from data using convex optimization. *SIAM Journal on Control and Optimization*, 58(5):2871–2899, 2020.

[25] Milan Korda, Didier Henrion, and Colin N. Jones. Convergence rates of moment-sum-of-squares hierarchies for optimal control problems. *Systems & Control Letters*, 100:1–5, 2017.

[26] Jean B. Lasserre, Didier Henrion, Christophe Prieur, and Emmanuel Trélat. Nonlinear optimal control via occupation measures
<table><caption>Table 4: Generalized Brockett integrator: lower and upper bounds, and computational times for various degrees of the SIP hierarchy ($D_d$)</caption><thead><tr><th>d &isin; N</th><th>Estimated value of (D<sub>d</sub>)</th><th>Certified LB for (D<sub>d</sub>)</th><th>Value feedback control (CL<sub>V</sub>)</th><th>Solution time (in s)</th><th>Iterations number</th><th>Certification time (in s)</th></tr></thead><tbody><tr><td>2</td><td>1.071</td><td>0.763</td><td>1.071</td><td>220</td><td>28</td><td>73</td></tr><tr><td>3</td><td>1.072</td><td>0.000</td><td>1.071</td><td>5 630</td><td>133</td><td>TL</td></tr><tr><td>4</td><td>1.072</td><td>0.000</td><td>1.070</td><td>147 400</td><td>319</td><td>TL</td></tr></tbody></table>

trajectory. Note that the estimated values of ($D_d$) computed by Algorithm 3.1 with the (inexact) sampling oracle are slightly larger than the value of the best trajectory we computed: thus, these estimates are not valid lower bounds, but only estimates of the value of the minimum time control problem. Regarding the certification of lower bounds, the global optimization solver `SCIP 8` fails to produce tight upper and lower bounds on $\phi(\theta)$, the infeasibility of the solution $\theta$ returned by Algorithm 3.1. Therefore, the resulting certified lower bounds are not tight either. In summary, we obtain a certified optimization gap of 29% for this minimal time control problem.

**6. Discussion.** We apply the dual approach in minimal time control, that consists in searching for maximal subsolutions of the HJB equation, to nonlinear, even non-polynomial, Lipschitz controlled systems. The basis functions used to generate these subsolutions are polynomials, that are subject to semi-infinite constraints. We prove the theoretical convergence of the resulting hierarchy of semi-infinite linear programs, and our numerical tests on three different systems show good convergence properties in practice. These results show that the use of a random sampling oracle allows a good approximation of the value of the control problem. For small systems, it is even possible to obtain tight and certified lower bounds, based on a global optimization solver. Finally, the numerical experiments also show that the computed subsolutions of the HJB equations help to recover near-optimal controls in a closed-loop form. As illustrated in these numerical experiments, the advantage of our approach based on semi-infinite programming, compared to the sum-of-squares approach, is the ability to handle non-polynomial systems. In the numerical example where a polynomial reformulation of the system was possible, the sum-of-squares approach was, however, faster.

A promising avenue for continuing this work is to investigate the use of a other basis of functions to search for an approximate value function, resulting in other semi-infinite programming hierarchies with convergence guarantees. In particular, it would be relevant to use non-differentiable functions in the basis to improve approximation capabilities for non-differentiable value functions. Another avenue of research is to extend the approach and theoretical results to a general optimal control problem, beyond the setting of minimal time control.

**Acknowledgments.** The authors would like to thank Maxime Dupuy, Leo Liberti and Claudia D'Ambrosio for fruitful discussions and advice. The authors would also like to thank two anonymous reviewers for their thoughtful comments.

REFERENCES

[1] Jean-Pierre Aubin. *Viability Theory*. Systems & control: foundations & applications. Springer Science+Business Media, 1991.

[2] Jean-Pierre Aubin and Hélène Frankowska. *Set-Valued Analysis*. Birkhäuser, 2009.

[3] T. Bagby, L. Bos, and N. Levenberg. Multivariate simultaneous approximation. *Constructive Approximation*, 18(4):569-577, 2002.

[4] Richard Bellman. Dynamic programming. *Science*, 153(3731):34-37, 1966. Publisher: American Association for the Advance-ment of Science.

[5] Eloïse Berthier, Justin Carpentier, Alessandro Rudi, and Francis Bach. Infinite-dimensional sums-of-squares for optimal control. In *2022 IEEE 61st Conference on Decision and Control (CDC)*, pages 577-582. IEEE, 2022.

[6] Ksenia Bestuzheva, Mathieu Besançon, Wei-Kun Chen, Antonia Chmiela, Tim Donkiewicz, Jasper van Doornmalen, Leon Eifler, Oliver Gaul, Gerald Gamrath, Ambros Gleixner, Leona Gottwald, Christoph Graczyk, Katrin Halbig, Alexander Hoen, Christopher Hojny, Rolf van der Hulst, Thorsten Koch, Marco Lübbecke, Stephen J. Maher, Frederic Matter, Erik Mühmer, Benjamin Müller, Marc E. Pfetsch, Daniel Rehfeldt, Steffan Schlein, Franziska Schlösser, Felipe Serrano, Yuji Shinano, Boro Sofranac, Mark Turner, Stefan Vigerske, Fabian Wegscheider, Philipp Wellner, Dieter Weninger, and Jakob Witzig. The `SCIP` optimization suite 8.0. Technical report, Optimization Online, 2021.

[7] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B. Shah. Julia: A fresh approach to numerical computing. *SIAM Review*, 59(1):65-98, 2017.
and LMI-relaxations. *SIAM Journal on Control and Optimization*, 47(4):1643-1666, 2008.

[27] Mirko Leomanni, Gabriele Costante, and Francesco Ferrante. Time-optimal control of a multidimensional integrator chain with applications. *IEEE Control Systems Letters*, 6:2371-2376, 2022. Conference Name: IEEE Control Systems Letters.

[28] Daniel Liberzon. *Calculus of Variations and Optimal Control Theory - A Concise Introduction*. Princeton University Press, 2012.

[29] Jérôme Lohéac and Jean-François Scheid. Time-optimal control for a perturbed brockett integrator. In *Variational methods*, pages 454-472. De Gruyter, 2017.

[30] Gianandrea Mannarini, Deepak N. Subramani, Pierre F. J. Lermusiaux, and Nadia Pinardi. Graph-search and differential equations for time-optimal vessel route planning in dynamic ocean waves. *IEEE Transactions on Intelligent Transportation Systems*, 21(8):3581-3593, 2020.

[31] G. J. Olsder. Time-optimal control of multivariable systems near the origin. *Journal of optimization theory and applications*, 16:497-517, 1975. Publisher: Springer.

[32] Antoine Oustry. *Global optimization of nonlinear semi-infinite programming problems: Applications in power systems and control*. phdthesis, Institut Polytechnique de Paris, 2023.

[33] Antoine Oustry, Matteo Tacchi, and Didier Henrion. Inner approximations of the maximal positively invariant set for polynomial dynamical systems. *IEEE Control Systems Letters*, 3(3):733-738, 2019. Conference Name: IEEE Control Systems Letters.

[34] Hans Josef Pesch. A practical guide to the solution of real-life optimal control problems. *Control and Cybernetics*, 23, 1996.

[35] L. S. Pontryagin. *Mathematical Theory of Optimal Processes*. CRC Press, 1987. Google-Books-ID: kwzq0F4cBVAC.

[36] Sebastian Sager, Mathieu Claeys, and Frédéric Messine. Efficient upper and lower bounds for global mixed-integer optimal control. *Journal of Global Optimization*, 61(4):721-743, 2015.

[37] Halil Mete Soner. Optimal control with state-space constraint i. *SIAM Journal on Control and Optimization*, 24(3):552-561, 1986. Publisher: Society for Industrial and Applied Mathematics.

[38] E. Trélat. Optimal control and applications to aerospace: Some results and challenges. *Journal of Optimization Theory and Applications*, 154(3):713-758, 2012.

[39] Robin Verschueren, Hans Joachim Ferreau, Alessandro Zanarini, Mehmet Mercangöz, and Moritz Diehl. A stabilizing nonlinear model predictive control scheme for time-optimal point-to-point motions. In *2017 IEEE 56th Annual Conference on Decision and Control (CDC)*, pages 2525-2530, 2017.

[40] Richard Vinter. Convex duality and nonlinear optimal control. *SIAM Journal on Control and Optimization*, 31(2):518-538, 1993. Publisher: Society for Industrial and Applied Mathematics.

[41] O. von Stryk and R. Bulirsch. Direct and indirect methods for trajectory optimization. *Annals of Operations Research*, 37(1):357-373, 1992.

[42] Saïd Zabi, Isabelle Queinnec, Germain Garcia, and Michel Mazerolles. Time-optimal control for the induction phase of anesthesia. In *20th IFAC World Congress*, page 12708, 2017.

## Appendix A. Technical lemmata.

**LEMMA A.1.** We consider a compact set $Z \subset \mathbb{R}^p$, and the family of compact sets $Z_\delta = \{z \in \mathbb{R}^p : d(z, Z) \le \delta\}$ for any $\delta \ge 0$ and a continuous function $\psi \in C(\mathbb{R}^p)$. Then, the function $\Psi(\delta) = \min_{z \in Z_\delta} \psi(z)$ is continuous at 0.

*Proof.* First of all, we notice that the function $\delta \mapsto \min_{z \in Z_\delta} \psi(z)$ is well-defined, since $\psi$ is continuous and $Z_\delta$ is compact. As $Z_{\delta_1} \subset Z_{\delta_2}$ for any $\delta_1 \le \delta_2$, the function $\Psi$ is non-increasing, which proves that the following limit exists:

$$ (A.1) \qquad \lim_{\delta \to 0^+} \Psi(\delta) = \Psi(0^+) \le \Psi(0). $$

We take a positive sequence $(\delta_k) \in \mathbb{R}_{++}^N$ such that $\delta_k \to 0$. Hence, $\Psi(\delta_k) \to \Psi(0^+)$ by definition of the right-limit. For any $k \in \mathbb{N}$, we define $z_k \in Z_{\delta_k}$ such that $\psi(z_k) = \Psi(\delta_k)$. The sequence $(\delta_k)$ being bounded, we can introduce an upper bound $\bar{\delta}$. Hence, any element of the sequence $(z_n)$ belongs to the compact set $Z_{\bar{\delta}}$, and up to the extraction of a subsequence, converges to a point $z$ being such that $\psi(z) = \Psi(0^+)$ by continuity of $\psi$ and uniqueness of the limit. As $d(z_k, Z)$, the distance between $z_k$ and the compact set $Z$, is bounded above by $\delta_k$ and is non-negative, it converges to 0. By continuity of the distance, we know that $d(z, Z) = 0$ and, thus, $\Psi(0^+) = \psi(z) \ge \Psi(0)$. Together with Eq. (A.1), this yields $\Psi(0^+) = \Psi(0)$. $\square$

**LEMMA A.2.** Let $Q \in C^1(\mathbb{R}^p)$ be a continuously differentiable function, with a locally Lipschitz gradient. Let $Z \subset \mathbb{R}^p$ be a compact set. Then, there exists a constant $A > 0$, and a sequence of polynomials $(w_d(x))_{d \in \mathbb{N}^*}$ such that for all $d \in \mathbb{N}^*$, $w_d \in \mathbb{R}_d[x_1, \dots, x_p]$ and

$$ (A.2) \qquad \sup_{x \in Z} |w_d(x) - Q(x)| \le \frac{A}{d} $$

$$ (A.3) \qquad \sup_{x \in Z} \|\nabla w_d(x) - \nabla Q(x)\| \le \frac{A}{d}. $$
[8] Jérôme Bolte and Edouard Pauwels. Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning. *Mathematical Programming*, 188(1):19–51, 2021.

[9] Loïc Bourdin and Emmanuel Trélat. Pontryagin maximum principle for optimal sampled-data control problems. *IFAC-PapersOnLine*, 48(25):80–84, 2015.

[10] Eduardo F. Camacho and Carlos Bordons Alba. *Model predictive control*. Springer science & business media, 2013.

[11] I. Capuzzo-Dolcetta and P.-L. Lions. Hamilton-jacobi equations with state constraints. *Transactions of the American Mathematical Society*, 318(2):643–683, 1990.

[12] Martina Cerulli, Antoine Oustry, Claudia D'Ambrosio, and Leo Liberti. Convergent algorithms for a class of convex semi-infinite programs. *SIAM Journal on Optimization*, 32(4):2493–2526, 2022. Publisher: Society for Industrial and Applied Mathematics.

[13] Frank H. Clarke. Generalized gradients and applications. *Transactions of the American Mathematical Society*, 205:247–262, 1975.

[14] Frank H. Clarke and Richard B. Vinter. The relationship between the maximum principle and dynamic programming. *SIAM Journal on Control and Optimization*, 25(5):1291–1311, 1987. Publisher: Society for Industrial and Applied Mathematics.

[15] Michael G. Crandall and Pierre-Louis Lions. Viscosity solutions of hamilton-jacobi equations. *Transactions of the American Mathematical Society*, 277(1):1–42, 1983. Publisher: American Mathematical Society.

[16] Hatim Djelassi, Alexander Mitsos, and Oliver Stein. Recent advances in nonconvex semi-infinite programming: Applications and algorithms. *EURO Journal on Computational Optimization*, 9:100006, 2021.

[17] Halina Frankowska. Optimal trajectories associated with a solution of the contingent hamilton-jacobi equation. *Applied Mathematics and Optimization*, 19(1):291–311, 1989.

[18] Hélène Frankowska and Franco Rampazzo. Filippov's and filippov-ważewski's theorems on closed domains. *Journal of Differential Equations*, 161(2):449–478, 2000.

[19] Kurt Otto Friedrichs. The identity of weak and strong extensions of differential operators. *Transactions of the American Mathematical Society*, 55(1):132–151, 1944.

[20] Ronald F. Gariepy and Lawrence Craig Evans. *Measure Theory and Fine Properties of Functions, Revised Edition*. Chapman and Hall/CRC, 2015.

[21] Enrico Giusti. *Minimal Surfaces and Functions of Bounded Variation*. Number 80 in Monographs in Mathematics. Birkhauser Verlag, 1984.

[22] M. A. Goberna and M. A. López. Recent contributions to linear semi-infinite optimization: an update. *Annals of Operations Research*, 271(1):237–278, 2018.

[23] LLC Gurobi Optimization. Gurobi optimizer reference manual, 2021.

[24] Didier Henrion and Milan Korda. Convex computation of the region of attraction of polynomial control systems. *IEEE Transactions on Automatic Control*, 59(2):297–312, 2014. IEEE Transactions on Automatic Control.

[25] Didier Henrion, Jean B. Lasserre, and Carlo Savorgnan. Nonlinear optimal control synthesis via occupation measures. In *2008 47th IEEE Conference on Decision and Control*, pages 4749–4754, 2008. ISSN: 0191-2216.

[26] Didier Henrion and Edouard Pauwels. Linear conic optimization for nonlinear optimal control, July 2014. arXiv:1407.1650 [math].

[27] Daniel Hernández-Hernández, Onésimo Hernández-Lerma, and Michael Taksar. The linear programming approach to deterministic optimal control problems. *Applications Mathematicae*, 24(1):17–33, 1996.

[28] R. Hettich and K. O. Kortanek. Semi-infinite programming: Theory, methods, and applications. *SIAM Review*, 35(3):380–429, 1993.

[29] Lars Hörmander. *The Analysis of Linear Partial Differential Operators I*. Classics in Mathematics. Springer, 2003.

[30] Reza Jazar. *Theory of Applied Robotics: Kinematics, Dynamics, and Control (2nd Edition)*. Springer, 2010.

[31] Morgan Jones and Matthew M. Peet. Polynomial approximation of value functions and nonlinear controller design with performance bounds.

[32] Milan Korda. Computing controlled invariant sets from data using convex optimization. *SIAM Journal on Control and Optimization*, 58(5):2871–2899, 2020.

[33] Milan Korda, Didier Henrion, and Colin N. Jones. Convergence rates of moment-sum-of-squares hierarchies for optimal control problems. *Systems & Control Letters*, 100:1–5, 2017.

[34] Jean B. Lasserre, Didier Henrion, Christophe Prieur, and Emmanuel Trélat. Nonlinear optimal control via occupation measures and LMI-relaxations. *SIAM Journal on Control and Optimization*, 47(4):1643–1666, 2008.

[35] Mirko Leomanni, Gabriele Costante, and Francesco Ferrante. Time-optimal control of a multidimensional integrator chain with applications. *IEEE Control Systems Letters*, 6:2371–2376, 2022. Conference Name: IEEE Control Systems Letters.

[36] Daniel Liberzon. *Calculus of Variations and Optimal Control Theor – A Concise Introduction*. Princeton University Press, 2012.

[37] Jérôme Lohéac and Jean-François Scheid. Time-optimal control for a perturbed brockett integrator. In *Variational methods*, pages 454–472. De Gruyter, 2017.

[38] Marco López and Georg Still. Semi-infinite programming. *European Journal of Operational Research*, 180(2):491–518, 2007.

[39] Gianandrea Mannarini, Deepak N. Subramani, Pierre F. J. Lermusiaux, and Nadia Pinardi. Graph-search and differential equations for time-optimal vessel route planning in dynamic ocean waves. *IEEE Transactions on Intelligent Transportation Systems*, 21(8):3581–3593, 2020.

[40] G. J. Olsder. Time-optimal control of multivariable systems near the origin. *Journal of optimization theory and applications*, 16:497–517, 1975. Publisher: Springer.

[41] Antoine Oustry. *Global optimization of nonlinear semi-infinite programming problems: Applications in power systems and control*. phdthesis, Institut Polytechnique de Paris, 2023.

[42] Antoine Oustry and Martina Cerulli. Convex semi-infinite programming algorithms with inexact separation oracles. *Optimization Letters*, 2024.

[43] Antoine Oustry, Matteo Tacchi, and Didier Henrion. Inner approximations of the maximal positively invariant set for polynomial dynamical systems. *IEEE Control Systems Letters*, 3(3):733–738, 2019. Conference Name: IEEE Control Systems Letters.
We underline that the constant A implicitly depends on p, Q and Z, but not on the polynomial $w_d(x)$, nor on its degree d.

**Proof.** We introduce a constant $R > 0$ such that $Z \subset B(0, R)$, and the function $\tilde{\omega} = \omega * 1_{B(0, R+1)}$, where $\omega$ is the mollifier introduced in the proof of Th. 2.5. We notice that $\tilde{\omega} \in C^\infty(\mathbb{R}^p)$ is supported on $\tilde{Z} = B(0, R+2)$ and constant equal to 1 over $B(0, R)$. We define $\tilde{Q}(x) = Q(x)\tilde{\omega}(x)$, and we notice that (i) $\tilde{Q}$ is supported on the compact set $\tilde{Z}$, which contains $Z$, (ii) for all $x \in Z$, $\tilde{Q}(x) = Q(x)$ and $\nabla \tilde{Q}(x) = \nabla Q(x)$. Applying [3, Th. 1] to the function $\tilde{Q}$, that has a compact support, we know that there exists a constant C such that for any $d \ge 1$, there exists a polynomial $w_d(x)$ of degree at most d such that

$$ (A.4) \qquad \sup_{x \in Z} |w_d(x) - \tilde{Q}(x)| \le \frac{C}{d} \kappa(\frac{1}{d}) \le C \kappa(\frac{1}{d}), $$

$$ (A.5) \qquad \sup_{x \in Z} |\partial_i (w_d - \tilde{Q})(x)| \le C \kappa(\frac{1}{d}) $$

where $\kappa(\delta) = \sup_{1 \le i \le p} \sup_{\substack{x,y \in \mathbb{R}^p \\ |x-y| \le \delta}} |\partial_i \tilde{Q}(x) - \partial_i \tilde{Q}(y)|$. We define $\tilde{Z} = \{x \in \mathbb{R}^p : d(x, Z) \le 1\}$. Since $\partial \tilde{Q}$ is uniformly null outside $\tilde{Z}$, and assuming that $\delta \le 1$, we notice that

$$ \kappa(\delta) = \sup_{1 \le i \le p} \sup_{\substack{x,y \in \tilde{Z} \times \mathbb{R}^p \\ |x-y| \le \delta}} |\partial_i \tilde{Q}(x) - \partial_i \tilde{Q}(y)| = \sup_{1 \le i \le p} \sup_{\substack{x,y \in \tilde{Z} \times \tilde{Z} \\ |x-y| \le \delta}} |\partial_i \tilde{Q}(x) - \partial_i \tilde{Q}(y)|. $$

We note that $\partial_i \tilde{Q}(x) = \tilde{\omega}(x)\partial_i Q(x) + Q(x)\partial_i \tilde{\omega}(x)$, and therefore, $|\partial_i \tilde{Q}(x) - \partial_i \tilde{Q}(y)| = |\tilde{\omega}(x)(\partial_i Q(x) - \partial_i Q(y)) + \partial_i Q(y)(\tilde{\omega}(x) - \tilde{\omega}(y)) + Q(x)(\partial_i \tilde{\omega}(x) - \partial_i \tilde{\omega}(y)) + \partial_i \tilde{\omega}(y)(Q(x) - Q(y))|$. We use, then, the triangle inequality and the facts that (i) $\tilde{\omega}$ is $C^\infty$, therefore bounded, Lipschitz continuous, and with a Lipschitz-continuous gradient over $\tilde{Z}$ and (ii) $Q$ is continuously differentiable, therefore bounded, and Lipschitz continuous over $\tilde{Z}$; by assumption it has a Lipschitz continuous gradient over the compact set $\tilde{Z}$. We deduce that $\partial_i \tilde{Q}$ is Lipschitz continuous over $\tilde{Z}$: there exists $L_i > 0$ such that $\sup_{\substack{x,y \in \tilde{Z} \times \tilde{Z} \\ |x-y| \le \delta}} |\partial_i \tilde{Q}(x) - \partial_i \tilde{Q}(y)| \le L_i \delta$, for all $\delta \in [0, 1]$. Defining $L = \max_i L_i$, we deduce

$\kappa(\delta) \le L\delta$. Then Eq. (A.4) reads $\sup_{x \in Z} |w_d(x) - \tilde{Q}(x)| \le \frac{CL}{d}$, and Eq. (A.5) reads $\sup_{x \in Z} |\partial_i (w_d - \tilde{Q})(x)| \le \frac{CL}{d}$ for all $i \in \{1, \dots, p\}$. We also deduce that $\sup_{x \in Z} \|\nabla w_d(x) - \nabla \tilde{Q}(x)\| \le \frac{CL_p}{d}$. Defining $A = CL_p$, and noticing that for all $x \in Z$, $\tilde{Q}(x) = Q(x)$ and $\nabla \tilde{Q}(x) = \nabla Q(x)$, one obtains the claimed statement. $\square$

**LEMMA A.3.** Under Assumption 1 and Assumption 2, we consider an admissible trajectory $(x(\cdot), u(\cdot))$ over $[0, t_1]$ of the minimal time control problem (2.1)-(2.4) starting from $(0, x_0)$. Then, for almost all $t \in [0, t_1]$, $f(t, x(t), u(t)) \in T_X(x(t))$.

**Proof.** We reduce to a time-invariant controlled system: we define, for any $y = (t, x) \in \mathbb{R}^{n+1}$ and $u \in \mathbb{R}^m$, $\tilde{f}(y, u) = \begin{pmatrix} 1 \\ f(y, u) \end{pmatrix}$, and the constant set-valued map $\tilde{U}(y) = U$. The control system $(\tilde{f}, \tilde{U})$ is a Marchaud control system [1, Def. 6.1.3], since: (i) the graph of $\tilde{U}$ is closed (ii) $\tilde{f}$ is continuous (iii) the velocity subsets $\{\tilde{f}(y, u) : u \in \tilde{U}(y)\}$ are convex due to Assumption 1, and (iv) the function $f$ has a linear growth since it is Lipschitz continuous, and the set-valued map are bounded, thus also has a linear growth. We define the set $C = \mathbb{R}_+ \times X$ and notice that the control $u(\cdot)$ regulates a trajectory $y(t) = \begin{pmatrix} t \\ x(t) \end{pmatrix}$ that remains in $C$, therefore according to [1, Th. 6.1.4], for all most all $t \in [0, t_1]$, $u(t) \in \{u \in \tilde{U}(y(t)): \tilde{f}(y(t), u(t)) \in T_C(y(t))\}$. We notice that $T_C(y(t)) \subset \mathbb{R} \times T_X(x(t))$, implying that $f(y(t), u(t)) \in T_X(x(t))$. $\square$

**LEMMA A.4.** For any locally Lipschitz continuous function $F: \mathbb{R}^p \to \mathbb{R}$, and for any Lipschitz continuous curve $y: [0, T] \to \mathbb{R}^p$, the function $t \mapsto F(y(t))$ is differentiable a.e. and satisfies

$$ (A.6) \qquad \min_{g \in \partial^c F(y(t))} g^\top \dot{y}(t) \le \frac{d}{dt} (F(y(t))) \le \max_{g \in \partial^c F(y(t))} g^\top \dot{y}(t), $$

for almost all $t \in [0, T]$.

The particular functions $F$ for which these three quantities are equal are called path-differentiable in [7].
[44] Hans Josef Pesch. A practical guide to the solution of real-life optimal control problems. *Control and Cybernetics*, 23, 1996.

[45] L. S. Pontryagin. *Mathematical Theory of Optimal Processes*. CRC Press, 1987. Google-Books-ID: kwzq0F4cBVAC.

[46] Sebastian Sager, Mathieu Claeys, and Frédéric Messine. Efficient upper and lower bounds for global mixed-integer optimal control. *Journal of Global Optimization*, 61(4):721–743, 2015.

[47] Halil Mete Soner. Optimal control with state-space constraint i. *SIAM Journal on Control and Optimization*, 24(3):552–561, 1986. Publisher: Society for Industrial and Applied Mathematics.

[48] E. Trélat. Optimal control and applications to aerospace: Some results and challenges. *Journal of Optimization Theory and Applications*, 154(3):713–758, 2012.

[49] Robin Verschueren, Hans Joachim Ferreau, Alessandro Zanarini, Mehmet Mercangöz, and Moritz Diehl. A stabilizing nonlinear model predictive control scheme for time-optimal point-to-point motions. In *2017 IEEE 56th Annual Conference on Decision and Control (CDC)*, pages 2525–2530, 2017.

[50] Richard Vinter. Convex duality and nonlinear optimal control. *SIAM Journal on Control and Optimization*, 31(2):518–538, 1993. Publisher: Society for Industrial and Applied Mathematics.

[51] O. von Stryk and R. Bulirsch. Direct and indirect methods for trajectory optimization. *Annals of Operations Research*, 37(1):357–373, 1992.

[52] Saïd Zabi, Isabelle Queinnec, Germain Garcia, and Michel Mazerolles. Time-optimal control for the induction phase of anesthesia. In *20th IFAC World Congress*, page 12708, 2017.

## Appendix A. Technical lemmata.

**LEMMA A.1.** We consider a compact set $Z \subset \mathbb{R}^p$, and the family of compact sets $Z_\delta = \{z \in \mathbb{R}^p : d(z, Z) \le \delta\}$ for any $\delta \ge 0$ and a continuous function $\psi \in C(\mathbb{R}^p)$. Then, the function $\Psi(\delta) = \min_{z \in Z_\delta} \psi(z)$ is continuous at 0.

*Proof.* First of all, we notice that the function $\delta \mapsto \min_{z \in Z_\delta} \psi(z)$ is well-defined, since $\psi$ is continuous and $Z_\delta$ is compact. As $Z_{\delta_1} \subset Z_{\delta_2}$ for any $\delta_1 \le \delta_2$, the function $\Psi$ is non-increasing, which proves that the following limit exists:

$$ (A.1) \qquad \lim_{\delta \to 0^+} \Psi(\delta) = \Psi(0^+) \le \Psi(0). $$

We take a positive sequence $(\delta_k) \in \mathbb{R}_{++}^N$ such that $\delta_k \to 0$. Hence, $\Psi(\delta_k) \to \Psi(0^+)$ by definition of the right-limit. For any $k \in \mathbb{N}$, we define $z_k \in Z_{\delta_k}$ such that $\psi(z_k) = \Psi(\delta_k)$. The sequence $(\delta_k)$ being bounded, we can introduce an upper bound $\bar{\delta}$. Hence, any element of the sequence $(z_n)$ belongs to the compact set $Z_{\bar{\delta}}$, and up to the extraction of a subsequence, converges to a point $z$ being such that $\psi(z) = \Psi(0^+)$ by continuity of $\psi$ and uniqueness of the limit. As $d(z_k, Z)$, the distance between $z_k$ and the compact set $Z$, is bounded above by $\delta_k$ and is non-negative, it converges to 0. By continuity of the distance, we know that $d(z, Z) = 0$ and, thus, $\Psi(0^+) = \psi(z) \ge \Psi(0)$. Together with Eq. (A.1), this yields $\Psi(0^+) = \Psi(0)$. $\square$

**LEMMA A.2.** Let $Q \in C^1(\mathbb{R}^p)$ be a continuously differentiable function, with a locally Lipschitz gradient. Let $Z \subset \mathbb{R}^p$ be a compact set. Then, there exists a constant $A > 0$, and a sequence of polynomials $(w_d(x))_{d \in \mathbb{N}^*}$ such that for all $d \in \mathbb{N}^*$, $w_d \in \mathbb{R}_d[x_1, \dots, x_p]$ and

$$ (A.2) \qquad \sup_{x \in Z} |w_d(x) - Q(x)| \le \frac{A}{d} $$

$$ (A.3) \qquad \sup_{x \in Z} \|\nabla w_d(x) - \nabla Q(x)\| \le \frac{A}{d}. $$

We underline that the constant $A$ implicitly depends on $p$, $Q$ and $Z$, but not on the polynomial $w_d(x)$, nor on its degree $d$.

*Proof.* We introduce a constant $R > 0$ such that $Z \subset B(0, R)$, and the function $\tilde{\omega} = \omega \star 1_{B(0, R+1)}$, where $\omega$ is the mollifier introduced in the proof of Th. 2.5. We notice that $\tilde{\omega} \in C^\infty(\mathbb{R}^p)$ is supported on $\tilde{Z} = B(0, R+2)$ and constant equal to 1 over $B(0, R)$. We define $\tilde{Q}(x) = Q(x)\tilde{\omega}(x)$, and we notice that (i) $\tilde{Q}$ is supported on the compact set $\tilde{Z}$, which contains $Z$, (ii) for all $x \in Z$, $\tilde{Q}(x) = Q(x)$ and $\nabla \tilde{Q}(x) = \nabla Q(x)$. Applying [3, Th. 1] to the function $\tilde{Q}$, that has a compact support, we know that there exists a constant $C$ such that for any $d \ge 1$, there exists a polynomial $w_d(x)$ of degree at most $d$ such that

$$ (A.4) \qquad \sup_{x \in Z} |w_d(x) - \tilde{Q}(x)| \le \frac{C}{d} \kappa(\frac{1}{d}) \le C\kappa(\frac{1}{d}), $$

$$ (A.5) \qquad \sup_{x \in Z} |\partial_t(w_d - \tilde{Q})(x)| \le C\kappa(\frac{1}{d}) $$
Proof. First, we notice that the functions $t \mapsto y(t)$ and $t \mapsto F(y(t))$ are Lipschitz continuous, therefore differentiable a.e. on $[0, T]$ due to the Rademacher theorem [16]. Hence, for almost all $t \in [0, T]$, both $y(t)$ and $F(y(t))$ are differentiable at $t$. We consider such a $t$, and we show that (A.6) holds for this particular $t$, which we prove the Lemma. Since $y$ is differentiable at $t$, $r(h) = y(t+h) - y(t) - h\dot{y}(t)$ is in $o_{h \to 0}(h)$.
Since $s \mapsto F(y(s))$ is differentiable at $t$, the following holds

$$ (A.7) \qquad \frac{d}{dt}(F(y(t))) = \lim_{h \to 0, h > 0} \frac{F(y(t+h)) - F(y(t))}{h} $$

$$ (A.8) \qquad = \lim_{h \to 0, h > 0} \frac{F(y(t) + h\dot{y}(t) + r(h)) - F(y(t))}{h} $$

Since $r(h) = o_{h \to 0}(h)$ and $F$ is locally Lipschitz, we know that $\lim_{h \to 0, h > 0} \frac{F(y(t))-F(y(t)+r(h))}{h} = 0$. Summing this with Eq. (A.8), we deduce that

$$ (A.9) \qquad \frac{d}{dt}(F(y(t))) = \lim_{h \to 0, h > 0} \frac{F(y(t) + r(h) + h\dot{y}(t)) - F(y(t) + r(h))}{h} $$

$$ (A.10) \qquad \leq \limsup_{\substack{y' \to y(t) \\ h \to 0, h > 0}} \frac{F(y' + h\dot{y}(t)) - F(y')}{h} = F^{\circ}(y(t), \dot{y}(t)), $$

where $F^{\circ}(y; v) = \limsup_{\substack{y' \to y \\ h \to 0, h > 0}} \frac{F(y'+hv)-F(y')}{h}$ is the $F^{\circ}(y; h)$ Clarke's directional derivative at $y \in \mathbb{R}^p$ in the direction

$v \in \mathbb{R}^p$. The inequality follows from the fact that $y(t) + r(h) \to y(t)$. By property of the Clarke subdifferential [12], we also know that $F^{\circ}(y; v) = \max_{g \in \partial^c F(y)} g^{\top}v$. Hence, in particular,

$$ (A.11) \qquad \frac{d}{dt}(F(y(t))) \leq \max_{g \in \partial^c F(y(t))} g^{\top} \dot{y}(t). $$

The reasoning that proved Eq. (A.11) is also applicable to $-F$, that is also locally Lipschitz, and such that $s \mapsto (-F)(s)$ is differentiable at $t$. Therefore,

$$ (A.12) \qquad \frac{d}{dt}(-F(y(t))) \leq \max_{g \in \partial^c(-F)(y(t))} g^{\top} \dot{y}(t). $$

As $\partial^c(-F)(y(t)) = -\partial^c F(y(t))$ by property of the Clarke subdifferential, we deduce that

$$ (A.13) \qquad -\frac{d}{dt}(F(y(t))) \leq \max_{g \in \partial^c F(y(t))} -g^{\top}\dot{y}(t) = -\min_{g \in \partial^c F(y(t))} g^{\top}\dot{y}(t), $$

and therefore $\frac{d}{dt}(F(y(t))) \geq \min_{g \in \partial^c F(y(t))} g^{\top}\dot{y}(t)$. $\square$

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