R ÉSUM É

La digestion anaérobique est un processus biologique au cours duquel des microorganismes décomposent de la matière organique pour produire du biogaz (dioxyde de carbone et methane) qui peut être utilisé comme source d'énergie renouvelable. Cette thèse porte sur l'élaboration de stratégies de contrôle et la conception de bioréacteurs qui maximisent la production de biogaz. La première partie se concentre sur le problème de contrôle optimal de la maximisation de la production de biogaz dans un chemostat avec un modèle à une réaction, en contrôlant le taux de dilution. Pour le problème à horizon fini, nous étudions des commandes type feedback, similaires à ceux utilisés en pratique et consistant à conduire le réacteur vers un niveau de substrat donné et à le maintenir à ce niveau. Notre approche repose sur une estimation de la fonction valeur inconnue en considérant différentes fonctions de coût pour lesquelles la solution optimale admet un feedback optimal explicite et autonome. En particulier, cette technique fournit une estimation de la sous-optimalité des régulateurs étudiés pour une large classe de fonctions de croissance dépendant du substrat et de la biomasse. À l'aide de simulations numériques, on montre que le choix du meilleur feedback dépend de l'horizon de temps et de la condition initiale. Ensuite, nous examinons le problème sur un horizon infini, pour les coûts moyen et actualisé. On montre que lorsque le taux d'actualisation tends vers à 0, la fonction valeur du problème actualisé converge vers la fonction valeur pour le coût moyen. On identifie un ensemble de solutions optimales pour le problème avec coût moyen comme étant les contrôles qui conduisent le système vers un état qui maximise le débit de biogaz sur un ensemble invariant. Nous revenons ensuite au problème à horizon fini fixe et avec le Principe du Maximum de Pontryagin, on montre que le contrôle optimal a une structure bang arc singulier. On construit une famille de contrôles extrémaux qui dépendent de la valeur constante du Hamiltonien. En utilisant l'équation de Hamilton-Jacobi-Bellman, on identifie le contrôle optimal comme étant celui associé à la valeur du Hamiltonien qui satisfait une équation de point fixe. On propose ensuite un algorithme pour calculer la commande optimale en résolvant cette équation de point fixe. On illustre enfin cette méthode avec les deux principales types de fonctions de croissance de Monod et Haldane. Dans la deuxième partie, on modélise et on étudie l'impact de l'hétérogénéité du milieu réactionnel sur la production de biogaz. Pour cela, on introduit un modèle de bioréacteur pilote qui décrit les caractéristiques spatiales. Ce modèle tire parti de la géométrie du réacteur pour réduire la dimension spatiale de la section contenant un lit fixe et, dans les autres sections, on considère les équations 3D de Navier-Stokes en régime permanent pour la dynamique des fluides. Pour représenter l'activité biologique, on utilise un modèle à deux réactions et pour les substrats, des équations advection-diffusion-réaction. On considère seulement les biomasses qui sont attachées au lit fixe et on modélise leur croissance avec une fonction densité dépendante. On montre que ce modèle peut reproduire le gradient spatial de données expérimentales et permet de mieux comprendre la dynamique interne du réacteur. En particulier, les simulations numériques indiquent qu'en mélangant moins, le réacteur est plus efficace, élimine plus de matières organiques et produit plus de biogaz.

i RESUMEN La digestión anaeróbica es un proceso biológico en el cual diferentes poblaciones microbianas transforman compuestos orgánicos en biogás (dióxido de carbono y metano), el cual puede ser luego utilizado como fuente de energía renovable. Esta tesis analiza distintas estrategias de control y diseño de bio-reactores que maximicen la producción de biogás. La primera parte se enfoca en el problema de control óptimo para maximizar la producción de biogás en un bio-reactor continuo o Quimiostato. Se considera el modelo de una reacción y la tasa de dilución es la variable de control. Para el problema con un horizonte finito, se estudia controles retro-alimentados (o tipo feedback), similares a los utilizados en la práctica, y que consisten en llevar el reactor hacia un nivel de sustrato determinado y mantenerlo allí. Nuestro enfoque se basa en establecer límites de la función valor considerando diferentes funciones de costo para las cuales la solución óptima admite una forma explícita, del tipo feedback, independiente del tiempo. En particular, esta técnica proporciona límites explícitos para la sub-optimalidad de los controles estudiados para una amplia clase de funciones de crecimiento dependientes de sustratos y biomasa. A continuación, consideramos el problema con horizonte infinito, tanto para un costo promedio como para uno descontado. Cuando la tasa de descuento tiende a cero, probamos que la función valor del problema descontado converge y que el límite es igual a la función valor para el costo promedio. Luego, se muestra que los controles óptimos para el problema con costo promedio son los que llevan al sistema a un estado que maximiza el flujo de biogás en un conjunto invariante. Posteriormente, volvemos al problema con horizonte finito dado y, usando el Principio Máximo de Pontryagin, demostramos que el control óptimo tiene una estructura bang -arco singular y somos capaces de construir una familia de controles parametrizadas por el valor constante del Hamiltoniano. Usando la ecuación de Hamilton-Jacobi-Bellman, el control óptimo se identifica como el asociado con el valor del Hamiltoniano que satisface una ecuación de punto fijo. A continuación, se propone un algoritmo para determinar el control óptimo mediante la resolución de esta ecuación de punto fijo. En la segunda parte se estudia el impacto de la heterogeneidad del medio en la producción de biogás. Este bioreactor se divide en tres secciones, siendo sólo la intermedia la que contiene biomasa. En dicha sección, el modelo matemático propuesta da cuenta de la geometría del reactor y reduce la dimensión espacial a una sola. Por otro lado, en las otras secciones, las ecuaciones 3D de Navier-Stokes son utilizadas para modelar la dinámica de fluidos. Para representar la actividad biológica se utiliza un modelo de dos reacciones y para los sustratos se utilizan ecuaciones de adveccióndifusión-reacción. Como ya establecido, sólo consideramos la biomasa que está fijada en la sección intermedia y modelamos su crecimiento con una función densidad dependiente. Hemos demostrado que nuestro modelo para este bio-reactor reproduce adecuadamente el gradiente espacial de datos experimentales y proporciona una mejor comprensión de la dinámica interna del reactor. En particular, las simulaciones numéricas indican que al mezclar menos, el reactor es más eficiente y produce más biogás. ABSTRACT Anaerobic digestion is a biological process in which organic compounds are degraded by different microbial populations into biogas (carbon dioxyde and methane), which can be used as a renewable energy source. This thesis works towards developing control strategies and bioreactor designs that maximize biogas production. The first part focuses on the optimal control problem of maximizing biogas production in a chemostat in several directions. We consider the single reaction model and the dilution rate is the controlled variable. For the finite horizon problem, we study feedback controllers similar to those used in practice and consisting in driving the reactor towards a given substrate level and maintaining it there. Our approach relies on establishing bounds of the unknown value function by considering different rewards for which the optimal solution has an explicit optimal feedback that is time-independent. In particular, this technique provides explicit bounds on the sub-optimality of the studied controllers for a broad class of substrate and biomass dependent growth rate functions. With numerical simulations, we show that the choice of the best feedback depends on the time horizon and initial condition. Next, we consider the problem over an infinite horizon, for averaged and discounted rewards. We show that, when the discount rate goes to 0, the value function of the discounted problem converges and that the limit is equal to the value function for the averaged reward. We identify a set of optimal solutions for averaged problems as the controls that drive the system towards a state that maximizes the biogas flow rate on an special invariant set. We then return to the problem over a fixed finite horizon and with the Pontryagin Maximum Principle, we show that the optimal control has a bang singular arc structure. We construct a one parameter family of extremal controls that depend on the constant value of the Hamiltonian. Using the Hamilton-Jacobi-Bellman equation, we identify the optimal control as the extremal associated with the value of the Hamiltonian which satisfies a fixed point equation. We then propose a numerical algorithm to compute the optimal control by solving this fixed point equation. We illustrate this method with the two major types of growth functions of Monod and Haldane. In the second part, we investigate the impact of mixing the reacting medium on biogas production. For this we introduce a model of a pilot scale upflow fixed bed bioreactor that offers a representation of spatial features. This model takes advantage of reactor geometry to reduce the spatial dimension of the section containing the fixed bed and in other sections, we consider the 3D steady-state Navier-Stokes equations for the fluid dynamics. To represent the biological activity, we use a 2 step model and for the substrates, advection-diffusionreaction equations. We only consider the biomasses that are attached in the fixed bed section and we model their growth with a density dependent function. We show that this model can reproduce the spatial gradient of experimental data and helps to better understand the internal dynamics of the reactor. In particular, numerical simulations indicate that with less mixing, the reactor is more efficient, removing more organic matter and producing more biogas.

Chapter 1 Introduction 1.1 Anaerobic Digestion and Biogas Production

Anaerobic Digestion (AD) is a natural process of degradation of organic matter by microorganisms in the absence of oxygen. It is the result of successive steps performed by different groups of microorganisms, progressively degrading matter from large organic polymers into simple monomers (sugars, amino acids and fatty acids) and further into biogas (methane and carbone dioxyde) [START_REF] Gujer | Conversion processes in anaerobic digestion[END_REF].

The interest in this process is that it allows the re-valorization of waste as a renewable energy source since the methane in biogas can be used as fuel for heating or producing electricity. This reduces greenhouse gas emissions by producing an alternative to fossil fuels and by capturing the emissions of biodegrading matter that would otherwise be released to the atmosphere. Furthermore, biogas presents a potential for electrical grid balancing since it can be stored and therefore the power generation is controllable, unlike other renewable energy sources, such as solar and wind. An additional environmental benefit is the production of digestate, the matter not consumed by microorganisms (minerals, certain organic compounds such as structural plant matter including lignin and cellulose) and this can be used as fertilizer.

This process is implemented in anaerobic digesters, a type of bioreactor, and these devices are designed to offer the best conditions to support microorganisms and maintain the process. Dedicated AD plants have been developed to process agricultural waste, manure and energy crops, but biogas can also be recovered from landfills or in wastewater treatment, during the removal of organic matter [START_REF] Metcalf | Wastewater engineering; treatment and reuse[END_REF].

The current installations in Europe and North America are mainly large scale electricity and heat biogas plants but in Asia and Africa, many small domestic scale digesters can be found in rural areas [START_REF]Renewable Energy Policy Network for the 21st Century, Renewables 2016 global status report[END_REF]. For the European Union, studies [START_REF] Scarlat | Biogas: Developments and perspectives in europe[END_REF] have underlined that despite a high potential, production remains low. Biogas represented 4.4% of natural gas consumption in 2015, however, this average hides a heterogeneous situation : this figure can be as high as [START_REF] Crespo | Modeling and optimization of activated sludge bioreactors for wastewater treatment taking into account spatial inhomogeneities[END_REF].2% in Sweden while in France it is only 1.5%. In Chile, again a high capacity is expected due to the importance of the agricultural sector and an estimated 4.1 TWh could be produced [START_REF] Rodríguez-Monroy | Electricity generation in chile using non-conventional renewable energy sources -a focus on biomass[END_REF] for a total electricity consumption in 2016 of 67 TWh [22], but only 2% of the potential is achieved.

Various factors could explain this situation, despite the fact that the necessary technology has been available for several decades now. First, a high starting investment is required as large AD bioreactors are needed to process the important volumes in wastewater treatement or to produce a significant amount of biogas. Secondly, there is a serious difficulty in assessing the economic viability of a project due to the volatility of energy markets and, in addition, estimating the production of a bioreactor in advance remains a challenge. Finally, there has been reports of poor profits due to high costs and low yields that have lead to the closing of biogas plants [START_REF] Beddoes | An analysis of energy production costs from anaerobic digestion systems on us livestock production facilities[END_REF].

For wastewater treatment, the use of AD for the removal of organic matter is still uncommon and the alternative method, aerobic digestion, has prevailed for a long time. Nonetheless, AD presents several advantages compared to aerobic digestion: less sludge is produced, treatment of water with higher organic loadings is possible and the production of biogas reduces the energetic cost of the process [START_REF] Riffat | Fundamentals of wastewater treatment and engineering[END_REF]. However, in addition to the problems previously mentioned, AD requires a longer time to start up, needs to be maintained at high temperature and is susceptible to toxic substances. This last point, and in general the instability of the process is an important aspect of AD when considering its implementation. Indeed, this process can be inhibited by wide range of substances, either that enter the reactor such as pathogens, or that are intermediate products of the process, specifically volatile fatty acids (VFA) that can significantly increase the pH [START_REF] Chen | Inhibition of anaerobic digestion process: A review[END_REF]. Another reason explaining the complexity of the process is the important differences between microorganisms producing acids and those producing methane [START_REF] Pohland | Developments in anaerobic stabilization of organic wastes -the two-phase concept[END_REF]. The consequence of this instability and the complexity of microbial ecosystems and bioprocesses, is that when developing implementations of AD, the stability of the process is considered first and foremost, often neglecting optimization of production.

Mathematical Models of Bioreactors

Bioreactors are used for science and industry, for the study of microorganisms, cell and tissue culture or production of derivatives and end products for chemical processes. The first models were introduced in the 1950s when the chemostat device was invented independently by Jacques Monod [START_REF]La technique de culture continue: theorie et applications[END_REF] and Aaron Novick and Leo Szilard [START_REF] Novick | Experiments with the chemostat on spontaneous mutations of bacteria[END_REF]. Since then, a large variety of models have been developed [START_REF] Vanrolleghem | Bioprocess model identification, in Advanced instrumentation, data interpretation, and control of biotechnological processes[END_REF][START_REF] Wade | Perspectives in mathematical modelling for microbial ecology[END_REF]. The purpose of this section is not to review the wide range of models but rather to present the key properties of the models used in this thesis.

The key aspect of microbial ecosystems is the important diversity of microorganisms, in terms of physiology, nutritional needs, growth kinetics, and sensitivity to environmental conditions. As a consequence, population models become rapidly complicated and their mathematical analysis extremely difficult. However, we are interested in the macroscopic behaviour of a bioreactor and in particular how matter is transformed. From this point of view, a bioreactor can be modelled as the combination of a biological process, the result of the microbial activity, and a physical process, the dynamics of the substances inside the reactor and how they enter and leave. This approach relies on characterizing microorganisms not as a population but instead representing the biological process essentially as a bio-chemical reaction. The basis of the model is then mass conservation and modelling the biological process corresponds to accounting how matter is transformed from one substance into another. The considered variables are the concentrations of reacting matter, in units of mass per volume. This includes microorganisms and they are characterized as biomass, which can be measured although it is a complicated task, in addition to being costly and imprecise. When the bioprocess is the result of a complex ecosystem, it can be viewed as the sum of several reactions, grouping microorganisms by functional role and how they transform matter.

An important property of these bio-reactions is that the microorganisms are at the same time a catalyst and a product : they both consume reactants and grow as a result. The usual hypothesis of biological activity is that consumption and production are proportional to growth, so that the reaction rate of each substance is obtained from a growth rate and a yield coefficient.

In the simplest cases, when it can be considered that microorganisms have equal access to the substrate they consume, the growth rate is considered linear with respect to biomass concentration. The dependence on the substrate consumed is however more complex, non linear and represented by a specific growth rate function, the rate of growth per unit of biomass. It generally depends on a single limiting substrate but eventually other factors affecting the microorganisms can be taken into account (inhibition by other substances, impact of physical or chemical properties of medium).

A wide range of specific growth rate functions have been considered and the Monod function [START_REF] Monod | The growth of bacterial cultures[END_REF] was one of the first to be proposed. It models the fact that the reaction can take place up to a maximum rate µ max . Here, s denotes the concentration of limiting substrate and K s the half-saturation constant,

µ M (s) = µ max s K s + s .
The Haldane function can be used to represent an inhibition of the reaction by the substrate, with K i the inhibition parameter,

µ H (s) = μ s K s + s + s 2 K i .
To represent crowding effects, for example when microorganisms do not have equal access to the substrate, the Contois function can be considered, which also depends on the biomass concentration

x, µ C (s, x) = µ max s K x x + s .
It is important to note that this remains an empirical approach, and that the parameters of growth functions and yield coefficients are obtained by fitting the model to measurements. In particular, for AD this needs to be done for each bioreactor, as the microbial community adapts itself to the feedstock. This has lead to the development of many different models of anaerobic digestion. The current standard, the ADM1 model [START_REF] Batstone | The IWA anaerobic digestion model no 1 (ADM1)[END_REF], was developed by the International Water Association (IWA) Anaerobic Digestion Modelling Task Group, and is composed of 19 reactions and at least 32 dynamic concentration state variables. Such large dimensions make this model impractical for analytical optimization or control, so simplified models have been proposed. A widely used model [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF], is composed of two biochemical processes to take into account the inhibition of the methane production by intermediate products, namely volatile fatty acids. Further simplification is possible, as it has been shown that a single step model can reproduce the qualitative behaviour of the anaerobic digestion process [START_REF] Bernard | Can we assess the model complexity for a bioprocess: theory and example of the anaerobic digestion process[END_REF]. Indeed, for the two-step model, the second reaction is the most limiting due to inhibition by the substrate and then a one-step model can be used to focus on the second reaction. In particular, a common assumption is to consider that the first step is fast and then the two reactions can be reduced to a single one with a slow-fast approximation, in which case it provides a good representation of the biogas production.

Monod Haldane

Modelling the physical part of the process consists in taking into account the design features of a bioreactor that have an impact on the biological activity. The first aspect to consider is how the reactor is fed and how products are retrieved. The 3 main modes of operation are :

• Batch : the contents are introduced at the beginning of the process, after which nothing is added or removed.

• Fed-batch : similar to the batch mode, although susbtrate can be added during process operation, but nothing is removed until the end.

• Continuously fed : during all of the process operation contents are added and removed from the reactor at the same rate such that the volume remains constant.

The last case is the most common for wastewater treatment and anaerobic digesters and the one we will consider in this thesis.

An important question is then if the reacting medium is homogeneous or heterogeneous. The first case occurs when the reacting contents are mixed and concentrations can be considered constant throughout the bioreactor such that the mass balance corresponds to a system of Ordinary Differential Equations (ODE). The basic example is the chemostat, a continuously fed, well mixed bioreactor [START_REF] Harmand | The Chemostat: Mathematical Theory of Microorganisms Cultures[END_REF][START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF]. For a single step biological processes, denoting s the substrate concentration and x the biomass, the model equations are

ṡ = D(s in -s) - 1 Y µ(s)x, ẋ = µ(s)x -Dx,
with s in the input concentration, µ(•) the specific growth rate, Y the yield coefficient, and D the dilution rate which is equal to the feeding rate divided by the constant volume.

However, there is a diversity of designs that take advantage of a spatial gradient in the concentrations, such as plug-flow reactors. In this case, bioreactors can be modeled using Partial Differential Equations (PDE) [START_REF] Samstag | Cfd for wastewater treatment: an overview[END_REF] but this leads to complex models, especially when authors try to account for all phenomenons (turbulence, different phases of matter, sedimentation,...). In particular, anaerobic digesters are generally very large since AD is a slow process and important volumes are treated, so that simulating an industrial scale bioreactor requires substantial computational resources. An alternative is compartment models, where the reactor is seen as a network of interconnected homogeneous zones, each of which is modeled as a chemostat, resulting in a system of ODE, although of large dimension.

Optimization of Biogas Production

The topic of this thesis is the maximization of biogas production, through the study of mathematical models and optimization problems. The first part will focus on bioreactor operation and the problem considered is the optimal control of the dilution rate for the chemostat. In a second part, a question of bioreactor design is investigated to understand the impact of heterogeneity of the reacting medium on biogas production.

It is important to point out that the aspect of stability of the bioprocess will not be directly considered in the problems addressed. Since biogas is a final product, maximizing it guarantees, although indirectly, that the process is kept in a healthy state. However, this reasoning is limited and the extreme nature of an optimum could result in running the process close to an unstable state, with the risk of tipping over inevitable due to the variability of bioprocesses.

In addition to discarding the issue of stability, we are looking at theses questions of bioreactor operation and design as mathematical problems and it is important to remember that a mathematical model is only an imperfect representation. There are also important sources of model errors, as we have seen in the previous section, rendering the models essentially qualitative rather than predictive.

The consequence is that we can not expect our results to be directly applicable. Instead, the objective here is to find a qualitative description of the optimal solution in order to understand the levers that can help increase biogas production. For instance, with optimal control problems, the aim is not to obtain controls to be used in practice but rather to find what characterizes the optimal control, to help practitioners improve their strategies in terms of biogas production. Furthermore, we are working with simple macroscopic models of an extremely variable microbial ecosystem and a key weakness of this modelling approach is the choice of a growth function. To remedy this, we want to obtain results for a general class of growth functions, characterized by general properties of microbial dynamics (for example, to generalize the Monod function, consider all functions that are increasing and bounded).

Ideally, for this we want to find the explicit or analytical expression of the solution of the optimization problems considered. However, the analytical resolution of optimization problems poses a serious challenge. For example, a 2 stage model of AD will have at least 4 dynamic variables and therefore carrying out the computations to obtain an analytical expression of an optimal control rapidly become complex.

On the other hand, numerical resolution of optimization problems is now well established and the progress in computer performance in the last decades mean that it is possible to solve large problems numerically. For the optimal control problems in consideration here, solvers such as BOCOP [START_REF] Bonnans | Bocop -A collection of examples[END_REF] can easily and rapidly compute an optimal control. However, it is obtained in open loop form whereas we aim at finding feedbacks, which are more robust in terms of possible time delays and measurement errors. More generally, the problem in using numerical optimization is that the computed solution is valid only for one set of parameters of a given growth function.

The consequence is that neither of these approaches alone, purely analytical or purely numerical, is adequate and this is a major challenge in the application of mathematical optimization for the development of bioprocesses. In this thesis, we propose several new methods that address this problem.

In Section 1.3.1, corresponding to Chapters 2 and 3, we present the main contributions of this thesis for the optimal control of the dilution rate. We work with a simple model to develop new techniques for the analysis of optimal control problems, that could be used with more complex models. We will also see how to carefully combine both analytical and numerical approaches. This will allow us to obtain expressions of optimal and sub-optimal feedback controls and practical ways to compute them and analyze their performance. In Section 1.3.2, corresponding to Chapter 4, dealing with bioreactor design, we show that it is not necessary to develop overly complex models to study spatial heterogeneity in bioreactors and we will see how numerical simulations can be used to obtain qualitative results.

Optimal Control of Biogas Production

Control of bioreactors is an important and active research topic and more specifically, control of AD has been the focus of many studies [START_REF] Dochain | Bioprocess control[END_REF][START_REF] Gaida | Feed control of anaerobic digestion processes for renewable energy production: A review[END_REF][START_REF] Holubar | Advanced controlling of anaerobic digestion by means of hierarchical neural networks[END_REF][START_REF] Jimenez | Instrumentation and control of anaerobic digestion processes: a review and some research challenges[END_REF][START_REF] Junker | Bioprocess monitoring and computer control: key roots of the current pat initiative[END_REF][START_REF] Nguyen | Automatic process control in anaerobic digestion technology: A critical review[END_REF]. The feeding rate is typically considered as the variable input, as it has a rapid and significant impact on the process. Other inputs have been considered, such as pH or alkalinity, but they are less cost effective as they require adding substances (such as a concentrated acid or base solution) [START_REF] Pind | Monitoring and control of anaerobic reactors[END_REF].

The various strategies developed in process control have been applied for the control of AD and they can be divided in 2 categories : model based and knowledge based. The latter are mainly proposed by experts, such as bioprocess engineers, who construct strategies based on practical knowledge. On the other hand, model based controllers are developed by the automatic control community and the main advantage of this approach is the theoretical properties that they can guarantee, such as robustness or performance [START_REF] Bonnans | Commande et optimisation de systemes dynamiques[END_REF].

