Je Voudrais

Pedro Gajardo

Joseph Frédéric Bonnans

Keywords: Optimal control, minimum time, stochastic control, bioremediation, chemostat. v vi

This thesis consists of two parts. In the first part we study minimal time strategies for the treatment of pollution in large water volumes, such as lakes or natural reservoirs, using a single continuous bioreactor that operates in a quasi-steady state. The control consists of feeding the bioreactor from the resource, with clean output returning to the resource with the same flow rate. We drop the hypothesis of homogeneity of the pollutant concentration in the water resource by proposing three spatially structured models. The first model considers two zones connected to each other by diffusion and only one of them treated by the bioreactor. With the help of the Pontryagin Maximum Principle, we show that the optimal state feedback depends only on the measurements of pollution in the treated zone, with no influence of volume, diffusion parameter, or pollutant concentration in the untreated zone. We show that the effect of a recirculation pump that helps to mix the two zones is beneficial if operated at full speed. We prove that the family of minimal time functions depending on the diffusion parameter is decreasing. The second model consists of two zones connected to each other by diffusion and each of them connected to the bioreactor. This is a problem with a non convex velocity set for which it is not possible to directly prove the existence of its solutions. We overcome this difficulty and fully solve the studied problem applying Pontryagin's principle to the associated problem with relaxed controls, obtaining a feedback control that treats the most polluted zone up to the homogenization of the two concentrations. We also obtain explicit bounds on its value function via Hamilton-Jacobi-Bellman techniques. We prove that the minimal time function is nonmonotone as a function of the diffusion parameter. The third model consists of a system of two zones connected to the bioreactor in series, and a recirculation pump between them. The control set depends on the state variable; we show that this constraint is active from some time up to the final time. We show that the optimal control consists of waiting up to a time from which it is optimal the mixing at maximum speed, and then to repollute the second zone with the concentration of the first zone. This is a non intuitive result. Numerical simulations illustrate the theoretical results, and the obtained optimal strategies are tested in hydrodynamic models, showing to be good approximations of the solution of the inhomogeneous problem. The second part consists of the development and study of a stochastic model of sequencing batch reactor. We obtain the model as a limit of birth and death processes. We establish the existence and uniqueness of solutions of the controlled equation that does not satisfy the usual assumptions. We prove that with any control law the probability of extinction is positive, which is a non classical result. We study the problem of the maximization of the probability of attaining a target pollution level, with the reactor at maximum capacity, prior to extinction. This problem does not satisfy any of the usual assumptions (non Lipschitz dynamics, degenerate locally Hölder diffusion parameter, restricted state space, intersecting reach and avoid sets), so the problem must be studied in two stages: first, we prove the continuity of the uncontrolled cost function for initial conditions with maximum volume, and then we develop a dynamic programming principle for a modification of the problem as an optimal control problem with final cost and without state constraint.

RÉSUMÉ

Cette thèse se compose de deux parties. Dans la première partie, nous étudions les stratégies de temps minimum pour le traitement de la pollution dans de grandes ressources en eau, par exemple des lacs ou réservoirs naturels, à l'aide d'un bioréacteur continu qui fonctionne à un état quasi stationnaire. On contrôle le débit d'entrée d'eau au bioréacteur, dont la sortie revient à la ressource avec le même débit. Nous disposons de l'hypothèse d'homogénéité de la concentration de polluant dans la ressource en proposant trois modèles spatialement structurés. Le premier modèle considère deux zones connectées l'une à l'autre par diffusion et seulement une d'entre elles connectée au bioréacteur. Avec l'aide du Principe du Maximum de Pontryagin, nous montrons que le contrôle optimal en boucle fermée dépend seulement des mesures de pollution dans la zone traitée, sans influence des paramètres de volume, diffusion, ou la concentration dans la zone non traitée. Nous montrons que l'effet d'une pompe de recirculation qui aide à homogénéiser les deux zones est avantageux si opérée à vitesse maximale. Nous prouvons que la famille de fonctions de temps minimal en fonction du paramètre de diffusion est décroissante. Le deuxième modèle consiste en deux zones connectées l'une à l'autre par diffusion et les deux connectées au bioréacteur. Ceci est un problème dont l'ensemble des vitesses est non convexe, pour lequel il n'est pas possible de prouver directement l'existence des solutions. Nous surmontons cette difficulté et résolvons entièrement le problème étudié en appliquant le principe de Pontryagin au problème de contrôle relaxé associé, obtenant un contrôle en boucle fermée qui traite la zone la plus polluée jusqu'au l'homogénéisation des deux concentrations. Nous obtenons des limites explicites sur la fonction valeur via des techniques de Hamilton-Jacobi-Bellman. Nous prouvons que la fonction de temps minimal est non monotone par rapport au paramètre de diffusion. Le troisième modèle consiste en deux zones connectées au bioréacteur en série et une pompe de recirculation entre elles. L'ensemble des contrôles dépend de l'état, et nous montrons que la contrainte est active à partir d'un temps jusqu'à la fin du processus. Nous montrons que le contrôle optimal consiste à l'atteinte d'un temps à partir duquel il est optimal de recirculer à vitesse maximale et ensuite ré-polluer la deuxième zone avec la concentration de la première. Ce résultat est non intuitif. Les stratégies optimales obtenues sont testées sur des modèles hydrodynamiques, en montrant qu'elles sont de bonnes approximations de la solution du problème inhomogène. La deuxième partie consiste au développement et l'étude d'un modèle stochastique de réacteur biologique séquentiel. Le modèle est obtenu comme une limite des processus de naissance et de mort. Nous établissons l'existence et l'unicité des solutions de l'équation contrôlée qui ne satisfait pas les hypothèses habituelles. Nous prouvons que pour n'importe quelle loi de contrôle la probabilité d'extinction de la biomasse est positive. Nous étudions le problème de la maximisation de la probabilité d'atteindre un niveau de pollution cible, avec le réacteur à sa capacité maximale, avant l'extinction. Ce problème ne satisfait aucune des suppositions habituelles, donc le problème doit être étudié dans deux étapes: en premier lieu, nous prouvons la continuité de la fonction de coût non contrôlée pour les conditions initiales avec le volume maximal et ensuite nous développons un principe de programmation dynamique pour une modification du problème original comme un problème de contrôle optimal avec coût final sans contrainte sur l'état.

Mots clés: Contrôle optimal, temps minimal, contrôle stochastique, biorestauration, chemostat.

iii RESUMEN La tesis de compone de dos partes. En la primera parte estudiamos estrategias de tiempo mínimo para el tratamiento de la contaminación en recursos acuíferos de gran volumen, tales como lagos o reservas naturales, mediante el uso de un biorreactor continuo que opera en un estado cuasi-estacionario. El control consiste en la alimentación del biorreactor desde el recurso, con un efluente más limpio siendo devuelto al recurso con el mismo caudal. Eliminamos la hipótesis de homogeneidad en la concentración del contaminante en el recurso proponiendo tres modelos espacialmente estructurados. El primer modelo considera dos zonas conectadas entre ellas mediante difusión y donde sólo una de ellas es tratada por el biorreactor. Con la ayuda el Principio del Máximo de Pontryagin probamos que el control retroalimentado óptimo depende sólo de las mediciones del contaminante en la zona tratada, sin dependencia del volumen, de la difusión, o de la concentración del contaminante en la zona no tratada. Mostramos que el efecto de añadir una bomba de recirculación que ayuda a mezclar ambas zonas es benéfico si ésta se opera a su máxima velocidad. El segundo modelo consiste en dos zonas conectadas entre sí por difusión y cada una de ellas conectada al biorreactor. Este es un problema donde el conjunto de velocidades es no convexo y para el cual no es posible probar directamente la existencia de soluciones. Superamos esta dificultad y resolvemos completamente el problema estudiad aplicando el principio de Pontryagin al problema asociado con controles relajados, obteniendo un control retroalimentado que trata le zona más contaminada hasta la homogeneización de ambas zonas. También obtenemos cotas explícitas sobre la función valor mediante técnicas de Hamilton-Jacobi-Bellman. Probamos que la función de tiempo mínimo es no-monótona como función del parámetro de difusión. El tercer modelo consiste en un sistema de dos zonas conectadas al biorreactor en serie, y una bomba de recirculación entre ellas. El conjunto de controles depende de la variable de estado; mostramos que esta restricción es activa a partir de cierto instante de tiempo hasta el final del proceso. Este es un resultado no intuitivo. Simulaciones numéricas ilustran los resultados teóricos, y las estrategias obtenidas son testeadas en modelos hidrodinámicos, mostrando ser buenas aproximaciones de la solución del problema no homogéneo. La segunda parte consiste en el desarrollo y estudio de un modelo estocástico de biorreactor secuencial por lotes. Obtenemos el modelo como un límite de procesos de nacimiento y muerte. Establecemos la existencia y unicidad de soluciones de la ecuación controlada que no satisface las hipótesis usuales. Probamos que para cualquier control, la probabilidad de extinción es positiva, resultado que no es clásico. Estudiamos el problema de la maximización de la probabilidad de llegar al nivel deseado de contaminación, con el reactor lleno, antes de la extinción. Este problema no satisface ninguna de las hipótesis usuales (dinámica no Lipschitz, coeficiente de difusión degenerado localmente Hölder, restricciones de espacio de estado, conjuntos objetivo y absorbente se intersectan), por lo que el problema debe ser estudiado en dos etapas: primero, probamos la continuidad de la función de costo sin control para condiciones iniciales con volumen máximo, y luego desarrollamos un principio de programación dinámica para una modificación del problema como un problema de control óptimo con costo final y sin restricciones de estado.

