Jorge Amaya

Cristopher Hermosilla

Giorgia Allen

Chourouk, Alexiane, Alexandre, Remi, Matthieu Antoine Anouk

Po-Yi Ludovic

Eugenio David

Dando

Abuelita Menche

Tía Sonia

Bea Tío Pez

Watón Raúl

APPLICATION OF OPTIMAL CONTROL TECHNIQUES TO NATURAL SYSTEMS MANAGEMENT

Keywords: optimal control, natural system management, mining, epidemiology Résumé contrôle optimal, gestion du système naturel, exploitation minière, épidémiologie

Muchas cosas han

Chapter 1 Introduction

For decades optimal control theory has been an important tool to apply in real world problems and engineering. A main motivation to develop this theory came from aerospace engineering in the context of the aerospace race in the Cold War, to solve problems such as the minimization of the fuel or time spent by an aircraft going from an initial position to another one. Many other examples exist in domains of aeronautics, chemical processes, vehicles and transportation, energy, biological systems, etc. Some of them are showed in [START_REF] Lee | Foundations of optimal control theory[END_REF][START_REF] Sethi | Optimal Control Theory[END_REF][START_REF] Trélat | Contrôle optimal: théorie & applications[END_REF][START_REF] Vinter | Optimal control[END_REF].

This work consists of applying analytical as well as numerical optimal control techniques to two specific problems. The first one corresponds to the Open Pit problem, well-known in context of mining, but using a new continuous formulation in contrast to the classical binary formulation. The second one arises during the pandemic motivated by overcrowded hospitals due to many infected individuals at the same time, therefore, the idea was to use optimal control tools in order to minimize the peak of infections.

This introduction is devoted to presenting the main optimal control techniques used and at the end, a brief explanation of each chapter. We refer to Chapters 2 and 5 as introductions to Parts I and II respectively.

General background in optimal control

Consider a time interval [0, T ], and (x(t), u(t)) ∈ R n × R m ∀t ∈ [0, T ] denoting x for the state and u for the control. An optimal control problem in a Bolza form is written as follows:

(P 0 )                Minimize ϕ(x(T )) + T 0 f 0 (t, x(t), u(t))dt

over all

x ∈ A C n [0, T ] and u ∈ U satisfying ẋ(t) = f (t, x(t), u(t)) for a.e.

x(0) = x 0 , x(T ) ∈ E

where the set of admissible controls is usually U = {u : [0, T ] → U, u ∈ L ∞ }. In this formulation, boundary conditions are represented by a fixed initial state x(0) = x 0 and a final state x(T ) belonging to a set E which typically correspond to {y ∈ R n |c 1 (y) = 0, c 2 (y) ≤ 0} where c 1 and c 2 are continuously differentiable functions with values in R k 1 and R k 2 respectively. Additionally, we can consider optimal control problems with constraints in the state written as (t, x(t)) ∈ A.

We then say that it is an optimal control problem with state constraints. Typically, state constraints are represented by g(x(t)) ≤ 0,t ∈ [0, T ] with g regular enough. In this same way, we define an optimal control problem with mixed constraints when we add (P 0 ) to the expression:

(t, x(t), u(t)) ∈ A

We do not focus on existence results for our problems, however those interested can refer to [START_REF] Cesari | Optimization-theory and applications: problems with ordinary differential equations[END_REF][START_REF] Clarke | Functional analysis, calculus of variations and optimal control[END_REF][START_REF] Lee | Foundations of optimal control theory[END_REF][START_REF] Trélat | Contrôle optimal: théorie & applications[END_REF]. In the following sections we will briefly recall the principle of 3 main classes of methods in optimal control, furthermore, we present Bocop solver, used in numerical experiments of this thesis.

Pontryagin Maximum Principle and indirect shooting methods

After its discovery (see [START_REF] Gamkrelidze | Discovery of the maximum principle[END_REF] for a review of the history), many versions of Maximum Principle can be found in the literature. Here we present one of them, but other versions are shown throughout this document when necessary. Firstly, we introduce the costate p, of the same dimension as the state x and define the Hamiltonian associated to (P 0 ) as H(t, x, p, p 0 , u) = p 0 f 0 (t, x, p) + p T f (t, x, u).

Let f 0 , f , ϕ be C 1 class and u ∈ L ∞ , the Maximum Principle is stated as follow Theorem 1.1.1. Under previous assumptions, if the pair (x(t), u(t)) is optimal for (P 0 ) then there exists an absolutely continuous function p : [0, T ] → R n and a real value p 0 ≤ 0 such that:

1. (p(•), p 0 ) ̸ = 0.

2. ẋ(t) = ∂ p H(t, x(t), p(t), p 0 , u(t)), ṗ(t) = ∂ x H(t, x(t), p(t), p 0 , u(t))

3. H(t, x(t), p(t), p 0 , u(t)) = max v∈U H(t, x(t), p(t), p 0 , v) [START_REF] Amaya | Optimality conditions for the continuous model of the final open pit problem[END_REF]. p(T ) ∈ p 0 ∇ϕ(x(T )) + N E (x(T ))

5. d dt H(t, x(t), p(t), p 0 , u(t)) = ∂ t H(t, x(t), p(t), p 0 , u(t))

A proof of this theorem can be seen in [START_REF] Cesari | Optimization-theory and applications: problems with ordinary differential equations[END_REF][START_REF] Lee | Foundations of optimal control theory[END_REF]. Note that this base version does not consider state constraints although these will appear in the Final Open Pit problem, therefore, we will present more general versions where appropriate.

Shooting method is based on the Pontryagin Maximum Principle and it is part of the indirect methods (for more detail see [START_REF] Ascher | Numerical solution of boundary value problems for ordinary differential equations[END_REF][START_REF] Bulirsch | Introduction to numerical analysis[END_REF][START_REF] Cesari | Optimization-theory and applications: problems with ordinary differential equations[END_REF]). To explain this method, assume the normal case, i.e., p 0 ̸ = 0. and free final conditions on x, that means, E = R n . The idea of the shooting method is to solve the boundary value problem derived from the Maximum Principle:

(BV P)          ẋ(t) = f (t, x(t), ū(t)) ṗ(t) = -∂ x H(t, x (
t), p(t), p 0 , ū(t)) x(0) = x 0 p(T ) = ∇ϕ(x(T )) using the control ū obtained in point 3 of theorem 1.1.1 when it can be written depending explicitly on x(t), p(t), i.e. ū(t) = Ψ(x(t), p(t)) for a certain function Ψ. We introduce the shooting function associated to this system as the functional that returns p(T ) -∇ϕ( x(T )) for a value z, where ( x, p) are solutions of the following Initial value problem

(IV P z )         
ẋ(t) = f (t, x(t), Ψ(x(t), p(t))) ṗ(t) = -∂ x H(t, x(t), p(t), p 0 , Ψ(x(t), p(t))) x(0) = x 0 p(0) = z .

If we call S : z → p(T ) -∇ϕ( x(T )) the shooting function, then, the shooting method consists in solving S(z) = 0, which is equivalent to solving the system (BVP) and therefore, obtaining a pair (x, p) fulfilling the conditions of theorem 1.1.1.

This idea can be generalized to a different set E and a more general optimal control problem, for instance, when x 0 is not fix but belongs to a set F ⊆ R n . In a more general case state constraints can be considered but the generalization is far more difficult.

Direct transcription approach

The so-called direct approach transforms the infinite dimensional optimal control problem (P 0 ) into a finite dimensional optimization one (typically non linear). The idea is to discretize the time interval [0, T ] in {t 0 = 0, . . . ,t N = t f } and apply it to the state and control variables, obtaining the finite dimensional variable X = {x 0 , . . . , x N , u 0 , . . . , u N-1 }.

Several options exist to discretize the Ordinary Differential Equation (ODE) of x, for instance, using Euler explicit form and considering an equidistant discretization on time, where we obtain

x i+1 = x i + h f (x i , u i ).
Finally, we get the following nonlinear programming problem:

(NLP)            min ϕ(x N ) + N-1 ∑ i=0 h f 0 (t i , x i , u i ) s.a x i+1 = x i + h f (x i , u i ) ∀i ∈ {0, .., N -1} u i ∈ U ∀i ∈ {0, .., N -1} x 0 = x0 , x N ∈ E
State and mixed constraints can be incorporated directly in each discretization point. These methods are widely used in industrial applications because they are more straightforward to apply than indirect methods. We refer the reader to [START_REF] Betts | Practical methods for optimal control using nonlinear programming[END_REF] and [START_REF] Nocedal | Numerical optimization[END_REF] for more details on direct transcription methods and NLP algorithms.

Hamilton-Jaccobi-Bellman equation

To introduce this results it is necessary to define the value function associated to (P 0 ). Then, for t ∈ [0, T ] and ξ ∈ R n , the value function V (t, ξ ) is defined by: V (t, ξ ) = min u∈U ϕ(x(T )) + T t f 0 (s, x(s), u(s))ds | ẋ(s) = f (s, x(s), u(s)) ∀s ∈ [t, T ]; x(t) = ξ .

It corresponds to the optimal value of (P 0 ) taking the initial time and position as t and ξ respectively. The main result is then Theorem 1.1.2. The value function V (t, ξ ) is a solution (in some specific sense) of the equation

∂ t V + min v∈U H(t, ξ , ∂ ξ V, -1, v) = 0, ∀(t, ξ ) ∈ (0, T ) × R n v(T, ξ ) = ϕ(ξ )
An important tool in finding the value function is the Dynamic Programming Principle which corresponds to the following expression V (t, ξ ) = min u∈U τ t f 0 (s, u(s), y(s))ds +V (y(τ), τ) .

(1.1)

Many methods mix the HJB equation and dynamic principle giving a global solution unlike direct transcription which is local, but the computational effort is normally higher.

For a complete introduction to the subject see [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF][START_REF] Sontag | Mathematical control theory: deterministic finite dimensional systems[END_REF], considering state-constraints [START_REF] Altarovici | A general Hamilton-Jacobi framework for nonlinear state-constrained control problems[END_REF][START_REF] Capuzzo-Dolcetta | Hamilton-Jacobi equations with state constraints[END_REF] and the link with Maximum Principle [START_REF] Clarke | The relationship between the maximum principle and dynamic programming[END_REF]. This approach is not mainly used in this thesis except for numerical simulations using the HJB version of Bocop, the solver that we are going to present now.

Bocop toolbox

In most of the numerical simulation of this thesis we use Bocop, an open source toolbox for optimal control problems, developed by Inria and which can be downloaded directly from https://www.bocop. org/. On the official web-page there is the option to choose between two different packages.

The first and original package implements a local optimization method. The optimal control problem is approximated by a finite dimensional optimization problem (NLP) using a time discretization (the direct transcription approach, see section 1.1.2). The NLP problem is solved by the well known software Ipopt [START_REF] Waechter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF], using sparse exact derivatives computed by CppAD [START_REF] Bell | CppAD: a package for C++ algorithmic differentiation[END_REF] 1 . In Figure 1.1 its interface is shown and it is possible to observe several options and functionalities of Bocop.

The second package BocopHJB implements a global optimization method. Similarly to the Dynamic Programming approach, the optimal control problem is solved in two steps. First we solve the Hamilton-Jacobi-Bellman equation satisfied by the value function of the problem. Then we simulate the optimal trajectory from any chosen initial condition. The computational effort is essentially taken by the first step, whose result, the value function, can be stored for subsequent trajectory simulations. 2 . Figure 1.2 shows interface of BocopHJB.

Bocop runs under Linux, Mac and Windows, besides, in both interfaces it is possible to export solutions in a specified data format allowing its treatment in Matlab or Python. We also take advantage of their interpolation tools, mainly in the mining context, where we have discretized data to be incorporated in the continuous optimal control model. A well detailed user guide for both packages can be found in https://www.bocop.org/download/.

Dissertation Outline

In this section we present how the manuscript is organized and the main subjects of each chapter.

Part I of this thesis begins with Chapter 2 where we introduce the Open Pit problem in mining. The focus of this chapter is to show the classical binary programming formulation introduced by Johnson in 1968 and the newer formulation introduced by Alvarez et al. in 2011. The last one, which uses continuous functions to retrieve the profile of a mine, will be the seminal model for Chapters 3 and 4.

In Chapter 3 we reformulate the continuous model of Alvarez et al, looking for optimal profiles among absolutely continuous functions rather than continuous ones. This allowed us to write the Liptchitz Figure 1.1: Bocop interface and some options modulus L p (x), introduced in equation (2.1), as a derivative of p. Using this slight reduction in the functional space, we can use optimal control tools to obtain optimality conditions for profiles in 2D and 3D. This is the main result of this chapter and corresponds to theorems 3.2.1 and 3.3.1.

On the other hand, Chapter 4 focuses on a numerical study of the Open Pit problem and its Sequential version. In this chapter we propose a new semi-continuous formulation, presented in section 4.2.2, of the Sequential Open Pit problem. This formulation is based on a discretization of the space domain and then the profile is represented by a finite set of variables at the discretization nodes. The explicit models (SOP) 2D S C and (SOP) 3D S C are shown in section 4.2. We present in this chapter several numerical experiments using local and global methods (using Bocop and BocopHJB) for the two dimensional case and for the first time in literature, a numerical solution for a 3 dimensional case in a continuous framework using this new semi-continuous model (see section 4.4.3).

Moving to Part II, in Chapter 5 we introduce two compartmental models for epidemiological purposes. The first one corresponds to the classical SIR model, which we used in Chapters 6 and 7. The second one is a small contribution developed in the context of a OPS (Organización Panamericana de la salud) project where we set a compartmental model that considers vaccination ans booster status (see section 5.2).

In Chapter 6 we work with optimal control problems consisted of minimizing the maximum of a state. We provide four alternative formulations to this problem in the Mayer form, two of which use state constraints, with one using mixed constraints and the other using differential inclusion. For the last one we formulated a family of more regular optimal control problems that approximate from below the optimal value. We compare the numerical performance of each one over both, an academic example and the minimization of peak of infectious in a SIR model with a L 1 constraint on the control. We summarize advantages and drawbacks of the different formulations for numerical computations in Table 6.8 section 6.6.

Chapter 7 shows proof that the optimal control that minimizes the peak of I, with a L 1 constraints in Figure 1.2: BocopHJB interface and some options the control, has the structure null-singular-null showed in equation (6.5.2). We also compare this solution with the one obtained by Morris et al that impose just one interval of time where the intervention can happen.

Finally, in Chapter 8 we show conclusion and perspectives for each part of this manuscript.

Part I

Mining context and the Open Pit problem

Chapter 2

The Open Pit problem

Mining is among the main Chilean economic activities, according to the technical report Anuario de la mineria de Chile prepared by Sernageomin [7], the governmental agency of geology and mining in Chile. It represented 12.5% of the CDP in 2020, mostly thanks to copper exports, of which Chile is the world's leading producer. Because of this, there exists in Chile a very active community of mining research which motives the first part of this work. The Final Open Pit problem (FOP) consists of finding the optimal shape of an open pit mine (in Figure 2.1 an example), in order to maximize the extraction profit while complying with a slope constraint which ensures the mine does not collapse. From this base problem, many other more complex ones can be set, for example, if we add a capacity in the number of blocks, the problem is called Capacitated Final Open Pit (CFOP). Another well studied problem considers the time periods for which we would like to obtain the shape of the mine, thus, obtaining an excavation plan. This problem is called Dynamical Final Open Pit (DFOP) which we also call Sequential Open Pit (SOP). We will present in the next section the classical formulation of the FOP, an integer programming prob-9 lem based on a block model introduced in 1968 by Johnson in [START_REF] Johnson | Optimum open pit mine production scheduling[END_REF], and then the continuous framework introduced in 2011 by Alvarez et al. [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF], which is the seminal paper of our work.

Integer programming formulation

Let B be a set of blocks where for every i ∈ B we know its benefit of extraction called b i . To extract a block i it is necessary to take other blocks in order to respect the slope constraint. This information is stored as an arc in a direct graph (B, A), i.e., (i, j) ∈ A if and only if to extract block i you have to take block j, in which case j is a predecessor of i. This is shown in Figure 2.2. To set the optimization problem we define the following variables for each block i ∈ B :

x i = 1 if block i is chosen for extraction 0 if not .
therefore, the binary optimization problem is:

(FOP 0 ) ∑ i∈B b i x i x j -x i ≥ 0 ∀(i, j) ∈ A x i ∈ {0, 1} ∀i ∈ B .
adding a maximal capacity C of extraction in terms of total mass, the model can be modified as

(CFOP 0 ) max ∑ i∈B b i x i x j -x i ≥ 0 ∀(i, j) ∈ A ∑ i∈N p i x i ≤ C x i ∈ {0, 1} ∀i ∈ B
where p i is the density of block i. To pose the dynamical problem we consider T time periods. For each time t ∈ {1, .., T } and each block i ∈ B we define the variable:

x t i = 1 if block i is extracted at time t 0 if not .
So, for a discount factor α ∈ (0, 1) the Dynamical Final Open Pit, named also by us as Sequential Open Pit (SOP), is written as follow

(SOP 0 ) max T ∑ t=1 ∑ i∈B b i (1 + α) t-1 x t i T ∑ t=1 x t i ≤ 1 ∀i ∈ B t ∑ l=1 x l j -x t i ≥ 0 ∀(i, j) ∈ A, t = 1, ..., T ∑ i∈B p i x t i ≤ C t t = 1, ..., T x t i ∈ {0, 1} ∀i ∈ B,t = 1, ..., T .
These kinds of formulations have been studied for several years and many others constraints and conditions can be added (see for example [START_REF] Espejo | A methodology for automatic ramp design in open pit mines[END_REF][START_REF] Jélvez | A new model for automated pushback selection[END_REF][START_REF] Saavedra-Rosas | Optimizing open-pit block scheduling with exposed ore reserve[END_REF]). A complete review of theses models can be found in [START_REF] Newman | A review of operations research in mine planning[END_REF]. In the following section, we will present a newer and less exploited model, which is based of the first part of this thesis.

Continuous framework

The formulation presented in this section comes from [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF] and we start by introducing the notation used.

Let Ω be the region of interest in R 1 or R 2 , supposed to be open and bounded. The border of a pit shall be determined by a continuous function p : Ω → R that is called profile and where p(x) represents the depth of the pit at point x ∈ Ω.

The slope of the pit at point x ∈ Ω will be computed as the Lipschitz modulus L p (x) defined as follows

L p (x) := lim sup x→x← x |p( x) -p( x)| || x -x|| (2.1)
and it will be bounded by a function w writing L p (x) ≤ w(x, p(x)), ∀x ∈ Ω. The natural shape of the ground will be represented by a continuous function p 0 , therefore, each feasible profile has to be deeper than this initial one and it has to connect with p 0 in the borders. That condition correspond to the following constraints:

p(x) -p 0 (x) ≥ 0, ∀x ∈ Ω p(x) -p 0 (x) = 0, ∀x ∈ ∂ Ω.
In Figure 2.3 it is possible to observe the previous description. To finally set the optimization problem, the effort and gain density functions are defined as e(x, z) ≥ e o > 0, g(x, z) ∈ R, ∀(x, z) ∈ Ω×Z, supposed uniformly bounded, and the functionals

G[p, q] := Ω q(x) p(x) g(x, z)dzdx, E[p, q] := Ω q(x) p(x) e(x, z)dzdx G[q] := G[p 0 , q] E[q] := E[p 0 , q] Figure 2.3: Explication image from [3].
The Final Open Pit problem in this continuous framework is then posed as:

(FOP c ) max G[p] p(x) -p 0 (x) ≥ 0, ∀x ∈ Ω p(x) -p 0 (x) = 0, ∀x ∈ ∂ Ω L p (x) ≤ w(x, p(x)) ∀x ∈ Ω p ∈ C (Ω)
and after adding the capacity constraint, the Capacitated Final Ope Pit correspond to:

(CFOP c ) max G[p] p(x) -p 0 (x) ≥ 0, ∀x ∈ Ω p(x) -p 0 (x) = 0, ∀x ∈ ∂ Ω L p (x) ≤ w(x, p(x)) ∀x ∈ Ω E[p] ≤ C p ∈ C (Ω)
.

To model the dynamical problem in [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF] they consider a continuous interval of time [0, T ] and for each time

t ∈ [0, T ], p(t, •) is a feasible profile of (FOP c ). Let also c ∈ L ∞ (0, T ), c(t) ≥ 0 be the capacity at time t and C(s,t) = t s c(τ)dτ the total capacity in the interval [s, t] ⊆ [0, T ].
Introducing a monotonically decreasing function ϕ ∈ C 1 [0, T ] representing the discount factor, such that ϕ(0) = 1 and 0 < ϕ(T ) < 1 (typically ϕ(t) = e -δt , for a given δ > 0), and noting P(t)(x) := p(t, x), the present value of the gain is then:

T 0 ϕ(t) Ω g(x, p(t, x))dxdP(t) = Ω T 0 ϕ(t)g(x, p(t, x)) ∂ p ∂t (t, x)dtdx.
Finally, the dynamical problem in the continuous framework can be written as:

(SOP c ) max Ω T 0 ϕ(t)g(x, p(t, x)) ∂ p ∂t (t, x)dtdx p(t, x) = p 0 (x) ∀x ∈ ∂ Ω, t ∈ [0, T ] L p(t,•) (x) ≤ w(x, p(t, x)) ∀x ∈ Ω, t ∈ [0, T ] p 0 = p(0, •) ≤ p(s, •) ≤ p(t, •) ∀s,t ∈ [0, t], s ≤ t Ω p(t,x) p(s,x) e(x, z)dzdx ≤ C(s,t) ∀s,t ∈ [0, t], s ≤ t p(t, •) ∈ C 1 (Ω)
Remark 2.2.1. The objective function of (SOP c ) suggest that p(t, x) must be differentiable with respect to t, so, to avoid that requirement, it can be written in the following form (expression obtained integrating by parts):

