# Equivalent Formulations of Optimal Control Problems with Maximum Cost and Applications

Emilio Molina, Alain Rapaport, Héctor Ramírez

► To cite this version:

Emilio Molina, Alain Rapaport, Héctor Ramírez. Equivalent Formulations of Optimal Control Problems with Maximum Cost and Applications. Journal of Optimization Theory and Applications, 2022, 195, pp.953-975. 10.1007/s10957-022-02094-z. hal-03746480

HAL Id: hal-03746480

https://hal.inrae.fr/hal-03746480

Submitted on 5 Aug 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Equivalent Formulations of Optimal Control Problems with Maximum Cost and Applications

Emilio Molina · Alain Rapaport · Héctor Ramírez

**Abstract** We revisit the optimal control problem with maximum cost with the objective to provide different equivalent reformulations suitable to numerical methods. We propose two reformulations in terms of extended Mayer problems with state constraints, and another one in terms of a differential inclusion with upper-semicontinuous right member without state constraint. For the latter we also propose a scheme that approximates from below the optimal value. These approaches are illustrated and discussed in several examples.

**Keywords** Optimal control · maximum cost · Mayer problem · state constraint · differential inclusion · numerical schemes · SIR model

**Mathematics Subject Classification (2000)** 49N90 · 49J35 · 49J45 · 65K10 · 90-08

## 1 Introduction

We consider the optimal control problem which consists in minimizing the maximum of a scalar function over a time interval

$$ \inf_{u(\cdot)} \sup_{t \in [t_0, T]} y(t) $$

where $y(t) = \theta(t, \xi(t))$ and $\xi(\cdot)$ is the solution of a controlled dynamics $\dot{\xi} = \phi(\xi, u)$, $\xi(t_0) = \xi_0$. This problem is not in the usual Mayer, Lagrange or Bolza forms of the optimal control theory, and therefore is not suitable to use the classical necessary optimality conditions of Pontryagin Maximum Principle or existing solving algorithms (based on direct method, shooting or Hamilton-Jacobi-Bellman equation). However, this problem falls into the class of optimal control with $L^\infty$ criterion, for which several characterizations of the value function have been proposed in the literature [3, 4, 12]. Typically, the value function is solution, in a general sense, of partial differential equation of the form

$$ \min \left( \partial_t V + \inf_u \partial_\xi V \cdot \phi(\xi, u) , V - \theta \right) = 0 $$

Emilio Molina

Department of Mathematical Engineering & Center for Mathematical Modeling

Universidad de Chile, Santiago, Chile

and LJLL, Sorbonne Université, CNRS, Paris, France

and CaGe, INRIA, Paris, France

emolina@dim.uchile.cl

Alain Rapaport, Corresponding author

MISTEA, Université Montpellier, INRAE, Institut Agro

Montpellier, France

alain.rapaport@inrae.fr

Héctor Ramírez

Department of Mathematical Engineering & Center for Mathematical Modeling

Universidad de Chile, Santiago, Chile

hramirez@dim.uchile.cl
without boundary condition. Nevertheless, although necessary optimality conditions and numerical procedures have been formulated [2, 8, 9, 11], there is no practical numerical tool to solve such problems as it exists for Mayer problems, to the best of our knowledge. The aim of the present work is to study different reformulations of this problem into Mayer form in higher dimension with possibly state or mixed constraints, for which existing numerical methods can be used. Indeed, it has already been underlined in the literature that discrete-time optimal control problems with maximum cost do not satisfy the Principle of Optimality but can be transformed into problems of higher dimension with additively separable objective functions [14, 15]. We pursue here this idea but in the continuous time framework, which faces the lack of differentiability of the max function.

This manuscript is organized as follows. In Section 2, we establish the setup and the hypotheses of this article, and define the problem. In Section 3, we provide equivalent formulations of the studied problem in the form of two Mayer problems with fixed initial condition, and under state or mixed constraint. In Section 4, we propose another formulation in terms of differential inclusion but without constraints, and then we show how the optimal value can be approximated from below by a sequence of more regular Mayer problems. Section 5 is devoted to numerical illustrations. We consider two problems for which the optimal solution can be determined explicitly (one borrowed from epidemiology), which allows to estimate and compare the numerical performances of the different formulations. These problems have been chosen linear with respect to the control variable in order to present discontinuous optimal controls, which are known to be numerically more sensitive. More precisely, the optimal solution of the first problem has pure bang-bang controls, while the second one possesses a singular arc. We discuss the issues arising in the numerical implementations of the different formulations, and compare numerically with $L^p$ approximations. Finally, we discuss in Section 6 about the potential merits of the different formulations as practical methods to compute optimal solution of $L^\infty$ control problems.

## 2 Problem and Hypotheses

We shall consider autonomous dynamical systems defined on an invariant domain $\mathcal{D}$ of $\mathbb{R}^{n+1}$ of the form

$$ \begin{cases} \dot{x} = f(x, y, u), \\ \dot{y} = g(x, y, u), \end{cases} \qquad (1) $$

(where $g$ is a scalar function) where the values of the control $u(\cdot)$ belong to a given set $U \subset \mathbb{R}^p$. More specifically, throughout the paper, we shall assume that the following properties are fulfilled.

### Assumption 1

i. $U$ is a compact set.

ii. The maps $f$ and $g$ are $C^1$ on $\mathcal{D} \times U$.

iii. The maps $f$ and $g$ have linear growth, that is there exists a number $C > 0$ such that

$$ \|f(x, y, u)\| + |g(x, y, u)| \le C(1 + \|x\| + \|y\|), \quad (x, y) \in \mathcal{D}, u \in U. $$

For instance, $y(\cdot)$ can be a smooth output of a dynamics

$$ \dot{x} = f(x, u), \quad y = h(x), $$

which can be rewritten as

$$ \begin{cases} \dot{x} = f(x, u), \\ \dot{y} = g(x, u) := \nabla h(x)^T \cdot f(x, u). \end{cases} $$

Let $\mathcal{U}$ be the set of measurable functions $u(\cdot) : [0, T] \mapsto U$ and consider $(x_0, y_0) \in \mathcal{D}, T > 0$. Under the usual arguments of the theory of ordinary differential equations, Assumption 1 ensures that for any $u(\cdot) \in \mathcal{U}$ there exists a unique absolutely continuous solution $(x(\cdot), y(\cdot))$ of (1) on $[0, T]$ for the initial condition $(x(0), y(0)) = (x_0, y_0)$ (see, for instance, [10]). Define then the solutions set

$$ \mathcal{S} := \{(x(\cdot), y(\cdot)) \in \mathcal{AC}([0, T], \mathbb{R}^{n+1}), \text{ sol. of (1) for } u(\cdot) \in \mathcal{U} \text{ with } (x(0), y(0)) = (x_0, y_0)\}. $$

We consider then the optimal control problem which consists in minimizing the "peak" of the function $y(\cdot)$

$$ \mathcal{P}: \inf_{u(\cdot) \in \mathcal{U}} \left( \max_{t \in [0, T]} y(t) \right) = \inf_{(x(\cdot), y(\cdot)) \in \mathcal{S}} \left( \max_{t \in [0, T]} y(t) \right). $$
## 3 Formulations With Constraint