Palabras claves: Control óptimo, tiempo mínimo, control estocástico, bioremediacion, quimiostato. In green, the washout equilibrium; in red, the nontrivial equilibrium. . . . . . . . 1.6 Equilibrium points for the Monod growth function. We see that D 1 < µ max , and then the concentration of substrate at equilibrium is s † r = λ 1 . For D 2 , the equilibrium is the washout. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Equilibrium points for the Haldane growth function. We see that for D 1 < µ(s † ) there exist two equilibrium points λ 1 (locally asymptotically stable) and λ 2 (unstable) besides the washout (unstable). For D 2 = µ(s † ) the equilibrium concentration is s ⋆ r = s † (unstable), and the washout is unstable. For D 3 , the equilibrium is the washout (stable). . . . . . . . . . . . . . . . . . . . . . . . . 1.8 Relation between fictitious time τ and real time t . . . . . . . . . . . . . . . . 1.9 Behavior of the inhomogeneous representation of a lake. From the homogeneous initial distribution of pollutant (on the left), the system evolves up to a point in which two zones are clearly differentiated (on the right). Taken from [

Bioprocesses and bioremediation

Bioremediation is the process that uses living organisms (usually microorganisms) or microbial processes to produce molecular transformations or degradations on environmental contaminants (hazardous to soils, groundwater, sediments, surface water, or air) into products of less toxic form. When these changes occur naturally without human intervention, the process is called natural attenuation. Nevertheless, the speed of such changes is slow. By means of appropriate control techniques these biological systems can be used to enhance the speed of the changes or degradations, as well as to use them in places with high pollutant concentrations.

Microorganisms have the ability to biodegrade most of the organic contaminants and many inorganic contaminats, for instance, hydrocarbon, pesticides, herbicides, petroleum, gasoil, heavy metals among others. Biological treatments of organic contaminations are based on the degradative abilities of the microorganisms [START_REF] Moreno | Biological treatments for contaminated soils: hydrocarbon contamination. fungal applications in bioremediation treatment[END_REF].

Bioremediation technologies can be broadly classified as ex situ and in situ [START_REF] Boopathy | Factors limiting bioremediation technologies[END_REF]. Ex situ technologies are those treatments which involve the physical removal of the contaminated material for treatment process. In situ techniques involve treatment of the contaminated material in place. Examples of in situ and ex situ bioremediation are

• Land farming: Solid-phase treatment system for contaminated soils: may be done in situ or ex situ.

• Composting: Aerobic, thermophilic treatment process in which contaminated material is mixed with a bulking agent; can be done using static piles or aerated piles.

• Bioreactors: Biodegradation in a container or reactor; may be used to treat liquids or slurries.

• Bioventing: Method of treating contaminated soils by drawing oxygen through the soil to stimulate microbial activity.

• Biofilters: Use of microbial stripping columns to treat air emissions.

• Bioaugmentation: Addition of bacterial cultures to a contaminated medium; frequently used in both in situ and ex situ systems.

• Biostimulation: Stimulation of indigenous microbial populations in soils or ground water by providing necessary nutrients.

• Intrinsic bioremediation: Unassisted bioremediation of contaminant; only regular monitoring is done.

• Pump and treat: Pumping ground water to the surface, treating, and reinjecting.

In the present work we are interested in the bioremediation of natural water resources with ex situ technologies, in particular, using biorectors. A bioreactor is a device that supports a biologically active environment, where chemical reactions (usually microbial fermentation or biotransformation) take place as a result of microbial metabolism, and can be used for bioremediation, cell or tissue culture, or the generation of derivatives or end products of interest of chemical processes, such as enzymes. The scheme of treatment of a water resource consists roughly of three steps: primary, secondary and tertiary treatments.

• Primary or mechanical treatment is designed to remove gross, suspended and floating solids from the effluent. It includes screening to trap solid objects and sedimentation by gravity to remove suspended solids. This level is sometimes referred to as mechanical treatment although chemicals are often used to accelerate the sedimentation process. Primary treatment can reduce the biochemical oxygen demand of the incoming wastewater by 20-30% and the total suspended solids by some 50-60%. Primary treatment is usually the first stage of wastewater treatment.

• Secondary or biological treatment removes the dissolved organic matter that escapes primary treatment. This is achieved by microorganisms consuming the organic matter as food, and converting it to carbon dioxide, water, and energy for their own growth and reproduction. The biological process is performed in bioreactor tanks, then followed by additional settling tanks (secondary sedimentation) to remove more of the suspended solids, or the excess of microorganisms. About 85% of the suspended solids and biochemical oxygen demand can be removed by a plant with secondary treatment. Secondary treatment technologies include the basic activated sludge process, the variants of pond and constructed wetland systems, trickling filters and other forms of treatment which use biological activity to break down organic matter.

• Tertiary treatment is additional treatment beyond secondary. Tertiary treatment can remove more than 99% of all the impurities from sewage, producing an effluent of almost drinking-water quality. The related technology can be very expensive, requiring a high level of technical know-how and well trained treatment plant operators, a steady energy supply, and chemicals and specific equipment which may not be readily available. An example of a typical tertiary treatment process is the modification of a conventional secondary treatment plant to remove additional phosphorus and nitrogen.

(See [START_REF] Ramalho | Introduction to Wastewater Treatment Processes[END_REF] for an extensive description of these processes). The scheme of the water treatment process is depicted in Figure 1.1. In the present thesis we focus on the optimization of the operation of the secondary treatment.

Tertiary treatment Primary treatment

Water resource From the point of view of mathematical modelling, biological reactors can be divided into two major classes [START_REF] Dochain | Automatic control of bioprocesses[END_REF]:

• stirred tank reactors, such as chemostats, sequencing batch reactors, etc., for which the reacting medium is homogeneous (this is achieved with an agitator) and the reaction is described by ordinary differential equations;

• reactors with a spatial concentration gradient, such as fixed beds, fluidized beds, air lifts, settlers, etc., for which the reactor is described by partial differential equations.

From the point of view of its operation, the biological reactors can be classified according the way in which the liquid exchange is performed. We can distinguish three main modes.

1. Discontinuous or batch, 2. Semi-continuous, sequencing-batch, or fed-batch, 3. Continuous or chemostat.

The batch reactor consists of an operation mode in which the culture medium with the nutrient is introduced to the tank at the begining of the process. After that point, nothing is added to or removed from the system, up to the end of the process, which is performed at a constant volume.

The sequencing-batch reactor (SBR) consists of a tank which is fed with a supply of nutritive elements at the beginning and during the process. This type of system is widely used in the industrial and municipal wastewater treatment plants. Basic fill-and-draw treatment systems such as the SBR have been used since the 19th century. According to [START_REF] Irvine | Sequencing batch biological reactors: An overview[END_REF][START_REF] Irvine | Sequencing batch reactors for biological waste-water treatment[END_REF], the cycle of operation of this type of systems consists of five basic stages:

• fill: water with nutrient is received in the tank,

• react: the desired reactions take place,

• settle: the microorganisms are separated from the treated water

• draw: the treated effluent is discharged,

• idle: period comprised between the discharge of the tank and its refilling.

The chemostat mode is the most widely used in wastewater treatment and bioremediation of water resources. The main characteristic is that the culture volume is constant, water with nutrient being added continuously to the system, and treated water being removed at the same rate. The reactor is operated in such a way to keep the system in a steady state by adding a constant nutrient concentration at a fixed rate, although the inflow rate can be controlled to change the desired steady state. The chemostat was invented independently by Monod, and Novick and Szilard in 1950, as a mean to study the steady behavior of microorganisms by continuously adding a constant nutrient concentration, and since its invention it has been a key tool to study microbial dynamics [START_REF] Monod | La technique de la culture continue, théorie et applications[END_REF][START_REF] Smith | The theory of the chemostat[END_REF].

Mathematical models of bioreactors and classical results

Mathematical models

Ever since their conception, bioreactors have been widely used to study microbial population dynamics [START_REF] Monod | Recherches sur la croissance des cultures bacteriennes[END_REF][START_REF] Monod | La technique de la culture continue, théorie et applications[END_REF][START_REF] Novick | Experiments with the chemostat on spontaneous mutations of bacteria[END_REF], and the use of bioreactors for wastewater treatment and bioremediation has been widely studied for the last 40 years [3,[START_REF] D'ans | Optimal control of bacterial growth[END_REF][START_REF] D'ans | Time-optimal control for a model of bacterial growth[END_REF][START_REF] Gajardo | Minimal time sequential batch reactors with bounded and impulse controls for one or more species[END_REF][START_REF] Van Impe | Optimal adaptive control of fedbatch fermentation processes[END_REF][START_REF] Van Impe | Optimal adaptive control of biotechnological processes[END_REF][START_REF] Kittisupakorn | Comparison of optimisation based control techniques for the control of a cstr[END_REF][START_REF] Moreno | Optimal time control of bioreactors for the wastewater treatment[END_REF][START_REF] Smets | Optimal adaptive control of (bio)chemical reactors: past, present and future[END_REF][START_REF] Smets | Optimal control of (bio-)chemical reactors: generic properties of time and space dependent optimization[END_REF][START_REF] Soukkou | Optimal control of a CSTR process[END_REF][START_REF] Srinivasan | Dynamic optimization of batch processes: I. characterization of the nominal solution[END_REF]. Typically, there are three types of models for the analysis of the system bioreactor-water resource: the simplest from the mathematical point of view, that relies on the hypotheses of homogeneity of the concentrations in the resource and/or in the bioreactor as well as instant mixing (among others), is a system of nonlinear ordinary differential equations (ODE) [START_REF] D'ans | Optimal control of bacterial growth[END_REF][START_REF] D'ans | Time-optimal control for a model of bacterial growth[END_REF][START_REF] Smith | The theory of the chemostat[END_REF]. Later, models with partial differential equations (PDE) that take into account fluid dynamics, mass conservation and pollutant diffusion have been introduced to address the inhomogeneity of the pollutant in the resource [4] or the bioreactor vessel [START_REF] Diehl | Control of an ideal activated sludge process in wastewater treatment via an ode-pde model[END_REF][START_REF] Diehl | A reduced-order ODE-PDE model for the activated sludge process in wastewater treatment: classification and stability of steady states[END_REF]; this effect naturally appears due to the speed of the reactions in the bioreactor and the slow diffusion speed as well as the slow mixing in large water resources; this is typically the situation of large scale reactors or large scale water resources, fluidized beds, settlers among others. The third kind of models are the stochastic models, that take into account the uncertainty of different types of variables [START_REF] Campillo | Weak convergence of a mass-structured individual-based model[END_REF][START_REF] Campillo | Stochastic modeling of the chemostat[END_REF][START_REF] Collet | Stochastic models for a chemostat and long-time behavior[END_REF][START_REF] Imhof | Exclusion and persistence in deterministic and stochastic chemostat models[END_REF].

Establishing a mathematical model for the dynamics in a bioreactor is not a simple task because of the large number of interconnected variables involved in the process, for instance, pH, temperature, aereation, nutrient, different bacterial species and nutrient substances, end products, etc. Thus, the choice of the number of reactions to be considered and components which intervene in these reactions is very important for modeling. It will be carried out based on the knowledge that we have on the process and measurements which could have been carried out. The reaction scheme conditions the structure of the model. It will thus have to be chosen with parsimony, bearing in mind the objectives of the model and the precision which is expected. The required number of reactions and the reaction scheme can be determined directly from a set of available experimental data [START_REF] Dochain | Automatic control of bioprocesses[END_REF].

The mathematical model of the dynamics in the stirred tank reactor rely on the main concept of mass balance [START_REF] Dochain | Automatic control of bioprocesses[END_REF][START_REF] Dochain | Dynamical modelling and estimation in wastewater treatment processes[END_REF][START_REF] Smith | The theory of the chemostat[END_REF], which can be broadly expressed as (

Suppose that X denotes a microbial species that will degrade a substrate denoted by S. This reaction is schematically described in the following form:

S r(•) -→ X. (1.2)
in which X plays the role of an autocatalyst, i.e., it is both a product and a catalyst. The rate at which the process occurs is related to enzyme kinetics, and depends on the concentration of substrate and biomass. Typically, the reaction rate r(•) is linear with respect to the biomass X, and it is assumed that microorganisms have uniform access to the substrate. With these assumptions, the reaction rate has the form r(S, X) = µ(S)X, where the function µ(•) is called the specific growth rate [START_REF] Dochain | Dynamical modelling and estimation in wastewater treatment processes[END_REF]. In practice, growth rate functions are obtained experimentally in laboratories.

There are two widely used expressions for the growth function µ(•). The first expression, due experimental works of Monod [START_REF] Monod | Recherches sur la croissance des cultures bacteriennes[END_REF][START_REF] Monod | La technique de la culture continue, théorie et applications[END_REF], is also called Michaelis-Menten formulation (from enzyme dynamics), has the form

µ(s) = µ max s K S + s . (1.3)
In (1.3), µ max is the maximum specific growth rate (in units of [1/h]) and K S is the halfsaturation constant (in units of [g/l]). This formula models the saturation or limited growth with respect to the substrate.

The second expression for the growth function takes into account the effect of inhibition of the growth of the microbial specie with respect to the excess of substrate. This expression is called the Haldane growth law:

µ(s) = μs K S + s + s 2 K l .
(1.4)

In (1.4) the term K l is the inhibition constant. This function has a unique maximum s † := √ K s K l , and is increasing in [0, s † ) and decreasing in the interval (s † , ∞), as Figure 1.4 shows. The main assumptions for the stirred bioreactors are

s [g/L] [1/h] [1/h] s [g/L]
• the culture vessel is well stirred, allowing to consider homogeneity of the nutrient and micro-organisms distributions inside the reactor,

• all the other significant parameters (e.g., temperature) affecting growth are kept constant,

• the reaction scheme summarizes the distribution of mass and flows between various reactions intervening in the process

• reproduction is proportional to nutrient uptake.

Let s r (t), x r (t), and v r (t) denote the concentration of nutrient, concentration of microorganisms, and the culture volume in the culture vessel at time t (in units of [g/l] and [l] respectively). Thus s r (t)v r (t) denotes the amount of nutrient in the vessel at time t; analogously, x r (t)v r (t) denotes the amount of biomass at time t. Denote by Q in and Q out the input and output flow rates respectively. The mathematical equations that model the time evolution of the quantities of interest in the system can be obtained by a mass balance of the form (1.1), whose expression is

     d dt (s r v r ) = Q in s in -Q out s r - 1 Y µ(s r )x r v r , d dt (x r v r ) = µ(s r )x r v r -Q out x r . (1.5)
The term Y is the yield coefficient that reflects the conversion of nutrient to microorganism, that is, the quantity of biomass which is produced when one unit of substrate is consumed by the reaction (1.2); it is assumed to be constant. Along with the mass balance equations we consider the variation of culture volume given by the equation

d dt v r = Q in -Q out . (1.6)
Finally, the expression of the equations of the bioreactor is

           ṡr = - 1 Y µ(s r )x r + Q in v r (s in -s r ), ẋr = µ(s r )x r - Q in v r x r , vr = Q in -Q out .
(1.7)

According to the operation mode, the inflow and outflow rates are

• Batch: Q in = Q out = 0.
• SBR: Q out = 0.

• Chemostat:

Q in = Q out .

Classical results on chemostats

Chemostats have been the object of extensive studies since the early works of Monod in the 40's [START_REF] Monod | Recherches sur la croissance des cultures bacteriennes[END_REF] and its work La technique de la culture continue [START_REF] Monod | La technique de la culture continue, théorie et applications[END_REF], the contributions of Novick and Szilard [START_REF] Novick | Experiments with the chemostat on spontaneous mutations of bacteria[END_REF] in 1950, and Herbert et al [START_REF] Herbert | The continuous culture of bacteria; a theoretical and experimental study[END_REF] in 1956. In [START_REF] Smith | The theory of the chemostat[END_REF] there is an extensive study of the chemostat with the main mathematical results up to 1995.

Let us remind that for the chemostat the culture volume is constant, so the input flow rate Q in and output flow rate Q out are equal at every time instant. Replacing in (1.2.2), we obtain the simple model of chemostat

       ẋr = µ(s r ) - Q in v r x r , ṡr = - 1 Y µ(s r )x r + Q in v r (s in -s r ).
(1.8)

The quantity D := Q in /v r is called dilution rate or washout rate; it has units of 1/t. Under an appropriate change of variable (x r = x r /Y , but we choose to keep x r instead of xr as variable), the equations of the chemostat take the form ẋr = (µ(s r ) -D) x r , ṡr = -µ(s r )x r + D(s in -s r ).

(1.9)

A basic mathematical result on the chemostat concerns the asymptotic behavior of the system. Depending on the behavior of the growth function and the dilution rate, there may exists a unique or several equilibrium points. Suppose that the chemostat is operated at a constant dilution rate D > 0. The trivial equilibrium point E 0 := (0, s in ) is called the washout. This equilibrium, in which there are no microorganisms in the reactor, is not desirable from the point of view of the depollution process. The other equilibrium points are the solutions of the system of equations

µ(s ⋆ r ) = D, x ⋆ r = s in -s ⋆ r .
(1.10)

We notice that all the equilibrium points rely on the straight line {x r + s r = s in }. Indeed, this is stated by the second equation of (1.10), provided that D > 0. In Figure 1.5 we show the phase portrait of the single species chemostat model with a Monod growth function with µ max = 1, The following proposition summarizes the properties of the equilibrium points.

K S = 1, s in = 4, for different dilution coefficients D ∈ [0, 1]. 0 0.5 1 1.5 2
Proposition 1.1 Suppose that the constant dilution rate D > 0 is such that the equation µ(s ⋆ r ) = D has a solution s ⋆ r ∈ (0, s in ). • The equilibrium point E 0 always exists (independently of D) and it is unstable.

• For all s ⋆ r < s in solution of µ(s ⋆ r ) = D there exists an equilibrium point E 1 := (x ⋆ r , s ⋆ r ) = (s in -s ⋆ r , s ⋆ r ). Any of these equilibrium points is locally asymptotically stable if µ ′ (s ⋆ r ) > 0; if µ ′ (s ⋆ r ) ≤ 0 the equilibrium point is unstable.
The proof of this proposition relies on the eigenvalues of the Jacobian matrix of the system (1.9). Notice that we need to impose the condition s ⋆ r ∈ [0, s in ] not just for mathematical reasons but for the interpretation of the variables; indeed, the concentration of biomass at the equilibrium is x ⋆ r = s in -s ⋆ r , which only has a biological meaning if it is a nonnegative quantity.

For different types of uptake functions µ(•), Proposition 1.1 gives different qualitative re-sults.

• For an increasing growth function (like the Monod's law), if µ(s in ) < D, the unique locally asymptotically stable equilibrium point is the washout. On the other hand, for D < µ(s in ), there exits a unique solution to the equation µ(s ⋆ r ) = D in the interval (0, s in ); the washout becomes unstable, and the point

(x ⋆ r , s ⋆ r ) = (s in -µ -1 (D), µ -1 (D)
) is locally asymptotically stable (see Figure 1.6).

s [g/l] s in [1/h] D 1 D 2 (s in ) Figure 1.6.
Equilibrium points for the Monod growth function. We see that D 1 < µmax, and then the concentration of substrate at equilibrium is s † r = λ 1 . For D 2 , the equilibrium is the washout.

• For an increasing-decreasing growth function (like the Haldane's law), there may exist 0, 1 or 2 solutions to the equation µ(s ⋆ r ) = D. Denote by s † the concentration of maximum efficiency (the maximum of the Haldane growth function), and suppose that s † < s in .

-If D > µ(s † ), there is no solution to the equation µ(s ⋆ r ) = D. In such case, the washout is the only equilibrium point, and is is locally asymptotically stable.

-If D = µ(s † ), there exist two equilibrium points: the washout (unstable), and the point

E 1 = (s in -s † , s † ) (unstable).
-If D < µ(s † ), there exist two solutions λ 1 < s † and λ 2 > s † to the equation µ(s ⋆ r ) = D. If λ 2 < s in , then there are two nontrivial equilibrium points E 1 = (s in -λ 1 , λ 1 ) (locally asymptotically stable) and

E 2 = (s in -λ 2 , λ 2 ) (unstable); If λ 2 > s in , the only nontrivial equilibrium point is E 1 . The washout is unstable.
The model can be extended in a natural way (obtained via the mass balance equations) to a model with several species. The way to model the system depending on the interactions between the different species in the chemostat. Suppose that the only interaction between a system of n species is due to the competition for a single nutrient. The equations are

     ẋr,i = (µ i (s r ) -D) x r,i , i = 1, . . . , n ṡr = - n j=1 1 Y j µ j (s r )x r,j + D(s in -s r ). (1.11)
The constants Y i correspond to the yield coefficients of the transformation of substrate into biomass of specie i. Suppose that the growth functions are increasing or increasing-decreasing

[1/h] s [g/l] D 3 D 2 D 1 s in (s in ) Figure 1.7.
Equilibrium points for the Haldane growth function. We see that for D 1 < µ(s † ) there exist two equilibrium points λ 1 (locally asymptotically stable) and λ 2 (unstable) besides the washout (unstable). For D 2 = µ(s † ) the equilibrium concentration is s ⋆ r = s † (unstable), and the washout is unstable. For D 3 , the equilibrium is the washout (stable).

(like Monod's or Haldane's law). Define, for each i ∈ {1, . . . , n}, the set Λ i (D) := {s r > 0 | µ i (s r ) > D}. Because of the assumptions on the growth functions, if Λ i (D) is not empty, then it is an interval Λ i (D) = (λ i , ν i ), where ν i can take the value +∞. Suppose the particular situation in which that all the λ i 's that are finite and different, and without loss of generality, that the numbers λ 1 , . . . , λ n are ordered such that 0

< λ 1 < • • • < λ m < s in < λ m+1 < • • • < λ n
for some m ∈ {1, . . . , n}. The equilibrium points of this system are

• the washout: E 0 = (s ⋆ r , x ⋆ r,1 , . . . , x ⋆ r,n ) = (s in , 0, . . . , 0); • E i = (s ⋆ r , x ⋆ r,1 , . . . , x ⋆ r,n ) = (λ i , 0, . . . , Y i (s in -λ i ), .
. . , 0) for i = 1, . . . , m (in this case, there is only substrate and specie i; this point exists because λ i < s in ); for i > m, E i does not exist.

• E i = (s ⋆ r , x ⋆ r,1 , . . . , x ⋆ r,n ) = (ν i , 0, . . . , Y i (s in -ν i ), .
. . , 0) for i = 1, . . . , m ; (only substrate and specie i), that exists only if ν i < s in and in such case it is unstable. For i > m, E i does not exist.

A well known result on the coexistence of several species in chemostats is the competitive exclusion principle. This principle states that in a chemostat with constant operation parameters D and s in , several species cannot coexist [START_REF] Hardin | The competitive exclusion principle[END_REF][START_REF] Smith | The theory of the chemostat[END_REF] (in a more general way, it states that when two species compete for the same critical resources within an environment, one of them will eventually outcompete and displace the other). This principle has a mathematical sustent in the work of [START_REF] Butler | A mathematical model of the chemostat with a general class of functions describing nutrient uptake[END_REF] (as stated in [START_REF] Smith | The theory of the chemostat[END_REF]).

Theorem 1.2 ([12]

) Define I = ∪ i=1,...,m (λ i , ν i ), with λ i , ν i as before. Suppose that I is a non-empty open interval. Then I = (λ i , ν j ), for some j.

1. If λ 1 < s in < ν j , then E 1 attracts all solutions with x r,i (0) > 0. The first point of Theorem 1.2 corresponds to the competitive exclusion principle. The second point of Theorem 1.2 states another phenomenon: the possibility that an excess of nutrient can lead to the washout of all populations. Theorem 1.2 does not hold if some of the λ i 's coincide.

Control problems arise naturally in the contex of biochemical processes, biotechnology, and wastewater treatment. Bioreactor control provides special challenges due to significant process variability, the complexity of biological systems, the need, in many cases, to operate in a sterile environment, and the relatively few real-time direct measurements available that help define the state of the culture. The study of control of the bioreactor is based on suitable manipulation of its operation parameters, the dilution rate, the input nutrient concentration, or the drop of the well-mixed hypothesis. The first works on the control of chemostats address the problem of coexistence of different species under periodic changes on the dilution rate, when a fraction of the biomass and growing medium are periodically harvested, and when both the dilution rate and the concentration of the substrate in the feed are varied simultaneously and in a periodic manner [START_REF] Butler | A mathematical model of the chemostat with periodic washout rate[END_REF][START_REF] Stephanopoulos | The growth of competing microbial populations in a cstr with periodically varying inputs[END_REF].

Since the work by [START_REF] D'ans | Optimal control of bacterial growth[END_REF], the optimization of bioreactor operation has received great attention in the literature, see [START_REF] Alford | Bioprocess control: Advances and challenges[END_REF]3,[START_REF] Rani | Control of fermenters: a review[END_REF] for reviews of the different optimization techniques that have been used in bioprocesses. Among them, the theory of optimal control has proven to be a generic tool for deriving practical optimal rules [START_REF] Van Impe | Optimal adaptive control of biotechnological processes[END_REF][START_REF] Smets | Optimal adaptive control of (bio)chemical reactors: past, present and future[END_REF][START_REF] Smets | Optimal control of (bio-)chemical reactors: generic properties of time and space dependent optimization[END_REF].

Typically, the optimal control of continuous processes usually involves a two-step procedure. First, the optimal steady state is determined as a nominal set point that maximizes a criterion [START_REF] Soukkou | Optimal control of a CSTR process[END_REF][START_REF] Spitzer | Maximization of steady-state bacterial production in a chemostat with ph and substrate control[END_REF]. The benefit of operating a periodic control about the nominal point can be analyzed [START_REF] Abulesz | Periodic optimization of continuous microbial growth processes[END_REF][START_REF] Ruan | Comparison of several periodic oscillations of a continuous fermentation process[END_REF]. Then, a control strategy that drives the state about the nominal set point from any initial condition is searched for [START_REF] Kittisupakorn | Comparison of optimisation based control techniques for the control of a cstr[END_REF], possibly in the presence of model uncertainty using extremum seeking techniques [6,[START_REF] Marcos | Adaptive extremum seeking control of a continuous stirred tank bioreactor with Haldane's kinetics[END_REF][START_REF] Wang | Optimizing bioreactors by extremum seeking[END_REF][START_REF] Zhang | Adaptive extremum seeking control of continuous stirred-tank bioreactors[END_REF].

Classical results on SBRs

SBRs are typically used to treat municipal and industrial wastewaters, particularly in areas characterized by low or varying flow patterns. Improvements in equipment and technology, especially in aeration devices and computer control systems, have made SBRs a viable choice over the conventional activated-sludge system. These plants are very practical for a number of reasons. For instance, in areas where there is a limited amount of space, treatment takes place in a single basin instead of multiple basins, allowing for a smaller footprint. The treatment cycle can be adjusted to undergo aerobic, anaerobic, and anoxic conditions in order to achieve biological nutrient removal, including nitrification, denitrification, and some phosphorus removal. SBRs offer a cost-effective way to achieve lower effluent limits.

Flow-paced batch operation is generally preferable to time-paced batch or continuous inflow systems. Under a flow-paced batch system, a plant receives the same volumetric loading and approximately the same organic loading during every cycle. Under a time-paced mode, each basin receives different volumetric and organic loading during every cycle, and the plant is not utilizing the full potential of this treatment method which is the ability to handle variable waste streams. Time-paced operation can lead to under-treated effluent if the cycle time is not adjusted. For an SBR to be effective, the plant must have proper monitoring, allow operators to adjust the cycle time, and have knowledgeable operators who are properly trained to make the necessary adjustments to the cycle [START_REF]Sequencing batch reactor design and operational considerations[END_REF].

Since the SBR is a time oriented system, the control objectives are usually to optimize trajectories to attain a prescribed target in finite time [START_REF] Hong | Optimal substrate feeding policy for fed batch fermentation with substrate and product inhibition kinetics[END_REF][START_REF] Irvine | Sequencing batch reactors for biological waste-water treatment[END_REF][START_REF] Johnson | The control of fed-batch fermentation processes: a survey[END_REF][START_REF] Kurtanjek | Optimal nonsingular control of fed-batch fermentation[END_REF][START_REF] Lim | Computational algorithms for optimal feed rates for a class of fed-batch fermentation: numerical results for penicillin and cell production[END_REF][START_REF] Shioya | Optimization and control in fed-batch bioreactors[END_REF][START_REF] Tsoneva | Augmented lagrange decomposition method for optimal control calculation of batch fermentation processes[END_REF] or to maximize production at a given time [START_REF] Gajardo | Minimal time sequential batch reactors with bounded and impulse controls for one or more species[END_REF][START_REF] Mazouni | Optimal time switching control for multi-reaction batch process[END_REF][START_REF] Moreno | Optimal time control of bioreactors for the wastewater treatment[END_REF][START_REF] Srinivasan | Dynamic optimization of batch processes: I. characterization of the nominal solution[END_REF]. The single species model of SBR is

           ṡr = - 1 Y µ(s r )x r + u v r (s in -s r ), ẋr = µ(s r )x r - u v r x r , vr = u, (1.12)
where u is the inflow rate, that is usually the control variable and is bounded between 0 and u max . This system has the constraint that v r ≤ v max .

In [START_REF] Moreno | Optimal time control of bioreactors for the wastewater treatment[END_REF] the problem of attaining a prescribed pollution level s out with the tank at its maximum capacity in minimal time is studied, for the single species case. The author considers a varying inflow substrate concentration s in to be a function of time. By means of Green theorem the optimal strategy is characterized in two cases depending on the behavior of the growth function. In the case that the growth function is monotone (Monod type), the optimal control consists on filling the tank at maximum speed u = u max until the tank is full v = v max , and then wait until s r ≤ s out . In the case that the growth function is nonmonotone with one maximum s † , the optimal strategy consists on a feedback control such that brings the pollutant concentration s r as fast as possible to s † and to keep the system in that state up to the time when v r = v max , and then wait.

In [START_REF] Gajardo | Minimal time sequential batch reactors with bounded and impulse controls for one or more species[END_REF], the authors extend the previous work to the case of several species and allowing impulsional controls, that is, adding arbitrarily large amounts of water in arbitrarily small time instants. The SBR with impulse is modeled by a system of equations without impulse in a fictitious time τ that allows to extend the real time t at the points where an impulse is made.

                 d dτ s r = -r n j=1 1 Y j µ j (s r )x r,j + u v r (s in -s r ), d dτ x r,i = rµ i (s r )x r,i - u v r x r,i , i = 1, . . . , n, d dτ v r = u, (1.13)
Here, u is the control variable that represents the inflow rate, and r is a control variable that takes the value 0 when there is an impulse (that lasts for the fictitious time intervals [τ 0 , τ 1 ] such that τ 1 τ 0 u(τ ) = v + -v -explains the difference of volumes before and after the impulse) and takes the value 1 otherwise. The behavior of the real time t as a function of τ is shown in Figure 1.8 The authors define the immediate one impulse (IOI) strategy as the control that consists of performing an impulse that fills the tank in time zero, and then to wait; also, the singular arc (SA) strategy to a level s ⋆ is introduced for the case in which the growth function has a nonmonotone behavior, consisting on a feedback control such that brings the pollutant concentration s r to s ⋆ in an impulse, then keeping the system in that state up to the time when v r = v max , and then wait, unless s ⋆ < s out , in which case there is a final impulse. With the Pontryagin's Maximum Principle and Hamilton-Jacobi-Bellman techniques, the authors characterize the cost of the one impulse strategy, and then prove that for the single species case with monotonic growth function the IOI strategy is optimal; for the nonmonotonic growth function, the SA strategy to the level s ⋆ = s † (the maximum of the growth function) is optimal. For the two species case, the ratio of the derivatives of the growth functions is strictly monotone, the optimal control is also characterized as the IOI strategy or the SA strategy to a level s ⋆ ∈ (s out , s in ).

Stochastic models of bioreactors

Stochastic models for a bioreactor are those that incorporate the effect of uncertainty in the behavior of the corresponding variables. The uncertainty in a model is due to a combination of

• uncertainty in input variables,

• uncertainty in parameter values,

• uncertainty on model structure.

Typically the models that can be found in the literature address the first type of incertainty that typically correspond to random variations in the inflow rate Q in and/or the input pollutant concentration s in , or to the second type of uncertainty, related to the uncertainty on the growth function parameters, leading to studies on parameter estimation. Lately, the uncertainty on model structure in the chemostat has called the attention of some authors [START_REF] Campillo | Weak convergence of a mass-structured individual-based model[END_REF][START_REF] Collet | Stochastic models for a chemostat and long-time behavior[END_REF][START_REF] Imhof | Exclusion and persistence in deterministic and stochastic chemostat models[END_REF]. The uncertainty on model structure is related to the microscopic mechanism of replication of microorganisms.

In [START_REF] Caraballo | Chemostats with random inputs and wall growth[END_REF], the authors introduce a novelty approach to tackle the problem of random inputs in the model of chemostat, particularly the random nutrient supplying rate or the random input nutrient concentration, with or without wall growth, from the mathematical point of view of random dynamical systems and using the concept of random attractors. The authors obtain results of existence of uniformly bounded non-negative solutions, existence of random attractors, and geometric details of random attractors for different values of parameters.

In [START_REF] Imhof | Exclusion and persistence in deterministic and stochastic chemostat models[END_REF], the authors model the influence of random fluctuations by setting up and analyzing a stochastic differential equation, and show that random effects may lead to extinction in scenarios where the deterministic model predicts persistence, and establish some stochastic persistence results. The single species model presented in that article is

dX 0 = (r -δX 0 -a(X 0 , X 1 ))dt + σ 0 X 0 dW 0 (t), dX 1 = (a(X 0 , X 1 ) -s(X 1 ))dt + σ 1 X 1 dW 1 (t), (1.14)
where X 0 stands for substrate concentration, X 1 denotes the biomass concentration of a microbe species feeding on the substrate; the substrate inflow rate r and the relative substrate outflow rate δ are positive constants; the substrate uptake rate a(X 0 , X 1 ) = µ(X 0 )X 1 , which is equal to the microbe growth rate, is non-negative and strictly increasing in both variables, with a(X 0 , 0) = a(0, X 1 ) = 0 for all X 0 , X 1 (this means that substrate uptake occurs only when substrate and microbes are present); the microbe removal rate s(x 1 ) is non-negative and strictly increasing; and σ 0 , σ 1 ≥ reflect the size of the stochastic effects. The derivation of said model is performed as a limit of a family of discrete Markov Chains whose random effects depend linearly on the respective variable. The authors prove that under suitable assumptions for every initial condition the equation (1.14) has a strong solution defined for every time instant, pathwise uniqueness holds, and with probability one the process stays in the interior of the positive orthant.

In [START_REF] Campillo | Weak convergence of a mass-structured individual-based model[END_REF], the authors propose a model of chemostat where the bacterial population is individuallybased, each bacterium is explicitly represented and has a mass evolving continuously over time. The substrate concentration is represented as a conventional ordinary differential equation. These two components are coupled with the bacterial consumption. Mechanisms acting on the bacteria such as growth, division and washout, are explicitly described, and bacteria interact via consumption. The authors prove the convergence of this process to the solution of an integro-differential equation when the population size tends to infinity. The equations of the model are

               Y t = 1 v r mmax 0 xp t (x)dx, Ṡt = D(s in -S t ) - k v r mmax 0 ρ(S t , x)p t (x)dx, ∂ ∂t p t (x) + ∂ ∂x (ρ(S t , x)p t (x)) + (λ(S t , x) + D)p t (x) = 2 mmax 0 λ(S t , z) z q x z p t (z)dz, for x ∈ [0, m max ].
Here p t (x) is the density of population with respect to its mass at a time t. Hence m 0 p t (x)dx is the number of cells which mass is between m 0 and m 1 and the average biomass concentration at time t is Y t ; ρ(s, x) and λ(s, x) are respectively the growth function and the division rate of a bacterium of mass x with a substrate concentration s (typically ρ(s, x) = µ(s)x), the mass distribution of the daughter cells is represented by the probability density function q(α) on [0, 1].

In [START_REF] Collet | Stochastic models for a chemostat and long-time behavior[END_REF], the authors introduce two stochastic chemostat models consisting of a coupled population-nutrient process reflecting the interaction between the nutrient and the bacteria in the chemostat with finite volume, where the nutrient concentration evolves continuously but depends on the population size, while the population size is a birth-and-death process with coefficients depending on time through the nutrient concentration. The nutrient is shared by the bacteria and creates a regulation of the bacterial population size. The latter and the fluctuations due to the random births and deaths of individuals make the population go almost surely to extinction. The authors study the long-time behavior of the bacterial population conditioned to non-extinction, prove the global existence of the process and its almost-sure extinction; the existence of quasi-stationary distributions is obtained based on a general fixed-point argument. The authors also prove the absolute continuity of the nutrient distribution when conditioned to a fixed number of individuals and the smoothness of the corresponding densities. The equation of the substrate concentration is given by

d dt S(t) = D(s in -S(t)) -b(S(t))N (t),
linked to the process of the population size by the infinitesimal generator [START_REF] Ethier | Markov processes[END_REF] of the process Z = (Z(t) := (N (t), S(t)), t ≥ 0) given by

L f (n, s) = b(s)nf (n + 1, s) + (D + d(s))nf (n -1, s) -(b(s) + D + d(s))nf (n, s) + (D(s in -s) -nb(s)) ∂ ∂s f (n, s).
In this model b(d) is the birthe rate per individual (typically the Monod growth function) and d(s) is the background death rate per individual, that depends on the nutrient. In this model it is important to remark the fact that there is almost sure extinction of the biomass, which is a result that does not hold for the deterministic models.

Model of inhomogeneous lake

This thesis studies the bioremediation of natural water reservoirs, such as lakes, ponds, lagoons, underground waters, water tables, etc., with the use of a bioreactor. From a mathematical point of view, this means that we couple the dynamics of a bioreactor with the dynamics of the pollutant concentration in the water resource. We suppose that there is a settler that perfectly separates the biomass from the effluent of the bioreactor in negligible time before returning the treated water to the resource.

Since typically water resources have a large volume and spatial constraints, inhomogeneity of the distribution of the pollutant in the resource naturally arises. The connection between the bioreactor and the water reservoir induces water movements in the latter. We model the water resource as a bounded open domain Ω ⊆ R n , with n = 2 for lakes of small depth, or n = 3 for a deep lake. We suposse that the velocity field v(t, x) depends only of the inflow/outflow discharges from the bioreactor, the water viscosity is constant and homogeneous, and the fluid is incompressible. Denote Γ in ⊆ ∂Ω is the area that receives the effluent of the reactor and Γ out ⊆ ∂Ω is the area from which water is taken to the bioreactor. Denote p the pressure of the fluid. The equations that model the velocity field in the resource are given by the Navier-Stokes equations of an incompressible fluid [4]:

   ∂ ∂t v + v • ∇ v + ∇p -ν v ∆v = 0, (t, x) ∈ (0, ∞) × Ω, ∇ • v = 0, (t, x) ∈ (0, ∞) × Ω.
(1.15)

With this equations we provide the initial condition

v(t = 0, x) = 0, x ∈ Ω,
and the boundary conditions for Γ in , Γ out , and

Γ 0 := ∂Ω\(Γ in ∪ Γ out ):      v(t, x) = Q(t) v in (x), x ∈ Γ in , v(t, x) = Q(t) v out (x), x ∈ Γ out , v(t, x) = 0, x ∈ Γ 0 .
The functions v in and v out are unitary parabolic vector fields that describe the velocity profile on Γ in and Γ out parallel to the outwards normal n and satisfy

- Γ in v in • ndS = Γout v out • ndS = 1.
The previous model allows to characterize the velocity field in the water resource as a function of the pumping speed Q(t). Numerical simulations can be performed to compute the velocity field independently of the pollutant concentration (since we have assumed that changes in the pollutant concentration do not affect the viscosity of the fluid). If the resource has a large depth, a 3 dimensional model can be considered. This is numerically costly because of the size as well as the geometry of the resource. If the depth is small, a 2 dimensional model gives good results. Now, denote by s l (t, x) the concentration of pollutant in the resource at the point x ∈ Ω in the time instant t. Suppose that the diffusivity coefficient ν s of the pollutant is constant and homogeneous. The equation that models the time evolution of the concentration of pollutant in the resource are composed by an equation of mass conservation of the pollutant (given the velocity field v previously computed) are

∂ ∂t s l -ν s ∆s l + v • ∇s l = 0, (t, x) ∈ (0, ∞) × Ω, (1.16)
with the initial condition

s l (t = 0, x) = s 0 l (x), x ∈ Ω,
and boundary conditions

               ν s ∂s l ∂n (t, x) = Q(t) |Γ in | v in (x) • n(s l (t, x) -s out r (t)), x ∈ Γ in , ∂s l ∂n (t, x) = 0, x ∈ Γ out , ∂s l ∂n (t, x) = 0, x ∈ Γ 0 .
Numerical simulations of this model in [4] show that the inhomogeneity can be considered as if there were two different zones: a first zone that is being actively treated by the bioreactor in which there could be considered a gradient of concentrations, and another zone that is being depolluted mainly by diffusion with the first zone. This amounts to consider simpler models for treating the inhomogeneity of the pollutant, for instance, compartimental models that consider the resource splitted into two or more zones, each of them with homogeneous concentration, and connected between each other by advection or diffusion coefficients according to geometric specifications.

Figure 1.9. Behavior of the inhomogeneous representation of a lake. From the homogeneous initial distribution of pollutant (on the left), the system evolves up to a point in which two zones are clearly differentiated (on the right). Taken from [4] The motivation of this part of the thesis is to establish optimal control rules for the depollution of water resources in minimal time, with the use of simpler ODE models as an approximation of the inhomogeneous model. In this regard, there exists a tradeoff between the pumping speed and the quality of depollution. On the one hand, if the pumping speed Q is too low, the quality of the depollution will be good but the treatment time will be long since it will take too much time to recirculate the water through the bioreactor. On the other hand, if the pumping speed is too high, the quality of the treatment in the reactor will be poor because the microorganisms do not have enough time to process the pollutant, and it could even lead to the washout in the reactor. Finding an optimal control for minimal time treatment has both a biological and an economical sense. The problem of the depollution of water resources consists in treating the resource by means of a bioreactor by controlling the inflow rate Q and making the pollutant concentration in the resource to decrease under a certain level s considered safe.

In this thesis, we consider three spatially distributed configurations to model the inhomogeneity. The first configuration, that we call the active-dead zones model. This configuration considers the resource split into two zones that are connected to each other, and only one of them is connected to the bioreactor, the other zone being depolluted only by diffusion with the first one. This model applies to confinement in lakes due to geometry, or resources of large proportions in which only a specific area is being treated, among other situations (see Figure 1.10). The third configuration models the gradient of concentrations in the resource by means of considering two separated zones connected in a series configuration. Examples of this configuration are two lakes situated at different heights connected by a water stream. Water is taken from the first zone to the bioreactor, where it is treated; cleaner water from the output of the bioreactor is injected in the second zone, which connects with the first zone (see Figure 1.12). In [4] the authors give a justification of the possibility of considering simpler models of two homogeneous zones as a mean to model the inhomogeneous behavior of the pollutant concentration; the authors perform a parameter estimation of the active-dead zones model to fit the output pollutant concentration of the inhomogeneous lake, showing that a two zones model explains well the inhomogeneity, and justify that this could be extended to several homogeneous zones with different configurations.

We could also consider inhomogeneous models of the chemostat, but in this work it will not be necessary, because of the different time scales that naturally arise in this kind of systems (see Section 1.4). For further references on models of unstirred chemostat, see [START_REF] Baxley | Differential equations and computational simulations ii coexistence in the unstirred chemostat[END_REF][START_REF] Hsu | Dynamics of competition in the unstirred chemostat[END_REF][START_REF] Smith | The theory of the chemostat[END_REF].

Singular perturbations

Coupled systems of compartimental water resources and bioreactors naturally present two time scales. Since the bioreactor vessels are typically much smaller than the water resources that they treat, the changes in the substrate concentration inside the bioreactors are much faster than those in the water resources (see [START_REF] Dochain | Automatic control of bioprocesses[END_REF][START_REF] Dochain | Dynamical modelling and estimation in wastewater treatment processes[END_REF]). This kind of systems in which there is a parameter that introduces such behavior is called singular perturbation model. Singular perturbations cause a multitime-scale behavior of dynamical systems characterized by the presence of slow and fast transients in the system's response.

Let us consider, as an example, the system homogeneous resource-bioreactor [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF]. This system consists of a water resource being depolluted by means of a continuous stirred reactor operated at a pumping speed Q(t) and a settler that performs a perfect separation of the biomass from the treated water at the effluent of the bioreactor, such that there is no biomass in the overflow of the settler. An important assumption that we make is that the separator operates in a negligible time. This assumption can be made because the different time scales between the reactors and the water resource, and would not be valid if the water resource, the bioreactor, and the settler had a similar volume scale. For a more realistic model of the coupling we refer to [START_REF] Diehl | Control of an ideal activated sludge process in wastewater treatment via an ode-pde model[END_REF][START_REF] Diehl | A reduced-order ODE-PDE model for the activated sludge process in wastewater treatment: classification and stability of steady states[END_REF]. Denote by s l and v l the pollutant concentration and water volume in the resource, respectively. The equations that model the dynamics of the system are the chemostat equations coupled with the bass balance for the substrate concentration in the resource (which is the input concentration of the chemostat):

               ẋr = µ(s r ) - Q v r x r , ṡr = -µ(s r )x r + Q v r (s l -s r ), ṡl = Q v l (s r -s l ).
(1.17)

The reasonable assumption that the volume of the reactor is much smaller than the volume of the resource generates two different time scales. Consider the time change τ = t/v l . Under this new time scale, (1.17) becomes

             ǫ d dτ x r = (µ(s r )v r -Q) x r , ǫ d dτ s r = -µ(s r )x r v r + Q(s l -s r ), d dτ s l = Q(s r -s l ). (1.18)
where ǫ := v r /v l is a parameter that indicates the difference of scales bioreactor/resource. If ǫ is small enough, in the natural time scale of the resource the changes in the reactor are almost instantaneous. This effect is called a slow-fast dynamic. It is a well known effect in chemical kynetics [START_REF] Van Breusegem | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF]. In the case that the variables (x r , s r ) reach an equilibrium, this equilibrium will depend on Q and s l . Such an equilibrium is called a quasi-steady-state. For (1.18), the equilibrium concentrations in the bioreactor (supposing that we can discard the washout, which is always possible by choosing an appropriate value of Q) are

s ⋆ r (s l , Q) = µ -1 (Q/v r ), x ⋆ r (s l , Q) = s l -s ⋆ r (s l , Q). (1.19)
We assume here that the value s ⋆ r (s l , Q) defined as a solution of equation (1.19) is well defined and it is unique, for each s l > 0 and appropriate Q > 0. These assumptions are satisfied if the growth function µ(•) is increasing in the interval [0, s l ] and the bioreactor is operated in the range of inflow rates Q such that Q/v r < µ(s l ). For instance, this is satisfied by the Monod's law.

A reduction of order of the system (1.18) can be performed by replacing (1.19) in the dynamic of s r in the time scale of the resource, as if the reactor continuously operates at a quasisteady state. The reduced model is

d dτ sl = Q(s ⋆ r (s l , Q) -sl ).
(1.20)

The reduced model (1.20) is much simpler in mathematical terms than the original model (1.18) and can be seen as an approximation of the latter. The remaining question is relative to the behavior of the approximated system with respect to the original one: Is it true that when ǫ approaches to 0, the solution of (1.18) approaches to the trajectory (x r , sr , sl ), where sl (t) is solution of (1.20) and

xr (t) = x ⋆ r (s l (t), Q), sr (t) = s ⋆ r (s l (t), Q)?
The answer to that question is given by Tychonnof theorem [START_REF] Khalil | Nonlinear Systems[END_REF]Theorem 11.1]. This theorem states that the solution of a system of the form

ẋ = f (t, x, z, ǫ), ǫ ż = g(t, x, z, ǫ), (1.21)
can be approximated by the solution of the reduced system

ẋ = f (t, x, h(t, x), 0), z = h(t, x), (1.22)
when the perturbation ǫ converges to 0, where h(t, x) is a solution of g(t, x, h(t, x), 0) = 0. More precisely, Tychonnof's theorem states the following: consider any time interval [t 0 , t 1 ], the singular perturbation problem (1.21), and let z = h(t, x) be an isolated root of g(t, x, h(t, x), 0) = 0. Assume that the following conditions are satisfied for all (t, x, z -

h(t, x), ǫ) ∈ [t 0 , t 1 ] × D x × D y × [0, ǫ 0 ], for some domains D x ⊆ R n convex and D y ⊆ R m
that contains the origin:

• the functions f , g, their first partial derivatives with respect to (x, z, ǫ), and the first partial derivative of g with respect to t are continuous; the functions (t, x) → h(t, x) and (t, x, z) → ∂g(t, x, z, 0)/∂z are differentiable with continuous derivatives;

• the reduced problem (1.22) has a unique solution x(t) ∈ S, t ∈ [t 0 , t 1 ], with S compact subset of D x ;

• the origin is an exponentially stable equilibrium of the equation

dy dt (τ ) = g(t, x, y(τ ) + h(t, x), 0). (1.23)
uniformly in (t, x); let R y ⊆ D y be the region of attraction of (1.23), and Ω y be a compact subset of R y .

Then, there exists a positive constant ǫ ⋆ > 0 such that for all initial condition z(0) -h(t 0 , x 0 ) in Ω y and 0 < ǫ < ǫ ⋆ , the singular problem (1.21) has a unique solution x(t, ǫ), z(t, ǫ) on [t 0 , t 1 ] and

x(t, ǫ) -x(t) = O(ǫ), z(t, ǫ) -h(t, x(t)) -ŷ(t/ǫ) = O(ǫ) (1.24)
holds uniformly for t ∈ [t 0 , t 1 ], where ŷ(τ ) is solution of (1.23) with initial condition ŷ(t 0 ) = z(t 0 ) -h(t 0 , x(t 0 )). Moreover, given any t b > t 0 there exists another

ǫ ⋆⋆ ≤ ǫ ⋆ such that uniformly for t ∈ [t b , t 1 ] z(t, ǫ) -h(t, x(t)) = O(ǫ) (1.25)
whenever ǫ < ǫ ⋆⋆ .

In Figure 1.13 (Figure 3.6 in Chapter 3) we present some numerical simulation of the trajectories obtained by applying the optimal feedback obtained from the reduced dynamics to the unreduced dynamics of the second spatially distributed model studied in this thesis, for different values of ǫ. We see that the trajectories converge to the trajectory of the reduced model if ǫ converges to zero. An important remark for this kind of systems, that couple the large water resource with a bioreactor, is that the equilibrium substrate concentration in the bioreactor only depends on the inflow rate Q (as we will see in Chapters 2-4), this quantities being linked by the formula

s ⋆ r (Q) = µ -1 (Q/v r ). (1.26)
This allows to consider a change of control variable from Q to s ⋆ r (Q) (if the relation is univoque), and then to control the equilibrium concentration in the bioreactor rather than the speed of the pump. The second part of the equation (1.19) imposes a state constraint on the new control variable:

s ⋆ r (Q) ∈ [0, s l ].
This constraint, if natural, poses problems when applying optimization and optimal control tools to the problem, given that the constraint can be active or not in a not negligible interval of time. This state constraint transforms the problem into a non standard problem.

Contributions of the thesis 1.5.1 Deterministic optimal control for continuous bioremediation processes

In the first part of this thesis we study minimal time strategies for the treatment of pollution in large water volumes, such as lakes or natural reservoirs, using a single continuous bioreactor that operates in a quasi-steady state. The control consists of feeding the bioreactor from the resource, with clean output returning to the resource with the same flow rate. We drop the hypothesis of homogeneity of the pollutant concentration in the water resource by proposing three spatially structured models, each consisting of configurations of two homogeneous zones where the nature of the connections between them as well as the differentiation of the zones may depend on the structure of the problem; for instance, if most of the water stream is passing from one zone to another, this effect can be seen as a transport term between two zones interconnected in a series configuration; this type of configuration is also suitable to model gradients of concentration. Other type of connection that models the effect of (mostly) pure diffusion between two zones is a configuration in parallel. The confinement of one part of the resource due to geography or reactor design also defines clearly differentiated zones that lead to strong inhomogeneity. Naturally, in a highly inomogeneous resource, in which many zones can be differentiated, a compartimental model can be implemented as an approximation of the dynamics of the system. As it was mentioned in Section 1.3, the problem of the depollution of water resources consists in treating the resource by means of a bioreactor leading the pollutant concentration in the resource to decrease under a certain level s > 0 considered safe in environmental terms. In any of the three studied spatial configurations the target of the process will be leading the substrate concentration in both zones under s. If we define as s = (s 1 , s 2 ) ∈ R 2 + the spatial variable that denotes the pollutant concentrations in the first and second zones respectively, the target set will be

T := {s = (s 1 , s 2 ) ∈ R 2 + | s 1 ≤ s, s 2 ≤ s}. (1.27)
In Chapter 2, we study the first configuration, that we call the active-dead zones. This configuration considers the resource split into two zones that are connected to each other, and only one of them is connected to the bioreactor, the other zone being depolluted only by diffusion with the first one. This model applies to confinement in lakes due to geometry, or resources of large proportions in which only a specific area is being treated, among other situations (see Figure 1.10). In this model the effect of adding a second pump which pumps water from one zone to the other and viceversa, enhancing the mixing of the water in the two zones and enhancing the diffusion and mixing of the pollutant through the total resource. The mathematical equations, considering that the bioreactor operates at equilibrium, are

       ṡ1 = v r v 1 µ(s ⋆ r )(s ⋆ r -s 1 ) + D v 1 + Q R v 1 (s 2 -s 1 ), ṡ2 = D v 2 + Q R v 2 (s 1 -s 2 ), (1.28)
where the controls are the output concentration of pollutant in the reactor s ⋆ r and the recirculation speed Q R . Here s ⋆ r is connected with the inflow rate Q by the equation

µ(s ⋆ r ) = Q v r ,
and must satisfy the state constraint s ⋆ r (t) ∈ [0, s 1 (t)] for every time instant t up to the final time. We first recall recent characterizations of optimal policies among feedback controls under the assumption of a uniform concentration in the resource, in which the optimal flow rate decreases with time. Then, we introduce a dead zone in the resource, considering two measurement points. In the first part of this work we study the problem without the effect of recirculation (equivalently, Q R = 0). We prove that the domain {s 2 ≥ s 1 } is positively invariant for the dynamics, allowing us to consider a simpler target set

T = s = (s 1 , s 2 ) ∈ R 2 + | s 2 ≤ s .
With the help of the Pontryagin Maximum Principle, we show that the variant of the dead zone has no influence on the optimal control law. More precisely, the main result is Proposition 1.3 (Chapter 2, Proposition 2.8) The optimal flow rate Q ⋆ , which solves the depollution problem is given by Q

⋆ (s 1 ) = v r µ(ŝ ⋆ r (s 1 ))
, where s ⋆ r (•) is the unique solution of the optimal control problem (P). This solution is characterized by the feedback control defined by ŝ⋆ r (s 1 ) := arg max

s ⋆ r ∈(0,s 1 ) µ(s ⋆ r )(s 1 -s ⋆ r ). (1.29)
In particular, the optimal feedback control s ⋆ r satisfies that s ⋆ r (t) ∈ (0, s 1 (t)), for all t ∈ [0, T opt ], where T opt is the optimal treatment time.

It is important to remark that the obtained optimal state-feedback control has the same form that the optimal control for the homogeneous problem as in [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF] where there is no dead zone, and only depends of the measurements on the active zone.

We also prove that when considering a non-null recirculation speed Q R , the effect of improving the mixing between the two zones by adding a new pump that connects the two zones results in a lower time of treatment, being the ideal case to have a high mixing speed. Proposition 1.4 (Chapter 2, Proposition 2.9) Suppose that (s ⋆ r , Q R ) is an optimal control for the problem with recirculation. Then, for all t ∈ [0, T opt ], the control s ⋆ r is given by the feedback control defined in (1.29), and we have that

Q R (t) = Q max R .
We also prove that the minimal time function, as a function of the diffusion parameter, is decreasing, and this result is independent of the optimal bioreactor's inflow rate Q ⋆ previously obtained. We perform numerical simulations that confirm the theoretical results. Proposition 1.3 and other results of Chapter 2 are part of the published article [START_REF] Gajardo | Bioremediation of natural water resources via optimal control techniques[END_REF].

In Chapter 3, we model the inhomogeneity of the concentration of pollutant in the resource as two clearly differentiated connected patches. The resource is treated by means of a continuously stirred bioreactor connected to each of the patches by pumps. The speed of each pump is controlled independently, which is equivalent to control the total pumping speed and the proportion of the effluent directed towards each zone (see Figure 1.11). The equations that model the dynamics of the pollutant concentration in the two zones, under the assumption that the bioreactor operates at equilibrium, are

     ṡ1 = α r µ(s ⋆ r )(s ⋆ r -s 1 ) + d r (s 2 -s 1 ) , ṡ2 = 1 -α 1 -r µ(s ⋆ r )(s ⋆ r -s 2 ) + d 1 -r (s 1 -s 2 ), (1.30)
with r = v 1 /(v 1 + v 2 ) and d = D/v r , and where α := q 1 /(q 1 + q 2 ) is the proportion of the pump speed directed towards zone 1. In this problem, the target set is

T = s = (s 1 , s 2 ) ∈ R 2 + | max{s 1 , s 2 } ≤ s .
The set where the control variable u = (α, s ⋆ r ) takes values is

U (s) = {(α, s ⋆ r ) | α ∈ [0, 1], s ⋆ r ∈ [0, αs 1 + (1 -α)s 2 ]} , (1.31)
which is a state-dependent control set. Under this constraint, the velocity set is not everywhere convex. Indeed, the dynamic (1.30) can be written as ṡ = F (s, u) + dG(s), and for each

s = (s 1 , s 2 ) such that s 1 = s 2 , the set F (s, U (s)) := {F (s, u) | u ∈ U (s)} is non-convex: the segment {θF (s, (1, ŝ⋆ r (s 1 ))) + (1 -θ)F (s, (0, ŝ⋆ r (s 2 ))) | θ ∈ (0, 1)} does not belong to F (s, U (s)) (where ŝ⋆ r (•) is defined in (1.33); see Figure 1.14).
We study the relaxed problem with a convexification of the control set for which we know there exists a solution, and we prove that the solution of the relaxed control is a solution of the original problem, and it can be characterized as a feedback control. The main results are the following:

s 1 s 2 Velocity set
u ⋆ [s] = (1, ŝ⋆ r (s 1 )) when s 1 > s 2 , (r, ŝ⋆ r (s 1 )) = (r, ŝ⋆ r (s 2 )) when s 1 = s 2 , (0, ŝ⋆ r (s 2 )) when s 1 < s 2 .
(1.32)

where ŝ⋆ r (σ) is the unique solution to

ŝ⋆ r (σ) = arg max s ⋆ r ≥0 µ(s ⋆ r )(σ -s ⋆ r ) (1.33)
In particular, we see that the optimal control satisfies that

s opt r (t) ∈ (0, α(t)s 1 (t) + (1 -α(t))s 2 (t)),
for all t, therefore never attaining an extreme value. We also study the behavior of the optimal treatment time as a function of the diffusion parameter by means of Hamilton-Jacobi-Bellman techniques; it turns out that for some initial conditions in which one of the concentrations is small enough a high diffusion is better than a low diffusion; for other initial conditions, a low diffusion is desirable.

Proposition 1.5 and other results of Chapter 3 are part of the published article [START_REF] Ramírez | Optimal Feedback Synthesis and Minimal Time Function for the Bioremediation of Water Resources with Two Patches[END_REF].

In Chapter 4, we model the gradient of concentrations in the resource by means of considering two separated zones connected in a series configuration. Examples of this configuration are two lakes situated at different heights connected by a water stream. Water is taken from the first zone to the bioreactor, where it is treated; cleaner water from the output of the bioreactor is injected in the second zone, which connects with the first zone. We control the equilibrium concentration of substrate in the bioreactor s r and the recirculation pump speed Q R ∈ [0, QR ] (see Figure 1.12). The mathematical model is given by the equations

       ṡ1 = v r v 1 µ(s r ) + Q R v 1 (s 2 -s 1 ), ṡ2 = v r v 2 µ(s r )(s r -s 2 ) + Q R v 2 (s 1 -s 2 ).
(1.34)

In [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF] the authors consider this problem without recirculation (or equivalently QR = 0) and with a state constraint for the concentration at the output of the reactor given by s r (t) ∈ [0, s 2 (t)]. This is considered as a way to not repollute the water resource while treating it. In this work we allow as a natural state constraint s r (t) ∈ [0, s 1 (t)], allowing the reactor concentration to be lower than its input concentration instead of being lower than the concentration in the zone at its effluent. The fact that the control variable can take an extreme value depending on the state variable poses an extra difficulty, as we show that before arriving to the target set this constraint must be active. We first prove that under the control constraints the set {s 1 > s 2 } is positively invariant, allowing us to consider a simpler target set:

T = s = (s 1 , s 2 ) ∈ R 2 + | s 1 ≤ s .
We prove the existence of functions γ (nonnegative, depending on the adjoint states) and ψ (depending on s) that characterize the behavior of the optimal control. We prove the following proposition.

Proposition 1.6 (Chapter 4, Proposition 4.8) Suppose that s opt is an optimal trajectory, associated to the optimal control u opt = (s opt r , Q opt R ). Then, there exist time instants

0 ≤ t 1 < t 2 < T opt such that 1. Q opt R (t) = 0 and s opt r (t) = ŝ⋆ r (s opt (t), γ(t)), where ŝ⋆ r (s, γ) satisfies ŝ⋆ r (s, γ) = arg max sr∈(0,s 1 ) µ(s r ) 1 + γ s 2 -s r s 1 -s 2 , ( 1
3. Q opt R (t) = QR and s opt r (t) = s opt 1 (t) a.e. t ∈ [t 2 , T opt ].
We notice the different qualitative behavior of the optimal control law in this case with respect to the parallel configurations. We perform numerical simulations in which we show the benefits of considering an enlarged control set, as well as introducing the recirculation. We show that the treatment time decreases with a larger recirculation upper limit.

In the three previous cases the study gives as a result an explicit formula for the optimal control problem as a feedback function that depends on the pollutant concentrations of each zone. These results allow decision makers to use simpler models to control their water treatment systems that approximate relatively well the more complex problem of inhomogeneity.

The results of Chapter 4 are part of an ongoing work.

Study of stochastic modeling of sequencing batch reactors

The second part of this thesis concerns the modelling and study of a stochastic model of sequencing batch reactor. This model is obtained as a limit of a sequence of continuous time Markov processes whose jump rates and jump laws are prescribed following the rules of replication and death of microorganisms, input and consumption of substrate, and input flow of water to the tank, that is, processes whose jump law depend only on the current state of the system. The obtained model is given by the controlled stochastic differential equation

             dx t = µ(s t ) - u t v t -β x t dt + γ x t v t dW t , x(0) = y, ds t = -µ(s t )x t + u t v t (s in -s t ) dt, s(0) = z, dv t = u t dt, v(0) = w, (1.36)
where u t ∈ [0, u max ] is a control process that satisfies a state constraint:

u t = 0 if v(t) ≥ v max , W = (W t ) t≥0 is a one-dimensional
= (y, z, w) ∈ (0, ∞) × [0, s in ] × [v min , v max ]
and admissible control u with respect to the brownian W , that is, progressively measurable with respect to the Brownian filtration, there exists a solution of equation (1.36) up to the extinction time τ E = inf{t ≥ 0 | x t = 0}. We prove the following result:

Proposition 1.7 (Chapter 5, Proposition 5.7) Let u = (u t ) t≥0 be an admissible control with respect to W = (W t ) t≥0 and X u = (x u , s u , v u ) the solution of (1.36). The probability that x u = (x u t ) t≥0 hits 0 at some time instant is positive.

The typical problem of depollution consists in reaching the target set

C := {(x, s, v) ∈ R 3 + | x > 0, s ≤ s out , v = v max },
where s out < s in is an acceptable level of pollution. With Proposition 1.7 in mind, the problem of depollution of water in minimal time is not well posed, because once extinction occurs the depollution process stops. We address the problem of the maximization of the probability of hitting the target set before extinction of biomass:

V (ξ) = max u P ξ [τ C ≤ τ E ] .
In this problem, none of the usual assumptions are satisfied, i.e, there is no Lipschitzianity of the coefficients of the equation, the diffusion coefficient is degenerate, the target set C and the extinction set

E := {(x, s, v) | x = 0, s ∈ [0, s in ], v ∈ [v min , v max ]} have intersecting closures,
the property of interior cone for the target set does not hold, the extinction set is absorbent, and the cost function is discontinuous. We split the study of this problem in two parts.

In the first part, we characterize the value function ṽ(ξ) for initial conditions ξ in the set V := (0, ∞) × [s out , s in ] × {v max } where the control is constrained to be null, and we prove that this function is continuous in V.

Later, we reformulate the original state-constrained problem by removing the constraint v ≤ v max and considering the cost in V as a final cost:

V (t, ξ) = max u∈Ut E t,ξ [✶ C (X τ u ) + ṽ(X τ u )✶ V (X τ u )] .
where

τ u := τ u C ∪ τ u V ∪ τ u E
is the first time that the process hits C ∪ V ∪ E. We conjecture a dynamic programming principle, and give the reasons why this is a non-standard and difficult problem. We present numerical simulations that show the behavior of the function ṽ(ξ) on the set V. This is an ongoing work. Section 2.4 is an additional material which has not been still submitted.

Introduction

Depollution of water resources is an important ecological issue in the fields of eutrophication's prevention and wastewater treatment. This problem has been often approached with the help of fed-batch or continuous bioreactors, receiving an important attention in the literature of the last 40 years [START_REF] D'ans | Optimal control of bacterial growth[END_REF][START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF][START_REF] Gajardo | Minimal time sequential batch reactors with bounded and impulse controls for one or more species[END_REF][START_REF] Hong | Optimal substrate feeding policy for fed batch fermentation with substrate and product inhibition kinetics[END_REF][START_REF] Moreno | Optimal time control of bioreactors for the wastewater treatment[END_REF]. Roughly speaking, in those bioremedation processes, biotic agents process the pollutant until its concentration decrease to acceptable levels.

In this work, we consider a natural water resource (for instance, a lake, a creek, or a water table) of volume v polluted with a toxic substrate. Our objective is to make the concentration of this substrate decrease, as fast as possible, to a prescribed value s. This process is driven by a continuous stirred bioreactor of volume v r . The bioreactor is fed from the resource with a flow rate Q, and its output returns to the resource with the same flow rate Q, after a separation of biomass and substrate due to a settler. The settler avoids the presence of excessive biomass used for treatment in the natural resource, which could result in undesirable sludge and possibly leads to an increase of eutrophication. We assume that during the entire treatment the volume v of the resource does not change. This volume is also assumed to be much larger than the bioreactor volume. As a first approach, one can assume that the toxic substrate is homogeneously distributed in the resource and is instantaneously mixing. This case was studied in Gajardo et al. [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF]. Therein, optimal strategies among constant controls and feedback controls were characterized. The last result is recalled in the next section. In this work we break the assumption of homogenous distribution of the pollutant. Our approach consists in splitting the lake into two zones, with volumes v 1 and v 2 and loaded with substrate concentrations s 1 and s 2 , respectively. Both zones are assumed to be homogeneous and instantaneously mixed. Moreover, those zones are connected one to each other via a fixed diffusion rate D > 0. The first zone, called the active zone, is directly connected to the reactor, while the second one, called dead zone, is only connected to the first zone (see Figure 3.1). This situation models a depollution process where one can clearly identify a (dead) zone of difficult access. This is typically due to the location of the bioreactor and to the large size of the water resource to be treated. A similar situation was studied in Gajardo et al. [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF], where a division of the water resource composed by two homogenous zones were also considered. However, both zones were considered active therein, which led to very different results. In this work we characterize the optimal strategy which solves the situation described above, and we compare it with the solution of the homogeneous case. We prove the somehow surprising result that both optimal strategies coincide.

The outline of the article is the following. In Section 2.2 we introduced our mathematical model. For this, we split this section into three subsections. In the first one, we present the dynamics of the bioprocess involved in our study. In the second one, we prove the first results on those dynamics. These results are needed in order to establish our minimal time optimal control problem in the third subsection. In Section 2.3, the Pontryagin maximum principle is applied to our optimal control problem. This allows to state the optimality conditions and the adjoint system. A fundamental analysis of the adjoint variables is carried out; at the end of this section we state our main result, which is the characterization of the optimal feedback strategy for our bioremediation problem. In Section 2.4 we study the effect of adding a pump to the system in order to increase the diffusion between the two zones, which leads to the study the behavior of the treatment time with respect to the diffusion parameter. Finally, in Section 2.5, some numerical simulations are presented.

Mathematical model 2.2.1 Description of the dynamics

We consider the usual chemostat model to describe the dynamics of the bioreactor [START_REF] Smith | The theory of the chemostat[END_REF]:

       ẋr = µ(s r )x r - Q v r x r , ṡr = -µ(s r )x r + Q v r (s 1 -s r ).
Here, s r and x r indicate the concentrations of substrate and biomass, respectively, and µ(•) denotes the growth function of x r .

The time evolution of the concentrations s i of the pollutants in the two zones (i = 1 for the active zone and i = 2 for the dead zone) are given by the equations

       ṡ1 = Q v 1 (s r -s 1 ) + D v 1 (s 2 -s 1 ), ṡ2 = D v 2 (s 1 -s 2 ), (2.1)
where the volumes v 1 and v 2 are assumed to be constant and much larger than v r (each), and D denotes the diffusion rate coefficient between the active and dead zones. The control variable is the flow rate Q of the pumps that take water from the active zone (with a pollutant concentration s 1 ) and reinject water to the active zone (with a smaller concentration s r ).

Throughout this article, we assume that µ(•) fulfills the following conditions:

Assumption 2.1 The function µ(•) is well defined and nonnegative on R + := [0, +∞), is continuously differentiable on R ++ := (0, +∞), is strictly increasing, concave and satisfies that µ(0) = 0.

Typical examples of functions satisfying Assumption 2.1 are:

• Linear growth function µ(s) = As, for all s, where A > 0.

• Monod's law growth function:

µ(s) = µ max s K + s , for all s, (2.2)
where µ max is called the maximum growth rate and K the half-saturation constant, both are positive.

Since we have supposed the reasonable fact that v 1 and v 2 are considerably larger than v r , we can suppose that our process has a slow-fast dynamic in which the slow variables are (s 1 , s 2 ) and the fast variables are (x r , s r ). This implies that the substrate and biomass concentrations present in the bioreactor can be approximated in the slow time scale by a quasi-stationary state

(x ⋆ r , s ⋆ r ) = (x ⋆ r (s 1 , s 2 , Q), s ⋆ r (s 1 , s 2 , Q))
satisfying the following relations:

µ(s ⋆ r ) = Q v r , x ⋆ r = s 1 -s ⋆ r . ( 2
ṡ1 = αµ(s ⋆ r )(s ⋆ r -s 1 ) + β 1 (s 2 -s 1 ), ṡ2 = β 2 (s 1 -s 2 ), (2.4)
where

α := v r v 1 , β 1 := D v 1 , and
β 2 := D v 2 .
Remark 2.1 Notice that condition s ⋆ r ∈ [0, s 1 ] has both a mathematical purpose as well as a physical meaning. Indeed, on the one hand, due to the adopted slow-fast approximation, this condition implies that the biomass concentration x ⋆ r = s 1 -s ⋆ r is nonnegative, which is required from a mathematical point of view. And, on the other hand, since the bioreactor is used here as a cleaning device, it should maintain its substrate concentration lower than the substrate concentration s 1 of the treated water resource.

We define the set of admissible controls as follows: Of course, we need to operate our bioprocess under that condition in order to work with the reduced dynamic (2.4). However, we will try to impose this condition from a mathematical point of view via the results of Section 2.2.2.

A = {s ⋆ r : [0, ∞) → R | s ⋆ r measurable}. ( 2
We end this section by recalling the characterization of optimal strategies for the homogeneous case obtained in Gajardo et al. [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF]. Notice that this case is simply a particular case of (2.1) when dynamic of v 2 is omitted (for simplicity, we say v 2 = 0) and D = 0. Proposition 2.2 Consider the homogeneous case (v 2 = 0 and D = 0) and the minimal time problem for which one drives s 1 to reach, as fast as possible, a prescribed level s. Then, the optimal feedback fulfills Q opt (s 1 ) = v r µ(s opt (s 1 )) with

s opt (s 1 ) ∈ arg max s∈(0,s 1 ) µ(s)(s 1 -s).
(2.6)

The main purpose of this chapter is to show that this feedback strategy is also optimal for our case with an active and a dead zone.

Preliminary results on the dynamics

Throughout this article we assume the following relation on the initial conditions:

Assumption 2.3 Denote z = (z 1 , z 2 ) ∈ R 2
+ the initial condition of the dynamic (2.4). It holds that 0 < z 1 ≤ z 2 , that is, the dead zone is more polluted than the active zone at the beginning of our depollution process.

Notice that this assumption can be made with out loss of generality. Indeed, if z 1 > z 2 , when we start the depollution process we should expect that, at some time t, we have s 1 ( t) ≤ s 2 ( t). Then, it is enough to redefine the starting time of the dynamics as t.

Proposition 2.4 For any admissible control s

⋆ r satisfying s ⋆ r (t) ∈ [0, s 1 (t)
] for all t, and any s 1 and s 2 generated by the dynamics (2.4), the set

I = {s 1 ≤ s 2 } is invariant.
Proof. It suffices to subtract the dynamics of s 2 and s 1 :

ṡ2 -ṡ1 = -(s 2 -s 1 )[β 1 + β 2 ] + αµ(s ⋆ r )(s 1 -s ⋆ r ).
If the set {s 1 = s 2 } is reached, then the difference s 2 -s 1 will be nondecreasing and, consecuently, s 2 remains greater than or equal to s 1 .

Proposition 2.5 For any admissible control s

⋆ r satisfying s ⋆ r (t) ∈ [0, s 1 (t)]
, for all t, and any s 1 and s 2 generated by the dynamics (2.4), there exists K (depending on z, v 1 and v 2 ) such that s i (t) < K for all t, for i = 1, 2. Moreover, the target set is attained from any initial condition z ∈ I.

Proof. Let us consider the dynamics of the total substrate mass present in the water resource. Its derivative satisfies:

d dt (v 1 s 1 + v 2 s 2 ) = v 1 (αµ(s ⋆ r )(s ⋆ r -s 1 ) + β 1 (s 2 -s 1 )) + v 2 β 2 (s 1 -s 2 ) = -v r µ(s ⋆ r )(s 1 -s ⋆ r ) So, since s 1 -s ⋆
r ≥ 0, the total substrate mass present in the resource is a nonincreasing function. Then, one has

v i s i (t) < v 1 s 1 (t) + v 2 s 2 (t) ≤ v 1 z 1 + v 2 z 2 , which implies that s i (t) < v 1 z 1 +v 2 z 2 v i := K i for i = 1, 2. We conclude by setting K := max{K 1 , K 2 }.
Now, in order to see that the target set can be attained from any initial condition, consider a constant control s ⋆ r (t) = s ⋆ to be determined. In this case, (2.4) is a linear system in s = (s 1 , s 2 ) that can be written as

ṡ(t) = A(s ⋆ )s(t) + b(s ⋆ ), s(0) = z, (2.7)
where

A(s ⋆ ) = -αµ(s ⋆ ) -β 1 β 1 β 2 -β 2 , b(s ⋆ ) = αµ(s ⋆ )s ⋆ 0.
From the analysis of the trace and the determinant of the matrix A(s ⋆ ), we conclude that this is an invertible matrix with two different real negative eigenvalues. Then, the solution of (2.7) converges to the equilibrium point s ∞ = -A(s ⋆ ) -1 b(s ⋆ ) = (s ⋆ , s ⋆ ), concluding that s 2 (t) converges to s ⋆ as t goes to infinity. We choose for the control the constant values s ⋆ = s/2 if s 1 (0) ≥ s or s ⋆ ∈ (0, s 1 (0)) if s 1 (0) < s. Notice that the set I ⋆ := {s 2 ≥ s 1 ≥ s ⋆ } is positively invariant. Indeed, for any t ≥ 0 such that s 2 (t) > s 1 (t) = s ⋆ , we have ṡ1 (t) > 0, and for any t ≥ 0 such that s 2 (t) = s 1 (t) > s ⋆ we can apply Proposition 2.4. Thus, the control

s ⋆ r (t) = s ⋆ satisfies s ⋆ r (t) ∈ [0, s 1 (t)
] for all t up to the final time, and there exists a finite time t † > 0 such that s 2 (t † ) = s, concluding that this control leads the system to the target set in finite time.

Minimal time optimal control problem

In this section we describe our problem as a minimal time optimal control problem. For this, first note that, thanks to Assumption 2.3 and Proposition 2.4, to know the time of depollution it is enough to measure the dead zone pollutant's concentration s 2 . Consequently, our target set can be defined as

T = {(s 1 , s 2 ) ∈ R 2 + | s 2 ≤ s} = R + × [0, s]. (2.8)
On the other hand, due to Proposition 2.5, we expect to deal with admissible bounded controls. Since we are interested in controls s ⋆ r (•) ∈ A that satisfy s ⋆ r (t) ∈ [0, s 1 (t)] for all t ≥ 0, our admissible control set will be restricted to

A U := {s ⋆ r ∈ A | s ⋆ r (t) ∈ [0, s 1 (t)]
, for all t ≥ 0}. Hence, in this work, we consider the following optimization problem:

(P) inf s ⋆ r ∈A U {T | s(0) = z, s(T ) ∈ T }.
In what follows will show that the optimal solution s ⋆ r of (P) is given by a feedback control that coincides with the solution of the problem with the homogeneous zone.

We have the following proposition: Proposition 2.6 For any initial condition z ∈ I there exists a unique solution for Problem (P).

Proof. The dynamics (2.4) are continuous with respect to (s, s ⋆ r ). Then, for fixed s ∈ R 2 + , the velocity set generated as the image of the compact and connected set [0, s 1 ] by dynamics (2.4) is a compact and connected set and, since the control variable s ⋆ r only acts in the first component, the velocity set is a segment and thus convex. This result, along with Propositions 2.4 and 2.5, assure that hypotheses of Filippov's theorem [17, Theorem 9.2.i] are satisfied, thus proving the existence of an optimal control for problem (P).

Application of Pontryagin maximum principle

The Hamiltonian of problem (P) is the mapping

H : R 2 + × R 2 × R + → R given by H(s, λ; s ⋆ r ) = -1 + λ 1 αµ(s ⋆ r )(s ⋆ r -s 1 ) + (-λ 1 β 1 + λ 2 β 2 )(s 1 -s 2 ).
According to the Pontryagin maximum principle [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF], if s ⋆ r (•) is a solution of (P) with corresponding optimal trajectory s(•) = (s 1 (•), s 2 (•)), then there exists an adjoint state λ(•) = (λ 1 (•), λ 2 (•)), with no both components equal to zero at the same time, such that

λ1 = (αµ(s ⋆ r ) + β 1 )λ 1 -β 2 λ 2 , λ2 = -β 1 λ 1 + β 2 λ 2 , (2.9)
and

H(s(t), λ(t); s ⋆ r (t)) = max sr∈[0,s 1 (t)]
H(s(t), λ(t); s r ).

(2.10)

With these dynamics we can associate the transversality conditions

λ 1 (T ) = 0 and λ 2 (T ) < 0, (2.11)
where T denotes the optimal time (this notation is consistent with the definition of problem (P)).

Remark 2.2 It is possible to write the maximum principle even though the control set depends on the state s 1 (t), at each time instant t ≥ 0. Indeed, it suffices to do a change of variable u(t) ∈ [0, 1] to have a control u(•) whose control set is compact and convex and whose associated s ⋆ r (t) = u(t)s 1 (t) is an admissible control, rewrite the dynamics accordingly, and restate the maximum principle. The result is equivalent to our approach.

Condition (2.10) states that, for a fixed time t, s ⋆ r (t) maximizes the function

φ(t, s ⋆ r ) = -λ 1 (t)µ(s ⋆ r )(s 1 (t) -s ⋆ r ), (2.12)
Concerning the behavior of the adjoint state λ 1 (•), we have the following proposition.

Proposition 2.7 The adjoint state λ(•) associated with an optimal solution of problem (P) satisfies that λ 1 (t) < 0 for all t ∈ [0, T ).

Proof. If one writes the adjoint equations (2.9) as λi = φ i (t, λ 1 , λ 2 ) (i = 1, 2), one can notice that the partial derivatives ∂ j φ i (i = j) are non-positive. From the theory of monotone dynamical systems (see for instance [START_REF] Smith | Monotone dynamical systems[END_REF]), the dynamics (2.9) are thus competitive or, equivalently, cooperative in backward time. As the transversality conditions (2.11) give λ i (T ) ≤ 0 (i = 1, 2), we deduce by the property of monotone dynamics that one should have λ i (t) ≤ 0 (i = 1, 2) for any t ≤ T . Now, if there exists t < T such that λ 1 (t) = 0, then one should have λ 2 (t) < 0. This implies λ1 (t) > 0 (because β 1 > 0), thus obtaining that the trajectory never returns to the third quadrant, leading to a contradiction with the trasversality condition λ 1 (T ) = 0.

Let us define, for s 1 > 0, the value ŝ⋆ r = ŝ⋆ r (s 1 ) as the solution s ⋆ r in (0, s 1 ) of the problem ŝ⋆ r (s 1 ) := arg max

s ⋆ r ∈(0,s 1 ) µ(s ⋆ r )(s 1 -s ⋆ r ), (2.13)
that is given by the unique solution of the equation

µ(s ⋆ r ) = µ ′ (s ⋆ r )(s 1 -s ⋆ r ). (2.14)
For each s 1 > 0 the solution of is unique. Indeed, if we define

β(s 1 , s ⋆ r ) := µ(s ⋆ r )(s 1 -s ⋆ r ), thanks to Assumption 2.1, for fixed s 1 > 0 the function s ⋆ r → β(s 1 , s ⋆ r ) is positive in the interval (0, s 1 ), β(s 1 , 0) = β(s 1 , s 1 ) = 0,
and its second derivative is strictly negative, concluding that it is strictly concave in (0, s 1 ).

We also define the function

γ(s 1 ) := max s ⋆ r ∈(0,s 1 ) µ(s ⋆ r )(s 1 -s ⋆ r ) = µ(ŝ ⋆ r (s 1 ))(s 1 -ŝ⋆ r (s 1 )). (2.15)
We are now in position to establish the main result of our article.

Proposition 2.8

The optimal flow rate Q ⋆ , which solves the depollution problem stated in Section 1, is given by

Q ⋆ (s 1 ) = v r µ(s ⋆ r (s 1 ))
, where s ⋆ r (•) is the unique solution of the optimal control problem (P). This solution is characterized by the feedback control defined in (2.3). In particular, the optimal feedback control s ⋆ r satisfies that s

⋆ r ∈ [0, s 1 (t)], for all t ∈ [0, T ].
Proof. Let φ be the function defined in (2.12). Since λ 1 (t) < 0, for all t ∈ [0, T ) (by Proposition 2.7), we notice s ⋆ r (t) maximizes the Hamiltonian with respect to s r , or equivalently, maximizes ψ(t, s r ) = µ(s r )(s 1 (t) -s r ) with respect to s r on [0, s 1 (t)]. Then, s ⋆ r has to have the form s ⋆ r (t) = ŝ⋆ r (s 1 (t)), or equivalently, s ⋆ r is characterized by the feedback rule (2.3). Moreover, at its maximum with respect to s ⋆ r , ψ satisfies Fermat's theorem on minima and maxima: ∂ ∂s ⋆ r ψ(t, s ⋆ r ) = 0, which exactly coincides with (2.14).

Remark 2.3 Proposition 2.8 shows us that the optimal feedback policy does not depend on s 2 . However, the optimal time does.

Remark 2.4 In Remark 2.2 it was mentioned that the change of control s ⋆ r (t) = u(t)s 1 (t) allows to deal with controls u(•) whose set is independent of the state variable s 1 . The complete expression of the adjoint equations is

λ1 = (αµ(us 1 ) + β 1 )λ 1 -β 2 λ 2 + λ 1 αu(µ ′ (us 1 )s 1 (1 -u) -µ(us 1 )), λ2 = -β 1 λ 1 + β 2 λ 2 , (2.16)
These equations, along with the maximum principle, allow to prove that the adjoint state λ 1 is strictly negative in the whole interval [0, T ) in the same manner than in Proposition 2.7 and, in particular, to conclude the equality (µ ′ (u(t)s 1 (t))s 1 (t)(1 -u(t)) -µ(u(t)s 1 (t))) = 0 along an optimal pair (s(•), u(•)). For the sake of interpretation of the variables, we continued working with s ⋆ r instead of u.

As it was pointed out in Gajardo et al. [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF], condition (2.14) gives explicit expressions for the optimal policy when some explicit growth functions µ(•) are chosen. For instance:

1. For any linear growth function, we have that

s ⋆ r (s 1 ) = s 1 2 .
2. For the Monod law described in (2.2) the unique maximum is computed as follows

s ⋆ r (s 1 ) = K 2 + Ks 1 -K.

The effect of recirculation

We consider now the previous problem, but now with the addition of two pumps that helps the mixing of the two zones, as a procedure to help with the homogeneization of the substrate concentration in the two zones (Figure 3.2). Suppose that the speed of this new pumps is

Q R ∈ [Q min R , Q max R ], with 0 ≤ Q min R ≤ Q max R . Q s 2 s 1 1 v 2 v r v s , x r
Q Q Q R Q R D Figure 2.2.
Connections between active zone, dead zone and bioreactor, with the recirculation pumps.

The time evolution of the concentrations s i of the pollutants in the two zones, after the reduction of the system, is given by the equations

       ṡ1 = αµ(s ⋆ r )(s ⋆ r -s 1 ) + β 1 + Q R v 1 (s 2 -s 1 ) ṡ2 = β 2 + Q R v 2 (s 1 -s 2 ), (2.17)
where we recall that α = v r /v 1 , β 1 = D/v 1 and β 2 = D/v 2 . In this case, the Hamiltonian of problem is

H(s, λ; s ⋆ r , Q R ) = -1 + λ 1 αµ(s ⋆ r )(s ⋆ r -s 1 ) + (-λ 1 β 1 + λ 2 β 2 )(s 1 -s 2 ) + Q R - λ 1 v 1 + λ 2 v 2 (s 1 -s 2 ). (2.18)
We say that a control

(s ⋆ r , Q R ) is admissible if it satisfies for all t ≥ 0 that s ⋆ r (t) ∈ [0, s 1 (t)] and Q R (t) ∈ [Q min R , Q max R ].
We consider the same problem of minimum-time depollution. Propositions 2.4 and 2.5 are still valid. The maximum principle states that if (s ⋆ r , Q R ) is a solution of the problem, then there exist an adjoint state λ(•) non-null at every time that solves the adjoint equations and transversality conditions

       λ1 = αµ(s ⋆ r )λ 1 -(D + Q R ) - λ 1 v 1 + λ 2 v 2 , λ 1 (T ) = 0, λ2 = (D + Q R ) - λ 1 v 1 + λ 2 v 2 , λ 2 (T ) < 0, (2.19)
and such that for each time t ∈ [0, T ] the Hamiltonian is maximized at the control value

(s ⋆ r , Q R (t)).
Notice that this is equivalent to the independent maximization of the functions φ(t, s r ) with respect to s r ∈ [0, s 1 (t)] and of the function

Q R → Q R - λ 1 v 1 + λ 2 v 2 (s 1 -s 2 ).
Proposition 2.9 Suppose that (s ⋆ r , Q R ) is an optimal control for the problem with recirculation. Then, for all t ∈ [0, T ], the control s ⋆ r is given by the feedback control defined in (2.3), and we have that

Q R (t) = Q max R .
Proof. The proof that s ⋆ r (t) = ŝ⋆ r (s 1 (t)) is the same as in Proposition 2.8; it relies in the facts that the Hamiltonian is decoupled in the variables s ⋆ r and Q R , and in the nonpositivity of the adjoint state λ 1 , whose proof is the same as in Proposition 2.7.

Let us define the function

η(t) = λ 1 (t) v 1 - λ 2 (t) v 2 .
(2.20)

Notice that the maximimization of the Hamiltonian along the optimal trajectory and the invariance of the set

I implies that Q R (t) maximizes the function Q R → Q R η(t), for all t ∈ [0, T ].
We have that the time derivative of η(•) is

η = λ 1 α 1 v 1 µ(s ⋆ r ) + (D + Q R ) 1 v 1 + 1 v 2 η,
with final value η(T ) = -λ 2 (T )/v 1 > 0. We prove that this function is positive for all t ∈ [0, T ). Suppose that for some time t † < T we have η(t † ) = 0. Then, we have that η(t † ) = α 1 µ(s ⋆ r (t † ))λ 1 (t † ) < 0, and thus η(t) ≤ 0 for all t ∈ [t † , T ]. This contradicts the fact that at the final time η(T ) > 0, concluding the statement.

Since η is positive along the optimal trajectory, then for all t ∈ [0, T ] the control value that maximizes

Q R → Q R η(t) is Q R = Q max R .
This last result states that if it is possible to add a recirculation to the process, it is better to use the maximum speed of pumping during all the process. It seems natural to eliminate Q R as a control parameter (since it will not change along the optimal trajectory), and to consider instead its maximum value Q max R as a modification of the diffusion term. We can rename D + Q max R simply as D (and, in consequence, return to the system (2.4)), and to study the behavior of the minimal-time function with respect to this parameter. It seems natural to ask if the best possible scenario is the homogeneous one, that is, the case with infinite diffusion.

Let us denote the minimal-time function by V D (•), indexed by the value of the diffusion parameter D:

V D (z) = inf s ⋆ r ∈A U {T ≥ 0 | s z,s ⋆ r ,D (T ) ∈ T },
where s z,s ⋆ r ,D (•) denotes the solution of with initial condition z = (z 1 , z 2 ) ∈ I and admissible control s ⋆ r . We have the following result.

Proposition 2.10 For any z ∈ I, the function D → V D (z) is decreasing.

Proof. From Proposition 2.8, we know that for every D > 0 and initial condition z ∈ I, the optimal control s ⋆ r (•) is given by s ⋆ r (t) = ŝ⋆ r (s z,s ⋆ r ,D 1 (t)) with ŝ⋆ r (•) the feedback given in (2.3) and s z,s ⋆ r ,D (•) the solution of (2.4) with initial condition z and control s ⋆ r . We split the proof in two parts: i. Notice that under the optimal feedback, s z,s ⋆ r ,D 2 (•) is decreasing; thus, we can define the time change

τ (t) = z 2 -s z,s ⋆ r ,D 2 (t).
The time τ ranges from 0 to τ f := z 2 -s (independent of D); from (2.4) we get the formula for the final time (depending on D)

T = T D = v 2 D τ f 0 dτ s 2 (τ, D) -s 1 (τ, D) . Now, s 1 (•) is solution of the nonautonomous equation in time τ d dτ s 1 = f D (τ, s 1 ) := - αv 2 D γ(s 1 ) z 2 -τ -s 1 + v 2 v 1 . (2.21)
Notice that for

D 1 < D 2 we have f D 1 (τ, s 1 ) < f D 2 (τ, s 1 ), which implies that for all τ ∈ [0, τ f ], s 1 (τ, D 1 ) < s 1 (τ, D 2 ). We conclude that s z,s ⋆ r ,D 1 (T ) is increasing with D. ii. Now, the variable m(t) = v 1 s z,s ⋆ r ,D 1 (t) + v 2 s z,s ⋆ r ,D 2 (t) is also decreasing in time, with ṁ(t) = -v r γ(s z,s ⋆ r ,D 1
(t)). We define the time change

θ(t) = m 0 -m(t), m 0 = v 1 z 1 + v 2 z 2 , which ranges from θ = 0 to θ D f := m 0 -(v 1 s z,s ⋆ r ,D 1 (T ) + v 2 s).
From the point i. we see that θ D f is decreasing with D. Again, from (2.4), we obtain an expression for the final time:

T D = 1 v r θ D f 0 dθ γ(s 1 (θ, D)) . (2.22) Now, s 1 (•) is solution of the nonautonomous equation in time θ d dθ s 1 = fD (θ, s 1 ) := - 1 v 1 + D v 1 v 2 v r m 0 -(v 1 + v 2 )s 1 -θ γ(s 1 ) . (2.23)
For D 1 < D 2 we have fD 1 (θ, s 1 ) < fD 2 (θ, s 1 ), which implies that for all θ ∈ [0,

θ D f ], s 1 (θ, D 1 ) < s 1 (θ, D 2 ).
The function s → γ(s) is increasing. Indeed, if s > s ′ , we have s -s r > s ′ -s r , and then µ(s r )(s -s r ) > µ(s r )(s ′ -s r ). Taking supremum on s r we conclude that s → γ(s) is increasing. This result proves that the integrand in (2.22) is decreasing with D and this, along with the fact that θ D f is decreasing with D, allows to conclude that the optimal time V D (z) = T D is a decreasing function of D.

Numerical Simulations

Simulations have been run for the Monod law

µ max = 1[h -1 ], K = 10[l -1 ], a diffusion rate D = 0.1[l • h -1 ], and volumes v 1 = v 2 = 100[l] and v r = 1[l].
We perform the simulations for three different values of initial substrate concentrations in the active zone (z 1 ) and in dead zone (z 2 ). For each of these initial conditions, we chose three different values of the prescribed final level of substrate; s = 1, 3 and 5[l -1 ]. For each of them, we compare the optimal time obtained by using feedback strategies (denoted by T feedback ; see (2.3)) with optimal time obtained by a constant strategy sought in the fixed interval s ⋆ r ∈ [0, 5] (denoted by T constant ). The results are shown in Tables 2 It can be seen that the gain obtained by feedback-type control strongly increases when the value of s decreases. This shows that the importance of using feedback controls becomes more relevant when we expect to depollute better our resource.

In Tables 2.4-2.6 we show the behavior of the value function or equivalently the optimal treatment time (obtained using the feedback control ŝ⋆ r ) with respect to the diffusion parameter D ∈ {0.1, 1, 10} and different initial conditions. We can clearly see that the increasing the diffusion parameter results in a much shorter treatment time, regardless of the target level s. This shows the benefits of adding a pump to the system to help the mixing between the two zones at maximum speed, resulting in a much shorter treatment time.

Conclusion

In this work, we have modeled the inhomogeneity of the pollutant concentration as two separated zones, one of them connected to the treating device, called the active zone, and the other one being depolluted by diffusion with the first zone, called dead zone. We obtained an optimal state-feedback control for the problem of minimal time depollution, and this feedback turns out to have the same form that the optimal control for the homogeneous case as stated in [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF] where there is no dead zone, depending only of the measurements on the active zone. This has some interesting features for practitioners, for instance it does not require the knowledge of the volumes of the two zones or the diffusion parameter.

We have also shown that the optimal treatment time is decreasing with the diffusion parameter. We show this behavior by means of numerical simulations.

Introduction

Today, the decontamination of water resources and reservoirs in natural environments (lakes, lagoons, etc.) and in industrial frameworks (basin, pools, etc.) is of prime importance. Due to the availability of drinking water becoming scarce on earth, efforts have to be made to reuse water and to preserve aquatic resources. To this end, biological treatment is a convenient way to extract organic or soluble matter from water. The basic principle is to use biotic agents (generally micro-organisms) that convert the pollutant until the concentration in the reservoir decreases to an acceptable level. Typically, the treatment is performed with the help of continuously stirred or fed-batch bioreactors. Numerous studies have been devoted to this subject over the past 40 years (see, for instance, [3,[START_REF] D'ans | Optimal control of bacterial growth[END_REF][START_REF] D'ans | Time-optimal control for a model of bacterial growth[END_REF][START_REF] Gajardo | Minimal time sequential batch reactors with bounded and impulse controls for one or more species[END_REF][START_REF] Van Impe | Optimal adaptive control of fedbatch fermentation processes[END_REF][START_REF] Van Impe | Optimal adaptive control of biotechnological processes[END_REF][START_REF] Kittisupakorn | Comparison of optimisation based control techniques for the control of a cstr[END_REF][START_REF] Moreno | Optimal time control of bioreactors for the wastewater treatment[END_REF][START_REF] Smets | Optimal adaptive control of (bio)chemical reactors: past, present and future[END_REF][START_REF] Smets | Optimal control of (bio-)chemical reactors: generic properties of time and space dependent optimization[END_REF][START_REF] Soukkou | Optimal control of a CSTR process[END_REF][START_REF] Srinivasan | Dynamic optimization of batch processes: I. characterization of the nominal solution[END_REF]).

The following main types of procedure are usually considered:

• The direct introduction of the biotic agents to the reservoir. This solution could lead to the eutrophication of the resource.

• The draining of the reservoir to a dedicated bioreactor and the filling back of the water after treatment. This solution attempts to eradicate various forms of life supported by the water resource, that cannot survive without water (such as fish, algae, etc.).

Alternatively, one can consider a side bioreactor that continuously treats the water pumped from the reservoir and that injects it back with the same flow rate so that the volume of the reservoir remains constant at all time. At the output of the bioreactor, a settler separates biomass from the water so that no biomass is introduced in the resource. Such an operating procedure is typically used for water purification of culture basins in aquaculture [START_REF] Crab | Nitrogen removal techniques in aquaculture for a sustainable production[END_REF][START_REF] Eding | Design and operation of nitrifying trickling filters in recirculating aquaculture: A review[END_REF][START_REF] Piedrahita | Reducing the potential environmental impact of tank aquaculture effluents through intensification and recirculation[END_REF].

The choice of the flow rate presents a trade-off between the speed at which the water is treated and the quality of decontamination. Recently, minimal-time control problems with simple spatial representations have been formulated and addressed [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF]. Under the assumption that the resource is perfectly mixed, an optimal state-feedback that depends on the characteristics of the micro-organisms and on the on-line measurement of the pollutant concentration has been derived. Later, an extension with a more realistic spatial representation was proposed in [START_REF] Gajardo | Bioremediation of natural water resources via optimal control techniques[END_REF] that considers two perfectly-mixed zones in the resource: an "active" zone, where the treatment of the pollutant is the most effective, and a more confined or "dead" zone that communicates with the active zone by diffusion of the pollutant. It has been shown that the optimal feedback obtained for the perfectly mixed case is also optimal when one applies it on the pollutant concentration in the active zone only. The fact that this controller does not require knowledge of the size of the dead zone or of the value of the diffusion parameter, neither of the online measurement of the pollutant in the dead zone, is a remarkable property. Nevertheless, the minimal time is impacted by the characteristics of the confinement. In the present work, we consider that the treatment of the water resource can be split into two zones i.e. the water is extracted from the resource at two different points (instead of one), and the treated water returns to the resource (with the same flows) at two different locations. A diffusion makes connection between the zones (see Fig. 3.1). Such a division into two patches can represent real situations such as:

• natural environments where water tables or lagoons are connected together by a small communication path (this modeling covers also the particular case of a null diffusion when one has to treat two independent volumes),

• resource hydrodynamics that reveal influence zones for each pumping devices, depending on the locations of the extraction and return points,

• accidental pollution as an homogeneous strain diffusing into the complementary part of the resource.

The control problem consists in choosing dynamically the total flow rate q and the flow distribution q 1 , q 2 between the two patches, with the objective of having both of them decontaminated in minimal time. Notice that a particular strategy consists in having all the time a flow distribution entirely with one zone, which amounts to the former problem with active and dead zones mentioned above. We study here the benefit of switching dynamically the treatment to the other patch or treating simultaneously both patches. The associated minimal-time problem is significantly more complex, because there are two controls and the velocity set is non-convex (this is shown in the next Section). Indeed, it is necessary to use different techniques to address the cases of non-null diffusion between the two zones and the limiting case of null diffusion between the two zones.

The paper is organized as follows. In the next section, definitions and assumptions are presented. In Section 3.3, properties of the optimization problem with relaxed controls and non-null diffusion are investigated. In Section 3.4, the optimal control strategy for the original problem with non-null diffusion is given and proven. In Section 3.5, we address the particular case of null diffusion and we provide explicit bounds on the minimal-time function. Finally, we show numerical computations that illustrate the theoretical results, and give concluding remarks.

Definitions and preliminaries

In what follows, we denote by R the set of real numbers, R + and R ⋆ + the sets of non-negative and positive real numbers respectively. Analogously, R -and R ⋆ -are the sets of non-positive and negative real numbers respectively. We set also

R 2 + = R + × R + and R 2 -= R -× R -.
The time evolution of the concentrations s i (i = 1, 2) of pollutants in the two patches are given by the equations

       ds 1 dt = q 1 v 1 (s r -s 1 ) + D v 1 (s 2 -s 1 ) , ds 2 dt = q 2 v 2 (s r -s 2 ) + D v 2 (s 1 -s 2 ) , (3.1)
where the volumes v i (i = 1, 2) are assumed to be constant and D denotes the diffusion coefficient of the pollutant between the two zones. The control variables are the flow rates q i of the pumps in each zone, which bring water with a low pollutant concentration s r from the bioreactor and remove water with a pollutant concentration s i from each zone i, with the same flow rates q i .

The concentration s r at the output of the bioreactor is linked to the total flow rate q = q 1 + q 2 by the usual chemostat model:

       ds r dt = -µ(s r )x r + q v r (s in -s r ) , dx r dt = µ(s r )x r - q v r x r , (3.2)
where x r is the biomass concentration, v r is the volume of the bioreactor and µ(•) is the specific growth rate of the bacteria (without a loss of generality we assume that the yield coefficient is equal to one). These equations describe the dynamics of a bacterial growth consuming a substrate that is constantly fed in a tank of constant volume (see for instance [START_REF] Smith | The theory of the chemostat[END_REF]). The input concentration s in is given here by the combination of the concentrations of the water extracted from the two zones:

s in = q 1 s 1 + q 2 s 2 q 1 + q 2 . (3.3)
We assume that the output of the bioreactor is filtered by a settler, that we assume to be perfect, so that the water that returns to the resource is biomass free (see [START_REF] Diehl | Control of an ideal activated sludge process in wastewater treatment via an ode-pde model[END_REF][START_REF] Diehl | A reduced-order ODE-PDE model for the activated sludge process in wastewater treatment: classification and stability of steady states[END_REF] for considerations of settler modeling and conditions that ensure the stability of the desired steady-state of the settler).

The target to be reached in the minimal time is defined by a threshold s > 0 of the pollutant concentrations, that is

T = s = (s 1 , s 2 ) ∈ R 2 + | max{s 1 , s 2 } ≤ s . (3.4)
In the paper, we shall denote t f as the first time at which a trajectory reaches the target (when it exists).

We make the usual assumptions on the growth function µ(•) in absence of inhibition.

Assumption 3.1 µ(•) is a C 1 increasing concave function defined on R + with µ(0) = 0.
Under this last assumption, we recall that under a constant s in , the dynamics (3.2) admit a unique positive equilibrium (s ⋆ r , x ⋆ r ) that is globally asymptotically stable on the domain R + × R ⋆ + provided that the condition q/v r ≤ µ(s in ) is satisfied (see, for instance, [START_REF] Smith | The theory of the chemostat[END_REF]). Then, s ⋆ r is defined as the unique solution of µ(s ⋆ r ) = q/v r and x ⋆ r = s in -s ⋆ r . Consequently, considering expression (3.3), the controls q 1 and q 2 are chosen such that

q 1 + q 2 ≤ v r µ q 1 s 1 + q 2 s 2 q 1 + q 2 . (3.5)
We assume that the resource to be treated is very large. This amounts to considering that the bioreactor is small compared to both zones of the resource. Assumption 3.2 v 1 and v 2 are large compared to v r .

Let us define α = q 1 /q, r = v 1 /(v 1 + v 2 ), d = D/v r , and ǫ = v r /(v 1 + v 2 ). Then, the coupled dynamics (3.1)-(3.2) with (3.3) can be written in the slow-fast form

                       ds r dt = -µ(s r )x r + q v r (αs 1 + (1 -α)s 2 -s r ) , dx r dt = µ(s r )x r - q v r x r , ds 1 dt = ǫ α r q v r (s r -s 1 ) + d r (s 2 -s 1 ) , ds 2 dt = ǫ 1 -α 1 -r q v r (s r -s 2 ) + d 1 -r (s 1 -s 2 ) . (3.6)
Provided that the initial conditions of the variables (s r , x r ) belong to R + × R ⋆ + , applying Tychonnof's Theorem (see for instance [START_REF] Khalil | Nonlinear Systems[END_REF]), the dynamics of the slow variables (s 1 , s 2 ) can be approached using the reduced dynamics

     ṡ1 = ds 1 dτ = α r µ(s ⋆ r )(s ⋆ r -s 1 ) + d r (s 2 -s 1 ) , ṡ2 = ds 2 dτ = 1 -α 1 -r µ(s ⋆ r )(s ⋆ r -s 2 ) + d 1 -r (s 1 -s 2 ) (3.7)
in the time scale τ = ǫt. In this formulation, the quasi-steady-state concentration s ⋆ r of the bioreactor can be considered as a control variable that takes values in [0, αs 1 +(1-α)s 2 ], which is equivalent to choosing q ∈ [0, v r µ(αs 1 + (1 -α)s 2 )] when Assumption 3.1 is satisfied. In the following, we shall consider the optimal control for the reduced dynamics only. Nevertheless, we give some properties of the optimal feedback for the reduced dynamics when applied to the un-reduced one, in Section 3.4 (Remark 3.2) and Appendix.

Notice that the control problem can be reformulated with the controls u = (α, s ⋆ r ) that belong to the state-dependent control set

U (s) = {(α, s ⋆ r ) | α ∈ [0, 1], s ⋆ r ∈ [0, αs 1 + (1 -α)s 2 ]} (3.8)
equivalently to controls q 1 and q 2 . In what follows, a measurable function u(•) such that u(t) ∈ U (s(t)) for all t is called an admissible control.

Lemma 3.3 The domain R 2 + is positively invariant by the dynamics (3.7) for any admissible controls u(•), and any trajectory is bounded. Furthermore, the target T is reachable in a finite time from any initial condition in R 2 + .

Proof. For s 1 = 0 and s 2 ≥ 0, one has ṡ1 ≥ 0. Similarly, one has ṡ2 ≥ 0 when s 1 ≥ 0 and s 2 = 0. By the uniqueness of the solutions of (3.7) for measurable controls u(•), we deduce that R 2 + is invariant. From equations (3.7), one can write

r ṡ1 + (1 -r) ṡ2 = µ(s ⋆ r )(s ⋆ r -(αs 1 + (1 -α)s 2 )
) ≤ 0 for any admissible controls. One then deduces

rs 1 (t) + (1 -r)s 2 (t) ≤ M 0 = rs 1 (0) + (1 -r)s 2 (0), ∀t ≥ 0,
which provides the boundedness of the trajectories.

Consider the feedback strategy

α = r , s ⋆ r = rs 1 + (1 -r)s 2 2 ,
and we write the dynamics of m = rs 1 + (1 -r)s 2 as follows:

ṁ = -µ m 2 m 2 < 0 , ∀m > 0 .
Then, from any initial condition in R 2 + , the solution m(t) tends to 0 when t tends to infinity. Therefore, m(•) reaches the set [0, min(r, 1 -r)s] in a finite time, which guarantees that s = (s 1 , s 2 ) belongs to T at that time.

For simplicity, we define the function

β(σ, s ⋆ r ) = µ(s ⋆ r )(σ -s ⋆ r ) (3.9)
so that the dynamics (3.7) can be written in the more compact form

ṡ = F (s, u) + dG(s) (3.10)
where F (•) and G(•) are defined as follows:

F (s, (α, s ⋆ r )) = -    α r β(s 1 , s ⋆ r ) 1 -α 1 -r β(s 2 , s ⋆ r )    , G(s) =    s 2 -s 1 r s 1 -s 2 1 -r    .
The dynamics can be equivalently expressed in terms of controls v = (α, ζ) that belong to the state-independent set V = [0, 1] 2 with the dynamics where the relaxed controls ũ = (u a , u b , p) = (α a , s ⋆ ra , α b , s ⋆ rb , p) belong to the set

ṡ = F (s, (α, ζ(αs 1 + (1 -α)s 2 ))) + dG(s) (3
Ũ (s) = U (s) 2 × [0, 1].
Thanks to Lemma 3.3, the hypotheses of Filippov's existence theorem [17, Theorem 9.2.i] are satisfied for the relaxed dynamics (that is, with relaxed controls) and thus, for any initial condition there exists a solution of the relaxed problem. In the next section, we show that the relaxed problem admits an optimal solution that is also a solution of the original (non-relaxed) problem.

Study of the relaxed problem

Throughout this section, we assume that the parameter d is positive. The particular case of d = 0 will be considered later in Section 3.5. Let us write the Hamiltonian of the relaxed problem

H(s, λ, (α a , s ⋆ ra , α b , s ⋆ rb , p)) = -1 + p Q(s, λ, (α a , s ⋆ ra )) + (1 -p) Q(s, λ, (α b , s ⋆ rb )) + d(s 2 -s 1 ) λ 1 r - λ 2 1 -r (3.14)
which is to be maximized w.r.t. (α a , s ⋆ ra , α b , s ⋆ rb , p) ∈ Ũ (s), where λ = (λ 1 , λ 2 ), and we have defined, for convenience, the function

Q(s, λ, (α, s ⋆ r )) = -α λ 1 r β(s 1 , s ⋆ r ) + (1 -α) λ 2 1 -r β(s 2 , s ⋆ r ) . (3.15)
The adjoint equations are

       λ1 = λ 1 p α a r µ(s ⋆ ra ) + (1 -p) α b r µ(s ⋆ rb ) + d r -λ 2 d 1 -r , λ2 = -λ 1 d r + λ 2 p 1 -α a 1 -r µ(s ⋆ ra ) + (1 -p) 1 -α b 1 -r µ(s ⋆ rb ) + d 1 -r , (3.16)
with the following transversality conditions

   s 1 (t f ) < s, s 2 (t f ) = s ⇒ λ 1 (t f ) = 0, λ 2 (t f ) < 0 , s 1 (t f ) = s, s 2 (t f ) < s ⇒ λ 1 (t f ) < 0, λ 2 (t f ) = 0 , s 1 (t f ) = s, s 2 (t f ) = s ⇒ λ 1 (t f ) ≤ 0, λ 2 (t f ) ≤ 0 with λ(t f ) = 0 . (3.17)
As usual, a triple (s(•), λ(•), ũ⋆ (•)) satisfying (3.12), (3.16), (4.7), and

H(s(t), λ(t), ũ⋆ (t)) = max ũ∈ Ũ (s(t)) H(s(t), λ(t), ũ) (3.18)
is called an admissible extremal.

Lemma 3.4 Along any admissible extremal, one has λ i (t) < 0 (i = 1, 2) for any t < t f .

Proof. If one writes the adjoint equations (3.16) as λi = φ i (t, λ 1 , λ 2 ) (i = 1, 2), one can notice that the partial derivatives ∂ j φ i (i = j) are non-positive. From the theory of monotone dynamical systems (see for instance [START_REF] Smith | Monotone dynamical systems[END_REF]), the dynamics (3.16) are thus competitive or, equivalently, cooperative in backward time. As the transversality conditions (4.7) gives λ i (t f ) ≤ 0 (i = 1, 2), we deduce by the property of monotone dynamics that one should have λ i (t) ≤ 0 (i = 1, 2) for any t ≤ t f . Moreover, λ = 0 is an equilibrium of (3.16) and λ(t f ) has to be different from 0 at any time t ≤ t f . Then, λ i (t) (i = 1, 2) cannot be simultaneously equal to zero. If there exists t < t f and i ∈ {1, 2} such that λ i (t) = 0, then one should have λ j (t) < 0 for j = i. However, d > 0 implies λi (t) > 0, thus obtaining a contradiction with λ i ≤ 0 for any time from t to the final time t f .

For the following, we consider the function

γ(σ) = max s ⋆ r ≥0 β(σ, s ⋆ r ), σ > 0 , (3.19)
which satisfies the following property:

Lemma 3.5 Under Assumption 3.1, for any σ > 0, there exists a unique ŝ⋆ r (σ) ∈ (0, σ) that realizes the maximum in (3.19). Furthermore, the function γ(•) is differentiable and increasing with γ

′ (σ) = µ(ŝ ⋆ r (σ)) . (3.20)
Proof. Consider the function ϕ : (σ, w) ∈ R + × [0, 1] → β(σ, wσ) and the partial function ϕ σ : w ∈ [0, 1] → ϕ(σ, w) for fixed σ > 0. Notice that ϕ σ (0) = ϕ σ (1) = 0 and that ϕ σ (w) > 0 for w ∈ (0, 1). Simple calculation gives ϕ ′′ σ (w) = µ ′′ (wσ)(1 -w)σ 3 -2µ ′ (wσ)σ 2 , which is negative. Therefore, ϕ σ (•) is a strictly concave function on [0, 1] and consequently admits a unique maximum w ⋆ σ on [0, 1]. We conclude that w ⋆ σ belongs to (0, 1) or, equivalently, that the maximum of s ⋆ r → β(σ, s ⋆ r ) is realized for a unique ŝ⋆ r (σ) = w ⋆ σ σ in (0, σ).

Furthermore, one has ϕ ′ σ (w) = σµ ′ (wσ)(σ -wσ) -σµ(wσ), and the necessary optimality condition ϕ ′ σ (w ⋆ σ ) = 0 gives the equality

µ(ŝ ⋆ r (σ)) = µ ′ (ŝ ⋆ r (σ))(σ -ŝ⋆ r (σ)) . (3.21)
Simple calculation shows that for each w ∈ [0, 1], the function σ → ϕ(σ, w) is convex. Because the maximizer w ⋆ σ of ϕ σ (•) is unique for any σ, one can apply the rules of differentiability of pointwise maxima (see, for instance, [19, Chap. 2.8]), which state that the function γ(σ) = max w∈[0,1] ϕ(σ, w) is differentiable with

γ ′ (σ) = ∂ϕ ∂σ (σ, w ⋆ σ σ) = w ⋆ σ µ ′ (w ⋆ σ σ)(σ -w ⋆ σ σ) + µ(w ⋆ σ σ)(1 -w ⋆ σ ) .
Equation (3.21) provides the simpler expression (3.20), which shows that γ(•) is increasing.

We now consider the variable

η = -λ 1 r γ(s 1 ) - -λ 2 1 -r γ(s 2 ) (3.22)
which will play the role of a switching function. Notice that this is not the usual switching function of problems with linear dynamics w.r.t. a scalar control because our problem has two controls α and s ⋆ r that cannot be separated, and the second control acts non-linearly in the dynamics.

Lemma 3.6 For fixed (s, λ) ∈ R 2 + × R 2 -, the pairs u ⋆ = (α, s ⋆ r ) ∈ U (s) that maximize the function Q(s, λ, •) are the following: i. when η > 0: u ⋆ = (1, ŝ⋆ r (s 1 )), ii. when η < 0: u ⋆ = (0, ŝ⋆ r (s 2 )), iii. when η = 0 and s 1 = s 2 : u ⋆ ∈ [0, 1] × {ŝ ⋆
r } where ŝ⋆ r = ŝ⋆ r (s 1 ) = ŝ⋆ r (s 2 ), iv. when η = 0 and s 1 = s 2 : u ⋆ = (1, ŝ⋆ r (s 1 )) or u ⋆ = (0, ŝ⋆ r (s 2 )).

Proof. When η > 0, one can write, using Lemma 3.5 and λ 1 , λ 2 < 0,

Q(s, λ, (1, ŝ⋆ r (s 1 ))) = -λ 1 r γ(s 1 ) > α -λ 1 r γ(s 1 ) + (1 -α) -λ 2 1 -r γ(s 2 ) , ∀α ∈ [0, 1) ≥ α -λ 1 r β(s 1 , s ⋆ r ) + (1 -α) -λ 2 1 -r β(s 2 , s ⋆ r ) , ∀α ∈ [0, 1), ∀s ⋆ r ∈ [0, αs 1 + (1 -α)s 2 ], ≥ Q(s, λ, (α, s ⋆ r )) , ∀α ∈ [0, 1), ∀s ⋆ r ∈ [0, αs 1 + (1 -α)s 2 ]
, and for α = 1, one has Q(s, λ, (1, ŝ⋆ r (s 1 ))) > Q(s, λ, (1, s ⋆ r )) , ∀s ⋆ r = ŝ⋆ r (s 1 ). Therefore, the maximum of Q(s, λ, •) is reached for the unique pair (α, s ⋆ r ) = (1, s ⋆ r (s 1 )).

Similarly, when η < 0, one can show that the unique maximum is (α, s ⋆ r ) = (0, s ⋆ r (s 2 )).

When η = 0, one has

-λ 1 r γ(s 1 ) = -λ 2 1 -r γ(s 2 ) > Q(s, λ, (α, s ⋆ r )) , ∀α ∈ [0, 1] , ∀s ⋆ r / ∈ {ŝ ⋆ r (s 1 ), ŝ⋆ r (s 2 )}.
If s 1 = s 2 , one necessarily has λ 1 /r = λ 2 /(1 -r) = 0, and thus,

Q(s, λ, (α, s ⋆ r )) = -λ 1 r β(s 1 , s ⋆ r ) < -λ 1 r γ(s 1 ) = Q(s, λ, (α, ŝ⋆ r (s 1 ))) , ∀s ⋆ 1 = ŝ⋆ r (s 1 ), for any α ∈ [0, 1]. The optimal s ⋆ r is necessarily equal to ŝ⋆ r (s 1 ) = ŝ⋆ r (s 2 ).
If s 1 = s 2 , one has s⋆ r (s 1 ) = ŝ⋆ r (s 2 ), and consequently, using Lemma 3.5 and the fact that λ 1 and λ 2 are both negative,

Q(s, λ, (α, ŝ⋆ r (s 1 ))) = α -λ 1 r γ(s 1 ) + (1 -α) -λ 2 1 -r β(s 2 , ŝ⋆ r (s 1 )) < -λ 1 r γ(s 1 ) , ∀α ∈ [0, 1) Q(s, λ, (α, ŝ⋆ r (s 2 ))) = α -λ 1 r β(s 1 , ŝ⋆ r (s 2 )) + (1 -α) -λ 2 1 -r γ(s 2 ) < -λ 2 1 -r γ(s 2 ) , ∀α ∈ (0, 1]
Then, (α, s ⋆ r ) = (1, ŝ⋆ r (s 1 )) and (α, s ⋆ r ) = (0, ŝ⋆ r (s 2 )) are the only two pairs that maximize Q(s, λ, •). Proposition 3.7 At almost any time, an optimal control ũ⋆ of the relaxed problem satisfies the following property:

1. when η = 0 or s 1 = s 2 , one has F (s, ũ⋆ ) = F (s, u ⋆ ), where u ⋆ is given by Lemma 3.6 i.-ii.-iii.

2. when η = 0 and s 1 = s 2 , one has

ũ⋆ ∈ {(1, ŝ⋆ r (s 1 )), (0, ŝ⋆ r (s 2 ))} × U (s) × {1} ∪ U (s) × {(1, ŝ⋆ r (s 1 )), (0, ŝ⋆ r (s 2 ))} × {0} ∪ {(1, ŝ⋆ r (s 1 ), 0, ŝ⋆ r (s 2 ))} × [0, 1] ∪ {(0, ŝ⋆ r (s 2 ), 1, ŝ⋆ r (s 1 ))} × [0, 1]. (3.23)
Proof. According to Pontryagin's Maximum Principle, an optimal control ũ = (u a , u b , p) has to maximize for a.e. time the Hamiltonian H given in (3.14) or, equivalently, the quantity

(u a , u b , p) → Q(s, λ, (u a , u b , p)) = pQ(s, λ, u a ) + (1 -p)Q(s, λ, u b )
where λ 1 and λ 2 are negative (from Lemma 3.4). Let us consider the maximization of the function Q(s, λ, •) characterized by Lemma 3.6.

In cases i and ii, the function Q(s, λ, •) admits a unique maximizer u ⋆ . Thus, Q(s, λ, •) is maximized for u a = u ⋆ with p = 1 independent of u b (or, symmetrically, for u b = u ⋆ with p = 0 independent of u a ) or for u a = u b = u ⋆ independent of p ∈ [0, 1]. In any case, one has F (s, ũ⋆ ) = F (s, u ⋆ ).

In case iii, the function Q(s, λ, •) is maximized for a unique value of s ⋆ r = ŝ⋆ r (s 1 ) = ŝ⋆ r (s 2 ) independent of α. Thus, Q(s, λ, •) is maximized when s ⋆ ra is equal to this value with p = 1 independent of u b (and, symmetrically, when s ⋆ rb is equal to this value with p = 0 independent of u a ) or when both s ⋆ ra and s ⋆ rb are equal to this value independent of α a , α b and p. In any case, one has F (s, ũ⋆ ) = F (s, u ⋆ ), where u ⋆ ∈ [0, 1] × {s ⋆ r }.

In case iv, the function Q(s, λ, •) admits two possible maximizers. Thus, Q(s, λ, •) is maximized when u a is equal to one of these maximizers with p = 1 independent of u b , when, symmetrically, u b is equal to one of these maximizers with p = 0 independent of u a , or when u a and u b are equal to the two different maximizers independent of p. All these cases appear in the set-membership (3.23).

Remark 3.1 In case 2 of Proposition 3.7, a relaxed control ũ⋆ with p ∈ (0, 1) can be approximated by a high-frequency switching between non-relaxed controls u = (1, ŝ⋆ r (s 1 )) and u = (0, ŝ⋆ r (s 2 )) (see the "chattering control" in [START_REF] Borisov | Theory of chattering control. Systems & Control: Foundations & Applications[END_REF]). In practice, such a high-frequency switching between the two pumps is not desired.

The following Lemma will be crucial later at several places.

Lemma 3.8 Along any extremal trajectory, one has at almost any time

η = d γ(s 1 ) r + γ(s 2 ) 1 -r λ 2 1 -r - λ 1 r + d λ 1 r 2 µ(ŝ ⋆ r (s 1 )) + λ 2 (1 -r) 2 µ(ŝ ⋆ r (s 2 )) (s 1 -s 2 ) . (3.24)
Proof. Let us write the time derivatives of the products λ 1 γ(s 1 ) and λ 2 γ(s 2 ) that appear in the expression of the function η using expressions (3.12), (3.16) and (3.20):

d dt [λ 1 γ(s 1 )] = λ 1 r δ 1 + dγ(s 1 ) λ 1 r - λ 2 1 -r + d λ 1 r µ(ŝ ⋆ r (s 1 ))(s 2 -s 1 )
where we put

δ 1 = pα a [µ(s ⋆ ra )γ(s 1 ) -µ(ŝ ⋆ r (s 1 ))β(s 1 , s ⋆ ra )] + (1 -p)α b [µ(s ⋆ rb )γ(s 1 ) -µ(ŝ ⋆ r (s 1 ))β(s 1 , s ⋆ rb )]
. One can easily check that for any optimal control ũ⋆ given by Proposition 3.7, one has δ 1 = 0. Similarly, one can write

d dt [λ 2 γ(s 2 )] = λ 2 1 -r δ 2 + dγ(s 2 ) λ 2 1 -r - λ 1 r + d λ 2 1 -r µ(ŝ ⋆ r (s 2 ))(s 1 -s 2 )
where

δ 2 = p(1-α a ) [µ(s ⋆ ra )γ(s 2 ) -µ(ŝ ⋆ r (s 2 ))β(s 2 , s ⋆ ra )]+(1-p)(1-α b ) [µ(s ⋆ rb )γ(s 2 ) -µ(ŝ ⋆ r (s 2 ))β(s 2 , s ⋆ rb )]
, with δ 2 = 0 for any optimal control ũ⋆ given by Proposition 3.7.

Then, one obtains the equality (3.24).

We now prove that the non-relaxed problem admits an optimal solution that is also optimal for the relaxed problem. Proposition 3.9 The optimal trajectories for the problem with the convexified dynamics (3.12) are admissible optimal trajectories for the original dynamics (3.10). Furthermore, the optimal control u ⋆ (•) satisfies the following property

s 1 (t) = s 2 (t) =⇒ u ⋆ (t) = (1, ŝ⋆ r (s 1 )) or u ⋆ (t) = (0, ŝ⋆ r (s 2 )), for a.e. t ∈ [0, t f ] .
Proof. We will prove that the set of times whereby the optimal relaxed strategy generates a velocity that belongs to the convexified velocity set but not to the original velocity set has Lebesgue measure zero. For this, consider s 1 > s 2 and η = 0. Because γ(•) is increasing (see Lemma 3.5), γ(s 1 ) > γ(s 2 ). Additionally, η = 0 implies that λ 1 /r > λ 2 /(1 -r). From equation (3.24) of Lemma 3.8, we deduce the inequality η < 0 (where λ 1 and λ 2 are negative by Lemma 3.4). Similarly, to consider s 2 > s 1 and η = 0 implies that η > 0. We conclude that case 2 of Proposition 3.7 can only occur at times in a set of null measure, from which the statement follows.

Now, because the optimal strategy of the convexified problem is (at almost any time) an admissible extremal for the original problem, and because the optimal time of the convexified problem is less than or equal to the optimal time of the original problem, the original problem has a solution, and it is characterized by point 1 of Proposition 3.7.

The last statement of the proposition follows from point 1 of Proposition 3.7.

Synthesis of the optimal strategy

According to Proposition 3.9, we can now consider optimal trajectories of the original (nonrelaxed) problem, knowing that the optimal strategy is "bang-bang" except on a possible singular arc that belongs to the diagonal set

∆ := {s ∈ R 2 + s.t. s 1 = s 2 }.
Proposition 3.10 For d > 0, the following feedback control drives any initial state in R 2 + \ T to the target T in minimal time:

u ⋆ [s] = (1, ŝ⋆ r (s 1 )) when s 1 > s 2 , (r, ŝ⋆ r (s 1 )) = (r, ŝ⋆ r (s 2 )) when s 1 = s 2 , (0, ŝ⋆ r (s 2 )) when s 1 < s 2 .
(3.25)

Proof. From Pontryagin's Maximum Principle, a necessary optimality condition for an admissible trajectory is the existence of a solution to the adjoint system

       λ1 = λ 1 α r µ(s ⋆ r ) + d λ 1 r - λ 2 1 -r , λ2 = λ 2 1 -α 1 -r µ(s ⋆ r ) + d λ 2 1 -r - λ 1 r , (3.26)
with the transversality conditions (4.7) and where u ⋆ = (α, s ⋆ r ) maximizes the Hamiltonian (3.24), one obtains the property s 1 > s 2 and η < 0 ⇒ η < 0 using the facts that λ i (i = 1, 2) are negative (Lemma 3.4) and that γ(•) is increasing (Lemma 3.5). When η < 0, one has u ⋆ = (0, ŝ⋆ r (s 2 )) from Lemma 3.6, and it is possible to write

H(s, λ, u) = -1 + Q(s, λ, u) + d(s 2 -s 1 ) λ 1 r - λ 2 1 -r w.r.t. u.

Consider the set

I -= (s, η) ∈ (R 2 + \ T ) × R s.t. s 1 > s 2 and η < 0 . From expression
ṡ1 -ṡ2 = - d r(1 -r) (s 1 -s 2 ) + γ(s 2 ) 1 -r ,
which shows that s 1 -s 2 remains positive for any future time. Thus, the set I -is positively invariant by the dynamics defined by systems (3.7) and 3.26). We deduce that the existence of a time t < t f such that (s(t), η(t)) ∈ I -implies (s(t f ), η(t f )) ∈ I -, and from the transversality condition (4.7), one obtains λ 1 (t f ) < λ 2 (t f ) = 0. Then, one should have η(t f ) = -λ 1 (t f )γ(s 1 (t f ))/r > 0, thus obtaining a contradiction. Similarly, one can show that the set

I + = (s, η) ∈ (R 2 + \ T ) × R s.t. s 1 < s 2
and η > 0 is positively invariant and that the transversality condition implies that (s, η) never belongs to I + along an optimal trajectory. Because ∆ is the only possible locus of a singular arc, we can form a conclusion about the optimality of (3.25) outside ∆.

Now, consider the function

L(s) = 1 2 (s 1 -s 2 ) 2
and write its time derivative along an admissible trajectory s(•) as follows:

L = ∇L, ṡ = - α r β(s 1 , s ⋆ r ) + 1 -α 1 -r β(s 2 , s ⋆ r ) (s 1 -s 2 ) - 2d r(1 -r) L .
Along an optimal trajectory, one has

L + 2d r(1 -r) L = - γ(s 1 ) r (s 1 -s 2 ) when s 1 > s 2 , γ(s 2 ) 1 -r (s 1 -s 2 ) when s 1 < s 2 ,
and deduces that the inequality L + 2d r(1-r) L ≤ 0 is satisfied. Consequently, the set ∆ ⊂ L -1 (0) is positively invariant by the optimal dynamics. On ∆, the maximization of Q(s, λ, •) gives the unique s * r = ŝ⋆ r (s 1 ) = ŝ⋆ r (s 2 ) because λ 1 , λ 2 are both negative (see Lemmas 3.4, 3.5 and 3.6). Finally, the only (non-relaxed) control that leaves ∆ invariant is such that α = r. Remark 3.2 The feedback (3.25) has been proved to be optimal for the reduced dynamics (3.7). In the Appendix, we prove that this feedback drives the state of the un-reduced dynamics (3.6) to the target in finite time, whatever is ǫ > 0. In Section 3.6, we show on numerical simulations how the time to reach the target is close from the minimal time of the reduced dynamics when ǫ is small.

Study of the minimal-time function

Define the function

T (σ) = max(0, T (σ)) with T (σ) = σ s dξ γ(ξ)
, σ > 0 .

Lemma 3.11 T (•) is strictly concave on [s, +∞).

Proof. Lemma 3.5 allows one to claim that T (•) is twice differentiable for any σ > 0 and that one has

T ′′ (σ) = - γ ′ (σ) γ(σ) 2 < 0 , ∀σ > 0 .
The function T (•) is strictly concave on R + , and because T (•) coincides with T (•) on [s, +∞), we conclude that T (•) is strictly concave on this interval.

Let us denote the minimal-time function by V d (•), indexed by the value of the parameter d:

V d (x) = inf u(•) {t > 0 | s(x, u, d, t) ∈ T } ,
where s(x, u, d, •) denotes the solution of (3.10) with the initial condition s(0) = x = (x 1 , x 2 ), the admissible control u(•) and the parameter value d. Lemma 3.3 ensures that these functions are well defined on R 2 + .

Proposition 3.12 The value functions V d (•) satisfy the following properties.

i. For any

d ≥ 0, V d (•) is Lipschitz continuous on R 2 + . ii. For d = 0, one has V 0 (x) = rT (x 1 ) + (1 -r)T (x 2 ) for any x ∈ R 2
+ , and the feedback (3.25) is optimal for both relaxed and non-relaxed problems.

Proof. On the boundary ∂ + T of the target that lies in the interior of the (positively) invariant domain R 2 + , the set N (•) of unitary external normals is

N (s) = 0 1 when s 1 < s and s 2 = s, cos θ sin θ θ∈[0,π/2]
when s 1 = s 2 = s, 1 0 when s 1 = s and s 2 < s .

At any s ∈ ∂ + T , one has inf u∈U (s) inf ν∈N (s) F (s, u) + dG(s), ν ≤ inf u∈U (s) inf ν∈N (s) F (s, u), ν = -γ(s) < 0 .
Furthermore, the maps

s → F (s, u) + dG(s)
are Lipschitz continuous w.r.t. s ∈ R 2 + uniformly in u. According to [5, Sect 1. and 4, Chap. IV], the target satisfies then the small time locally controllable property, and the value functions

V d (•) are Lipschitz continuous on R 2 + .
When d = 0, the feedback (3.25) provides the following dynamics

ṡ = -1 r γ(s 1 ) 0 when s 1 > max(s 2 , s) , - γ(s 1 ) γ(s 2 ) when s 1 = s 2 > s , -1 1-r 0 γ(s 2 )
when s 2 > max(s 1 , s) , and one can explicitly calculate the time to go to the target for any initial condition x ∈ R 2 + , which we denote as W 0 (x):

W 0 (x) = -r max(x 2 ,s) x 1 ds 1 γ(s 1 ) - s max(x 2 ,s) ds 1 γ(s 1 ) when x 1 ≥ x 2 , -(1 -r) max(x 1 ,s) x 2 ds 2 γ(s 2 ) - s max(x 1 ,s) ds 2 γ(s 2 ) when x 1 ≤ x 2 .
One can check that W 0 is Lipschitz continuous and that it can be written as W 0 (x) = rT (x 1 ) + (1 -r)T (x 2 ). We now show that W 0 is a viscosity solution of the Hamilton-Jacobi-Bellman equation associated to the relaxed problem

H(x, ∇W 0 (x)) = -1+ max (ua,u b ,p)∈ Ũ (x) pQ(x, -∇W 0 (x), u a )+(1-p)Q(x, -∇W 0 (x), u b ) = 0, x / ∈ T , (3.27)
(where Q is defined in (3.15)) with the boundary condition W 0 (x) = 0, x ∈ T .

(3.28)

Consider the C 1 functions W 0,1 (x) = rT (x 1 ) , W 0,2 (x) = (1 -r)T (x 2 ) and W 0 (x) = W 0,1 (x) + W 0,2 (x) defined on R 2 + . One has ∇W 0,1 (x) =   r γ(x 1 ) 0   and ∇W 0,2 (x) =    0 1 -r γ(x 2 )   
which are non-negative vectors. One can then use Lemma 3.6 to obtain the property

H(x, ∇W 0,1 (x)) = H(x, ∇W 0,2 (x)) = H(x, ∇W 0 (x)) = 0, x ∈ R 2 + ,
which shows that W 0,1 , W 0,2 and W 0 are solutions of (3.27) in the classical sense.

At x / ∈ T with x i = s (i = 1, 2), W 0 is C 1 and locally coincides with W 0 . Then, it satisfies equation (3.27) in the classical sense.

At x /

∈ T with x 1 = s or x 2 = s, W 0 is not differentiable but locally coincides with max(W 0 , W 0,2 ) or max(W 0 , W 0,1 ). From the properties of viscosity solutions (see, for instance, [5, Prop 2.1, Chap. II]), one must simply check that W 0 is a super-solution of (3.27). At such points, the Fréchet sub-differential of W 0 is

∂ -W 0 (x) = 0, r γ(s) × 1 -r γ(x 2 ) when x 1 = s , r γ(x 1 ) × 0, 1 -r γ(s) when x 2 = s .
Because any sub-gradient δ -∈ ∂ -W 0 (x) is a non-negative vector, one can again use Lemma 3.6 and obtain

H(x, δ -) = 0, ∀δ -∈ ∂ -W 0 (x) ,
which proves that W 0 is a viscosity solution of (3.27). Moreover, W 0 satisfies the boundary condition (3.28). Finally, we use the characterization of the minimal-time function as the unique viscosity solution of (3.27) in the class of Lipschitz continuous functions with boundary conditions (3.28) (see [5, Th. 2.6, Chap IV]) to conclude that W 0 is the value function of the relaxed problem. Because the time W 0 (x) to reach the target from an initial condition x / ∈ T is obtained with the non-relaxed control (3.25), we also deduce that V 0 and W 0 are equal. Remark 3.3 In the case d = 0, the control given by (3.25) is optimal but not the unique solution of the problem. Indeed, in Proposition 3.12, we proved that V 0 (•) is the unique viscosity solution to equation (3.27), where one of the possible maximizers of the Hamiltonian given in (3.27) is given by (3.25), but on the set (s, ∞) 2 \ ∆ there are more choices for u; for instance,

u ⋆ [s] = (1, ŝ⋆ r (s 1 )) when s 2 ≤ s < s 1 , (0, ŝ⋆ r (s 2 ))
when s 1 > s 2 > s , (r, ŝ⋆ r (s 1 )) = (r, ŝ⋆ r (s 2 )) when s 1 = s 2 , (1, ŝ⋆ r (s 1 )) when s < s 1 < s 2 , (0, ŝ⋆ r (s 2 )) when s 1 ≤ s < s 2 satisfies (3.27).

Proposition 3.13

The functions V d (•) satisfy the following properties:

i. V d (x) = T (x 1 ) = T (x 2 ) for any x ∈ ∆ and d ≥ 0, ii. V ∞ (x) = lim d→+∞ V d (x) = T (rx 1 + (1 -r)x 2 ) for any x ∈ R 2 + , and iii. d → V d (x) is increasing for any x ∈ (s, +∞) 2 \ ∆.
Proof. Consider an initial condition x in ∆ \ T . The optimal synthesis given in Proposition 3.10 shows that the set ∆ is invariant by the optimal flow and that the dynamics on ∆ are ṡi = -γ(s i ), i = 1, 2 independent of d. We then conclude that V d (x) = T (x i ) for i = 1, 2.

Consider d > 0 and x /

∈ ∆ ∪ T . Denote for simplicity s(•) as the solution s(x, u ⋆ , d, •) with the feedback control u ⋆ given in Proposition 3.10, and t f = V d (x). Define t ∆ as the first time t such that s(t) ∈ ∆ (here, we allow the solution s(•) to possibly enter the target T before reaching ∆).

From equation (3.10) with control (3.25), one can easily check that the following inequalities are satisfied

x 1 > x 2 ⇒ x 1 > s 1 (t) ≥ s 2 (t) > x 2 , ∀t ∈ [0, t ∆ ] , x 1 < x 2 ⇒ x 1 < s 1 (t) ≤ s 2 (t) < x 2 , ∀t ∈ [0, t ∆ ] .
Then, because the function γ(•) is increasing (Lemma 3.5), one can write, if the state s has not yet reached ∆,

- d r(1 -r) |s 1 -s 2 |-M + ≤ d dt |s 1 -s 2 |≤ - d r(1 -r) |s 1 -s 2 |-M - (3.29)
with M -= min(γ(x 2 )/r, γ(x 1 )/(1 -r)) and M + = max(γ(x 1 )/r, γ(x 2 )/(1 -r)). Then, we obtain an upper bound on the time

t ∆ t ∆ ≤ r(1 -r) d log 1 + d |x 1 -x 2 | M -r(1 -r) (3.30)
which tends to zero when d tends to infinity. From (3.29), we can also write

|x 1 -x 2 |-M + t ∆ ≤ d r(1 -r) t ∆ 0 |s 2 (τ ) -s 1 (τ )|dτ ≤ |x 1 -x 2 |-M -t ∆
and finally, one obtains from (3.10) the following bounds on s i (t ∆ ) (i = 1, 2):

rx 1 + (1 -r)x 2 -max(r, (1 -r))M + t ∆ ≤ s i (t ∆ ) ≤ rx 1 + (1 -r)x 2 -min(r, (1 -r))M -t ∆ .
(3.31) Therefore, s 1 (t ∆ ) = s 2 (t ∆ ) converges to rx 1 + (1 -r)x 2 when d tends to +∞. Furthermore, one has

t f = t ∆ + T (s(t ∆ )) when s(t ∆ ) / ∈ T , t f < t ∆ when s(t ∆ ) ∈ T . Because t ∆ → 0 and because T (•) is continuous with T (rx 1 + (1 -r)x 2 ) = 0 when rx 1 + (1 - r)x 2 ≤ s, we obtain the convergence V ∞ (x) = lim d→+∞ V d (x) = T (rx 1 + (1 -r)x 2 ) .

Now, consider the domain

D + = {s ∈ R 2 + | s 1 ≥ s 2 > s},
and let us show that any trajectory of the optimal flow leaves D + at (s, s) with the help of this simple argumentation on the boundaries of the domain:

s 2 = s ⇒ ṡ2 = d 1 -r (s 1 -s) ≥ 0 , s 1 = s 2 ⇒ ṡ1 = ṡ2 .
It is convenient to consider the variable s = rs 1 + (1 -r)s 2 , whose optimal dynamics in D + are simply

ṡ(t) = -γ(s 1 (t)) , t ∈ [0, t f ] . (3.32)
Because s(•) is strictly decreasing with time, an optimal trajectory in D + can be parameterized by the fictitious time

τ (t) = rx 1 + (1 -r)x 2 -s(t) , t ∈ [0, t f ] (3.33)
(where x is an initial condition in D + ). The variable s 1 (•) is then a solution of the scalar non-autonomous dynamics

ds 1 dτ = f d (τ, s 1 ) = - 1 r -d s 1 + τ -(rx 1 + (1 -r)x 2 ) r(1 -r)γ(s 1 ) when s 1 + τ > rx 1 + (1 -r)x 2 , -1 when s 1 + τ = rx 1 + (1 -r)x 2 ,
with the terminal fictitious time

τ f = rx 1 + (1 -r)x 2 -s .
Notice that τ f is independent of d. One then deduces the inequalities

d 1 > d 2 and s 1 + τ > rx 1 + (1 -r)x 2 =⇒ f d 1 (τ, s 1 ) < f d 2 (τ, s 1 )
and thus,

d 1 > d 2 and x ∈ D + \ ∆ =⇒ s 1 (x, u ⋆ , d 1 , τ ) < s 1 (x, u ⋆ , d 2 , τ ) , ∀τ ∈ [0, τ f ] . (3.34)
Finally, from equations (3.32) and (3.33), the time to reach the target can be expressed as

t f = τ f 0 dτ γ(s 1 (τ )) . (3.35)
Because the function γ(•) is increasing and because τ f is independent of d, one can conclude from (3.34) and (3.35) that

d 1 > d 2 and x ∈ D + \ ∆ =⇒ V d 1 (x) > V d 2 (x) .

The case of initial conditions in

D -\ ∆, with D -= {s ∈ R 2 + | s 2 ≥ s 1 > s}, is symmetric. Remark 3.4
The tightness V ∞ -V 0 of the bounds on the value function V d on (s, +∞) 2 \ ∆ is related to the concavity of the function T (•) on (s, +∞) (the less the concavity max σ∈[s,+∞) |T ′′ (σ)| is, the tighter the bounds are).

The bounds V 0 ≤ V d < V ∞ that are satisfied on the set (s, +∞) 2 \ ∆ are not necessarily satisfied outside this set: for x outside the target but such that rx 1 + (1 -r)x 2 < s, one has V ∞ (x) = 0 and V 0 (x) > 0. Therefore, we conclude that a large diffusion negatively impacts the time to treat the resource when both zones are initially polluted; however, when one of the two zones is initially under the pollution threshold, a large diffusion could positively impact the duration of the treatment.

Numerical illustrations

We consider the Monod (or Michaelis-Menten) growth function, which is quite popular in bioprocesses and which satisfies Assumption 3.1:

µ(s) = µ max s K s + s , with the parameters µ max = 1[h -1 ] and K s = 1[gl -1 ].
The corresponding function γ(•) is depicted in Fig. 3.2. The threshold that defines the target has been chosen as Several optimal trajectories in the phase portrait are drawn in Fig. 3.3 for small and large values of the parameter d. Finally, level sets of the value functions V 0 and V ∞ are represented in Fig. 3.4. One can make the following observations concerning the influence of the diffusion on the treatment duration, that we consider to be valuable from a practical viewpoint.

s = 1[gl -1 ]. s [g/L] µ 1 0.8 [1/h] 1 s [g/L] [g/L/h]
• When pollution is homogeneous, the best is to maintain it homogeneous, and the treatment time is then independent of the diffusion.

• A high diffusion is favorable for having fast treatments when initial concentrations are strongly different for the two zones. Typically, when the pollutant concentration is below the threshold in one patch, a high diffusion can reduce significantly the treatment time compared to a small diffusion.

• When initial concentrations in the two patches are close, a small diffusion leads to faster treatment than a large diffusion.

For various initial condition s(0), we have also performed numerical comparisons of the minimal time V d (s(0)) given by the feedback strategy (3.25) against two other non-optimal control strategies:

1. the best constant control that gives the smallest time T ⋆ cst to reach the target among constant controls, 2. the optimal one-pump feedback strategy obtained in the former work [START_REF] Gajardo | Bioremediation of natural water resources via optimal control techniques[END_REF]. This last control strategy considers that only one patch can be treated (that we called the "active zone"). The problem amounts then to consider the same dynamics (3.7) but one seeks the feedback s ⋆ r (•) that gives the minimal time T ⋆ one when α is imposed to be constantly equal to 1 (or 0 depending which patch is treated). In [START_REF] Gajardo | Bioremediation of natural water resources via optimal control techniques[END_REF], it has been proved that the feedback s 1 → ŝ⋆ r (s 1 ) is optimal. The results presented in Table 3.1 show first that the benefit of using the optimal feedback strategy over the other strategies increases with the level of initial pollution. The simulations also demonstrate the gain of using two pumps instead of one: for large concentrations of pollutant at initial time, one can see on the tables that a constant two-pumps strategy can be even better that the optimal feedback strategy restricted to the use of one pump only. This kind of situations typically occurs when diffusion is low and the time required by the optimal strategy for using simultaneously the two pumps is large compared to the overall duration. This is particularly noticeable when the initial pollution is homogeneous and the use of two pumps allows to maintain the levels of concentrations equal in both patches. We conclude that, for small diffusion, treating only one patch without the possibility to allocate the treatment in both patches could be quite penalizing. Figure 3.5 illustrates the time history of the two feedback controllers.

Furthermore, the Table 3.1 illustrates the effect of diffusion on the treatment times. One can first notice that the relative effect of the diffusion parameter d on the optimal time V d is decreasing with the threshold s. This can be explained by the fact that the proportion of the time spent on the set s 1 = s 2 , that is independent of the parameter d, is larger when one begins further away from the target. One can also see that the differences between strategies decrease when the diffusion increases. Intuitively, a high diffusion makes the resource behave quickly close to a perfectly mixed resource with one patch, leading consequently to less benefit of using more than one pump. Nevertheless, one can see that considering feedback controls remain quite efficient compared to constant ones when initial pollution is high.

Finally, we illustrate on Fig. 3.6 the effect of approximating the original dynamics (3.6) by the reduced one (3.7), when applying the feedback (3.25).

As proven in the Appendix, the feedback (3.25) drives the state to the target in finite time for any ǫ > 0.

Conclusion

In this work, we have shown that although the velocity set of the control problem is not convex, there exists an optimal solution with ordinary controls that is also optimal among relaxed controls. The optimal strategy consists in the most rapid approach to the homogenized concentration of pollutant in both patches. For the particular case of null diffusion, the most rapid approach path is not the unique solution of the problem. This optimal state-feedback has some Thus rs 1 + ǫs r is decreasing and has a limit when t tends to +∞. Since the trajectories are bounded, r ṡ1 + ǫ ṡr is uniformly continuous, and we conclude by Barbalat's Lemma (see for instance [START_REF] Khalil | Nonlinear Systems[END_REF]) that r ṡ1 +ǫ ṡr converges to 0, which implies that the positive quantities s 1 -s 2 and µ(s r )x r have to converge also to 0. Notice that s r = 0 implies ṡr = µ(s ⋆ r )s 1 > µ(ŝ ⋆ r (s))s > 0. So s r cannot tend to 0 and x r has necessarily to converge to 0. Write now the dynamics

d dt (s 1 -s r ) = -1 + ǫ r µ(s ⋆ r )(s 1 -s r ) -d ǫ r (s 1 -s 2 ) -µ(s r )x r ,
where µ(s ⋆ r ) > µ(ŝ ⋆ r (s)) > 0 and d ǫ r (s 1 -s 2 ) -µ(s r )x r tends to 0. Thus, there exists a time T > 0 large enough such that

s r (t) > s 1 (t) -η ≥ s ⋆ r (t) , ∀t > T ,
which implies to have µ(s r ) -µ(s ⋆ r ) > 0 for large t, thus a contradiction with the convergence of x r to 0.

Clearly the feedback (3.25) leaves the set {s 1 = s 2 } invariant. Denote for simplicity s l = s 1 = s 2 , and write ṡl + ǫ ṡr + ǫ ẋr = -ǫµ(s ⋆ r )x r < 0 . Trajectories are thus bounded, and by Barbalat's Lemma one obtains that µ(s ⋆ r )x r tends to 0. We prove now that s l has to reach s in finite time. If not, µ(s ⋆ r (t)) > µ(s) for any time and x r tends to zero. Write the dynamics

d dt (s l -s r ) = -(1 + ǫ)µ(s ⋆ r )(s l -s r ) + µ(s r )x r .
As before, we deduce that there exits a time T ′ > 0 such that

s r (t) > s l (t) -η ≥ s ⋆ r (t) , ∀t > T ′ ,
leading to a contradiction with the convergence of x r to 0.

Chapter 4

Minimal-time bioremediation of natural water resources with gradient of pollutant

This chapter corresponds to an ongoing work.

Introduction

The fight against eutrophication of lakes and natural reservoirs (that is, the excessive development of phytoplankton associated with an excess of nutrients) constitutes a major challenge. This ecological issue has given rise to many studies over the last 30 years (see, for instance, the surveys [START_REF] Gulati | Lakes in the Netherlands, their origin, eutrophication and restoration: state of the art review[END_REF] or [START_REF] Sondergaard | Lake restoration: successes, failures and long term effects[END_REF] and the references herein). To remediate eutrophication, various techniques such as bio-manipulation or ecological control have been proposed for mitigation. A common point across the proposed remediation approaches is that they are usually based on "biotic" actions on the lake trophic chain dedicated to the restoration of equilibrium in local ecosystems. To do so, most studies are based on empirical knowledge. However, since the 1970s, the use of eutrophication models (ranging from heuristic data-based models at a steady state to more recent dynamic mass-balance-based models) together with optimal control techniques have been proposed ( [START_REF] Estrada | Addressing the control problem of algae growth in water reservoirs with advanced dynamic optimization approaches[END_REF] and references herein).

Since the pioneering work by [START_REF] D'ans | Optimal control of bacterial growth[END_REF], the optimization of bioreactor operation has received great attention in the literature; see [START_REF] Alford | Bioprocess control: Advances and challenges[END_REF]3,[START_REF] Rani | Control of fermenters: a review[END_REF] for reviews of the different optimization techniques that have been used in bioprocesses. Among them, the theory of optimal control has proven to be a generic tool for deriving practical optimal rules [START_REF] Van Impe | Optimal adaptive control of biotechnological processes[END_REF][START_REF] Smets | Optimal adaptive control of (bio)chemical reactors: past, present and future[END_REF][START_REF] Smets | Optimal control of (bio-)chemical reactors: generic properties of time and space dependent optimization[END_REF]. The optimal control of continuous processes usually involves a two-step procedure. First, the optimal steady state is determined as a nominal set point that maximizes a criterion [START_REF] Soukkou | Optimal control of a CSTR process[END_REF][START_REF] Spitzer | Maximization of steady-state bacterial production in a chemostat with ph and substrate control[END_REF]. The benefit of operating a periodic control about the nominal point can be analyzed [START_REF] Abulesz | Periodic optimization of continuous microbial growth processes[END_REF][START_REF] Ruan | Comparison of several periodic oscillations of a continuous fermentation process[END_REF]. Then, a control strategy that drives the state about the nominal set point from any initial condition is searched for [START_REF] Kittisupakorn | Comparison of optimisation based control techniques for the control of a cstr[END_REF], possibly in the presence of model uncertainty using extremum seeking techniques [6,[START_REF] Marcos | Adaptive extremum seeking control of a continuous stirred tank bioreactor with Haldane's kinetics[END_REF][START_REF] Wang | Optimizing bioreactors by extremum seeking[END_REF][START_REF] Zhang | Adaptive extremum seeking control of continuous stirred-tank bioreactors[END_REF].

It is assumed that a small continuously stirred bioreactor is available to treat the polluted water by removing a substrate considered to be in excess. Particularly, we consider a natural water resource of volume v polluted with a substrate of concentration s l . As underlined above, typical examples of natural water resources in need of treatment are lakes or water tables that have been contaminated with diffused pollutants, such as organic matter or nutrients. The objective of the treatment is to make the concentration of the pollutant s l decrease as quickly as possible to a prescribed value s, with the help of a continuous stirred bioreactor of volume v r . The reactor is fed from the resource with a flow rate Q, and its output returns to the resource with the same flow rate Q after separation of biomass and substrate in a settler. The settler avoids the presence of excessive biomass used for treatment in the natural resource, which could result in undesirable sludge and possibly lead to an increase of eutrophication. We assume that during the entire treatment, the volume V of the resource does not change.

Q r v s , x r r biomass Q Q bioreactor settler s 2 1 v 2 v Q Q R s 1 Q R Figure 4.1.
Modeling scheme of the nonhomogeneity of concentrations in the resource as two zones, and a recirculation pump.

In [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF] the authors address the problem of the nonhomogeneity of the pollutant concentration in the resource, by means of a model of two compartments of volumes v 1 , v 2 such that v = v 1 + v 2 , connected in series to the bioreactor. The connection in series represents a gradient of concentration in the resource (see Figure 4.1). They assume the volumes to be large with respect to v r . Water is pumped from the first one while the clean one is rejected in the second one. One of the key assumption is that every admissible control law should satisfy that the pollutant concentration that is reinjected in the second zone must be smaller than the present pollutant concentration in said zone, in order to not repollute it. Under the previous constraint, the authors prove that the optimal control law consists on the maximization of a certain nonlinear function up to the moment when the substrate concentration in the reactor reaches the concentration in the second zone, moment from which is optimal to keep them even, and from the model dynamics it is concluded that in this stage the optimal control is constant.

In this paper we revisit the previously mentioned work, but considering two main changes in the model. The first change is that we drop the constraint of the output concentration of the reactor being smaller than the concentration in the second zone. The second change is that we add to the system a pump that introduces a recirculation of water from the first zone to the second zone, helping the system to mix. Two extreme cases are remarkable: if the speed of this pump is small, we should recover the results shown in [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF] for the series configuration; on the other hand, if the speed of this pump is too high, we should recover the homogeneous case (also studied in [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF]). We consider the speed of the pump Q R as a new control parameter. This paper is organized as follows. In Section 4.2, we introduce the mathematical model, and definitions and assumptions are presented. In Section 4.3 the minimal time optimal control problem is introduced and the optimal control is charaterized. Section 4.4 is devoted to numerical simulations and discussion.

Definitions and preliminaries

In what follows, we denote by R the set of real numbers, R + and R ⋆ + the sets of nonnegative and positive real numbers respectively. Analogously, R -and R ⋆ -are the sets of nonpositive and negative real numbers respectively.

For the time evolution of the concentrations of pollutant s i , i = 1, 2, in the two zones of the resource, by means of a mass balance we obtain the equations:

       ṡ1 = Q + Q R v 1 (s 2 -s 1 ), ṡ2 = Q v 2 (s r -s 2 ) + Q R v 2 (s 1 -s 2 ), (4.1)
where the volumes v i , i = 1, 2 are assumed to be constant, Q is the flow rate of the pump that connects the first zone with the bioreactor, and Q R is the flow rate of the pump that reinjects water from the zone 1 to the zone 2 (see Figure 4.1).

We consider the usual chemostat model to describe the dynamics of the bioreactor [START_REF] Smith | The theory of the chemostat[END_REF]:

       ṡr = -µ(s r )x r + Q v r (s 1 -s r ), ẋr = µ(s r )x r - Q v r x r (4.2)
where s r and x r indicate the concentrations of substrate and biomass inside the bioreactor, respectively.

The growth rate function µ(•) fulfills the following properties.

Assumption 4.1 The growth function µ(•) is defined on R + , continuously differentiable on R ⋆ + , increasing, concave, and satisfies µ(0) = 0.

Since we have supposed that v 1 and v 2 are considerably larger than v r , the process has a slowfast dynamic in which the slow variables are (s 1 , s 2 ) and the fast variables are (x r , s r ). This implies that the substrate and biomass concentrations present in the bioreactor can be approximated in the slow time scale by a quasi-stationary state

(x ⋆ r , s ⋆ r ) = (x ⋆ r (s 1 , s 2 , Q), s ⋆ r (s 1 , s 2 , Q))
The optimization problem consists of driving down the concentration of the resource to a prescribed value s > 0 from an initial condition z ∈ D in minimal time:

(P) min u∈U {t ≥ 0 | s(t) ∈ T , s(0) = z} .
The Hamiltonian of the problem is

H(s, λ; u) = -1 + µ(ws 1 ) [α 1 λ 1 (s 2 -s 1 ) + α 2 λ 2 (ws 1 -s 2 )] + Q R (s 1 -s 2 ) - λ 1 v 1 + λ 2 v 2 ,
that is maximized along the optimal trajectory; if

u opt = (w opt , Q opt R ) is a solution of problem (P), the adjoint state λ = (λ 1 , λ 2 ) satisfies                λ1 = λ 1 α 1 µ(w opt s 1 ) + Q opt R v 1 + α 1 w opt µ ′ (w opt s 1 )(s 1 -s 2 ) -λ 2 Q opt R v 2 + α 2 w opt µ(w opt s 1 ) + α 2 w opt µ ′ (w opt s 1 )(w opt s 1 -s 2 ) , λ2 = -λ 1 α 1 µ(w opt s 1 ) + Q opt R v 1 + λ 2 α 2 µ(w opt s 1 ) + Q opt R v 2 , (4.6)
with the transversality conditions

λ 1 (T opt ) < 0 , λ 2 (T opt ) = 0. (4.7)
Regarding the behavior of the adjoint states, we have the following lemma.

Lemma 4.3 Along any admissible extremal, one has λ 1 (t) < 0 for any t ≤ T opt , and λ 2 (t) < 0 for any t < T opt .

Proof. If one writes the adjoint equations (4.6) as λi = φ i (t, λ 1 , λ 2 ) (i = 1, 2), one can notice that the partial derivatives ∂ j φ i (i = j) are non-positive. From the theory of monotone dynamical systems (see for instance [START_REF] Smith | Monotone dynamical systems[END_REF]), the dynamics (4.6) is thus competitive or, equivalently, cooperative in backward time. As the transversality conditions (4.7) gives λ i (T opt ) ≤ 0 (i = 1, 2), we deduce by the property of monotone dynamics that one should have λ i (t) ≤ 0 (i = 1, 2) for any t ≤ T opt . Moreover, λ = 0 is an equilibrium of (4.6) and λ(T opt ) has to be different from 0 at any time t ≤ T opt . Then, λ i (t) (i = 1, 2) cannot be simultaneously equal to zero. If there exists t < T opt and i ∈ {1, 2} such that λ i (t) = 0, then one should have λ j (t) < 0 for j = i. However, µ(•) being positive and Q R ≥ 0 implies λi (t) ≥ 0; moreover, if λ 1 (t) = 0, then λ1 (t) ≥ 0 and we contradict the transversality condition λ 1 (T opt ) < 0; in the same way, if λ 2 (t) = 0, then λ2 (t) > 0 and we contradict the transversality condition λ 2 (T opt ) = 0.

Consider the following function

γ(t) = α 2 λ 2 (t) α 1 λ 1 (t) , (4.8)
whose time derivative γ fulfills

γ = µ(w opt s 1 ) ((α 2 -α 1 )γ -α 2 ) + Q opt R v 1 α 2 α 1 + γ (-1 + γ) + w opt α 1 γ -µ ′ (w opt s 1 )(s 1 -s 2 ) + γµ(w opt s 1 ) + γµ ′ (w opt s 1 )(w opt s 1 -s 2 ) , (4.9)
with final condition γ(T opt ) = 0. Also note that γ = 0 implies γ < 0 and therefore γ(t) > 0 for any t ∈ [0, T opt ). With this new variable, we can rewrite the Hamiltonian as

H(s, λ 1 , γ; u) = -1 -α 1 λ 1 (s 1 -s 2 )µ(ws 1 ) 1 + γ s 2 -ws 1 s 1 -s 2 + Q R (s 1 -s 2 ) -λ 1 v 1 [1 -γ] , (4.10
) For s 1 > s 2 > 0, w ∈ [0, 1], and γ > 0, we define the following functions

φ(s 1 , s 2 , γ, w) = µ(ws 1 ) 1 + γ s 2 -ws 1 s 1 -s 2 , (4.11)
ψ(s 1 , s 2 , γ) = µ ′ (s 1 )(1 -γ) -γ µ(s 1 ) s 1 -s 2 . (4.12) Notice that ψ(s 1 , s 2 , γ) = 1 s 1 ∂ ∂w φ(s 1 , s 2 , γ, w)| w=1 . When s 1 > s 2 and γ > 0, the max- imization of the Hamiltonian with respect to (w, Q R ) ∈ [0, 1] × [0, QR ] is equivalent to the maximization of the functions w → φ(s 1 , s 2 , γ, w) with w ∈ [0, 1] and Q R → (1 -γ)Q R with Q R ∈ [0, QR ].
Such maximum w opt is unique, as we see in Lemma 4.4. Lemma 4.4 Under Assumption 4.1, for s 1 > s 2 > 0 and γ > 0, the function φ(s 1 , s 2 , γ, •) is strictly concave on the interval 0, s 1 -s 2 γs 1 + s 2 s 1 .

Proof. Notice that the function φ can be written as

φ(s 1 , s 2 , γ, w) = γ s 1 -s 2 µ(ws 1 ) (c -ws 1 ) , c := s 1 -s 2 γ + s 2 . (4.13)
Define the function β(s, s r ) := µ(s r )(s -s r ), (s, s r ) ∈ R 2 + . This is a twice differentiable function, and for fixed s > 0, the function s r → β(s, s r ) is positive on the interval (0, s), negative on (s, ∞), and it is strictly concave in the interval [0, s]. Indeed,

∂ 2 ∂s 2 r β(s, s r ) = µ ′′ (s r )(s -s r ) -2µ ′ (s r ), (4.14)
which thanks to Assumption 4.1 in strictly negative. Taking s = c = s 1 -s 2 γ + s 2 , we see that β(c, •) is strictly concave in the set [0, c], and using that s r = ws 1 is a linear increasing function of w we conclude that the function w → β(c, ws 1 ) is strictly concave in [0, c/s 1 ].

The previous lemma allows to perform an analysis of the necesary conditions for the optimal control values. In the following lemma we characterize the possible values for an optimal control u ⋆ dependig on the functions ψ and γ.

Lemma 4.5 For fixed (s, λ 1 , γ) ∈ R 2 + × R -× R + , with s 1 > s 2 > 0 and γ > 0, the value u ⋆ = (w ⋆ , Q ⋆ R )
that maximizes the function u → H(s, λ 1 , γ; u) (with H as in (4.10)) are the following:

1. when γ ≥ 1, Q ⋆ R = 0 and w ⋆ is the unique solution of ∂ ∂w φ(s 1 , s 2 , γ, w ⋆ ) = 0, (4.15)
2. when γ < 1 and ψ(s, γ) < 0, Q ⋆ R = QR and w ⋆ is the solution of (4.15), 3. when γ < 1 and ψ(s, γ) ≥ 0, Q ⋆ R = QR and w ⋆ = 1.

Proof. Consider again the function β defined in Lemma 4.4. The maximization of φ(s 1 , s 2 , γ, •) is equivalent to the maximization of w → β(c, ws 1 ) in [0, 1], with c as in (4.13).

Consider first the case γ ≥ 1. This implies that s 1 ≥ c, and then w ⋆ belongs to the interval (0, c/s 1 ) (because in the interval (c/s 1 , 1] the function w → β(c, ws 1 ) is negative). In this interval the function w → β(c, ws 1 ) is differentiable and strictly concave (Lemma 4.4) and then the maximizer w ⋆ is unique and satisfies w ⋆ < 1. We conclude that w ⋆ satisfies (4.15). Now, we look at the case γ < 1. This implies that s 1 < c, and we cannot straightforwardly ensure that w ⋆ < 1; nevertheless, depending on the sign of ψ(s 1 , s 2 , γ) we can distinguish two situations: if ψ(s 1 , s 2 , γ) ≤ 0 we have that φ(s 1 , s 2 , γ, •) is decreasing in the interval [1, c/s 1 ] and consequently its maximum w ⋆ belongs to the interval (0, 1); in particular w ⋆ is unique and satisfies (4.15). On the other hand, if ψ(s 1 , s 2 , γ) > 0, that means that in the interval [0, 1] the function φ(s 1 , s 2 , γ, •) is increasing, being maximized in the extreme point w ⋆ = 1 (which is again unique).

The variable Q ⋆ R maximizes a linear function which only depends on the sign of 1 -γ (because s 1 > s 2 and by Lemma 4.3 we have λ 1 < 0). If γ > 1, then 1 -γ is negative and the maximum is the extreme point Q ⋆ R = 0; if γ < 1, then 1 -γ is positive and the maximum is the extreme point Q ⋆ R = QR . Notice that singular arcs with respect to Q R do not exist, because γ = 1 and γ = 0 imply s ⋆ r = 0, which can only happen if s 1 = 0.

Remark 4.1 The previous lemma can be expressed in terms of the equilibrium substrate concentration at the output of the bioreactor s r as

1. when γ ≥ 1, Q ⋆ R = 0 and s ⋆ r is the unique solution of ∂ ∂s r H(s 1 , s 2 , γ, s ⋆ r ) = 0, (4.16)
2. when γ < 1 and ψ(s, γ) < 0, Q ⋆ R = QR and s ⋆ r is the solution of (4.16), 3. when γ < 1 and ψ(s, γ) ≥ 0, Q ⋆ R = QR and s ⋆ r = s 1 .

In order to characterize the behavior of the optimal control u opt = (w opt , Q opt R ) along time, we need then to study the behavior of the function ψ as a function of time. The time derivative of ψ is

ψ = µ ′′ (s 1 )(1 -γ) -µ ′ (s 1 ) γ s 1 -s 2 ṡ1 -µ(s 1 ) d dt γ s 1 -s 2 -µ ′ (s 1 ) γ,
where a straightforward calculation provides

d dt γ s 1 -s 2 = α 2 µ(w opt s 1 ) s 1 -s 2 γ 1 - s 2 -w opt s 1 s 1 -s 2 -1 + Q opt R v 1 (s 1 -s 2 ) γ 2 + 2 α 2 α 1 γ - α 2 α 1 + γ α 1 w opt s 1 -s 2 -µ ′ (w opt s 1 )(s 1 -s 2 ) + γµ(w opt s 1 ) + γµ ′ (w opt s 1 )(w opt s 1 -s 2 ) .
The full expression for ψ can be written as ψ = I 1 + I 2 + I 3 + I 4 , where

I 1 = -µ ′′ (s 1 )(1 -γ) + µ ′ (s 1 ) γ s 1 -s 2 α 1 µ(w opt s 1 ) + Q opt R v 1 (s 1 -s 2 ), I 2 = -µ(s 1 ) α 2 µ(w opt s 1 ) s 1 -s 2 γ 1 - s 2 -w opt s 1 s 1 -s 2 -1 + Q opt R v 1 (s 1 -s 2 ) γ 2 + 2 α 2 α 1 γ - α 2 α 1 , I 3 = -µ ′ (s 1 ) µ(w opt s 1 ) ((α 2 -α 1 )γ -α 2 ) + Q opt R v 1 γ 2 + α 2 α 1 -1 γ - α 2 α 1 , I 4 = -γ α 1 w opt s 1 -s 2 -µ ′ (w opt s 1 )(s 1 -s 2 ) + γµ(w opt s 1 ) + γµ ′ (w opt s 1 )(w opt s 1 -s 2 ) × (µ(s 1 ) + µ ′ (s 1 )(s 1 -s 2 )). (4.17)
At the final time we have γ(T opt ) = 0 (thus Q opt R (T opt ) = QR ) and ψ(T opt ) = µ ′ (s) > 0 (thus w opt (T opt ) = 1); so there exists a time t ⋆ < T opt such that for all t ∈ (t ⋆ , T opt ], we have γ(t) < 1 and ψ(t) > 0. In terms of the optimal control u opt , by Lemma 4.5 this translates into Q opt R (t) = QR and w opt (t) = 1 (equivalently s opt r (t) = s 1 (t)) for all t in a time interval (t ⋆ , T opt ] before the final time.

Remark 4.2

The control Q opt R takes its maximum value before w opt can reach its maximum value. Indeed, if for some t > 0 we have γ(t) = 1, then ψ(t) < 0, so w opt (t) still belongs to the interval (0, 1) (Lemma 4.5). In term of the substrate concentration at the output of the reactor, s opt r (t) still belongs to the interval (0, s 1 (t)).

We have the two following lemmas concerning the behavior of the switching functions γ and ψ. Lemma 4.6 Suppose that at a certain time t † < T opt we have γ(t † ) = 1. Then, for all t ∈ (t † , T opt ] we have γ(t) < 1.

Proof. Suppose that γ = 1. Then, evaluating in the expression of γ given in (4.9), we have γ = -µ(s opt r )α 1 < 0. We conclude that once the set {γ ≤ 1} is reached, it is never left.

Lemma 4.7 Suppose that at a certain time t † < T opt we have ψ(t † ) = 0. Then, for all t ∈ (t † , T opt ] we have ψ(t) > 0.

Proof. Suppose ψ ≥ 0. Then γ belongs to the set [0, γ max (s)] where γ max (s) := µ ′ (s 1 )(s 1 -s 2 ) µ ′ (s 1 )(s 1 -s 2 )+µ(s 1 ) < 1. According to the control rule given in Lemma 4.5, w opt = 1 and Q opt R = QR . In this case, the formula for ψ is

ψ = -µ ′′ (s 1 )(1 -γ) + µ ′ (s 1 ) γ s 1 -s 2 α 1 µ(s 1 ) + QR v 1 (s 1 -s 2 ) -µ(s 1 ) α 2 µ(s 1 ) s 1 -s 2 (2γ -1) + µ ′ (s 1 ) ((α 2 -α 1 )γ -α 2 ) - QR v 1 µ(s 1 ) s 1 -s 2 γ 2 + 2 α 2 α 1 γ - α 2 α 1 + µ ′ (s 1 ) γ 2 + α 2 α 1 -1 γ - α 2 α 1 -γ α 1 s 1 -s 2 (µ ′ (s 1 )(s 1 -s 2 )(-1 + γ) + γµ(s 1 )) (µ(s 1 ) + µ ′ (s 1 )(s 1 -s 2 )) (4.18)
The previous expression corresponds to a concave function of second degree as function of γ; let us denote it by h s (γ). So, if we prove that h s is positive for γ = 0 and γ = γ max (s), we conclude that ψ ≥ 0. First, we take γ = 0, and compute

h s (0) = -µ ′′ (s 1 ) α 1 µ(s 1 ) + QR v 1 (s 1 -s 2 ) + α 2 α 1 µ(s 1 ) s 1 -s 2 + µ ′ (s 1 ) α 1 µ(s 1 ) + QR v 1 ,
which is positive. Now, we evaluate γ = γ max (s) in the formula of h s given by (4.18) (for ease we write γ max (s) as γ max ), which gives as result h s (γ max ) = I 1 + I 2 + I 3 + I 4 , with

I 1 = -µ ′′ (s 1 )(1 -γ max ) + µ ′ (s 1 ) γ max s 1 -s 2 α 1 µ(s 1 ) + QR v 1 (s 1 -s 2 ), I 2 = -µ(s 1 ) α 2 µ(s 1 ) s 1 -s 2 (2γ max -1) + µ ′ (s 1 ) ((α 2 -α 1 )γ max -α 2 ) , I 3 = - QR v 1 µ(s 1 ) s 1 -s 2 γ 2 max + 2 α 2 α 1 γ max - α 2 α 1 + µ ′ (s 1 ) γ 2 max + α 2 α 1 -1 γ max - α 2 α 1 I 4 = -γ max α 1 s 1 -s 2 (µ ′ (s 1 )(s 1 -s 2 )(-1 + γ max ) + γ max µ(s 1 )) (µ(s 1 ) + µ ′ (s 1 )(s 1 -s 2 )).
Since γ max < 1, we conclude that I 1 > 0. Now, for the analysis of I 2 , we notice that

2γ max -1 = 2 µ ′ (s 1 )(s 1 -s 2 ) µ ′ (s 1 )(s 1 -s 2 ) + µ(s 1 ) -1 = µ ′ (s 1 )(s 1 -s 2 ) -µ(s 1 ) µ ′ (s 1 )(s 1 -s 2 ) + µ(s 1 ) < 0,
the last inequality obtained from the fact that µ(•) is strictly concave and positive. We also have

(α 2 -α 1 )γ max -α 2 = α 2 (γ max -1) -α 1 < 0,
concluding that I 2 > 0. For I 3 we develop the quadratic expressions

γ 2 max + 2 α 2 α 1 γ max - α 2 α 1 = 1 α 1 (α 1 + α 2 )µ ′ (s 1 ) 2 (s 1 -s 2 ) 2 -α 2 µ(s 1 ) 2 (µ ′ (s 1 )(s 1 -s 2 ) + µ(s 1 )) 2 ,
and

γ 2 max + α 2 α 1 -1 γ max - α 2 α 1 = - µ(s 1 ) α 1 α 2 µ(s 1 ) + (α 1 + α 2 )µ ′ (s 1 )(s 1 -s 2 ) (µ ′ (s 1 )(s 1 -s 2 ) + µ(s 1 )) 2 ,
concluding that

I 3 = QR α 2 v 1 α 1 1 s 1 -s 2 µ(s 1 ) 2 µ ′ (s 1 )(s 1 -s 2 ) + µ(s 1 )
, is positive. Now, I 4 can be written as

I 4 = -γ max α 1 s 1 -s 2 (µ(s 1 ) + µ ′ (s 1 )(s 1 -s 2 ))J, with J = µ ′ (s 1 )(s 1 -s 2 )(-1 + γ max ) + γ max µ(s 1 ) = µ ′ (s 1 )(s 1 -s 2 ) -1 + µ ′ (s 1 )(s 1 -s 2 ) µ ′ (s 1 )(s 1 -s 2 ) + µ(s 1 ) + µ ′ (s 1 )(s 1 -s 2 ) µ ′ (s 1 )(s 1 -s 2 ) + µ(s 1 )
µ(s 1 ) = 0.

Thus, we conclude that if γ ∈ [0, γ max ], then ψ > 0.

Lemmas 4.6 and 4.7 summarize the behavior of any optimal control strategy u opt = (w opt , Q opt R ) in terms of the state variable s opt , the adjoint variable γ, and the function ψ. We state this result in the following proposition. Proposition 4.8 Suppose that s opt is an optimal trajectory, associated to the optimal control u opt = (w opt , Q opt R ). Then, there exist time instants 0 ≤ t 1 < t

and Q R → (1 -γ(t))Q R with Q R ∈ [0, QR ].
Suppose γ(0) > 1 (otherwise, t 1 = 0). Since γ is continuous with γ(T opt ) = 0, there exist a time instant t 1 ≥ 0 such that γ(t

Numerical simulations

We consider the Monod growth function, widely used in bioprocesses and satisfies Assumption 4.1

µ(s) = µ max s K s + s , (4.19)
with parameters µ max = 1[h -1 ] and K s = 1[gl -1 ], and volumes

v 1 = v 2 = 500[l], v r = 1[l].
The threshold level is set to s = 1[gl -1 ]. In Table 4.1 we compare the treatment time of the strategy presented in this work, that is, with the constraint s r ∈ [0, s 1 ] and without recirculation, with the strategy presented in [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF], with the constraint s r ∈ [0, s 2 ] and without recirculation. We see that the effect of allowing the second zone to repollute is benefical over the treatment time (in units of [h]), and its effect is even more notorious for initial conditions such that the second zone has a small concentration with respect to the threshold. This effect can be attributed to the fact that under the constraint s r ∈ [0, s 2 ] the variable s 2 can only be decreasing, providing a much more restricted control set than the constraint s r ∈ [0, s 1 ], notoriously for the initial conditions s 2 (0) close to 0. For a more mixed resource, this effect is not too notorious.

Initial condition Optimal time in [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF] T opt Improving

(2.5, 0. In Table 4.2 we show the effect of the recirculation in the time of treatment, with the optimal strategy described in Proposition 4.8. We compare the improving of the time with respect to the time obtained without recirculation. We can see that increasing the recirculation upper limit is benefical for the treatment time, and again the better results are obtained for initial conditions such that the pollution in the second zone is low compared to the threshold. 4.4, obtained by backward integration. We represent in blue the part of the trajectory in which there is no saturation of any of the controls, related to the point 1. of Proposition 4.8; in green we represent the part in which Q R takes its maximum value, related to the point 2. of Proposition 4.8; in red, Q R takes its maximum value and s r is equal to s 1 , related to the point 3. of Proposition 4.8.

Initial condition

T opt T opt T opt (z 1 , z 2 ) ( QR = 0) ( QR = 1) ( QR = 10) (2.

]

Phase portrait of optimal trajectories, Q R = 0

]

Phase portrait of optimal trajectories, Q R = 1

]

Phase portrait of optimal trajectories, Q R = 10

Conclusions

In this work we extend the results of [START_REF] Gajardo | Minimal time bioremediation of natural water resources[END_REF] by adding two changes to the proposed model. First, we consider that the concentration at the effluent of the bioreactor can be larger that the concentration in the second zone, allowing it to be repolluted; it is bounded by the concentration from where the bioreactor is fed. We prove that this enlargement of the control set is positive for the treatment, which is not an intuitive result. The second change is the introduction of a recirculation pump that helps the concentrations in the two zones to mix. We prove that every optimal control is composed of at most three steps, that involve the knowledge of an adjoint variable. In the first part of the treatment the recirculation pump is not active, up to a moment in which it is activated up to the final time. After this, at some moment the constraint s r ≤ s 1 becomes active, up to the final time. In this final period the pollution in zone 2 increases, which can be interpreted as the fact that an overdepollution of the second zone is performed to finish with a strong diffusion between both concentrations. For such a method to be appicable it is necessary to perform measurements at different places of the resource and of the initial speed of the variation in concentration at two remote locations in the resource for the identification of the inhomogeneity parameter of the model. Also, the control do not take the form of a feedback strategy, but with a proper estimation of the parameters of the model it is not difficult to perform the backwards integration to have the knowledge of the adjoint function and the switching times.

the input flow rate it is possible to control the time of treatment, in particular, the duration time of the fill mode. The tank has a maximum volume, and it is reasonable to establish as one of the goals to treat the maximum possible amount of wastewater in each cycle. This problem was studied in [START_REF] Moreno | Optimal time control of bioreactors for the wastewater treatment[END_REF] for the single-species case with monotonic and nonmonotonic kinetics using techniques in the plane such as Green's theorem. For the case of monotonic growth functions as the Monod law (see [START_REF] Smith | The theory of the chemostat[END_REF]) the optimal solution is the strategy called the most rapid approach strategy, which consists in filling the tank as fast as possible up to its maximum capacity and then wait until the target level of pollution is reached. For the case of nonmonotonic growth functions, the optimal solution consist in leading the system as fast as possible to a singular arc, that consists of maintaining the pollutant concentration in the level that maximizes the efficiency of the bacterial growth, and remain in the singular arc for as long as possible. In [START_REF] Gajardo | Minimal time sequential batch reactors with bounded and impulse controls for one or more species[END_REF] the authors study the problem of depollution in minimal time with several species of microorganism and monotonic and nonmonotonic growth functions and allowing impulse strategies. In the mentioned work the authors proved with techniques such as Pontryagin's maximum principle and Hamilton-Jacobi-Bellman equation that the previously mentioned strategies are optimal, but now with an impulse at the beginning of the process that takes the system instantaneously to a state in which the tank is full for the case of monotonic growth function, and in which the substrate concentration coincides with the maximum of the growth function in the case of nonmonotonic growth function.

Stochastic models for water treatment have been introduced for the continuously stirred bioreactors [START_REF] Campillo | Stochastic modeling of the chemostat[END_REF][START_REF] Collet | Stochastic models for a chemostat and long-time behavior[END_REF][START_REF] Imhof | Exclusion and persistence in deterministic and stochastic chemostat models[END_REF]. The models in these articles are obtained as an answer to the problem that the standard assumptions on the behavior of microbial population at the macroscopic time scale do not hold when considering a small number of individuals. Then, an analysis of the behavior of the population must be performed at a microscopic scale considering the demographic and environmental sources of randomness into the model. To the best or our knowledge, there is no literature for stochastic SBR models.

The outline of the present work is as follows. In Section 2 we develop a stochastic model of SBR as a limit of multidimensional birth-death processes. In Section 3 we study the existence of solutions of the controlled stochastic differential equation obtained. In Section 4 we discuss the problem of the maximization of the probability of reaching the target before extinction.

Stochastic SBR model

In what follows, we denote by R the set of real numbers, N and Z the sets of natural and integer numbers respectively, R + and R ⋆ + the sets of non-negative and positive real numbers respectively. For a set A and n ∈ N we denote A n := {x = (x 1 , . . . , x n )|x i ∈ A, i = 1, . . . , n}. For two real numbers a, b we define a ∧ b := min{a, b}. For A ⊆ R m , we denote by

C k (A; R n ) (C k b (A; R n ))
the set of k-times differentiable (bounded) functions f : A → R n with continuous k-th derivatives, and

C k (A) = C k (A; R) (C k b (A) = C k b (A; R)); D([0, ∞); R n
) denotes the set of right continuous functions f : [0, ∞) → R n with left limits everywhere. For a topological space (X, τ ) we denote by B(X) its Borel σ-algebra, and for a measurable space (S, S) we define P(S, S) (or simply P(S)) the set of probability measures on (S, S).

Denote by x, s, and v the biomass and substrate concentration, and culture volume inside the bioreactor vessel. We suppose that water with a constant concentration of nutrient s in and without bacteria is poured into the bioreactor at an inflow rate u ∈ [0, u max ]. We adopt the usual assumptions that the growth rate of microorganisms is proportional to the mass of microorganisms and it depends of the substrate concentration by means of the uptake function µ(•) [START_REF] Smith | The theory of the chemostat[END_REF]. We introduce an individual microbial death rate β ≥ 0. We suppose that the yield coefficient of the reaction, which is the amount of biomass produced by the bacterial specie when one unit of substrate is consumed by the reaction, is a constant Y > 0. We suppose that the culture volume is bounded between a minimum volume v min assumed to be the lower level to which the tank is emptied during the draw mode, and a maximum volume v max given by the maximum operative capacity of the tank. For the growth function µ(•) we make the following assumption:

Assumption 5.1 The growth function µ(•) is defined in [0, ∞), µ(0) = 0, is non-negative, bounded by above by a constant µ max > 0, and is Lipschitz continuous.

Our aim in this section is to develop a stochastic population process describing the dynamics of a bacterial specie, substrate, and volume inside the SBR. For this, we consider an individual-based tridimensional birth and death process that represents the discretized total number of microbial cells, substrate molecules, and the water molecules, scaled by scale parameters K x , K s , K v > 0 that represent the change of units from number of molecules into grams for a large population. We suppose that two of the following events cannot occur at the same time: the division of a microbial cell, the death of a cell, the entry of a unit of substrate into the tank (and at the same time the entry of unit of water), and the consumption of a unit of substrate by a cell. Let us define xK t the amount of cells composing the biomass scaled by the parameter K x , ŝK t the amount of substrate molecules scaled by the parameter K s , and vK t the amount of water molecules scaled by the parameter K v at the time instant t. At the same time, we introduce demographic randomness by introducing a perturbation parameter γ ≥ 0 in the birth and death rates of biomass, following [START_REF] Méléard | Quasi-stationary distributions and population processes[END_REF]. This procedure generates a pure jump Markov process η K t := (x K t , ŝK t , vK t ) that takes values in the state space D 3 K := (Z/K x ) × (Z/K s ) × (Z/K v ). Denote by ξ = (x, ŝ, v) a generic element in D 3 K and e i (i = 1, 2, 3) the components of the canonical basis of R 3 . The scheme of jumps of the process η K can be summarized as follows:

• ξ → ξ + e 1 K x with rate xK x µ ŝ v + γK x , • ξ → ξ - e 1 K x with rate xK x [β + γK x ], • ξ → ξ + e 2 K s + e 3 K v with rate uK v , • ξ → ξ - e 2 K s with rate xK x µ ŝ v .
This process is characterized by its infinitesimal generator L K defined as follows: for ξ = (x, ŝ, v) ∈ D 3 K and φ ∈ C (R 3 ),

L K φ( ξ) = lim hց0 E(φ(η K h ) | η K 0 = ξ) -φ( ξ) h ,
that in this case has the explicit formula

L K φ( ξ) := xK x µ ŝ v + γK x φ ξ + e 1 K x -φ( ξ) + xK K x [β + γK x ] φ ξ - e 1 K x -φ( ξ) + K v u φ ξ + e 2 K s + e 3 K v -φ( ξ) + K x xµ ŝ v φ ξ - e 2 K s -φ( ξ) , (5.1)
An explicit pathwise representation of η K as a semimartingale process can be obtained by applying L K to the projections p x ( ξ) = x, p s ( ξ) = ŝ and p v ( ξ) = v [START_REF] Joffe | Weak convergence of sequences of semimartingales with applications to multitype branching processes[END_REF]:

                 xK t = xK 0 + t 0 xK r µ ŝK r vK r -xK r β dr + M K, (x) t
,

ŝK t = ŝK 0 + t 0 K v K s u - K x K s xK r µ ŝK r vK r dr + M K, (s) t
,

vK t = vK 0 + t 0 udr + M K,(v) t , (5.2)
where M K = (M K,(x) , M K,(s) , M K,(v) ) is a discontinuous local martingale associated to the process η K whose quadratic variation terms can be explicitly computed:

M K,(x) t = t 0 2γ xK r dr + 1 K x t 0 xK r µ ŝK r vK r + β dr, M K,(s) t = 1 K s t 0 K v K s u + K x K s xK r µ ŝK r vK r dr, M K,(v) t = 1 K v t 0 udr, M K,(s,v) t = 1 K s t 0 udr, M K,(x,s) t = M K,(x,v) t = 0.
(5.3) Notice that K s /K x = 1/Kx 1/Ks denotes the amount of biomass generated by the consumption of one unit of substrate, which is by definition the yield coefficient Y ; in the same manner, K v /K s = 1/Ks 1/Kv denotes the increase of substrate due to the addition of one unit of water to the system, which is s in by definition. The scale parameters K s and K v can be expressed in terms of just one scale parameter K := K x and two parameters s in and Y as K x = K, K s = Y K, K v = Y s in K. The system (5.2) can be rewritten as

                 xK t = xK 0 + t 0 xK r µ ŝK r vK r -xK r β dr + M K,(x) t , ŝK t = ŝK 0 + t 0 us in - 1 Y xK r µ ŝK r vK r dr + M K,(s) t , vK t = vK 0 + t 0 udr + M K,(v) t , (5.4)
and we can replace the scale parameters in (5.3) accordingly. The only term that is not divided by the scale parameter K is the first term of M K,(x) . The model of bioreactor (5.4) is well defined, is non explosive, and when K tends to infinity, the limit process η := (x t , ŝt , vt ) t≥0 is a weak solution to the stochastic differential equation

             dx t = µ ŝt vt -β xt dt + 2γ xt dW t , dŝ t = - 1 Y µ ŝt vt xt + us in dt, dv t = udt, (5.5)
where (W t ) t≥0 is an standard Brownian motion in R. This result is stated in the next proposition:

Proposition 5.2 Suppose Assumption 5.1 holds. The sequence of processes (η K ) K∈N with infinitesimal generator (5.1) and deterministic initial condition η K 0 = ξ0 = (ŷ, ẑ, ŵ) is tight in P(D([0, ∞); R 3 )) and converges in distribution to a weak solution η of the stochastic differential equation (5.5) with the same initial condition η 0 = ξ0 .

Proof. The proof consists of various steps. First of all, the rate function λ K (•) and jump function ν K (•) are defined for ξ = (x, ŝ, v)

∈ D 3 K as                λ K x + ( ξ) = xK x µ ŝ v + γK x ; ν K x + ( ξ) = 1 Kx e 1 , λ K x -( ξ) = xK x [β + γK x ] ; ν K x -( ξ) = -1 Kx e 1 , λ K sv + ( ξ) = uK v ; ν K sv + ( ξ) = 1 Ks e 2 + 1 Kv e 3 , λ K s -( ξ) = xK x µ ŝ v ; ν K s -( ξ) = -1 Ks e 2 .
(5.6)

For any K > 0, the components of the process η K only take nonnegative values. Indeed, the underlying Markov chain takes values on D 3 K ∩ R 3 + , and it is not difficult to see for any point of D 3 K ∩ R 3 + with at least one component equal to zero, the probability of jumping to a state with a negative component is zero.

For any K > 0 the processes are well defined. This can be seen because for any point ξ = (x, ŝ, v), the jump rate λ K ( ξ) = θ λ K θ ( ξ) is positive, and

j λ K j ( ξ)|ν K j ( ξ)| 2 ≤ x2 + 1 K µ max 1 + 1 Y 2 + β 2 + u Y K s in + 1 s in .
So, according to [35, Corollary 2], η K t has finite moments of order two, and then it is well defined and it is non explosive. To prove the tension of the family (P K ) K of the laws of the processes η K on D([0, ∞); R 3 ) it is enough to check hypotheses (H 1 ), (H 2 ), and (H 3 ) of [START_REF] Joffe | Weak convergence of sequences of semimartingales with applications to multitype branching processes[END_REF]Section 3.2]. For this, define (according to notation of [START_REF] Joffe | Weak convergence of sequences of semimartingales with applications to multitype branching processes[END_REF]

) A K (t) = t, C = C b (R 3 ), L K (φ, ξ, t, ω) = L K φ( ξ), and b K ( ξ, t, ω) =   xµ(ŝ/v) -xβ us in -1 Y xµ(ŝ/v) u   and a K ( ξ, t, ω) =     2γ x + 1 K x µ ŝ v + β 0 0 0 1 Y K us in + x Y µ ŝ v u Y K 0 u Y K u Y s in K     .
• (H 1 ) The point i. of this hypothesis is satisfied. Indeed, we have

|b K ( ξ)| 2 = xµ ŝ v -xβ 2 + us in - 1 Y xµ ŝ v 2 + u 2 ≤ u 2 (1 + 2s 2 in ) + 2x 2 µ 2 max (1 + 1 Y 2 ) + β 2 , trace(a K ( ξ)) = 2γ x + 1 K x µ ŝ v + β + 1 Y K us in + 1 Y xµ ŝ v + 1 Y s in K u ≤ γ 2 + x2 1 + 1 2K + 1 2K µ max 1 + 1 Y 2 + β 2 + u Y K s in + 1 s in . Defining C K 1 = γ 2 + 1 2K µ max 1 + 1 Y 2 + β 2 + u Y K s in + 1 s in + u 2 (1 + 2s 2 in ), C K 2 = 1 + 1 2K + 2 µ 2 max (1 + 1 Y 2 ) + β 2 , then |b K ( ξ)| 2 + trace(a K ( ξ)) ≤ C K 1 + x2 C K 2 ≤ C K 2 C K 1 C K 2 + | ξ| 2 .
We can find uniform bounds on

C K 2 y D K := C K 1 /C K 2 (independent of K): C K 2 ≤ 2 1 + µ 2 max (1 + 1 Y 2 ) + β 2 := C, and
D K ≤ γ 2 + µ max 1 + 1 Y 2 + β 2 + u Y s in + 1 s in + u 2 (1 + 2s 2 in ) 1 + 2 µ 2 max (1 + 1 Y 2 ) + β 2 := D.
Finally, for all K ∈ N,

|b K ( ξ)| 2 + trace(a K ( ξ)) ≤ C D + | ξ| 2 ,
where D is a deterministic constant independent of t, which proves point ii. of this hypothesis.

• (H 2 ) The initial condition of each process η K is deterministic and independent of K, so (H 2 ) is satisfied.

• (H 3 ) Taking α(t) = t y ρ K = 0, (H 3 ) is satisfied.

Thus, for each T > 0 there exists a constant C T and K 0 such that for all K ≥ K 0 ,

E K sup t≤T |η K t | 2 ≤ C T (1 + E K |η K 0 | 2 ), (5.7)
(from [46, Lemma 3.2.2]), and the family (η K ) K∈N is tight (from [46, Proposition 3.2.3]), and then there exists a limit probability law P on D([0, ∞); R 3 + ). Also, for any T > 0,

E sup t≤T |η K t -η K t -|≤ 1/K holds, concluding that P is supported on C ([0, ∞); R 3 + ).
It remains to show that P solves the martingale problem associated to equation (5.5) (for the martingale problem, see [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF]).

Consider Ω = C ([0, ∞); R 3 + ), the σ-field F = B(Ω), and the system (Ω, F , (F t ) t≥0 , P ), where

F t = σ ω(s), 0 ≤ s ≤ t, ω ∈ C ([0, ∞); R 3 + ) . We define for φ ∈ C 2 b (R 3 ) y ξ = (x, ŝ, v), L φ( ξ) = ∂ ∂x φ( ξ) µ ŝ v -β x+ ∂ ∂s φ( ξ) - 1 Y µ ŝ v x + us in + ∂ ∂v φ( ξ)u+ 1 2 ∂ 2 ∂x 2 φ( ξ)2γ x;
(5.8) for ω ∈ Ω, t ≤ 0, the function M t φ : Ω → R 3 + defined by

M t φ(ω) = φ(ω(t)) -φ(ω(0)) - t 0 L φ(ω(r))dr,
and the canonical process on Ω given by η t (ω) = ω(t). We have to prove that ∀φ ∈ C 2 b (R 3 ), the process (M t φ(η)) t≥0 is an (F t )-local martingale under P . From (5.1) and (5.8),

|L K φ( ξ) -L φ( ξ)| ≤ γK 2 x φ ξ + e 1 K + φ ξ - e 1 K -2φ( ξ) - 1 K 2 ∂ 2 ∂x 2 φ( ξ) +K xµ ŝ v φ ξ + e 1 K -φ( ξ) - 1 K ∂ ∂x φ( ξ) +Kβ x φ ξ - e 1 K -φ( ξ) + 1 K ∂ ∂x φ( ξ) +Y s in Ku φ ξ + e 2 Y K + e 3 Y s in K -φ( ξ) - 1 Y K ∂ ∂s φ( ξ) - 1 Y Ks in ∂ ∂v φ( ξ) +K xµ ŝ v φ ξ - e 2 Y K -φ( ξ) + 1 Y K ∂ ∂s φ( ξ) .
Considering Taylor expansions of first and second order on φ in the previous formula,

|L K φ( ξ) -L φ( ξ)| ≤ γK 2 xo 1 K 2 + K xµ max o 1 K + Kβ xo 1 K + uY s in Ko 1 K + K xµ max o 1 K ≤ C(1 + | ξ| 2 )O 1 K . Define for φ ∈ C 2 b (R 3 ), y ∈ D([0, ∞); R 3 ), n ∈ N, 0 < s 1 ≤ s 2 ≤ . . . ≤ s n ≤ s ≤ t, y F s 1 , . . . , F sn ∈ C b , M K t φ(y) =φ(y(t)) -φ(y(0)) - t 0 L K φ(y(s))ds, G K t,s φ(y) =M K t φ(y) -M K s φ(y), G t,s φ(y) =M t φ(y) -M s φ(y), H n (y) =F 1 (y(s 1 )) . . . F n (y(s n )).
We have to prove that under P , (M t φ(η)) t≥0 is a martingale, this is,

E (G t,s φ(η)F 1 (η s 1 ) . . . F n (η sn )) = 0.
Notice that for fixed K,

E K G t,s φ(η K )H n (η K ) = E K [G t,s φ(η K ) -G K t,s φ(η K )]H n (η K ) +E K G K t,s φ(η K )H n (η K ) , (5.9
) where the second term of the right-hand side of (5.9) is zero because (M K t φ(η K )) t≥0 is a martingale under P K . It is easy to see that for 0 < s < t,

G t,s φ(η K ) -G K t,s φ(η K ) = [M t φ(η K ) -M s φ(η K )] -[M K t φ(η K ) -M K s φ(η K )] = t s L K φ(η K r )dr - t s L φ(η K r )dr.
Then,

E K |G t,s φ(η K ) -G K t,s φ(η K )|H n (η K ) = E K t s |L K φ(η K r ) -L φ(η K r )|dr H n (η K ) ≤ O 1 K E K t s (1 + |η K r | 2 )drH n (η K ) .
Since η K t has bounded moments of second order, uniformly bounded on finite intervals (from (5.10)), the previous expression converges to zero when K → ∞. On the other side, since for any t the function y(•) → L φ(y(t)) is continuous, the weak convergence of P n to P allows to say

E K G t,s φ(η K )F 1 (η K s 1 ) . . . F n (η K sn ) K→∞ -→ E (G t,s φ(η)F 1 (η s 1 ) . . . F n (η sn )) .
Then, taking limits when K → ∞ in equation (5.9), we conclude that

E (G t,s φ(η)H n (η)) = 0, (5.10)
which means that L satisfies the martingale problem under the law P . Now, with the previous result it is possible to prove that the process η is a solution of (5.5). The process η t = (x t , ŝt , vt ) has the integral representation

                 xt = x0 + t 0 µ ŝr vr xr -β xr dr + M (x) t , ŝt = ŝ0 + t 0 - 1 Y µ ŝr vr xr + us in dr + M (s) t , vt = v0 + t 0 udr + M (v) t ,
where M (x) , M (s) , M (v) are local martingales with quadratic variation

M (x) t = t 0 2γ xr dr, M (s) t = M (v) t = 0, M (i) , M (j) t = 0 i = j.
Consider (Ω ′ , F ′ , (F ′ t ) t≥0 , P ′ ) a probability space on which is defined an (F ′ t ) t≥0 -Brownian motion W ′ independent of M (x) . Define Ω * = Ω × Ω ′ , F * = F ⊗ F ′ , P * = P ⊗ P ′ , and on (Ω * , F * , P * ), the process

W t (ω, ω ′ ) = t 0 1 √ 2γ xr ✶ {xr =0} dM (x) r (ω) + t 0 ✶ {xr=0} dW ′ r (ω ′ ).
We state that W is a Brownian motion. Indeed, is a continuous process and its quadratic variation is

W t = t 0 1 2γ xr ✶ {xr =0} d M (x) r + t 0 ✶ {xr=0} d W ′ r + t 0 2 √ 2γ xr ✶ {xr =0} ✶ {xr=0} d M (x) , W ′ r = t 0 1 2γ xr ✶ {xr =0} 2γ xr dr + t 0 ✶ {xr=0} dr = t,
concluding thanks to Lévy's Theorem [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF]Theorem 3.3.16]. Now, notice that

E * M (x) t - t 0 2γ xr dW r 2 = E * M (x) t - t 0 ✶ {xr =0} dM (x) r + 2γ xr ✶ {xr=0} dW ′ r 2 = E * t 0 ✶ {xr=0} dM (x) r 2 = E * t 0 ✶ {xr=0} d M (x) r = E * t 0 ✶ {xr=0} 2γ xr dr = 0.
Then, M (x) has the brownian representation

M (x) t = t 0 2γ xr dW r , M (s) t = 0, M (v) t = 0.
Define γ = √ 2γ. We define the limit process of concentrations (X t ) t≥0 = (x t , s t , v t ) t≥0 , where x t := xt /v t , s t := ŝt /v t and v t := vt . Using Itô's formula [85, Theorem 1.5.5], we obtain that the dynamics of (X t ) t≥0 are characterized by the following SDE:

             dx t = µ(s t ) -β - u v t x t dt + γ x t v t dW t , ds t = - 1 Y µ(s t )x t + u v t (s in -s t ) dt, dv t = udt.
(5.11)

Remark 5.1 If we consider a null demographic stochasticity parameter γ = 0 and a null mortality rate β = 0, the model given in (5.11) coincides with the deterministic model (1.12). Then, (5.11) extends the deterministic model of SBR. We can also see that the diffusion coefficient in (5.11) shows that the amplitude of the variations of the biomass decrease with the increase of the volume. This is an interesting behavior that cannot be observed in the chemostat models, since in the chemostat the culture volume is maintained constant.

Existence of solutions of the controlled stochastic model

From now on, we consider u as a control parameter in the model instead of a constant parameter.

In this section, we assume that the death rate β is null and, without loss of generality, the yield coefficient will be assumed Y = 1. Let (Ω, F , P) be a complete probability space that supports a one-dimensional Brownian motion W = (W t ) t≥0 . The set where the control variable takes values is U = [0, u max ], and let u be an U -valued measurable proces defined on (Ω, F , P), which we call control process. We say that the control u = (u t ) t≥0 is admissible with respect to W = (W t ) t≥0 (or simply admissible) if u(•) is adapted with respect to the natural filtration of the Brownian W , that we denote by (F t ) t≥0 , and u(t, ω) ∈ U a.a (t, ω). The purpose of this section is to study the existence and uniqueness of solutions of the controlled stochastic differential equation

             dx t = µ(s t ) - u t v t x t dt + γ x t v t dW t , ds t = -µ(s t )x t + u t v t (s in -s t ) dt, dv t = u t dt.
(5.12)

for u admissible control. Define the coefficients of the stochastic differential equation (5.12):

b(x, s, v; u) =   µ(s)x -u v x -µ(s)x + u v (s in -s) u   and σ(x, s, v) =   γ √ x √ v 0 0   . (5.13)
These coefficients do not satisfy the Lipschitz condition (not even locally Lipschitz condition) and in consequence, existence of solutions of (5.12) is not guaranteed by the usual theorems (for instance, [START_REF] Yong | of Applications of Mathematics[END_REF]Theorem 1.6.16]). In order to study existence of solutions to (5.12) we consider a regularization and a truncation on the coefficients b(•) and σ(•). Let us define, for n ∈ N, the functions p n , q n by p

n (x) = -n✶ (-∞,n) (x) + x✶ [-n,n] (x) + n✶ (n,∞) (x), q n (x) = √ n|x|✶ [0,1/n) (|x|) + |x|✶ [1/n,n] (|x|) + √ n✶ (n,∞) (|x|),
and the following stochastic differential equations:

             dx n t = µ(s n t ) - u t v t p n (x n t )dt + γ q n (x n t ) √ v t dW t , ds n t = -µ(s n t )p n (x n t ) + u t v t (s in -s n t ) dt, dv t = u t dt, (5.14)
whose coefficients we denote by b n (x, s, v; u) and σ n (x, s, v), and solution X n . The parameter n controls the regularization of the square root near zero in the diffusion coefficients, and sets a bound on the variable x, turning the coefficients of (5.14) bounded and Lipschitz continuous and with sublinear growth with respect to (x, s) under some assumptions on the growth function µ(•). Under Assumption 5.1, the coefficients of the equation (5.14) satisfy the Lipschitz condition and the sublinear growth condition in (x, s), but not in (x, s, v). Nevertheless, we have the following lemma:

Lemma 5.3 For any admissible control u, for each n ∈ N there exists a unique solution of the equation (5.14).

Proof. Using Itô's formula for Z n , we get

Z n t = Z n 0 + t 0 [x n r -p n (x n r )] u r dr + t 0 γq n (x n r ) √ v r dW r , with Z n 0 = (x n 0 + s n 0 -s in )v 0 = (y + z -s in )w.
In the set {r ≤ τ n } we have p n (x n r ) = x n r and q n (x n r ) = √ x n r ; then, we can write

Z n t∧τ n = Z n 0 + γ t∧τ n 0 x n r v r dW r , (5.19)
from which the result follows.

Lemma 5.5 For all T ≥ 0, we have

sup n>0 E sup t≤T |X n t∧τ n | < ∞.
Proof. We analyze, as a first step, the behavior of s n t . The process s n t has paths of finite variation, so it can be studied in terms of its infinitesimal variations. It is non-negative, because if at some time instant t > 0 we have s n t = 0, the equation of s n reads

ds n t = u t v t s in dt ≥ 0,
so s n will be increasing. The process s n t is also bounded. Indeed, thanks to pathwise uniqueness, x n t is always nonnegative. If at some time instant t ≥ 0 we have that s n t ≥ s in , then

ds n t = -µ(s n t )p n (x n t ) + u t v t (s in -s n t ) dt ≤ 0,
so s n t will be decreasing. Now, consider the process Z n defined by (5.18). Thanks to (5.19), we can write

[x n t∧τ n + s n t∧τ n -s in ] v t∧τ n = [y + z -s in ] w + t∧τ n 0 γ x n r v r dW r . (5.20)
The process s n = (s n t ) t≥0 is always non-negative, and (v t ) t≥0 is non-decreasing and bounded from above by v max , so

x n t v 0 = x n t v t + s n t v t -s in v t + s in v t ≤ [x n t + s n t -s in ] v t + s in v max .
( Using the Burkholder-Davis-Gundy inequality [START_REF] Karatzas | Brownian motion and stochastic calculus[END_REF]Theorem 3.3.28], we obtain the existence of a constant C > 0 independent of n such that

E sup t≤T t∧τ n 0 γ x n r v r dW r ≤ CE T ∧τ n 0 γ2 x n r v r dr 1 2
, hence

E sup t≤T x n t∧τ n ≤ κ w + 1 w CE T ∧τ n 0 γ2 x n r v r dr 1 2 ≤ 1 w κ + 1 w C γE 1 + T ∧τ n 0 x n r v r dr ≤ 1 w [κ + C γ] + γC v max w E T 0 sup r ′ ≤r
x n r ′ ∧τ n dr .

Using Gronwall's inequality [48, Problem 5.2.7], and taking supremum over n > 0, we conclude that

sup n>0 E sup t≤T x n t∧τ n < ∞.
We define the process X = (x, s, v) in the following way: for n ≥ 1, n ∈ N,

X t = X n t , in the set {t ≤ τ n }.
It is easy to see that P-a.s., for all t ≥ 0,

X t∧τ n = ξ + t∧τ n 0 b(X r , u r )dr + t∧τ n 0 σ(X r )dW r , (5.22)
where b(•) and σ(•) are the coefficients defined in (5.13). The process X is defined up to the final time τ , and then, equation (5.12) in its integral form (5.22) has a solution defined in the time interval [0, τ ]. It remains to prove that τ is the first time that x t reaches 0 (possibly infinite).

Proposition 5.6 P-almost surely, τ = τ E := inf{t ≥ 0 : x t = 0}.

Proof. Suppose that P(τ = τ E ) > 0, and define

A = {ω ∈ Ω | τ (ω) = τ E (ω)}. Notice that if τ = ∞, then τ E = ∞, so necessarily A ⊆ {τ < ∞}.
Then, for almost all ω ∈ A there exists an increasing sequence

(n k (ω)) k∈N ր ∞ such that x τ n k (ω) (ω) = n k (ω) ր ∞.
Notice that there must exist a T > 0 such that P(A ∩ {τ ≤ T }) > 0. Indeed, if such T does not exist, for every N ∈ N we have P(A ∩ {τ ≤ N }) = 0. Notice that the sequence of sets A ∩ {τ ≤ N } is increasing with N and converges to A ∩ {τ < ∞}, concluding that P(A ∩ {τ < ∞}) = 0, which contradicts the initial assumption. Now, consider the T such that P(A ∩ {τ ≤ T }) > 0. For almost each ω ∈ A ∩ {τ ≤ T } we have that sup t∈[0,T ] x t∧τ n k (ω) (ω) = n k (ω), and then

sup n∈N E sup t∈[0,T ] x t∧τ n ≥ sup k∈N E sup t∈[0,T ] x t∧τ n k (ω) (ω)✶ A∩{τ ≤T } = ∞,
which contradicts Lemma 5.5.

The last proposition shows that the equation (5.12) has a solution up to τ , and that pathwise uniqueness holds. Now, we present a result that shows the behavior of the process with respect to extinction of the biomass. In the deterministic case the extinction time is infinite unless the initial condition is such that x(0) = 0, but in the stochastic framework the possibility of extinction must be considered (see [START_REF] Campillo | Stochastic modeling of the chemostat[END_REF][START_REF] Imhof | Exclusion and persistence in deterministic and stochastic chemostat models[END_REF]). This is a serious issue in terms of the depollution process, because extinction of the biomass before reaching the target stops the treatment, making it impossible to reach the target without considering other types of strategies such as adding biomass to the system when a critical level is reached. For equation (5.12), we have the following result: Proposition 5.7 Let u = (u t ) t≥0 be an admissible control with respect to W = (W t ) t≥0 and X u = (x u , s u , v u ) the solution of (5.12). The probability that x u = (x u t ) t≥0 hits 0 at some time instant is positive.

Proof. Consider the continuous local martingale

M u t = t∧τ 0 1 √ v u r dW r ,
whose quadratic variation is

M u t = t∧τ 0 1 v u r dr.
Since (v u t ) t≥0 es positive, non-decreasing and bounded up to τ , then S := lim t→∞ M u t = τ 0 1 v u r dr, P -a.s. Let us define, for r ∈ [0, ∞), the stopping time

T (r) = inf{t ≥ 0 | M u t > r}, r < S, ∞, r ≥ S.
According to [40, Theorem II.7.2'], defining B r = M u T (r) , G r = F T (r) , it is possible to construct an extension of (Ω, F , (F t ) t≥0 , P), namely ( Ω, F , ( Ft ), P) on which there is a Brownian motion B ′ independent of (N r , G r ) such that the process

Br = B r + B ′ r -B ′
r∧S is a standard unidimensional Brownian motion defined in ( Ω, F , ( Ft ), P). Then, X u (the solution of (5.12)) can be written in terms of the new time r and new Brownian motion ( Br ) r≥0 (in ( Ω, F , ( Ft ), P)) as follows: renaming xr = x T (r) , sr = s T (r) , ṽr = v T (r) and ũr = u T (r) , X ũ = (x ũ, sũ , ṽũ ) solves      dx r = [µ(s r )ṽ r -ũr ] xr dr + γ xr d Br , ds r = [-µ(s r )ṽ r xr + ũr (s in -sr )] dr, dṽ r = ũr ṽr dr,

with initial condition ξ = (y, z, w), for 0 ≤ r ≤ S (the process is stopped at r = S). We notice that

P ξ (∃t ∈ [0, ∞) : x u t = 0) = Pξ ∃r ∈ [0, ∞) : xũ r = 0 . Now,
r = (µ(s r )ṽ r -ũr )x r , β

r = µ max v max xr , b 1 (r, x) = b 2 (r, x) = µ max v max x), we get Pξ x ũ r ≤ xr , ∀r ∈ [0, ∞) = 1. (2)
Comparing the CIR model in [52, Proposition 6.2.4] and (5.24) (with a = 0, b = -µ max v max y σ = γ), the probability that x hits zero belongs to (0, 1). We conclude then

P ξ (∃t ∈ [0, ∞) : x u t = 0) = Pξ ∃r ∈ [0, ∞) : xũ r = 0 ≥ Py (∃r ∈ [0, ∞) : xr = 0) > 0.
Thus, the probability of extinction of biomass is positive.

The optimal reach-avoid problem

Consider the domain D

:= [0, ∞)×[0, s in ]×[v min , v max ], the target set C := (0, ∞)×[0, s out ]× {v max }, the full tank set V := (0, ∞)×[s out , s in ]×{v max }, the extinction set E := {0}×[0, s in ]× [v min , v max ],
and an initial condition ξ ∈ D. Consider a probability space (Ω, F , P) on which a one-dimensional Brownian motion W = (W t ) t≥0 is defined, and consider the Brownian natural filtration F := (F t ) t≥0 . We define the set of admissible controls U = {u : [0, ∞) × Ω → [0, u max ] | u is F -progresively measurable} , and denote X ξ,u the solution of the controlled stochastic differential equation (5.12) with control u ∈ U and initial condition ξ. This process takes values on the set D. For a Borel measurable set A ⊆ R 3 define the hitting times τ u A = inf {t ≥ 0 | X u t ∈ A}. Consider, for ξ ∈ D and u ∈ U , the cost functional given by the probability of reaching the target set before extinction of the biomass: J(ξ; u)

:= P ξ [τ u C ≤ τ u E ] , (5.25)
and the value function (P) V (ξ) := sup u∈U J(ξ; u).

(5.26)

Remark 5.2 Notice that the value function is bounded between 0 and 1 since it is defined as the supremum of probabilities of events, and it is discontinuous at the boundary of the domain. Indeed, for ξ ∈ C the value function takes the value V (ξ) = 1, for ξ ∈ E we have V (ξ) = 0, and cl(C) ∩ E = φ. The problem is that the reach set C and the avoid set E cannot be separated by nonintersecting open sets. The fact that the control on the set V is forced to take the value 0 poses another problem with respect to the standard formulation of the problem as a reach-avoid problem, since we have a state constraint. Nevertheless, for initial conditions ξ 0 = (x 0 , s 0 , v max ) in the set V it is possible to explicitly compute the probability of hitting the target before extinction.

Indeed, in the set V the constraints of the problem impose u ≡ 0, and then the process satisfies the SDE    dx t = µ(s t )x t dt + γ √ v max √ x t dW t , t ∈ [0, τ E ∧ τ C ],

ds t = -µ(s t )x t dt,

(5.27)

Notice that t → s t is decreasing up to τ E ∧ τ C P-a.s., thus we can define the time-change functions θ t := s 0 -s t = t 0 µ(s r )x r dr, T r := inf{t ≥ 0 : θ t > r},

(that depend of (x 0 , s 0 )) where the new time will be denoted by r. Denoting xr = x Tr , Wr = N

Tr , with N We make another time-change (that depends of s 0 ) in order to absorb the diffusion coefficient in (5.28). Define Proof. We study formula (5.30). First, the function s 0 → ϕ s 0 r is continuous: Then, s 0 → T s 0 r and s 0 → T s 0 = ϕ s 0 s 0 -sout are continuous, which proves the continuity of the first term in (5.30). For the continuity of the second term of the right-hand side, define (5.31)

|ϕ s 1 r -ϕ s 2 r | ≤ 1 µ(s 1 -r)µ(s 2 -r)
where the last inequality comes from Cauchy-Schwartz inequality and the fact that F (x, s) ≤ 1.

Notice that the exponent in the definition of G is a Gaussian variable with null mean and variance equal to From the Itô isometry, we have where we know that the terms inside the integrand converge to zero pointwise and are bounded; the term inside the expectation also converges to zero pointwise and is bounded. We conclude by the dominated convergence theorem that ∆ converges to zero in L 2 if s → s 0 .

E |∆| 2 ≤ 2E
not lead any initial condition with volume v < v(ξ 0 ) to the set C ∪ V because the volume will never reach its maximum value v max . The assumptions of works as [START_REF] Bouchard | Weak dynamic programming for generalized state constraints[END_REF][START_REF] Bouchard | Weak dynamic programming principle for viscosity solutions[END_REF][START_REF] Esfahani | The stochastic reach-avoid problem and set characterization for diffusions[END_REF] are not satisfied in this problem. However, we were able to prove the point 3. provided that the point 2. is fulfilled:

Lemma 5.10 Suppose that for every fixed control u ∈ U , the function (t, ξ) → X t,ξ,u τ t,ξ,u is continuous in probability. Then, the function (t, ξ) → J(t, ξ; u) is lower semicontinuous.

Proof. Consider (t 0 , ξ 0 ) ∈ [0, T ] × D and a sequence (t n , ξ n ) n∈N ⊆ [0, T ] × D converging to (t 0 , ξ 0 ) that realizes the lim inf, this is, l(t 0 , ξ 0 ) = lim n→∞ J(t n , ξ n , u) = lim inf (t,ξ)→(t 0 ,ξ 0 ) J(t, ξ, u). Denote τ n := τ tn,ξn,u and X τ n := X tn,ξn,u τ n . Since X τ n converges in probability to X τ 0 , then there exists a subsequence (t n k , ξ n k ) k∈N such that X τ n k ∈ Γ converges almost surely to X τ 0 ∈ Γ. Since g : Γ → [0, 1] is lower semicontinuous, and using Fatou's lemma,

J(t 0 , ξ 0 ) = E [g(X τ 0 )] = E lim inf k→∞ g (X τ n k ) ≤ lim k→∞ E tn k ,ξn k [g(X τ n k )]
= lim n→∞ E tn,ξn [g(X τ n )] = l(t 0 , ξ 0 ) = lim inf (t ′ ,ξ ′ )→(t 0 ,ξ 0 ) J(t ′ , ξ ′ , u).

Numerical simulations and conclusions

We consider the Monod (or Michaelis-Menten) growth function, which is quite popular in bioprocesses and which satisfies Assumption We compute numerically the function ṽ defined in Proposition 5.8, which is the probability of reaching the target before extinction occurs starting from an initial configuration in which the tank is full (v = v max ), when it is not possible to control the system. The results are shown in We see that the probability of hitting the target before extinction decreases with the parameter γ. For really γ, the probability of hitting the target before extinction is almost one, for initial conditions (x 0 , s 0 ) such that x 0 is far from 0. For large values of γ the effect of the demographic stochasticity is notorious, and the probability of extinction becomes large.

We also notice that ṽ(•) is discontinuous on the set {x = 0}, and it is non-smooth on the set {s = s out }. We see then the necessity to work with the lower and upper semicontinuous envelopes of the value function V (•).

In this work we propose a stochastic model of sequencing batch reactor. This model is obtained as the limit of an individual-based pure jump Markov processes, and it extends the deterministic model in the sense that when the parameter of demographic stochasticity is null we recover the deterministic dynamic. The obtained model is given by a degenerate stochastic differential equation whose diffusion coefficient only acts directly on the biomass component. We obtain a result that states that the probability of extinction of the biomass is positive, which gives an important message to practitioners that may work with this model. The extinction of biomass does not allow to define the problem of depollution of water in minimal time, and then the maximization of the probability of reaching the target before the extinction is a natural problem to study. This problem does not satisfy any of the usual assumptions, giving as a Chapter 6

Conclusions and perspectives

The purpose of this thesis is the mathematical study of optimal control problems of bioremediation of water resources. The main contributions of this thesis are:

• Chapter 2 -When there exists an active-dead zones scheme, the optimal control is a feedback law that depends only of the pollutant concentration of the active zone, that is, the zone that is treated directly by the reactor, and of the parameters of the growth function; it does not depend on the diffusion parameter or the volumes of the two zones. Nevertheless, the optimal treatment time depends of the diffusion parameter and the concentration of pollutant in the dead zone.

-The addition of a recirculation pump as a mean of enhancing the diffusion between the two zones is proved to be beneficial for the treatment.

• Chapter 3

-When we consider a scheme of two patches depolluted by the same bioreactor, the optimal control is again a feedback law that depends of the pollutant concentration of the most polluted zone. When pollution os homogeneous, the best is to maintain it homogeneous up to the final time, and the treatment time is independent of the diffusion.

-A high diffusion is favorable for the treatment when the initial concentrations are strongly different for the two zones. This is the case when one of the two concentrations is small compared to the threshold.

-For initial conditions in the two zones close to each other, a small diffusion parameter leads to faster treatment than a large diffusion.

-The possibility of treating two zones of the resource instead of one (as in the activedead model) increases the efficiency of the depolluting process, allowing to treat the most polluted zone as opposed to the active-dead zone that treates the less polluted one.

• Chapter 4

-The consideration of a gradient of concentrations by means of a configuration of two zones in series shows that the depollution process must be carried out in a way to over depollute the zone at the effluent of the bioreactor (which is the less polluted zone), in order to repollute it at some point of the process near the end of the process.

-A recirculation pump that helps to homogeneize the two zones is good. This recirculation pump is always optimally activated before the end of the process, and makes the repollution of the second zone to occur faster.

-A good estimation of the parameters of the system must be carried out in order to apply the optimal control, that in this case is not a feedback control, contrary to what happens in Chapters 2 and 3.

• Chapter 5

-We present a new stochastic model of sequencing batch reactor, which is a kind of model that has not been explored in the literature.

-We present a result of extinction of biomass regardless of the feeding strategy. This is an important message for practitioners to devise different alternatives to address the stop of the depollution process such as adding biomass to the tank when the biomass level is too low.

-We derive a dynamic programming principle for the problem of reaching the target concentration before that extinction of biomass occurs. This is an ongoing work.

The perspectives of future work of this thesis are the following.

• In Chapter 2 we proved that the optimal control law does not depend on the diffusion parameter D or the repartition of volumes v 1 , v 2 , for the case when these parameters are constant. We would like to study the effect of considering a diffusion parameter and a repartition of volumes that depends of the time. Our guess is that the main result remains true.

• The study of the effect of splitting the resource into several zones connected in different configurations. In Chapter 3 we proved that when considering two patches connected to the bioreactor, it is optimal to treat the most polluted zone at the optimal flow rate up to the moment in which the two concentrations are even (if this happens before reaching the desired pollutant concentration), and from that moment it is optimal to keep them equal up to the time in which the pollutant concentration reaches the threshold. If we consider an extension of several zones interconnected to each other by diffusion and some (or all) of them connected to the bioreactor, a natural candidate to the optimal control is to consider that the most polluted zone will be treated until its concentration reaches the same level of the second most polluted zone, and then treat both of them keeping the concentrations even until the level of pollution of the third most polluted zone is reached, and so on.

Consider a simple extension that consists of adding a third zone in the two patches model studied in Chapter 3. This third zone is adjacent to the other two zones that have access to the bioreactor and is connected by diffusion with them, as Figure 6.1 shows. In this problem, we see that there is no direct way to depollute the intermediate zone, it must be depolluted in an indirect way through the control of the pollution in zones 1 and 3. This is a problem whose velocity set is not convex, bringing the same issues that the problem studied in Chapter 3; it also has controllability issues, in the sense that depending on the parameters of the system, the configuration s 1 = s 2 , s 1 = s 3 , s 2 = s 3 cannot be maintained. It is also not clear what to do when the intermediate zone is the most polluted one. We would like to extend the results of Chapters 2-4 to several zones connected in different configurations, and to explore the effects of combining the parallel and series configurations to generate a more accurate model of inhomogeneity in the resource.

• In Chapter 3 we proved that despite of the velocity set being not convex, the solution of the relaxed problem is a non-convexified control. This is due to the particular structure of the dynamic that depends of two control parameters α and s ⋆ r , where α acts linearly and s ⋆ r not.

whenever s 1 = s 2 , the velocity set is not convex. Nevertheless, its extreme points are given by points parametrized by values α = 0 and s ⋆ r = ŝ⋆ r (s 2 ), or α = 1 and s ⋆ r = ŝ⋆ r (s 1 ), or s ⋆ r = 0 and α ∈ [0, 1]. This last option is not admissible. -when s 1 = s 2 , the velocity set becomes convex, so in this case there is no need of convexification.

An interesting possible future work can be to obtain results of Chapter 3 for systems whose dynamics are governed by equations such that the velocity set is non-convex but satisfy the previous property.

• Another natural extension of the work presented in this thesis concerns the study of the depollution of a water resource polluted by two substrates. The water is treated by using a continuous stirred bioreactor consisting of two zones. In the first zone the substrate of In this problem the velocity set is non-convex, and the Hamiltonian is non-smooth. Then, a relaxation of the problem must be performed, and tools such as the extended maximum principle of Pontryagin with subdifferentials must be used. We would like to study in a first step the minimal time optimal control problem with constant volumes, and then, as a second step, to extend the result to obtain the best initial repartition of volumes depending of the initial state to perform the treatment in minimal time.

• In Chapter 5 we study and derive a dynamic programming principle for the problem of the maximization of the probability of attaining a target prior to extinction of biomass.

2 . 1

 21 vii viii L'apprentissage ne se termine jamais. ix x tant 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Definitions and preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Optimal control problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .5 Stochastic modelling of sequencing batch reactors for wastewater treatment 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Stochastic SBR model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Existence of solutions of the controlled stochastic model . . . . . . . . . . . . 5.4 The optimal reach-avoid problem . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Numerical simulations and conclusions . . . . . . . . . . . . . . . . . . . . . Comparison between optimal treatment time and the time given by the best constant control, from initial condition z = (5, 8) . . . . . . . . . . . . . . . . 2.2 Comparison between optimal treatment time and the time given by the best constant control, from initial condition z = (6, 7) . . . . . . . . . . . . . . . . 2.3 Comparison between optimal treatment time and the time given by the best constant control, from initial condition z = (2, 10) . . . . . . . . . . . . . . . 2.4 Comparison of optimal times with respect to diffusion parameter D from initial condition z = (5, 8) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Comparison of optimal times with respect to diffusion parameter D from initial condition z = (6, 7) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Comparison of optimal times with respect to diffusion parameter D from initial condition z = (2, 10) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Time comparisons (in hours) for r = 0.3 and target value s = 1 (top), 0.1 (bottom) [g/l] (initial condition s(0) and diffusion parameter d are given in [g/l] and [1/h], respectively) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Comparison of treatment times for the optimal strategies . . . . . . . . . . . . 4.2 Comparison of optimal times with respect to maximum recirculation parameter QR for different initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . xv List of Figures 1.1 Scheme of water treatment . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Industrial bioreactor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Scheme bioreactor-settler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Typical growth functions. On the left, the Monod uptake function; on the right, the Haldane uptake function. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Trajectories on the phase portrait, depending on the diffusion parameter D.

4 ] 1 .

 41 10 First model of inhomogeneity: the active-dead zones configuration . . . . . . . 1.11 Second model of inhomogeneity: the model with two patches . . . . . . . . . . 1.12 Third model of inhomogeneity: the configuration in series with recirculation . . 1.13 Total pollutant concentration in the resource of the full dynamics (3.6) with the strategy (3.25), for different values of ǫ. . . . . . . . . . . . . . . . . . . . . . 1.14 Velocity set of the dynamic (1.30) when s 1 = s 2 . . . . . . . . . . . . . . . . . 2.1 Connections between active zone, dead zone and bioreactor. . . . . . . . . . . 2.2 Connections between active zone, dead zone and bioreactor, with the recirculation pumps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Modeling scheme of the treatment of two interconnected patches (definitions of the variables and parameters are given in Section 3.2) . . . . . . . . . . . . 3.2 Graphs of µ(•) and corresponding γ(•). . . . . . . . . . . . . . . . . . . . . . . 3.3 Optimal paths for d = 0.1[h -1 ] (left) and d = 10[h -1 ] (right) with r = 0.3. . . . 3.4 Level sets (in hours) of V 0 (left) and V ∞ (right) for r = 0.3. . . . . . . . . . . .

Figure 1 . 1 .Figure 1 . 2 .

 1112 Figure 1.1. Scheme of water treatment

Figure 1 . 4 .

 14 Figure 1.4. Typical growth functions. On the left, the Monod uptake function; on the right, the Haldane uptake function.

Figure 1 . 5 .

 15 Figure 1.5. Trajectories on the phase portrait, depending on the diffusion parameter D. In green, the washout equilibrium; in red, the nontrivial equilibrium.

2 .

 2 If s in > ν j , then E 0 and E 1 are local attractors, their basins of attractions are nonempty open sets, and the complement of the union of the two basins of attraction has zero Lebesgue measure.

Figure 1 . 8 .

 18 Figure 1.8. Relation between fictitious time τ and real time t

m 1

 1

Figure 1 . 10 .Figure 1 . 11 .

 110111 Figure 1.10. First model of inhomogeneity: the active-dead zones configuration The second model consists of two clearly differentiated connected patches, similar to the active-dead zone. The difference is that the resource is treated by means of a continuously stirred bioreactor connected to each of the patches by independently operated pumps (see Figure 1.11)

Figure 1 . 12 .

 112 Figure 1.12. Third model of inhomogeneity: the configuration in series with recirculation

Figure 1 . 13 .

 113 Figure 1.13. Total pollutant concentration in the resource of the full dynamics (3.6) with the strategy (3.25), for different values of ǫ.

Figure 1 . 14 . 2 Proposition 1 . 5 (

 114215 Figure 1.14. Velocity set of the dynamic (1.30) when s 1 = s 2

Figure 2 . 1 .

 21 Figure 2.1. Connections between active zone, dead zone and bioreactor.

. 5 )

 5 Note that we have not yet imposed the condition s⋆ r ∈ [0, s 1 ].

3 Table 2 . 4 . 3 Table 2 . 5 .

 324325 3204 × 10 3 5.6069 × 10 2 4.7702 × 10 2 3 3.1402 × 10 3 1.8718 × 10 3 1.7450 × 10 3 1 1.0064 × 10 4 7.7433 × 10 3 3.6877 × 10 Comparison of optimal times with respect to diffusion parameter D from initial condition z = (51680 × 10 3 5.4748 × 10 2 4.7585 × 10 2 3 2.9992 × 10 3 1.8585 × 10 3 1.7452 × 10 3 1 9.9235 × 10 3 7.7265 × 10 3 3.5543 × 10 Comparison of optimal times with respect to diffusion parameter D from initial condition z = (6370 × 10 2 4.5386 × 10 2 3.4634 × 10 2 3 2.6976 × 10 3 1.7657 × 10 3 1.6173 × 10 3 1 9.6207 × 10 3 7.1613 × 10 3 3.5902 × 10 3

Chapter 3 Optimal

 3 feedback synthesis and minimal time function for the bioremediation of water resources with two patches This chapter corresponds to the published article Ramírez C., H.; Rapaport, A.; Riquelme, V. Optimal Feedback Synthesis and Minimal Time Function for the Bioremediation of Water Resources with Two Patches. SIAM J. Control Optim. 54 (2016), no. 3, 1697-1718.

Figure 3 . 1 .

 31 Figure 3.1. Modeling scheme of the treatment of two interconnected patches (definitions of the variables and parameters are given in Section 3.2)

γ 1 Figure 3 . 2 .

 132 Figure 3.2. Graphs of µ(•) and corresponding γ(•).

Figure 3 . 3 .

 33 Figure 3.3. Optimal paths for d = 0.1[h -1 ] (left) and d = 10[h -1 ] (right) with r = 0.3.

Figure 3 . 4 .

 34 Figure 3.4. Level sets (in hours) of V 0 (left) and V∞ (right) for r = 0.3.

1 ]Figure 3 . 5 .Figure 3 . 6 .

 13536 Figure 3.5. Trajectories and controls generated by the two-pumps and one-pump optimal feedback, for r=0.3, s=1 [g/l] and s(0)=(3,10) [g/l]. On the left d=0.1[1/h], and on the right d=10[1/h].

  -

Figure 4 . 2 .

 42 Figure 4.2. Backward integration of the extremals, with QR = 0.

Figure 4 . 3 .

 43 Figure 4.3. Backward integration of the extremals, with QR = 1.

Figure 4 . 4 .

 44 Figure 4.4. Backward integration of the extremals, with QR = 10.

. 21 )

 21 Replacing(5.20) in(5.21), with v 0 = w,x n t∧τ n w ≤ [y + z -s in ] w + s in v max + t∧τ n 0 γ x n r v r dW r .Define κ := [y + z -s in ] w + s in v max . Taking supremum on [0, T ] and expectation, E w sup t≤T x n t∧τ n ≤ κ + E sup t≤T t∧τ n 0 γ x n r v r dW r .

1 µ

 1 r )x r dW r , the law of this process satisfies the following SDE independent of s t :dx r = dr + γ √ v max 1 µ(s 0 -r) d Wr , x0 = x 0 . (5.28) for r ∈ [0, θ τ ], with θ τ = θ τ C ∧ θ τ E .In this new time scale, in the event {τ C < τ E } the time to the target is θ τ C = s 0 -s out (constant, depending only of s 0 ), and the time to extinction is θ τ E . The coefficient that accompanies the brownian part in equation (5.28) does not explode, and this equation provides an explicit expression for xr :xr = x 0 + r ∧ θ τ E ∧ (s 0 -s out ) + γ √ v max r∧θτ E ∧(s 0 -sout) 0 (s 0 -h) d Wh .

d

  := inf{r ≥ 0 : ϕ r > h}, In the time h, defining xh = x Th and Wh = N Wh , xh solves dx h = v max γ2 µ(s 0 -Th )dh + d Wh , x0 = x 0 . (5.29) for h ∈ [0, ϕ θτ ]. Define for simplicity, T s 0 = T := ϕ s 0 θ C , TE = ϕ s 0 θτ E. Notice that in the set{τ C < τ E } we have T s 0 = ϕ s 0 -sout . Define also A s 0 h = A h := vmax γ2 h 0 µ(s 0 -Tα )dα. For a process M denote M ⋆ t = max s∈[0,t] M s .In these terms, hitting the target before extinction corresponds to the event { TE > T, x T > 0}.P ξ 0 [τ C < τ E ] = P TE > T, x T > 0 = P min h∈[0, T ] x 0 + A h + Wh > 0, TE > T = P min h∈[0, T ] A h + Wh > -x 0 .Consider the probability law Q under which the process Wh -A h is an Fh -Brownian motion (where Fh = F TT h), which has density with respect to P given by M 0 -Tα ) 2 dα .Then,P ξ 0 [τ C < τ E ] = P ∀h ∈ [0, T ], Wh > -x 0 -A h = Q ∀h ∈ [0, T ], Wh -A h > -x 0 -A h = E Q ✶ {min h∈[0, T ] Wh >-x 0 } = E ✶ {min h∈[0, T ] Wh >-x 0 } M T = E ✶ {min h∈[0, T ] Wh >-x 0 } exp 0 -Tα ) 2 dαFinally, we get the expression for the cost on V:P ξ 0 [τ C < τ E ] = exp -0 -Tα ) 2 dα × E ✶ { W ⋆ T <x 0 } exp -v max γ2T 0 µ(s 0 -Tα )d Wα .