For knowledge based controllers, recent works have incorporated the aspect of optimizing performance, but avoiding failure of the process is still prioritized and only when the process is stable does the controller attempt to push the system towards higher biogas production [START_REF] Steyer | Advanced control of anaerobic digestion processes through disturbances monitoring[END_REF]. For example, the feedback controller developed in [START_REF] Rodrıguez | A hydrogenbased variable-gain controller for anaerobic digestion processes[END_REF] uses a set-point for the hydrogen outflow concentration, since it is very sensitive to process destabilization and therefore guarantees a fast response of the controller. To optimize the biogas production the gain is dependent on the methane concentration and allows the controller to push the system to higher dilution rates as long as hydrogen remains close to the set-point level. The PID (Proportional Integral Derivative) controller developed in [START_REF] García-Diéguez | Multi-objective cascade controller for an anaerobic digester[END_REF] uses a cascade strategy consisting of 2 control loops: the first is used to control the VFA concentration, to avoid the acidification of the reactor and guarantee process stability. The second loop then attempts to optimize production by driving the system towards a set-point that maximizes the theoretical methane outflow. These knowledge based controllers have the advantage to be simple, robust and effective although they are dependent on the gain parameters and finding the optimal values can be challenging.

With model based control, maximization of biogas production has been taken into account more directly but the optimization is generally of static nature in the sense that controls are designed to drive the process towards a steady state, computed from a model, that maximizes biogas production. For example, the equilibriums of the 2 step model have been computed and the dilution rate corresponding to the most productive steady state can be found by solving an optimization problem [START_REF] Bayen | On the steady state optimization of the biogas production in a two-stage anaerobic digestion model[END_REF] or extremum seeking algorithms can be used to reach it [START_REF]Nonlinear adaptative stabilizing control of an anaerobic digestion model with unknown kinetics[END_REF][START_REF] Lara-Cisneros | On the dynamic optimization of methane production in anaerobic digestion via extremum-seeking control approach[END_REF]. The few applications of optimal control theory also follow this direction: the minimal time problem of reaching an optimal steady state has been solved on a 2 dimensional invariant set of a 2 stage model [START_REF] Bayen | Analysis of an optimal control problem related to the anaerobic digestion process[END_REF].

In contrast, there has been much less work considering the dynamic optimization problem over the transients. In [START_REF] Sbarciog | On the optimization of biogas production in anaerobic digestion systems[END_REF], the authors study the problem of driving the system to a neighborhood of an optimal steady state of a 2 step model, while maximizing the biogas outflow rate but with a penalization of the control. The Pontryagin Maximum Principle is used to show that the optimal control is of bang-bang type and although the switching is not explicit, this has lead to the development of a control with a heuristic switching strategy based on methane measurements [START_REF] Sbarciog | A biogas-based switching control policy for anaerobic digestion systems[END_REF].

We will consider here the optimal control problem of maximizing the total methane produced for the single step model of the chemostat. For a fixed final time

T ∈ R, initial data ξ := (t 0 , s 0 , z 0 ) ∈ D := (-∞, T ) × [0, s in ) × (0, ∞)
and maximum dilution rate

D max > 0, the problem is Maximize J(ξ, D(•)) := T t 0 µ(s(t))x(t) dt over all D : [t 0 , T ] → [0, D max ] measurable such that ṡ = D(s in -s) -µ(s)x, s(t 0 ) = s 0 , ẋ = µ(s)x -Dx, x(t 0 ) = x 0 . (P bio )
This problem was first stated over 20 years ago [START_REF] Stamatelatou | Optimal and suboptimal control of anaerobic digesters[END_REF] and using the Pontryagin Maximum Principle, it was shown that the optimal control is bang-bang singular arc but the optimal synthesis was not given and remains unknown today.

Ghouali et al. [START_REF] Ghouali | Maximizing biogas production from the anaerobic digestion[END_REF] have made progress by solving (P bio ) for initial conditions in the invariant set I = {x + s = s in }, for which the dynamics reduce to a single scalar equation This is achieved for a large class of growth functions, only requiring the existence of a unique maximum on the set I of the growth rate µ(s)(s in -s), or equivalently, monotonicity on either side of the maximum. This allows to use a comparison result for scalar ODEs to compare trajectories associated with different controls and then establish which achieves the greatest production, since the biogas flow rate is proportional to the growth rate. The optimal control is a most rapid approach path (MRAP) feedback to the unique maximizer s

ṡ = D -µ(s) (s in -s).
ψ s(s) = 0 if s > s, µ(s) if s = s, D max if s < s. (1.1)
Optimal control problems over a fixed time horizon, in general, possess a time-dependent optimal synthesis, while the duration of process operation is often poorly known. The reduced problem exhibits the remarkable feature of having an optimal synthesis independent of the terminal time, which makes it attractive from an application point of view. However, this will not be the case for the general problem (P bio ) and as we shall see later on, the optimal control will depend on the time horizon and initial condition considered. This makes it considerably more difficult to solve and in addition, as the dynamics are 2 dimensional, we can not use the same comparison technique.

The feedback (1.1) has the additional benefit that it can be implemented without complex biomass measurements, and only requires estimations of substrate concentration and biogas flow rate. Finally, this control strategy of driving the process to an optimal substrate level resembles those used in practice [START_REF] García-Diéguez | Multi-objective cascade controller for an anaerobic digester[END_REF][START_REF] Steyer | Advanced control of anaerobic digestion processes through disturbances monitoring[END_REF]. For these reasons, we have developed a method of estimating the sub-optimality for all initial conditions of similar MRAP feedbacks to any substrate level s * ∈ [0, s in ], and identifying the best one.

Sub-Optimal Feedback Controls with Guaranteed Value

The starting point of this work is the observation of monotonicity properties arising from different sources. The first comes from natural assumptions on the growth rate and we will consider here substrate and biomass dependent functions. Assumption 1.1 We suppose that µ : R + × R + → R + is a Lipschitz continuous function that satisfies, for all x > 0 µ(0, x) = 0 and µ(s, x) > 0 for s > 0.

We suppose as well that x → µ(s, x) is non increasing, which models crowding effects, and x → µ(s, x)x is non decreasing, which models the fact that having more biomass provides at least the same growth.

A typical instance of this class is the Contois growth function, but this class of functions also contains growth functions that depend only on the substrate concentration, such as the Monod and the Haldane functions.

The next monotonic behaviour can be observed in the chemostat model, in the way trajectories are attracted to the invariant domain I, independently of the control. This is immediately apparent with the change of variables

z = x s in -s , ż = µ s, (s in -s)z (1 -z)z
The consequence is that we can use the fact that z(•) is monotonic for all controls, in conjunction with the properties of the growth rate to establish a relation between the value functions of (P bio ) and of an auxiliary problem with an objective function depending only on s(•). To simplify notations, we denote φ(s, z) = µ s, (s in -s)z (s in -s) and the growth rate is then φ(s, z)z. We now consider the problem, for z 1 0,

Maximize J z 1 (ξ, D(•)) := T t 0 φ(s(t), z 1 ) dt over all D : [t 0 , T ] → [0, D max ] measurable such that ṡ = D(s in -s) -φ(s, z)z, s(t 0 ) = s 0 , ż = µ s, (s in -s)z (1 -z)z, z(t 0 ) = z 0 . (P z 1 )
We denote the value function of (P bio ) as V (•) and of (P z 1 ) as W z 1 (•), that we can relate in the following manner. Proposition 1.2 For any initial data ξ ∈ D and any z 1 ∈ [min(z 0 , 1), max(z 0 , 1)], we have the following frame for the value function V of the original problem

min(z 0 , 1)W z 1 (ξ) V (ξ) max(z 0 , 1)W z 1 (ξ). (1.2)
Then, any optimal control u z 1 (•) of the auxiliary problem (P z 1 ) guarantees a (sub-optimal) value for the original criterion J(ξ, •) that satisfies

min(z 0 , 1)W z 1 (ξ) J(ξ, u z 1 (•)) max(z 0 , 1)W z 1 (ξ)
and we have the following sub-optimality estimation

V (ξ) -J(ξ, u z 1 (•)) |1 -z 0 |W z 1 (ξ). (1.
3)

The first frame (1.2) is remarkable as it gives an estimation of the value function, without having to solve the associated problem.

On the other hand, the sub-optimality estimation (1.3) will be the basis for assessing the performance of MRAP feedbacks. It turns out that the optimal controls of (P z 1 ) are similar to the solution (1.1) of the reduced problem. Indeed, the objective function J z 1 (•) of (P z 1 ) does not depends directly on z(•) and since this variable is monotone, the impact of a control is seen through the s(•) variable. Therefore, we can employ a similar comparison technique, as used for the reduced problem, to compare the rewards associated with different controls and prove the optimality of an MRAP feedback. For this, we need again to assume the existence of a unique maximum, but this time for the growth rate on the set {z = z 1 }. Assumption 1.3 For each z 1 0, the function s → φ(s, z 1 ) admits a unique maximum on (0, s in ), and we denote the substrate level at which this maximum is attained as s(z 1 ) = arg max s∈(0,s in ) φ(s, z 1 ).

We then have the following result.

Proposition 1.4 For all initial data ξ ∈ D, the MRAP feedback to s(z 1 ),

ψ s(z 1 ) (s, z) = 0 if s > s(z 1 ), µ(s(z 1 ), (s in -s(z 1 ))z) z if s = s(z 1 ), D max if s < s(z 1 ), (1.4)
is optimal for the auxiliary problem (P z 1 ).

In the case of a substrate only dependent growth function, these feedbacks all coincide with the MRAP to s(1), the maximum on the invariant set I = {z = 1}, since φ(s, z) = µ s (s in -s). However, for the substrate and biomass dependent case, these controls can differ and various factors impact which feedback is the best.

We can expect the initial condition to have an influence on the performance of (1.4), as can been seen with the sub-optimality estimation (1.3). In particular, due to the the term |1 -z 0 |, the distance to the set I will be important. However, through W z 1 (•), the impact of the initial condition will differ for each growth function. Indeed, in Figure 1.3 we can observe that, for the Contois growth function, W z 1 (•) varies significantly with the initial biomass and this can be attributed to the dependence of the Contois growth function on biomass concentration.

On the other hand, the impact of the time horizon is more straightforward. Since the trajectories are attracted to the invariant set, we can expect ψ s(1) (•) to be the best when the horizon is sufficiently large. However, when the horizon is small, the feedback ψ s(z 0 ) would seem to be the best option since this strategy consists in remaining close to the maximum of the biogas flow rate corresponding to the initial condition, whereas another feedback could drive the system away, towards another maximizing state but that can not be reached in time.

To examine this we proceed with numerical simulations by computing the reward for a range of values of z 1 ∈ [min(z 0 , 1), max(z 0 , 1)] and of final times for fixed initial data. In order to identify the maximum of J(ξ, ψ s(z 1 ) (•)) with respect to z 1 for different final times, we normalize the reward by computing

J N (T, z 1 ) = J(ξ, ψ s(z 1 ) (•)) -min y J(ξ, ψ s(y) (•)) max y J(ξ, ψ s(y) (•)) -min y J(ξ, ψ s(y) (•))
where the minimum and maximum are taken for y ∈ [min(z 0 , 1), max(z 0 , 1)]. Hence, for each final time T , the maximum reward is achieved for z 1 such that J N (T, z 1 ) = 1 and the minimum when J N (T, z 1 ) = 0. Figure 1.4 clearly shows that there is a monotonous relation between the final time and the best feedback: the greater the time horizon, the closer z 1 must be chosen to the invariant set I. This has lead us to examine the case of a infinite time horizon.

Infinite Horizon Problems

Stating the optimal control problem (P bio ) on a finite horizon raises a number of issues. In practical applications, it can be difficult or even impossible to specify a final time in advance, especially considering that, in general, solutions of finite horizon problems depend on the given time interval and therefore any change mid-course of the planning horizon will result in loss of optimality. For a long time now, researchers working on optimization related to economics have dealt with these difficulties by considering problems over an infinite horizon [START_REF] Kamien | Dynamic optimization: the calculus of variations and optimal control in economics and management[END_REF][START_REF] Seierstad | Optimal control theory with economic applications[END_REF]. Such a formulation of optimal control problems also reflects the need for preserving the viability of a system indefinitely.

As the reward is unbounded on a infinite horizon, we have to choose a concept of optimality [START_REF] Carlson | Infinite horizon optimal control: deterministic and stochastic systems[END_REF]. In most cases, the process is operated for a very long duration and the performance expected from the practitioners is to maintain a high average value over time, so a natural choice is to consider the limit of the averaged reward

J T (ξ, D(•)) = 1 T T 0 µ s(t), x(t) x(t) dt. (1.5)
However, we have to consider the inferior and superior limit as it is possible to construct controls for which the average reward does not converge. We thus consider two optimal control problems and denote their value functions

V ∞ (ξ) = sup lim inf T →∞ J T (ξ, D(•)) : D(•) ∈ L ∞ ((0, ∞), [0, D max ]) , (1.6)
V ∞ (ξ) = sup lim sup T →∞ J T (ξ, D(•)) : D(•) ∈ L ∞ ((0, ∞), [0, D max ]) . (1.7)
Another choice of objective function, often used in problems related to economics, is the discounted reward, for a discount rate δ > 0,

J δ (ξ, D(•)) = ∞ 0 δe -δt µ s(t), x(t) x(t) dt. (1.8)
The term e -δt represents a discount rate or a preference for earlier rather than later production. For a positive discount rate, the optimal control for this reward could be very different from the solution of problems (1.6) and (1.7) but, when the discount rate δ goes to 0, the average and discounted rewards are in fact related. This is the reason for rescaling the integral in (1.8) with the discount rate δ, in order to guarantee that the limit remains finite. We thus consider the following value function

V δ (ξ) = sup J δ (ξ, D(•)) : D(•) ∈ L ∞ ((0, ∞), [0, D max ]) .
(1.9)

The relation between the value functions (1.6), (1.7) and the limit of (1.9) as δ goes to 0 has been studied by Grüne [START_REF]On the relation between discounted and average optimal value functions[END_REF] and the basis of his work is the following result [41, Lemma 3.1].

Lemma 1.5 Let q : R → R be a measurable and bounded function. If the average

1 T T 0 q(t) dt converges when T → ∞ then lim T →∞ 1 T T 0 q(t) dt = lim δ→0 ∞ 0 δe -δt q(t) dt.
A similar result can be established for objective and value functions of optimal control problems for general control systems and we transcribe it here in our setting.

Lemma 1.6 For all initial data ξ ∈ D, and controls

D(•) ∈ L ∞ ((0, ∞), [0, D max ]), lim inf T →∞ J T (ξ, D(•)) lim inf δ→0 J δ (ξ, D(•)) lim sup δ→0 J δ (ξ, D(•)) lim sup T →∞ J T (ξ, D(•)), and V ∞ (ξ) lim inf δ→0 V δ (ξ) lim sup δ→0 V δ (ξ) V ∞ (ξ).
The immediate consequence is that equality holds in both these frames if the limit of the average reward exists. The main result of Grüne shows that the average and discounted value functions are not only equal but also piecewise constant for general affine control systems satisfying a controllability assumption. However, the work of Grüne gives no information on the optimal controls, whether the solutions of the discounted problems for a positive discount rate might converge to a solution of the limit problem.

In our case, we prove the existence of optimal controls for the discounted problem, which is a first step in showing the convergence of the optimal controls thanks to the concept of Γ-limit, the convergence notion guaranteeing that optimal solutions converge to a maximizer of the limit problem.

Proposition 1.7 For all ξ ∈ D and for all δ > 0, the suprema are attained,

V δ (ξ) = max D(•) J δ (ξ, D(•)).
If the Γ-limit of J δ (•) exists as δ goes to 0,

J 0 (ξ, D(•)) := Γ -lim δ→0 J δ (ξ, D(•)),
then the maxima converge, pointwise in ξ, to the maximum of the limit,

V 0 (ξ) := lim δ→0 V δ (ξ) = max D(•) J 0 (ξ, D(•)).
(1.10)

Furthermore, if D δ (•) is an optimal control for (1.9), i.e. if V δ (ξ) = J δ (ξ, D δ (•)) and if D δ (•) converges to D 0 (•) in L ∞ ((0, ∞), [0, D max ]) then D 0 (•) is an optimal control for (1.10) and V 0 (ξ) = J 0 (ξ, D 0 (•)) = lim δ→0 J δ (ξ, D δ (•)).
Although our problem does not verify the controlability assumption of Grüne, Lemma 1.6 gives valuable information and allows us to show the equality of the value functions (1.6), (1.7) and (1.10). The intuition here is that, under suitable conditions, the limit of the average of a function is its value at infinity. In our case, we can show that the average reward associated with a converging trajectory is equal to the biogas flowrate reached at infinity. Moreover, we can exploit a key property of the chemostat model : the invariant domain I is attractive for all controls 1 and therefore the maximum of the biogas flowrate on the set I is an upper bound of the value functions. We denote s = s(1), a substrate level, not necessarily unique this time, at which the maximum biogas flowrate is attained. Then the average reward of a control that drives the system to s is equal to the upper bound and thus such a control is optimal. Proposition 1.8 For any initial data ξ ∈ D, any D(•) ∈ L ∞ ((0, ∞), [0, D max ]) that drives the system asymptotically to the state (s, s in -s) is optimal for problems (1.6), (1.7) and (1.10). We then have

V ∞ (ξ) = V 0 (ξ) = V ∞ (ξ) = µ(s, s in -s)(s in -s).
This result allows us to prove that the MRAP feedback (1.4) to s is optimal for the average reward problems since this control brings the system towards s and maintains it there as it reaches the set I. However, this is not the only optimal control and, for example, in the case of a growth function that depends only on the substrate and that is monotone (such as the Monod growth function), the constant control D = µ(s) can also drive the system to the state (s, s in -s).

More generally, Proposition 1.8 offers a simple characterization of the optimal controls and makes this approach valuable for understanding what constitutes a control strategy that maximizes biogas production. In addition, posing the problem on an infinite horizon is interesting when considering more complex models due to the relative simplicity of resolution of these types of problems, compared to the finite horizon problem for instance.

Fixed Point Algorithm

We now consider the full problem (P bio ) on a finite horizon. As mentioned previously, the optimal synthesis for this type of problem is in general time dependent and a key difficulty here is understanding the impact of the time horizon on the optimal control. Indeed, when this problem was first studied [START_REF] Stamatelatou | Optimal and suboptimal control of anaerobic digesters[END_REF], although it was shown that the optimal control is bang-bang singular, the dependence of the singular arcs on the initial data was not given. It is in this aspect that we seek to advance and we focus here on proposing a candidate to optimal control in feedback form and giving a practical way to compute it.

To gain information on the singular arcs, we use the Pontryagin Maximum Principle (PMP) [START_REF] Clarke | Functional analysis, calculus of variations and optimal control[END_REF], which states that an optimal control maximizes the Hamiltonian, a function of the state variables and associated adjoint states. An important fact here, is that the Hamiltonian does not depend explicitly on time and thus is constant and equal to some h = h(ξ) ∈ R. Combining this with the maximum condition, we can establish an equation that is valid on a singular arc. Here, in the case of a substrate only dependent growth function, we can view this as the fact that extremal trajectories during the singular arc remain in the graph of

s → x h (s) := h µ (s)(s in -s) µ(s) 2 .
(1.11)

1 except a certain set of controls but that can not be optimal The consequence is that if we know the value of h, we can identify the singular arc and if we know how the control switches from bang to bang and from bang to singular arc, we can construct an admissible extremal control.

We assume here, from the knowledge of the solution on the invariant set I, that the optimal trajectories follow a most rapid approach path to the singular arc, and remain on it. Althought computing a feedback D h (s) depending only on s and h, to stay on the graph of (1.11) is straightforward, there are several challenges to fully construct a control for a large class of growth functions. In particular, it might not be possible to remain on the singular arc with an admissible control and then the graph of (1.11) could be divided in several disjoint admissible sets. Another key point is that bang-singular arc trajectories might not cover all of the state space and therefore a switching curve for bangbang-singular arc trajectories must be determined. However, for the Monod and Haldane growth functions, thanks to further assumptions, we can determine the admissible section G h of the singular arc and identify the set G 0 h (resp. G max h ) on which the control is 0 (resp. D max ). The result is the following feedback, illustrated in Figure 1.5,

ψ h (x, s) =      0, if (x, s) ∈ G 0 h , D max , if (x, s) ∈ G max h , D h (s), if (x, s) ∈ G h . (1.12)
In order to make this a suitable candidate to optimal control, we can identify the value of the Hamiltonian h for a given initial data thanks to the link between the PMP and the Hamilton-Jacobi-Bellman (HJB) equation.

(-h, p x (t), p s (t)) ∈ ∂ + V (t, x(t), s(t)), a.e. on [t 0 , T ].
where ∂ + V is the viscosity superdifferential of the value function of (P bio ) and p x (•), p s (•) are the adjoint states of the PMP.

This Lemma implies that whenever the value function is differentiable at ξ = (t 0 , x 0 , s 0 ), we should have that h

= -∂ t 0 V (ξ).
If ψ h is actually an optimal control, we can write the value function as the cost of the control ψ h , that is, V (ξ) = J(ξ, ψ h ) and therefore

h = -∂ t 0 J(ξ, ψ h ), (1.13)
In other words, h is a fixed point of the mapping η → -∂ t 0 J(ξ, ψ η ) and this is the key point we use to construct an algorithm to identify h.

In addition, equation (1.13) could be seen as a certificate of optimality for the feedback (1.12) if we could prove that it is an extremal control. Indeed, in this case, if the algorithm converges to a fixed point, then the computed feedback is a good approximation of an optimal control, because the reward associated with the feedback ψ h is an approximated solution to the HJB equation.

To solve (1.13), we consider the classical iterative scheme for finding a fixed point of a function by repeatedly computing the image of the previous iterate. A particularity here is that the function is composed of a partial derivative and therefore, to approximate it with a finite difference, we work with a range of initial times and use the fact that J(t 0 = T, x 0 , s 0 , ψ h ) = 0 to start, running through the initial times backwards until reaching the desired starting time. Testing this algorithm with parameters from published works, we found that it performs well and converges in only a few iterations, with the relative error decreasing rapidly. This method is a promising illustration of an interesting combination of analytical and numerical approaches. First, it allows to gain an understanding of the problem by obtaining a sufficiently explicit expression of the solution. Furthermore, the computed solutions show the complexity of the problem, with trajectories from different initial data having different singular arcs. This could be used to further justify the use of the simpler and easier to implement sub-optimal controllers (1.4) previously studied, as these have a similar form to (1.12) and they can be seen as a way of approximating the singular arcs. On the other hand, this work shows that numerical computations could be used to confirm the validity of a candidate, in addition to giving a practical way of computing a precise solution for a particular case.

Spatially Heterogeneous Bioreactor Modelling

The most common type of bioreactor in use today for anaerobic digestion is the continuously stirred tank reactor (CSTR), which uses a mixing system to homogenize reactor contents. The advantages of such a design is that it avoids the accumulation of toxic substances in one area of the tank and, in addition, mixing can prevent dead zones, with little substrate and thus a low reaction rate. However, for homogeneous reactors, the output is equal to the inner concentration, so that, to have a high level of organic matter removal, the bioreactor must be maintained at a low concentration. This means that CSTR type devices need to be operated with low dilution rates, either with a low feeding rate or by designing large reactors. Furthermore, a major disadvantage of the CSTR is that mixing represents a substantial portion of the energy required to run an anaerobic digester and can thus offset the benefits of producing biogas.

An alternative design is a un-mixed reactor that has a gradient of concentration from input to output, such as the plug flow or tubular reactor. This would potentially allow a high reaction rate at the beginning, and thus high biogas production, but with a low output concentration. However, few full scale bioreactors have been implemented for AD, as a better understanding of the impact of heterogeneity of the reacting medium on performance and stability is still required.

A variety of models have been developed to study questions related to mixing efficiency and heterogeneity. On one hand, several studies have worked with very complex models, considering multi-phase and turbulent flow, and generally using comercial computer fluid dynamics (CFD) software to run simultions [START_REF] Leonzio | Study of mixing systems and geometric configurations for anaerobic digesters using cfd analysis[END_REF][START_REF] Terashima | Cfd simulation of mixing in anaerobic digesters[END_REF][START_REF] Vesvikar | Flow pattern visualization in a mimic anaerobic digester using cfd[END_REF][START_REF] Wu | Cfd simulation of mixing in egg-shaped anaerobic digesters[END_REF]. However, the complexity of these models means that they are computationally intensive to simulate and therefore cannot be used to study optimization problems.