ϕ(T ) Ω p(T,x) p(0,x) g(x, z)dzdx + Ω T 0 -ϕ ′ (t) p(t,x) p(0,x) g(x, z)dzdtdx
In this work the authors proved existence results for previous problems (see propositions 6 and 11 in [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF]) but neither was any characterization of solutions nor numerical experiments shown. These two subjects are the novelty of Part I of this thesis.

Chapter 3

Analysis of optimality conditions for 2D and 3D Final Open Pit

This chapter corresponds to the article [START_REF] Amaya | Optimality conditions for the continuous model of the final open pit problem[END_REF] entitled "Optimality conditions for the continuous model of the final open pit problem"

Introduction

The long term planning of a mine operation consists of defining a sequence for the extraction of material from the mine in order to maximize profit. As a first step in this process, decision-makers usually must decide the final pit limit, which corresponds to the identification of a maximum value on the total mass to be extracted from the site, which enables an upper bound on the discounted value of the profit over several periods to be defined. This first step is called the Final Open Pit or Ultimate Open Pit problem. A very early contribution to the practical resolution of this problem was proposed by Lerchs and Grossman [START_REF] Lerchs | Optimum Design of Open-Pit Mines[END_REF] and, since then, a great variety of models and algorithms have been proposed. See Hustrulid et al [START_REF] Hustrulid | Open pit mine planning and design, two volume set & CD-ROM pack[END_REF] and Newman et al [START_REF] Newman | A review of operations research in mine planning[END_REF] for a more thorough introduction to open pit mine planning. The first effort to formally describe a practical mathematical model to solve this problem in an integrated way seems to be the work by Johnson [START_REF] Johnson | Optimum open pit mine production scheduling[END_REF].

Three different problems are usually considered for the economic valuation, design and planning of open pit mines. The first is the Final Open Pit (FOP) problem, which aims at finding the region of maximal economic value under geotechnical stability constraints. Another more realistic problem is what we call here the Capacity Final Open Pit (CFOP), which adds an additional constraint on the total capacity for extraction. The third problem is a multi-period version of the latter, which we call the Capacity Dynamic Open Pit (CDOP) problem, with the goal of finding an optimal sequence of volumes to be extracted with bounded capacities during each period.

The usual formulation of these problems consists of describing an ore reserve as a three-dimensional block model. Each block corresponds to a unitary volume of extraction, characterized by several physical and economic attributes, most of which are estimated from experimental sampling. Block models can be represented as directed graphs where nodes represent the blocks and arcs determineblock precedence (order of extraction). Block precedence is essentially induced by operational constraints, such as those derived from slope stability. This discrete approach usually gives rise to huge, combinatorially large-scale instances of Integer Programming, such as that presented by Cacetta [START_REF] Caccetta | Application of optimisation techniques in open pit mining[END_REF]. A great number of publications dealing with discrete block modeling for open pit mines have been published over the last 60 years. The seminal methodology for obtaining the ultimate pit limit, introduced by Lerchs and Grossman [START_REF] Lerchs | Optimum Design of Open-Pit Mines[END_REF], has been extensively applied in real mines for many years. The capacity dynamic problem is more difficult to solve and many methods using discrete optimization techniques have been proposed by Boland et al [START_REF] Boland | A strengthened formulation for the open pit mine production scheduled problem[END_REF], Cacetta and Hill [START_REF] Caccetta | An application of branch and cut to open pit mine scheduling[END_REF] and Hochbaum and Chen [START_REF] Hochbaum | Performance analysis and best implementations of old and new algorithms for the open-pit mining problem[END_REF]. This problem is beyond the scope of this chapter, but we can mention some dynamic programming formulations, for instance, Johnson and Sharp [START_REF] Johnson | A Three-dimensional dynamic programming method for optimal ultimate open pit design[END_REF] and Wright [START_REF] Wright | The use of dynamic programming for open pit mine design: some practical implications[END_REF]. Metaheuristic and evolutionary algorithms have also been extensively tested by Denby and Schofield [START_REF] Denby | Open-pit design and scheduling by use of genetic algorithms[END_REF] and Ferland et al [START_REF] Ferland | Application of a particle swarm algorithm to the capacitated open pit mining problem[END_REF].

In this chapter we use an alternative approach to the above mentioned (CFOP) problem based on a continuous framework, proposed by Alvarez et al [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF]. The basic idea is to describe pit contours by a continuous real-valued function, which maps each pair of horizontal coordinates to the corresponding vertical depth. Slope stability is ensured by means of a spatially distributed constraint on the local Lipschitz constant of the profile function. The maximal feasible local slope may vary throughout the site, depending on the geotechnical properties of the mineral deposit. The extraction capacity and operational costs are described by a possibly discontinuous effort density, a scalar function defined on the three-dimensional mining site. Concerning the continuous approach, we mention here the contribution by Ekeland and Queyranne [START_REF] Ekeland | Optimal pits and optimal transportation[END_REF], who proposed an alternative approach based on determining an optimum pit from an optimum dual solution of a particular transportation problem. Additionally, in [START_REF] Griewank | Duality results for stationary problems of open pit mine planning in a continuous function framework[END_REF], the authors derive duality results for the stationary open pit problem in the continuous framework, employing an additional condition called convex-likeness. The same authors, in [START_REF] Strogies | A PDE constraint formulation of open pit mine planning problems[END_REF], propose a partial differential equation model and show that, under suitable assumptions, the physically stable excavation path is the solution of a certain Hamilton-Jacobi equation.

The economic value of the blocks is given by a gain density defined on the deposit, which can also be a discontinuous function. Our goal here is to extend the existence results develop by Alvarez et al [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF] to the qualitative properties of the optimal solutions. This qualitative characterization is derived from the optimality conditions in the calculus of variations and control theory.

The chapter is organized as follows. In Section 3.2 we describe the stationary problem in terms of continuous profile functions and we establish the basis of our approach, in the context of a "2D-mine", which permits to give a simple motivation of the real 3D problem and to derive relevant results that can be generalized to the real case. Section 3.3 is devoted to the study of the realistic 3D instance, extending the main results of the previous section. By using tools from the calculus of variations, we derive an operational characterization of the optimal profile, particulary to show that the gain function must take the value zero along the border of the optimal profile, unless the capacity or slope constraints are active. In Section 3.4 we briefly summarize the main contributions of this chapter and indicate some avenues for future research.

The 2D open pit problem

To fix ideas, we begin by considering the idealized case of an open pit on the plane, that is, the framework where the profiles are modeled using a continuous function that depends only on a single space variable (denoted x for simplicity).

Statement of the problem

For the sake of notation, we assume that the depth of a profile is always positive. In this framework, an admissible profile is a function p : [a, b] → R that must satisfy some conditions, the first one being as follows: given an initial profile p 0 : [a, b] → R an admissible profile has to satisfy

p 0 (x) ⩽ p(x), ∀x ∈ [a, b]. Ground level x a b x p(x) F grav
L p (x) := lim sup x→x← x |p( x) -p( x)| | x -x| .
Due to the risk of landslides, the slope of a profile cannot be too steep. Note that the maximal slope allowed may change depending on the position and depth in the pit. This constraint is then represented via the condition

L p (x) ⩽ κ(x, p(x)), ∀x ∈ [a, b],
where κ(x, z) represents the maximal slope at the point (x, z) allowed for a profile p : [a, b] → R to be admissible. Note that if the profile is continuously differentiable on (a, b), then the slope agrees with the absolute value of the profile's derivative (see [42, Section 1D] ), that is,

L p (x) = | ṗ(x)|, ∀x ∈ (a, b).
However, in our setting, working with smooth functions is too restrictive. For this reason we choose to work with a broader class of functions, namely, the collection of continuous functions whose derivatives exist almost everywhere on [a, b] and which satisfy

p(x) = p(a) + x a ṗ(s)ds, ∀x ∈ [a, b].
This class of functions is the so-called set of absolutely continuous functions, which we denote by A C [a, b]. It turns out that absolutely continuous functions are well behaved with respect to the slope, in the sense that the slope agrees almost everywhere with the derivative of an absolutely continuous profile.

Lemma 3.2.1. Let p ∈ A C [a, b], then L p (x) = | ṗ(x)| almost everywhere on [a, b].
On the other hand, due to physical or economic constraints, the capacity of extraction is indeed limited. Given a position x, the effort associated with extracting a block at depth z ⩾ p 0 (x) can be represented by a nonnegative quantity e(x, z). Thus, given a maximal budget c max > 0, the capacity constraints associated with a profile p : We are now in a position to formally state the 2D Final Open Pit problem:

                     Maximize b a p(x) p 0 (x) g(x, z)dzdx over all p ∈ A C [a, b] subject to p(a) = p 0 (a), p(b) = p 0 (b) p 0 (x) ⩽ p(x), for all x ∈ [a, b], | ṗ(x)| ⩽ κ(x, p(x)) for a.e. x ∈ [a, b] b a p(x) p 0 (x)
e(x, z)dzdx ⩽ c max .

(P 2D )

Standing Assumptions

Throughout the remainder of this section, unless otherwise stated, we will assume -∞ < a < b < +∞ and c max > 0 are fixed parameters of the problem. The initial profile p 0 : [a, b] → R is a given continuously differentiable function.

The profit objective function g : [a, b] × R is assumed to be a nonnegative, bounded and piecewise continuous function. The marginal cost of extraction e : [a, b] × R → R is assumed to be a nonnegative, bounded and continuous function. Also, the maximal slope allowed κ : [a, b] × R → R is assumed to be continuous, nonnegative and bounded with p → κ(x, p) being continuously differentiable for any x ∈ [a, b] fixed and such that (x, q) → ∇ q κ(x, q) is bounded on [a, b] × R n .

Under these assumptions, the existence of an optimal profile is ensured, as proved by Alvarez et al [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF]. Moreover, the fact that an optimal profile is absolutely continuous (Lipschitz continuous actually) is enforced by the boundedness and continuity of the maximal slope κ. This existence result concerns as well the 3D case studied in §3.3. Remark 3.2.1. In the light of Alvarez et al [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF]Lemma 1], the feasible set of (P 2D ) without the capacity constraint, is convex provided z → κ(x, z) is concave for any x ∈ [a, b] fixed. Moreover, by [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF]Proposition 5], if z → e(x, z) is monotonically increasing and z → g(x, z) is monotonically decreasing (for x ∈ [a, b] fixed), then the problem (P 2D ) turns out to be a convex one. The previous paragraph can also be applied to the 3D case. However, under the assumptions we have done so far, these hypotheses cannot be assured. As a matter of fact, the problem (P 2D ) may have several local minima, which are not necessarily global. It worths to mention then that in general setting of this manuscript we deal with non-convex problems. The previous comment also applies to the 3D case.

Basics on state constrained optimal control

Let us point out that the formulation of (P 2D ) is slightly more restrictive than what has been treated in Alvarez et al [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF]. Essentially, we restrict our analysis to a small class of functions, those that are absolutely continuous. The main advantage of doing so is that now the Final Open Pit problem can be treated as an optimal control problem with state constraints, and optimality conditions can be derived by fairly standard methods.

For the sake of completeness, we state the main tool from optimal control theory we are going to use in the analysis provided in this section. Let us consider a general Mayer optimal control problem on R n :

                   Minimize ϕ(q(b)) over all q ∈ A C n [a, b
] and measurable functions u satisfying q(x) = f (x, q(x), u(x)) , for a.e. x ∈ [a, b], u(x) ∈ U, for a.e. x ∈ [a, b], h(x, q(x)) ⩽ 0, for any x ∈ [a, b], (q(a), q(b)) ∈ E.

(P M ) Here, q ∈ A C n [a, b] means that q : [a, b] → R n and if q = (q 1 , . . . , q n ), then each component q i ∈ A C [a, b].
Furthermore, for the purposes of our analysis, we only need to consider the case in which:

• ϕ : R n → R is a continuously differentiable function, • f : [a, b] × R n × R m → R n is such that, x → f (x, q, û) is measurable, q → f ( x, q, û) is Lipschitz continuous (uniformly with respect to ( x, û)) and u → f ( x, q, u) is continuous for any ( x, q, û) ∈ [a, b] × R n ×U fixed, • h : [a, b] × R n → R is continuous, with q → h( x, q) being differentiable for any x ∈ [a, b] fixed and such that (x, q) → ∇ q h(x, q) is continuous on [a, b] × R n , • U ⊆ R m is a given nonempty compact set, • E ⊆ R n × R n is a nonempty closed convex set.
It is worth recalling that the (convex) normal cone to a set S ⊆ R k is defined by

N S (s) := {η ∈ R k | ⟨η, s -s⟩ ⩽ 0, ∀ s ∈ S}, ∀s ∈ S.
In particular, given s 0 ∈ R, we have

N (-∞,s 0 ] (s) = {0} if s < s 0 [0, +∞) if s = s 0 ,
and N {s 0 } (s 0 ) = R.

Definition 3.2.1. An arc q ∈ A C n [a, b] admissible for (P M ) is said to be a weak local minimizer of the problem (related to an optimal control ū) if there is ε > 0 such that

q ∈ A C n [a, b] is admissible for (P M ) and ∥q -q∥ W 1,1 ⩽ ε =⇒ ϕ ( q(b)) ⩽ ϕ (q(b)) .
Here ∥ • ∥ W 1,1 stands for the usual norm of the Sobolev space W

: [a, b] → R n satisfying γ(x) = ∇ q h(x, q(x)), for µ-a.e. x ∈ [a, b], such that 1. (λ , µ, η) ̸ = (0, 0, 0); 2. -λ (x) ∈ ∂ C q H (x, q(x), ξ (x), ū(x)) for a.e. x ∈ [a, b]; 3. (λ (a), -ξ (b)) ∈ {0} × {η∇ϕ( q(b))} + N E ( q (a), q(b));
4. H (x, q(x), ξ (x), ū(x)) = max u∈U H (x, q(x), ξ (x), u) for a.e. x ∈ [a, b];

5. supp(µ) ⊆ {x ∈ [a, b] | h(x, q(x)) = 0}.
Here H (x, q, ξ , u) = ⟨ξ , f (x, q, u)⟩ and