(5. 30 ) 5 . 8

 3058 Proposition The function V ∋ ξ → ṽ(ξ) := P ξ [τ u C ≤ τ u E ] is continuous.

  1 -h) -µ(s 2 -h)|dh ≤ r µ(s 1 -r)µ(s 2 -r) max h∈[0,r] |µ(s 1 -h) -µ(s 2 -h)|

F 1 2+

 1 (x 0 , s 0 ) = ✶ W ⋆ T <x 0 , G(s 0 ) = exp -v max γ2 T s 0 0 µ(s 0 -T s 0 α )d Wα .Then, for ξ = (x, s, v max ) ∈ V,E|F (x 0 , s 0 )G(s 0 ) -F (x, s)G(s)| ≤ E [|F (x 0 , s 0 ) -F (x, s)|G(s 0 )] + E [F (x, s)|G(s 0 ) -G(s)|] ≤ E |F (x 0 , s 0 ) -F (x, s)| 2 1 2 E G(s 0 ) 2 E [|G(s 0 ) -G(s)|] .

≤ exp p 2 2 v

 2 0 -T s 0 α ) 2 dα, so G(s 0 ) is finite (depending on s 0 ). Indeed, with an appropriate change of variables, we can write for any p ≥ 1E [G(s 0 ) p ] = exp p max γ2 (s 0 -s out )µ(s 0 ) < ∞. (5.32)From the previous expression, we can see that for any s ≤ s in it holdsE [G(s) p ] ≤ E [G(s in ) p ] <∞, and we conclude that {G(s)} s∈[sout,s in ] is uniformly integrable. To prove the convergence of G(s) to G(s 0 ) in L 1 when s → s 0 it is enough to prove that 0 -T s 0 α )d Wα → 0 in L 2 if s → s 0 .

  T s α ) -µ(s 0 -T s 0 α )| 2 dα + 2E T s T s 0 |µ(s 0 -T s 0 α )| 2 dα ≤ 4L 2 µ E T s 0 (|s -s 0 | 2 +| T s α -T s 0 α | 2 )dα + 2µ 2 max E| T s 0 -T s |,(5.33)