At the other end of the complexity spectrum, there is compartment models, which represent a reactor as a network of interconnected well-mixed zones and thus use systems of ordinary differential equations. The advantage of reducing complexity is that it allows a more in depth analysis and for example, studies with this type of models have found that the impact of heterogeneity depends on the graph of interconnections [START_REF] Bello-Mendoza | Modelling the effects of imperfect mixing on the performance of anaerobic reactors for sewage sludge treatment[END_REF][START_REF] Bensmann | Reactor configurations for biogas plants-a model based analysis[END_REF][START_REF] Donoso-Bravo | Modelling of an anaerobic plug-flow reactor. process analysis and evaluation approaches with non-ideal mixing considerations[END_REF][START_REF] Haidar | Effects of spatial structure and diffusion on the performances of the chemostat[END_REF][START_REF] Rapaport | Some non-intuitive properties of simple extensions of the chemostat model[END_REF].

A third type of model attempts to find a compromise between model complexity and physical accuracy. They represent a reactor in 1 or 2 spatial dimensions, often assuming that fluid velocity is constant in space to focus on the biological activity [START_REF] Dochain | Dynamical modelling & estimation in wastewater treatment processes[END_REF][START_REF] Kung | The growth of pure and simple microbial competitors in a moving distributed medium[END_REF][START_REF] Mu | Anaerobic digestion model no. 1-based distributed parameter model of an anaerobic reactor: I. model development[END_REF][START_REF] Winkin | Dynamical analysis of distributed parameter tubular reactors[END_REF]. Recently, 2D models coupling bio-reactions with fluid dynamics have been developed and this has allowed the consideration of optimization problems [START_REF] Barbier | Modelling of biological decontamination of a water resource in natural environment and related feedback strategies[END_REF][START_REF] Crespo | Modeling and optimization of activated sludge bioreactors for wastewater treatment taking into account spatial inhomogeneities[END_REF][START_REF]Shape optimization of spatial chemostat models[END_REF][START_REF] Moya | Mathematical modeling and optimization of bioreactors and liquid crystals[END_REF].

We introduce here a similar type of reduced complexity model for a real pilot scale bioreactor. This device is particularly adapted to develop a spatially heterogeneous model as experimental data was gathered by collecting substrate at different points along the main axis of the reactor and a spatial gradient was observed. A compartment model has already been developed for this bioreactor [START_REF] Hmissi | Evaluation of alkalinity spatial distribution in an up-flow fixed bed anaerobic digester[END_REF], considering two interconnected homogeneous zones and the two reaction model of [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF]. The experimental data was used to fit the biological parameters and this model was able to reproduce the spatial gradient roughly. The present work aims at improving the modelling of spatial features, by representing more faithfully the details of the reactors geometry. This device can be divided in three sections: the liquid enters at the bottom, goes up through a fixed bed, which helps to fix the biomass, and the output of the reactor is at the top (Figure 1.7). The contents of the reactor can be mixed by recirculating liquid from the output back into the input or by an auxiliary system that pumps liquid from the very bottom back into the tank at the same height. This makes this bioreactor particularly interesting for the study of heterogeneity as it can either be operated as a CSTR or as a tubular reactor by changing the recirculation flow rate.

The configuration of the input and output flows, means that we need to consider a 3D model for the lower and upper parts of the reactor. However, as the fixed bed is made of narrow PVC tubes, the middle section can be modelled as an array of parallel 1D tubular reactors. This considerably reduces the complexity of the model, and in particular of the fluid dynamics in this section. For the bottom and top sections, we will model the liquid as an incompressible viscous fluid and since the timescale of the dynamics of fluids is much shorter than the biological timescales, we consider the steady state Navier-Stokes equations, denoting U = (u x , u y , u z ) the fluid velocity, p the pressure and g the acceleration due to gravity,

U • ∇U -ν∆U + ∇p = g, ∇ • U = 0.
For the biological activity, we use the two reaction model of [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF] and therefore need to compute the spatial distribution of the concentrations of two substrates and two biomasses. However, to further reduce the complexity of the model, we will consider that the bio-reactions take place only in the middle section since it has been observed that most of the active biomass is attached to the fixed-bed. Then, in the bottom and top sections, the substrate concentrations S k , k = 1, 2, satisfy advection-diffusion equations

∂ t S k + U • ∇S k -D k ∆S k = 0
where D k are the diffusion coefficients. For the middle section, in each tube, we consider advection-diffusion-reaction equations,

∂ t S k + u z ∂ z S k -D k ∂ 2 zz S k = f k ,
where f k is the reaction term which depends on the concentrations of substrate and biomass through the growth rate functions.

As we suppose that the micro-organisms are fixed, the biomass concentrations B k , k = 1, 2 are neither transported nor diffused and only react. However, if we use a reaction term only based on a classical growth rate function, which depends linearly on the biomass concentration, then it will result in unbounded exponential growth of the micro-organisms.

To remedy this, we will consider a death rate τ k and to take into account crowding effects, we add a density dependent term to the growth function, taken from [START_REF] Lobry | A new hypothesis to explain the coexistence of n species in the presence of a single resource[END_REF],

g k (B k ) = 1 1 + c k √ B k .
Then, the biomasses satisfy a distributed ODE, denoting µ k the specific growth rates,

∂ t B k = µ k (S k )g k (B k )B k -τ k B k ,
The model equations can be straightforwardly solved with the Finite Element Method, but the difference of physical and biological timescales make this type of model difficult to simulate efficiently. Indeed, we are primarily interested in observing the biological activity of the reactor over the course of several days, but to correctly resolve the physical processes, we need to use a time step of the order of seconds for the bottom section due to the high fluid velocities caused by the mixing system. In addition, with small diffusion coefficients, fine meshes are required and using the reported values for D k results in long computations. However, the solution retains the same macroscopic behaviour even if we take a value for D k two orders of magnitudes greater and thus we can use a coarser mesh.

The detailed representation of the processes of advection and diffusion makes this model a valuable tool to study the impact of physical operational parameters, such as the recirculation and mixing flowrates. We first explore this question with numerical simulations of tracer experiments, which consist in releasing a short pulse of an inert substance in the inflow and measuring the concentration at the output (Figure 1.8). These reveal that the influence of the bottom mixing system is mainly seen during the transients and its efficiency is governed by the ratio between the recirculation and mixing flowrates.

Concerning the biological activity, the computationally cost of simulations make estimating the model parameters accurately unpractical. It is important to point out that we are not trying to construct a predictive model here, but instead we want to be able to reproduce the qualitative behaviour of the bioreactor. Then we can use the growth function parameters found for the compartment model of [START_REF] Hmissi | Evaluation of alkalinity spatial distribution in an up-flow fixed bed anaerobic digester[END_REF] and choose the extra parameters (death rates τ k and the density dependence parameter c k ) so that model reproduces experimental data. The results are satisfactory, considering that measurements of bioprocesses have high error margins. In particular, the ODE model has a tendency to under estimate the spatial variation of substrate concentration whereas we can reproduce it more accurately (Figure 1.9). However, the transient behaviour is poorly captured but this can be attributed to the use of a simple two reaction model and therefore, we will mainly focus on reactor performance at steady state.

To study the impact of heterogeneity on the biological activity, we run simulations for a range of recirculation flow rates. We observe that the reactor operates more efficiently with lower recirculation, removing more organic matter and producing more biogas. For small recirculation flow rates, the dilution rate is effectively lower, so that the substrate concentration entering the reactor is much higher than standard operating conditions. In this case, our simulations showed a very strong gradient of substrate concentrations, with high levels of biomass at the very beginning of the fixed bed. This indicates that the inhibition phenomenons are currently not well captured with this model as the biological parameters where obtained with data that only had low concentrations of substrate. These preliminary results of ongoing work must therefore be taken with caution.

Chapter 2 Optimal and Sub-Optimal Feedback Controls for Biogas Production

This chapter corresponds to the published article Haddon, A., Ramírez, H., and Rapaport, A. Optimal and Sub-optimal Feedback Controls for Biogas Production.

J Optim Theory Appl (2019) 183:642. https://doi.org/10.1007/s10957-019-01570-3

Introduction

Anaerobic digestion is a biological process in which organic matter is transformed by microbial species into biogas (methane and carbon dioxide). Such transformations have been used for a long time in waste water-treatment plants to purify water [START_REF] Russell | Practical wastewater treatment[END_REF]. Valorizing biogas production while treating wastewater has received recently great attention, as a way of producing valuable energy and limiting the carbon footprint of the process [START_REF] Rehl | CO2 abatement costs of greenhouse gas (GHG) mitigation by different biogas conversion pathways[END_REF]. As a final product of the biological reaction, the total production of biogas measures the performances of the biological transformation. Therefore, there is a strong interest in determining control strategies maximizing biogas production.

With continuous-stirred bioreactors, two kinds of anaerobic models are usually considered for control purposes in the literature: the one-step model, which corresponds to the classical chemostat model [START_REF] Harmand | The Chemostat: Mathematical Theory of Microorganisms Cultures[END_REF], and the two-step model that has been proposed by Bernard et al. [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF]. Although these models only have few dynamic variables, it has been shown that they are capable of reproducing the qualitative behavior of the anaerobic digestion process [START_REF] Bernard | Can we assess the model complexity for a bioprocess: theory and example of the anaerobic digestion process[END_REF]. Furthermore, in the two-step model, the second reaction is the most limiting due to inhibition by the substrate and we can then consider that a one-step model can be used to focus on the second reaction. In particular, a common assumption is to consider that the first step is fast and then the two reactions can be reduced to a single one with a slow-fast approximation and in this case, the one-step model provides a good representation of the biogas production.

The control variable is typically the input flow rate (or equivalently the dilution rate, since the volume of the reactor is constant in continuous operating mode). Several works have already considered the static optimization problem of maximizing the output flow rate of biogas at steady state, and various control strategies have been proposed to stabilize the processes at these nominal states (see for instance [START_REF] Steyer | Advanced control of anaerobic digestion processes through disturbances monitoring[END_REF][START_REF] Rodrıguez | A hydrogenbased variable-gain controller for anaerobic digestion processes[END_REF][START_REF] Dimitrova | Nonlinear stabilizing control of an uncertain bioprocess[END_REF][START_REF]Nonlinear adaptive control of a model of an uncertain fermentation process[END_REF][START_REF]Nonlinear adaptative stabilizing control of an anaerobic digestion model with unknown kinetics[END_REF][START_REF] Sbarciog | A biogas-based switching control policy for anaerobic digestion systems[END_REF][START_REF] Vargas | On-line maximization of biogas production in an anaerobic reactor using a pseudo-super-twisting controller[END_REF]).