A first approach considers the family of constrained sets of solutions

$$S_z := \{(x, y) \in S, y(t) \le z, t \in [0, T]\}, \quad (z \in \mathbb{R})$$

and to look for the optimization problem

$$\inf\{z; S_z \neq \emptyset\}.$$

This problem can be reformulated as a Mayer problem

$$P_0: \inf_{u(\cdot) \in U} z(T)$$

for the extended dynamics in $\mathcal{D} \times \mathbb{R}$

$$\left\{ \begin{array}{l} \dot{x} = f(x, y, u), \\ \dot{y} = g(x, y, u), \\ \dot{z} = 0, \end{array} \right.$$

under the state constraint

$$C: \quad z(t) - y(t) \ge 0, \ t \in [0, T]$$

where $z(0)$ is free. Direct methods can be used for such a problem. However, as $z(0)$ is free, solutions are not sought among solutions of a Cauchy problem, which prevents using other methods based on dynamic programming such as the Hamilton-Jacobi-Bellman equation.

We propose another extended dynamics in $\mathcal{D} \times \mathbb{R}$ with an additional control $v(\cdot)$ with values in $[0, 1]$

$$\left\{ \begin{array}{l} \dot{x} = f(x, y, u), \\ \dot{y} = g(x, y, u), \\ \dot{z} = \max(g(x, y, u), 0)(1 - v). \end{array} \right. \qquad (2)$$

Let $\mathcal{V}$ be the set of measurable functions $v: [0, T] \mapsto [0, 1]$. Note that under Assumption 1, for any $(x_0, y_0, z_0) \in \mathcal{D} \times \mathbb{R}$ and $(u, v) \in \mathcal{U} \times \mathcal{V}$, there exists an unique absolutely solution $(x(\cdot), y(\cdot), z(\cdot))$ of (2) on $[0, T]$ for the initial condition $(x(0), y(0), z(0)) = (x_0, y_0, z_0)$. Here, we fix the initial condition with $z_0 = y_0$ and consider the Mayer problem

$$P_1: \quad \inf_{(u(\cdot), v(\cdot)) \in \mathcal{U} \times \mathcal{V}} z(T) \quad \text{under the constraint } C$$

and shows its equivalence with problem $\mathcal{P}$. We first consider fixed controls $u(\cdot)$.

**Proposition 3.1** For any control $u(\cdot) \in \mathcal{U}$, the optimal control problem

$$\inf_{v \in \mathcal{V}} z(T) \quad \text{under the constraint } C \qquad (3)$$

admits an optimal solution. Moreover, an optimal solution verifies

$$z(T) = \max_{t \in [0, T]} y(t) \qquad (4)$$

and is reached for a control $v(\cdot)$ that takes values in $\{0, 1\}$.

*Proof* From equations (2), one get that any solution $z(\cdot)$ is non decreasing, and as $z$ satisfies the constraint $z \ge y$, we deduce that one has

$$z(T) \ge \max_{t \in [0, T]} y(t) \qquad (5)$$

for any solution of (2), and thus

$$\max_{t \in [0, T]} y(t) \le \inf_{v \in \mathcal{V}} z(T) \quad \text{under the constraint } z(t) \ge y(t), \ t \in [0, T].$$

Let $x(\cdot), y(\cdot)$ be the solution of (1) for the control $u(\cdot)$ and let $I$ be the set of invisible points from the left of $y$, that is

$$I := \{t \in (0, T), y(t') > y(t) \text{ for some } t' < t\}.$$
Consider then the control

$$v(t) = \begin{cases} 1, & t \in \text{int } I, \\ 0, & t \notin \text{int } I. \end{cases} \qquad (6)$$

When $I$ is empty, $y(\cdot)$ is a non decreasing function and, when $v(t) = 0$ for all $t \in [0, T]$, one has $z(t) = y(t)$ for any $t \in [0, T]$. Therefore, one has

$$z(T) = y(T) = \max_{t \in [0,T]} y(t).$$

When $I$ is non empty, there exists, from the sun rising Lemma [17], a countable set of disjoint non-empty intervals $I_n = (a_n, b_n)$ of $[0, T]$ such that

- the interior of $I$ is the union of the intervals $I_n$,

- one has $y(a_n) = y(b_n)$ if $b_n \neq T$,

- if $b_n = T$, then $y(a_n) \ge y(b_n)$.

Note that when $t \notin \text{int } I$, one has $y(t) \ge y(t')$ for any $t' \le t$. Therefore, the solution $z$ with control (6) verifies

$$z(t) = \begin{cases} y(t), & t \notin \text{int } I, \\ y(a_n), & t \in I_n \text{ for some } n, \end{cases}$$

(see Figure 1 as an illustration). Let $\bar{t}$ be the largest time in $[0, T]$ such that

$$y(\bar{t}) = \max_{t \in [0,T]} y(t),$$

which implies that any point $t' > \bar{t}$ in $[0, T]$ is invisible from the left. Therefore, one has $z(T) = z(\bar{t})$. As for any $t$ $z(t)$ is equal to $y(\tau)$ for some $\tau$, one has necessarily $z(\bar{t}) \le y(\bar{t})$. Thus, from (5) we obtain

$$\max_{t \in [0,T]} y(t) = z(T)$$

and deduce

$$\max_{t \in [0,T]} y(t) = \inf_{v(\cdot) \in \mathcal{V}} z(T) \text{ under the constraint } \mathcal{C}.$$

<figure><img src="image_2.png" /><figcaption>Fig. 1 Illustration of the function $z$ (in red) corresponding to a function $y$ (in blue) with the control given by expression (6).</figcaption></figure>
**Remark 3.1** The proof of Proposition 3.2 gives an optimal construction of $z(\cdot)$ which is the lower envelope of non decreasing continuous functions above the function $y(\cdot)$, as depicted on Figure 1. However, there is no uniqueness of the optimal control $v(\cdot)$. Any admissible solution $z(\cdot)$ that is above $y(\cdot)$ and such that $z(t) = \hat{y}$ for $t \ge \hat{t} = \min\{t \in (0, T], y(t) = \hat{y}\}$, where $\hat{y} := \max_{s \in [0, T]} y(s)$, is also optimal.

We then obtain the equivalence between problems $\mathcal{P}_1$ and $\mathcal{P}$ in the following sense.

**Proposition 3.2** If $(u^*(\cdot), v^*(\cdot))$ is optimal for Problem $\mathcal{P}_1$, then $u^*(\cdot)$ is optimal for Problem $\mathcal{P}$. Conversely, if $u^*(\cdot)$ is optimal for Problem $\mathcal{P}$, then $(u^*(\cdot), v^*(\cdot))$ is optimal for Problem $\mathcal{P}_1$ where $v^*(\cdot)$ is optimal for the problem (3) for the fixed control $u^*(\cdot)$.

Let us give another equivalent Mayer problem but with a mixed constraint (this will be useful in the next section). We consider again the extended dynamics (2), with a control $v(\cdot)$ which values belong to $[0, 1]$ and the initial conditions $(x(0), y(0), z(0)) = (x_0, y_0, y_0)$. Define then the mixed constraint

$$ \mathcal{C}_m: \max(y(t) - z(t), 0)(1 - v(t)) + z(t) - y(t) \ge 0, \quad \text{a.e. } t \in [0, T] $$

and the optimal control problem

$$ \mathcal{P}_2: \inf_{(u(\cdot), v(\cdot)) \in \mathcal{U} \times \mathcal{V}} z(T) \quad \text{under the constraint } \mathcal{C}_m. $$

**Proposition 3.3** Problems $\mathcal{P}_1$ and $\mathcal{P}_2$ are equivalent.

*Proof* One can immediately see that for any admissible solution that satisfies constraint $\mathcal{C}$, the constraint $\mathcal{C}_m$ is necessarily fulfilled as $\max(y - z, 0)$ is identically null.

Conversely, fix an admissible control $u(\cdot)$ and consider a control $v(\cdot)$ that satisfies $\mathcal{C}_m$. We show that this implies that the solution $(y(\cdot), z(\cdot))$ verifies necessarily $z(t) \ge y(t)$ for any $t \in [0, T]$. If not, consider the non-empty set

$$ E := \{t \in [0, T]; z(t) - y(t) < 0\}. $$

which is open as $z - y$ is continuous. Note that one has $\dot{z}(t) - \dot{y}(t) \ge 0$ for a.e. $t \in E$. Therefore $z - y$ is non decreasing in $E$ and we deduce that for any $t \in E$, the interval $[0, t]$ is necessarily included in $E$, which then contradicts the initial condition $z(0) = y(0)$.

# 4 Formulation Without State Constraint

We posit $\Pi = (x, y, z) \in \mathcal{D} \times \mathbb{R}$ and consider the differential inclusion

$$ \dot{\Pi} \in F(\Pi) := \bigcup_{(u,v) \in \mathcal{U} \times [0,1]} \begin{bmatrix} f(x, y, u) \\ g(x, y, u) \\ h(x, y, z, u, v) \end{bmatrix} \quad (7) $$

with

$$ h(x, y, z, u, v) = \max(g(x, y, u), 0)(1 - v1_{\mathbb{R}^+}(z - y)) $$

where $1_{\mathbb{R}^+}$ is the indicator function

$$ 1_{\mathbb{R}^+}(\zeta) = \begin{cases} 1, & \zeta \ge 0, \\ 0, & \zeta < 0. \end{cases} $$

Let $\Pi_0 = (x_0, y_0, y_0)$ and denote by $S_\ell$ the set of absolutely continuous solutions of (7) with $\Pi(0) = \Pi_0 \in \mathcal{D} \times \mathbb{R}$. We consider the Mayer problem

$$ \mathcal{P}_3: \inf_{\Pi(\cdot) \in S_\ell} z(T). $$

# Assumption 2

$$ \forall (x, y) \in \mathcal{D}, \quad G(x, y) := \bigcup_{u \in \mathcal{U}} \begin{bmatrix} f(x, y, u) \\ g(x, y, u) \end{bmatrix} \quad \text{is convex.} $$
**Proposition 4.1** Under Assumption 2, problem P3 admits an optimal solution. Moreover, any optimal solution $\Pi(\cdot) = (x(\cdot), y(\cdot), z(\cdot))$ verifies

$$z(T) = \max_{t \in [0,T]} y(t)$$

with $(x(\cdot), y(\cdot))$ solution of (1) for some control $u(\cdot) \in U$, that in turn is optimal for problem P.

*Proof* We fix the initial condition $\Pi(0) = \Pi_0$ and consider the augmented dynamics

$$\dot{\Pi} \in F^{\dagger}(\Pi) := \bigcup_{(u,v,\alpha) \in U \times [0,1]^2} \begin{bmatrix} f(x, y, u) \\ g(x, y, u) \\ h^{\dagger}(x, y, z, u, v, \alpha) \end{bmatrix} \quad (8)$$

with

$$h^{\dagger}(x, y, z, u, v, \alpha) = (1 - \alpha)h(x, y, z, u, v) + \alpha \max_{w \in U} h(x, y, z, w, 0).$$

Under Assumption 2, the values of $F^{\dagger}$ are convex compact. One can straightforwardly check that the set-valued map $F^{\dagger}$ is upper semi-continuous<fnref n="1" /> with linear growth. Therefore, the reachable set $S_{\ell}^{\dagger}(T)$ (where $S_{\ell}^{\dagger}$ denotes the set of absolutely continuous solutions of (8) with $\Pi(0) = \Pi_0$) is compact (see for instance [1, Proposition 3.5.5]). Then, there exists a solution $\Pi^{*}(\cdot) = (x^{*}(\cdot), y^{*}(\cdot), z^{*}(\cdot))$ of (8) which minimizes $z(T)$.

Note that any admissible solution $(x(\cdot), y(\cdot), z(\cdot))$ of system (2) that satisfies the constraint $C_m$ belongs to $S_{\ell} \subset S_{\ell}^{\dagger}$. We then get the inequality

$$z^*(T) \le \inf\{z(T); (x(\cdot), y(\cdot), z(\cdot)) \text{ sol. of (2) with } C_m\}. \quad (9)$$

Let us show that any solution $\Pi(\cdot) = (x(\cdot), y(\cdot), z(\cdot))$ in $S_{\ell}$ verifies

$$z(T) \ge \max_{t \in [0,T]} y(t). \quad (10)$$

We show that one has $z(t) \ge y(t)$ for any $t \in [0, T]$. We proceed by contradiction, as in the proof of Proposition 3.3. If the set $E = \{t \in (0, T); z(t) - y(t) < 0\}$ is non-empty, one has $\dot{z}(t) - \dot{y}(t) \ge 0$ for a.e. $t \in E$ which implies, by continuity, that one has $z(0) - y(0) < 0$ which contradicts the initial condition $z(0) = y(0)$. Moreover, as the map $h$ is non-negative, $z(\cdot)$ is non decreasing and we conclude that (10) is verified.

On another hand, thanks to Assumptions 1 and 2, we can apply Filippov's Lemma to the set-valued map $G$, which asserts that $(x(\cdot), y(\cdot))$ is solution of (1) for a certain $u(\cdot) \in U$. With(10), we obtain

$$z^*(T) \ge \max_{t \in [0,T]} y^*(t) \ge \inf_{u \in U} \left\{ \max_{t \in [0,T]} y(t); (x(\cdot), y(\cdot)) \text{ sol. of (1)} \right\} \quad (11)$$

where $(x^*(\cdot), y^*(\cdot))$ is solution of (1) for a certain $u^*(\cdot) \in U$.

Finally, inequalities (9) and (11) with Propositions 3.2 and 3.3 show that $z^*(T)$ is reached by a solution of (2) under the constraint $C_m$, and that $u^*(\cdot)$ is optimal for problem P. We also conclude that the optimal value $z^*(T)$ is reached by a solution in $S_{\ell}$, which is thus optimal for problem P3.

**Remark 4.1** Let us stress that the function $h$ is not continuous, which does not allow to use Filippov's Lemma for the set valued map $F$. This means that one cannot guarantee a priori that an absolutely continuous solution $\Pi(\cdot) = (x(\cdot), y(\cdot), z(\cdot))$ can be synthesized by a measurable control $(u(\cdot), v(\cdot))$. Proposition 4.1 shows that $(x(\cdot), y(\cdot))$ is indeed a solution of system (1) for a measurable control $u(\cdot)$, but one cannot guarantee a priori that $z(\cdot)$ can be generated by a measurable control $v(\cdot)$, which is irrelevant for our purpose.

<footnote n="1">¹ A set-valued map $F: \mathcal{X} \rightsquigarrow \mathcal{X}$ is upper semi-continuous at $\xi \in \mathcal{X}$ if and only if for any neighborhood $\mathcal{N}$ of $F(\xi)$, there exists $\eta > 0$ such that for any $\xi' \in B_{\mathcal{X}}(\xi, \eta)$ one has $F(\xi') \subset \mathcal{N}$ (see for instance [1]).</footnote>
We end this section by exhibiting an approximation scheme from below of the optimal cost. These approaches are of major interest for minimization problems because, since upper bounds are commonly obtained via any sub-optimal control of problem $P_0, P_1, P_2$ or $P_3$ (provided typically by a numerical scheme), they are useful to frame the optimal value of the problem. This will be illustrated in Section 5.

Let us consider the family of dynamics parameterized by $\theta > 0$

$$ \begin{cases} \dot{x} = f(x, y, u), \\ \dot{y} = g(x, y, u), \\ \dot{z} = h_{\theta}(x, y, z, u, v), \end{cases} \qquad (12) $$

with

$$ h_{\theta}(x, y, z, u, v) = \max(g(x, y, u), 0)(1 - v e^{-\theta \max(y-z, 0)}). $$

Here the expression $e^{-\theta \max(y-z, 0)}$ plays the role of an approximation of $1_{\mathbb{R}^+}(z - y)$ when $\theta$ tends to $+\infty$.

We then define the family of Mayer problems

$$ \mathcal{P}_3^{\theta}: \inf_{\Pi(\cdot) \in \mathcal{S}_{\theta}} z(T) $$

where $\mathcal{S}_{\theta}$ denotes the set of absolutely continuous solutions $\Pi(\cdot) = (x(\cdot), y(\cdot), z(\cdot))$ of (12) for the initial condition $\Pi(0) = \Pi_0$. Let us underline that, for problems with Lipschitz dynamics and without state constraints, necessary conditions based on Pontryagin Maximum Principle can be derived, leading to shooting methods that are known to be very accurate. They can be initialized from numerical solutions of problems $P_1$ or $P_2$, that in turn can be obtained, for instance, through direct methods.

**Proposition 4.2** Under Assumption 2, consider any increasing sequence of numbers $\theta_n$ ($n \in \mathbb{N}$) that tends to $+\infty$. For each $n \in \mathbb{N}$, the problem $\mathcal{P}_3^{\theta_n}$ admits an optimal solution, and for any sequence of optimal solutions $(x_n(\cdot), y_n(\cdot), z_n(\cdot))$ of $\mathcal{P}_3^{\theta_n}$, the sequence $(x_n(\cdot), y_n(\cdot))$ converges, up to sub-sequence, uniformly to an optimal solution $(x^*(\cdot), y^*(\cdot))$ of Problem $\mathcal{P}$, and its derivatives weakly to $(\dot{x}^*(\cdot), \dot{y}^*(\cdot))$ in $L^2$. Moreover, $z_n(T)$ is an increasing sequence that converges to $\max_{t \in [0,T]} y^*(t)$.

*Proof* As in the proof of Proposition 4.1, we consider for any $\theta > 0$ the convexified dynamics

$$ \begin{cases} \dot{x} = f(x, y, u), \\ \dot{y} = g(x, y, u), \\ \dot{z} = h_{\theta}^{\dagger}(x, y, z, u, v, \alpha) := (1 - \alpha)h_{\theta}(x, y, z, u, v) + \alpha \max_{w \in U} h_{\theta}(x, y, z, w, 0), \end{cases} $$

where $\alpha \in [0, 1]$. Then, there exists an absolutely continuous solution $(x_{\theta}^{*}(\cdot), y_{\theta}^{*}(\cdot), z_{\theta}^{*}(\cdot))$ and a measurable control $(u_{\theta}^{*}(\cdot), v_{\theta}^{*}(\cdot), 0)$ that minimize $z(T)$. For the control $(u_{\theta}^{*}(\cdot), v_{\theta}^{*}(\cdot), 0)$, the solution is given by $(x_{\theta}^{*}(\cdot), y_{\theta}^{*}(\cdot), \tilde{z}_{\theta}^{*}(\cdot))$ where $\tilde{z}_{\theta}^{*}(\cdot)$ is solution of the Cauchy problem

$$ \dot{z} = \tilde{l}_{\theta}(t, z) := h_{\theta}^{\dagger}(x_{\theta}^{*}(t), y_{\theta}^{*}(t), z; u_{\theta}^{*}(t), v_{\theta}^{*}(t), 0), \quad z(0) = y(0) $$

while $z_{\theta}^{*}(\cdot)$ is solution of

$$ \dot{z} = l_{\theta}(t, z) := h_{\theta}^{\dagger}(x_{\theta}^{*}(t), y_{\theta}^{*}(t), z, u_{\theta}^{*}(t), v_{\theta}^{*}(t), \alpha_{\theta}^{*}(t)), \quad z(0) = y(0). $$

One can check that the inequality

$$ \tilde{l}_{\theta}(t, z) \le l_{\theta}(t, z), \quad t \in [0, T], \quad z \in \mathbb{R} $$

is fulfilled, which gives by comparison of solutions of scalar ordinary differential equations (see for instance [18]) the inequality

$$ \tilde{z}_{\theta}^{*}(t) \le z_{\theta}^{*}(t), \quad t \in [0, T]. $$

We deduce that $(x_{\theta}^{*}(\cdot), y_{\theta}^{*}(\cdot), z_{\theta}^{*}(\cdot))$ is necessarily a solution of (12).

Let

$$ \bar{y} := \inf_{u \in \mathcal{U}} \left\{ \max_{t \in [0, T]} y(t); (x(\cdot), y(\cdot)) \text{ sol. of (1)} \right\}. $$
By Proposition 4.1, we know that there exists an optimal solution $(x(\cdot), y(\cdot), z(\cdot))$ of problem $\mathcal{P}_3$ such that $z(T) = \bar{y}$. Clearly, this solution belongs to $S_\theta$ for any $\theta$, and we thus get

$$z_{\theta}^{*}(T) \leq \bar{y}. \qquad (13)$$

Let

$$F_{\theta}(\Pi) := \bigcup_{(u,v) \in U \times [0,1]} \begin{bmatrix} f(x, y, u) \\ g(x, y, u) \\ h_{\theta}(x, y, z, u, v) \end{bmatrix}$$

and note that one has

$$\lim_{\theta \to +\infty} d(F_{\theta}(\Pi), F(\Pi)) = 0, \quad \Pi \in \mathcal{D} \times \mathbb{R}. \qquad (14)$$

Consider an increasing sequence of numbers $\theta_n$ ($n \in \mathbb{N}$), and denote $\Pi_n(\cdot) = (x_n(\cdot), y_n(\cdot), z_n(\cdot))$ an optimal solution of problem $\mathcal{P}_3^{\theta_n}$. Note that one has

$$S_{\theta_{n+1}} \subset S_{\theta_n} \cdots \subset S_{\theta_0}. \qquad (15)$$

Therefore, the sequence $\dot{\Pi}_n(\cdot)$ is bounded, and $\Pi_n(\cdot)$ as well. As $F$ is upper semi-continuous, we obtain that $\Pi_n(\cdot)$ converges uniformly on $[0, T]$, up to a sub-sequence, to a certain $\Pi^*(\cdot) = (x^*(\cdot), y^*(\cdot), z^*(\cdot))$ which belongs to $S_l$ (see for instance [7, Th. 3.1.7]). From property(15), we obtain that $z_n(T)$ is a non decreasing sequence that converges to $z^*(T)$, and from (13), we get passing at the limit

$$z^*(T) \leq \bar{y}.$$

On another hand, $(x^*(\cdot), y^*(\cdot), z^*(\cdot))$ belongs to $S_l$ and we get from Proposition 4.1 the inequality

$$z^*(T) \geq \bar{y}.$$

Therefore, one has $z^*(T) = \bar{y}$ and $(x^*(\cdot), y^*(\cdot), z^*(\cdot))$ is then an optimal solution of problem $\mathcal{P}_3$. From Proposition 4.1, we obtain that one has necessarily

$$z^*(T) = \max_{t \in [0,T]} y^*(t).$$

Finally, the sequence $(\dot{x}_n(\cdot), \dot{y}_n(\cdot))$ being bounded, it converges, up to a sub-sequence, weakly to $(\dot{x}^*(\cdot), \dot{y}^*(\cdot))$ in $L^2$ thanks to Alaoglu's Theorem.

# 5 Numerical Illustrations

Our aim is to illustrate the different formulations on problems for which the optimal solution is known. We have considered here only direct and Hamilton-Jacobi-Bellman numerical methods, and not indirect methods. As underlined for instance in [6], computational schemes based on the Maximum Principle of Pontryagin (such as the shooting method), although usually quite accurate, require an a priori knowledge of the structure of the solution as well as a good approximation of the adjoint state. Here, our objective is to show that solving numerically problems with maximum cost considering the equivalent formulations we propose can be done straightforwardly with existing software based on direct or Hamilton-Jacobi-Bellman methods, without having to provide the structure of the optimal solution.

## 5.1 A Particular Class of Dynamics

We consider dynamics of the form

$$(\Sigma): \begin{cases} \dot{x} = f(x), \\ \dot{y} = g(x, u). \end{cases} \quad x \in \mathbb{R}^n, u \in U,$$

**Proposition 5.1** A feedback control $x \mapsto \phi^*(x)$ such that

$$g(x, \phi^*(x)) = \min_{u \in U} g(x, u), \quad x \in \mathbb{R}^n$$

is optimal for problem $\mathcal{P}$.
Proof For a given $x_0$ in $\mathbb{R}^n$, let $x(\cdot)$ be the solution of $\dot{x} = f(x)$, $x(0) = x_0$ independently to the control $u(\cdot)$. For any control $u(\cdot)$, the corresponding solution $y(\cdot)$ verifies

$$y(t) = y(0) + \int_{0}^{t} g(x(\tau), u(\tau)) d\tau \geq y(0) + \int_{0}^{t} \min_{v \in U} g(x(\tau), v) d\tau, \quad t \geq 0. \quad (16)$$

Consider the (measurable) map

$$\vartheta(t, u) := g(x(t), u), \quad (t, u) \in [0, T] \times U.$$

From the existence of measurable selection of extrema of measurable maps (see, for instance, Theorem 2 in [5]), there exists a measurable control $u^*(\cdot)$ such that

$$\vartheta(t, u^*(t)) = \min_{v \in U} \vartheta(t, v), \quad t \in [0, T].$$

Let $y^*(\cdot)$ be defined as

$$y^*(t) := y(0) + \int_{0}^{t} g(x(\tau), u^*(\tau)) d\tau = y(0) + \int_{0}^{t} \min_{v \in U} g(x(\tau), v) d\tau, \quad t \geq 0.$$

Clearly, $y^*(\cdot)$ is solution of $\Sigma$ for the control $u^*(\cdot)$ and one gets from (16)

$$\max_{t} y(t) \geq \max_{t} y^{*}(t).$$

We conclude that $y^*(\cdot)$ is an optimal trajectory of problem $\mathcal{P}$ for the control generated by the feedback $\phi^*$.

As a toy example, we have considered the system

$$\begin{cases} \dot{x} = 1, & x(0) = 0, \\ \dot{y} = (1-x)(2-x)(4-x)(1+u/2), & y(0) = 0, \end{cases} \quad u \in [-1, 1]$$

for which

$$\phi^*(x) = -\operatorname{sign}((1-x)(2-x)(4-x))$$

is an optimal control which minimizes $\max_{t \in [0,T]} y(t)$. The optimal control is thus pure bang-bang. Remark that this problem can be equivalently written with a scalar non-autonomous dynamics

$$\dot{y} = (1 - t)(2 - t)(4 - t)(1 + u/2)$$

for which the open-loop control

$$u^*(t) = -\operatorname{sign}((1-t)(2-t)(4-t))$$

is optimal.

For $T = 5$, we have first computed the exact optimal solution of problem $\mathcal{P}$ with the open-loop $u^*(\cdot)$, by integrating the dynamics with Scipy in Python software (see Figure 2). Effects of perturbations on the switching times on the criterion are presented in Table 1, which show a quite high sensitivity of the optimal control for this problem (as it often the case for bang-bang controls). Then, we have solved numerically problems $P_0$ to $P_2$ with a direct method

<table><caption>**Table 1** Sensitivity to the optimal switching.</caption><thead><tr><th>disturbance</th><th>max<sub>t∈[0,T]</sub> y(t)</th><th>error</th></tr></thead><tbody><tr><td>0</td><td>2.24985</td><td>0</td></tr><tr><td>0.001%</td><td>2.24985</td><td>4.10<sup>-6%</sup></td></tr><tr><td>0.01%</td><td>2.25010</td><td>0.01%</td></tr><tr><td>0.1%</td><td>2.69457</td><td>20%</td></tr></tbody></table>

(Bocop software using Gauss II integration scheme) for 500 time steps and an optimization relative tolerance equal to
<figure><img src="image_3.png" /><figcaption>Fig. 2 Optimal solution: $y(\cdot)$ on the left, $u^*(\cdot)$ on the right</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3 Comparisons of the three methods on $y(\cdot)$, $u(\cdot)$, $z(\cdot)$ and $v(\cdot)$.</figcaption></figure>

$10^{-10}$. For problem $\mathcal{P}_3$, as the dynamics is not continuous, direct methods do not work well and we have used instead a numerical scheme based on dynamic programming (BocopHJB software) with 500 time steps and a discretization of 200 × 200 points of the state space. For the additional control $v$, we have considered only two possible values 0
<table><caption>**Table 2** Comparison of the numerical results.</caption><thead><tr><th>problem</th><th>max <i>y</i>(<i>t</i>)<br><i>t</i>∈[0,<i>T</i>]</th><th>error</th><th>computation time</th></tr></thead><tbody><tr><td><i>P</i></td><td>2.24705</td><td>0</td><td>—</td></tr><tr><td><i>P</i><sub>0</sub></td><td>2.249888</td><td>0.126%</td><td>0.5 <i>s</i></td></tr><tr><td><i>P</i><sub>1</sub></td><td>2.24998</td><td>0.130%</td><td>1.8 <i>s</i></td></tr><tr><td><i>P</i><sub>2</sub></td><td>2.249941</td><td>0.129%</td><td>3.8 <i>s</i></td></tr><tr><td><i>P</i><sub>3</sub></td><td>2.26778</td><td>0.8%</td><td>248 <i>s</i></td></tr></tbody></table>

and 1 as we know that the optimal solution is reached for $v \in \{0, 1\}$ (see Proposition 3.1). The numerical results and computation times are summarized in Table 2, while Figure 3 presents the corresponding trajectories.

We note that the direct method give very accurate results, and the computation time for problem $\mathcal{P}_0$ is the lowest because it has only one control. The computation time for problem $\mathcal{P}_2$ is slightly higher than for $\mathcal{P}_1$ because the mixed constraint $\mathcal{C}_m$ is heavier to evaluate. The numerical method for problem $\mathcal{P}_3$ is of completely different nature as it computes the optimal solution for all the initial conditions on the grid, which explains a much longer computation time. The accuracy of the results is also directly related to the size of the discretization grid and can be improved by increasing this size but at the price of a longer computation time.

It can be noticed in Figure 3 some differences between the obtained trajectories. Let us underline that after the peak of $y(\cdot)$, there is no longer uniqueness of the optimal control (indeed any value of the control does not modify the value of the peak). Hence, we may think about considering an additional cost after the peak to avoid this multiplicity. However, since neither the value of the peak, nor the time when it is reached, are known, this does not seem easy to implement.

## 5.2 Application to an Epidemiological Model

The SIR model is one of the most basic transmission model in epidemiology for a directly transmitted infectious disease (see [19] for a complete introduction) and it retakes great importance nowadays due to Covid-19 epidemic.

Consider on a time horizon $[0, T]$ variables $S(t)$, $I(t)$ and $R(t)$ representing the fraction of susceptible, infected and recovery individuals at time $t \in [0, T]$, so that one has $S(t) + I(t) + R(t) = 1$ with $S(t)$, $I(t)$, $R(t) \ge 0$. Let $\beta > 0$ be the rate of transmission and $\gamma > 0$ the recovery rate. Interventions as lock-downs and curfew are modeled as a factor in rate transmission that we denote $u$ and which represents our control variable taking values in $[0, u_{max}]$ with $u_{max} \in (0, 1)$, where $u = 0$ means no intervention and $u = u_{max}$ the most restrictive one which reduces as much as possible contacts among population. The SIR dynamics including the control is then given by the following equations:

$$ \dot{S} = -(1-u)\beta SI, \qquad (17) $$

$$ \dot{I} = (1-u)\beta SI - \gamma I, \qquad (18) $$

$$ \dot{R} = \gamma I. \qquad (19) $$

When the reproduction number $R_0 = \beta/\gamma$ is above one and the initial proportion of susceptible is above the herd immunity threshold $S_h = R_0^{-1}$, it is well known that there is an epidemic outbreak. Then, the objective is to minimize the peak of the prevalence

$$ \max_{t \in [0, T]} I(t) $$

with respect to control $u(\cdot)$ subject to a $L^1$ budget

$$ \int_{0}^{T} u(t) \le Q \qquad (20) $$

on a given time interval $[0, T]$ where $T$ is in chosen large enough to ensure the herd immunity of the population is reached at date $T$. Note that one can drop the $R$ dynamics to study this problem. If the constraint (20) were not imposed, then the optimal solution would be the trivial control $u(t) = u_{max}$, $t \in [0, T]$, which is in general unrealistic from a operational point of view. A similar problem has been considered in [16] but under the constraint
that intervention occurs only once on a time interval of given length, that we relax here. Note that the constraint (20) can be reformulated as a target condition, considering the augmented dynamics

$$ \dot{S} = -(1-u)\beta SI, \qquad (21) $$

$$ \dot{I} = (1-u)\beta SI - \gamma I, \qquad (22) $$

$$ \dot{C} = -u(t), \qquad (23) $$

with initial condition $C(0) = Q$ and target $\{C \ge 0\}$. Extension of the results of Sections 3 and 4 to problems with target do not present any particular difficulty, and is left to the reader.

For initial conditions $I_0 = I(0) > 0$ and $S_0 = S(0) > S_h$, the optimal solution has been determined in [13] as the feedback control

$$ \psi(I, S) := \begin{cases} 1 - \frac{S_h}{S} & \text{if } I = \bar{I} \text{ and } S > S_h, \\ 0 & \text{otherwise,} \end{cases} $$

where

$$ \bar{I} := \frac{I_0 + S_0 - S_h - S_h \log \left( \frac{S_0}{S_h} \right)}{Q \beta S_h + 1} $$

is the optimal value of the peak. The proof of the optimality of this feedback is out of the scope of the present paper and can be found in [13]. This control strategy consists in three phases:

1. no intervention until the prevalence $I$ reaches $\bar{I}$ (null control),

2. maintain the prevalence $I$ equal to $\bar{I}$ until $S$ reaches $S_h$ (singular control),

3. no longer intervention when $S > S_h$ (null control),

as illustrated in Figure 4 for the parameters given in Table 3. Note that differently to the previous example, this

<figure><img src="image_5.png" /><figcaption>Fig. 4 The optimal solution for the SIR problem.</figcaption></figure>

control strategy is intrinsically robust with respect to a bad choice of $\bar{I}$: the maximum value of $I$ is always guaranteed to be equal to $\bar{I}$. However, a mischoice of $\bar{I}$ has an impact on the budget (see [13] for more details).

<table><caption>Table 3 Parameters used in numerical computations and optimal value of the peak.</caption><thead><tr><th>β</th><th>γ</th><th>T</th><th>Q</th><th>S(0)</th><th>I(0)</th><th>II</th></tr></thead><tbody><tr><td>0.21</td><td>0.07</td><td>300</td><td>28</td><td>1 - 10<sup>-6</sup></td><td>10<sup>-6</sup></td><td>0.105</td></tr></tbody></table>

Adding the z-variable, we end up with a dynamics in dimension four, which is numerically heavier than for the previous example. In particular, methods based on the value function are too time consuming to obtain accurate results for refined grids in a reasonable computation time. So we have considered direct methods only. We do not consider here problem $P_3$, but instead its regular approximations $P_3^\theta$ suitable to direct methods. For direct methods that use
algebraic differentiation of the dynamics, convergence and accuracy are much better if one provides differentiable dynamics. This is why we have approximated the max(·, 0) operator for problems $P_1$ and $P_2$ by the Laplace formula

$$ \frac{\log(e^{\lambda\xi} + 1)}{\lambda} \underset{\lambda \to +\infty}{\longrightarrow} \max(\xi, 0), \quad \xi \in \mathbb{R} $$

with $\lambda = 100$ for the numerical experiments. For problem $P_3^\theta$, one has to be careful about the interplay between the approximations of $\max(\cdot, 0)$ and the sequence $\theta_n \to +\infty$, to provide approximations from below of the optimal value. The function $h_\theta$ is thus approximated by the expression

$$ h_{\theta}(x, y, z, u, v) \approx \frac{\log(e^{\lambda_1 g(x, y, u)} + 1)}{\lambda_1} \left(1 - v e^{\frac{\theta}{\lambda_2} \log(e^{\lambda_2(y-z)} + 1)}\right) $$

which depends on three parameters $\lambda_1$, $\lambda_2$ and $\theta$. Posit for convenience

$$ \alpha := \frac{\theta}{\lambda_2} $$

and consider the function

$$ \omega_{\alpha, \lambda_2}(\xi) := e^{-\alpha \log(e^{-\lambda_2 \xi} + 1)}, \quad \xi \in \mathbb{R} $$

which approximates the indicator function $1_{\mathbb{R}^+}$.

**Lemma 5.1** One has the following properties.

1. For any positive numbers $\alpha$, $\lambda_2$, the function $\omega_{\alpha, \lambda_2}$ is increasing with

$$ \lim_{\xi \to -\infty} \omega_{\alpha, \lambda_2}(\xi) = 0, \quad \lim_{\xi \to +\infty} \omega_{\alpha, \lambda_2}(\xi) = 1. $$

2. For any $\varepsilon \in (0, 1)$, one has $\omega_{\alpha, \lambda_2}(-\varepsilon^2) = \varepsilon$ and $\omega_{\alpha, \lambda_2}(0) = 1 - \varepsilon$ exactly for

$$ \alpha = -\frac{\log(1 - \varepsilon)}{\log(2)}, \quad \lambda_2 = \frac{\log(\varepsilon^{-\frac{1}{\alpha}} - 1)}{\varepsilon^2}. \qquad (24) $$

*Proof* One has first

$$ \omega'_{\alpha, \lambda_2}(\xi) = \lambda_2 \alpha \frac{e^{-\lambda_2 x}}{e^{-\lambda_2 \xi} + 1} \omega_{\alpha, \lambda_2}(\xi) > 0 $$

and the function $\omega_{\alpha, \lambda_2}(\cdot)$ is thus increasing. From

$$ \lim_{\xi \to -\infty} -\alpha \log(e^{-\lambda_2 \xi} + 1) = -\infty $$

one get

$$ \lim_{\xi \to -\infty} \omega_{\alpha, \lambda_2}(\xi) = 0 $$

and similarly

$$ \lim_{\xi \to +\infty} -\alpha \log(e^{-\lambda_2 \xi} + 1) = 0 $$

implies

$$ \lim_{\xi \to +\infty} \omega_{\alpha, \lambda_2}(\xi) = 1. $$

Finally, with simple algebraic manipulation of the conditions $\omega_{\alpha, \lambda_2}(-\varepsilon^2) = \varepsilon$ and $\omega_{\alpha, \lambda_2}(0) = 1 - \varepsilon$, one obtains straightforwardly the expressions (24).
<table><caption>**Table 4** Values of parameters α, λ₂ for different ε.</caption><thead><tr><th>ε</th><th>α</th><th>λ<sub>2</sub></th></tr></thead><tbody><tr><td>0.2</td><td>0.32</td><td>124</td></tr><tr><td>0.15</td><td>0.234</td><td>360</td></tr><tr><td>0.1</td><td>0.152</td><td>1514</td></tr><tr><td>0.075</td><td>0.112</td><td>4094</td></tr><tr><td>0.05</td><td>0.074</td><td>16193</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Approximation of the indicator function with different values of ε (zoom on the abscissa axis on the right).</figcaption></figure>

<table><caption>**Table 5** Comparison of performances for problems P₀, P₁, P₂.</caption><thead><tr><th>problem</th><th>max<sub>t∈[0,T]</sub> y(t)</th><th>computation time</th></tr></thead><tbody><tr><td>P<sub>0</sub></td><td>0.1015</td><td>10 s</td></tr><tr><td>P<sub>1</sub></td><td>0.1015</td><td>12 s</td></tr><tr><td>P<sub>2</sub></td><td>0.1015</td><td>13 s</td></tr></tbody></table>

<table><caption>**Table 6** Comparison of performances for problem P<sub>3</sub><sup>θ</sup>.</caption><thead><tr><th>ε</th><th>z(T)</th><th>max<sub>t∈[0,T]</sub> y(t)</th><th>computation time</th></tr></thead><tbody><tr><td>0.2</td><td>0.0684</td><td>0.1038</td><td>80 s</td></tr><tr><td>0.15</td><td>0.0823</td><td>0.1038</td><td>65 s</td></tr><tr><td>0.1</td><td>0.0954</td><td>0.1037</td><td>51 s</td></tr><tr><td>0.075</td><td>0.0993</td><td>0.1050</td><td>83 s</td></tr><tr><td>0.05</td><td>0.1010</td><td>0.1036</td><td>97 s</td></tr></tbody></table>

We have taken λ₁ = 5000 and considered a sequence of approximations of the indicator function for the values given in Table 4 according to expressions (24) of Lemma 5.1 (see Figure 5).

Computations have been performed with Bocop software on a standard laptop computer (with a Gauss II integration scheme, 600 time steps and relative tolerance 10<sup>-10</sup>). As one can see in Figure 6 and Table 5 problems P₀, P₁, P₂ give the peak values with a very good accuracy, and present similar performances in terms of computation time. In Figure 7 and Table 6, the numerical solutions of P<sub>3</sub><sup>θ</sup> are illustrated for the values of α and λ₂ given in Table 4. As expected, the numerical computation of the family of problems P<sub>3</sub><sup>θ</sup> provides an increasing sequence of approximation from below of the optimal value and thus complements the computation of problems P₀, P₁ or P₂. From Figures of Tables 5 and 6, one can safely guarantee that the optimal value belongs to the interval [0.1010, 0.1015]. However, the trajectories found for P<sub>3</sub><sup>θ</sup> are not as closed as the ones of problems P₀, P₁ or P₂. This can be explained by the fact that problems P<sub>3</sub><sup>θ</sup> are not subject to the constraint z(t) ≥ y(t) and thus provides trajectories for which z(T) is indeed below max<sub>t</sub> y(t).
<figure><img src="image_7.png" /><figcaption>Fig. 6 Comparisons of numerical results for the methods $P_0, P_1, P_2$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7 Comparison of the numerical results for problem $P_3^\theta$.</figcaption></figure>

Finally, we have compared our approximation technique with the classical approximation of the $L^\infty$ criterion by $L^p$ norms

$$ \mathcal{P}_{L^p}: \inf_{u(\cdot) \in \mathcal{U}} \|y(t)\|_p $$

with the same direct method. To speed up the convergence, we have used the Bocop facility which allows a batch mode which consists in initializing the search from a solution found for a former value of $p$, that have been taken $p \in \{2, 5, 10, 15\}$ (see Figure 8). Besides, to ensure convergence it was necessary take 1200 time step instead of 600 as in previous simulations. The total time of the process is 78s after summing computation times given in Table 7. However, one can see that the trajectory found for $p = 15$ is quite far to give a peak value close from the other
<figure><img src="image_9.png" /><figcaption>Fig. 8 Numerical solutions for problems $P_{LP}$.</figcaption></figure>

<table><caption>Table 7 Comparison of the numerical results with the $L^p$ approximation.</caption><thead><tr><th>p</th><th>max<sub>t∈[0,T]</sub> y(t)</th><th>||y(t)||<sub>p</sub></th><th>computation time</th></tr></thead><tbody><tr><td>2</td><td>0.119653</td><td>1.0222</td><td>34 s</td></tr><tr><td>5</td><td>0.105244</td><td>0.2474</td><td>14 s</td></tr><tr><td>10</td><td>0.105375</td><td>0.15678</td><td>13 s</td></tr><tr><td>15</td><td>0.105170</td><td>0.13549</td><td>17 s</td></tr></tbody></table>

methods. Moreover, the same method for $p = 15$ but initialized from the solution found for $p = 2$ gives poor results for a computation time of 50s (see Figure 9). We conclude that the $L^p$ approximation is not practically reliable for this kind of problems.

<figure><img src="image_10.png" /><figcaption>Fig. 9 Numerical solution for $P_{L15}$ without batch iteration (computation time 50s).</figcaption></figure>

## 6 Discussion and Conclusions

In this work, we have presented different reformulations of optimal control problems with maximum cost in terms of extended Mayer problems, and tested them numerically on two examples whose optimal solution has bang-bang
controls and singular arcs. We have proposed two kinds of formulations: one with state or mixed constraints suitable to direct methods, and another one without any constraint but less regular and suitable to dynamical programming type methods. Moreover, for the latter one, we have proposed an approximation scheme generated by a sequence of regular Mayer unconstrained problems, which performs better than approximations based on $L^p$ norms. However, although this second approach requires larger computation time, it complements the first one providing approximations of the optimal value from above.

This first work puts in perspective the study of necessary optimality conditions for the maximum cost problems with the help of these formulations, which will be the matter of a future work.

Finally, we summarize advantages and drawbacks of the different formulations for numerical computations in Table 8, that could help practitioners in the choice of the method.

<table><caption>**Table 8** Comparison of the different formulations.</caption><thead><tr><th>Formulation</th><th>P<sub>0</sub></th><th>P<sub>1</sub> or P<sub>2</sub></th><th>P<sub>3</sub></th><th>P<sub>3</sub><sup>θ</sup></th></tr></thead><tbody><tr><td>suitable to direct methods</td><td>yes</td><td>yes</td><td>no</td><td>yes</td></tr><tr><td>suitable to Hamilton-Jacobi-Bellman methods</td><td>no</td><td>yes</td><td>yes</td><td>yes</td></tr><tr><td>suitable to shooting methods without constraint</td><td>no</td><td>no</td><td>no</td><td>yes</td></tr><tr><td>provides approximations from below</td><td>no</td><td>no</td><td>no</td><td>yes</td></tr></tbody></table>

**Acknowledgements** The authors are grateful to Pierre Martinon for fruitful discussions and advices. This work was partially supported by ANID-PFCHA/Doctorado Nacional/2018-21180348, FONDECYT grant 1201982 and Centro de Modelamiento Matemático (CMM) FB210005 BASAL fund for center of excellence, all of them from ANID (Chile).

## References

1. Aubin, J.-P., Viability Theory, Springer (2009).

2. Barron, E. N., The Pontryagin maximum principle for minimax problems of optimal control, Nonlinear Analysis: Theory, Methods & Applications 15, 1155-1165 (1990).

3. Barron, E.N. and Ishii, H. The Bellman equation for minimizing the maximum cost, Nonlinear Analysis: Theory, Methods & Applications, 13(9), 1067-1090 (1989).

4. Barron, E. N., Jensen, R. R. and Liu, W., The $L^\infty$ control problem with continuous control functions, Nonlinear Analysis: Theory, Methods & Applications 32, 1-14 (1998).

5. Brown, L. and Purves, R., Measurable selections of extrema, The Annals of Statistics, 1(5), 902-912 (1973).

6. Caillau, J.B., Ferretti, R., Trélat, E. and Zidani, H., Numerics for finite-dimensional optimal control problems. Preprint hal-03707475 (2022).

7. Clarke, F., Optimization and Nonsmooth Analysis, SIAM Classics in Applied Mathematics (1990).

8. Di Marco, A. and Gonzalez, R. L. V., A numerical procedure for minimizing the maximum cost, in: System Modelling and Optimization (Prague, 1995), Chapman & Hall, London, 285-291 (1996).

9. Di Marco, A. and Gonzalez, R. L. V., Minimax optimal control problems. Numerical analysis of the finite horizon case, ESAIM: Mathematical Modelling and Numerical Analysis 33(1), 23-54 (1999).

10. Fleming, W. and Rishel, R., Deterministic and Stochastic Optimal Control, Springer-Verlag, New-York (1975).

11. Gianatti, J.; Aragone, L.; Lotito, P. and Parente, L., Solving minimax control problems via nonsmooth optimization, Operations Research Letters 44, 680-686 (2016).

12. Gonzalez, R.L.V. and Aragone, L., A Bellman's equation for minimizing the maximum cost, Indian Journal of Pure & Applied Mathematics 31(12), 1621-1632 (2000).

13. Molina, E. and Rapaport, A., An optimal feedback control that minimizes the epidemic peak in the SIR model under a budget constraint, Automatica, in press (2022).

14. Morgan, J. and Peet, M., Extensions of the Dynamic Programming Framework: Battery Scheduling, Demand Charges, and Renewable Integration, IEEE Transactions on Automatic Control, 66(4), 1602-1617 (2021).

15. Morgan, J. and Peet, M., A generalization of Bellman's equation with application to path planning, obstacle avoidance and invariant set estimation, Automatica, 127, n. 109510 (2021).

16. Morris, D., Rossine, F., Plotkin, J. and Levin, S., Optimal, near-optimal, and robust epidemic control, Communications Physics 4(78) (2021).

17. Tao, T., An introduction to measure theory, Graduate Studies in Mathematics 126, American Mathematical Society (2011).

18. Walter, W., Ordinary Differential Equations, Springer New-York (1998).

19. Weiss, H., The SIR model and the foundations of public health, MATerials MATemàtics, 2013(3) (2013).