3 . 1 :

 31 µ(s) = µ max s K s + s , with the parameters µ max = 2[h -1 ] and K s = 1[gl -1 ]. The threshold that defines the target has been chosen as s out = 1[gl -1 ]. The input concentration of substrate is s in = 5[gl -1 ].

Figures

Figure 5 . 1 . 1 0 5 . 2 . 5 0 5 . 3 .

 51152553 Figure 5.1. Plot of the function ṽ for γ = 1

1 q = q 1 +q 2 3 v 1 vFigure 6 . 1 .

 13161 Figure 6.1. First extension: model with three patches and two pumps

Figure 6 . 2 . 1 Q v 1 , 2 Q v 2

 621122 Figure 6.2. Second extension: problem with two pollutantsIn this problem it is also possible to suppose that the volume of the resource is large compared to the volume of the reactor, which leads to a quasi-steady-state approximation of the dynamics. Here it is convenient to keep Q as a control variable instead of replacing it by the concentration of equilibrium in each zone of the bioreactor (as it was done in Chapters 2-4). The equations that model the time evolution of this problem are given byṡa = Q (s ⋆ a1 (Q, s a ) -s a ) , ṡb = Q (s ⋆ b2 (Q, s a , s b ) -s b ) ,where s ⋆ a1 (Q, s a ) and s ⋆ b2 (Q, s a , s b ) are the corresponding concentrations of equilibrium inside each zone of the bioreactor, given by the formulas

  [t 1 , t 2 ], where t 2 is the first time in which ψ(t) ≥ 0;