There has been comparatively much less work considering the dynamic optimization problem over the transients, while bio-processes are often not initialized at their optimal nominal state. Although the optimal control problem, which consists in maximizing biogas production over a given time interval, has been posed a long time ago [START_REF] Stamatelatou | Optimal and suboptimal control of anaerobic digesters[END_REF], it is still unsolved today (even for the one-step model). Let us mention two attempts to solve approximately or partially this problem. Sbarciog et al. [START_REF] Sbarciog | An optimizing start-up strategy for a bio-methanator[END_REF] have considered the two-step anaerobic model and proposed a strategy for maximizing biogas production as an optimal control to drive the system in finite time in a neighborhood of the optimal steady state, with additive penalty terms in the criterion. In [START_REF] Ghouali | Maximizing biogas production from the anaerobic digestion[END_REF], Ghouali et al. give a complete solution of the original optimal control problem for the one-step model, but for a particular subset of initial conditions which belong to an invariant manifold of the system (see also [START_REF]Optimal Control in Bioprocesses: Pontryagin's Maximum Principle in Practice[END_REF]). The dynamics can be then reduced to a scalar one and the authors show that the optimal solution exhibits a singular arc with a "most rapid approach path" optimal strategy. Let us underline that optimal control problems over a fixed time horizon possess generally a time-dependent optimal synthesis, while the duration of process operation is often poorly known. However, the scalar reduced problem exhibits the remarkable feature of having an optimal synthesis independent of the terminal time, which makes it quite attractive from an application view point.

The purpose of the present article is to propose new control strategies for the one-step model, as time-independent feedbacks for general initial conditions

• either considering an infinite horizon,

• either considering sub-optimal controllers for the finite horizon.

For the infinite horizon (see for instance the book [START_REF] Carlson | Infinite horizon optimal control: deterministic and stochastic systems[END_REF]), we consider the limit of the discounted criterion (when the discount factor tends to zero) and the average cost. We study optimal strategies and compare their related optimal costs. This study extends the preliminary results presented in the conference paper [START_REF] Haddon | First results of optimal control of average biogas production for the chemostat over an infinite horizon[END_REF] and considers a large class of growth functions, that can be in particular density-dependent (such as the Contois law) or not (such as the Monod or Haldane law). Our work for the finite horizon exploits and extends an approximation technique presented in [START_REF] Haddon | Guaranteed value strategy for the optimal control of biogas production in continuous bio-reactors[END_REF]. This consists, for a given initial condition, in framing the optimal solution by considering a different reward for which the optimal solution can be determined exactly and that possess the property of having a time-independent optimal synthesis (i.e. whatever is the time horizon, finite or infinite). This technique has moreover the advantage of providing bounds on the sub-optimality of the controllers. The results are again obtained for a large class of growth functions and we show that density dependent growth functions lead to more sophisticated feedback laws.

The paper is organized as follows. Section 2.2 specifies dynamics, control, criterion and hypotheses, and gives some preliminary results about controllability and asymptotic behavior of solutions. Sections 2.3 and 2.4 study the optimal solutions, respectively for the infinite and finite time horizons. Finally, Section 2.5 illustrates our results on various growth functions.

Preliminaries

In this work, we consider the classical chemostat model [START_REF] Harmand | The Chemostat: Mathematical Theory of Microorganisms Cultures[END_REF]. This represents a well-mixed continuously fed bioreactor in which a substrate of concentration s is treated (and then transformed into biogas) by a population of microorganisms of concentration x

ṡ = u(s in -s) - 1 Y µ(s, x)x, (2.1)
ẋ = µ(s, x)x -ux. (2.2)
We denote s in > 0 the inflow concentration of substrate, Y the yield coefficient, µ(•, •) the specific growth rate and u the dilution rate, which is the control.

The biogas flowrate is assumed proportional to the growth rate so that the biogas produced during a time interval [t 0 , T ] is proportional to

T t 0 µ(s(t), x(t))x(t) dt
and, without loss of generality, we will suppose that the proportionality coefficient as well as the yield coefficient are equal to 1.

We will consider the following class of growth functions : Assumption 2.1 We suppose that µ : R + × R + → R + is a Lipschitz continuous function that satisfies, for all x > 0 µ(0, x) = 0 and µ(s, x) > 0 for s > 0.

We suppose as well that x → µ(s, x) is non increasing, which models crowding effects, and x → µ(s, x)x is non decreasing, which models the fact that having more biomass provides at least the same growth.

A typical instance of this class is the Contois growth function, defined later in (2.37), but note that this class of functions also contains growth functions that depend only on the substrate concentration, such as the Monod (2.35) and the Haldane (2.36) functions.

We will study the problem of maximizing the accumulated biogas for controls in the following set of admissible controls

U([t 0 , T ]) = u(•) ∈ L ∞ (t 0 , T ; R) : u(t) ∈ [0, u max ] for t ∈ [t 0 , T ]
with t 0 ∈ R and T ∈ R ∪ {+∞}, and where u max > 0 is a given parameter that represents the maximal dilution rate. We will consider initial conditions taken in the invariant set

D := [0, s in [×]0, ∞[
which corresponds to the most common operating conditions. Notice that for initial conditions in D, any solution of (2.1)-(2.2) cannot reach s = s in in finite time and stays non negative. Therefore the set D is (forward) invariant.

Properties of the Dynamics

On the invariant domain D, we introduce the change of variables

ζ = (s, z) with z = x s in -s ,
under which the dynamics become

ζ = ṡ ż = f (ζ, u) := u -µ s, (s in -s)z z (s in -s) µ s, (s in -s)z (1 -z)z . (2.3)
We will denote s t 0 ,ξ,u (•) and z t 0 ,ξ,u (•) the solution of (2.3), with initial condition ξ = (s 0 , z 0 ) = (s(t 0 ), z(t 0 )) ∈ D and control u(•) ∈ U([t 0 , T ]). The cumulated biogas production becomes

T t 0 φ s t 0 ,ξ,u (t), z t 0 ,ξ,u (t) z t 0 ,ξ,u (t) dt (2.4) with φ(s, z) = µ s, (s in -s)z (s in -s) (2.5)
and we will denote

φ(z) = max s∈]0,s in [ φ(s, z). (2.6)
We can now establish an important property of the controlled dynamics.

Lemma 2.2

The trajectories of the system (2.3) for a given initial condition ξ = (s 0 , z 0 ) ∈ D, for all admissible controls, remain in the set

L(ξ) = [0, s in ] × [min(z 0 , 1), max(z 0 , 1)].
(2.7)

Proof. From Assumption 2.1 we have that µ(•, •) 0 and since the solutions z(•) satisfy (2.3), we then have the following min(z 0 , 1) z t 0 ,ξ,u (t) max(z 0 , 1) for all t 0, for any admissible control u(•).

In the following, we consider initial conditions that guarantee the controllability of the s variable.

Assumption 2.3 We suppose that the initial condition ξ ∈ D is such that max (s,z)∈L(ξ) µ s, (s in -s)z z < u max .
In practice, for a given initial condition it possible to choose u max such that the previous inequality is satisfied.

We now define a class of feedbacks, that will play an important role, and that are based on the notion of most rapid approach path, a well known concept in the theory of optimal control; see, for example, [START_REF] Rapaport | Competition between most rapid approach paths: necessary and sufficient conditions[END_REF][START_REF] Hartl | A new sufficient condition for most rapid approach paths[END_REF].

Definition 2.4 For (s, z) ∈ L(ξ), we define the most rapid approach feedback to a given substrate level s * ∈ [0, s in [, as

ψ s * (s, z) = 0 if s > s * , µ(s * , (s in -s * )z) z if s = s * , u max if s < s * .
(2.8)

Clearly, with Assumption 2.3 this feedback is well defined, so that, associated with this control, for every initial condition ξ ∈ D, there exists a unique absolutely continuous solution for the dynamics (2.3).

Lemma 2.5 For any ξ ∈ D satisfying Assumption 2.3, a given substrate level s * ∈]0, s in [ is reachable in finite time with the feedback ψ s * .

Proof. First, using the monotonicity properties of µ(•, •) of Assumption 2.1, it is clear that ψ s * is admissible provided Assumption 2.3 is satisfied.

To show that s * is reachable in finite time, it is enough to note that when s t 0 ,ξ,ψ s * (t) > s * , for t in a given open interval I, we have

ṡt 0 ,ξ,ψ s * (t) = -µ s, (s in -s)z z(s in -s) k -< 0, ∀ t ∈ I with k -= -min s∈]s * ,s in [ µ s, (s in -s) min(z 0 , 1) min(z 0 , 1)(s in -s * ). This insures that s * is always reachable in finite time from s 0 > s * . Analogously, if s t 0 ,ξ,ψ s * (t) < s * , for t ∈ I, we have from Assumption 2.3 ṡt 0 ,ξ,ψ s * (t) = u max -µ s, (s in -s)z z (s in -s) k + > 0, ∀ t ∈ I with k + = u max -max s∈]0,s * [
µ s, (s in -s) max(z 0 , 1) max(z 0 , 1) (s in -s * ). Then s * is reachable from s 0 < s * , again in finite time.

Remark It should be pointed out that there is a similarity with the turnpike property [START_REF] Zaslavski | Turnpike properties in the calculus of variations and optimal control[END_REF][START_REF] Trélat | The turnpike property in finite-dimensional nonlinear optimal control[END_REF] when using the controller (2.8). The turnpike property has received great attention in the literature (see, for instance, [START_REF] Grüne | On the relation between turnpike properties for finite and infinite horizon optimal control problems[END_REF][START_REF] Rapaport | Competition between most rapid approach paths: necessary and sufficient conditions[END_REF][START_REF] Hartl | A new sufficient condition for most rapid approach paths[END_REF][START_REF]Nonturnpike optimal solutions and their approximations in infinite horizon[END_REF]), and recent results give sufficient optimality conditions [START_REF] Faulwasser | On turnpike and disspativity properties of continuous-time optimal control problems[END_REF][START_REF] Faulwasser | Exact turnpike properties and economic nmpc[END_REF]. However, we shall show in the next sections that the value s * , which determines the turnpike, has to depend on the initial condition (excepted for the very particular case when the initial condition belongs to the invariant set {z = 1} that has been solved in [START_REF] Ghouali | Maximizing biogas production from the anaerobic digestion[END_REF]). So, we are not in the usual framework of a single turnpike [START_REF] Faulwasser | On turnpike and disspativity properties of continuous-time optimal control problems[END_REF][START_REF] Faulwasser | Exact turnpike properties and economic nmpc[END_REF] or isolated turnpikes [START_REF] Rapaport | Turnpike theorems by a value function approach[END_REF], and the results of the literature do not apply.

For the problem on an infinite horizon, we will consider persistently exciting controls, which are defined as satisfying

T t 0 u(t) dt -→ T →∞ ∞.
As the next Lemma shows, the trajectories associated with these controls are such that z t 0 ,ξ,u (t) converges to 1, which is essential in our approach. Furthermore, for non persistently exciting controls, s t 0 ,ξ,u (t) converges to 0 and thus the biogas production also converges to 0. As a consequence, the controls that maximize biogas production are necessarily persistently exciting controls. Lemma 2.6 For all initial conditions ξ ∈ D and for all persistently exciting controls u

(•) ∈ U([0, ∞[), we have lim t→∞ z 0,ξ,u (t) = 1
and

lim δ→0 ∞ 0 δe -δt z 0,ξ,u (t) dt = lim T →∞ 1 T T 0 z 0,ξ,u (t) dt = 1.
Moreover, for non persistently exciting controls, we have

lim t→+∞ s 0,ξ,u (t) = 0.
Proof. From equation (2.3), the solution z(•) = z 0,ξ,u (•) can be written as follows

z(t) = z 0 + e t t 0 µ(s(τ ),x(τ )) dτ 1 + z 0 e t t 0 µ(s(τ ),x(τ )) dτ -1 (2.9)
where

s(•) = s 0,ξ,u (•), x(•) = x 0,ξ,u (•). From equation (2.2), the solution x(•) is such that x(t) = x(t 0 )e t t 0 µ(s(τ ),x(τ ))-u(τ ) dτ .
Therefore, if the integral function

t → t t 0 µ(s(τ ), x(τ )) dτ, t ≥ t 0 (2.10)
is bounded, then x(t) must converge asymptotically to 0 when t goes to +∞ and u(•) is a persistently exciting control. Moreover, from equations (2.1), (2.2) we have

d dt s(t) + x(t) = u(t) s in -s(t) + x(t) so that s(t) + x(t) = s in + (s(t 0 ) + x(t 0 ) -s in e -t t 0 u(τ ) dτ
and then s(t) must converge to s in when t goes to +∞. Consequently, by continuity of the function µ, there exists T > t 0 such that µ(s(t), x(t))) > µ(s in , 0)/2 > 0 for any t > T , which implies that the integral defined in (2.10) goes to +∞ when t goes to +∞, which is a contradiction. We deduce that this integral cannot be bounded and from equation (2.9) that z(t) converges to 1 when t goes to +∞.

A proof of the equality of limits of the integrals

lim δ→0 ∞ 0 δe -δt z 0,ξ,u (t) dt = lim T →∞ 1 T T 0 z 0,ξ,u (t) dt
can be found in [START_REF] Grüne | Asymptotic controllability and exponential stabilization of nonlinear control systems at singular points[END_REF]Lemma 3.5]. For the value of the limits we use the fact that z 0,ξ,u (t) converges to 1 : for all ε > 0, there exits a time t ε such that, for all t t ε,

|z 0,ξ,u (t) -1| < ε.
Then, for all T max(t ε, t ε/ε)

1 T T 0 z 0,ξ,u (t) dt -1 1 T t ε 0 |z 0,ξ,u (t) -1| dt + 1 T T t ε |z 0,ξ,u (t) -1| dt < t ε T |z 0 -1| + 1 - t ε T ε < ε (|z 0 -1| + 1) .
With this, for all ε > 0, we can take ε = ε/(|z 0 -1| + 1) and then we have, for

T max(t ε, t ε/ε) 1 T T 0 z 0,ξ,u (t) dt -1 < ε.
Finally, we prove that for non persistently exciting controls, s 0,ξ,u (t) converges to 0. Therefore, suppose that u(•) is an admissible control with a finite integral and we define, for all t 0,

I(t) := t 0 u(τ ) dτ < ∞ and ϕ(t) := (s in -s 0,ξ,u (t))e I(t) .
Then ϕ (t) = φ s 0,ξ,u (t), z 0,ξ,u (t) z 0,ξ,u (t)e I(t) 0 and since ϕ(t) is bounded, we can deduce that ϕ(t) converges as t goes to infinity. Note as well that ϕ is absolutely continuous and thus uniformly continuous. We can therefore use Barbalat's Lemma [START_REF] Khalil | Nonlinear systems[END_REF]Lemma 4.2] to get that ϕ (t) converges to 0. Then, as z 0,ξ,u (t) cannot reach 0 (Lemma 2.2), we have that φ s 0,ξ,u (t), z 0,ξ,u (t) must converge to 0 and by continuity we conclude that s 0,ξ,u (t) converges to 0.

Infinite Horizon and Average Reward

In this section, we study the problem of maximizing biogas production over an infinite horizon. Since the dynamics (2.3) are autonomous, without loss of generality, we can assume here that t 0 = 0 and we will then denote s ξ,u (•) and z ξ,u (•) solutions of (2.3).

We start by defining the average biogas production during a time interval [0, T ] as

J T (ξ, u(•)) = 1 T T 0 φ s ξ,u (t), z ξ,u (t) z ξ,u (t) dt (2.11)
and we consider the inferior and superior limits as T goes to infinity

J ∞ (ξ, u(•)) = lim inf T →∞ J T (ξ, u(•)), (2.12
)

J ∞ (ξ, u(•)) = lim sup T →∞ J T (ξ, u(•)). (2.13)
The optimal control problems in consideration here consist in maximizing these functionals with respect to the dilution rate u(•) ∈ U([0, ∞[), for any initial condition ξ ∈ D. More precisely, the value functions of these optimal control problems are

V ∞ (ξ) = sup J ∞ (ξ, u(•)) : u(•) ∈ U([0, ∞[) , (2.14) V ∞ (ξ) = sup J ∞ (ξ, u(•)) : u(•) ∈ U([0, ∞[) .
(2.15)

We need to consider the inferior and superior limits here as there exists controls for which the rewards (2.12) and (2.13) may differ. Indeed, this is the case for certain oscillating controls as can be seen in the example in the Appendix. Nevertheless, we will show that the value functions (2.14) and (2.15) are in fact equal. Moreover, we will connect these problems to the problem with a discounted reward when the discount factor goes to 0, as in [START_REF]On the relation between discounted and average optimal value functions[END_REF].

To this end, we now define the following discounted reward, for a discount rate δ > 0

J δ (ξ, u(•)) = δ ∞ 0 e -δt φ s ξ,u (t), z ξ,u (t) z ξ,u (t) dt.
(2.16)

This type of cost function is often used in problems related to economics for which the term e -δt represents a discount rate or a preference for the present [START_REF] Carlson | Infinite horizon optimal control: deterministic and stochastic systems[END_REF]. In our setting, the use of this discounted reward can be seen as a preference for earlier rather than later production. Here, the integral is rescaled with the discount factor δ in order to guarantee that, when we take the limit as δ goes to 0, the reward remains finite.

The value function of the optimal control problem for a given δ is then

V δ (ξ) = sup J δ (ξ, u(•)) : u(•) ∈ U([0, ∞[) .
(2.17)

Note that both average rewards (2.12) and (2.13), as well as the discounted reward (2.16), are well defined as the following Lemma shows.

Lemma 2.7 For all ξ ∈ D, for all admissible controls u(•) ∈ U([0, ∞[) and for all δ > 0, the rewards J ∞ (ξ, u(•)), J ∞ (ξ, u(•)) and J δ (ξ, u(•)) are uniformly bounded.

Proof. From the monotonicity properties of Assumption 2.1, we have that the function z → φ(s, z) is non increasing. for all s > 0. Thus, for all t 0

φ(s ξ,u (t), z ξ,u (t)) φ(0).
The uniform boundedness of the rewards then follows from Lemma 2.2.

Solution of Optimal Control Problems for the Average Rewards

We now solve the optimal control problems (2.14) and (2.15). We start by determining an upper bound for the value functions and then we will exhibit controls that attain this bound.

Proposition 2.8 For all initial conditions ξ ∈ D

V ∞ (ξ) V ∞ (ξ) max s∈]0,s in [ φ(s, 1).
Proof. With the monotonicity properties of µ(•, •) of Assumption 2.1, we have that z → φ(s, z) is non increasing and z → φ(s, z)z is non decreasing. This implies that φ(s, max(z 0 , 1)) φ(s, z) φ(s, min(z 0 , 1)) (2.18) and φ(s, min(z 0 , 1)) min(z 0 , 1) φ(s, z)z φ(s, max(z 0 , 1)) max(z 0 , 1).

(2.19)

First, we consider the case when z 0 1. For any control u(•), we have

J T (ξ, u(•)) 1 T T 0 φ s(t), max(z 0 , 1) max(z 0 , 1) dt max s∈]0,s in [ φ(s, 1) = φ(1).
Taking the lower and upper limit as T goes to infinity and the supremum with respect to u(•) we get the result.

Next, for z 0 1, we have

J T (ξ, u(•)) 1 T T 0 φ s(t), min(z 0 , 1) z(t) dt max s∈]0,s in [ φ(s, 1) 1 T T 0 z(t) dt.
Using Lemma 2.6, we get that J ∞ (ξ, u(•)) J ∞ (ξ, u(•)) φ(1) and we conclude taking the supremum with respect to u(•).

Note that the existence of a maximum of s → φ(s, 1) = µ(s, s in -s)(s in -s) on ]0, s in [ follows from Assumption 2.1. We will denote a substrate level at which such a maximum is attained as s = arg max

s∈]0,s in [ φ(s, 1)
Proposition 2.9 For any initial condition ξ ∈ D, any control u(•) ∈ U([0, ∞[) that drives the system asymptotically to the state (s, 1) is optimal for problems (2.14) and (2.15). We then have

V ∞ (ξ) = V ∞ (ξ) = φ(s, 1) = φ(1). ( 2

.20)

Proof. The continuity of φ implies that for all ε > 0, there exists a time t ε 0 such that, for all t t ε , φ s ξ,ū (t), z ξ,ū (t) z ξ,ū (t) -φ s, 1 < ε.

(2.21)

Since s ξ,ū (•) and z ξ,ū (•) take values in the compact set L(ξ) (2.7), there is a constant M ξ > 0 such that, for all t 0,

φ s ξ,ū (t), z ξ,ū (t) z ξ,ū (t) < M ξ . (2.22)
Then, for all T t ε , from (2.21) and (2.22)

J T (ξ, u(•)) -φ s, 1 1 T tε 0 φ s ξ,ū (t), z ξ,ū (t) z ξ,ū (t) -φ s, 1 dt + 1 T T tε φ s ξ,ū (t), z ξ,ū (t) z ξ,ū (t) -φ s, 1 dt < 2M ξ t ε T + 1 - t ε T ε and we have J ∞ (ξ, u(•)) = J ∞ (ξ, u(•)) = φ(s, 1).
Using Proposition 2.8, we get the equality of value functions (2.20) and deduce the optimality of u(•) for both average biogas production problems (2.14) and (2.15).

With Lemma 2.6, we know that all persistently exciting admissible controls make z(•) converge to 1, and from Lemma 2.5, we know that the feedback ψ s * defined in (2.8) with s * = s guarantees that s(•) reaches s. Then, from the previous Proposition we have the following result.

Proposition 2.10 For any initial condition ξ ∈ D satisfying Assumption 2.3, the most rapid approach feedback to s, defined in (2.8) and denoted ψ s, is optimal for both average production problems (2.14) and (2.15). Clearly, there is not a unique optimal control for the infinite horizon problems that we have considered. For example, in the case of a growth function that depends only on the substrate and that is monotone (such as the Monod growth function), the constant control u = µ(s) can also drive the system to the state (s, 1). Nonetheless, for the control ψ s, we are able to state an estimation of the sub-optimality for the finite horizon problem.

Relation between Average and Discounted Biogas Production Problems

We now discuss the relation between the average and discounted biogas production problems. We first show that the value function of the discounted problem converges when the discount factor δ goes to 0.

Proposition 2.11 For all ξ ∈ D, the value function of the discounted problem (2.17) converge as δ goes to 0 to the value functions of the average problems (2.14) and (2.15),

V 0 (ξ) := lim δ→0 V δ (ξ) = V ∞ (ξ) = V ∞ (ξ). ( 2

.23)

Proof. This is a consequence of [START_REF]On the relation between discounted and average optimal value functions[END_REF] on the relation between average and discounted functionals. First, [START_REF]On the relation between discounted and average optimal value functions[END_REF]Lemma 3.3] gives

sup u(•) lim inf T →∞ J T (ξ, u(•)) = lim T →∞ sup u(•) inf τ T J τ (ξ, u(•)) and sup u(•) lim sup T →∞ J T (ξ, u(•)) = lim T →∞ sup u(•) sup τ T J τ (ξ, u(•))
Now, denote M ξ the upper bound of φ(•, •) on the compact set L(ξ) (2.7), then [START_REF]On the relation between discounted and average optimal value functions[END_REF]Lemma 3.4] states that for all T > 0, ε > 0 and δ < ε/2T M ξ and all admissible controls, inf

τ T J τ (ξ, u(•)) -ε J δ (ξ, u(•)) sup τ T J τ (ξ, u(•)) + ε. (2.24)
Therefore, taking the supremum over all admissible controls and the limit as T goes to infinity, the result is obtained thanks to the equality of the average value functions of Proposition 2.9.

Now the question is whether the limit (2.23) is the value function corresponding to the problem for limit of the discounted reward when δ goes to 0. For this we would first need to show the convergence of the discounted rewards but this remains an open question for a general control. Notice nonetheless, that [START_REF]On the relation between discounted and average optimal value functions[END_REF]Lemma 3.4] actually gives us the following estimation,

J ∞ (ξ, u(•)) lim inf δ→0 J δ (ξ, u(•)) lim sup δ→0 J δ (ξ, u(•)) J ∞ (ξ, u(•))
Then, for a given control, if the average cost converges when T goes to infinity, i.e.

J ∞ (ξ, u(•)) = J ∞ (ξ, u(•))
, then the discounted reward also converges. In particular, this means that for the MRAP feedback to s, defined in (2.8), the discounted reward converges to the limit of the value function,

lim δ→0 J δ (ξ, ψ s) = lim δ→0 V δ (ξ). (2.25)
Another interesting question is to determine if the optimal controls converge when δ goes to 0. For this, we could show the Γ-convergence of the rewards and the existence of optimal controls for δ > 0. We specify this last point in Proposition 2.14 and for this, we will consider the discounted reward (2.16) as a function of the trajectory ζ(•) = s ξ,u (•), z ξ,u (•) instead of the control and with a slight abuse of notation, we will denote it as J δ (ζ(•)).

Define the set valued map

F (ζ) := u∈[0,umax] f (ζ, u)
and consider the set of all forward trajectories of (2.3) with initial condition ξ,

S(ξ) := ζ(•) ∈ AC([0, ∞[, L(ξ)) : ζ(0) = ξ, ξ(t) ∈ F (ξ(t)) a.e. t ∈ [0, ∞[ ,
where AC([0, ∞[, L(ξ)) denotes the set of absolutely continuous functions from [0, ∞[ to L(ξ). We recall from the Filippov Selection Theorem (see for instance [START_REF] Vinter | Optimal Control[END_REF]) that the optimal control problem (2.17) is equivalent to the optimization problem on S(ξ),

V δ (ξ) = sup J δ (ζ(•)) : ζ(•) ∈ S(ξ) .
We now specify the topology that we will use. Definition 2.12 For b > 0, we denote by L 1 0, ∞; R 2 , e -bt dt the weighted Lebesgue space of measurable functions y(

•) from [0, ∞[ to R 2 such that ∞ 0 ||y(t)||e -bt dt < ∞
and we denote W 1,1 0, ∞; R 2 , e -bt dt the weighted Sobolev space of measurable functions y(•) satisfying

y(•) ∈ L 1 0, ∞; R 2 , e -bt dt and ẏ(•) ∈ L 1 0, ∞; R 2 , e -bt dt .
We consider the topology on W 1,1 0, ∞; R 2 , e -bt dt for which a sequence y n (•) converges to y(•), if and only if,

-y n (•) converges uniformly to y(•) on compact intervals, -ẏn (•) converges weakly to ẏ(•) in L 1 0, ∞; R 2 , e -bt dt .
Next, we define the notion of Γ-limit in our context (see [START_REF] Maso | An introduction to Γ-convergence[END_REF] for further details). Definition 2.13 For a given initial condition ξ ∈ D and trajectory ζ(•) ∈ S(ξ), the Γ-lower limit and Γ-upper limit of J δ (•) are

Γ -lim inf δ→0 J δ (ζ(•)) = sup V∈N (ζ(•)) lim inf δ→0 inf η(•)∈V J δ (η(•)) Γ -lim sup δ→0 J δ (ζ(•)) = sup V∈N (ζ(•)) lim sup δ→0 inf η(•)∈V J δ (η(•)).
Here, N (ζ(•)) denotes the set of all open neighbourhoods of ζ(•), for the topology of the space W 1,1 0, ∞; R 2 , e -bt dt given in Definition 2.12. If both of these limits coincide, then the Γ-limit of J δ (•) is

Γ -lim δ→0 J δ (ζ(•)) = Γ -lim inf δ→0 J δ (ζ(•)) = Γ -lim sup δ→0 J δ (ζ(•)).
We now have the following result.

Proposition 2.14 For all ξ ∈ D and for all δ > 0, the suprema are attained,

V δ (ξ) = max ζ(•) J δ (ζ(•)).
If the Γ-limit of J δ (•) exists as δ goes to 0,

J 0 (ζ(•)) := Γ -lim δ→0 J δ (ζ(•)),
then the maxima converge, pointwise in ξ, to the maximum of the limit,

V 0 (ξ) := lim δ→0 V δ (ξ) = max ζ(•) J 0 (ζ(•)).
(2.26)

Furthermore, if ζ δ (•) is an optimal trajectory for (2.17), i.e. if V δ (ξ) = J δ (ζ δ (•)), and if ζ δ (•) converges to ζ 0 (•) in S(ξ), then ζ 0 (•)
is an optimal trajectory for (2.26) and

V 0 (ξ) = J 0 (ζ 0 (•)) = lim δ→0 J δ (ζ δ (•)).
Proof. To show that the suprema are attained we show that the the set of all forward trajectories of (2.3) with initial condition ξ is compact for the topology on W 1,1 0, ∞; R 2 , e -bt dt given in Definition 2.12.

For each ξ ∈ D we set

F ξ (ζ) := F P L(ξ) (ζ)
where P L(ξ) is the projection on the convex set L(ξ). Then F ξ has linear growth, so that we can define

c = sup ζ∈Dom(F ξ ) ||F ξ (ζ)|| ||ζ|| + 1 ,
where

||F ξ (ζ)|| := sup η∈F ξ (ζ) ||η||.
Note that F is upper semi-continuous and has compact non-empty convex images (such a map is known as a Marchaud map [START_REF] Aubin | Viability theory[END_REF]). With this, the set S(ξ) is the set of absolutely continuous solutions of the differential inclusion

ζ(t) ∈ F ξ (ζ(t)), ζ(0) = ξ.
We can therefore use [1, Theorem 3.5.2] to establish that S(ξ) is compact for the topology of W 1,1 0, ∞; R 2 , e -bt dt for b > c, thereby proving the existence of optimal trajectories in S(ξ).

In addition, this allows us to show that the maxima converge to to the maximum of the limit. Indeed, when the rewards Γ-converge, it is sufficient to show that there exists a countably compact set on which the suprema are attained for all δ [25, Theorem 7.4]. The set S(ξ) is clearly independent of δ and countably compact, since it is compact. Finally, the convergence of optimal trajectories can be shown with [START_REF] Maso | An introduction to Γ-convergence[END_REF]Corollary 7.20].

Remark Notice that estimation (2.24) is obtained uniformly in the controls, so that the discounted reward converges uniformly over the class of controls for which the average cost converges. Uniform convergence implies Γ-convergence ([25, Proposition 5.2]), so the consequence is that if we restrict the problem to this class of controls then the limit problem is well defined. Then (2.25) can be written as

J 0 (ξ, ψ s) = V 0 (ξ)
which proves the optimality of the MRAP feedback ψ s for the limit problem restricted to the aforementioned class of controls.

Finite Horizon and Sub-optimal Controls

We now examine the problem of maximizing biogas production over a finite horizon for a time interval [t 0 , T ] where T is fixed. For this we consider the following reward

J(t 0 , ξ, u(•)) = T t 0 φ s t 0 ,ξ,u (t), z t 0 ,ξ,u (t) z t 0 ,ξ,u (t) dt (2.27)
where we recall that s t 0 ,ξ,u (•), z t 0 ,ξ,u (•) is the solution of (2.3) with control u(•) ∈ U([t 0 , T ]) and initial condition ξ ∈ D. The optimal control problem consists in maximizing this functional with respect to the dilution rate, so that the associated value function is

V (t 0 , ξ) = sup J(t 0 , ξ, u(•)) : u(•) ∈ U([t 0 , T ]) .
(2.28)

We also consider auxiliary optimal control problems, which consist in maximizing the cost, for a given z 1 ∈ [min(z 0 , 1), max(z 0 , 1)],

J z 1 (t 0 , ξ, u(•)) = T t 0 φ(s t 0 ,ξ,u (t), z 1 ) dt (2.29)
for the same dynamics (2.3). The value functions of these auxiliary problems are then defined as

W z 1 (t 0 , ξ) = sup J z 1 (t 0 , ξ, u(•)) : u(•) ∈ U([t 0 , T ]) . (2.30)
The resolution of these auxiliary problems will be presented in Section 2.4.1.

We now show that the value functions of the original problem (2.28) and the auxiliary problems (2.30) are related.

Proposition 2.15 For all ξ ∈ D, t 0 < T and any z 1 ∈ [min(z 0 , 1), max(z 0 , 1)], we have the following frame for the value function V of the original problem

min(z 0 , 1)W z 1 (t 0 , ξ) V (t 0 , ξ) max(z 0 , 1)W z 1 (t 0 , ξ).
(2.31)

Proof. We start with the case z 0 1. For a given control u(•) ∈ U([t 0 , T ]), we define the following time

t 1 = inf {t t 0 : z t 0 ,ξ,u (t) = z 1 } ∧ T
which it is well defined since z t 0 ,ξ,u (•) is monotonous. Then, for t 0 t t 1 we have z 0 z t 0 ,ξ,u (t) z 1 1 and with the monotonicity properties of µ(•, •) of Assumption 2.1 we have φ(s t 0 ,ξ,u (t), z 1 )z 0 φ(s t 0 ,ξ,u (t), z t 0 ,ξ,u (t))z t 0 ,ξ,u (t) φ(s t 0 ,ξ,u (t), z 1 )z 1 .

Next, for t 1 t T we have z 0 z 1 z t 0 ,ξ,u (t) 1 and φ(s t 0 ,ξ,u (t), z 1 )z 1 φ(s t 0 ,ξ,u (t), z t 0 ,ξ,u (t))z t 0 ,ξ,u (t) φ(s t 0 ,ξ,u (t), z 1 ).

Combining these inequalities we get

t 1 t 0 φ(s t 0 ,ξ,u (t), z 1 )z 0 dt+ T t 1 φ(s t 0 ,ξ,u (t), z 1 )z 1 dt J(t 0 , ξ, u(•)) t 1 t 0 φ(s t 0 ,ξ,u (t), z 1 )z 1 dt + T t 1 φ(s t 0 ,ξ,u (t), z 1 ) dt. Now, since z 0 z 1 1 we have z 0 J z 1 (t 0 , ξ, u(•)) J(t 0 , ξ, u(•)) J z 1 (t 0 , ξ, u(•)).
For the case z 0 1, we proceed in a similar way to get

J z 1 (t 0 , ξ, u(•)) J(t 0 , ξ, u(•)) z 0 J z 1 (t 0 , ξ, u(•)).
We conclude by taking the supremum over all admissible controls.

The interest of the previous frames on the value functions is that it allows to find controls for which we have an estimation of sub-optimality for the original problem. Proposition 2.16 For all ξ ∈ D and all t 0 < T , any optimal control u z 1 (•) for the reward J z 1 (t 0 , ξ, •) guarantees a (sub-optimal) value for the original criterion J(t 0 , ξ, •) that satisfies min(z 0 , 1)W z 1 (t 0 , ξ) J(t 0 , ξ, u z 1 (•)) max(z 0 , 1)W z 1 (t 0 , ξ)

(2.32)

and we have the following estimation of the value function V

V (t 0 , ξ) -J(t 0 , ξ, u z 1 (•)) |1 -z 0 |W z 1 (t 0 , ξ). ( 2

.33)

Proof. From the proof of Proposition 2.15, for any control u(•) ∈ U([t 0 , T ]), we have min(z 0 , 1)J z 1 (t 0 , ξ, u(•)) J(t 0 , ξ, u(•)) max(z 0 , 1)J z 1 (t 0 , ξ, u(•)).

Evaluating this for any optimal control u z 1 (•) for the reward J z 1 (t 0 , ξ, •) gives the suboptimality frame (2.32). The sub-optimality estimation (2.33) then follows from (2.31) and (2.32).

Resolution of Auxiliary Problems

In order to obtain sub-optimal controls for problem (2.28) we now need to solve the auxiliary problem (2.30) for a given z 1 ∈ [min(z 0 , 1), max(z 0 , 1)]. The optimal control of this auxiliary problem is an autonomous feedback, even though the horizon is fixed and finite. It is similar to the optimal feedback for the infinite horizon problem ψ s, defined in (2.8), and it drives the system towards a maximizer of s → φ(s, z 1 ) but now, this maximizing substrate level depends on z 1 . We first need an assumption on the uniqueness of a maximum of φ(•, z 1 ).

Assumption 2.17 For each z 1 0, the function s → φ(s, z 1 ) admits a unique maximum on ]0, s in [, and we denote the substrate level at which this maximum is attained as

s(z 1 ) = arg max s∈]0,s in [ φ(s, z 1 ). ( 2

.34)

Note that implies that s → φ(s, z 1 ) is increasing on ]0, s(z 1 )] and decreasing on [s(z 1 ), s in [. Proposition 2.18 For all ξ ∈ D satisfying Assumption 2.3 and all t 0 < T , the most rapid approach feedback to s(z 1 ), defined in (2.8) and denoted ψ s(z 1 ) , is optimal for the auxiliary problem (2.30).

Proof. We start with the case s 0 s(z 1 ). With the control u = 0, the solution of (2.3) is such that s t 0 ,ξ,0 (•) is monotonic and non increasing. Therefore there exists a time t min , possibly larger than T , such that s t 0 ,ξ,0 (t min ) = s(z 1 ) and then the solution with the feedback (2.8) is, with t * = min(t min , T )

s t 0 ,ξ,ψ s(z 1 ) (t) = s t 0 ,ξ,0 (t), if t 0 t < t * , s(z 1 ), if t * t T.
Next, for all u ∈ [0, u max ] and for all (s, z) ∈ L(ξ),

-µ(s, (s in -s)z)(s in -s)z (s in -s)u -µ(s, (s in -s)z)(s in -s)z.
By the theorem of comparison of solutions of scalar differential equations, this implies that s t 0 ,ξ,0 (t) s t 0 ,ξ,u (t), up to time t * , for all controls u(•) ∈ U([t 0 , T ]). Since s → φ(s, z 1 ) is decreasing on [s(z 1 ), s in [, we have φ(s t 0 ,ξ,0 (t), z 1 ) φ(s t 0 ,ξ,u (t), z 1 ).

Finally, as s → φ(s, z 1 ) reaches its maximum at s(z 1 ) we get

J z 1 (t 0 , ξ, ψ s(z 1 ) ) = t * t 0 φ(s t 0 ,ξ,0 (t), z 1 )dt + T t * φ(s(z 1 ), z 1 )dt T t 0 φ(s t 0 ,ξ,u (t), z 1 )dt = J z 1 (t 0 , ξ, u).
We now consider s 0 < s. From Assumption 2.3, the feedback is admissible and we have

u max µ(s, (s in -s)z)z for all (s, z) ∈ L(ξ)
Thus, with the control u = u max , the solution of (2.3) is such that s t 0 ,ξ,umax (•) is monotone and non decreasing. Therefore, there exists a time t max , possibly larger than T , such that s t 0 ,ξ,umax (t max ) = s(z 1 ) and then the solution with the feedback (2.8) is, with t * = min(t max , T )

s t 0 ,ξ,ψ s(z 1 ) (t) = s t 0 ,ξ,umax (t), if t 0 t < t * , s(z 1 ), if t * t T.
Next, for all u ∈ [0, u max ] and for all (s, z) ∈ L(ξ)

(s in -s)(u max -µ(s, (s in -s)z)z) (s in -s)(u -µ(s, (s in -s)z)z)
and this implies that s t 0 ,ξ,umax (t) s t 0 ,ξ,u (t), up to time t * , for all controls u(•) ∈ U([t 0 , T ]). Since s → φ(s, z 1 ) is increasing on ]0, s(z 1 )], we have φ(s t 0 ,ξ,umax (t), z 1 ) φ(s t 0 ,ξ,u (t), z 1 ).

Finally, since s → φ(s, z 1 ) reaches its maximum at s(z 1 ), we get

J z 1 (t 0 , ξ, ψ s(z 1 ) ) = t * t 0 φ(s t 0 ,ξ,umax (t), z 1 )dt + T t * φ(s(z 1 ), z 1 )dt T t 0 φ(s t 0 ,ξ,u (t), z 1 )dt = J z 1 (t 0 , ξ, u).

Application to Particular Growth Functions

The controls that we have considered up to now are all most rapid approach feedbacks to s(z 1 ), with z 1 ∈ [min(z 0 , 1), max(z 0 , 1)], and this leads to the question of which is best in terms of biogas production. It turns out that it depends on the initial conditions and the horizon considered.

Indeed, we know that for an infinite horizon, the feedback ψ s(z 1 ) with z 1 = 1 is optimal and we can then expect that when the horizon is large, the best of the considered feedbacks would be for z 1 close to 1. On the other hand, when the horizon is small, the feedback ψ s(z 0 ) would seem to be the best option since this strategy consists in remaining close to the maximum of the biogas flow rate corresponding to the initial condition, whereas another feedback could drive the system away, towards another maximizing state but that can not be reached in time.

In this section, we apply our main results to the most common growth functions and explore with numerical simulations the question of determining the best feedback ψ s(z 1 ) for a given initial condition and final time. In particular, we will work with the Monod function

µ M (s) = µ max s K s + s (2.35) the Haldane function µ H (s) = μs K s + s + s 2 K i (2.36)
and the Contois function

µ C (s, x) = µ max s K s x + s (2.37)
where µ max , μ, K s and K i are positive numbers. We shall see later that these functions satisfy our assumptions (Lemma 2.19).

First, note that the Monod and Haldane functions only depend on the substrate, so that in this case, the maximizers s(z 1 ), defined in (2.34), are all equal to s(1) = s, for all z 1 ∈ [min(z 0 , 1), max(z 0 , 1)]. We illustrate the associated feedback ψ s for a Haldane function with a graph of the state space trajectories in Figure 2.1. The case of a Monod function leads to a similar dynamical behavior and the only major difference is the value of s.

From now on we will only consider the Contois growth function, for which we plot the trajectories in state space obtained with the feedback ψ s(z 0 ) in Figure 2.2.

To determine which of the feedbacks ψ s(z 1 ) is the best, we now compute the associated reward for a range of values of z 1 ∈ [min(z 0 , 1), max(z 0 , 1)] and of final times for a given initial condition. In order to easily identify the maximum of J(ξ, ψ s(z 1 ) (•)) with respect to z 1 , we normalize the average reward (2.11) by computing

J N (T, z 1 ) = J T (ξ, ψ s(z 1 ) (•)) -min y J T (ξ, ψ s(y) (•)) max y J T (ξ, ψ s(y) (•)) -min y J T (ξ, ψ s(y) (•))
where the minimum and maximum are taken for y ∈ [min(z 0 , 1), max(z 0 , 1)]. Hence, for each final time T , the maximum reward is achieved for z 1 such that J N (T, z 1 ) = 1 and the minimum when J N (T, z 1 ) = 0. We can see clearly that for small final times, the maximum is attained for a value of z 1 close to z 0 and that for z 1 = 1 the reward is the smallest. However, as the final time increases, the value of z 1 for which the reward is maximum approaches 1, and with the feedback ψ s(z 0 ) the reward is the smallest. In particular, we can see that the best of the feedbacks ψ s(z 1 ) depends on the final time. This leads us to consider a new feedback that keeps the system in the set of maximizers S = (s, z) ∈ D : s = s(z) .

(2.38)

We therefore introduce the following most rapid approach feedback to S

ψ S (s, z) = 0, if s > s(z), ū(s, z), if s = s(z), u max , if s < s(z), (2.39)
where ū(s, z) is the feedback that keeps the system in the set S, that we compute by differentiating with respect to time the equation s(t) = s(z(t)).

We first illustrate this feedback in Figure 2.5 where we show the states as functions of time and the open loop realizations of the feedbacks ψ s(z 0 ) , ψ s(1) and ψ S . Next, in Figure 2.6 we compare the reward of the feedback ψ S to the others and we can notice that the reward associated with the feedback ψ S is always one of the best, although for any given final time it is possible to do better with a feedback ψ s(z 1 ) for the right z 1 .

Note also that the feedback ψ S will drive the system asymptotically towards the state (s, z) = (s, 1) so that it is also optimal for the infinite horizon problems (2.14), (2.15) and (2.23).

In Figure 2.8, we show the difference between the rewards of the feedbacks ψ s(1) and ψ s(z 0 ) as a function of the initial condition for various final times. From this, we see that Difference between rewards associated to the feedbacks ψ s(1) and ψ s(z 0 ) as functions of the initial condition and for various final times : (x 0 , s 0 ) → J(0, x 0 , s 0 , ψ s(1) (•)) -J(0, x 0 , s 0 , ψ s(z 0 ) (•)). Contois growth function (µ max = 0.74, K s = 1) with s in = 100, u max = 1.5. the feedback that is best changes, depending on the initial condition and the horizon considered.

The sub-optimality estimation (2.32) is affected similarly, as this bound depends on the initial condition and in particular, the distance to the set {z = 1} has a major impact on the sub-optimality of the considered feedbacks. In addition, the growth function has an influence on our estimation, through W z 1 (•), and we illustrate this in Figure 2.9 by plotting this value function for the Haldane and the Contois growth function. Observe that, for the Contois growth function, W z 1 (•) varies significantly with the initial biomass and thus the sub-optimality bound as well. This can be attributed to the dependence of the Contois growth function on biomass concentration and this effect is not seen with the Haldane growth function, which depends only on the substrate.

We finish this section with a Lemma that shows that the considered growth functions satisfy our assumptions. Lemma 2.19 For all positive µ max , μ, K s and K i the Monod, Haldane and Contois growth functions satisfy Assumptions 2.1 and 2.17.

Proof. Notice that the function φ with the Monod or Haldane function does not depend on z. Let us show that the function µ M is increasing and strictly concave

µ M (s) = µ max K s (K s + s) 2 > 0, µ M (s) = -2 µ max K s (K s + s) 3 < 0.
Now, since the function φ(•, 1) is non-negative on [0, s in ] and vanishes at 0 and s in it admits a maximum on ]0, s in [. One has

d ds φ(s, 1) = µ M (s)(s in -s) -µ M (s), (2.40)
d 2 ds 2 φ(s, 1) = µ M (s)(s in -s) -2µ M (s). (2.41)
The function φ(•, 1) is thus strictly concave on ]0, s in [, which provides the uniqueness of its maximum.

For the Haldane function, we have

d ds φ(s, 1) = μs in K s -2K s s -s 2 (1 + s in K i ) (K s + s + s 2 K i ) 2 such that d
ds φ]0, 1) > 0 and d ds φ(s in , 1) < 0 and since d ds φ(•, 1) is continuous it must have an odd number of zeroes in the interval ]0, s in [. But notice that the equation d ds φ(s, 1) = 0 admits at most 2 solutions and φ]0, 1) = φ(s in , 1) = 0 and therefore φ(•, 1) has a unique maximum.

For the Contois function, notice that µ C (s, x) = µ M (s/x) so that,

φ(s, z 1 ) = µ M s (s in -s)z 1 (s in -s),
for z 1 ∈ [min(z 0 , 1), max(z 0 , 1)], and since s → s (s in -s)z 1 is an increasing function, φ(•, z 1 ) is also strictly concave.

Conclusions

In this work, we have proposed a novel approach to obtain autonomous sub-optimal feedbacks for the open problem of maximizing biogas production in the chemostat model out of equilibrium. These controllers generalize the "most-rapid approach path" feedback control that is known to be optimal when the initial condition belongs to a certain manifold. Indeed, we obtain a family of feedback controls of similar structure, for which we are able to give bounds on the sub-optimality. This last point merits to be underlined as it usually difficult to evaluate a priori the performances of sub-optimality without having to determine or compute the optimal solution. This choice gives also flexibility for the practitioners to choose a controller depending on the time horizon or simply to pick one when the finite horizon is poorly known (as each controller guarantees a sub-optimality bound), or to adjust it when the horizon is changed. For infinite horizon we show that each controller guarantees the same optimal averaged cost. This methodology, based on a framing of the dynamics, could be investigated for a larger class of dynamics, such as the two-step model, and be the matter of future work.

Appendix: A Particular Example

We construct here a control u(•) for which the average rewards (2.12) and (2.13) do not coincide. For this, let us consider an initial condition ξ = (s 0 , z 0 ) = (ε, 1), with ε ∈]0, s in [ fixed. The set {(s, 1) ∈ R 2 + : s ∈ [0, s in ]} is clearly invariant for the dynamics (2.3) and therefore the chosen initial condition ensures that trajectories (s ξ,u (•), z ξ,u (•)) remains in this set.

Now consider the 2 following paths :

(A) Starting at ξ := (ε, 1), use the control u = u max to reach a prescribed level of substrate s * ∈ (ε, s in [ in finite time. Then, apply the control u = 0 to return to ξ in finite time, which is possible by Assumption 2.3. Denote this control by u * , and let t * be the (finite) time necessary to follow this path and I * be the biogas produced by this path.

(B) Starting at ξ := (ε, 1), use u = µ(ε, s in -ε) to stay at (s = ε, z = 1) for any time interval.

Then, define control u(•) as follows:

• For t ∈ [0, t * ], set u(t) = µ(ε, s in -ε) so that the biogas production for this period is I ε := t * φ(ε, 1).

• For t ∈]2 2k t * , 2 2k+1 t * ], with k ∈ N, set u = u * in order to follow the path (A) repeatedly 2 2k times. For each of these intervals the biogas production is 2 2k I * .

• For t ∈]2 2k+1 t * , 2 2k+2 t * ], with k ∈ N, set u = µ(ε, s in -ε).
For each of these intervals the biogas production is 2 2k+1 I ε .

Thus, when we apply control u(•) up to a time 2 2N t * , for a given N 1, the average biogas production is computed as follows

K N = 1 2 2N t * 2 2N t * 0 φ(s ξ,u (t), 1) dt = 1 2 2N t * I ε + N -1 k=0 2 2k I * + N -1 k=0 2 2k+1 I ε = I * + 2I ε t * N j=1 2 -2j + I ε 2 2N t * which yields K N -→ K ∞ := I * + 2I ε 3t * as N → +∞.
We have used here the fact that the sum s N = N j=1 2 -2j converges to 1/3. Indeed, this follows from the identity

4s N = N j=1 2 2(-j+1) = N -1 i=0 2 -2i = 1 + s N -2 -2N .
However, for the same control u(•), the average biogas production is, up to time 2 2N +1 t * , computed as follows

L N = 1 2 2N +1 t * 2 2N +1 t * 0 φ(s ξ,u (t), 1) dt = 1 2 2N +1 t * 2 2N t * K N + 2 2N I * = 1 2 K N + I * t * which yields L N -→ L ∞ := 2I * + I ε 3t * as N → +∞.
Since s * > ε, it follows that I * > I ε , and consequently, L ∞ > K ∞ . We thus obtain that

J ∞ (ξ, u(•)) L ∞ > K ∞ J ∞ (ξ, u(•)).
Chapter 3

An Algorithm for Maximizing the Biogas Production in a Chemostat

This chapter corresponds to the published article Haddon, A. and Hermosilla, C. An Algorithm for Maximizing the Biogas Production in a Chemostat.

J Optim Theory Appl (2019) 182:1150. https://doi.org/10.1007/s10957-019-01522-x

Introduction

Biogas is a product of the anaerobic digestion process, in which several populations of microorganisms break down organic matter in the absence of oxygen. This process is an interesting technology for the treatment of liquid and solid waste since the collected biogas is mainly composed of methane and therefore can be used as a renewable energy source [START_REF] Rehl | CO2 abatement costs of greenhouse gas (GHG) mitigation by different biogas conversion pathways[END_REF]. In this context, it is relevant to develop control strategies that maximize methane production, in order to increase the efficiency and sustainability of waste treatment. As a matter of fact, a major reason that has been reported for the closing of anaerobic digestion plants, is the insufficient profits associated with poor biogas production [START_REF] Beddoes | An analysis of energy production costs from anaerobic digestion systems on us livestock production facilities[END_REF]. Substantial expertise is needed to operate the anaerobic digestion process properly as it is a complex non-linear and unstable process. Although it is possible to use various inputs for control, such as pH or alkalinity [START_REF] Guwy | Neural network and onoff control of bicarbonate alkalinity in a fluidised-bed anaerobic digester[END_REF], the dilution rate (also called feeding rate) is considered in general as the variable input. It is important to note that most studies on the control of anaerobic digestion have focused primarily on process stability [START_REF] Nguyen | Automatic process control in anaerobic digestion technology: A critical review[END_REF]. However, recently, some works have incorporated the aspect of optimizing performance and, among these, a wide range of control strategies have been used: PID controllers [START_REF] García-Diéguez | Multi-objective cascade controller for an anaerobic digester[END_REF], expert systems [START_REF] Rodrıguez | A hydrogenbased variable-gain controller for anaerobic digestion processes[END_REF], fuzzy logic [START_REF] Djatkov | Method for assessing and improving the efficiency of agricultural biogas plants based on fuzzy logic and expert systems[END_REF] and adaptive control [START_REF]Nonlinear adaptative stabilizing control of an anaerobic digestion model with unknown kinetics[END_REF], to mention a few strategies.

In this work, we address the problem of biogas production from an optimal control point of view. We focus our attention on the one reaction model in a chemostat. We are particularly interested in providing a practical method to determine an optimal control in feedback form for maximizing the production of biogas. The numerical scheme we propose for such purpose (Algorithm 1) has been obtained by combining the two major techniques in optimal control, namely, the Pontryagin Maximum Principle (PMP) and the Hamilton-Jacobi-Bellman (HJB) equation. The PMP allows us to describe the structure of the optimal synthesis and the singular curve in terms of a given parameter, while the HJB equation gives a practical way to compute such a parameter.

As far as we are aware of, there are few works dealing with the dynamic optimization problem of biogas production. Actually, due to the complexity of the problem, only models with one or two bio-reactions have been considered [START_REF] Sbarciog | An optimizing start-up strategy for a bio-methanator[END_REF][START_REF] Bayen | Analysis of an optimal control problem related to the anaerobic digestion process[END_REF]. In addition, only problems on well-mixed continuously stirred tanks have been studied, since the nonlinearities and the high dimension of a more complex bio-reactor model make the analysis of the associated optimal control problem hard to handle. It is worth mentioning that models with only few dynamic variables are capable of describing the qualitative behaviour of the anaerobic digestion process [START_REF] Bernard | Can we assess the model complexity for a bioprocess: theory and example of the anaerobic digestion process[END_REF]. The tradeoff between practical solvability and qualitative description justifies the use of these simplified models, and in particular the one reaction model we study in this paper.

The problem for a one reaction model was first considered in [START_REF] Stamatelatou | Optimal and suboptimal control of anaerobic digesters[END_REF] and later solved for a special set of initial conditions for which the model reduces to a one dimensional problem [START_REF] Ghouali | Maximizing biogas production from the anaerobic digestion[END_REF]. More recently, the general one reaction model has been revisited to propose a suboptimal control for which there is an estimation of sub-optimality [START_REF] Haddon | Guaranteed value strategy for the optimal control of biogas production in continuous bio-reactors[END_REF]. Let us mention that the problem has also been considered in the infinite horizon case [START_REF] Haddon | First results of optimal control of average biogas production for the chemostat over an infinite horizon[END_REF]. To the best of our knowledge, a complete synthesis for the problem of maximizing biogas production over a fixed finite horizon has not been addressed before, even for the single reaction model. This work contributes in this direction, by proposing a candidate to optimal synthesis and giving a practical way to compute it.

Notice that for general optimal control problems a wide range of algorithms have been studied and implemented as open source software. Either based on Shooting methods, Dynamic Programming or Discretize-then-optimize methods such as Nonlinear Model Predictive Control ; see for example [START_REF]Bocop: an open source toolbox for optimal control[END_REF][START_REF] Gerdts | Optimal Control and Parameter Identification with Differential-Algebraic Equations of index 1: User's Guide[END_REF][START_REF] Bokanowski | a C++ Library for Solving HJ equations[END_REF]114]. The problem we study in this paper can in principle be solved numerically with any of these methods, provided that one knows exactly the parameters of the model (which are hard to estimate in practice). However, as pointed out earlier, the goal of studying simplified models is to provide a good picture of how an optimal synthesis may look like (qualitative description) rather than giving a specific solution for the maximization of biogas production problem. This paper is organized as follows. In Section 2 we describe the problem at hand. In Section 3 we analyze the optimality conditions and we identify a class of extremal controls. In Section 4, we explain the algorithm we propose and we provide some numerical simulations in Section 5. Finally, in the appendix, we give an analytic proof (based on the HJB approach) for the optimality of the feedback law we found in Section 3 for a special set of initial conditions.

Problem Statement

In this work, we consider a single reaction model of the anaerobic digestion process where a substrate of concentration s is transformed by a microbial population of concentration x into biogas. The bioreactor is assumed to be continuously-fed and well-mixed, for which the mass balance equations are the classical chemostat equations [START_REF] Lobry | The Chemostat: Mathematical Theory of Microorganism Cultures[END_REF] :

ẋ = µ(s)x -Dx, ṡ = D(s in -s) -µ(s)x, (3.1)
where s in > 0 is the substrate inflow concentration and D is the dilution rate, which will be the controlled variable (it is assumed to be a measurable function of time). We suppose here, without loss of generality, that the units are chosen such that the yield coefficient of the reaction is equal to 1.

The specific growth rate of the microorganisms µ(•) is usually chosen of Monod type (µ M ) or of Haldane type (µ H ): given

µ max , K, K i > 0 µ M (s) := µ max s K + s , µ H (s) := µ max s K + s + s 2 /K i . (3.2)
However, at first we will study the optimal control problem for a rather general class of functions, which in particular covers the Monod and Haldane cases.

Standing Assumptions: The growth rate of the microorganisms µ(•) is a twice continuously differentiable function on [0, +∞[ such that µ(0) = 0, µ(s) > 0 and d ds

µ (s) µ(s) 2 = 0, ∀s > 0.
In the Monod and Haldane cases we have that

µ M (s) µ M (s) 2 = K µ max s 2 and µ H (s) µ H (s) 2 = K µ max s 2 - 1 µ max K i
respectively. Thus in particular, they satisfy our Standing Assumptions.

The biogas flow-rate is assumed proportional to the growth rate of the microorganisms [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF] and therefore the biogas production during a time interval [t 0 , T ] for a given substrate concentration s(•) and a given microbial population concentration x(•) is

T t 0 µ(s(t))x(t)dt.
The goal of the problem we deal with here is to maximize the biogas production over a finite horizon [t 0 , T ] for a given initial condition x 0 , s 0 > 0 by controlling the dilution rate t → D(t) of the bioreactor under the constraint that D(t) ∈ [0, D max ], where D max > 0 is the maximal dilution rate allowed.

In summary, the optimal control problem that we will study is the following

                   Maximize T t 0 µ(s(t))x(t)dt over all D : [t 0 , T ] → [0, D max ] measurable such that ẋ = µ(s)x -Dx, x(t 0 ) = x 0 , ṡ = D(s in -s) -µ(s)x, s(t 0 ) = s 0 , 0 ≤ s(t) ≤ s in and 0 ≤ x(t), ∀t ∈ [t 0 , T ]. (P bio )
With a slight abuse of notation, we may sometimes write

J(t 0 , x 0 , s 0 , ψ) := T t 0 µ(s(t))x(t)dt
for the biogas production associated with a control ψ : R 2 → [0, D max ] in feedback form and some initial conditions (x 0 , s 0 ) = (x(t 0 ), s(t 0 )). Under these circumstances, the functions s(•) and x(•) denote a solution of control system (3.1) in closed-loop form associated with these data.

Remark Since feedback controls ψ : R 2 → [0, D max ] are not necessarily continuous functions of the state variables, the classical theory of ordinary differential equations (ODEs) cannot be evoked for ensuring the existence and uniqueness of solutions to the ODEs system:

ẋ = µ(s)x -ψ(s, x)x, ṡ = ψ(s, x)(s in -s) -µ(s)x.
In our setting, the feedback controls are going to be regular enough to ensure the well-posedness (existence and uniqueness) of the control system (3.1) in closed-loop form. This is due to the fact that the feedback controls considered later on have an underlying stratified structure and so they can be handled with a tailored ODEs theory; see for instance [START_REF] Hermosilla | Stratified discontinuous differential equations and sufficient conditions for robustness[END_REF].

About the State-Constraints

Let us point out that in the formulation of the problem we have included state-constraints over the system, described by a set K := [0, +∞[×[0, s in ]. In the rest of the paper this restriction will be disregarded. The main reason for doing so is that the system (3.1) is invariant on K (see for example [START_REF] Clarke | Nonsmooth analysis and control theory[END_REF]Theorem 4.3.8]). Indeed, the set-valued map

F (x, s) := {(µ(s)x -Dx, D(s in -s) -µ(s)x) : D ∈ [0, D max ]}
is locally Lipschitz continuous, has linear growth, has nonempty compact and convex images and satisfies the invariance condition F (x, s) ⊆ T K (x, s) for any (x, s) ∈ K; where T K stands for the Contingent Cone. The last affirmation comes from the fact that F (0, s) ⊆ {0} × [0, +∞[ for any s ∈ [0, s in ] and

F (x, 0) ⊆] -∞, 0] × [0, +∞[, F (x, s in ) ⊆ R×] -∞, 0[, ∀x > 0.
By similar arguments we can see that the set {0}×[0, s in ] is also invariant, which means that no trajectory of the system that starts from x(t 0 ) = x 0 > 0 and s(t 0 ) = s 0 ∈]0, s in [ will reach that set. Moreover, the fact that F (x, s in ) is contained in R×] -∞, 0[ for any x > 0 implies that no trajectory can reach the level s = s in provided that s(t 0 ) < s in . Also, note that for any x ≥ 0 we can find s x > 0 small enough such that D(s in -s x )-µ(s x )x ≥ 0 for any D ∈]0, D max ]. Moreover, a trajectory of (3.1) associated with D = 0 cannot reach the level s = 0 in finite time, otherwise there would be two backward solutions to the corresponding ODE starting from the same point. In practice, this means that whenever the initial conditions are taken such that s 0 ∈]0, s in [ and x 0 > 0, we will have that the condition over the states of the system holds, and is even stronger, in the sense that we will also have that 0 < s(t) < s in and 0 < x(t), ∀t ∈ [t 0 , T ].

(3.3)

Optimality Conditions

The preceding discussion implies in particular that admissible trajectories exist for the optimal control problem (P bio ). Furthermore, since the objective function to be maximized does not depend explicitly on the control function D(•), standard assumptions that guarantee the existence of optimal control can be evoked (for example [START_REF] Clarke | Functional analysis, calculus of variations and optimal control[END_REF]Theorem 23.11]). Thus in the rest of the paper we might assume that optimal trajectories for the maximization of biogas production problem exists and focus on optimality condition to understand and approximate such solutions. Also, since we are mainly interested in the case that the initial conditions are such that s 0 ∈]0, s in [ and x 0 > 0, we will assume, unless otherwise stated, that optimal trajectories satisfy (3.3).

Pontryagin Maximum Principle

We begin our study of problem (P bio ) by establishing necessary conditions of optimality with the Pontryagin Maximum Principle (PMP) [START_REF] Clarke | Functional analysis, calculus of variations and optimal control[END_REF]Corollary 22.3]. For this we set the Hamiltonian H :

R 2 × R 2 × [0, D max ] → R as H(x, s, p x , p s , D) := µ(s)x + p s (D(s in -s) -µ(s)x) + p x (µ(s)x -Dx). (3.4)
We consider an optimal control D(•) of (P bio ) and its associated states x(•) and s(•), solution of (3.1) with initial condition (x 0 , s 0 ) = (x(t 0 ), s(t 0 )). Then, the PMP states that there exist adjoint states p s , p x : [t 0 , T ] → R satisfying, for almost every t ∈ [t 0 , T ], the adjoint equations

ṗx = Dp x -µ(s)(1 + p x -p s ), ṗs = Dp s -µ (s)x(1 + p x -p s ), (3.5)
the transversality condition p x (T ) = p s (T ) = 0 and the maximum condition, for almost every t ∈ [t 0 , T ],

H(x(t), s(t), p x (t), p s (t), D(t)) = max d∈[0,Dmax] H(x(t), s(t), p x (t), p s (t), d). (3.6)
In addition, since the Hamiltonian does not depend explicitly on time, it is constant, which means that for some c = c(t 0 , x 0 , s 0 ) ∈ R we have Since the Hamiltonian is affine in the control variable, an extremal control will depend on the sign of the commutation function

H(x(t), s(t), p x (t), p s (t), D(t)) = c, a.e. on [t 0 , T ]. ( 3
φ(t) := ∂ ∂D H(x(t), s(t), p x (t), p s (t), D(t)) = p s (t)(s in -s(t)) -p x (t)x(t).
We then have that D(t) = 0 if φ(t) < 0 and D(t) = D max if φ(t) > 0, while no information can be directly obtained from the PMP in the case φ(t) = 0.

We recall that a singular arc is a time interval during which we have φ(t) = 0 and since this equation is valid along a singular arc, we also have d dt φ(t) = 0. Therefore, during a singular arc the state variables and the adjoint states satisfy the following equations

(s in -s)p s -xp x = 0, µ (s)(s in -s)(1 + p x -p s ) = µ(s). (3.8)
With this, we can get an equation that the state variables satisfy during a singular arc, that depends only on the constant value c of the Hamiltonian. Indeed, when the commutation function vanishes, we have c = µ(s)x(1 + p x -p s ) and using (3.8) we get

c µ (s)(s in -s) = µ(s) 2 x. (3.9)
We now define the following function, for h > 0 given

x h (s) := h µ (s)(s in -s) µ(s) 2 , 0 < s < s in .
Then, from (3.9), we have that the extremal state trajectories during the singular arc remain in the graph of s → x c (s). This means that if we knew the value of c then we would be able to construct the singular arc and construct an admissible extremal control in feedback form for the optimal control for problem (P bio ).

Remark In the Monod and Haldane cases, we have that the curve described above has, respectively, the form

x M h (s) := hK(s in -s) µ max s 2 , x H h (s) := hK(s in -s) µ max s 2 - h(s in -s) µ max K i .
Remark Note that in general the curve x h (s) → 0 when s → s in , and x h (s) → +∞ if s → 0 provided that µ (s) µ(s) 2 → +∞ as s → 0; this is for instance the case of the Monod and Haldane growth rate functions.

Construction of Extremal Controls

To identify the extremal controls, we start by constructing explicitly a control that drives the system to a singular arc associated with a given h > 0.

We first need to compute the control D h that keeps the system on the singular curve {(x h (s), s) : 0 < s < s in }. For this, we differentiate with respect to time the relation

x(t) = x h (s(t)) to get µ(s)x h -D h x h = D h (s in -s) -µ(s)x h ∂ s x h
and we then have the following expression for the control on the singular arc

D h (s) = µ(s)x h (s)(1 + ∂ s x h (s)) x h (s) + (s in -s)∂ s x h (s) , (3.10)
with

∂ s x h (s) = h d ds µ (s) µ(s) 2 (s in -s) - µ (s) µ(s) 2 , ∀s ∈]0, s in [. (3.11)
Our Standing assumptions, in particular the fact that d ds µ

= 0 implies that D h (s) is well defined (as a real valued function) for any s ∈]0, s in [. The control D h is not necessarily an admissible control for the problem at hand. For some s ∈]0, s in [ it could happen that D h (s) ∈ [0, D max ]. We assume that the singular control D h (s) is admissible only on a bounded interval I s where the bounds s 0 and s max are defined as the solutions of

D h (s 0 ) = 0 and D h (s max ) = D max . (3.12)
Remark It is straightforward to see that the control on the singular curve associated with a Monod growth rate function is given by

D M h (s) = hK(2s in -s) -µ max s 3 2s(K + s)(s in -s) , ∀s ∈]0, s in [. It follows that D M h (s) → +∞ if s → 0.
However, the behavior of D M h (s) when s → s in depends on the data of the problem. As a matter of fact lim

s→s in D M h (s) =      +∞, if hK > µ max s 2 in , -∞, if hK < µ max s 2 in , 2µmaxs in K+s in , if hK = µ max s 2 in .
This means that, depending on the data of the problem, singular optimal trajectories may not occur at all; for instance if D M h (s) > D max for every s ∈]0, s in [. We plan to study this issue in more details and for general growth rate functions elsewhere. For (s, x) ∈ I s × x h (I s ), we extend the singular curve such that the control to stay on that curve is equal to 0 or D max . For this we integrate the dynamics backwards with D = 0 (respectively D = D max ) starting from s 0 (respectively s max ). We therefore have the following singular curve :

G h :=    x h (s), s : s ∈ I s x(τ, x h (s 0 ), 0), s(τ, s 0 , 0) : τ 0 x(τ, x h (s max ), D max ), s(τ, s max , D max ) : τ 0    (3.13)
where we denote s(τ, s max , D max ) the value at time τ of the solution with control D max starting at s max at time τ = 0 and similarly for x(τ, x h (s max ), D max ), x(τ, x h (s 0 ), 0) and x(τ, x h (s 0 ), 0). Note that with D = 0 we have ẋ+ ṡ = 0 so that the trajectory { x(τ, x h (s 0 ), 0), s(τ, s 0 , 0) : τ 0} corresponds to the graph of the mapping s → -s+s 0 +x h (s 0 ). This is a decreasing function of s and with D = 0 the trajectories are such that s(t) is also decreasing and therefore { x(τ, x h (s 0 ), 0), s(τ, s 0 , 0) : τ 0} corresponds to the set

-s + s 0 + x h (s 0 ), s : s 0 < s < s in , -s + s 0 + x h (s 0 ) > 0 . (3.14)
The singular curve G h divides the state space in 2 sets on which the control must be either 0 or D max and we thus denote G 0 h (respectively G max h ) the set on which the control is 0 (respectively D max ). To distinguish these sets we use again the fact with D = 0 the trajectories are such that s(t) is decreasing and therefore either the trajectory reaches the singular curve or approaches asymptotically the set {(x, 0) : x > 0}. This corresponds to determining whether there exists a point (-s + x 0 + s 0 , s) that belongs to G h for s ∈]0, s 0 ]. We then have the following expression for G 0

h G 0 h = {(x, s) : ∃ s s such that (-s + x + s, s) ∈ G h } and then G max h := R 2 + \ (G 0 h ∪ G h ).
With this, we have the following family of feedback controls

ψ h (x, s) =      0, if (x, s) ∈ G 0 h , D max , if (x, s) ∈ G max h , D h (s), if (x, s) ∈ G h . (3.15)
An example of the trajectoires obtained with this feedback is show in Figure 3.1. Note that, because of the way the feedback control ψ h has been constructed, solution of related the closed-loop system do exist and are uniquely determined by the initial data.

Hamilton-Jacobi-Bellman Equation

In order to make the feedback law (3.15) a suitable candidate to optimal control, we now need to identify the value of the Hamiltonian for a given initial condition and initial time, that is, we need to calculate or approximate c. For this purpose, we use the Hamilton-Jacobi-Bellman (HJB) equation, motivated by the fact that the cost associated with an optimal control, seen as a function of the initial data, can be completely characterized by an appropriate HJB equation.

The value function for the production of biogas problem (without state constraints) is

V (t 0 , x 0 , s 0 ) := sup D(•) T t 0 µ(s(t))x(t)dt : ẋ = µ(s)x -Dx, x(t 0 ) = x 0 , ṡ = D(s in -s) -µ(s)x, s(t 0 ) = s 0
where the maximum is taken over all D : [t 0 , T ] → [0, D max ] measurable. It is not difficult to see that, thanks to the continuity of the trajectories of the control system (3.1) with respect the initial data, the value function (t 0 , x 0 , s 0 ) → V (t 0 , x 0 , s 0 ) is continuous. We have already discussed that optimal controls do exist, and then the supremum is actually a maximum. Furthermore, because of the invariance of the set [0, +∞[×[0, s in ] with respect to the control system (3.1), this value function agrees with the value function of the original problem (P bio ) with state constraints. Let us mention that problems with state constraints are considerably harder to deal with and so the fact stated above simplifies considerably the ensuing analysis (cf. [START_REF] Hermosilla | Infinite horizon problems on stratifiable stateconstraints sets[END_REF]).

The HJB equation for the problem we are dealing with is

∂ t u + sup D∈[0,Dmax] H(x, s, ∂ x u, ∂ s u, D) = 0, in ] -∞, T [×R 2 , (3.16)
where H is the Hamiltonian given in (3.4). Existence and uniqueness of solutions in the viscosity sense for HJB equations is a well-known and studied fact, see for instance [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF].

As a matter of fact, the value function V is the unique viscosity solution to (3.16) that satisfies the terminal condition

u(T, x, s) = 0, ∀x, s ∈ R.
Furthermore, the HJB equation and the PMP are related via the following lemma, which links the derivatives of the value function with the adjoint arcs.

(c, p x (t), p s (t)) ∈ ∂ + V (t, x(t), s(t)), a.e. on [t 0 , T ].
where c = H(x(t), s(t), p x (t), p s (t), D(t)) for a.e. t ∈ [t 0 , T ] and

∂ + u(z) := q ∈ R n : lim sup y→z u(y) -u(z) -q (y -z) |y -z| ≤ 0
stands for the viscosity superdifferential of a function u : R n → R.

The preceding lemma implies that whenever the value function is differentiable at (t 0 , x 0 , s 0 ), we should have that (p x (t 0 ), p s (t 0 )) = ∇ (x 0 ,s 0 ) V (t 0 , x 0 , s 0 ) and c = -∂ t 0 V (t 0 , x 0 , s 0 ). This fact is the key point we use for proposing an algorithm for solving the production of biogas problem. Indeed, we have seen that the control that maximizes the Hamiltonian is ψ h given by (3.15), where h > 0 is the (constant) value of the Hamiltonian. The value of the Hamiltonian can be obtained, for example, by evaluating at initial time:

h = H(x 0 , s 0 , p x (t 0 ), p s (t 0 ), ψ h (x 0 , s 0 ))
From the previous section, we can deduce that ψ c is the optimal control associated with optimal singular trajectories (it is the unique candidate to be an extremal control in this case). Also, in the appendix we show that ψ c is the optimal control for a particular choice of initial conditions1 . If ψ c is actually an optimal control, we can write the value function as the cost of the control ψ c , that is, V (t 0 , x 0 , s 0 ) = J(t 0 , x 0 , s 0 , ψ c ). With Lemma 3.1 we then get c = H(x 0 , s 0 , ∇ (x 0 ,s 0 ) J(t 0 , x 0 , s 0 , ψ c ), ψ c (x 0 , s 0 ))

Thus, using the HJB equation (3.16), we also have

c = -∂ t 0 J(t 0 , x 0 , s 0 , ψ c ) (3.17)
In other words, c is a fixed point of the mapping h → -∂ t 0 J(t 0 , x 0 , s 0 , ψ h ). Hence, if we are able to compute or approximate a fixed point of the mapping h → -∂ t 0 J(t 0 , x 0 , s 0 , ψ h ) we will be able to reconstruct an optimal synthesis for the production of biogas problem.

In the next section we will present an algorithm, one of the main contribution of this work, based on a classical iterative scheme for finding a fixed point of a function by repeatedly computing the image of the previous iterate.

Let us point out that the HJB equation is valid regardless of the structure of the optimal control. This means that this equation can be seen as a certificate of optimality for the feedback control (3.15), in the sense that if the algorithm converges, then the proposed feedback control is a good approximation of an optimal control, because the value function obtained with the feedback ψ c is an approximated solution to the HJB equation.

An Algorithm for Maximizing the Production of Biogas

We present now a way to compute the extremal feedback control ψ c by solving the fixed point equation (3.17) numerically, in order to get the value of the Hamiltonian c for any initial condition (x 0 , s 0 ) ∈]0, +∞[×]0, s in [.

HJB Fixed Point Algorithm

The algorithm we propose is based on a classical iterative scheme for finding a fixed point of a function by repeatedly computing the image of the previous iterate. More precisely, if the equation to be solved is F (h) = h for some given mapping F : R d → R d , then starting from an initial guess h 0 , the numerical scheme consists in computing h n+1 = F (h n ) for n = 0, 1, 2, ... The algorithm is then considered to have converged to a fixed point when the iterates stabilize to a given tolerance ε, specifically when ||h n+1 -h n || < ε.

In our case, the function for which we need to compute a fixed point, -∂ t 0 J(t 0 , x 0 , s 0 , ψ h ), is composed of a partial derivative and therefore to estimate it numericaly with a finite stoping condition that guarantees the convergence of the last fixed point to the desired tolerance.

On the other hand, by computing the whole vector h n+1 at every iteration, we can stop the algorithm when the vector has converged for the maximum norm, that is, when

max k |h k n+1 -h k n | < ε.
This guarantees that all components of the vector have converged to a desired tolerance. In addition, the algorithm will keep iterating on the first components (h N -1 n , h N -2 n , ...), which converge to the desired tolerance first, but as such it will keep on reducing the errors automatically to achieve convergence of the last components (h 1 n , h 2 n , ...).

Initial Guess

Concerning the initial guess h 0 , recall that with the optimal control we have c = µ(s(T ))x(T ).

In most cases, we can broadly approximate this by taking h k 0 = µ(s 0 )x 0 , for all k, and Algorithm 1 will converge.

However, in the most difficult cases, this is not sufficient and the accumulation of errors that were previously mentioned can cause the algorithm to diverge. In fact, the only problems we encountered were when the solution has a bang-bang-singular arc control with a switch from D = 0 to D = D max before reaching the singular arc. To deal with these cases, we propose to first to identify an extremal candidate by solving the fixed point equation c = µ s(T, t 0 , x 0 , s 0 , ψ c ) x(T, t 0 , x 0 , s 0 , ψ c ) Then using the obtained fixed point as an initial guess for Algorithm 1, we can check the optimality of the associated extremal candidate.

In this case, as we do not need to compute the partial derivative of the cost with respect to initial time we do not need to compute the fixed points for a range of initial times simultaneously. Other than this, the algorithm to solve this equation is similar to the previous and is shown below as Algorithm 2.

Algorithm 2:

Input:

MaxIterations, ε h 0 ← µ(s 0 )x 0 for n < MaxIterations do Solve D hn (s 0 ) = 0 and D hmax (s max ) = D max Compute singular curve G hn Compute x(•, t 0 , x 0 , s 0 , ψ hn ), s(•, t 0 , x 0 , s 0 , ψ hn ) h n+1 ← µ(s(T, t 0 , x 0 , s 0 , ψ hn ))x(T, t 0 , x 0 , s 0 , ψ hn ) if |h n+1 -h n | < ε then return h n+1

Numerical Simulations

In this section, we illustrate the fixed point algorithm with the growth functions of Monod and Haldane (3.2) with parameter values from [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF].

In Figures 3.2 and 3.3, we show examples of optimal trajectories in state space for various initial conditions but with the same initial and final times for each growth function. The solutions are similar for both growth functions and there is both bang-singular arc and bang-bang-singular arc optimal solutions. Note that, the singular curve varies for each initial condition but that all trajectories that reach a singular arc finish with

s(T ) = s * = arg max s∈[0,s in ] µ(s)(s in -s).
This is expected since c = µ(s(T ))x(T ) and using the expression for the singular curve (3.9) evaluated at final time we get µ (s(T ))(s in -s(T )) = µ(s(T )) and we recognize this as a necessary condition for maximizing s → µ(s)(s in -s).

Next, Figures 3.4 initial times but for fixed initial conditions and final time. We can see that the singular curve varies for different initial times and the strategy can also change. For instance, in Figure 3.5, we can see that for t 0 = 1.875 the optimal trajectory corresponds to the control D = 0 and as the initial time decreases, the optimal control switchs to a bangbang-singular arc with first D = D max and then D = 0 before reaching the singular arc.

Alongside each set of trajectories is also shown the corresponding values of the Hamiltonian as a function of the initial time. Although the function t 0 → c(t 0 ) appears to be continuous it is clearly not continuously differentiable everywhere and the points at which this function is not smooth correspond to initial times when there is a change in the type of control. For example, in Figure 3.4 we can see on the state space trajectories graph that the optimal control for t 0 = 0.5 is a single bang D = 0 whereas for t 0 = 0.25 it is bang-bang with a switch from D = 0 to D = D max and on the graph of t 0 → c(t 0 ) there is indeed a point of irregularity near t 0 = 0.4. We illustrate the performance and convergence of Algorithm 1 and 2 in Figure 3.7 with graphs of error as function of iterations :

n → |h k n+1 -h k n |.
For Algorithm 1, we can see that for initial times close to the final time, the convergence is very fast. However, as the horizon increases, not only is convergence slower but there is a limit for the errors and eventually they stop decreasing. This is likely due to the numerical errors when computing the finite difference approximation of ∂ t 0 J(t 0 , x 0 , s 0 , ψ h ) since this behaviour is not seen for Algorithm 2 which does not need the computation of ∂ t 0 J(t 0 , x 0 , s 0 , ψ h ). Notice however, that Algorithm 2 requires more iterations but that the convergence accelerates at the end and reaches machine precision.

Comparison with Bocop

Finally, we compare our feedback to the control obtained with the open source toolbox for optimal control Bocop [START_REF]Bocop: an open source toolbox for optimal control[END_REF][START_REF] Bonnans | Bocop -A collection of examples[END_REF]. This package implements a direct method that approximates the optimal control problem by a finite dimensional optimization problem using a time discretization. Table 3.1 presents a performance comparison by looking at the biogas production of each control and the relative difference. We can see that our feedback achieves nearly as much as Bocop and that the difference is greater for the last 2 rows which correspond to trajectories that are bang-bang-singular arc as in Figures 3.4 Next, in Table 3.2 we show some computational times associated with our feedback (for various error tolerances ε) and Bocop. We can see that the time necessary to compute our feedback is similar to the time reported by Bocop, although it is important to note that Bocop only computes the control for a single initial time, where as our algorithm for a range of initial times.

Conclusions

In this work, we have given an algorithm to compute an extremal control for the problem of maximizing biogas production for the classical model of the chemostat for a fixed finite horizon. The extremal control is obtained in state feedback form which has advantages in terms of robustness with respect to pertubations on the initial data. In order to achieve this we first studied necessary optimality conditions thereby obtaining an analytical expression of a family of extremal feedbacks. Then we use a sufficient optimality condition (the HJB equation) to single out one of the extremal feedbacks as a candidate to optimal control. The resulting algorithm is fast and converges rapidly in practice. As pointed out before, the HJB equation can be seen as a test of optimality for the proposed feedback control (3.15), in the sense that if the algorithm converges, then the proposed feedback control is a good approximation of an optimal control. This fact has also been corroborated with the numerical examples we have exhibited and the comparison done with Bocop.

Let us finally mention that the technique we have introduced in this paper is well suited for the one reaction model. Some extensions to more general cases, such as two reactions models, should be possible. This is work in progress.

Appendix: Reduced Model

In this final, part we provide a HJB proof for the optimal synthesis for the reduced model, that is, the case where the initial data satisfy s in = x 0 +s 0 . In particular, we show how the fixed point characterization of the optimal control can be used analytically in a special case when the dynamics reduces to a single equation. A well known property of the chemostat model is that the set I := {(x, s) ∈ R : x + s = s in } is invariant for the dynamics (3.1) and thus, for initial conditions in I, the dynamics reduce to ṡ = Dµ(s) (s in -s). This special case was solved in [START_REF] Ghouali | Maximizing biogas production from the anaerobic digestion[END_REF], with the following assumptions (H1) The function s → µ(s)(s in -s) has a unique maximizer s * on [0, s in ].

(H2) The upper bound on the controls is such that D max > µ(s * ).

The optimal control is then D * (s) = 0 if s > s * , D * (s) = D max if s < s * and D * (s) = µ(s * ) if s = s * . Here, we will give another proof of the optimality of this control, by using the fixed point characterization. First, we can identify the control D * as a control of the type (3.15) where the singular arc is reduced to s = s * . In other words, it corresponds to the control ψ h * where h * satisfies equation (3.9) for the singular arc with s = s * , which in this case is h * µ (s * ) = µ(s * ) 2 . Next, since s * is a maximizer we have µ (s * )(s in -s * )-µ(s * ) = 0 and therefore h * = µ(s * )(s in -s * ).

To prove the optimality of ψ h * , we must now show that h * is a fixed point of the mapping h → -∂ t 0 J(t 0 , x 0 , s 0 , ψ h ). For this we first study the trajectories obtained with the feedback control ψ h * . We denote in the remainder of the section the right-hand side of the differential equation for s(•) with control ψ h * as f (s) := (ψ h * (s) -µ(s))(s in -s).

Notice that for s > s * we have f (s) = -µ(s)(s in -s) < 0 and for s < s * we have f (s) = (D max -µ(s))(s in -s) > 0 from assumption (H2). Thus, s * is reachable from any initial condition in I with control ψ h * . We define the time t * when s * is reached, from a given initial condition s 0 ∈ [0, s in ] and initial time t 0 with control ψ h * , that is, t * := inf {t t 0 : s(t, t 0 , s 0 , ψ h * ) = s * } . Finally, note that with control D = µ(s * ) the point s = s * becomes a steady state. Therefore the trajectories with control ψ h * are s(t) = s(t, t 0 , s 0 , ψ h * ), for t 0 t min(t * , T ), s * , for min(t * , T ) t T.

We can now compute ∂ t 0 J(t 0 , x 0 , s 0 , ψ h * ) and for this we need the following. We are now in a position to prove the optimality of the feedback control proposed earlier.

Proposition 3.3 For any initial condition (x 0 , s 0 ) ∈ I and for any initial time t 0 such that s * is reachable, that is when t * T , we have ∂ t 0 J(t 0 , x 0 , s 0 , ψ h * ) = -µ(s * )(s in -s * ), so that ψ h * is the optimal control.

Proof. We start by writing the cost as

J(t 0 , x 0 , s 0 , ψ h * ) = t * t 0 µ(s(t))(s in -s(t)) dt + (T -t * )µ(s * )(s in -s * )
differentiating with respect to t 0 we get

∂ t 0 J(t 0 , x 0 , s 0 , ψ h * ) = -µ(s 0 )(s in -s 0 ) + t * t 0 ∂ s µ(s(t))(s in -s(t)) ∂ t 0 s(t) dt.
Note that the terms with ∂ t 0 t * cancel out because s(t * ) = s * . Now, using Lemma 3.2 we get

∂ t 0 J(t 0 , x 0 , s 0 , ψ h * ) = -µ(s 0 )(s in -s 0 ) - t * t 0 ∂ s µ(s(t))(s in -s(t)) ṡ(t) dt = -µ(s 0 )(s in -s 0 ) - t * t 0 d dt µ(s(t))(s in -s(t)) dt = -µ(s * )(s in -s * ).
We conclude by recalling that h * = µ(s * )(s in -s * ) and therefore h * is a fixed point of h → -∂ t 0 J(t 0 , x 0 , s 0 , ψ h ).

Chapter 4

Spatially Heterogeneous Modelling of an Upflow Fixed-bed Bioreactor

This chapter corresponds to ongoing work.

Introduction

Anaerobic digestion (AD) is a commonly used process for the removal of organic matter in the treatment of wastewater and agricultural waste. Through a series of biochemical reactions, organic compounds are degraded by microorganisms into biogas, a mixture of methane and carbon dioxyde, which can be used as a renewable energy source, thereby reducing the energetic cost of the process.

The design and operation of AD bioreactors raises a number of challenges, since anaerobic digestion is known to be a complex, nonlinear and unstable process. The most common type of reactor in use for AD is the continuously stirred tank reactor (CSTR), which uses a mixing system to homogenize the reacting medium. Increasingly, reactors with heterogeneous contents are being developed and implemented. However, the benefits of a spatial gradient of concentration are unclear: while some experimental studies report better performance or process stability with mixing [START_REF] Borole | Methane production in a 100-l upflow bioreactor by anaerobic digestion of farm waste[END_REF][START_REF] Kalia | Effect of mixing digested slurry on the rate of biogas production from dairy manure in batch fermenter[END_REF], others have observed the opposite [START_REF] Gómez | Anaerobic codigestion of primary sludge and the fruit and vegetable fraction of the municipal solid wastes: Conditions for mixing and evaluation of the organic loading rate[END_REF][START_REF] Stroot | Anaerobic codigestion of municipal solid waste and biosolids under various mixing condi-tions-i. digester performance[END_REF] or that mixing has little impact on biogas production [START_REF] Karim | Anaerobic digestion of animal waste: Effect of mixing[END_REF].

Mathematical modeling of bioreactors has been recognized as an important tool for the analysis, control and optimization of the process [START_REF] Batstone | Mathematical modelling of anaerobic digestion processes: applications and future needs[END_REF] and in addition to giving a deeper insight into the process, modeling allows to evaluate different bioreactor designs and operational scenarios. Therefore, models have been developed that are capable of representing spatial variations of reactor contents in order to study the impact of heterogeneity. These differ from the classical models that focus on the biochemical kinetics and generally assume that the reactor is perfectly mixed.

The first examples of such models, known as compartment models or gradostat [START_REF] Smith | The theory of the chemostat: dynamics of microbial competition[END_REF], represent a reactor as a network of interconnected well-mixed zones and thus use systems of ordinary differential equations, which facilitates analysis and simulation. Studies with this type of models have found that the impact of heterogeneity depends on the graph of interconnections [START_REF] Haidar | Effects of spatial structure and diffusion on the performances of the chemostat[END_REF][START_REF] Rapaport | Some non-intuitive properties of simple extensions of the chemostat model[END_REF] and that mixing is either be beneficial [START_REF] Bensmann | Reactor configurations for biogas plants-a model based analysis[END_REF], disadvantageous [START_REF] Bello-Mendoza | Modelling the effects of imperfect mixing on the performance of anaerobic reactors for sewage sludge treatment[END_REF] or has little impact [START_REF] Donoso-Bravo | Modelling of an anaerobic plug-flow reactor. process analysis and evaluation approaches with non-ideal mixing considerations[END_REF].

The emergence of computer fluid dynamics (CFD) has allowed the simulation of models of AD bioreactors that take into account complex physical phenomenons and although many studies only consider fluid dynamics without bioreactions, this leads to complex models. This has been widely used to study and optimize various mixing systems for CSTR [START_REF] Leonzio | Study of mixing systems and geometric configurations for anaerobic digesters using cfd analysis[END_REF][START_REF] Terashima | Cfd simulation of mixing in anaerobic digesters[END_REF][START_REF] Vesvikar | Flow pattern visualization in a mimic anaerobic digester using cfd[END_REF][START_REF] Wu | Cfd simulation of mixing in egg-shaped anaerobic digesters[END_REF]. Recently however, simulations have been performed that consider bioreactions [START_REF] Meister | Integrating hydrodynamics and biokinetics in wastewater treatment modelling by using smoothed particle hydrodynamics[END_REF][START_REF] Rezavand | A fully lagrangian computational model for the integration of mixing and biochemical reactions in anaerobic digestion[END_REF] and one study [START_REF]Integration of mixing, heat transfer, and biochemical reaction kinetics in anaerobic methane fermentation[END_REF] has concluded that mixing has little impact on methane yield.

A third type of model attempts to find a compromise between model complexity and physical accuracy. These models represent a reactor in 1 or 2 spatial dimensions with simple fluid dynamics or even assuming that fluid velocity is constant in space to focus on the bioreactions [START_REF] Dochain | Dynamical modelling & estimation in wastewater treatment processes[END_REF][START_REF] Kung | The growth of pure and simple microbial competitors in a moving distributed medium[END_REF][START_REF] Mu | Anaerobic digestion model no. 1-based distributed parameter model of an anaerobic reactor: I. model development[END_REF][START_REF] Winkin | Dynamical analysis of distributed parameter tubular reactors[END_REF]. The advantage of reducing complexity is that it allows a more in depth analysis or the consideration of optimization problems [START_REF] Barbier | Modelling of biological decontamination of a water resource in natural environment and related feedback strategies[END_REF][START_REF] Crespo | Modeling and optimization of activated sludge bioreactors for wastewater treatment taking into account spatial inhomogeneities[END_REF][START_REF]Shape optimization of spatial chemostat models[END_REF][START_REF] Moya | Mathematical modeling and optimization of bioreactors and liquid crystals[END_REF].

The first objective of the present work is to continue the development of this type of model that offers a more accurate representation of spatial features than compartment models but that is still tractable for optimization of bioreactor design and operation. In addition, we are interested in investigating the impact of heterogeneity on organic matter removal and biogas production.

For this, we model here a pilot scale reactor operated in Guadalajara-Jalisco (Mexico) that is used to treat diluted tequila vinasses by anaerobic digestion and produce biogas. This bioreactor is particularly interesting for the study of heterogeneity as it is half way in between a CSTR and tubular unmixed reactor. Indeed, there is no internal mixing system and instead, homogenization is achieved primarily by recirculating liquid from the output back into the input and therefore, by changing the recirculation flow rate, this reactor can either be operated as a CSTR or as a tubular reactor. For the moment, experiments were run only with a high recirculation, but nonetheless, a spatial gradient was observed as data was gathered by collecting substrate at different points along the main axis of the reactor. This was used to develop a compartment model consisting of two interconnected homogeneous zones [START_REF] Hmissi | Evaluation of alkalinity spatial distribution in an up-flow fixed bed anaerobic digester[END_REF]. This model was able to reproduce the spatial gradient roughly and the present work aims at improving the representation of spatial features.

In Section 2, we present our model and the challenges of its numerical simulation. Section 3 deals with parameter estimation and model verification. Finally, in Section 4, we study the impact of mixing on biological activity.

Model

The bioreactor that we model here is a vertical cylinder with the liquid to be treated entering at the bottom and exiting at the top (Figure 4.1). A portion of the liquid from the top is recirculated and mixed with the influent before re-entering the tank at the bottom. Mixing is also helped by pumping out liquid from the very bottom of the tank which is injected back into the tank at the same height. Furthermore, both the inflows at the bottom are such that the liquid enters nearly tangentially in order to induce a circular current in the bottom of the tank. In the middle, there is a fixed bed made of vertical PVC tubes with honeycomb structures (cloisonyl), which helps to fix the biomass.

The configuration of this reactor and in particular the arrangement of the input and output flows, means that we need to consider a 3D model and can not use symmetry to reduce the spatial dimensions for the lower and upper parts of the reactor. However, the PVC tubes in the middle section of the reactor are narrow and once the fluid has entered one of these tubes from the bottom, it remains in the same tube until it reaches the top section. We can therefore consider that the contents inside each tube are homogeneous in the horizontal directions and only model the spatial distribution of substances along the vertical dimension. Furthermore, observations of this reactor have reported that most of the active biomass is attached to the walls of the tubes and we will thus neglect the suspended biomass, considering that the bio-reactions take place only in the middle section.

The result is a model of the reactor in 3 parts : bottom and top sections in 3D, where we must consider the fluid dynamics to compute the spatial distribution of the substrates; and the middle section, modelled as an array of parallel 1D tubular reactors, thereby reducing considerably the complexity of the model is this section.

Fluid Dynamics

The vinasses that are treated are completely liquid and do not contain any solid material, so that we can consider that the contents of the reactor are only in liquid or gas phases. We will model the liquid phase as an incompressible viscous fluid, having similar physical properties to water. It is important to note that the timescale of the dynamics of fluids is much shorter than the biological timescales and therefore we can consider that the liquid flow remains in a steady state, as long as the input-output flows are constant. We therefore use the steady state Navier Stokes equations to compute the fluid velocity and pressure. The gases produced by the biological activity collect at the top of the reactor and the main impact gases can have on reactor dynamics is through pressure. Indeed, the solubility of certain gases can vary with the pressure, such as CO 2 for instance. This effect could be taken into account by adding a boundary condition for the pressure at liquid-gas interface at the top of the reactor, with the pressure of the gas phase obtained from the ideal gas law. However, we are primarily interested in the biological activity and we can estimate biogas production through the methane generation. Indeed, due to the very low solubility of methane, its molar flow rate is unaffected by pressure. Then, since the pressure in the Navier-Stokes equations is defined only up to a constant, we will neglect the dynamics of gases and consider that the upper boundary of the liquid is fixed and similar to a wall.

We now detail the model equations for the fluid dynamics in the three sections of the reactor.

Bottom This part is a cylinder Ω B ⊂ R 3 , with boundary ∂Ω B = Γ in ∪ Γ BM ∪ Γ mi ∪ Γ Bmo ∪ Γ Bw .
The fluid enters on the side through Γ in , leaves towards the middle section through Γ BM and Γ Bw is the wall of the tank, through which the fluid cannot pass. For the mixing system, the fluid is pumped out through Γ mo and renters through Γ mi .

We denote the fluid velocity U (t, x, y, z) = (u x (t, x, y, z), u y (t, x, y, z), u z (t, x, y, z)), and p(t, x, y, z) the pressure and g the acceleration due to gravity. The steady state incompressible Navier-Stokes equations in 3 dimensions are

U • ∇U -ν∆U + ∇p = g, ∇ • U = 0. ( 4.1)
Concerning boundary conditions, for the walls of the bioreactor we consider a no-slip condition, U (t, x, y, z) = 0 on Γ Bw .

For the inflow and the mixing system, we consider classical Poiseuille flow profiles: let U Bin , U mi and U mo be unitary parabolic functions such that

Γ Bin U Bin dσ = Γ mi U mi dσ = Γmo U mo dσ = 1.
We denote the total input flow rate as

Q tot = Q r + Q in , Q in the influent
-3 m 2 s -1 D s Substrate diffusion coefficient 10 -7 m 2 s -1 Q in Influent flow rate 5 L h -1 Q r Recirculation flow rate 150 L h -1 Q m Mixing flow rate 1000 L h -1
direction of injection and the outwards unit normal n, the boundary conditions are then

U (t, x, y, z) = Q tot U Bin (x, y, z) cos α, U Bin (x, y, z) sin α, 0 on Γ in , U (t, x, y, z) = Q m U mi (x, y, z) cos α, U mi (x, y, z) sin α, 0 on Γ mi , U (t, x, y, z) = Q m U mo (x, y, z) n on Γ mo .
For the boundary Γ BM between the bottom and middle section, where the vertical inner tubes begin, we will neglect a small interface layer and consider that the flow is completely in the vertical direction and homogeneous in the horizontal directions. Denoting |Γ BM | the surface area of the interface, we have

U (t, x, y, z) = Q tot |Γ BM | 0, 0, 1 on Γ BM .
Middle This is the section of the bioreactor with the honeycomb structure : there are large tubes each divided in smaller tubes, for a total of N t ≈ 100 tubes. We model these as 1D vertical tubular bioreactors, indexed by a horizontal position (x i , y i ) ∈ Γ BM , i ∈ {1, ..., N t }, with each tube Ω (x i ,y

i ) = [z 0 , z 1 ] × {(x i , y i )}.
The fluid velocity in each tube Ω (x i ,y i ) is supposed constant in space, equal to the velocity at the bottom interface Γ BM , so that we have

u x (t, x, y, z) = 0, u y (t, x, y, z) = 0, u z (t, x, y, z) = Q tot |Γ BM | .
Top The upper section, Ω U is similar to the bottom although the fluid enters from the bottom through Γ U M , leaves through Γ U o and does not pass through the walls Γ U w . The fluid velocity satisfies the same Navier-Stokes equations (4.1) as in the bottom section but the boundary conditions are

U (t, x, y, z) = Q b |Γ BM | 0, 0, 1 on Γ U M U (t, x, y, z) = 0 on Γ U w U (t, x, y, z) = Q in U U o (x, y, z) on Γ U o
where U U o is again a unitary Poiseuille flow.

Biokinetics

For the anaerobic digestion process, we consider the two reaction model of [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF]. The first reaction, acidogenesis, represents the degradation of organic matter S 1 into volatile fatty acids (VFA) S 2 by acidogenic microorganisms B 1 . Then methanogenes B 2 transform the VFA into methane.

k 1 S 1 r 1 -→ B 1 + k 2 S 2 , k 3 S 2 r 2 -→ B 2 + k 4 CH 4
Here reaction rate r 1 , r 2 are the reaction rates and k 1 , k 2 , k 3 , k 4 are the yield coefficients.

In [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF] the reaction rate were taken as r k = µ k (S k )B k , for k = 1, 2, where the specific growth rates µ 1 , µ 2 are the Monod and Haldane growth functions

µ 1 (s) = µ max 1 s K s 1 + s , µ 2 (s) = µ max 2 s K s 2 + s + s 2 /K i .
where S r k is the substrate recirculation concentration coming from the upper section of the bioreactor which is computed as, denoting Γ U o the section of the boundary where the liquid leaves the reactor,

S r k (t) = 1 Γ U o Γ U o S k (t, x, y, z)dσ
For the auxiliary mixing system, we take similar boundary conditions

(S k U + D k ∇S k ) • n = S mix k U • n on Γ mi ∇S k • n = 0 on Γ mo with as before S mix k (t) = 1 Γ Bmo Γ Bmo S k (t, x, y, z)dσ
At the interface with the middle section Γ BM , we could consider flux equality conditions

(S k U + D k ∇S k ) • n Ω B = (S k U + D k ∇S k ) • n Ω M on Γ BM
where Ω M is the middle section. However, it is well know that for advection-diffusion equations, when the diffusion coefficient is small, this type of interface conditions can be well approximated by using the following artificial transparent boundary condition [START_REF] Halpern | Artificial boundary conditions for the linear advection diffusion equation[END_REF] ∇S k • n = 0 on Γ BM .

Middle The various concentrations in consideration here are functions of (t, z) (time and the vertical direction) and are also indexed by a horizontal position (x i , y i ) ∈ Γ BM and we denote them S k (t, z; x i , y i ) for the substrate with k = 1, 2 and similarly for the biomasses.

The substrates satisfy the following equation, in each Ω (x i ,y i )

∂ t S k + u z ∂ z S k -D k ∂ 2 zz S k = f k (ξ)
, where f k (ξ) corresponds to the biological reaction with ξ = (S 1 , S 2 , B 1 , B 2 ),

f 1 (ξ) = -k 1 µ 1 (S 1 )g 1 (B 1 )B 1 f 2 (ξ) = k 2 µ 1 (S 1 )g 1 (B 1 )B 1 -k 3 µ 2 (S 2 )g 2 (B 2 )C 2 .
For boundary conditions we take, at the bottom interface

S k (t, z 0 ; x i , y i ) = S k (t, x i , y i , z 0 ) Ω B ,
and for the top interface we take again an artificial transparent boundary

∂S k ∂z (t, z 1 ; x i , y i ) = 0.
We suppose that the micro-organisms are fixed to the wall of the tubes so that the biomass concentrations are neither advected nor diffused, but we must take into account their death rate

τ k ∂ t B k = µ k (S k )g k (B k )B k -τ k B k ,
Note that this is a family of ordinary differential equations, so that no boundary conditions are required.

For the methane flow rate, as proposed in [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF], we will consider that it is proportional to the growth rate of the second biomass. To compute the total flow rate for the whole reactor, we need to integrate over all the middle section Top We again suppose that there is no biomass in this section, so that the substrates satisfy the same equations as in the bottom. The boundary condtions are

Q CH 4 = |Γ BM | Nt i=1 z 1 z 0 k 4 µ 2 S 2 g 2 B 2 B 2 dz.
s k (t, x, y, z) Ω U = s k (t, z 1 ; x, y) Ω M on Γ U M ∇s k • n = 0 on Γ U o ∪ Γ U w

Numerical Resolution

To solve the model equation we use the Finite Element method implemented in the open source library FreeFem++ [START_REF] Hecht | New development in freefem++[END_REF]. For the fluid dynamics, we solve the time dependent Navier-Stokes equations until reaching a steady state. We consider the well established combination of a characteristics method, to deal with the non-linear convection term [START_REF] Pironneau | On the transport-diffusion algorithm and its applications to the navier-stokes equations[END_REF], and P2-P1 type finite elements [START_REF] Glowinski | Finite element methods for navier-stokes equations[END_REF]. As expected, the fluid converges to a steady state rapidly and Figure 4.2 illustrates the computed flow for standard operating conditions.

For the substrate equations of the bottom and top sections, we also use a finite element method with an implicit time discretization, which presents the advantage of being stable for larger time steps and coarser meshes than explicit time discretizations. However, to correctly solve these equations it is necessary to take time steps that are very small compared to the biological timescale. Indeed, especially in the bottom with high fluid velocities due to the mixing system, we need to take time steps of the order of seconds, whereas the typical timescale of a bioprocess is of the order of a day.

Another complication comes from the very small diffusion coefficients of the substrates. For example, a typical component of VFAs, acectic acid is reported to have a diffusion coefficient of the order of 10 -9 m 2 /s [START_REF] Winkelmann | Diffusion coefficient of acetic acid into water and ethanol solution[END_REF]. The consequence is that the substrate concentrations can present sharp spatial variations and to resolve these correctly, fine meshes must be used. Then, for smaller diffusion coefficients, the computational cost of simulations is substantially increased, as for example with D = 10 -8 it is necessary to use a mesh with 4 times more points than for D = 10 -7 (Table 4.3).

To investigate the effect of lowering the diffusion coefficient on the solution, we have run simulations of the bottom section for different values of D. Starting from an homogeneous initial condition, our test case considers a 10% increase in input concentration. The substrate concentrations at the interface between the bottom and middle sections after 8 hours are shown in Figures 4.3 (with mixing) and 4.4 (without mixing). Note that the main change of behaviour occurs between D = 10 -6 and D = 10 -7 , when the spatial variations of concentration become sharper and the homogenization effect of diffusion less important. Therefore we will use D = 10 -7 as we are essentially interested in the macroscopic behavior of the reactor and the qualitative effects of heterogeneity.

Parameter Estimation and Model Verification

The computational cost of simulations for this type of model makes it unpractical to fit many parameters by minimizing the difference between model outputs and experimental data. Instead, we will use parameters values from [START_REF] Hmissi | Evaluation of alkalinity spatial distribution in an up-flow fixed bed anaerobic digester[END_REF], that were obtained by fitting a simpler compartment ODE model to experimental data of the reactor that we study here. This simpler model was also based on the 2 reactions model of [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF] and it was shown that it is capable of reproducing the input-output behavior of the reactor, so that we will use the same growth function parameters and yield coefficients.

However, this work did not consider density dependent growth functions or biomass death rates and therefore we must estimate these parameters. For the death rates, it is generally accepted that they are small compared to the maximum growth rate and therefore we will take, somewhat arbitrarily, τ 1 = 0.1 • µ max For the density dependence parameters (c 1 , c 2 ), we will use the same experimental data used to fit the model of [START_REF] Hmissi | Evaluation of alkalinity spatial distribution in an up-flow fixed bed anaerobic digester[END_REF], which contains records of VFA concentration. There is however no data for S 1 and therefore to get a first estimate of c 1 , we use the simpler model of [START_REF] Hmissi | Evaluation of alkalinity spatial distribution in an up-flow fixed bed anaerobic digester[END_REF]. Indeed, this model takes into account the physico-chemical dynamics of the reactor (including pH, alkalinity, strong ions, inorganic carbon and carbon dioxyde) for which experimental data was available and since these variables depend on S 1 , we can assume that the model represents correctly the input-output behavior of S 1 , especially when the reactor is in a steady state. Therefore, to estimate c 1 , we have run simulations for various values of this parameter and we then take the value for which our model converges to the same steady state as the model of [START_REF] Hmissi | Evaluation of alkalinity spatial distribution in an up-flow fixed bed anaerobic digester[END_REF]. Figure 4.5 illustrates the best fit. Now, using the experimental data, we can estimate c 2 , the parameter of the second growth function. The VFA concentration was measured from the output but also at several points within the reactor, at heights of 30 and 90 cm from the bottom. This will allow us to check how well our model can represent the spatial gradient of concentrations within the reactor. Figure 4.6 illustrates the best fit, along side the output of the simpler model. Notice that the 3D model captures well the position where most of the degradation takes place and is also capable of reproducing the vertical variation better than the ODE model.

Impact of Recirculation and Mixing

The value of this model is that it can reproduce the spatial distribution of the contents of the reactor and thus can help to understand the internal dynamics. In particular, the detailed representation of the physical processes of advection and diffusion makes this model a valuable tool to study the impact of physical operational parameters, such as the recirculation and mixing flowrates. In this section we explore this question with numerical simulations of tracer experiments and of the full model.

Tracer Pulse

The first simulations that we propose are numerical versions of tracer pulse experiments, which is a common method of characterizing and understanding the hydrodynamics of reactors. These consist in releasing a short pulse of an inert substance, the tracer, in the inflow and measuring the concentration at the output. For our simulations, we will therefore run the model with only one substrate and no bio-reactions. Starting from a homogeneous initial condition, we consider a short pulse of high influent concentration, with S in = 60 g/L for 1 min. Figure 4.7 shows the tracer concentration in the reactor during the first 5 hours. The first observation that we can make is that the bottom mixing system does not induce a circular current in the whole of the bottom and instead, mixes the tracer only in a small section. This can be attributed to the fact that the input and output of the mixing system are close to each other and therefore the fluid exiting this system is pulled back towards it directly. This suggests that the complete homogenization of the bottom is achieved not only through advection but that diffusion is also important. This simulation also shows that, although a portion of the tracer goes directly up and out, an important part of the tracer is pulled in by the mixing system and is delayed before going up. This is confirmed by comparing the output concentrations with and without mixing, i.e. for Q m = 1000 L/h and Q m = 0, as show in Figure 4.8. Indeed, the initial spike and the total mass exiting during the first days (proportional to the integrals of the curves) is much higher without mixing. This indicates that with mixing an important part of the tracer remains trapped in the bottom of the reactor and therefore the mixing system also acts as a delay. Now, to study the impact of the recirculation flow rate, we have run tracer simulations for Q r ∈ {50, 100, 150, 200} and we show the tracer output concentration in Figure 4.9. Notice that for Q r = 100, the initial spike is greatly diminished and seems to disappears completely for Q r = 50. In fact, for Q r = 50, the maximum output concetration is not reached at first, which indicates that the pulse of tracer has mainly gone through the mixing pump. More generally, the mixing system becomes proportionally stronger as the ratio Q m /Q r becomes higher. The result is that for a given mixing flowrate, better homogenization is achieved with smaller recirculation flow rates, in the sense that perturbations of the input, like the pulse of tracer, are absorbed by the mixing system.

Nonetheless, in all cases, after the initial spike the reactor behaves similarly to a well mixed chemostat, with the output concentration slowly decreasing. Indeed, on longer timescales, the diffusion process helps with the homogenization of the reactor, for all mixing and recirculation flowrates.

Impact on Biological Activity

Unlike the inert tracer, the substrates will not always tend towards an even distribution, since they are consumed through the bio-reactions. Therefore, we now look at the impact of the flowrates on the biological activity and we are primarily interested in reactor performance at steady state.

To study the influence of the recirculation flowrate Q r , we have run simulations for a range of flowrates from 50 to 200 L/h and we plot the steady state values for the output substrate concentration S 1 and S 2 and the methane flowrate in Figure 4.10. It appears that a smaller recirculation flowrate gives better results, in particular in terms of output concentration and for Q r = 50 L/h there is 7% less S 1 and 12% less S 2 than for Q r = 150L/h. For the biogas production, the difference is a lot smaller, with Q r = 50 L/h there is only 1.7 % more methane compared to standard operating conditions. 4.12 show the distributions of VFA and methanogenic biomass for Q r = 50 and 150 L/h. Observe that for the slower recirculation rate, the variations of substrate and biomass concentrations are a lot sharper and the reaction takes place mainly at the beginning of the tubes. This can be explained by the fact that, because of the lower recirculation rate, the influent concentration S if k (k = 1, 2) is less diluted as the substrate concentration entering the bottom is

S in k = Q r S r k + Q in S if k Q r + Q in .
Then, a higher substrate concentration arrives in the tubes, which translates into a higher growth rate. This leads to more biomass, which also increases the growth rate and this turns into a positive feedback, which is only balanced by the density dependence of the growth rate. The result is that the substrate is rapidly consumed until the concentration drop to levels at which the growth rate is very low and therefore little reaction occurs in the upper part of the tubes. This phenomenon is the reason why the biomass density dependence and substrate inhibition of the growth rate must be correctly taken into account. In our case, these results must be taken with caution, as the growth rate parameters where obtained with data that had low concentrations of substrate and thus the obtained values are not very accurate for high substrate concentrations. Notice in particular that the substrate inhibition constant K i of the second growth rate is very high, so that this function is very close to a Monod function [START_REF] Khedim | Effect of control parameters on biogas production during the anaerobic digestion of proteinrich substrates[END_REF] and therefore a lower value should be used.

Conclusion

In this work, a 3D model of an Upflow Fixed Bed bioreactor was developed and it was shown that it is capable of reproducing experimental data and in particular, the vertical variations of substrate concentration. Then, with numerical simulations, the impact of heterogeneity on biological activity was studied. A number of conclusions can be drawn from this work.

• First, we have shown the value of simpler models in terms of physical representation: without turbulence or multi-phase flow, this model can reproduce experimental data and help to better understand internal dynamics of the reactor.

• Despite the reduced complexity compared to other 3 dimensional models, the difference of timescales between the physical and biological process mean that simulations are computationally intensive. Serious work is needed to optimize the numerical method and implementation.

• The efficiency of the bottom mixing system is governed by the ratio Q m /Q r and its influence is mainly seen during the transients, by delaying or absorbing and smoothing out perturbations. The closeness of the input and output of this mixing system is detrimental and the system might be more efficient if they were further apart. • The bioreactor could be more efficient with a lower recirculation flow rate. However, this leads to higher substrate concentration at bottom which could destabilize reactor. Further experimental data is necessary to confirm this result and in order to have adequate representation of inhibition due to high substrate and biomass levels.

Chapter 5

Conclusions

The purpose of this thesis was the study of optimal control problems and mathematical models for the maximisation of biogas production. We present here the main contribution of this thesis and the perspectives of future work.

• Chapter 2

-We presented a method to obtain an estimation of the value function of the problem of maximization of biogas production in a chemostat on a finite horizon, thereby establishing sub-optimality bounds for controls used in practice, for a large class of substrate and biomass dependent growth functions. For the infinite horizon, we showed that the value function of the discounted problem converge when the discount rate goes to 0 and that the limit is equal to the value function of the average reward.

-The time horizon influences the choice of the best MRAP type feedback. For short horizons it is best to drive the chemostat to a substrate level that maximize the biogas flow rate on the set { x s in -s = x 0 s in -s 0 }. On the other hand, as the horizon becomes longer, it is best to drive the system to a maximizer on the set {x + s = s in }, and this becomes optimal, on average, for an infinite horizon.

-The methods developed could be extended to more complex models. For example, it is straightforward to show that for the two reaction model, there is a 2 dimensional invariant set that is attractive for persistently exciting controls. Then, it could be shown that any control driving the system to a maximizer on the invariant set is optimal for the average reward on an infinite horizon, although the existence of such a control is not guaranteed. On the other hand, with appropriate conditions on the growth function, we could establish an estimation of the value function and we could find sub-optimal control by studying the problem for initial conditions on the invariant set.

• Chapter 3

-We presented a method to transform the Hamilton-Jacobi-Bellman (HJB) equation into a fixed point equation by writing extremal controls as functions of the Hamiltonian. This allowed us to developed an algorithm to identify the singular arc of a candidate to optimal control. The computed control is obtained in feedback form and depends on the initial condition and time horizon considered.

-This method could be extended for more general problems to identify singular arcs. For example, for systems with dynamics that are affine in the control, it should be possible to write a expression of extremal controls that depends on the Hamiltonian and then our algorithm could be used to compute the singular arc.

• Chapter 4

-We have introduced a spatially heterogeneous model of a real bioreactor, capable of reproducing the spatial gradient of substrate concentration. It helps to better understand the internal dynamics of the reactor and shows that it is not always important to consider overly complex fluid dynamics.

-Our preliminary results indicate that a heterogeneous reactor is more efficient than a well mixed device. However, we must accurately take into account inhibition phenomenons to confirm our results. For example, we could consider a different density dependence term in the growth rate, such as

g(B) = 1 1 + cB 2
in order to have a reaction rate that goes to 0, when biomass concentrations become very large.

-The difference of timescales between the physical and biological process make it impossible to consider complex optimization problems due to the high computational cost of simulations. However, simulations have shown that although there are variations at reactor scale, concentrations are 'locally' homogeneous. This could be accurately approximated by a compartment model, although more complex than the 2 zone model of [START_REF] Hmissi | Evaluation of alkalinity spatial distribution in an up-flow fixed bed anaerobic digester[END_REF]. Such a model would be much faster to simulate and realistic optimization problems could be considered. Our model could be used to determine the configuration of the interconnections of the compartment model by minimizing the difference between both models.

Figure 1 . 1 :

 11 Figure 1.1: Examples of Monod and Haldane growth functions

Figure 1 . 2 :

 12 Figure 1.2: State space trajectories with feedback ψ s. The black line represents the invariant set {x + s = s in }. Haldane growth function (μ = 0.74, K s = 9.28, K i = 256) with s in = 100, D max = 3.

Figure 1 . 3 :Figure 1 . 4 :

 1314 Figure 1.3: Auxiliary value function (x 0 , s 0 ) → W z 1 (0, x 0 , s 0 ) with z 1 = 1. On the left, Contois growth function (µ max = 0.74, K s = 1, u max = 1.5) and on the right, Haldane growth function (μ = 0.74, K s = 9.28, K i = 256, u max = 3). In both cases, s in = 100 and T = 2.

Figure 1 . 5 :

 15 Figure 1.5: Example of trajectories obtained with feedback ψ h with h = 2 and Monod growth function (µ max = 1.2, K = 7.1), D max = 0.7 and s in = 10.

Lemma 1 . 9 ([ 3 ,

 193 Theorem III.3.42]) A measurable function D : [t 0 , T ] → [0, D max ] maximizes (P bio ), the production of biogas problem, if and only if the maximum condition of the PMP holds and

Figure 1 . 6 :

 16 Figure 1.6: Optimal trajectories in state space, for the Monod growth function (µ max = 1.2, K = 7.1) with t 0 = 0, T = 2.5, D max = 0.7 and s in = 10. The initial conditions are on the left (x 0 , s 0 ) = {(1, 3), (5, 1), (11, 2), (15.5, 4.5)} and on the right, (x 0 , s 0 ) = {(1, 5.5), (2, 6.5), (5, 7), (10, 7)}.

Figure 1 . 7 :

 17 Figure 1.7: Schematic view of the reactor

Figure 1 . 8 :Figure 1 . 9 :

 1819 Figure 1.8: Tracer concentration for standard conditions, at 1, 2, 3, 4 and 5 hours after pulse.

Figure 2 . 1 :

 21 Figure 2.1: State space trajectories with feedback ψ s. The black line represents the invariant set {(x, s) : x + s = s in }. Haldane growth function (μ = 0.74, K s = 9.28, K i = 256) with s in = 100, u max = 3.

Figure 2 . 2 :

 22 Figure 2.2: State space trajectories with feedback ψ s(z 0 ) for z 0 ∈ {0.2, 0.7, 1.5, 3} and s 0 ∈ {10, 60, 75}. The color and type of line indicates the value of z 0 . Contois growth function (µ max = 0.74, K s = 1) with s in = 100, u max = 1.5.

Figure 2 . 3 :

 23 Figure 2.3: Normalized average reward J N (T, z 1 ) as a function of z 1 ∈ [z 0 , 1] and T ∈ [0.5, 6] for the initial condition (x 0 , s 0 ) = (20, 20). Contois growth function (µ max = 0.74, K s = 1) with s in = 100, u max = 1.5.

Figure 2 .

 2 Figure 2.3 shows a case when z 0 < 1 and Figure 2.4 is an example of z 0 > 1.We can see clearly that for small final times, the maximum is attained for a value of z 1 close to z 0 and that for z 1 = 1 the reward is the smallest. However, as the final time increases, the value of z 1 for which the reward is maximum approaches 1, and with the feedback ψ s(z 0 ) the reward is the smallest. In particular, we can see that the best of the feedbacks ψ s(z 1 ) depends on the final time.

Figure 2 . 4 :

 24 Figure 2.4: Normalized average reward J N (T, z 1 ) as a function of z 1 ∈ [1, z 0 ] and T ∈ [0.5, 6] for the initial condition (x 0 , s 0 ) = (70, 60). Contois growth function (µ max = 0.74, K s = 1) with s in = 100, u max = 1.5.

Figure 2 . 5 :

 25 Figure 2.5: On the left, t → x(t) (solid lines) and t → s(t) (dashed lines) with feedbacks ψ ∈ {ψ s(z 0 ) , ψ s(1) , ψ S } and on the right, the corresponding open loop controls. Contois growth function (µ max = 0.74, K s = 1) with s in = 100, u max = 1.5 and initial condition (x 0 , s 0 ) = (30, 2).

Figure 2 . 6 :

 26 Figure 2.6: On the left, average reward as function of final time T→ J T (ξ, ψ(•)) with feedback ψ ∈ {ψ s(z 0 ) , ψ s(z 1 ) , ψ s(1), ψ S } with z 0 = 0.25 and z 1 = 0.625. On the right, the corresponding state space trajectories. Contois growth function (µ max = 0.74, K s = 1) with s in = 100, u max = 1.5 and initial condition (x 0 , s 0 ) =[START_REF] Clarke | Functional analysis, calculus of variations and optimal control[END_REF][START_REF] Clarke | Functional analysis, calculus of variations and optimal control[END_REF].

Figure 2 . 7 :

 27 Figure 2.7: On the left, average reward as function of final time T→ J T (ξ, ψ(•)) with feedback ψ ∈ {ψ s(z 0 ) , ψ s(z 1 ) , ψ s(1), ψ S } with z 0 = 1/3 and z 1 = 2/3. On the right, the corresponding state space trajectories. Contois growth function (µ max = 0.74, K s = 1) with s in = 100, u max = 1.5 and initial condition (x 0 , s 0 ) =[START_REF] Bello-Mendoza | Modelling the effects of imperfect mixing on the performance of anaerobic reactors for sewage sludge treatment[END_REF][START_REF] Meister | Integrating hydrodynamics and biokinetics in wastewater treatment modelling by using smoothed particle hydrodynamics[END_REF].

Figure 2 . 9 :

 29 Figure 2.9: Auxiliary value function (x 0 , s 0 ) → W z 1 (0, x 0 , s 0 ) with z 1 = 1. On the left, Contois growth function (µ max = 0.74, K s = 1, u max = 1.5) and on the right, Haldane growth function (μ = 0.74, K s = 9.28, K i = 256, u max = 3). In both cases, s in = 100 and T = 2.

. 7 )

 7 This, along with the transversality condition, yields c = µ(s(T ))x(T ) > 0. Let us call extremal trajectory and extremal control to any trajectory (x(•), s(•), p x (•), p s (•)) and control D(•) satisfying (3.1)-(3.5)-(3.6)-(3.7).

Figure 3 . 1 :

 31 Figure 3.1: Example of trajectories obtained with feedback ψ h with h = 2 and Monod growth function (µ max = 1.2, K = 7.1), D max = 0.7 and s in = 10. The thicker lines correspond to the singular curve G h .

Lemma 3 . 1 ([ 3 ,

 313 Theorem III.3.42]) Under the Standing Assumptions, a measurable functions D : [t 0 , T ] → [0, D max ] maximizes (P bio ), the production of biogas problem, if and only if the maximum condition (3.6) holds and

Figure 3 . 2 :

 32 Figure 3.2: Optimal trajectories in state space, for the Monod growth function (µ max = 1.2, K = 7.1) with t 0 = 0, T = 2.5, D max = 0.7 and s in = 10. The initial conditions are on the left (x 0 , s 0 ) = {(1, 3), (5, 1), (11, 2), (15.5, 4.5)} and on the right, (x 0 , s 0 ) = {(1, 5.5), (2, 6.5), (5, 7), (10, 7)}.

Figure 3 . 3 : 75 Figure 3 . 4 :

 337534 Figure 3.3: Optimal trajectories in state space, for the Haldane growth function (µ max = 0.74, K = 9.28, K i = 256) with t 0 = 0, T = 2, D max = 1 and s in = 100. The initial conditions are on the left (x 0 , s 0 ) = {(5, 10), (40, 5), (110, 1), (230, 23)} and on the right, (x 0 , s 0 ) = {(15, 40), (20, 60), (50, 70), (140, 70)}.

Figure 3 . 5 : 5 Figure 3 . 6 :

 35536 Figure 3.5: For the initial condition (x 0 , s 0 ) = (1, 4.5), on the left, optimal trajectories in state space and on the right, t 0 → c(t 0 ) value of the Hamiltonian as a function of initial time. Monod growth function (µ max = 1.2, K = 7.1) with T = 2.5, D max = 0.7 and s in = 10.

  and 3.5.

75 Figure 3 . 7 :

 7537 Figure 3.7: Error as a function of iterations : n → |h k n+1 -h k n | of Algorithm 1 on the left and of Algorithm 2 on the right. Monod growth function (µ max = 1.2, K = 7.1) with T = 1, D max = 0.7 and s in = 10.On the left initial condition (x 0 , s 0 ) = (2, 2) and on the left (x 0 , s 0 ) =[START_REF] Bayen | Analysis of an optimal control problem related to the anaerobic digestion process[END_REF][START_REF] Bayen | On the steady state optimization of the biogas production in a two-stage anaerobic digestion model[END_REF]

Figure 4 . 1 :

 41 Figure 4.1: Schematic view of the reactor

  flow rate, Q r the re-circulation flow rate and Q m the mixing flow rate. If α is the angle between the

Figure 4 . 2 :

 42 Figure 4.2: Streamlines of fluid flow, for standard operating conditions Q r = 150 L/h and Q mix = 1000 L/h. Bottom section (left) and top section (right).

Figure 4 . 3 :

 43 Figure 4.3: Bottom-Middle interface at 8h after a 10 % change of inflow concentration with mixing flow rate Q mix = 1000 L/h and different diffusion parameters. Left D = 10 -6 , Center D = 10 -7 , Right D = 10 -8 .

Figure 4 . 4 :

 44 Figure 4.4: Bottom-Middle interface at 8h after a 10 % change of inflow concentration with mixing flow rate Q mix = 0 L/h and different diffusion parameters. Left D = 10 -6 , Center D = 10 -7 , Right D = 10 -8 .

Figure 4 . 5 :

 45 Figure 4.5: Comparison with ODE model

1 and τ 2

 12 = 0.1 • µ max 2 .

9 S2

 9 [mmol/L] simulated 3D 30cm simulated 3D 90cm simulated 3D out VFA measured (30cm) VFA measured (90cm) VFA measured out btm (ODE) top (ODE)

Figure 4 . 6 :

 46 Figure 4.6: VFA (S 2 ) concentrations at different heights within the reactor (30 and 90 cm from the bottom) and at output, for the 3D model and experimental data. Values of the 2 zone ODE model of [56] are also shown.

Figure 4 . 7 :

 47 Figure 4.7: Tracer concentration for standard conditions Q r = 150 L/h and Q m =1000 L/h, at 1, 2, 3, 4 and 5 hours after pulse.

Figure 4 . 8 :

 48 Figure 4.8: Tracer output concentration with and without bottom mixing and Q r =150 L/h

Figures 4. 11

 11 Figures 4.11 and 4.12 show the distributions of VFA and methanogenic biomass for Q r = 50 and 150 L/h. Observe that for the slower recirculation rate, the variations of substrate and biomass concentrations are a lot sharper and the reaction takes place mainly at the beginning of the tubes. This can be explained by the fact that, because of the lower recirculation rate, the influent concentration S if k (k = 1, 2) is less diluted as the substrate concentration entering the bottom is

  Figures 4.11 and 4.12 show the distributions of VFA and methanogenic biomass for Q r = 50 and 150 L/h. Observe that for the slower recirculation rate, the variations of substrate and biomass concentrations are a lot sharper and the reaction takes place mainly at the beginning of the tubes. This can be explained by the fact that, because of the lower recirculation rate, the influent concentration S if k (k = 1, 2) is less diluted as the substrate concentration entering the bottom is

Q r = 50 QFigure 4 . 9 :

 5049 Figure 4.9: Tracer output concentration for various recirculation flowrates and Q m =1000 L/h

CH 4 [

 4 mmol/h]

Figure 4 . 10 :

 410 Figure 4.10: Steady state values for the output substrate concentration S 1 and S 2 and the methane flowrate for a range of recirculation flowrates Q r ∈ [50, 200] L/h.

Figure 4 . 11 :

 411 Figure 4.11: Distribution of S 2 and B 2 at steady state with Q r = 50 L/h.

Figure 4 . 12 :

 412 Figure 4.12: Distribution of S 2 and B 2 at steady state with Q r = 150 L/h.

Table 3 .

 3

		1: Performance comparaison with Bocop
	(x 0 , s 0 ) Biogas (ψ c ) Biogas (Bocop) Relative Difference
	(3, 2)	3.2232	3.2235	9 • 10 -5
	(3, 6)	5.3285	5.3290	9 • 10 -5
	(1, 4.5)	1.8725	1.8729	2 • 10 -4
	(7, 8)	4.7904	4.7933	6 • 10 -4

Table 3 .

 3 2: Computation time (in seconds) comparison with Bocop (x 0 , s 0 ) CPU time (ψ c ) ε = 10 -4 CPU time (ψ c ) ε = 10 -6 CPU time (Bocop)

	(3, 2)	2.32	3.60	2.06
	(3, 6)	1.37	1.74	1.74

  Lemma 3.2 For any initial condition (x 0 , s 0 ) ∈ I, for the trajectories with control ψ h * we have∂ t 0 s(t) = -f (s(t)) at time t ∈ [t 0 , t * ].

									t
	Proof. We can write the differential equation satisfied by s(•) as s(t) = s 0 +	f (s(τ )) dτ
							t		t 0
	and differentiating we get ∂ t 0 s(t) = -f (s 0 ) +	f (s(τ ))∂ t 0 s(τ ) dτ. This is a linear dif-
							t 0	t
	ferential equation and the solution is ∂ t 0 s(t) = -f (s 0 ) exp	f (s(τ )) dτ . Now, as
	t 0 f (s(t)) does not change sign for t ∈ [t 0 , t * ) and since f (s(t)) is the derivative of s(t) we
	have	t t 0	f (s(τ )) dτ =	s(t) s 0	f (s) f (s)	ds = ln	f (s(t)) f (s 0 )	.

Table 4 .

 4

		1: Physical parameters
	Parameter		Value Unit
	g	Acceleration due to gravity	9.81 m 2 s -1
	ν	Fluid viscosity	10

Table 4 .

 4

	3: Bottom section mesh parameters
	Diffusion Coefficient [m 2 /s]	Number of Points
	10 -7	61 882
	10 -8	244 860

The approach we have taken provides a new proof for the optimality of the synthesis already known for the reduced model, that is, the case where s in = x 0 + s 0 . To show consistency of our approach, the details for this case have been included in the appendix.

Acknowledgements

First and foremost, I would like to acknowledge my advisors, Héctor Ramírez and Alain Rapaport, for their guidance, inspiration, and support. I have learned a great deal from you about research, mathematics and otherwise. I thank you for the diversity of projects and opportunities that you offered me and the patience that you maintained while I tried to navigate my way through it all.

I have greatly benefited from working with Cristopher Hermosilla, Antoine Rousseau and Jérôme Harmand, I deeply enjoyed the time spent working with you and I thank you for many inspirational discussions.

Acknowledgments

Powered@NLHPC: This research/thesis was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02)

I would also like to thank my fellow students, who have helped make my time at university such a pleasure, in Chile and France.

Last, but certainly not least, I would like to thank my family, for your love and encouragement, and for supporting me throughout my studies and my life in general. This thesis was made possible thanks to the financial support from CONICYT-Chile through the doctoral fellowship CONICYT-PFCHA / Doctorado Nacional / 2017-21170249, FONDECYT Grants 1160567 and 1160204, grant REDES 150011 and Basal Program CMM-AFB 170001 from CONICYT-Chile.

difference approximation, we must work with a discrete range of initial times. For N ∈ N, we denote {t k 0 } k=1,..,N a set of initial times with constant step ∆t 0 = t k+1 0 -t k 0 . We will therefore compute a vector of fixed points c = (c k ) ∈ R N , where each c k will correspond to the value of the Hamiltonian for the initial time t k 0 .

To obtain h n+1 from the previous iterate h n = (h k n ) we start by computing for each t k 0 the trajectories with the control ψ h k n and the associated cost J(t k 0 , x 0 , s 0 , ψ h k n ) with standard numerical integration tools. For this the singular curve G h k n must be first determined by solving equations (3.12) to establish the admissible range [s max , s 0 ] and then integrating backwards to obtain x(•, x h (s max ), D max ), s(•, s max , D max ) .

We can then approximate the partial derivative of the cost as

and if we set t N 0 = T , we can use that J(T, x 0 , s 0 , ψ h ) = 0 to start the computations of these partial derivatives, running through the range of initial times backwards beginning with t N -1 0 and ending with t 1 0 .

In summary, for a fixed initial condition (x 0 , s 0 ) ∈]0, +∞[×]0, s in [ and final time T , the algorithm is shown below.

The main issue that can prevent the convergence of this algorithm is the accumulation of numerical errors that can propagate through the finite difference approximation of ∂ t 0 J(t 0 , x 0 , s 0 , ψ h ). Indeed, since we need J(t k+1 0 , x 0 , s 0 , ψ h k+1 n

) to compute h k n+1 , any errors made to get h k+1 n will propagate to h k n+1 and all following values h j n+1 for j < k. Another consequence of this inter-dependance is that to have h k n converge, h k+1 n must have already converged to a fixed point.

With these considerations in mind, it might seem unnecessary to compute the whole vector h n+1 at every iteration and instead computing one fixed point at a time would appear to be more efficient. An alternative algorithm would then consist in first iterating only on h N -1 n until convergence, which is possible because we only need J(T, x 0 , s 0 , ψ h ). Then using the obtained fixed point to get J(t N -1 0 , x 0 , s 0 , ψ c K-1 ) we could move on to computing c K-2 . Repeating this process, we can thus find all the fixed points c k until reaching the desired initial time. However, due to the accumulation of errors, to get the convergence of h j n with a certain tolerance it is necessary to get the convergence of h k n for j < k with a smaller tolerance. It is then complicated in practice to determine an efficient However in this work, in order to take into account the effects of crowding and the fact that a majority of the biomass is fixed to the walls of the tubes in the middle section, we will consider density dependent growth rates r k = µ k (S k )g k (B k )B k . Since we model the tubes in 1 dimension, we consider that every point of these tubes there is a 2D colony of microorganisms attached to the wall of the tube with the substrate passing through the center. Biomass inhibition functions have been proposed, for a 2D dimensional colony of microorganisms in [START_REF] Lobry | A new hypothesis to explain the coexistence of n species in the presence of a single resource[END_REF] and we will consider the following function

We now detail the model equation for the substrate and biomass concentrations in the bottom, middle and top sections.

Bottom We suppose that there is no biomass in this section of the reactor and therefore no reaction. The evolution of the substrate concentration S k (t, x, y, z) for k = 1, 2 is modeled by the following advection-diffusion equation

where D k is the diffusivity of the substrate S k . The boundary conditions for the inflow and the walls are

where n is the outwards unit normal and S in k is the substrate inflow concentration, for k = 1, 2