ξ (x) = λ (x) + [a,x[ γ(s)µ(ds) ∀x ∈ [a, b[ and ξ (b) = λ (b) + [a,b]
γ(s)µ(ds)

Optimality conditions for the 2D open pit problem

In this part of the chapter, we analyze the behavior of an optimal profile by using the tools from optimal control theory described earlier. The following result can in principle be stated for local optima as well. However, to keep the presentation of the chapter simple, we prefer to present it only for a global optimum. Proof. The proof of the result is based on a transformation of the Final Open pit problem (P 2D ) into a Mayer problem such as (P M ), for which p provides a weak local minimizer related to the optimal control ū(x) :=

   ṗ(x) κ(x, p(x)) if κ(x, p(x)) ̸ = 0 0 otherwise. (3.1)
We divide the proof into several parts for the sake of exposition.

1. First we show that (P 2D ) is an instance of the Mayer problem (P M ). The key points here are to interpret the slope condition as a controlled ordinary differential equation and to be able to handle the capacity constraints On the other hand, note that (3.2) is actually an isoperimetric inequality constraint. To deal with it, we introduce a new auxiliary state. Let q 1 : [a, b] → R be given by

q 1 (t) = t a p(x) p 0 (x) e(x, z)dzdx, ∀t ∈ [a, b].
Thus, it is clear that (3.2) can be written as q 1 (b) ⩽ c max . Furthermore, q 1 (a) = 0 and the velocity of q 1 is given by the expression q1 (x) = Also, by defining q 2 : [a, b] → R via the formula

q 2 (t) = t a p(x) p 0 (x) g(x, z)dzdx, ∀t ∈ [a, b],
it is clear that the total profit is given by q 2 (b), and that this new state satisfies q2 (x) = p(x) p 0 (x) g(x, z)dz, for a.e. x ∈ [a, b], with q 2 (a) = 0.

Therefore, setting q(x) = (q 1 (x), q 2 (x), p(x)) for any x ∈ [a, b], we see that (P 2D ) is an instance of the Mayer problem (P M ) with ϕ(q) = -q 2 , h(x, q) = p 0 (x)q 3 , U = [-1, 1],

f (x, q, u) = q 3 p 0 (x)
e(x, z)dz,

q 3 p 0 (x) g(x, z)dz, uκ(x, q 3 ) and E = (α, β ) ∈ R 3 × R 3 | α 1 = α 2 = 0, α 3 = p 0 (a), β 1 ⩽ c max and β 3 = p 0 (b) .
2. Now, since p is assumed to be an optimal solution of (P 2D ), it follows that p provides a weak local minimizer of (P M ), related to the optimal control defined in (3.1).

Moreover, the condition under which Lemma 3.2.2 has been stated are satisfied by the data provided in the preceding part; the Lipschitz continuity of q → f ( x, q, û) (uniformly with respect to ( x, û) ∈ [a, b] × U) comes from the fact that e and g are measurable bounded functions and κ is Lipschitz continuous in the second variable, uniformly with respect to the first one. Therefore, we can apply Lemma 3.2.2, and so, there exist λ

∈ A C 3 [a, b], η ∈ {0, 1}, a (positive) Radon measure µ on [a, b],
and a Borel measurable function γ : [a, b] → R 3 fulfilling the conditions in Lemma 3.2.2. Note first that the Hamiltonian does not depend on q 1 nor on q 2 , and also that ∇ q h(x, q) = (0, 0, -1).

Because of point 2 in Lemma 3.2.2 and the definition of x → ξ (x), we can deduce that there are λ1 , λ2 ∈ R such that

ξ 1 (x) = λ 1 (x) = λ1 and ξ 2 (x) = λ 2 (x) = λ2 , ∀x ∈ [a, b],
and

ξ 3 (x) = λ 3 (x) -µ([a, x[), ∀x ∈ [a, b[, and ξ 3 (b) = λ 3 (b) -µ([a, b]).
Note also that ∇ϕ(q) = (0, -1, 0) and Note that whenever κ(x, p(x)) ̸ = 0 ( a.e. on [a, b]) we have that

N E ( q(a), q(b)) = R 3 × β ∈ R 3 | β 1 ⩾ 0, β 1 ( q1 (b) -c max ) =
ξ 3 (x) ṗ(x) = |ξ 3 (x)|κ(x, p(x)).
This implies that whenever ṗ(x) < κ(x, p(x)), then necessarily ξ 3 (x) = 0, and so

ξ 3 (x) ( ṗ(x) -κ(x, p(x))) , for a.e. x ∈ [a, b].
Finally, note that since e is continuous and κ is continuously differentiable in the second variable, for any x ∈ [a, b] ξ ∈ R 3 and u ∈ [-1, 1] fixed, such that p → q(x, p) is continuous at p = q 3 , we have ∇ q H(x, q, ξ , u) = (0, 0, ξ 1 e(x, q 3 ) + ξ 2 g(x, q 3 ) + ξ 3 u∂ p κ(x, q 3 ))

In particular, by point 2 in Lemma 3.2.2, for a.e. x ∈ [a, b] such that p → q(x, p) is continuous at p = p(x) we have

-λ3 (x) = λ1 e(x, p(x)) + ηg(x, p(x)) + ξ 3 (x) ū(x)∂ p κ(x, p(x))
because in this case q → H(x, q, ξ (x), ū(x)) is continuously differentiable at

q = x a p( x) p 0 ( x) e( x, z)dzd x, x a p( x) p 0 ( x) g( x, z)dzd x, p(x) .
Since the functions g is piecewise continuous, for any x ∈ [a, b], if p → q(x, p) is not continuous at p = p(x) we have that

∂ C q H(x, q, ξ , u) = {(0, 0)} × (ξ 1 e(x, q 3 ) + ξ 2 G(x, q 3 ) + ξ 3 u∂ p κ(x, q 3 )) where G(x, p) = co{g(x, p -), g(x, p + )}, ∀x ∈ [a, b], p ∈ R.
Also, on the one hand, by Maximum Principle (point 4 in Lemma 3.2.2), for a.e. x ∈ [a, b] such that κ(x, p(x)) > 0 we must have that ξ 3 (x) ū(x) = |ξ 3 (x)|. On the other hand, if κ(x, p(x)) = 0, we must have that ∂ p κ(x, p(x)) = 0 because p = p(x) is a local minimum of p → κ(x, p). Combining these two issues we get

ξ 3 (x) ū(x)∂ p κ(x, p(x)) = |ξ 3 (x)|∂ p κ(x, p(x)), for a.e. x ∈ [a, b].
Therefore, setting λ = λ1 and ζ = λ 3 the conclusion follows.

We now state a direct consequence of the preceding theorem in the case when the slope condition is not active, and the state constraint is only active at the end-points.

The 3D open pit problem

We now turn into the more realistic case of an open pit in the 3D space. The profiles in this framework are modeled using a continuous function that depends on the two horizontal space variable (denoted x and y for simplicity). Generically, we denote a profile of an open pit by p : Ω → R where Ω ⊆ R 2 is the bounded domain in R 2 that represents the open pit and where p(x, y) represents the depth of the profile at the point (x, y) ∈ Ω.

Statement of the problem

As done for the 2D case, we assume that the depth of a profile is always positive. The final open pit problem in the 3D case has the same structure as in the 2D case. This means that for a given initial profile p 0 : Ω → R, the total profit and total extraction associated with an admissible profile p : Ω → R are given respectively by e(x, y, z)dzdxdy.

The maximal slope allowed is also considered to be bounded, and thus profiles are Lipschitz continuous mappings. The associated constraint is then represented via the condition

L p (x, y) := lim sup ( x, ȳ)→(x,y)←( x, ŷ) |p( x, ȳ) -p( x, ŷ)| | x -x| 2 + | ȳ -ŷ| 2 ⩽ κ(x, y, p(x, y)), ∀(x, y) ∈ Ω.
Therefore, the Final Open Pit problem in the 3D case is the following: e(x, y, z)dzdxdy ⩽ c max .

                           Maximize Ω p(x,
(P 3D )

Remark 3.3.1. A more general model that considers profiles having a time dependance has been studied by Álvarez et al in [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF]. The analysis of this problem, called Capacitated Dynamic Open Pit, becomes more difficult and we plan to study it in details elsewhere.

Remark 3.3.2. Similarly as for the 2D case, a control setting can be introduced to deal with the 3D case; see the proof of Theorem 3.2.1. This is certainly a suitable approach to handle Theorem 3.3.1, however, in this setting the control is distributed and also subject to constraints; see the discussion in §3.4. Moreover, optimality conditions for problems of this kind are known to be harder to handle and for this reason we take another path to prove Theorem 3.3.1 base on classical calculus of variations.

Standing Assumptions

Throughout the remainder of this section, unless otherwise stated, we will assume Ω ⊆ R 2 is an open bounded domain and c max > 0 is fixed parameter of the problem. The initial profile p 0 : Ω → R is a given continuously differentiable function.

The densities of gain and effort are now g : Ω × R and e : Ω × R → R, respectively. They are assumed to be bounded, measurable and the second one (the densities of effort) nonnegative. Also, the maximal slope allowed κ : Ω × R → R is assumed to be continuous, nonnegative and bounded with p → κ(x, p) being continuously differentiable for any x ∈ Ω fixed and such that (x, q) → ∂ p κ(x, p) is bounded on Ω × R n .

Optimality conditions

We now present some necessary optimality conditions that extend the one given for the 2D case. The conditions obtained in this case do not require the continuity of the gain function g nor the continuity of the effort e. However, because of the nonholonomic character of the slope constraints, the result we present is only valid for the case when optimal profiles do not saturate this condition (see assumption (3.4) below). Nonholonomic constraints are hard to handle in calculus of variations of multiple integral and require technical assumptions which may be too strong for the scope of this chapter. The main difficulty is that the construction of suitable variations is not always ensured; see for instance [START_REF] Giaquinta | Calculus of variations I[END_REF]Chapter 2]. It remains then as an open problem and future work to provide necessary optimality conditions for the general case where nonholonomic restriction may be active. Theorem 3.3.1. Let p ∈ Lip Ω be an optimal profile of (P 3D ). Assume that p ̸ = p 0 and let Ω 0 ⊆ Ω be the open domain of R 2 given by

Ω 0 = (x, y) ∈ R 2 | p 0 (x, y) < p(x, y) .
If the slope constraints is not active on Ω 0 , that is,

sup (x,y)∈Ω 0 L p(x, y) -κ(x, y, p(x, y)) < 0, (3.4)
then there is λ ⩽ 0 such that g(x, y, p(x, y)) + λ e(x, y, p(x, y)) = 0, ct p. in Ω 0 .

Furthermore, λ satisfies the following properties:

1. If Ω p(x,y) p 0 (x,y)
e(x, y, z)dzdxdy < c max then λ = 0.

λ can be taken to be any value

λ = - Ω g(x, y, p(x, y))ψ(x, y)dxdy, provided that ψ ∈ C ∞ 0 (Ω) is such that Ω e(x, y, p(x, y))ψ(x, y)dxdy = 1.
Proof. The proof follows rather standard arguments in calculus of variations, adapted to be able to handle the integral inequality constraint of isoperimetric type. Assume first that e(x, y, p(x, y)) is not identically zero in Ω 0 . Take some ψ ∈ C ∞ 0 (Ω 0 ) such that Ω 0 e(x, y, p(x, y))ψ(x, y)dxdy = 1.

Since e(x, y, p(x, y)) ⩾ 0 in (x, y) ∈ Ω 0 and it is not identically zero, the existence of such function ψ is guaranteed. Now take ϕ ∈ C ∞ 0 (Ω 0 ) arbitrary and define for s,t ∈ R the profile p s,t ∈ Lip Ω given by p s,t (x, y) = p(x, y) + sϕ(x, y) +tψ(x, y). Consider the functions (s,t) → f (s,t) and (s,t) → h(s,t) defined on R 2 via the formulas

f (s,t) := Ω 0 p s,t (x,y) p 0 (x,y) g(x, y, z)dzdxdy and h(s,t) := Ω 0 p s,t (x,y) p 0 (x,y) e(x, y, z)dzdxdy.
By the definition of Ω 0 , the continuity of p → κ(x, y, p) and (3.4), it follows then that there is δ > 0 such that for any (x, y) ∈ Ω 0 and any s,t ∈ (-δ , δ ) we have p s,t (x, y) > p 0 (x, y) and L p s,t (x, y) < κ(x, y, p s,t (x, y)).

Since p is an optimal profile for the final open pit problem (P 3D ), it follows that (0, 0) is a local maximum of the problem Maximize f (s,t) over all s,t ∈ R subject to h(s,t) ⩽ c max .

This nonlinear optimization problem satisfies the so-called Mangasarian-Fromovitz condition because

∂ t h(0, 0) = Ω 0 e(x, y, p(x, y))ψ(x, y)dxdy = 1.
Where the first equality is justified by the Dominated Convergence Theorem and the fact that e ∈ L ∞ (Ω × R). Therefore, by the Karush-Kuhn-Tucker theorem, there is λ ⩽ 0, which in principle depends on ψ and ϕ, such that ∇ f (0, 0) + λ ∇h(0, 0) = 0 and λ (h(0, 0)c max ) = 0.

Moreover, similarly as justified above, it is not difficult to see that

∂ s h(0, 0) = Ω 0 e(x, y, p(x, y))ϕ(x, y)dxdy, ∂ s f (0, 0) = Ω 0 g(x, y, p(x, y))ϕ(x, y)dxdy, ∂ t f (0, 0) = Ω 0 g(x, y, p(x, y))ψ(x, y)dxdy.
On the one hand, by the condition over the partial derivatives with respect to the t variable we get λ does not depend on ϕ because

Ω 0 g(x, y, p(x, y))ψ(x, y)dxdy = ∂ t f (0, 0) = -λ ∂ t h(0, 0) = -λ .
On the other hand, by the condition over the partial derivatives with respect to the s variable we get Ω 0 g(x, y, p(x, y)) + λ e(x, y, p(x, y)) ϕ(x, y)dxdy = 0.

But, since ϕ ∈ C ∞ 0 (Ω) is arbitrary and (x, y) → g(x, y, p(x, y)) + λ e(x, y, p(x, y)) belongs in particular to L 2 (Ω 0 ), by the fundamental lemma of the calculus of variations (cf. [66, Lemma 3.2.3]) the conclusion follows.

Finally, for the case that e(x, y, p(x, y)) is identically zero in Ω 0 it is enough to define s → f (s) on R via the formula

f (s) := Ω 0 p(x,y)+sϕ(x,y) p 0 (x,y) g(x, y, z)dzdxdy, with ϕ ∈ C ∞ 0 (Ω 0 )
arbitrary and check that s = 0 is a local minimum of f . Then the conclusion follows by using the Fermat rule ( ḟ (0) = 0) and fundamental lemma of the calculus of variations.

Remark 3.3.3. Note in particular that Theorem 3.3.1 says that if the marginal cost associated with extracting a block at any depth is zero (there is no capacity constraint), that is, e(x, y, z) = 0 for any (x, y) ∈ Ω and z ⩾ p 0 (x, y), then the marginal gain of extracting a block at any depth must be zero on the subsection Ω 0 , that is, g(x, y, p(x, y)) = 0, ct p. in Ω 0 .

Future work and final remarks

In this chapter we have provided necessary optimality conditions for a profile of an open pit mine to be an optimum for the final open pit problem in the 2D as well as on the 3D. Both settings involve isoperimetric restriction, which are hard to handle in general. Nonholonomic restrictions, such as the slope condition has been treated only for the 2D case. The 3D case remains as an open question that deserves some attention, and which we plan to address in a future work. Obtaining numerical solutions is still an open problem. For this purpose, a classical direct method or an indirect method using the results of this chapter can be implemented. Preliminary simulations have been done with the help of the INRIA solver for optimal control problems BOCOP [START_REF] Bonnans | Bocop -A collection of examples[END_REF]. This is an issue that need to be investigated in more details.

Finally, let us mention that, similarly as done for the 2D case, the maximal slope condition in the 3D can actually be subsumed by a control type condition

∇p(x, y) = κ(x, y, p(x, y)) (cos(θ ), sin(θ )) , for a.e. (x, y) ∈ Ω, θ ∈ [0, 2π).
Thus, a possible way to address the final open pit problem is to study the optimal control problem associated with this constraints. This issue needs to be investigated in details.

Chapter 4

Numerical study of 2D and 3D Final and Sequential Open Pit

This chapter corresponds to the submitted pre-print [START_REF] Molina | Optimal control approaches for Open Pit planning[END_REF] titled "Optimal control approaches for Open Pit planning"

Introduction

In long-term planning of mine operation, a common task consists in determining the profile of the total mass of material to be extracted from the site to optimally design an opencast mine. This so-called Final Open Pit problem was introduced in the early works Ref. [START_REF] Lerchs | Optimum Design of Open-Pit Mines[END_REF][START_REF] Johnson | Optimum open pit mine production scheduling[END_REF], with a more recent overview in Ref. [START_REF] Newman | A review of operations research in mine planning[END_REF].

The typical approach used to solve this problem is based on a discrete block model of the site, each block having an associated extraction cost and profit value, based on topographical and geological data. Using a graph of block precedence (i.e. order of extraction) allows to take into account slope constraints for the mine stability, and gives rise to large Integer Programming problems, see for instance Ref. [START_REF] Caccetta | Application of optimisation techniques in open pit mining[END_REF]. The dynamic programming approach was also investigated in this framework, see e.g. Ref. [START_REF] Wright | The use of dynamic programming for open pit mine design: some practical implications[END_REF]. Another approach presented in Ref. [START_REF] Strogies | A PDE constraint formulation of open pit mine planning problems[END_REF] uses a PDE formulation for time labeling functions.

The present chapter follows the continuous approach introduced in Ref. [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF] with the reformulation of the Open Pit using a calculus of variation framework, and then in Chapter 3 as an optimal control problem. The main contributions of the present work include the analysis of the Final Open Pit (FOP from now) with capacity, slope and initial profile constraints, using Pontryagin's Maximum Principle to extend the results previously obtained in chaper 3. Then we introduce a new semi-continuous formulation that can handle the Sequential Open Pit problem (SOP from now) ,i.e. optimization of the mine profile over a sequence of several time-frames, for a 2D space domain (3D mine profile). Finally, numerical simulations are provided for both the continuous and semi-continuous approaches, including global optimization for the 2D FOP case, and to our knowledge the first results for the 3D profile optimization as an optimal control problem. The outline of the chapter is as follows. After the introduction presenting context, Section II covers the SOP problem statement with the continuous approach, and introduces the semicontinuous formulation. Section III presents the FOP analysis using Pontryagin's Maximum Principle and in particular discusses the control structure in terms of bang, constrained and singular arcs. Section IV present the numerical simulations for three test cases: 2D FOP, 2D SOP and 3D SOP, and is followed by the conclusions.

Problem statement

For a given spatial domain Ω, we consider a continuous function p : Ω → R called profile that delimits the shape of the mine pit. The aim is to determine the profile that maximizes the gain from the excavated soil, while respecting some limits for the excavated capacity and maximal slope of the mine. We recall now the continuous approach for open pit planning, and introduce a new semi-continuous approach that can handle the 3D profile case. Both of these approaches lead to optimal control formulations of the problem.

Continuous formulation

The key idea in the so-called continuous formulation, originally introduced in [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF], is to use the distance (position along the x-axis) as independent variable, which allows to define the mine profile as a function of this new 'time'. Introducing a suitable dynamics for this function, with the associated control function, allows to formulate the open pit planning as an optimal control problem (OCP).

Final open pit planning problems (FOP)

For the 1-dimensional case, the domain Ω = [a, b] will correspond to the independent variable or 'time' of the optimal control problem. Consider the state variables P, c : [a, b] → R + for the depth profile of the pit and the excavated capacity of the mine. Let us also denote P 0 ∈ C 1 ≥ 0 the initial profile corresponding to the natural shape of the ground. We set the state constraint P(t) ≥ P 0 (t), ∀t ∈ [a, b], and the boundary conditions P(a) = P 0 (a), P(b) = P 0 (b). An additional final condition is that the total excavated capacity is limited, i.e. c(b) ≤ c max .

We also introduce κ : [a, b] × R → R * such that κ(t, z) is the maximal pit slope at position t and depth z. Instead of the original dynamics Ṗ = u, we choose to use a normalized control u : [a, b] → [-1, 1] which is a bit simpler than having the mixed state-control constraint |u(t)| ≤ κ(t, P(t)) for the maximal slope. As part of the soil characteristics, we also note g, e : [a, b] × R → R the densities of gain and effort for excavating at a given position and depth. The optimal control formulation of (FOP) is then as follows:

(FOP)                                        max b a P(t) P 0 (t) g(t, z)dzdt Ṗ(t) = u(t)κ(t, P(t)) ∀t ∈ [a, b] ċ(t) = P(t) P 0 (t) e(t, z)dz ∀t ∈ [a, b] u(t) ∈ [-1, 1] ∀t ∈ [a, b] P 0 (t) -P(t) ≤ 0 ∀t ∈ [a, b] P(a) = P 0 (a), P(b) = P 0 (b) c(a) = 0, c(b) ≤ c max Remark 4.2.1.
In the following we take the basic effort function e ≡ 1. The gain function g is typically defined by interpolation over tabular data, and has to be integrated numerically along the depth z.

Sequential open pit planning (SOP)

We introduce now an extended version of (FOP), in which we want to schedule an extraction program over N consecutive time-frames. This case is quite relevant in mine planning since mining companies divide the digging process into periods for operational purposes. We extend the notations of (FOP) to the multi-frame framework, and note P i the mine profile at time-frame i, with the associated control u i , while c i is the excavated capacity during time-frame i. Each mine profile has to be deeper than the previous one, i.e. the constraint P ≥ P 0 from (FOP) is generalized as P i ≥ P i-1 , i = 1 . . . N. The capacity limit c i max is now enforced at each individual time-frame. Finally, the objective function now takes into account a depreciation rate α > 0 over time, with the gains for the more distant time-frames being valued less than for the more immediate time-frames. This new optimal control problem reads as follows

(SOP)                                          max N ∑ i=1 b a P i (t) P i-1 (t) g(t, z) (1 + α) t-1 dzdt Ṗi (t) = u i κ(t, P i (t)) ∀t ∈ [a, b], i = 1, . . . , N ċi = P i (t) P i-1 (t) e(t, z)dz ∀t ∈ [a, b], i = 1, . . . , N u i (t) ∈ [-1, 1], ∀t ∈ [a, b], i = 1, . . . , N P i-1 (t) -P i (t) ⩽ 0 ∀t ∈ [a, b], i = 1, . . . , N P i (a) = P 0 (a), P i (b) = P 0 (b) i = 1, . . . , N c i (a) = 0, c i (b) ⩽ c i max i = 1, . . . , N Remark 4.2.2.
Note that (SOP) with N = 1 corresponds to (FOP). Numerically, the multi-process (SOP) can be reformulated by duplicating the state and control variables (as well as the constraints) for each time-frame. Adding the proper linking constraints between the final and initial conditions of the successive time-frames, we obtain a single process version of (SOP) that can be solved by standard methods. The overall problem dimension, however, is higher, therefore computationally expensive methods such as global optimization may be able to handle (FOP) but not (SOP), see section 4.4.

Remark 4.2.3. For the discrete (block) formulation, it is known (see for instance Ref. [START_REF] Caccetta | An application of branch and cut to open pit mine scheduling[END_REF] and [START_REF] Molina | Analytical Properties of the Feasible and Optimal Profiles in the Binary Programming Formulation of Open Pit[END_REF]) that each profile (or 'pit') which is solution of (SOP) is not deeper than the optimal pit of (FOP) with the same parameters and infinite capacity. A similar result has been obtained for the continuous framework in Ref. [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF].

Semi continuous formulation for SOP

The main limitation of the continuous approach is that using the independent variable to represent the position in space makes it difficult to handle the 3D profile case, both in terms of dynamics / controls and profile slopes. This is why we introduce a new approach called the semi-continuous formulation, based on an explicit discretization of the space domain Ω. The mine profile is therefore represented by a finite set of variables at the discretization nodes, as illustrated in Figure 4.1. Control variables are defined as the excavation effort at each discretization node. Slope constraints are modeled as state constraint linking each node with their neighbors. The independent variable is here standard time, expressed in time-frames such that the final time T is the total number of time-frames. Since SOP is a multi-phase problem, one standard way to formulate it is to normalize the time interval to [0, 1] and duplicate the variables for each time-frame. This approach yields another optimal control formulation of the Sequential Open Pit problem, for which extension from 1D to 2D space domain is rather straightforward, at the cost of an increase in overall problem dimension.

Notations. In the context of the semi-continuous approach, for functions of both space and time such as profile, controls and slopes, we will typically use subscripts for the space discretization node in Ω, and exponents for the time-frame of the multi-phase Sequential Open Pit. For instance, P k i will represent the profile depth at node i and time-frame k, and P k := (P k i ), i = 0, . . . , N refers to the mine profile at timeframe k. Similarly, U k i will denote the digging at node i and time-frame k, with U k := (P k i ), i = 0, . . . , N corresponding to the overall excavation effort over the domain Ω at time-frame k. We will also denote by

P k+1 P k
the integral of a function between the two mine profiles at time-frames k and k + 1; for 2D profiles this is a 2D integral along x and the depth z, and a 3D integral along x, y, z for 3D profiles. N , with N +1 discretization nodes (x i ), and note I = {0, . . . , N} the set of indices for the nodes. We define the state variables for the profile nodes (P i ) i∈I as functions of time. We also introduce the control variables at each node (U i ) i∈I ≥ 0, corresponding to the excavation effort, so that the profile variables follow the simple dynamics Ṗi

(t) = U i (t) , ∀i ∈ I , ∀t ∈ [0, T ]. (4.1)
Gain. The achieved gain during time-frame k is the integral of g between the current profile P k and the previous P k-1 . Taking into account the depreciation rate α introduced in 4.2.1, the overall gain to be maximized is

T ∑ k=1 P k P k-1 g(x, z) (1 + α) t k -1 dxdz. (4.2)
The computation of this objective is detailed in 4.6.

Slope. We denote S k i the slope at node i and time-frame k, which is a function of time. The maximal slope condition writes as

-1 ≤ S k i (t) κ x i , P k i (t) ≤ 1 , ∀i ∈ I , ∀k = 0 . . . T , ∀t ∈ [0, 1] (4.3)
In the 2D case we will use the simple slope formula

S k i = (P k i -P k i-1 )/∆x (4.4)
and the slope limits are state constraints.

Capacity. The excavation effort at each time-frame k corresponds to the integral of the effort E between the two consecutive profiles P k-1 and P k , and the capacity limit writes as

P k P k-1 e(x, z)dxdz ≤ C k , ∀k = 1, . . . , T. (4.5)
The computation of this integral is detailed in 4.6.

Initial profile. This is now a standard initial condition of the form

P 0 i (0) = p 0 (x i ) , ∀i ∈ I. (4.6)
We obtain the following multi-phase problem (SOP) 2D SC with Fig. 4.2 illustrating the profile discretization in the 2D case, with N = 7 and T = 2. Implementation details regarding the approximation of the various integrals are presented in 4.6

(SOP) 2D SC                                                  max T ∑ k=1 P k P k-1 g(x, z) (1 + α) k-1 dxdz Ṗk i (t) = U k i (t), i ∈ I, k = 1, . . . , T, t ∈ [0, 1] -1 ≤ S k i (t) κ(x i ,P k i (t) ≤ 1, i ∈ I, k = 1, . . . , T, t ∈ [0, 1] P k P k-1 e(x, z)dxdz ≤ C k , k = 1, . . . , T P 0 i (0) = p 0 (x i ), i ∈ I U k i (t) ≥ 0, i ∈ I, k = 1, . . . , T, t ∈ [0, 1]
Remark 4.2.4. Setting T = 1 corresponds to the Final Open Pit problem with a single time-frame.

Remark 4.2.5. The boundary condition P |∂ Ω = 0 is in practice built in directly in the problem formulation by eliminating the profile and control variables at the nodes corresponding to the boundary of the space domain.

Remark 4.2.6. Moreover, the constraint that each profile must be deeper than the previous one, which was a state constraint in the continuous formulation, is now simply enforced by the conditions U i ≥ 0.

Remark 4.2.7. The step size ∆x for the discretization of Ω in the semi-continuous approach can be seen as the analogue of the time step ∆t for the continuous approach, which uses distance as independent variable.

Ω a b

N intervals M respectively, and obtain a grid with (N + 1) × (M + 1) nodes. Noting J = {0, . . . , M}, we introduce the state variables (functions of time) (P i, j ) i, j∈I×J representing the mine depth at each node (x i , y j ) := (a+i∆x, c+ j∆y). The mine profile at time-frame k is now a surface represented by the set of points P k := (P k i, j (0)). We introduce the (N +1)×(M +1) non-negative controls U k i, j ≥ 0 , i, j ∈ I ×J, with the same dynamics Ṗk i, j = U k i, j .

P 0 P 1 P 2 X i P k i P k-1 i P k i+1 U k i S k i
Initial profile conditions are written as:

P 0 i, j (0) = p 0 (x i , y j ) , i, j ∈ I × J. (4.7)
The objective and capacity limit are similar to the 2D case, except that the integrals of g and e between two consecutive profiles are now in 3D instead of 2D. The relevant implementation details are provided in 4.6.

The main adjustment concerns the slope constraint: for each point P i, j of the profile we now choose to consider two slopes S i, j and T i, j , in the x-axis and y-axis directions respectively. Using the same basic forward finite differences as in 2D, we obtain the two sets of slope constraints at time-frame k:

-1 ≤ P k i+1, j (t) -P k i, j (t) κ x i , y j , P k i, j (t) ∆x ≤ 1, ∀i = 0, . . . , N -1, j = 0, . . . , M -1, ∀t ∈ [0, 1]. (4.8) -1 ≤ P k i, j+1 (t) -P k i, j (t) κ x i , y j , P k i, j (t) ∆y ≤ 1, ∀i = 0, . . . , N -1, j = 0, . . . , M -1, ∀t ∈ [0, 1]. (4.9)
Note that more sophisticated choices could be used for the slopes, such as centered differences formulas or increasing the number of slopes considered at each point. The formulation of (SOP) 3D SC is summarized below, with Fig. 4.3 illustrating the profile discretization in the 3D case with N = 5 and M = 3.

(SOP) 3D SC                                                            max T ∑ k=1 P k P k-1 g(x, y, z) (1 + α) k-1 dxdydz Ṗk i, j (t) = U k i, j (t), ∀(i, j) ∈ I × J, k = 1, . . . , T, t ∈ [0, 1] -1 ≤ S k i, j (t) κ x i ,y j ,P k i, j (t) ≤ 1, ∀(i, j) ∈ I × J, k = 1, . . . , T, t ∈ [0, 1] -1 ≤ T k i, j (t) κ x i ,y j ,P k i, j (t) ≤ 1, ∀(i, j) ∈ I × J, k = 1, . . . , T, t ∈ [0, 1] P k P k-1 e(x, y, z)dxdydz ≤ C k , k = 1, . . . , T P 0 i, j (0) = p 0 (x i , y j ), ∀(i, j) ∈ I × J Uk i, j (t) ≥ 0, ∀(i, j) ∈ I × J, k = 1, . . . , T, t ∈ [0, 1]

Analysis and optimality conditions for FOP

We study the final open pit problem in continuous formulation (FOP) by applying Pontryagin's Maximum Principle (Ref. [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF]), and look at the possible control structure of optimal profiles. Optimality conditions for (SOP) are not detailed here, and are more involved in particular due to the state constraint P ≤ P 0 being generalized over the sequence of time-frames, i.e. P i ≤ P i-1 , i = 1, . . . , N.

Applying Pontryagin's Maximum Principle

Following the formulation in Ref. [START_REF] Bonnans | Revisiting the analysis of optimal control problems with several state constraints[END_REF], we now state the PMP for (FOP). We denote y the state variables, p the associated costate variables, l the running cost, f the dynamics and h the state constraint. In all the following we assume the so-called normal case, i.e. the multiplier associated to the cost is nonzero and View is from 'above', with the state variable P k i, j giving the profile depth at node (x i , y j ) at time-frame k. Slopes S k i, j , T k i, j with neighbors along the x-axis and y-axis must be smaller than the maximal allowed slopes given by the function κ. The control U k i, j (along the z-axis) corresponds to the excavated depth from the same profile node at the previous time-frame P k-1 i, j . can be normalized to 1. Let us define the pre-Hamiltonian, omitting the argument t of functions u, y, p for clarity: Then for any local optimum ( ȳ, ū), there exists a non-trivial set of multipliers ( p, μ) such that the following relations are satisfied: i.e. μ is an non-decreasing function and is constant when the state constraint is not active.

Ω x y 0 N M P k i,j P k i+1,j P k i,j+1 S k i,j T k i,j U k i,j
H(t,
i) Adjoint equation d pP (t) = (g(t, P) -pP (t) ū(t)κ P (t, P) -pc (t)) dt + d μ(t) (4.
Remark 4.3.1. The state constraint is of order 1 since the control appears in its first time derivative ḣ = -uκ. We refer the reader to, for instance, Ref. [START_REF] Bonnans | Revisiting the analysis of optimal control problems with several state constraints[END_REF] for a more in-depth analysis of state constraints in the PMP framework, and especially the so-called "alternate adjoint" formulation.

Inactive case: bang/singular control

We start by studying the case when the state constraint is not active. Since per (4.16) the optimal control minimizes the pre-Hamiltonian which is linear in the control, solutions typically consist in a sequence of bang (saturated control) and/or singular control arcs. We introduce the switching function whose sign will determine the optimal control

ψ(t) := H u (t) = p P (t)κ(t, P(t)) (4.18)
As κ has strictly positive values we obtain the control law:

ū(t) =      1 if p P (t) < 0 -1 if p P (t) > 0 u s (t) if p P (t) = 0 over an interval (4.19)
The value of the singular control u s is traditionally determined from the fact that ψ and all its time derivatives vanish over a singular arc.

Over a singular arc, ψ vanishes and the first time derivative of the switching function can be reduced to ψ(t) = (g(t, P(t))p c )κ(t, P(t)) (4.20)

and similarly, by plugging ψ(t) = 0 in the second derivative and recalling ṗc = 0, we obtain ψ(t) = (g t (t, P(t)) + uκg P (t, P(t)))κ (4.21)

Solving ψ(t) = 0 for the singular control leads to

u s (t) = - g t (t, P(t)) κg P (t, P(t)) (4.22)
We can now derive the two following lemmas concerning singular arcs.

Lemma 4.3.1. A singular arc is not admissible when | g t (t,P) κg P (t,P) | > 1, and in particular when g P (t, P) = 0.

Proof. Immediate consequence of (4.22) and the control constraint u ∈ [-1, 1] Lemma 4.3.2. Let P be an optimal profile solution of (FOP), then, over a singular arc the curve (t, P(t)) follows the geodesic of g. Moreover, when maximal capacity is not reached, singular arcs follow more specifically the geodesics of null gain g = 0.

Proof. From (4.20), over a singular arc the equation ψ = 0 indicates that the derivative ġ := g t (t, P(t)) + Ṗg P (t, P(t)) vanishes, therefore the mine profile will follow the geodesics of g. If the maximal capacity constraint is not active, then the associated costate p c is zero (see (4.15)), and ψ = 0 then gives g = 0.

These lemmas expand the analysis of singular arcs obtained in Chapter 3 with calculus of variations techniques.

Active state constraint case

Over a constrained arc, the control u c is such that the constraint remains active, i.e. h = P 0 -P = 0, leading to the expression

u c (t) = Ṗ0 (t) κ(t, P(t)) (4.23)
Note that a constrained arc can only occur if the u c is admissible, i.e. | Ṗ0 (t)| ≤ κ(t, P(t)). This simply means that the initial profile must satisfy the maximal slope constraint.

Control structure summary

To summarize, the optimal profile, in terms of control structure, is a sequence of bang, singular and/or constrained arcs. Constrained arcs are where the profile is the same as the initial one, meaning there was no further excavation on these parts of the domain. Bang arcs correspond to the parts of the profile where the slope reaches its maximal allowed value, i.e. the digging is as steep as possible. Singular arcs, on the other hand, follow the geodesics of the gain function, meaning the gain is constant along these parts of the profile. Moreover, if the capacity limit is not reached, then this geodesic is more specifically the one of null gain, i.e. the digging stops where excavation is not profitable anymore.

Optimality conditions for the semi-continuous formulation are more involved and remain to be investigated thoroughly. The main complications arise from the spatial discretization of the profile, leading to maximal slope limits now being state constraints that involve adjacent nodes variables (including controls).

Numerical simulations

We present now the numerical simulations that illustrate the continuous and semi-continuous formulations of the Open Pit problem. After a brief description of the algorithms used for the global and local optimizations, we detail three test cases. First is the 2D FOP with limited capacity, that we solve with the continuous (both global and local optimization) and semi-continuous formulation (local optimization). The second example is the 2D SOP with limited capacity, for which we compare the results of both continuous and semi-continuous formulations (both with local optimization). Finally, we present a test case for the 3D SOP problem with the semi-continuous formulation. All simulations were carried out on a standard laptop, with numerical settings for all methods recalled in Table 4.1 p.39.

Numerical methods

Global optimization: FOP with continuous formulation

The Final Open Pit problem is low-dimensional, with only two state variables (not counting the running cost) and one control variable. Therefore it makes sense to try a global optimization method such as dynamic programming, or the so-called Hamilton-Jacobi-Bellman (HJB) approach. We use here the software BOCOPHJB [START_REF] Bonnans | Bocop -A collection of examples[END_REF], and refer to for instance Ref. [START_REF] Falcone | Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations[END_REF] for a detailed description of the HJB method. In this approach the value function of a fully discretized (time, state and control variables) version of the problem is computed, with the global optimum then being reconstructed from this information.

Local optimization: FOP and SOP with continuous and semi-continuous formulations

Since the numerical cost of the global method is too high for the Sequential Open Pit problem, we also use a local optimization method, namely the direct transcription approach. This method approximates the original (OCP) problem by a discretized reformulation as a nonlinear optimization problem (NLP), using a discretization of the time interval. We refer interested readers to for instance Ref. [START_REF] Betts | Practical methods for optimal control using nonlinear programming[END_REF] for a review of direct methods. We use here the software BOCOP [START_REF] Bonnans | BOCOP -A toolbox for optimal control problems[END_REF], based on the solver IPOPT [START_REF] Waechter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF] with sparse derivative computed by the automatic differentiation tool CPPAD [START_REF] Bell | CppAD: a package for C++ algorithmic differentiation[END_REF]. This local optimization method is also used for the semi-continuous formulation with explicit discretization of the space domain.

Numerical settings

In Table 4

2D FOP with global optimization for continuous approach

The solution obtained by the global optimization is displayed in Figure 4.4. At first glance, the control structure seems to be of the form Constrained-Bang-Bang-Constrained. On both sides the constraint P = P 0 is active, meaning there is no additional digging from the initial profile. In the middle, digging occurs with maximal slope, leading to the two bang arcs. Remark 4.4.3. It is worth noting that an estimate of the PMP costate can be derived from the gradient of the value function computed by the global method, see for instance Refs. [START_REF] Clarke | The relationship between the maximum principle and dynamic programming[END_REF][START_REF] Cristiani | Initialization of the shooting method via the Hamilton-Jacobi-Bellman approach[END_REF]. In the present case however, the gradient turns out to be quite noisy and of little practical use. This could be improved by increasing the discretizations, although the increase in computational times would not be competitive with respect to using a direct method.

2D FOP with local optimization for continuous approach

The solution from the local optimization is displayed in Figure 4.5 with the optimal profile and control as well as the PMP costate check. This solution is actually extremely close to the one in section 4.4.2, which indicates that the direct method actually found the global optimum as well, with the benefit of a more accurate solution. In particular, we can here clearly see that the first two arcs with nonzero control around x = 200 are not bang arcs since |u| < 1: they are actually part of the first constrained arc and correspond to the region where P 0 varies, thus the control u c (t) = Ṗ0 (t) κ(t,P(t)) from (4.23) is not just zero. Moreover, we can now check that the Constrained-Bang-Bang-Constrained control structure is consistent with the switching function and the path constraint. We observe a perfect match between the adjoint estimate from the discretized problem and the recomputed PMP costate. Figure 4.5 shows the value of the state constraint h = P 0 -P and its associated multiplier dµ. We retrieve dµ from the multiplier of the state constraint in the discretized problem (the correspondence can be inferred from comparing the expression of the PMP Hamiltonian and the Lagrangian of the NLP problem). In accordance with (4.17), the multiplier dµ is positive, and null when the constraint is not active. We also observe that the costate p P is continuous at the junctions between bang and singular arcs, while the control is discontinuous. Remark 4.4.4. In this particular case, the solution has no singular arcs, which is due to the capacity limit that prevents reaching the null gain region. The examples with unlimited capacity in 4.7 and 4.7.2 illustrate solutions with singular arcs where the optimal profile follows the geodesic g = 0. We finally present the solution obtained for the same problem using the semi-continuous formulation with a single phase (i.e T = 1). As can be seen in Figure 4.6 and Table 4.2, the solution is similar to the global and local optimizations using the continuous approach, with close values for the objective. CPU times are of the same order of magnitude for the two local optimizations with continuous and semi-continuous formulation, while global optimization is significantly slower (two orders).

Method

Sequential Open Pit (2D and 3D): local optimization

In this section we present solutions for the Sequential Open Pit. First we solve a 2D example using both the continuous and semi-continuous formulations. Then we show a solution for a more realistic 3D problem using the semi-continuous formulation. To our knowledge, this is the first attempt to tackle the 3D case in an optimal control framework. Continuous approach. We solve the 2D SOP problem for 12 time-frames, with a constant function κ = 1 , a rate α = 0.1 and a maximal capacity c max = 1e4 for each time-frame. The solution indicates that most of the excavation effort is concentrated in the high gain regions of the domain, which is not surprising. Semi continuous approach. We now solve the same SOP problem with the semi-continuous formulation. Both approaches give similar solutions, as can be seen in Figure 4.7. The objective values showed in Table 4.3 are quite close with a difference of 1.3%, while CPU times are in the same order of magnitude.

Method

Objective CPU Continuous 89939 31s Semi-continuous 91153 43s Table 4.3: Solutions for the 2D SOP problem.

3D SOP with semi continuous approach

For the 3D case we consider a domain Ω = [0, 1200] × [0, 400] and an analytical gain density function stated by

g(x, y, z) = 1000 -(x -600) 2 + (y -200) 2 + (z -350) 2 (4.24)
that reaches its maximal value in (600, 200, 350) and which decreases radially from this point. The initial profile used in this instance is showed in Figure 4.8. We set a discount factor α = 0.1 and solve the 3D SOP for different capacity limits and number of time-frames, using a 30 × 10 discretization of Ω. Figures 4.9 4. [START_REF] Amaya | Optimality conditions for the continuous model of the final open pit problem[END_REF] shows the objective values and CPU times. Results are consistent overall, with solutions trying to reach the region of highest gain as fast as allowed by the slope and capacity constraints. Increasing the capacity limit and / or the duration of the time interval both yield better objective values, as expected. CPU times are still reasonable, with the longest run at 139s. Capacity limit: 10

Conclusions

In the present chapter we focused on the Open Pit problem in an optimal control framework. We extended some previous results on the optimality conditions for the Final Open Pit, and introduced a new semi-continuous formulation that handles the 3D profile sequential optimization. Numerical simulations are provided for the continuous and semi-continuous approaches on several test cases. The 2D FOP case showed a good consistency between global and local optimization for the continuous approach, as well as local optimization for semi-continuous, and matched the optimality conditions from Pontryagin's Principle. Then the 2D SOP case again indicated a good match for the continuous and semi-continuous formulations. Finally we solved a 3D SOP test case, to our knowledge for the first time in an optimal control framework. Perspectives in the continuation of the present work include solving a more complete 3D SOP example using 3D interpolated data for the gain and maximal slope, as well as studying the optimality conditions for the semi-continuous approach. The latter could prepare for the use of indirect shooting methods such as HAMPATH [START_REF] Caillau | Differential pathfollowing for regular optimal control problems[END_REF], especially since the local optimization method used here can provide the knowledge of the optimal control structure and a costate approximation.

Implementation details for the semi-continuous approach

Time discretization. The Sequential Open Pit for the semi-continuous approach described in 4.2.2 is a multi phase problem. Instead of duplicating the variables for each time-frame, we use here in practice a more compact implementation, by using a time step ∆t of 1 time-frame, i.e. the time discretization t k = 0 . . . T is the sequence of time-frames. This choice makes sense from the operational point of view, since the sequential open pit planning precisely consists in determining the optimal mine profile at each time-frame. It also simplifies a lot the computation of the integrals of the gain and effort functions between two successive mine profiles. We choose an implicit Euler scheme for the time discretization, which gives the trivial discrete dynamics

P k+1 i = P k i +U k+1 i (4.25)
that easily gives the next / previous mine profile when needed in the computations.

Gain. An additional state variable g is added to represent the gain realized along the time-frames, whose dynamics can be written as

ġ(t k ) = 1 (1 + α) k-1 P k P k-1 g(x, z)dxdz, ∀k = 1, . . . , T (4.26)
The objective is then to maximize g(T ). For the 2D case, we approximate the 2-dimensional integral of g by trapezoidal rule over x then along z. In the 3D profile case, the 3D integral of g for the computation of the gain is approximated using a 2D trapezoidal rule along (x, y) then a standard trapezoidal rule along z. Capacity. At each time-frame, the integral of the excavation effort over the domain Ω can be approximated by

P k P k-1 e(x, z)dxdz ≈ N ∑ i=0 ∆x P k i P k-1 i e(x i , z)dz (4.27)
Since E = 1 and from the discrete dynamics P k i = P k-1 i +U k i , we can use the following formula

P k P k-1 e(x, z)dxdz ≈ ∆x N-1 ∑ i=0 U k i . (4.28)
Similarly, for the 3D profile case, the excavation effort at time-frame k is approximated as

P k P k-1 e(x, y, z)dxdydz ≈ ∆x∆y N-1 ∑ i=0 M-1 ∑ j=0 U k i, j . (4.29)

Additional examples for the final open pit -continuous formulation 4.7.1 FOP with infinite capacity and constant slope

We show here the basic example with unconstrained capacity, namely c max = ∞. Fig. 4.11 shows the solution obtained by the global method, and Fig. 4.12 shows the solution from the local method, and we observe that both solutions match. With infinite capacity, the solution, as expected, digs as much as possible with respect to the maximal slope, until it reaches negative gain. This corresponds to the observed Bang-Singular-Bang control structure (neglecting the two very small constrained arcs P = P 0 = 0 at the extremities). As stated in Lemma 2, the singular arc in the middle follows the geodesic g = 0. The corresponding control also matches the theoretical expression of the singular control (4.19), despite some oscillations at the junctions with the bang arcs.

FOP with infinite capacity and variable slope

Here we illustrate a case with a non-constant maximal slope κ. For this, we consider following arbitrary κ function:

κ(x, z) =      0.5 x ∈ [0, 330) 1 x ∈ [330, 960) 5 x ∈ [960, 1230] (4.30)
We chose a piece-wise constant function so that the assumption in section 4.3 is satisfied almost everywhere. The solutions from the global and local optimizations are shown in Fig. 4.13 and Fig. 4.14 respectively. This time we obtain a control structure that includes all possible types of arcs: Bang-Singular-Bang-Constrained. The main difference compared to the constant slope case is that the optimal profile digs less ground on the right side region where the gain is negative, as a steeper slope is allowed there.

As for the previous examples we observe that the singular control and costate from the solution closely match their formal expressions from the PMP.

Compartmental models

Since the Covid-19 outbreak, numerous research works have coupled optimal control techniques with epidemiological models for the study and management of the pandemic, such as reducing the number of infected people (see for example [START_REF] Bliman | How best can finite-time social distancing reduce epidemic final size?[END_REF][START_REF] Grigorieva | Optimal quarantine strategies for COVID-19 control models[END_REF][START_REF] Kantner | Beyond just "flattening the curve": Optimal control of epidemics with purely non-pharmaceutical interventions[END_REF][START_REF] Ketcheson | Optimal control of an SIR epidemic through finite-time non-pharmaceutical intervention[END_REF][START_REF] Kovacevic | A Distributed Optimal Control Model Applied to COVID-19 Pandemic[END_REF][START_REF] Palmer | Optimal control of COVID-19 infection rate with social costs[END_REF]). In the context of this PhD, I had the opportunity to contribute to the study of the so-called overcrowding problem, which is a significant issue for hospitals. More precisely, our main objective was to minimize the peak of infections in order to keep the health system below saturation levels.

We worked with the most widespread compartmental model in epidemiology called the SIR model, which we recall in the following section. More specifically we started from the work of Morris et al. [START_REF] Morris | Optimal, near-optimal, and robust epidemic control[END_REF] on peak minimization for the SIR model, which despite not using optimal control techniques per se, is closely related to control theory in its formulation.

After the classical SIR model we introduce a new compartmental model that takes vaccination into account. This model was developed in the framework of a Pan American Health Organization 1 project by a joint team of mathematicians and doctors devoted to the evaluation of the vaccination campaign of which I was a part of. We hope that this new model can be used in future works.

SIR model

The so-called SIR model introduced in 1927 by Kermack and McKedrick [69] is the most widespread model for the study of directly transmitted infectious diseases that spread through contact from an individual to another, such as, for instance, flu, tuberculosis, MERS-CoV and Covid-19.

The model derives its name from partitioning a given population into three compartments, namely Susceptible, Infected and Recovered individuals. We note respectively S(t), I(t) and R(t) the size at time t ∈ [0, T ] of each group. Many assumption can be made modifying the final model, we consider here the case of a constant population, i.e. without any natural births or deaths. Normalizing the total population, we have the relation

S(t) + I(t) + R(t) = 1, ∀t ≥ 0.
The dynamic of the three groups follows the ODE system:

Ṡ = -β SI İ = β SI -γI Ṙ = γI, (5.1)
with constants β > 0 and γ > 0 corresponding to transmission and recovery rates respectively. The so-called basic reproductive number R 0 , defined as the average number of secondary infection cases 1 In Spanish Organización Panamericana de la Salud (OPS) 51 caused by a single primary case in a susceptible population, is therefore given by R 0 = β /γ in the SIR model. An obvious property of the model is that a disease will spread among the population when R 0 > 1, making this case the most relevant for study. We refer readers interested in more properties of SIR models and basic reproductive number to for instance [START_REF] Anderson | Discussion: the Kermack-McKendrick epidemic threshold theorem[END_REF][START_REF] Jones | Notes on R0[END_REF][START_REF] Kermack | A contribution to the mathematical theory of epidemics[END_REF][START_REF] Weiss | The SIR model and the foundations of public health[END_REF].

In order to fight the pandemic, local governments have implemented a variety of measures to reduce the transmission of the disease. These efforts can typically be categorized as non-pharmaceutical, such as lock-downs and quarantines, and pharmaceutical such as vaccines and treatments. The first kind can be introduced in the SIR model as a control variable b(t) ∈ [0, 1] applied as a factor to the transmission rate beta:

Ṡ = -b(t)β SI İ = b(t)β SI -γI Ṙ = γI (5.2)
In this case, b = 1 corresponds to no intervention, while an identically null control b = 0 will completely contain the infection, however, this is a hardly applicable policy at large scale, since it would amount to cutting off contact between individuals completely. This is why a more realistic constraint is usually imposed on the control, for example, in Chapters 6 and 7 we considered a maximal budget limit Q imposed on the L 1 norm of the policy intensity 1b. This is somewhat less restrictive that the original constraint in [START_REF] Morris | Optimal, near-optimal, and robust epidemic control[END_REF] that limited the control effectively applied (i.e b < 1) to a certain time window [t i ,t i + τ] to be determined. With this constraint, the authors main result on the peak minimization, i.e. min b max t∈[0,T ] I(t), in [START_REF] Morris | Optimal, near-optimal, and robust epidemic control[END_REF] was:

Proposition 5.1.1. The optimal intervention which is a solution of min b max t∈[0,T ] I(t) is: b opt (t) =    γ β S(t) , t ∈ [t i ,t i + f τ) 0, t ∈ [t i + f τ,t i + τ] for a certain f ∈ [0, 1].
This was our starting point to study the peak minimization in the SIR model, and in Chapter 7 we present our results and compare them to this original work.

Regarding pharmaceutical measures, we studied how to include vaccine programs in compartmental models. The result of this research will be presented in the next section and it corresponds to an original model posed by our team.

A new compartmental model including vaccines

We introduce here an extension of the SIR model that includes the effect of vaccination campaigns. This model was developed in the framework of a project supported by the OPS and the MINSAL (Chile ministry of health). The aim of this model is to evaluate the effectiveness of different.

Many extensions of the original SIR model, typically with more compartments, have been proposed, see for instance [START_REF] Angelov | Optimal vaccination strategies using a distributed epidemiological model applied to COVID-19[END_REF][START_REF] Aronna | A model for COVID-19 with isolation, quarantine and testing as control measures[END_REF][START_REF] Giordano | Modeling vaccination rollouts, SARS-CoV-2 variants and the requirement for non-pharmaceutical interventions in Italy[END_REF][START_REF] Ramos | Modeling the impact of SARS-CoV-2 variants and vaccines on the spread of COVID-19[END_REF]. The main novelty of the model presented here is the partition of the individuals into three so-called fundamental parts, according to their vaccination status: (i) not fully vaccinated, (ii) fully vaccinated and (iii) fully vaccinated with additional booster injection. This division follows the idea that the disease dynamics are altered by vaccination. Each fundamental part is further modeled with 6 compartments (state variables): susceptible (S), asymptomatic infected (I a ), symptomatic infected (I s ), asymptomatic recovered (R a ), symptomatic recovered (R s ) and deceased (D). To distinguish variables from the three fundamental parts, we add the subscripts N (non fully vaccinated), V (fully vaccinated) and V r (vaccinated with booster). We suppose the population modeled by the base structure is homogeneous (meaning there is no age structure in the model), isolated, susceptible to only one variant circulating and constant (no natural birth or death). We denote the transmission rates by β k i for i ∈ {N,V,V R } and k ∈ {a, s}, using a for asymptomatic and s for symptomatic individuals. For example, β a N represents the transmission rate after contact with a non-vaccinated, asymptomatic, infected person. The overall transmission rate for the susceptible partition S i with i ∈ {N,V,V R } will noted Λ i defined as:

S I

a I s R a R s D q a R a q s R s φ Λ (1 -φ )Λ δ I s γ s I s γ a I a
Λ i (t) = (1 -f i (t)) ∑ j∈{N,V,V r } ((1 -u(t))β a i, j I a j (t) + (1 -µ)β s i, j I s j (t)) S i (t) i ∈ {N,V,V r }.
Here µ ∈ [0, 1] represents the portion (assumed constant) of infected, symptomatic population that is completely isolated. The factor u(t) ∈ [0, µ] corresponds to non pharmaceutical measures, similar to b(t) in the basic SIR model except it only affects now the asymptomatic individuals, since the symptomatic portion are assumed to be in maximum isolation. We also denote by f i (t) the vaccine effectiveness to reduce transmission, being f N (t) ≡ 0 for the nonfully vaccinated part. Noting e S N (τ), e R a N (τ), e R s N (τ) ∈ [0, 1] the effectiveness for the compartments (S N ), (R a N ) and (R s N ), respectively, τ days after becoming fully vaccinated, and v S N (t), v R a N (t) and v R s N (t) the vaccination rates for each compartment, then the overall effectiveness f i is modeled similarly as in [START_REF] Ramos | Modeling the impact of SARS-CoV-2 variants and vaccines on the spread of COVID-19[END_REF] by the equation:

f V (t) = t-t V 0 e S N (τ)v S N (t -τ) + e R a N (τ)v R a N (t -τ) + e R s N (τ)v R s N (t -τ) dτ t-t V 0 v S i (s) + v R a i (s) + v R s i (s)ds t ≥ t V ,
where t V is the start of the vaccination plan.

Similarly, for a booster we note r S V (τ), r R a V (τ), r R s V (τ) ∈ [0, 1] the effectiveness in (S V ), (R a V ) and (R s V ), respectively, τ days after applying the booster dose and v S V (t), v R a V (t) and v R s V (t) the booster injection rates for each relevant compartment, then the total effectiveness of a booster dose is given by:

f V r (t) = t-t V 0 r S V (τ)v S V (t -τ) + r R a V (τ)v R a V (t -τ) + r R s V (τ)v R s V (t -τ) dτ t-t V 0 vV (t)(s)ds t ≥ t V ,
Finally we describe the remaining parameters of the model:

• Recovery rates: We denote γ k i , with i ∈ {N,V,V r } and k ∈ {a, s} the recovery rates of infected individuals who are in the fundamental part i and could be symptomatic (s) or asymptomatic (a). A natural assumption is that recovery should be quicker for vaccinated people (as reported for example in [START_REF] Maclean | Why do breakthrough COVID-19 infections occur in the vaccinated?[END_REF]), leading to the following inequalities

γ k N ≤ γ k V ≤ γ k V r k ∈ {a, s}.
• Death rates: We denote δ i , with i ∈ {N,V,V r } the death rates in the infected population in fundamental part i. Once again assuming vaccination reduces the mortality (as the evidence show in [START_REF] Arbel | Effectiveness of a second BNT162b2 booster vaccine against hospitalization and death from COVID-19 in adults aged over 60 years[END_REF][START_REF] Lv | Mortality rate and characteristics of deaths following COVID-19 vaccination[END_REF][START_REF] Kennedy | Vaccine hesitancy: a growing concern[END_REF][START_REF] Roghani | The influence of COVID-19 vaccination on daily cases, hospitalization, and death rate in Tennessee, United States: Case Study[END_REF]), we have the inequalities

δ V r ≤ δ V ≤ δ N .
• Loss of immunity rates : The rate of loss of immunity for a recovered individual will be denoted by q k i with i ∈ {N,V,V r } and k ∈ {a, s}.

• Probability of being asymptomatic: The probability of a susceptible individual, belonging to the fundamental part i ∈ {N,V,V r }, to be asymptomatic after infection will be denoted by

φ i ∈ [0, 1]. These probabilities should fulfill φ N ≤ φ V ≤ φ V r ,
because the probability of being asymptomatic should increase as vaccination advances as it is shown in [START_REF] Aguiar | The role of mild and asymptomatic infections on COVID-19 vaccines performance: a modeling study[END_REF][START_REF] Tang | Asymptomatic and symptomatic SARS-CoV-2 infections after BNT162b2 vaccination in a routinely screened workforce[END_REF] .

Table 5.1 shows a summary of the parameters used in this model.

S N I a N I s N R a N R s N D N S V I a V I s V R a V R s V D V S V r I a V r I s V r R a V r R s V r D V r q a N R a N q s N R s N φ N Λ N (1 -φ N )Λ N δ N I s N γ s N I s N γ a N I a N q a V R a V q s V R s V φ V Λ V (1 -φ V )Λ V δ V I s V γ s V I s V γ a V I a V q a Vr R a Vr q s Vr R s Vr φ Vr Λ Vr (1 -φ Vr )Λ Vr δ Vr I s Vr γ s Vr I s Vr γ a Vr I a Vr v SN v R a N v R s N v SV v R a V v R s V Figure 5
.2: Base structure with 3 fundamental parts and 18 compartments in total. Green arcs indicate vaccination actions that link the three parts, red arcs correspond to loss of immunity.

Chapter 6

A Mayer formulation for peak minimization problems

This chapter corresponds to the accepted pre-print [START_REF] Molina | Equivalent formulations of optimal control problems with maximum cost and applications[END_REF] titled "Equivalent formulations of optimal control problems with maximum cost and applications" to appear in Journal of Optimization Theory and Applications.

Introduction

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

u(•) ess sup t∈[t 0 ,T ] y(t)
where y(t) = θ (t, ξ (t)) and ξ (•) is the solution of a controlled dynamics ξ = φ (ξ , u), ξ (t 0 ) = ξ 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-Bellman Jacobi equation).

However, this problem falls into the class of optimal control with L ∞ criterion, for which several characterizations of the value function have been proposed in the literature [START_REF] Barron | The Bellman equation for minimizing the maximum cost[END_REF][START_REF] Barron | The L ∞ control problem with continuous control functions[END_REF][START_REF] González | A Bellman's equation for minimizing the maximum cost[END_REF]. Typically, the value function is solution, in a general sense, of a variational inequality of the form

min ∂ t V + inf u ∂ ξ V.φ (x, u) , V -θ = 0
without boundary condition. Nevertheless, although necessary optimality conditions and numerical procedures have been formulated [START_REF] Barron | The Pontryagin maximum principle for minimax problems of optimal control[END_REF][START_REF] Marco | A numerical procedure for minimizing the maximum cost[END_REF][START_REF] Marco | Minimax optimal control problems. Numerical analysis of the finite horizon case[END_REF][START_REF] Gianatti | Solving minimax control problems via nonsmooth optimization[END_REF], 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 [START_REF] Jones | Extensions of the dynamic programming framework: Battery scheduling, demand charges, and renewable integration[END_REF][START_REF] Jones | A generalization of Bellman's equation with application to path planning, obstacle avoidance and invariant set estimation[END_REF]. 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 6.2, we establish the setup and the hypotheses of this article, and define the problem. In Section 6.3, we provide equivalent formulations of the stud-57 ied problem in the form of two Mayer problems with fixed initial condition, and under state or mixed constraint. In Section 6.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 6.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 also compare numerically with L p approximations. Finally, we discuss in Section 6.6 about the potential merits of the different formulations as practical methods to compute optimal solution of L ∞ control problems.

Problem and hypotheses

We shall consider autonomous dynamical systems defined on a invariant domain D of R n+1 of the form ẋ = f (x, y, u) ẏ = g(x, y, u) (6.1) (where g is a scalar function) where the values of the control u(•) belong to a given set U ⊂ R p . More specifically, throughout the chapter, 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 D ×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)| ⩽ C(1 + ||x|| + |y|), (x, y) ∈ D, u ∈ U
For instance, y(•) can be a smooth output of a dynamics

ẋ = f (x, u), y = h(x)
which can be rewritten as ẋ

= f (x, u) ẏ = g(x, u) := ∇h(x) T • f (x, u)
Let U be the set of measurable functions u(•) : [0, T ] → U and consider (x 0 , y 0 ) ∈ D, T > 0. Under the usual arguments of the theory of ordinary differential equations, Assumption 1 ensures that for any u(•) ∈ U there exists a unique absolutely continuous solution (x(•), y(•)) of (6.1) on [0, T ] for the initial condition (x(0), y(0)) = (x 0 , y 0 ) (see for instance [START_REF] Fleming | Deterministic and stochastic optimal control[END_REF]). Define then the solutions set

S := {(x(•), y(•)) ∈ A C ([0, T ], R n+1 ), sol. of (6.1) for u(•) ∈ U 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(•):

P : inf u(•)∈U max t∈[0,T ] y(t) = inf (x(•),y(•))∈S max t∈[0,T ]
y(t) .

Formulations with constraint

A first approach considers the family of constrained sets of solutions

S z := {(x, y) ∈ S , y(t) ⩽ z, t ∈ [0, T ]}, (z ∈ R)
and to look for the optimization problem inf{z; S z ̸ = / 0}.

This problem can be reformulated as a Mayer problem

P 0 : inf u(•)∈U z(T )
for the extended dynamics in

D × R    ẋ = f (x, y, u) ẏ = g(x, y, u) ż = 0 under the state constraint C : z(t) -y(t) ⩾ 0, t ∈ [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

D × R with an additional control v(•) with values in [0, 1]    ẋ = f (x, y, u) ẏ = g(x, y, u) ż = max(g(x, y, u), 0)(1 -v) (6.2)
Let V be the set of measurable functions v : [0, T ] → [0, 1]. Note that under Assumption 1, for any (x 0 , y 0 , z 0 ) ∈ D × R and (u, v) ∈ U × V , there exists an unique absolutely solution (x(•), y(•), z(•)) of (6.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 : inf (u(•),v(•))∈U ×V z(T ) under the constraint C
and shows its equivalence with problem P. We first consider fixed controls u(•). and is reached for a control v(•) that takes values in {0, 1}.

Proof. From equations (6.2), one get that any solution z(•) is non decreasing, and as z satisfies the constraint z ⩾ y, we deduce that one has z(T ) ⩾ max t∈[0,T ] y(t) (6.5)

for any solution of (6.2), and thus

max t∈[0,T ] y(t) ⩽ inf v∈V z(T ) under the constraint z(t) ⩾ y(t), t ∈ [0, T ].
Let x(•), y(•) be the solution of (6.1) for the control u(•) and let I be the set of invisible points from the left of y, that is

I := t ∈ (0, T ); y(t ′ ) > y(t) for some t ′ < t .
Consider then the control

v(t) = 1, t ∈ int I, 0, t / ∈ int I (6.6)
When I is empty, y(•) is a non decreasing function and, when v(t) = 0 for all t ∈ [0, T ], one has z(t) = y(t) for any t ∈ [0, T ]. Therefore, one has

z(T ) = y(T ) = max t∈[0,T ] y(t).
When I is non empty, there exists, from the sun rising Lemma [START_REF] Tao | An introduction to measure theory[END_REF], 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 ̸ = T , -if b n = T , then y(a n ) ⩾ y(b n ).
Note that when t / ∈ int I, one has y(t) ⩾ y(t ′ ) for any t ′ ⩽ t. Therefore, the solution z with control (6.6) verifies z(t) = y(t), t / ∈ int I y(a n ), t ∈ I n for some n (see Figure 6.1 as an illustration). Let t ∈ [0, T ] be such that

y(t) = max t∈[0,T ] y(t),
which implies that any point t ′ > t in [0, T ] is invisible from the left. Then, one has z(T ) = z(t) ⩽ y(t). Thus, from (6.5), we obtain max

t∈[0,T ] y(t) = z(T )
and deduce max However, there is no uniqueness of the optimal control v(•). Any admissible solution z(•) that is above y(•) and such that z(t) = ŷ for t ⩾ t = min{t ∈ (0, T ], y(t) = ŷ}, where ŷ := max s∈[0,T ] y(s), is also optimal.

t∈[0,T ] y(t) = inf v(•)∈V z(T ) under the constraint C .
We then obtain the equivalence between problems P 1 and P in the following sense.

Proposition 6.3.2. If (u ⋆ (•), v ⋆ (•)) is optimal for Problem P 1 , then u ⋆ (•) is optimal for Problem P. Conversely, if u ⋆ (•) is optimal for Problem P, then (u ⋆ (•), v ⋆ (•)) is optimal for Problem P 1 where v ⋆ (•)
is optimal for the problem (6.3) for the fixed control u ⋆ (.).

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 (6.2), with a control v(•) 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

C m : max(y(t) -z(t), 0)(1 -v(t)) + z(t) -y(t) ⩾ 0, a.e. t ∈ [0, T ]
and the optimal control problem

P 2 : inf (u(•),v(•))∈U ×V z(T ) under the constraint C m .
Proposition 6.3.3. Problems P 1 and P 2 are equivalent.

Proof. One can immediately see that for any admissible solution that satisfies constraint C , the constraint C m is necessarily fulfilled as max(yz, 0) is identically null. Conversely, fix an admissible control u(•) and consider a control v(•) that satisfies C m . We show that this implies that the solution (y(•), z(•)) verifies necessarily z(t) ⩾ y(t) for any t ∈ [0, T ]. If not, consider the non-empty set

E := {t ∈ [0, T ]; z(t) -y(t) < 0},
which is open as zy is continuous. Note that one has ż(t)ẏ(t) ⩾ 0 for a.e. t ∈ E. Therefore zy is non decreasing in E and we deduce that for any t ∈ E, the interval [0,t] is necessarily included in E, which then contradicts the initial condition z(0) = y(0).

Formulation without state constraints

We posit Π = (x, y, z) ∈ D × R and consider the differential inclusion

Π ∈ F(Π) := (u,v)∈U×[0,1]   f (x, y, u) g(x, y, u) h(x, y, z, u, v)   (6.7) with h(x, y, z, u, v) = max(g(x, y, u), 0)(1 -v⊮ R + (z -y))
where ⊮ R + is the indicator function

⊮ R + (ζ ) = 1, ζ ⩾ 0 0, ζ < 0
Let Π 0 = (x 0 , y 0 , y 0 ) and denote by S ℓ the set of absolutely continuous solutions of (6.7) with Π(0) = Π 0 ∈ D × R. We consider the Mayer problem with (x(•), y(•)) solution of (6.1) for some control u(•) ∈ U , that in turn is optimal for problem P.

P 3 : inf Π(•)∈S ℓ z(T ).
Proof. We fix the initial condition Π(0) = Π 0 and consider the augmented dynamics

Π ∈ F † (Π) := (u,v,α)∈U×[0,1] 2   f (x, y, u) g(x, y, u) h † (x, y, z, u, v, α)   (6.8) with h † (x, y, z, u, v, α) = (1 -α)h(x, y, z, u, v) + α max w∈U h(x, y, z, w, 0).
Under Assumption 2, the values of F † are convex compact. One can straightforwardly check that the set-valued map F † is upper semi-continuous1 with linear growth. Therefore, the reachable set S † ℓ (T ) (where S † ℓ denotes the set of absolutely continuous solutions of (6.8) with Π(0) = Π 0 ) is compact (see for instance [START_REF] Aubin | Viability Theory[END_REF]Proposition 3.5.5]). Then, there exists a solution Π ⋆ (•) = (x ⋆ (•), y ⋆ (•), z ⋆ (•)) of (6.8) which minimizes z(T ).

Note that any admissible solution (x(•), y(•), z(•)) of system (6.2) that satisfies the constraint C m belongs to S ℓ ⊂ S † ℓ . We then get the inequality z ⋆ (T ) ⩽ inf{z(T ); (x(•), y(•), z(•)) sol. of (6.2) with C m }. (6.9)

Let us show that any solution

Π(•) = (x(•), y(•), z(•)) in S ℓ verifies z(T ) ⩾ max t∈[0,T ] y(t) (6.10)
We show that one has z(t) ⩾ y(t) for any t ∈ [0, T ]. We proceed by contradiction, as in the proof of Proposition 6.3.3. If the set E = {t ∈ (0, T ); z(t)y(t) < 0} is non-empty, one has ż(t)ẏ(t) ⩾ 0 for a.e. t ∈ E which implies, by continuity, that z(s)y(s) < 0 ∀s ∈ (0,t),t ∈ E and then z(0)y(0) < 0 which contradicts the initial condition z(0) = y(0). Moreover, as the map h is non-negative, z(•) is non decreasing and we conclude that (6.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(•), y(•)) is solution of (6.1) for a certain u(•) ∈ U . With(6.10), we obtain

z ⋆ (T ) ⩾ max t∈[0,T ] y ⋆ (t) ⩾ inf u∈U max t∈[0,T ] y(t); (x(•), y(•)) sol. of (6.1) (6.11)
where (x ⋆ (•), y ⋆ (•)) is solution of (6.1) for a certain u ⋆ (•) ∈ U . Finally, inequalities (6.9) and (6.11) with Propositions 6.3.2 and 6.3.3 show that z ⋆ (T ) is reached by a solution of (6.2) under the constraint C m , and that u ⋆ (•) is optimal for problem P. We also conclude that the optimal value z ⋆ (T ) is reached by a solution in S ℓ , which is thus optimal for problem P 3 . Remark 6.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 Π(•) = (x(•), y(•), z(•)) can be synthesized by a measurable control (u(•), v(•)). Proposition 6.4.1 shows that (x(•), y(•)) is indeed a solution of system (6.1) for a measurable control u(•), but one cannot guarantee a priori that z(•) can be generated by a measurable control v(•), which is irrelevant for our purpose.

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 6.5.

Let us consider the family of dynamics parameterized by θ > 0

   ẋ = f (x, y, u) ẏ = g(x, y, u) ż = h θ (x, y, z, u, v) (6.12)
with

h θ (x, y, z, u, v) = max(g(x, y, u), 0)(1 -v e -θ max(y-z,0) ).
Here the expression e -θ max(y-z,0) plays the role of an approximation of ⊮ R + (zy) when θ tends to +∞.

We then define the family of Mayer problems

P θ 3 : inf Π(•)∈S θ z(T )
where S θ denotes the set of absolutely continuous solutions Π(•) = (x(•), y(•), z(•)) of (6.12) for the initial condition Π(0) = Π 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 6.4.2. Under Assumption 2, for any increasing sequence of numbers θ n (n ∈ N) that tends to +∞, the problem P θ n 3 admits an optimal solution, and for any sequence of optimal solutions (x n (•),

y n (•), z n )(•)) of P θ n
3 , the sequence (x n (•), y n (•)) converges, up to sub-sequence, uniformly to an optimal solution (x ⋆ (•), y ⋆ (•)) of Problem P, and its derivatives weakly to ( ẋ⋆ (•), ẏ⋆ (•)) in L 2 . Moreover, z n (T ) is an increasing sequence that converges to max t∈[0,T ] y ⋆ (t).

Proof. As in the proof of Proposition 6.4.1, we consider for any θ > 0 the convexified dynamics

  Let ȳ := inf u∈U max t∈[0,T ] y(t); (x(•), y(•)) sol. of (6.1)
By Proposition 6.4.1, we know that there exists an optimal solution (x(•), y(•), z(•)) of problem P 3 such that z(T ) = ȳ. Clearly, this solution belongs to S θ for any θ , and we thus get

z ⋆ θ (T ) ⩽ ȳ. (6.13) Let F θ (Π) := (u,v)∈U×[0,1]   f (x, y, u) g(x, y, u) h θ (x, y, z, u, v)  
and note that one has lim

θ →+∞ d (F θ (Π), F(Π)) = 0, Π ∈ D × R (6.14)
Consider an increasing sequence of numbers θ n (n ∈ N), and denote Π n (•) = (x n (•), y n (•), z n (•)) an optimal solution of problem P θ n 3 . Note that one has

S θ n+1 ⊂ S θ n • • • ⊂ S θ 0 . (6.15)
Therefore, the sequence Πn (•) is bounded, and Π n (•) as well. As F is upper semi-continuous, we obtain that Π n (•) converges uniformly on [0, T ], up to a sub-sequence, to a certain

Π ⋆ (•) = (x ⋆ (•), y ⋆ (•), z ⋆ (•))
which belongs to S l (see for instance [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF]Th. 3.1.7]). From property(6.15), we obtain that z n (T ) is a non decreasing sequence that converges to z ⋆ (T ), and from (6.13), we get passing at the limit

z ⋆ (T ) ⩽ ȳ On another hand, (x ⋆ (•), y ⋆ (•), z ⋆ (•)
) belongs to S l and we get from Proposition 6.4.1 the inequality

z ⋆ (T ) ⩾ ȳ Therefore, one has z ⋆ (T ) = ȳ and (x ⋆ (•), y ⋆ (•), z ⋆ (•)
) is then an optimal solution of problem P 3 . From Proposition 6.4.1, we obtain that one has necessarily

z ⋆ (T ) = max t∈[0,T ] y ⋆ (t).
Finally, the sequence ( ẋn (•), ẏn (•)) being bounded, it converges, up to a sub-sequence, weakly to ( ẋ⋆ (•), ẏ⋆ (•)) in L 2 thanks to Alaoglu's Theorem.

Numerical illustrations

Our aim is to illustrate the different formulations on problems for which the optimal solution is known.

A particular class of dynamics

We consider dynamics of the form

(Σ) : ẋ = f (x) ẏ = g(x, u) x ∈ R n , u ∈ U Proposition 6.5.1. A feedback control x → φ ⋆ (x) such that g(x, φ ⋆ (x)) = min u∈U g(x, u), x ∈ R n
is optimal for problem P.

Proof. For a given x 0 in R n , let x(•) be the solution of ẋ = f (x), x(0) = x 0 independently to the control u(•). Then, for any solution y(•), one has

y(t) = y(0) + t 0 g(x(τ), u(τ)) dτ ⩾ y(0) + t 0 min v∈U g(x(τ), v) dτ, t ⩾ 0 Let y ⋆ (•) be defined as y ⋆ (t) := y(0) + t 0 min v∈U g(x(τ), v) dτ, t ⩾ 0 Clearly, one has max t y(t) ⩾ max t y ⋆ (t)
where y ⋆ (•) is a solution of Σ for any measurable control u ⋆ (•) such that g(x(t), u ⋆ (t)) = min v∈V g(x(t), v), a.e. t ⩾ 0.

We conclude that y ⋆ (•) is an optimal trajectory of problem P for the control generated by the feedback φ ⋆ .

As a toy example, we have considered the system

ẋ = 1, x(0) = 0 ẏ = (1 -x)(2 -x)(4 -x)(1 + u/2), y(0) = 0 u ∈ [-1, 1]
for which

φ ⋆ (x) = -sign (1 -x)(2 -x)(4 -x)
is an optimal control which minimizes max t∈[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

ẏ = (1 -t)(2 -t)(4 -t)(1 + u/2)
for which the open-loop control

u ⋆ (t) = -sign (1 -t)(2 -t)(4 -t)
is optimal.

For T = 5, we have first computed the exact optimal solution of problem P with the open-loop u ⋆ (•), by integrating the dynamics with Scipy in Python software (see Figure 6.2). Effects of perturbations 6.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 (Bocop software using Gauss II integration disturbance max t∈[0,T ] y(t) error 0 2.24985 0 0.001% 2.24985 4.10 -6 % 0.01% 2.25010 0.01% 0.1% 2.69457 20% Table 6.1: Sensitivity to the optimal switching scheme) for 500 time steps and an optimization relative tolerance equal to 10 -10 . For problem 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 and 1 as we know that the optimal solution is reached for v ∈ {0, 1} (see Proposition 6.3.1).

0 1 2 3 4 5 -5 -4 -3 -2 -
The numerical results and computation times are summarized in Table 6.2, while Figure 6.3 presents the corresponding trajectories. We note that the direct method give very accurate results, and the computation time for problem P 0 is the lowest because it has only one control. The computation time for problem P 2 is slightly higher than for P 1 because the mixed constraint C m is heavier to evaluate. The numerical method for problem P 3 is of completely different nature as it computes the optimal solution for all the initial conditions on problem max 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.

On Figure 6.3, one may notice some difference between the obtained trajectories. Let us underline that after the peak of y(•), there is no longer uniqueness of the optimal control.

Application to an epidemiological model

The SIR model is one of the most basic transmission model in epidemiology for a directly transmitted infectious disease (for a complete introduction, see for instance [START_REF] Weiss | The SIR model and the foundations of public health[END_REF]) 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 ∈ [0, T ], so that one has S(t) + I(t) + R(t) = 1 with S(t), I(t), R(t) ≥ 0. Let β > 0 be the rate of transmission and γ > 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 ∈ (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:

Ṡ = -(1 -u)β SI (6.16) İ = (1 -u)β SI -γI (6.17) Ṙ = γI (6.18)
When the reproduction number R 0 = β /γ is above one and the initial proportion of susceptible is above the herd immunity threshold S h = R -1 0 , it is well known that there is an epidemic outbreak. Then, the objective is to minimize the peak of the prevalence max t∈[0,T ]

I(t) with respect to control u(•) subject to a L 1 budget T 0 u(t) ≤ Q (6.19)
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 (6.19) were not imposed, then the optimal solution would be the trivial control u(t) = u max , t ∈ [0, T ], which is in general unrealistic from a operational point of view. A similar problem has been considered in [START_REF] Morris | Optimal, near-optimal, and robust epidemic control[END_REF] but under the constraint that intervention occurs only once on a time interval of given length, that we relax here. Note that the constraint (6.19) can be reformulated as a target condition, considering the augmented dynamics

Ṡ = -(1 -u)β SI (6.20) İ = (1 -u)β SI -γI (6.21) Ċ = -u(t) (6.22)
with initial condition C(0) = Q and target {C ⩾ 0}. Extension of the results of Sections 6.3 and 6.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 Chapter 7 as the feedback control

ψ(I, S) := 1 -S h S if I = Ī and S > S h 0 otherwise
where

Ī := I 0 + S 0 -S h -S h log S 0 S h Qβ 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 chapter and can be found in Chapter 7. This control strategy consists in three phases:

1. no intervention until the prevalence I reaches Ī (null control), 2. maintain the prevalence I equal to Ī until S reaches S h (singular control), 3. no longer intervention when S > S h (null control) as illustrated in Figure 6.4 for the parameters given in Table 6.3. Note that differently to the previous The optimal solution for the SIR problem example, this control strategy is intrinsically robust with respect to a bad choice of Ī: the maximum value of I is always guaranteed to be equal to Ī. However, a mischoice of Ī has an impact on the budget (see Chapter 7 for more details).

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 suitable to direct methods. For direct methods that use algebraic differentiation of the dynamics, convergence and β γ T Q S(0) I(0) Ī 0.21 0.07 300 28 1 -10 -6 10 -6 0.105 Table 6.3: Parameters used in numerical computations and optimal value of the peak 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

log e λ ξ + 1 λ -→ λ →+∞ max(ξ , 0), ξ ∈ R
with λ = 100 for the numerical experiments. For problem P θ 3 , one has to be careful about the interplay between the approximations of max(•, 0) and the sequence θ n → +∞, to provide approximations from below of the optimal value. The function h θ is thus approximated by the expression

h θ (x, y, z, u, v) ≃ log e λ 1 g(x,y,u) + 1 λ 1 1 -ve θ λ 2
log(e λ 2 (y-z) +1)

which depends on three parameters λ 1 , λ 2 and θ . Posit for convenience

α := θ λ 2
and consider the function ω α,λ 2 (ξ ) := e -α log(e -λ 2 ξ +1) , ξ ∈ R which approximates the indicator function ⊮ R + . One has the following properties.

Lemma 6.5.1.

1. For any positive numbers α, λ 2 , the function ω α,λ 2 is increasing with

lim ξ →-∞ ω α,λ 2 (ξ ) = 0, lim ξ →+∞ ω α,λ 2 (ξ ) = 1 2. For any ε ∈ (0, 1), one has ω α,λ 2 -ε 2 = ε and ω α,λ 2 (0) = 1 -ε exactly for α = - log(1 -ε) log(2) , λ 2 = log(ε -1 α -1) ε 2 (6.23)
Proof. One has first

ω ′ α,λ 2 (ξ ) = λ 2 α e -λ 2 x e -λ 2 ξ + 1 ω α,λ 2 (ξ ) > 0 and the function ω α,λ 2 (•) is thus increasing. From lim ξ →-∞ -α log(e -λ 2 ξ + 1) = -∞ one get lim ξ →-∞ ω α,λ 2 (ξ ) = 0 and similarly lim ξ →+∞ -α log(e -λ 2 ξ + 1) = 0 implies lim ξ →+∞ ω α,λ 2 (ξ ) = 1
Finally, with simple algebraic manipulation of the conditions ω α,λ 2 -ε 2 = ε and ω α,λ 2 (0) = 1ε, one obtains straightforwardly the expressions (6.23).

We have taken λ 1 = 5000 and considered a sequence of approximations of the indicator function for the values given in Table 6.4 according to expressions (6.23) of Lemma 6.5.1 (see Figure 6 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 -10 ). As one can see in Figure 6.6 and Table 6.5 problems P 0 , P 1 , P 2 give the peak values with a very good accuracy, and present similar performances in terms of computation time. In Figure 6.7 and Table 6.6, the numerical solutions of P θ 3 are illustrated for the values of α and λ 2 given in Table 6.4. As expected, the numerical computation of the family of problems P θ 3 provides an increasing sequence of approximation from below of the optimal value and thus complements the computation of problems P 0 , P 1 or P 2 . From Figures of Tables 6.5 and 6.6, one can safely guarantee that the optimal value belongs to the interval [0.1010, 0.1015]. However, the trajectories found for P θ 3 are not as closed as the ones of problems P 0 , P 1 or P 2 . This can be explained by the fact that problems P θ 3 are not subject to the constraint z(t) ⩾ y(t) and thus provides trajectories for which z(T ) is indeed below max t y(t). is quite far to give a peak value close from the other 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 6.9). We conclude that the L p approximation is not practically reliable for this kind of problems.

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 6 A feedback strategy for peak minimization in SIR model

This chapter corresponds to the accepted pre-print [START_REF] Molina | An optimal feedback control that minimizes the epidemic peak in the SIR model under a budget constraint[END_REF] titled "An optimal feedback control that minimizes the epidemic peak in the SIR model under a budget constraint" to appear in Automatica journal.

Introduction

Since the pioneering work of Kermack and McKendrick [70], the SIR model has been very popular in epidemiology, as the basic model for infectious diseases with direct transmission (see for instance [START_REF] Weiss | The SIR model and the foundations of public health[END_REF][START_REF] Li | An introduction to mathematical modeling of infectious diseases[END_REF] as introductions on the subject). It retakes great importance nowadays due to the recent coronavirus pandemic. In face of a new pathogen, non-pharmaceutical interventions (such as reducing physical distance in the population) are often the first available means to reduce the propagation of the disease, but this has economic and social prices. In [START_REF] Morris | Optimal, near-optimal, and robust epidemic control[END_REF][START_REF] Lobry | Qu'est ce que le pic d'une épidémie et comment le contrôler[END_REF], the authors underline the need of control strategies for epidemic mitigation by "flattering the epidemic curve", rather than eradication of the disease that might be too costly. Several works have applied the optimal control theory considering interventions as a control variable that reduces the effective transmission rate of the SIR model, and studied optimal strategies with criteria based on running and terminal cost over fixed finite interval or infinite horizon [START_REF] Behncke | Optimal control of deterministic epidemics[END_REF][START_REF] Bolzoni | Time-optimal control strategies in SIR epidemic models[END_REF][START_REF] Bolzoni | Optimal control of epidemic size and duration with limited resources[END_REF][START_REF] Kantner | Beyond just "flattening the curve": Optimal control of epidemics with purely non-pharmaceutical interventions[END_REF][START_REF] Palmer | Optimal control of COVID-19 infection rate with social costs[END_REF][START_REF] Bliman | How best can finite-time social distancing reduce epidemic final size?[END_REF][START_REF] Caulkins | The optimal lockdown intensity for COVID-19[END_REF][START_REF] Freddi | Optimal control of the transmission rate in compartmental epidemics[END_REF][START_REF] Ketcheson | Optimal control of an SIR epidemic through finite-time non-pharmaceutical intervention[END_REF][START_REF] Bliman | Optimal immunity control and final size minimization by social distancing for the sir epidemic model[END_REF]. However, the highest peak of the epidemic appears to be the highly relevant criterion to be minimized (especially when there is an hospital pressure to save individuals with severe forms of the infection). In [START_REF] Morris | Optimal, near-optimal, and robust epidemic control[END_REF], the authors studied the minimization of the peak of the infected population under the constraint that interventions occur on a single time interval of given duration. In the present work, we consider the same criterion, but under a budget constraint on the control (as an integral cost) that we believe to be more relevant as it takes into account the strength of the interventions and does not impose an a priori single time interval of given length for the interventions to take place (although we have been able to prove that the optimal solution consists indeed in having interventions on a single time interval but with a control strategy different that the one obtained in [START_REF] Morris | Optimal, near-optimal, and robust epidemic control[END_REF]). Let us also mention a more recent work [START_REF] Avram | Optimal control of a SIR epidemic with ICU constraints and target objectives[END_REF] that considers a kind of "dual" problem, which consists in minimizing an integral cost of the control under the constraint that the epidemic stays below a prescribed value and an additional constraint on the state at a fixed time. The structure of the optimal strategy given by the authors in [START_REF] Avram | Optimal control of a SIR epidemic with ICU constraints and target objectives[END_REF] is similar to the one we obtained without having to fix a time horizon and a terminal constraint. All the cited works rely on numerical methods to provide the effective control. Here, we give an explicit analytical expression of the optimal control. Let us stress that optimal control problems with maximum cost are not in the usual Mayer, Lagrange or Bolza forms of the optimal control theory [START_REF] Cesari | Optimization-theory and applications: problems with ordinary differential equations[END_REF], for which the necessary optimality conditions of Pontryagin's Principle apply, but fall into the class of optimal control with L ∞ criterion, for which characterizations have been proposed in the literature mainly in terms of the value function (see for instance [START_REF] Barron | The Bellman equation for minimizing the maximum cost[END_REF]). Although necessary optimality conditions and numerical procedures have been derived from theses characterizations (see for instance [START_REF] Barron | The Pontryagin maximum principle for minimax problems of optimal control[END_REF][START_REF] Marco | Minimax optimal control problems. Numerical analysis of the finite horizon case[END_REF]), these approaches remain quite difficult and numerically heavy to be applied on concrete problems. On the other hand, for minimal time problems with planar dynamics linear with respect to the control variable, comparison tools based on the application of the Green's Theorem have shown that it is possible to dispense with the use of necessary conditions to prove the optimality of a candidate solution [START_REF] Hermes | Functional Analysis and Time Optimal Control[END_REF]. Although our criterion is of different nature, we show in the present work that it is also possible to implement this approach for our problem.

The chapter is organized as follows. In the next section, we posit the problem of peak minimization to be studied. In Section 7.3, we define a class of feedback strategies that we called "NSN", and give some preliminary properties. Section 7.4 proves the existence of an NSN strategy which is optimal for our problem, and makes it explicit. Finally, Section 7.5 illustrates the optimal solutions on numerical simulations and discusses about the optimal strategy.

Definitions and problem statement

We consider the SIR model

   Ṡ = -β SI(1 -u) İ = β SI(1 -u) -γI Ṙ = γI (7.1)
where S, I and R denotes respectively the proportion of susceptible, infected and recovered individuals in a population of constant size. The parameters β and γ are the transmission and recovery rates of the disease. The control u, which belongs to U := [0, 1], represents the efforts of interventions by reducing the effective transmission rate. For simplicity, we shall drop in the following the R dynamics. Throughout the chapter, we shall assume that the basic reproduction number R 0 is larger than one, so that an epidemic outbreak may occur.

Assumption 3.

R 0 := β γ > 1.
For a positive initial condition (S(0), I(0)) = (S 0 , I 0 ) with S 0 + I 0 ⩽ 1, we consider the optimal control problem which consists in minimizing the epidemic peak under a budget constraint inf

u(•)∈U max t⩾0 I(t), (7.2)
where U denotes the set of measurable functions u(•) that take values in U and satisfy the L 1 constraint

+∞ 0 u(t)dt ⩽ Q . Remark 7.2.1.
From equations (7.1), one can easily check that the solution I(t) tends to zero when t tends to +∞ whatever is the control u(•), so that the supreme of I(•) over [0, +∞) in (7.2) is reached.

2. maintain the prevalence I equal to Ī until S reaches S h (singular control), 3. no longer intervention when S < S h (null control)

Remark 7.3.1. There is no switch of the control between phases 2 and 3, because u(t) tends to zero when S(t) tends to S h , according to expression (7.6).

One can check straightforwardly the following properties are fulfilled.

Lemma 7.3.1. For any Ī ∈ [I 0 , I h ], the maximal value of the control u ψ Ī (•) is given by

u max ( Ī) := 1 - S h S < 1,
where S is solution of S -S h log S = S 0 + I 0 -S h log S 0 -Ī.

Moreover, any solution given by the NSN strategy verifies

max t⩾0 I(t) = Ī.

Optimal strategy

We first show that the function L can be made explicit.

Proposition 7.4.1. One has

L ( Ī) = I h - Ī β S h Ī , Ī ∈ [I 0 , I h ]. (7.7)
Proof. Note first that whatever is Ī, S(•) is decreasing with the control (7.6). One can then equivalently parameterize the solution I(•), C(•) by σ (t) := S 0 -S(t)

instead of t. Posit σ h := σ (t h ) = S 0 -S h .
As long as I < Ī, one has u = 0 which gives

     dI dσ = f (σ ) := 1 -S h S 0 -σ > 0 dC dσ = 0
Remind, from the definition of I h , that the solution I(•) with u = 0 reaches I h in finite time. Therefore, one can define the number σ := inf{σ ⩾ 0,

I(σ ) = Ī} ⩽ σ h , which verifies σ 0 f (σ ) dσ = Ī -I 0 . (7.8) For σ ∈ [ σ , σ h ], one has u = 1 -S h /S, that is      dI dσ = 0 dC dσ = -1 β S h Ī 1 -S h S 0 -σ = -f (σ ) β S h Ī < 0 One then obtains L ( Ī) = C(0) -C(σ h ) = 1 β S h Ī σ h σ ( Ī) f (σ ) dσ
and with (7.8) one can write

L ( Ī) = 1 β S h Ī σ h 0 f (σ ) dσ + I 0 -Ī .
On the other hand, one has

σ h 0 f (σ ) dσ = σ h + S h log S h S 0 = I h -I 0
which finally gives the expression (7.7).

Then, the best admissible NSN control can be given as follows.

Corollary 7.4.1.

When Q ⩽ I h -I 0 β S h I 0
, the smallest Ī ∈ [I 0 , I h ] for which the solution with the NSN strategy is admissible, is given by the value Ī⋆ (Q) := I h Qβ S h + 1 (7.9)

and one has

L ( Ī⋆ (Q)) = Q. ( 7.10)
We give now our main result that shows that the NSN strategy is optimal.

Proposition 7.4.2. Let Assumptions 3 and 4 be fulfilled. Then, the NSN feedback is optimal with

Ī = Ī⋆ (Q), Q < I h -I 0 β S h I 0 I 0 , Q ⩾ I h -I 0 β S h I 0
where Ī⋆ (Q) is defined in (7.9), and Ī is the optimal value of problem (7.2).

Proof.

When Q ⩾ I h -I 0 β S h I 0
, the NSN strategy is admissible and the corresponding solution verifies max t⩾0 I(t) = I 0 which is thus optimal. Consider now

Q < I h -I 0 β S h I 0 . Let (S ⋆ (•), I ⋆ (•),C ⋆ (•)
) be the solution generated by the NSN strategy with Ī = Ī⋆ (Q), and denote u ⋆ (•) the corresponding control. Let S := S ⋆ (t) where t = inf{t > 0, I ⋆ (t) = Ī} and t ⋆ h := inf{t > t, S ⋆ (t) = S h }. We consider in the (S, I) plane the curve

C ⋆ := {(S ⋆ (t), I ⋆ (t)); t ∈ [0,t ⋆ h ]}
For S ⩾ S, the control (7.6) is null and a upward normal to C ⋆ is given by the expression

⃗ n(S, I) = β SI -γI β SI , (S, I) ∈ C ⋆ with S ∈ [ S, S 0 ].
On the other hand, the vector field in the (S, I) plane of any admissible solution is

⃗ v(S, I, u) = -β SI(1 -u) β SI(1 -u) -γI .
Then One has thus I(t h ) < Ī. Define

T := t ⋆ h + 1 γ log Ī I(t h ) > t ⋆ h and consider the non-admissible solution ( S(•), Ĩ(•), C(•)) of (7.3) on [0, T ] defined by the control ũ(t) = u ⋆ (t), t ∈ [0,t ⋆ h ) 1, t ∈ [t ⋆ h , T ]
.

One can straightforwardly check with equations (7.3) that the solution ( S(t), Ĩ(t), C(t)) is

(S ⋆ (t), I ⋆ (t),C ⋆ (t)), t ∈ [0,t ⋆ h ) (S h , Ī exp(-γ(t -t ⋆ h )),C ⋆ (t ⋆ h ) + t ⋆ h -t), t ∈ [t ⋆ h , T ].
Remind, from Corollary 7.4.1, that one has C ⋆ (t ⋆ h ) = 0 by equation (7.10)). Clearly, one has ( S(T ), Ĩ(T )) = (S h , I(t h )) and C(T ) < 0. We consider now in the (S, I) plane the simple closed curve Γ which is the concatenation of the trajectory ( S(•), Ĩ(•)) on forward time with the trajectory (S(•), I(•)) in backward time:

Γ := {( S(τ), Ĩ(τ)), τ ∈ [0, T ]} ∪ {(S(T + t h -t), I(T + t h -t)), τ ∈ [T, T + t h ]} that is anticlockwise oriented by τ ∈ [0, T + t h ]. Then one has C(T ) -C(t h ) = Γ dC.
From equations (7.3), one gets By the Green's Theorem, one obtains

dC = -dS β SI -dt = -dS β SI + dS+dI γI = 1 -S h
C(T ) -C(t h ) = D ∂ Q ∂ S (S, I) - ∂ P ∂ I (S, I) dSdI = D 1 - S h S 1 γI 2 dSdI > 0
where D is the domain bounded by Γ (see Figure 7.1 as an illustration). This implies C(t h ) < C(T ) < 0 and thus a contradiction with the admissibility condition (7.4) of the solution (S(•),

I(•),C(•)). We conclude that (S ⋆ (•), I ⋆ (•),C ⋆ (•)) is optimal.

Numerical illustrations and discussion

We have considered the same parameters and initial condition as in [START_REF] Morris | Optimal, near-optimal, and robust epidemic control[END_REF] (see Table 7

R 0 = 3, S h = 1 3 , I h ≃ 0.3.
Figure 7.2 presents a simulation of the optimal solution for the budget Q = 28, as an example (the minimum peak is reached for Ī ≃ 0.1015). As a comparison, the optimal strategy obtained by Morris et al. in [START_REF] Morris | Optimal, near-optimal, and robust epidemic control[END_REF] for a fixed time duration of interventions without consideration of any budget is quite different (see Figure 7.3). It consists in four phases: no intervention, maintain I constant, apply the maximal control an NSN strategy appears thus less restrictive to be applied in practice. The strategy proposed by Morris et al. induces also a second peak: after the third phase, the prevalence I increases again up to a peak which has to be equal to the level maintained during the second phase if it is optimally chosen. But this second peak turns out to be non robust under a mischoice (or mistiming) of the second phase (see [START_REF] Morris | Optimal, near-optimal, and robust epidemic control[END_REF] for more details). Comparatively, the NSN is naturally robust with respect to a bad choice of Ī: the maximum value of I is always guaranteed to be equal to Ī. However, a mischoice of Ī has an impact on the budget of the NSN strategy, given by expression (7.7) and illustrated in Table 7.2 (for model parameters given in Table 7.1 and Q = 28).

In case of a new epidemic among a large population, one can consider that the initial number of infected individuals is very low, while all the remaining population is susceptible. Therefore, one has S 0 + I 0 = 1 with I 0 very small, and the optimal value of Ī can be well approximated by its limiting Ī -Ī⋆ -10% -5% -1% +5% +10%

L ( Ī) -Q +17% +8% +1.5% -7% -14% From property (7.5), one also gets an approximation of the value Sℓ of S when I reaches Īℓ with u = 0, as the solution of the equation Sℓ + Īℓ -S h log( Sℓ ) = 1, and then an approximation of the duration of the intervention is given by

d ℓ := S h -Sℓ γ Īℓ ,
(one can easily check that along the singular arc I = Ī, one has Ṡ = -γ Ī). For the parameters of The question of parameters estimation in the SIR model from data is out of the scope of the present work. However, while reaching I = Īℓ without intervention, one may expect refinement of the estimates and thus an adjustment of the value of Īℓ . Note that if it is rather the height of the peak Ī that is imposed, the corresponding effort can be determined with expression (7.11), that is

Q = 1 β S h 1 -S h Ī -1 ,
as well with the duration of the intervention.

To have a better insight of the impacts of the available budget Q on the course of the epidemic, we have considered four characteristics numbers:

• t i : the starting date of the intervention,

• d: the duration of the intervention,

• Ī: the height of the peak,

• u max : the maximal value of the control, of the optimal solution, depicted on Figure 7.4 as a function of Q for I 0 = 10 -6 and S 0 + I 0 = 1. Let us note that the maximal budget Q under which it is not possible to immediately slow down the progress of the epidemic is given, according to Proposition 7.4.2, by

Q max := I h -I 0 β S h I 0 ≃ 4.3 10 6 ,
which is quite high. Moreover, the maximal value of the control is bounded by the value far from the value 1 (that would consists in a total lockdown of the population). On Figure 7.4, one can see that the peak Ī can be drastically reduced under a reasonable budget, and that taking larger budgets slows down the decrease of the peak, while the duration of the intervention carries on increasing, almost linearly. Indeed, remind that one has d = ( S -S h )/(γ Ī) and for an optimal value of Ī, one has Q = (I h -Ī)/(γ Ī) from (7.9). Then one gets

d = S -S h I h - Ī Q,
but for large values of Q, Ī is small and S closed to one, which gives an approximation of d as the linear function of

Q d ≃ 1 -S h I h Q ≃ 2.194 Q.
This implies that for a long duration, fixing the budget Q or the duration d tends to be equivalent. Therefore, for a same large duration, the optimal peak gets closed from the optimal one of the strategy of Morris et al. which constraints the duration only, but the difference of the budgets of these two strategies gets increasing with always a lower one for the NSN strategy, as one can see on Figure 7.5.

Finally, this analysis highlights (as already mentioned in [START_REF] Morris | Optimal, near-optimal, and robust epidemic control[END_REF][START_REF] Lobry | Qu'est ce que le pic d'une épidémie et comment le contrôler[END_REF]) the importance to do not intervene too early (unless one has a very large budget) and to choose the "right" time to launch interventions. We believe that curves as in Figure 7.4 might be of some help for decision makers. Chapter 8

Conclusion and perspectives

The objective of this thesis was to apply optimal control techniques in resource management problems. In particular we worked in two contexts, mining and epidemiology. Now, we will present the main results obtained and open problems which we identified.

We begin with the mining context corresponding to Part 1 of this document.

• Conclusions: In Chapter 3, we proposed a formulation based on the continuous framework of [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF], where the only change was to reduce the functional space of profiles from continuous to absolutely continuous functions. Then, we characterized analytically optimal profiles of (FOP) in 2D (Theorem 3.2.1) and 3D (Theorem 3.3.1), both being, main results of that section.

In Chapter 4 we showed (FOP) and (SOP) versions in a 2D case assuming continuity in densities g and e. Then we recovered Theorem 3.2.1 for these hypotheses. Moreover, we proposed a new version of the sequential problem (SOP) that we called Semi continuous formulation consisting in discretizing the space domain Ω and in each point of this discretization take a continuous function characterizing the depth at that point. This formulation is applied on 2D ((SOP) 2D SC ) and 3D cases ((SOP) 3D SC ) and allows us to manipulate the 3D case numerically. We showed numerical experiments for the 2D case, comparing both formulations and methods of resolution (local optimization and HBJ). Besides, we present, for the first time in literature to our knowledge, numerical solutions of (SOP) in the 3D case.

• Perspectives: The hypothesis in Theorem 3.3.1 are a bit strong, therefore, the search for weaker

hypotheses remains an open problem. An option would be work with profiles parameterized in each coordinate by another variable, reducing then the dimension of the optimal control problem space.

In the analytical part, another open problem is to explore optimality conditions for ((SOP) 2D SC ) and ((SOP) 3D

SC ) which could give extra information of optimal profiles. With respect to the numerical part, there does not exists a density function g in literature built from real data. That information is found in blocks as it was presented in the Introduction. A challenge is to design an efficient interpolation tool to pass from a block model to a continuous function g, and then use the same optimization tools used in Chapter 4 to solve an academic example. Now, the epidemiology context corresponding to Part 2 of this document.

• Conclusions: In Chapter 6 we worked with general optimal control problems consisting of minimizing the maximum value of a state which appear in epidemiology to minimize the peak of infected individuals after an outbreak disease. We presented four different equivalent formulations 87

Abstract

Optimal control techniques have numerous applications in engineering and real world problems. This thesis is devoted to using these techniques in two contexts, mining and epidemiology, dividing this document in two respective parts.

In the first part related to mining, we work with the continuous formulation of the Open Pit Problem consisting of finding the optimal shape of an opencast mine representing its profile by a continuous function. Optimality in this context corresponds to maximizing the profit of mineral extraction. We introduce for the first time optimal control models of this problem. We present optimality conditions of solutions along with numerical experiments using local and global methods. Another relevant problem in this context corresponds to the Sequential version of the Open Pit Problem, which consists of scheduling an extraction program over consecutive time frames (for example, a profile each 6 months), finding nested profiles maximizing a discounted profit. We proposed a novel semi-continuous model to obtain solutions of the sequential problem and we use it to present for the first time, to the best of our knowledge, numerical solutions of a three dimensional case (a possible real world mine) including the original Open Pit Problem (case with a single time-frame).

In the second part we deal with optimal control problems minimizing the maximum value of a state. This problematic was inspired by Covid-19, where hospitals and ICU beds were overcrowded due to a high amount of simultaneous infections. We present four different reformulations of this kind of optimal control problem as a Mayer one, each one having its pros and cons. We present the numerical performance of each formulation in an academic example and in a more realistic SIR model where the problem corresponds to minimizing the peak of infectious compartment with integral constraint in the control. With respect to the latter problem, we prove analytically that the structure of the optimal control is null-singular-null and we used it to assess numerical solutions.

Keywords: optimal control, natural system management, mining, epidemiology

Résumé

Les techniques de contrôle optimale ont plusieurs applications dans le domaines de l'engenierie et les problèmes de la vie pratique. Cette thèse est consacrée a l'utilitation de ces techniques dans deux contextes, l'exploitation minière et l'épidémiologie, divisant ce document en deux parties respectives. Dans la première partie relative à l'exploitation minière, nous travaillons avec la formulation continue du problème Final Open Pit consistant à trouver la forme optimale d'une mine à ciel ouvert représentant son profil par une fonction continue. L'optimalité dans ce contexte correspond à maximiser le bénéfice de l'extraction minérale. Nous présentons pour la première fois des modèles de contrôle optimal pour ce problème. Nous présentons des conditions optimales de solutions avec des simulations numériques utilisant des méthodes locales et globales. La version séquentielle du Final Open Pit est également très pertinent dans le contexte de cette thèse. Le problème consiste à planifier un programme d'extraction sur des périodes de temps consécutives (par exemple, un profil tous les 6 mois), en trouvant des profils imbriqués maximisant un bénéfice actualisé. Nous avons formulé un nouveau modèle semi-continu permettant d'obtenir des solutions du problème séquentiel et que nous utilisons pour proposer des solutions numériques dans un cas tridimensionnelle ( une mine possible du monde réel) y compris le problème Open Pit original (cas avec une seule période de temps). Cela, à notre connaissance, n'avait jamais été fait auparavant. Dans la deuxième partie, nous traitons les problèmes de contrôle optimal minimisant la valeur maximale d'un état. Ce problème a été inspiré par le contexte sanitaire très difficile suite a la covid-19, où avoir de beaucoup infectés en même temps sature les hôpitaux et les lits de soin intensif. Nous présentons quatre reformulations différentes de ce type de problème comme un Mayer de contrôle optimal, ayant chacun ses avantages et ses inconvénients. Nous évaluons la performance numérique de chaque formulation dans un exemple académique et dans un modèle SIR plus réaliste où le problème correspond à minimiser le pic du compartiment infecté avec une contrainte intégrale dans le contrôle. En ce qui concerne ce dernier problème, nous prouvons analytiquement que la structure du contrôle optimal est null-singular-null et nous l'avons utilisé pour évaluer les solutions numériques.

Mots clés : contrôle optimal, gestion du système naturel, exploitation minière, épidémiologie Laboratoire Jacques-Louis Lions Sorbonne Université -Campus Pierre et Marie Curie -4 place Jussieu -75005 Paris -France

Figure 2 . 1 :

 21 Figure 2.1: Chuquicamata mine, the biggest open pit mine in the world located in Atacama region of Chile

Figure 2 . 2 :

 22 Figure 2.2: Block model example.

  Generically, we denote a profile of an open pit by p : [a, b] → R + where a and b are the extreme points of the open pit (there is no loss of generality in taking a < b) and where p(x) represents the depth of the profile at the point x ∈ [a, b]; see Figure 3.1.

Figure 3 . 1 :

 31 Figure 3.1: Profile of an open pit on the plane

  [a, b] → R can be expressed via the condition , z)dzdx ⩽ c max . Concerning optimality, the marginal profit at each x ∈ [a, b] of an admissible profile p : [a, b] → R is given by p(x) p 0 (x) g(x, z)dz where g(x, z) represents the profit earned (or gain) for carrying out extraction at the block (x, z) for any z ∈ [p 0 (x), p(x)]. Therefore, the total profit associated with an admissible profile p : [a, b] → R is given by

Theorem 3 . 2 . 1 .Furthermore

 321 Let p ∈ A C [a, b] be an optimal profile of the problem (P 2D ). Then there are ζ ∈ A C [a, b], η ∈ {0, 1}, λ ⩽ 0 and a (positive) Radon measure µ on [a, b], with at least one of them not equal to zero, such that ζ (x) ∈ ηG(x, p(x)) + λ e(x, p(x)) + |µ([a, x[)ζ (x)|∂ p κ(x, p(x)), a.e. on [a, b], with supp(µ) ⊆ {x ∈ [a, b] | p 0 (x) = p(x)}, and where G(x, p) = co g(x, p -), g(x, p + ) , ∀x ∈ [a, b], ∀p ∈ R. z)dzdxc max = 0 and (ζ (x)µ([a, x[))(| ¯p(x)|κ(x, p(x))) = 0 for a.e. x ∈ [a, b].

  -point constraint of an additional state. Let p ∈ A C [a, b] be a given profile. On the one hand, note that for any x ∈ [a, b] such that κ(x, p(x)) ̸ = 0, the condition | ṗ(x)| ⩽ κ(x, p(x)) is equivalent to -1 ⩽ u(x) := ṗ(x) κ(x, p(x)) ⩽ 1. This implies then that the condition | ṗ(x)| ⩽ κ(x, p(x)) is actually equivalent to ṗ(x) = u(x)κ(x, p(x)), with -1 ⩽ u(x) ⩽ 1, for a.e. x ∈ [a, b].

p(x) p 0

 0 (x) e(x, z)dz, for a.e. x ∈ [a, b].

Corollary 3 . 2 . 1 .. 3 )

 3213 Let p ∈ A C [a, b] be an optimal profile of the problem (P 2D ). Suppose that p 0 (x) < p(x), ∀x ∈]a, b[ and | ṗ(x)| < κ(x, p(x)), for a.e. x ∈ [a, b]. (3Then there are η ∈ {0, 1} and λ ⩽ 0, such that 0 ∈ ηG(x, p(x)) + λ e(x, p(x)), a.e. on [a, b]. z)dzdxc max = 0. In particular, if p → g(x, p) is continuous for any x ∈ [a, b] fixed, then the condition reduces to ηg(x, p(x)) + λ e(x, p(x)) = 0, ∀x ∈ [a, b]. Moreover, 1. if the marginal cost associated with extracting a block at any depth is zero (there is no capacity constraint), that is, e(x, z) = 0 for any x ∈ [a, b] and z ⩾ p 0 (x), then the marginal gain of extracting a block at any depth must be zero on the subsection [a, b], that is, g(x, p(x)) = 0, ∀x ∈ [a, b]. 2. if the marginal cost associated with extracting a block at any depth is positive (there is an effective capacity constraint), that is, e(x, z) > 0 for any x ∈ [a, b] and z ⩾ p 0 (x), then η = 1 and g(x, p(x)) + λ e(x, p(x)) = 0, ∀x ∈ [a, b]. Proof. It is enough to apply directly Theorem 3.2.1, and check that ζ (x) = µ([a, x[) for a.e. x ∈ [a, b] and note that µ([a, x[) = 0 for a.e. x ∈ [a, b] because supp(µ) ⊆ {a, b}.

Figure 4 . 1 :

 41 Figure 4.1: 1D / 2D Space discretization of the mine profile

Figure 4 . 2 :

 42 Figure 4.2: Illustration of the 2D profile model discretized w.r.t. space and time as a set of P ki state variables with i = 0, . . . , N profile nodes and k = 0, . . . , T time-frames. Controls U k i are the depths excavated from the previous time-frame at each node. Slopes S k i between neighbor nodes must be smaller than the local maximal slopes i.e. κ(X i , P k i ).

Figure 4 . 3 :

 43 Figure 4.3: Illustration of the 3D profile model with N = 5 and M = 3.View is from 'above', with the state variable P k i, j giving the profile depth at node (x i , y j ) at time-frame k. Slopes S k i, j , T k i, j with neighbors along the x-axis and y-axis must be smaller than the maximal allowed slopes given by the function κ. The control U k i, j (along the z-axis) corresponds to the excavated depth from the same profile node at the previous time-frame P k-1 i, j .

  13) d pc (t) = 0 (4.14) ii) Transversality conditions pP (a), pP (b), pc (a) are free; pc (b) ≥ 0 with pc (b) = 0 if c(b) < c max (4.15) iii) Hamiltonian minimization ū(t) ∈ argmin w H(t, w, ȳ(t), p(t)) a.e. on (a, b) (4.16) iv) State constraint complementary relations d μ(t) ≥ 0, b a (P 0 (t) -P(x))d μ(t) = 0 and μ(b) = 0 (4.17

  )

Remark 4 . 4 . 2 .

 442 The non zero control around x = 200 simply follows the existing initial profile P 0 , and is part of the first constrained arc. See 4.4.2 for more details.

Figure 4 . 4 :

 44 Figure 4.4: 2D profile with limited capacity -global optimization (HJB method)

Figure 4 . 5 :

 45 Figure 4.5: 2D profile with limited capacity -local optimization (direct method) and optimality conditions

Figure 4 . 6 :

 46 Figure 4.6: 2D profile, limited capacity -local optimization using semi-continuous formulation

Figure 4

 4 Figure 4.7: 2D SOP: continuous and semi-continuous formulations.

Figure 4 . 8 :

 48 Figure 4.8: Initial profile for the 3D SOP case

Figure 4 . 9 :Figure 4 . 10 :

 49410 Figure 4.9: 3D profile optimization with limited capacity C = 1e6, for T = 2, 3 and 6 time-frames.

Figure 4 . 11 :

 411 Figure 4.11: 2D profile with infinite capacity -global optimization (HJB method)

Figure 4 . 12 :Figure 4 . 13 :Figure 4 . 14 :

 412413414 Figure 4.12: 2D profile with infinite capacity -local optimization (direct method) and consistency with PMP optimality condition

Figure 5 . 1 :

 51 Figure 5.1: Schema of one fundamental part as a 6-compartment model: susceptible (S), asymptomatic infected (I a ), symptomatic infected (I s ), asymptomatic recovered (R a ), symptomatic recovered (R s ) and dead (D). Red arcs represent the loss of immunity after infection.

Figure 5 .

 5 Figure 5.1 shows the compartments and transitions inside of each fundamental part. The complete model that we will call base structure has 18 compartments (6 for each part) as shown in Figure5.2, with the connections between the three fundamental parts corresponding to vaccination actions.

Proposition 6 . 3 . 1 .

 631 For any control u(•) ∈ U , the optimal control problem inf v∈V z(T ) under the constraint C (6.3) admits an optimal solution. Moreover, an optimal solution verifies z(T ) = max t∈[0,T ]

Remark 6 . 3 . 1 .

 631 The proof of Proposition 6.3.2 gives an optimal construction of z(•) which is the lower envelope of non decreasing continuous functions above the function y(•), as depicted on Figure6.1.

Figure 6 . 1 :

 61 Figure 6.1: Illustration of the function z (in red) corresponding to a function y (in blue) with the control given by expression (6.6)

Assumption 2 .Proposition 6 . 4 . 1 .

 2641 ∀(x, y) ∈ D, G(x, y) := u∈U f (x, y, u) g(x, y, u) is convex. Under Assumption 2, problem P 3 admits an optimal solution. Moreover, any optimal solution Π(•) = (x(•), y(•), z(•)) verifies z(T ) = max t∈[0,T ] y(t)

Figure 6 . 2 :

 62 Figure 6.2: Optimal solution: y(•) on the left, u ⋆ (•) on the right on the switching times of the control are presented in Table6.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 (Bocop software using Gauss II integration

3 Figure 6 . 3 :

 363 Figure 6.3: Comparisons of the three methods on y(•), u(•), z(•) and v(•).

Figure 6 . 4 :

 64 Figure 6.4: The optimal solution for the SIR problem

Figure 6 . 5 :

 65 Figure 6.5: Approximation of the indicator function with different values of ε (zoom on the abscissa axis on the right)

15 Figure 6 . 9 :

 1569 Figure 6.9: Numerical solution for P L 15 without batch iteration (computation time 50s)

  , one has ⃗ n(S, I) •⃗ v(S, I, u) = -β γSI 2 u ⩽ 0, ∀(S, I) ∈ C ⋆ with S ∈ [ S, S 0 ], which shows that any admissible solution is below the curve C ⋆ in the (S, I) plane for S ∈ [ S, S 0 ]. For S ∈ [S h , S], the curve C ⋆ is an horizontal line with I = Ī. Therefore, if there exists an admissible solution (S(•), I(•),C(•)) with max t I(t) < Ī, its trajectory in the (S, I) plane has to be below the curve C ⋆ for any S ∈ [S h , S 0 ]. Let t h := inf{t > 0, S(t) = S h }.

Figure 7 . 1 :

 71 Figure 7.1: The closed curve Γ is composed of the trajectory (S ⋆ (•), I ⋆ (•)) in blue up to to the point (S h , Ī), the additional part ( S(•), Ĩ(•)) in red and the hypothetical better trajectory (S(•), I(•)) in backward time in green.

Figure 7 . 2 :Figure 7 . 3 :

 7273 Figure 7.2: Optimal solution for Q = 28.

Figure 7 . 4 :

 74 Figure 7.4: Characteristics numbers as functions of Q.

Figure 7 . 5 :

 75 Figure 7.5: Comparison of the performances of the optimal strategies with same duration (in abscissa).

1 B. 2 . 2

 122 Problem functions // +++ DRAFT +++ This class implements the OCP functions // It derives from the generic class bocop3OCPBase // OCP functions are defined with templates since they will be called // from both the NLP solver ( double arguments ) and AD tool ( ad_double arguments ) # include < OCP .h > // gloabl N , DX , T ? file generation for mine problem ( ' PDE formulation 'dim_state = (N -1) *( M -1) + 1 # interior grid points for profile + gain dim_control = (N -1) *(M -1) # digging effort for interior grid points dim_boundaryconditions = (N -1) *(M -1) + 2 # initial interior profile + initial gain + final gain dim_pathconstraints = 2* M*N + 1 # X /Y slopes for points except farther boundaries + capacity limit # later use default options ? with open ( filename , 'w ') as deffile : deffile . write ( '# Definition file \ n\n ') deffile . write ( '# Dimensions \n ')

  1,1 ([a, b]; R n ). Under the conditions stated above, if q ∈ A C n [a, b] is a weak local minimizer of (P M ) related to the optimal control ū, then there exist λ ∈ A C n [a, b], η ∈ {0, 1}, a (positive) Radon measure µ on [a, b], and a Borel measurable function γ

	It turns out that, in this setting, weak local minimizers of (P M ) satisfy Maximum Principle for State
	Constrained problems (Vinter [97, Theorem 9.3.1]).
	Lemma 3.2.2.

  0 and β 2 = 0 .

	By point 3 in Lemma 3.2.2, we have that λ2 = η ∈ {0, 1} and λ1 ⩽ 0 with
	λ1	a	b	p(x) p 0 (x)	e(x, z)dzdx -c max = 0.
	By point 4 in Lemma 3.2.2 we have, since κ is nonnegative, that
	ξ 3 (x) ū(x)κ(x, p(x)) = |ξ 3 (x)|κ(x, p(x)), for a.e. x ∈ [a, b].

  .1 are the settings for the different numerical methods used in the simulations.We start with the 2D FOP as first example, since it is the only one for which all formulations, including global optimization, are available. We set a maximal capacity c max = 20.000. The gain function g is interpolated from values found in the Marvin block model of MINELIB, a publicly available library of test problem instances for open pit mining problems (see[START_REF] Espinoza | MineLib: a library of open pit mining problems[END_REF]). Solutions for the unlimited capacity case, with a different control structure (i.e. singular arcs), are shown in 4.7, with both constant and variable maximal slope.

	Remark 4.4.1.

2D FOP (global) t : 123 steps; P : 50 steps, e : 210 steps; u : 100 steps 2D FOP/SOP (local) t : 123 steps; tol = 10 -10 ; maxiter = 10000 2D SOP SC (local) t : T steps; N = 123 nodes; tol = 10 -10 ;maxiter = 10000 3D SOP SC (local) t : T steps; 30 × 10 nodes; tol = 10 -6 ;maxiter = 10000 Table 4.1: Numerical settings for the continuous and semi continuous formulations 4.4.2 Final open pit (2D): global and local optimization

Table 4 .

 4 2: Solutions for the 2D FOP with limited capacity.

		Objective CPU
	2D FOP Global optim.	10846 369s
	2D FOP Local optim.	11093	3s
	2D FOP SC Local optim.	11100	2s

Table 4 . 4

 44

	6

6

Capacity limit : 5 • 10 : 3D SOP: solutions from the semi-continuous formulation, for different time intervals and capacity limit per time-frame.

Table 6 .

 6 2: Comparison of the numerical results

		t∈[0,T ]	y(t)	error	computation time
	P P 0	2.24705 2.249888 0.126% 0	-0.5 s
	P 1	2.24998	0.130%	1.8 s
	P 2	2.249941 0.129%	3.8 s
	P 3	2.26778	0.8%	248 s

Table 6 .

 6 .5). 4: Values of parameters α, λ 2 for different ε

	ε	α	λ 2
	0.2	0.32	124
	0.15	0.234 360
	0.1	0.152 1514
	0.075 0.112 4094
	0.05	0.074 16193

Table 6 .

 6 Figure 6.6: Comparisons of numerical results for the methods P 0 , P 1 , P 2 5: Comparison of performances for problems P 0 , P 1 , P 2 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 ∈ {2, 5, 10, 15} (see Figure6.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

	0.08 0.10						y 0 y 1 y 2		0.5 0.6						u 0 u 1 u 2
	0.06								0.4
									0.3
	0.04								0.2
	0.02								0.1
	0.00								0.0
	0	50	100	150	200	250	300		0	50	100	150	200	250	300
	0.10								1.0
	0.08								0.9
	0.06								0.8
	0.04								0.7
	0.00 0.02						z 0 z 1 z 2		0.6 0.5						v 1 v 2
	0	50	100	150	200	250	300		0	50	100	150	200	250	300
				problem max t∈[0,T ]	y(t) computation time
					P 0		0.1015		10 s
					P 1		0.1015		12 s
					P 2		0.1015		13 s
				ε		z(T )	max t∈[0,T ]	y(t) computation time
				0.2	0.0684	0.1038	80 s
				0.15 0.0823	0.1038	65 s
				0.1	0.0954	0.1037	51 s
				0.075 0.0993	0.1050	83 s
				0.05 0.1010	0.1036	97 s
			Table 6.6: Comparison of performances for problem P θ 3
	Finally, we have compared our approximation technique with the classical approximation of the L ∞

Table 6 .

 6 7: Comparison of the numerical results with the L p approximation

  .8, that could help practitioners in the choice of the method.

	Chapter 7
	Formulation	P 0 P 1 or P 2 P 3 P θ 3
	suitable to direct methods	yes	yes	no	yes
	suitable to Hamilton-Jacobi-Bellman methods	no	yes	yes yes
	suitable to shooting methods without constraint no	no	no	yes
	provides approximations from below	no	no	no	yes
	Table 6.8: Comparison of the different formulations

  .1). For these values,

	β	γ	S(0)	I(0)
	0.21 0.07 1 -10 -6 10 -6

Table 7

 7

	.1: SIR parameters and initial condition
	one computes

Table 7 .

 7 2: Variation of the control budget of the NSN strategy under a mischoice of Ī

	expression for I 0 = 0, that is	Īℓ :=	1 -S h + S h log(S h ) Qβ S h + 1	.	(7.11)

Table 7 .

 7 Table 7.1, one obtains the limiting values given in Table7.[START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF]. This means that depending on the budget Q only, one can determine the minimal peak and the optimal strategy to apply, without the knowledge of the initial size of the infected population, provided that parameters β and γ of the disease are known.

	Īℓ	Sℓ	d ℓ
	0.1015 0.8406 71.39

3:

The limiting optimal values for arbitrarily small I 0 (with Q = 28)

  template instanciation for template functions , with double and double_ad // +++ could be in an included separate file ? // but needs to be done for aux functions too ? APPARENTLY NOT ! template void OCP :: finalCost < double >( double initial_time , double final_time , const double * initial_state , const double * final_state , const double * parameters , const double * constants , double & final_cost ); template void OCP :: dynamics < double >( double time , const double * state , const double * control , const double * parameters , const double * constants , double * state_dynamics ); template void OCP :: boundaryConditions < double >( double initial_time , double final_time , const double * initial_state , const double * final_state , const double * parameters , const double * constants , double * boundary_conditions ); template void OCP :: pathConstraints < double >( double time , const double * state , const double * control , const double * parameters , const double * constants , double * path_constraints ); template void OCP :: finalCost < double_ad >( double initial_time , double final_time , const double_ad * initial_state , const double_ad * final_state , const double_ad * parameters , const double * constants , double_ad & final_cost ); template void OCP :: dynamics < double_ad >( double time , const double_ad * state , const double_ad * control , const double_ad * parameters , const double * constants , double_ad * state_dynamics ); template void OCP :: boundaryConditions < double_ad >( double initial_time , double final_time , const double_ad * initial_state , const double_ad * final_state , const double_ad * parameters , const double * constants , double_ad * boundary_conditions ); template void OCP :: pathConstraints < double_ad >( double time , const double_ad * state , const double_ad * control , const double_ad * parameters , const double * constants , double_ad * path_constraints );

		control . set uniform bcp :: readFileToVector ( " data / depth . data " , depth_grid );
		# uniform case bcp :: readCSVToMatrix (" data / gain . csv " , G_values , '; ' ,0 ,1) ;
		control .0. lowerbound -1 bcp :: readCSVToMatrix (" data / slope . csv " , K_values , '; ' ,0 ,1) ;
	}	control .0. upperbound 1
		control .0. steps 100
	# System modes ( >=1) and admissible transitions for switching between modes // ///////////////////////////////////////////////////////////////////
	# all_transitions ; user_file ; user_function # state jumps at transitions : none ; user_function system . modes 1 Appendix B // explicit
		admissible . transitions all_transitions
		state . jumps none
	# Value function save files : # Previous steps type : resume ; overwrite ; ask Codes local optimization (Bocop)
		# Output format type : text ; binary ; none
		# Output path : . --> here ; other repertoire
		valueFunction . previous . steps overwrite
		valueFunction . output . format text
	valueFunction . output . path valueFunction / B.1 1D SOP -continuous formulation
		# SimulatedTrajectory :
	# Computation : none ; after_valueFunction ; read_valueFunction # Noise : none ; gaussian ; user_function B.1.1 Problem definition
		# Starting mode : best_mode ; user_function ; value :[0 , nbmode -1]
	# Output path : . --> here ; other repertoire simulatedTrajectory . computation after_valueFunction # Definition file
	simulatedTrajectory . output . path trajectory / # Dimensions simulatedTrajectory . noise gaussian dim . state 3 simulatedTrajectory . starting . mode 0 dim . control 1 simulatedTrajectory . starting . state .0 0.001 dim . boundaryconditions 6 simulatedTrajectory . starting . state .1 0.1 dim . pathconstraints 1 simulatedTrajectory . other_output 0 dim . parameters 0
	# ProcessLaw : dim . constants 2
	# Computation : true ; false # Time interval # Output Path : . --> here ; other repertoire initial . time 0 processLaw . computation false processLaw . output . path processLaw / final . time 1230
	processLaw . initial . path processLaw / initialDistribution / # Constants
	constant .0 1 # slope type ( value or '0 ' for table ) HD . option false constant .1 0 # quadratic regularization . effect not clear ...
	# Time discretisation NB ; singular case , CV not necessarily better when increasing A.1.2 Objective steps from 1000 to 2000
	time . steps 1000
		// This code is published under the Eclipse Public License
	// Authors : Daphne Giorgi , Benjamin Heymann , Jinyan Liu , Pierre Martinon , Olivier # Bounds for constraints
	Tissot boundarycond .0. lowerbound 0
	// Inria Saclay and Cmap Ecole Polytechnique boundarycond .0. upperbound 0
	// 2014 -2017 boundarycond .1. lowerbound 0
	boundarycond .1. upperbound 0
	boundarycond .2. lowerbound 0
	// Function for the running cost boundarycond .2. upperbound 0
	// Input : boundarycond .3. lowerbound 0
	// time : current time t boundarycond .3. upperbound 0
	// initial_time : t0 boundarycond .4. lowerbound 0
	// final_time : tf boundarycond .5. upperbound 2 e4
	// state : vector of state variables x pathconstraint .0. lowerbound 0
		// control : vector of control variables u
	// mode : mode of the system i # Bounds for variables
	// constants : vector of constants state .0. upperbound 500
		97

  / +++ DRAFT +++ This class implements the OCP functions // It derives from the generic class bocop3OCPBase // OCP functions are defined with templates since they will be called // from both the NLP solver ( double arguments ) and AD tool ( ad_double arguments ) // # pragma once # include < OCP .h > // /////////////////////////////////////////////////////////////////// // space domain is discretized on a X/Y grid with N x M intervals // this means ( N +1) ( M +1) points on the grid including the boundaries // since profile is equal to 0 on the boundary , only the (N -1) (M -1) interior points are modeled as state variables // this matrix is stored by columns ie 1D index is l = j * (N -1) + iAPPLICATION OF OPTIMAL CONTROL TECHNIQUES TO NATURAL SYSTEMS MANAGEMENT

	deffile . write ( ' boundarycond . '+ str (i )+ '. upperbound 0\ n ')
	# final gain
	deffile . write ( ' boundarycond . '+ str ( dim_boundaryconditions -1) + '. lowerbound '+ str
	( objective_lowerbound )+ '\ n ')
	# profile slope constraints
	for i in range ( dim_pathconstraints -1) :
	deffile . write ( ' pathconstraint . '+ str (i)+ '. lowerbound -'+ str ( snorm )+ '\n ')
	deffile . write ( ' pathconstraint . '+ str (i)+ '. upperbound '+ str ( snorm )+ '\ n ')
	# capacity constraint
	deffile . write ( ' pathconstraint . '+ str ( dim_pathconstraints -1) + '. lowerbound 0\ n ')
	deffile . write ( '\n# Bounds for variables \n ')
	# profile bounds
	for i in range ( dim_state -1) :
	deffile . write ( ' state . '+ str (i)+ '. lowerbound 0\ n ')
	deffile . write ( ' state . '+ str (i)+ '. upperbound 500\ n ')
	# gain
	deffile . write ( ' state . '+ str ( dim_state -1) + '. lowerbound 0\ n ')
	# digging bounds
	for i in range ( dim_control ):
	deffile . write ( '\n# Constants \ n ')
	deffile . write ( ' constant .0 ' + str ( N)+ '\n ')
	deffile . write ( ' constant .1 ' + str ( M)+ '\n ')
	deffile . write ( ' constant .2 ' + str ( Cmax_per_time_unit * T / steps )+ '\n ')
	deffile . write ( ' constant .3 ' + str ( alpha ) + '\n ')
	deffile . write ( '\n# Time discretisation \ n ')
	deffile . write ( ' ode . discretization euler_implicit \n ')
	deffile . write ( ' time . steps ' + str ( steps ) + '\n ')
	deffile . write ( '\n# Bounds for constraints \n ') // constants : [N , M , Cmax_per_time_unit , alpha ]
	# initial profile and gain
	for i in range ( dim_boundaryconditions -1) : double X_length = 1230;
	deffile . write ( ' boundarycond . '+ str (i )+ '. lowerbound 0\ n ') double Y_length = 400;

deffile . write ( ' dim . state ' + str ( dim_state ) + '\n ') deffile . write ( ' dim . control ' + str ( dim_control ) + '\n ') deffile . write ( ' dim . boundaryconditions ' + str ( dim_boundaryconditions ) + '\n ') deffile . write ( ' dim . pathconstraints ' + str ( dim_pathconstraints ) + '\n ') deffile . write ( ' dim . parameters 0\ n ') deffile . write ( ' dim . constants 4\ n ') deffile . write ( '\n# Time interval \ n ') deffile . write ( ' initial . time 0\ n ') deffile . write ( ' final . time ' + str ( T)+ '\n ') deffile . write ( ' control . '+ str (i )+ '. lowerbound 0\ n ') deffile . write ( ' control . '+ str (i )+ '. upperbound 500\ n ') deffile . write ( '\n# Initialization for discretized problem \n ') for i in range ( dim_state ) : deffile . write ( ' state . '+ str (i)+ '. init 10\ n ') for i in range ( dim_control ): deffile . write ( ' control . '+ str (i )+ '. init 1\ n ') deffile . write ( '\n# Names \n\n# Ipopt \n ') deffile . write ( ' ipoptIntOption . print_level 5\ n ') deffile . write ( ' ipoptIntOption . max_iter '+ str ( maxiter )+ '\ n ') deffile . write ( ' ipoptStrOption . mu_strategy monotone \ n ') deffile . write ( ' ipoptNumOption . tol '+ str ( tol )+ '\n ') deffile . write ( '\n# Misc \n ') deffile . write ( ' ad . retape 1\ n ') B.3.2 Problem functions /

Explication obtained from Bocop user guide in https://www.bocop.org/download/

A set-valued map F : X ⇝ X is upper semi-continuous at ξ ∈ X if and only if for any neighborhood N of F(ξ ), there exists η > 0 such that for any ξ ′ ∈ B X (ξ , η) one has F(ξ ′ ) ⊂ N (see for instance[START_REF] Aubin | Viability Theory[END_REF]).

Acknowledgements

Pierre and me, helping to start the cotutelle. No quiero dejar de agradecer en esta parte

Notation Description

Portion of population which is infected, symptomatic and is completely isolated

Factor in transmission rate due to non pharmaceutical measures.

Effectiveness of vaccines in transmission, i ∈ {N,V,V R } γ k i Recovery rates of infected individuals, i ∈ {N,V,V r } and k ∈ {a, s} δ i Death rates, i ∈ {N,V,V r }. q k i Rate of loss of immunity for a recovered individual, i ∈ {N,V,V r } and k ∈ {a, s}. φ i Probability of being asymptomatic, i ∈ {N,V,V r } Table 5.1: Parameters used in the model.

Fundamental part (N), non vaccinated:

; Fundamental part (V ), vaccinated without booster:

t); Fundamental part (V r ), vaccinated with booster:

Equivalently, one can consider the extended dynamics.

with the initial condition (S(0), I(0),C(0)) = (S 0 , I 0 , Q) and the state constraint C(t) ⩾ 0, t ⩾ 0. (7.4)

A solution of (7.3) is admissible if the control u(•) takes its values in U and the condition (7.4) is fulfilled.

The NSN feedback

Let us denote the immunity threshold

Note that S(•) is a non increasing function and that one has İ ⩽ 0 when S ⩽ S h , whatever is the control. If S 0 ⩽ S h , the maximum of I(•) is thus equal to I 0 for any control u(•), which solves the optimal control problem. We shall now consider that the non-trivial case.

Assumption 4. S 0 > S h .

Under this assumption, we thus know that for any admissible solution, the maximum of I(•) is reached for S ⩾ S h . For the control u = 0, one can easily check that the following property is fulfilled S(t) + I(t) -S h log(S(t)) = S 0 + I 0 -S h log(S 0 ), t > 0, (7. [START_REF] Anderson | Discussion: the Kermack-McKendrick epidemic threshold theorem[END_REF] and the maximum of I(•) is then reached for the value I h := I 0 + S 0 -S h -S h log S 0 S h .

We define the "NSN" (for null-singular-null) strategy as follows. We denote the L 1 norm associated to the NSN control

where u ψ Ī (•) is the control generated by the feedback (7.6).

This control strategy consists in three phases:

1. no intervention until the prevalence I reaches Ī (null control), CHAPTER 8. Conclusion and perspectives as a Mayer optimal control problem with and without state constraints. These formulations let us apply classical optimal control tools.

We showed the numerical experiment in an academic example and in a more realistic SIR problem. We compared performance among each formulation and L p approximation, used to approximate L ∞ norm. The main advantages and disadvantages in a numerical sense are summarized in the following table. In Chapter 7 we worked the particular problem of minimizing the peak of infected individuals over a SIR model with a L 1 cost in the control. We found the analytical expression for the optimal strategy which we called NSN (null-singular-null) and which corresponds to the feedback control (6.5.2).

Finally we presented numerical experiment and a parametric analysis of solutions comparing the budget Q with the duration of the intervention and height of the optimal peak.

• Perspectives: Chapter 6 was focused on modeling and numerical experiments, therefore, there is still work to be done in finding optimality conditions using one of the four formulations.

The strategy NSN presented in Chapter 7 seems possible to extend to more general planar optimal control problems and it is a work in progress. Besides it would be interesting investigate the performance of this strategy with other related problems, for example, minimizing the final value of susceptible S(T ) as in [START_REF] Bliman | Optimal immunity control and final size minimization by social distancing for the sir epidemic model[END_REF].

The compartment model introduced in section 5.2 has not been exploited in an optimal control sense which would be the following step to take. For example, it would be interesting to study what happens with results of Chapters 6 and 7 using this new model rather than the SIR one or indeed consider the vaccination rate as a control to be optimize.

Appendix A

Codes global optimization (BocopHJB)

A. return noise ; } void other_outputs ( double time , vector < double > state , vector < double > control , vector < double > noise , vector < double > constants , vector < double >& output ) { } # include " dependencies . hpp " vector < double > distance_grid ; vector < double > depth_grid ; vector < double > p0_values ; vector < vector < double > > G_values , K_values ; # include < cmath > # include " publicTools . // full grid is (N +1) (M +1) points with i|j indices from 0 to N|M // interior is (N -1) (M -1) points with i| j indices from 1 to N -1| M -1 // summation for gain is done for indices i =0.. N -1 and j =0.. M -1 since we take [i , i +1] x [j ,j +1] cells // +++ use something more centered here , with loop on all interior points exactly ? // for ( int i = -1; i <=N -2; i ++) // for ( int j = -1; j <=M -2; j ++) for ( int i =0; i <=N -1; i ++) for ( int j =0; j <=M -1; j ++) int M = ( int ) constants [START_REF] Aguiar | The role of mild and asymptomatic infections on COVID-19 vaccines performance: a modeling study[END_REF]; // 2D (N -1) x(M -1) grid stored by columns ie 1D index is l = j * (N -1) + i // digging progress : p_ij ^k +1 -p_ij ^k = u_ij ^k ( not for boundaries which are fixed to 0) for ( int l =0; l <(N -1) *(M -1) ; l ++) state_dynamics [l] = control [l ]; // gain at timeframe TIME STEP SHOULD BE EQUAL TO 1 ( or gain will be scaled wrongly ) // NB . if using implicit euler the state will correctly be the profile at the end of the time step ( NOT if using midpoint !) double alpha = constants [START_REF] Alvarez | A continuous framework for open pit mine planning[END_REF];