.35) a.e. t ∈ [0, t 1 ], where t 1 is the first time that γ(t) ≤ 1 (it may occur that t 1 = 0), 2. Q opt R (t) = QR and s opt r (t) = ŝ⋆ r (s opt (t), γ(t)), where ŝ⋆ r (s, γ) satisfies (4.15), a.e. t ∈

Table 2 . 1 .

 21 .1-2.3. Comparison between optimal treatment time and the time given by the best constant control, from initial condition z =(5,[START_REF] Boopathy | Factors limiting bioremediation technologies[END_REF]

	s	T feedback	T constant	Gain
	5 1.3204 × 10 3 3 3.1402 × 10 3 1 1.00644 × 10 4 1.13642 × 10 4 11.44% 0.31% 1.3245 × 10 3 1.77% 3.1970 × 10 3
	s	T feedback	T constant	Gain
	0.72% 2.43% 1 9.9235 × 10 3 1.12667 × 10 4 11.92% 5 1.1680 × 10 3 1.1765 × 10 3 3 2.9992 × 10 3 3.0739 × 10 3

Table 2 . 2 .

 22 Comparison between optimal treatment time and the time given by the best constant control, from initial condition z = (6, 7) × 10 3 2.7086 × 10 3 0.41% 1 9.6207 × 10 3 1.05199 × 10 4 8.54%

	s	T feedback	T constant	Gain
	5 9.370 × 10 2 3 2.6976	9.3280 × 10 2 0.11%

Table 2 . 3 .

 23 Comparison between optimal treatment time and the time given by the best constant control, from initial condition z =[START_REF] Alford | Bioprocess control: Advances and challenges[END_REF][START_REF] Bouchard | Weak dynamic programming for generalized state constraints[END_REF]

Table 2 . 6 .

 26 Comparison of optimal times with respect to diffusion parameter D from initial condition z =[START_REF] Alford | Bioprocess control: Advances and challenges[END_REF][START_REF] Bouchard | Weak dynamic programming for generalized state constraints[END_REF]

Table 3 .

 3

	5, 0)	0.42	0.01	0.42	0.01	0.42	0.01
	Increase:			(+ 1.45 %) (+ 0.00 %)	(+ 0.00 %)	(+ 0.00 %)
	s(0) = (3, 0)	1.01	0.06	1.05	0.06	1.01	0.06
	Increase:			(+ 3.90 %) (+ 0.85 %)	(+ 0.00 %)	(+ 0.00 %)
	s(0) = (4, 0.5)	1.33	2.17	1.39	2.23	1.37	2.21
	Increase:			(+ 4.68 %) (+ 2.62 %)	(+ 2.73 %)	(+ 1.55 %)
	s(0) = (4, 1.5)	3.20	3.65	3.67	3.75	8.27	3.72
	Increase:			(+ 14.76 %) (+ 2.58 %) (+ 158.27 %) (+ 1.91 %)
	s(0) = (4, 4)	5.45	5.45	5.74	5.71	18.25	5.53
	Increase:			(+ 5.43 %) (+ 4.90 %) (+ 235.01 %) (+ 1.59 %)
		V d		T ⋆ cst		T ⋆ one
		d = 0.1 d = 10	d = 0.1	d = 10	d = 0.1	d = 10
	s(0) = (1.5, 0)	25.95	34.12	38.65	38.81	34.03	34.14
	Increase:			(+ 48.93 %) (+ 13.74 %) (+ 31.14 %) (+ 0.05 %)
	s(0) = (3, 0)	32.91	39.91	50.08	50.12	45.89	40.15
	Increase:			(+ 52.18 %) (+ 25.58 %) (+ 39.45 %) (+ 0.60 %)
	s(0) = (4, 0.5)	41.08	42.86	58.65	58.02	61.51	42.94
	Increase:			(+ 42.77 %) (+ 35.37 %) (+ 49.74 %) (+ 0.1 %)
	s(0) = (4, 1.5)	43.69	44.37	63.59	63.28	70.81	44.49
	Increase:			(+ 45.57 %) (+ 42.61 %) (+ 62.08 %) (+ 0.27 %)
	s(0) = (4, 4)	45.94	45.94	71.67	71.04	81.58	46.17
	Increase:			(+ 56.02 %) (+ 54.64 %) (+ 77.60 %) (+ 0.51 %)

1. Time comparisons (in hours) for r = 0.3 and target value s = 1 (top), 0.1 (bottom) [g/l] (initial condition s(0) and diffusion parameter d are given in [g/l] and [1/h], respectively)

  2 < T opt such that 1. Q opt R (t) = 0 and w opt (t) = ŵ(s opt (t), γ(t)) a.e. t ∈ [0, t 1 ], where ŵ(s, γ) satisfies (4.15), where t 1 is the first time that γ(t) ≤ 1 (it may occur that t 1 = 0), 2. Q opt R (t) = QR and w opt (t) = ŵ(s opt (t), γ(t)) a.e. t ∈ [t 1 , t 2 ], where ŵ(s, γ) satisfies (4.15), where t 2 is the first time in which ψ(t) ≥ 0; 3. Q opt R (t) = QR and w opt (t) = 1 a.e. t ∈ [t 2 , T opt ]. Proof. According to Pontryagin's Maximum Principle, an optimal control u opt has to maximize for a.e. t ∈ [0, T opt ] time the Hamiltonian H given in (4.10) or, equivalently, since λ(t) ≤ 0 and γ(t) ≥ 0 (Lemma 4.3), maximize the functions w → φ(s opt (t), γ(t), w) with w ∈ [0, 1]

  1 ) = 1. For a.e. t ∈ [0, t 1 ], point 1. of Lemma 4.5 is satisfied. For t > t 1 Lemma 4.6 states that γ(t) < 1 and then Lemma 4.5 states that for t ≥ t 1 , Q opt R . Now, Remark 4.2 states that ψ(t 1 ) < 0 and we know that ψ(T opt ) > 0, so there exists a time instant t 2 < T opt such that ψ(t 2 ) > 0 and t 1 < t 2 . Lemma 4.7 says that a.e. t ∈ [t 1 , t 2 ] point 2. of Lemma 4.5 holds, and a.e. t ∈ [t 2 , T opt ] the point 3. of Lemma 4.5 holds.

Table 4 . 1 .

 41 Comparison of treatment times for the optimal strategies

Table 4 .

 4 2. Comparison of optimal times with respect to maximum recirculation parameter QR for different initial conditions Several optimal trajectories in the phase portrait are shown in Figures 4.2

	5, 0.1)	739.25	545.00 (-26.0%)	336.47 (-54.5%)
	(2.5, 0.5)	951.33	761.63 (-19.9%)	570.49 (-40.0%)
	(2.5, 1.5)	1504.33 1327.03 (-11.8%) 1160.94 (-22.8%)

  consider the process x = (x r ) r≥0 defined as the solution of dx r = µ max v max xr dr + γ√ xr d Br ,

	x0 = y.	(5.24)
	Strong existence and pathwise uniqueness hold for the equation (5.24) (from [52, Theorem 6.2.3]),

and it is defined for r ∈ [0, ∞). Since µ(•) is bounded, then for any r ≥ 0, µ(s r ) ≤ µ max and ṽr ≤ v max . Then, thanks to comparison principles (for instance, see

[START_REF] Ikeda | Stochastic differential equations and diffusion processes[END_REF] Theorem VI.1.1]

, with β

Remerciements

Nacional Convocatoria 2013 folio 21130840. Cette thèse a été développée le projet BION-ATURE de CIRIC INRIA CHILE, et elle a été partiellement financé par les projets CONICYT REDES 130067 et MOMARE SticAmsud. Aussi, je remercie le soutien financier de CON-ICYT sous le projet ACT 10336, les projets FONDECYT 1080173, 1110888, 1160204, le projet BASAL du Centro de Modelamiento Matemático de l'Université du Chili, et le projet MathAmsud N • 15MATH-02. Je remercie le soutien financier du Départemento de Postgrado y Postítulo de la Vicerrectoría de Asuntos Académicos de l'Université du Chili et l'Institut Français (l'Ambassade de la France au Chili). Je remercie ma famille par son appui le long de ma vie et de mes études, et ma fiancée Pascale, qui m'a appuyée tout ce temps, même dans mes longs séjours à l'étranger sans sa compagnie. xi xii

interesting features for the practitioners and controllers:

1. it does not require knowledge of the diffusion parameter D to be implemented, and 2. if the ratio r of the volumes of the two patches is not known, the optimal trajectory can be approximated by a regularization of the bang-bang control about the neighborhood of the set s 1 = s 2 that keeps the trajectory in this neighborhood.

Furthermore, is has been shown in simulations that the benefit of using two pumps instead of one can be significant when the diffusion is low. We have also proposed explicit bounds on the minimal-time function, characterizing the extreme cases d = 0 and d = +∞. We have shown that a large diffusion rate increases the treatment time when the pollution concentration is above the desired threshold in both zones, while in contrast, it can be beneficial when the concentration in one of the two zones is below the desired threshold. This remarkable feature could serve practitioners in the choice of pump positioning in an originally clean water resource that is suddenly affected by a local pollution. Such an investigation could be the matter of future work.

Appendix

Proposition 3.14 For any ǫ > 0, the feedback strategy (3.25) applied to the full dynamics (3.6) with x r (0) > 0 drives the state to the target in finite time.

Proof. Without any loss of generality, we assume that s 1 (0) ≥ s 2 (0) (the proof is similar when s 1 (0) ≤ s 2 (0)).

If s 1 (0) > s 2 (0), we prove that s 1 = s 2 is reached in finite time. If not one, one should have s 1 (t) > s 2 (t) with s 1 (t) ≥ s for any t > 0. This implies to have α ⋆ (t) = 1 and s ⋆ r (t) = ŝ⋆ r (s 1 (t)) at any time t > 0 and one has from equations (3.6):

which implies that the trajectories are bounded. For any σ ≥ s, ŝ⋆ r (σ) being the unique maximizer of the function β(σ, •), one has

.

The function µ(•) being increasing and concave, one obtains the inequality

Furthermore, one can write

satisfying the following relations:

Under Assumption 4.1, system (4.3) has unique solution for s ⋆ r ∈ [0, s 1 ]. Notice that this imposes a constraint on s ⋆ r for the positivity of the biomass concentration in the bioreactor. For simplicity, we keep the notation s r instead of s ⋆ r for the concentration of equilibrium in the bioreactor. Denoting α i = v r /v i (i = 1, 2) results in the following dynamics:

From now on, we work with the dynamics (4.4). Notice that the problem can be reformulated with the controls (s r , Q R ) instead of (Q, Q R ), but this formulation has the inconvenient that the set where the control takes values depends of the state variable. Nevertheless, we can make a change of variable of the form s r = ws 1 , with w ∈ [0, 1] and take (w, Q R ) as the control variable instead of (Q, Q R ). We keep s r as variable when there is no confusion, always keeping in mind that s r = ws 1 . Define the set of admissible controls

Proof. Suppose s = (s 1 , s 2 ) ∈ D and a control u ∈ U. By the continuity of the solutions of (4.4), the only way to leave the set D is through s 1 = s 2 . The equation that s 1 -s 2 satisfy is

so the set D cannot be abandoned.

Optimal control problem

Due to Lemma 4.2 we may suppose that the initial condition for the process belongs to the set D. Consequently, the target set can be considered as

Chapter 5

Stochastic modelling of sequencing batch reactors for wastewater treatment

This chapter corresponds to an ongoing work.

Introduction

Sequencing batch reactors (SBR) are a biological water treatment system. They are widely used for industrial and municipal wastewater. They are a time-oriented system with flow, energy input, and tank volume varying according to some prescribed strategy. In its most basic form, the SBR system consists of one or many tanks that operate in a fill-and-draw basis: the tank is filled during a discrete period of time and then operated as a batch reactor. The essential difference between the SBR and a conventional continuous bioreactor is that inside each tank functions as equalization, aeration, and sedimentation can be performed in a time sequence rather than a spatial sequence. This kind of waste-water treatment system is well described in [START_REF] Irvine | Sequencing batch biological reactors: An overview[END_REF][START_REF] Irvine | Sequencing batch reactors for biological waste-water treatment[END_REF].

The cycle of operation for each tank in an SBR system consists of five basic operating modes performed in a time sequence: fill, react, settle, draw, and idle modes. During fill mode the wastewater is received in the tank; during react mode the desired reactions take place; in settle mode the microorganisms are separated from the treated water; in draw mode the treated effluent is discharged, and idle mode corresponds to the period comprised between the discharge of the tank and its refilling. The time of a cycle is the total time between the beginning of fill up to the end of idle. The necessary time to achieve the total cycle may be long, and this may be economically expensive. On the other hand, one should expect that a longer treatment time leads to a better depollution level.

The problem consists in the minimization of the time of the treatment in each cycle of operation, in which the time associated to draw and idle modes is fixed. Then, by changing Proof. Given u = (u t ) t≥0 admissible control, the volume process v = (v t ) t≥0 has the explicit representation v t (ω) = w + t 0 u r (ω)dr.

(5.15)

Define the process Xn = (x n t , s n t ) t≥0 as the solution of the following stochastic differential equation with random coefficients bn (t, x, s; ω) and σn (t, x, s; ω) given by

(5.16)

We notice that given the control u(•), for each (x, s) ∈ R 2 , t ≥ 0, and ω ∈ Ω,

hence, on (Ω, F , (F t ) t≥0 , P), with control u, Brownian motion W , and initial condition ξ 0 , the solutions of (5.14) and (5.16) are the same. The coefficients bn (•) and σn (•) are Lipschitz continuous in (x, s) and have sublinear growth in (x, s), uniformly in t:

where L µ is the Lipschitz constant of µ, and

Consequently, according to [85, Theorem 1.6.16], for all n ∈ N the equation (5.16) has a unique solution; then, for all n ∈ N, equation (5.14) has a controlled solution X n = (X n t ) t≥0 .

Let us define for each n ∈ N the stopping time

n}, which is the first time when the dynamics (5.14) are truncated, and τ := sup n∈N τ n . The stopping time τ is well defined, because it is the limit of a sequence of increasing stopping times, and thanks to the pathwise uniqueness, for m > n it is not difficult to see that X m t = X n t in the set {t ≤ τ n }.

Lemma 5.4 The process Z n

•∧τ n = (Z n t∧τ n ) t≥0 , where Z n t is defined by

is a continuous local martingale with respect to (F t ) t≥0 .

Then,

Wh < x 0 ∨ x .

(5.34)

When (x, s) → (x 0 , s 0 ), the first term of (5.34) converges to P max h∈[0, T s 0 ] Wh = x 0 , which is null because the running maximum of a Brownian motion has density. The second term of (5.34) has an explicit formula:

that converges to 0 as (x, s) → (x 0 , s 0 ), concluding the proof. Now, we restate the problem (P) as a problem without state constraints. For this purpose, we consider the cost function ṽ(•) on V introduced in Proposition 5.8, and we define the set Γ := C ∪ V ∪ E. For every t ≥ 0, define F t = (F t r ) r≥0 where F t r is the σ-field σ(W θ -W t : t ≤ θ ≤ r ∨ t). We define the set of admissible controls

and denote X t,ξ,u the solution of the controlled stochastic differential equation (5.12) with control u ∈ U t and initial condition ξ at initial time t. Define

and the new value function

Notice that the cost function can be written as

where g : Γ → R defined by

is a lower semicontinuous function in its domain. We introduce, as in [80, Section 3.2], the lower and upper semicontinuous envelopes of V ,

We are interested in developing a dynamic programming principle for the value function. Our conjecture is the following: Conjecture 5.9 (Dynamic Programming Principle) Let (t, ξ) ∈ D be fixed. For every F tstopping time h, the value function V (•) satisfies the dynamic programming inequalities:

We expect to prove Conjecture 5.9 by following the approach given in [START_REF] Bouchard | Weak dynamic programming principle for viscosity solutions[END_REF][START_REF] Esfahani | The stochastic reach-avoid problem and set characterization for diffusions[END_REF]. For this, we need to establish the following steps: 1. To prove that for every control u ∈ U , the function (t, ξ) → τ t,ξ;u is continuous in probability.

2. To prove that for every control u ∈ U , the function (t, ξ) → X t,ξ;u τ t,ξ,u is continuous in probability.

3. To prove that the function (t, ξ) → J(t, ξ; u) = E t,ξ [g(X u τ u )] is lower semicontinuous. 4. To prove a Markov property of the type: for every control u ∈ U and F t -stopping time h, E t,ξ [g(X u τ u )|F h∧τ u ] (ω) = E h∧τ u (ω),X t,ξ;u h∧τ u (ω) g(X u ω,h∧τ u τ u ω,h∧τ u ) , P -a.s., where a.a (ω, t), u ω,t (ω ′ , •) are the controls that coincide with u(ω, •) in the interval [0, t] (as in [START_REF] Claisse | A Pseudo-Markov Property for Controlled Diffusion Processes[END_REF]).

The proof of the previous steps is non-standard, because of the particularities of the problem. Indeed, the first difficulty is the treatment of the extinction, since once extinction occurs, the system is frozen in the state in which extinction occurred. We need to prove for a fixed control u ∈ U the continuity in probability of the hitting times τ t,ξ;u E and of the hitting times τ t,ξ;u C∪V with respect to the initial conditions (t, ξ). But since the set E is absorbent and its closure intersects the closure of the set C ∪ V, the analysis of the continuity of the hitting times is more complicated than the usual case where the diffusion parameter is nondegenerate. Also, the constraints of the problem as well as measurability issues impose that the convergence of initial conditions (t, η) to a prescribed initial condition (t 0 , ξ 0 ) must satisfy t < t 0 , t ր t 0 , and v(η) > v(ξ 0 ), v(η) ց v(ξ 0 ). Indeed, for any initial condition (t 0 , ξ 0 ) the control u defined as u(r) = u max ✶ [t,t+(vmax-v(ξ 0 ))/umax] (r) leads either to extinction or to the set C ∪ V, but it does result a discontinuous value function. We conjecture a Dynamic Programming Principle for this problem. This is an ongoing work. This is an ongoing work, and the next steps for this problem are to derive the Hamilton-Jacobi-Bellman (HJB) equation satified by the value function, to establish if the obtained equation characterizes the value function, that is, if the value function is the unique viscosity solution of the HJB equation. This is a non-standard problem, because we work in an unbounded domain, with a differential operator that is degenerate and not elliptic, and the value function is discontinuous at the boundary of the domain. We expect to be able to give some information about the optimal control and to solve numerically the Hamilton-Jacobi-Bellman equation of the problem. Also, we expect to be able to prove the existence of solutions of the optimal control problem (Filippov's like theorem).

• We intend to extend the study of the stochastic bioreactor model to the case of the chemostat. It will be interesting to compare a stochastic model of the chemostat obtained by means of the methods described in this thesis, taking into account demographic stochasticity, with those existing in the literature. We are also interested in the study of the optimal control problem of maximizing the probability of reaching a target level of depollution prior to the extinction of the biomass, and to derive a dynamic programming principle for this problem.