# An Algorithm for Maximizing the Biogas Production in a Chemostat

Antoine Haddon, Cristopher Hermosilla

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

## To cite this version:

Antoine Haddon, Cristopher Hermosilla. An Algorithm for Maximizing the Biogas Production in a Chemostat. 2019, 182 (3), pp.1150-1170. 10.1007/s10957-019-01522-x. hal-02617851

HAL Id: hal-02617851

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

Submitted on 13 Jun 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# An Algorithm for Maximizing the Biogas Production in a Chemostat

Antoine Haddon¹,² Cristopher Hermosilla³

¹ Mathematical Engineering Department and Center for Mathematical Modelling (CNRS UMI 2807), Universidad de Chile, Santiago, Chile (ahaddon@dim.uchile.cl)

² MISTEA, Université Montpellier, INRA, Montpellier SupAgro, Montpellier, France

³ Departamento de Matemática, Universidad Técnica Federico Santa María, Valparaíso, Chile

## Abstract

In this work, we deal with the optimal control problem of maximizing biogas production in a chemostat. The dilution rate is the controlled variable and we study the problem over a fixed finite horizon, for positive initial conditions. We consider the single reaction model and work with a broad class of growth rate functions. With the Pontryagin Maximum Principle, we construct a one parameter family of extremal controls of type bang-singular arc. The parameter of these extremal controls is 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.

Keywords: Optimal Control; Chemostat Model; Pontryagin Maximum Principle; Hamilton-Jacobi-Bellman Equations; Optimal Synthesis

## 1 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 [1]. 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 [2].

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 [3], 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 [4]. 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 [5], expert systems [6], fuzzy logic [7] and adaptive control [8], 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 [9, 10]. In addition, only problems on well-mixed continuously stirred tanks have been studied, since the non-linearities 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 [11]. 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 [12] and later solved for a special set of initial conditions for which the model reduces to a one dimensional problem [13]. More recently, the general one reaction model has been revisited to propose a sub-optimal control for which there is an estimation of sub-optimality [14]. Let us mention that the problem has also been considered in the infinite horizon case [15]. 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 [16, 17, 18, 19]. 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.

## 2 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 [20] :

$$ \dot{x} = \mu(s)x - Dx, \quad \dot{s} = D(s_{in} - s) - \mu(s)x, \quad (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 $\mu(\cdot)$ is usually chosen of Monod type ($\mu_M$) or of Haldane type ($\mu_H$): given $\mu_{max}$, $K$, $K_i > 0$

$$ \mu_M(s) := \mu_{max} \frac{s}{K+s}, \quad \mu_H(s) := \mu_{max} \frac{s}{K+s+s^2/K_i}. \qquad (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 $\mu(\cdot)$ is a twice continuously differentiable function on $[0, +\infty[$ such that

$$ \mu(0) = 0, \quad \mu(s) > 0 \quad \text{and} \quad \frac{d}{ds} \left( \frac{\mu'(s)}{\mu(s)^2} \right) \neq 0, \quad \forall s > 0. $$

In the Monod and Haldane cases we have that

$$ \frac{\mu'_{M}(s)}{\mu_{M}(s)^2} = \frac{K}{\mu_{max}s^2} \quad \text{and} \quad \frac{\mu'_{H}(s)}{\mu_{H}(s)^2} = \frac{K}{\mu_{max}s^2} - \frac{1}{\mu_{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 [21] and therefore the biogas production during a time interval $[t_0, T]$ for a given substrate concentration $s(\cdot)$ and a given microbial population concentration $x(\cdot)$ is

$$ \int_{t_0}^{T} \mu(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 \mapsto D(t)$ of the bioreactor under the constraint that $D(t) \in [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

$$ \left\{ \begin{array}{ll} \text{Maximize} & \displaystyle\int_{t_0}^{T} \mu(s(t))x(t)dt \\ \text{over all} & D: [t_0, T] \to [0, D_{\text{max}}] \text{ measurable} \\ \text{such that} & \dot{x} = \mu(s)x - Dx, \quad x(t_0) = x_0, \\ & \dot{s} = D(s_{\text{in}} - s) - \mu(s)x, \quad s(t_0) = s_0, \\ & 0 \le s(t) \le s_{\text{in}} \text{ and } 0 \le x(t), \quad \forall t \in [t_0, T]. \end{array} \right. \qquad (\text{P}_{\text{bio}}) $$

With a slight abuse of notation, we may sometimes write

$$ J(t_0, x_0, s_0, \psi) := \int_{t_0}^{T} \mu(s(t))x(t)dt $$

for the biogas production associated with a control $\psi: \mathbb{R}^2 \to [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(\cdot)$ and $x(\cdot)$ denote a solution of control system (1) in closed-loop form associated with these data.

**Remark 1.** Since feedback controls $\psi: \mathbb{R}^2 \to [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:

$$ \dot{x} = \mu(s)x - \psi(s, x)x, \quad \dot{s} = \psi(s, x)(s_{\text{in}} - s) - \mu(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 (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 [22].
## 2.1 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, +\infty] \times [0, s_{in}]$. In the rest of the paper this restriction will be disregarded. The main reason for doing so is that the system (1) is invariant on $K$ (see for example [23, Theorem 4.3.8]). Indeed, the set-valued map

$$F(x, s) := \{(\mu(s)x - Dx, D(s_{in} - s) - \mu(s)x) : D \in [0, D_{max}]\}$$

is locally Lipschitz continuous, has linear growth, has nonempty compact and convex images and satisfies the invariance condition $F(x, s) \subseteq T_K(x, s)$ for any $(x, s) \in K$; where $T_K$ stands for the Contingent Cone. The last affirmation comes from the fact that $F(0, s) \subseteq \{0\} \times [0, +\infty[$ for any $s \in [0, s_{in}]$ and

$$F(x, 0) \subseteq] -\infty, 0] \times [0, +\infty[, \quad F(x, s_{in}) \subseteq \mathbb{R} \times] -\infty, 0[, \quad \forall x > 0.$$

By similar arguments we can see that the set $\{0\} \times [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 \in ]0, s_{in}]$ will reach that set. Moreover, the fact that $F(x, s_{in})$ is contained in $\mathbb{R} \times] -\infty, 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 \ge 0$ we can find $s_x > 0$ small enough such that $D(s_{in} - s_x) - \mu(s_x)x \ge 0$ for any $D \in ]0, D_{max}]$. Moreover, a trajectory of (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 \in ]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} \text{ and } 0 < x(t), \quad \forall t \in [t_0, T]. \qquad (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(\cdot)$, standard assumptions that guarantee the existence of optimal control can be evoked (for example [24, 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 \in ]0, s_{in}]$ and $x_0 > 0$, we will assume, unless otherwise stated, that optimal trajectories satisfy (3).

## 3.1 Pontryagin Maximum Principle

We begin our study of problem ($P_{bio}$) by establishing necessary conditions of optimality with the Pontryagin Maximum Principle (PMP) [24, Corollary 22.3]. For this we set the Hamiltonian $H: \mathbb{R}^2 \times \mathbb{R}^2 \times [0, D_{max}] \to \mathbb{R}$ as

$$H(x, s, p_x, p_s, D) := \mu(s)x + p_s(D(s_{in} - s) - \mu(s)x) + p_x(\mu(s)x - Dx). \qquad (4)$$

We consider an optimal control $D(\cdot)$ of ($P_{bio}$) and its associated states $x(\cdot)$ and $s(\cdot)$, solution of (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] \to \mathbb{R}$ satisfying, for almost every $t \in [t_0, T]$, the adjoint equations

$$\dot{p}_x = Dp_x - \mu(s)(1 + p_x - p_s), \quad \dot{p}_s = Dp_s - \mu'(s)x(1 + p_x - p_s), \qquad (5)$$
the transversality condition $p_x(T) = p_s(T) = 0$ and the maximum condition, for almost every $t \in [t_0, T]$,

$$H(x(t), s(t), p_x(t), p_s(t), D(t)) = \max_{d \in [0, D_{\text{max}}]} H(x(t), s(t), p_x(t), p_s(t), d). \quad (6)$$

In addition, since the Hamiltonian does not depend explicitly on time, it is constant, which means that for some $\mathbf{c} = \mathbf{c}(t_0, x_0, s_0) \in \mathbb{R}$ we have

$$H(x(t), s(t), p_x(t), p_s(t), D(t)) = \mathbf{c}, \quad \text{a.e. on } [t_0, T]. \quad (7)$$

This, along with the transversality condition, yields $\mathbf{c} = \mu(s(T))x(T) > 0$.

Let us call extremal trajectory and extremal control to any trajectory $(x(\cdot), s(\cdot), p_x(\cdot), p_s(\cdot))$ and control $D(\cdot)$ satisfying (1)-(5)-(6)-(7).

Since the Hamiltonian is affine in the control variable, an extremal control will depend on the sign of the commutation function

$$\phi(t) := \frac{\partial}{\partial D} H(x(t), s(t), p_x(t), p_s(t), D(t)) = p_s(t)(s_{\text{in}} - s(t)) - p_x(t)x(t).$$

We then have that $D(t) = 0$ if $\phi(t) < 0$ and $D(t) = D_{\text{max}}$ if $\phi(t) > 0$, while no information can be directly obtained from the PMP in the case $\phi(t) = 0$.

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

$$(s_{\text{in}} - s)p_s - xp_x = 0, \quad \mu'(s)(s_{\text{in}} - s)(1 + p_x - p_s) = \mu(s). \quad (8)$$

With this, we can get an equation that the state variables satisfy during a singular arc, that depends only on the constant value $\mathbf{c}$ of the Hamiltonian. Indeed, when the commutation function vanishes, we have

$$\mathbf{c} = \mu(s)x(1 + p_x - p_s)$$

and using (8) we get

$$\mathbf{c}\mu'(s)(s_{\text{in}} - s) = \mu(s)^2x. \quad (9)$$

We now define the following function, for $h > 0$ given

$$x_h(s) := h \frac{\mu'(s)(s_{\text{in}} - s)}{\mu(s)^2}, \quad 0 < s < s_{\text{in}}.$$

Then, from (9), we have that the extremal state trajectories during the singular arc remain in the graph of $s \mapsto x_c(s)$. This means that if we knew the value of $\mathbf{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_{\text{bio}})$.

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

$$x_h^M(s) := \frac{hK(s_{\text{in}} - s)}{\mu_{\text{max}}s^2}, \quad x_h^H(s) := \frac{hK(s_{\text{in}} - s)}{\mu_{\text{max}}s^2} - \frac{h(s_{\text{in}} - s)}{\mu_{\text{max}}K_i}.$$

**Remark 3.** Note that in general the curve $x_h(s) \to 0$ when $s \to s_{\text{in}}$, and $x_h(s) \to +\infty$ if $s \to 0$ provided that $\frac{\mu'(s)}{\mu(s)^2} \to +\infty$ as $s \to 0$; this is for instance the case of the Monod and Haldane growth rate functions.
## 3.2 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

$$\mu(s)x_h - D_h x_h = (D_h(s_{in} - s) - \mu(s)x_h)\partial_s x_h$$

and we then have the following expression for the control on the singular arc

$$D_h(s) = \frac{\mu(s)x_h(s)(1 + \partial_s x_h(s))}{x_h(s) + (s_{in} - s)\partial_s x_h(s)}, \quad (10)$$

with

$$\partial_s x_h(s) = h \left( \frac{d}{ds} \left( \frac{\mu'(s)}{\mu(s)^2} \right) (s_{in} - s) - \frac{\mu'(s)}{\mu(s)^2} \right), \quad \forall s \in ]0, s_{in}]. \quad (11)$$

Our Standing assumptions, in particular the fact that $\frac{d}{ds} \left( \frac{\mu'(s)}{\mu(s)^2} \right) \neq 0$ implies that $D_h(s)$ is well defined (as a real valued function) for any $s \in ]0, s_{in}]$.

The control $D_h$ is not necessarily an admissible control for the problem at hand. For some $s \in ]0, s_{in}]$ it could happen that $D_h(s) \notin [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 \quad \text{and} \quad D_h(s^{max}) = D_{max}. \quad (12)$$

**Remark 4.** It is straightforward to see that the control on the singular curve associated with a Monod growth rate function is given by

$$D_h^M(s) = \frac{hK(2s_{in} - s) - \mu_{max}s^3}{2s(K + s)(s_{in} - s)}, \quad \forall s \in ]0, s_{in}].$$

It follows that $D_h^M(s) \to +\infty$ if $s \to 0$. However, the behavior of $D_h^M(s)$ when $s \to s_{in}$ depends on the data of the problem. As a matter of fact

$$\lim_{s \to s_{in}} D_h^M(s) = \begin{cases} +\infty, & \text{if } hK > \mu_{max}s_{in}^2, \\ -\infty, & \text{if } hK < \mu_{max}s_{in}^2, \\ \frac{2\mu_{max}s_{in}}{K+s_{in}}, & \text{if } hK = \mu_{max}s_{in}^2. \end{cases}$$

This means that, depending on the data of the problem, singular optimal trajectories may not occur at all; for instance if $D_h^M(s) > D_{max}$ for every $s \in ]0, s_{in}]$. We plan to study this issue in more details and for general growth rate functions elsewhere.

For $(s, x) \notin I_s \times 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 := \left\{ \begin{array}{ll} (x_h(s), s) & : s \in I_s \\ (x(\tau, x_h(s^0), 0), s(\tau, s^0, 0)) & : \tau \le 0 \\ (x(\tau, x_h(s^{max}), D_{max}), s(\tau, s^{max}, D_{max})) & : \tau \le 0 \end{array} \right\} \quad (13)$$

where we denote $s(\tau, s^{max}, D_{max})$ the value at time $\tau$ of the solution with control $D_{max}$ starting at $s^{max}$ at time $\tau = 0$ and similarly for $x(\tau, x_h(s^{max}), D_{max})$, $x(\tau, x_h(s^0), 0)$ and $x(\tau, x_h(s^0), 0)$.

Note that with $D = 0$ we have $\dot{x} + \dot{s} = 0$ so that the corresponding trajectory $\{(x(\tau, x_h(s^0), 0), s(\tau, s^0, 0)) : \tau \le 0\}$ corresponds to the graph of the mapping $s \mapsto -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(\tau, x_h(s^0), 0), s(\tau, s^0, 0)) : \tau \le 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\}. \quad (14)$$
<figure><img src="image_3.png" /><figcaption>Figure 1: Example of trajectories obtained with feedback $\psi_h$ with $h = 2$ and Monod growth function ($\mu_{max} = 1.2$, $K = 7.1$), $D_{max} = 0.7$ and $s_{in} = 10$. The thicker lines correspond to the singular curve $G_h$.</figcaption></figure>

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_h^0$ (respectively $G_h^{max}$) 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 \in ]0, s_0]$. We then have the following expression for $G_h^0$

$$G_h^0 = \{(x, s) : \exists \tilde{s} \le s \text{ such that } (-\tilde{s} + x + s, \tilde{s}) \in G_h\}$$

and then $G_h^{max} := \mathbb{R}_+^2 \setminus (G_h^0 \cup G_h)$.

With this, we have the following family of feedback controls

$$\psi_h(x, s) = \begin{cases} 0, & \text{if } (x, s) \in G_h^0, \\ D_{\max}, & \text{if } (x, s) \in G_h^{\max}, \\ D_h(s), & \text{if } (x, s) \in G_h. \end{cases} \qquad (15)$$

An example of the trajectories obtained with this feedback is shown in Figure 1. Note that, because of the way the feedback control $\psi_h$ has been constructed, solution of related the closed-loop system do exist and are uniquely determined by the initial data.

## 3.3 Hamilton-Jacobi-Bellman Equation

In order to make the feedback law (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(\cdot)} \left\{ \int_{t_0}^{T} \mu(s(t))x(t)dt : \begin{array}{l} \dot{x} = \mu(s)x - Dx, \quad x(t_0) = x_0, \\ \dot{s} = D(s_{in} - s) - \mu(s)x, \quad s(t_0) = s_0 \end{array} \right\}$$
where the maximum is taken over all $D: [t_0, T] \to [0, D_{\max}]$ measurable. It is not difficult to see that, thanks to the continuity of the trajectories of the control system (1) with respect the initial data, the value function $(t_0, x_0, s_0) \mapsto 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, +\infty[ \times [0, s_{in}]$ with respect to the control system (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. [25]).

The HJB equation for the problem we are dealing with is

$$ \partial_t u + \sup_{D \in [0, D_{\max}]} H(x, s, \partial_x u, \partial_s u, D) = 0, \quad \text{in } ]-\infty, T[ \times \mathbb{R}^2, \quad (16) $$

where $H$ is the Hamiltonian given in (4). Existence and uniqueness of solutions in the viscosity sense for HJB equations is a well-known and studied fact, see for instance [26]. As a matter of fact, the value function $V$ is the unique viscosity solution to (16) that satisfies the terminal condition

$$ u(T, x, s) = 0, \quad \forall x, s \in \mathbb{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.

**Lemma 1** ([26, Theorem III.3.42]). Under the Standing Assumptions, a measurable functions $D: [t_0, T] \to [0, D_{\max}]$ maximizes $(P_{bio})$, the production of biogas problem, if and only if the maximum condition (6) holds and

$$ (c, p_x(t), p_s(t)) \in \partial^+ V(t, x(t), s(t)), \quad \text{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 \in [t_0, T]$ and

$$ \partial^+ u(z) := \left\{ q \in \mathbb{R}^n : \limsup_{y \to z} \frac{u(y) - u(z) - q^\top(y-z)}{|y-z|} \le 0 \right\} $$

stands for the viscosity superdifferential of a function $u: \mathbb{R}^n \to \mathbb{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)) = \nabla_{(x_0, s_0)} V(t_0, x_0, s_0) \quad \text{and} \quad c = -\partial_{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 $\psi_h$ given by (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), \psi_h(x_0, s_0)) $$

From the previous section, we can deduce that $\psi_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 $\psi_c$ is the optimal control for a particular choice of initial conditions<fnref n="1" />. If $\psi_c$ is actually an optimal control, we can write the value function as the cost of the control $\psi_c$, that is, $V(t_0, x_0, s_0) = J(t_0, x_0, s_0, \psi_c)$. With Lemma 1 we then get

$$ c = H(x_0, s_0, \nabla_{(x_0, s_0)} J(t_0, x_0, s_0, \psi_c), \psi_c(x_0, s_0)) $$

<footnote n="1">¹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.</footnote>
Thus, using the HJB equation (16), we also have

$$ \mathbf{c} = -\partial_{t_0} J(t_0, x_0, s_0, \psi_c) \qquad (17) $$

In other words, c is a fixed point of the mapping $h \mapsto -\partial_{t_0}J(t_0, x_0, s_0, \psi_h)$. Hence, if we are able to compute or approximate a fixed point of the mapping $h \mapsto -\partial_{t_0}J(t_0, x_0, s_0, \psi_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 (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 $\psi_c$ is an approximated solution to the HJB equation.

# 4 An algorithm for Maximizing the Production of Biogas

We present now a way to compute the extremal feedback control $\psi_c$ by solving the fixed point equation (17) numerically, in order to get the value of the Hamiltonian c for any initial condition $(x_0, s_0) \in ]0, +\infty[ \times ]0, s_{in}[$.

## 4.1 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(\mathbf{h}) = \mathbf{h}$ for some given mapping $F: \mathbb{R}^d \to \mathbb{R}^d$, then starting from an initial guess $\mathbf{h}_0$, the numerical scheme consists in computing $\mathbf{h}_{n+1} = F(\mathbf{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 $\varepsilon$, specifically when $||\mathbf{h}_{n+1} - \mathbf{h}_n|| < \varepsilon$.

In our case, the function for which we need to compute a fixed point, $-\partial_{t_0}J(t_0, x_0, s_0, \psi_h)$, is composed of a partial derivative and therefore to estimate it numerically with a finite difference approximation, we must work with a discrete range of initial times. For $N \in \mathbb{N}$, we denote $\{t_0^k\}_{k=1,...,N}$ a set of initial times with constant step $\Delta t_0 = t_0^{k+1} - t_0^k$. We will therefore compute a vector of fixed points $\tilde{\mathbf{c}} = (\mathbf{c}^k) \in \mathbb{R}^N$, where each $\mathbf{c}^k$ will correspond to the value of the Hamiltonian for the initial time $t_0^k$.

To obtain $\mathbf{h}_{n+1}$ from the previous iterate $\mathbf{h}_n = (h_n^k)$ we start by computing for each $t_0^k$ the trajectories with the control $\psi_{h_n^k}$ and the associated cost $J(t_0^k, x_0, s_0, \psi_{h_n^k})$ with standard numerical integration tools. For this the singular curve $G_{h_n^k}$ must be first determined by solving equations (12) to establish the admissible range $[s^{max}, s^0]$ and then integrating backwards to obtain $(x(\cdot, x_h(s^{max}), D_{max}), s(\cdot, s^{max}, D_{max}))$.

We can then approximate the partial derivative of the cost as

$$ \partial_{t_0} J(t_0^k, x_0, s_0, \psi_{h_n^k}) \approx \frac{J(t_0^{k+1}, x_0, s_0, \psi_{h_n^{k+1}}) - J(t_0^k, x_0, s_0, \psi_{h_n^k})}{\Delta t_0} $$

and if we set $t_0^N = T$, we can use that $J(T, x_0, s_0, \psi_h) = 0$ to start the computations of these partial derivatives, running through the range of initial times backwards beginning with $t_0^{N-1}$ and ending with $t_0^1$.

In summary, for a fixed initial condition $(x_0, s_0) \in ]0, +\infty[ \times ]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 $\partial_{t_0}J(t_0, x_0, s_0, \psi_h)$.
**Algorithm 1:**

Input: $N$, MaxIterations, $\varepsilon$, $\mathbf{h}_0$

for $n < \text{MaxIterations}$ do

for $k = N - 1, \dots, 1$ do

Solve $D_{h_n^k}(s^0) = 0$ and $D_{h_{\max}^k}(s^{\max}) = D_{\max}$

Compute singular curve $G_{h_n^k}$

Compute $x(\cdot, t_0^k, x_0, s_0, \psi_{h_n^k})$, $s(\cdot, t_0^k, x_0, s_0, \psi_{h_n^k})$ and $J(t_0^k, x_0, s_0, \psi_{h_n^k})$

$h_{n+1}^k \leftarrow -(J(t_0^{k+1}, x_0, s_0, \psi_{h_n^{k+1}}) - J(t_0^k, x_0, s_0, \psi_{h_n^k}))/\Delta t_0$

if $||h_{n+1} - h_n|| < \varepsilon$ then

return $h_{n+1}$

Indeed, since we need $J(t_0^{k+1}, x_0, s_0, \psi_{h_n^{k+1}})$ to compute $h_{n+1}^k$, any errors made to get $h_n^{k+1}$ will propagate to $h_{n+1}^k$ and all following values $h_{n+1}^j$ for $j < k$. Another consequence of this inter-dependance is that to have $h_n^k$ converge, $h_n^{k+1}$ must have already converged to a fixed point.

With these considerations in mind, it might seem unnecessary to compute the whole vector $\mathbf{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^{N-1}$ until convergence, which is possible because we only need $J(T, x_0, s_0, \psi_h)$. Then using the obtained fixed point to get $J(t_0^{N-1}, x_0, s_0, \psi_{c^{K-1}})$ we could move on to computing $\mathbf{c}^{K-2}$. Repeating this process, we can thus find all the fixed points $\mathbf{c}^k$ until reaching the desired initial time. However, due to the accumulation of errors, to get the convergence of $h_n^j$ with a certain tolerance it is necessary to get the convergence of $h_n^k$ for $j < k$ with a smaller tolerance. It is then complicated in practice to determine an efficient stoping condition that guarantees the convergence of the last fixed point to the desired tolerance.

On the other hand, by computing the whole vector $\mathbf{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_{n+1}^k - h_n^k| < \varepsilon$. 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^{N-1}$, $h_n^{N-2}$, ...), 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_n^1$, $h_n^2$, ...).

## 4.2 Initial Guess

Concerning the initial guess $\mathbf{h}_0$, recall that with the optimal control we have $\mathbf{c} = \mu(s(T))x(T)$. In most cases, we can broadly approximate this by taking $h_0^k = \mu(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

$$ \mathbf{c} = \mu(s(T, t_0, x_0, s_0, \psi_c))x(T, t_0, x_0, s_0, \psi_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 \leftarrow \mu(s_0)x_0$

for $n < \text{MaxIterations}$ do

Solve $D_{h_n}(s^0) = 0$ and $D_{h_{\max}}(s^{\max}) = D_{\max}$

Compute singular curve $G_{h_n}$

Compute $x(\cdot, t_0, x_0, s_0, \psi_{h_n})$, $s(\cdot, t_0, x_0, s_0, \psi_{h_n})$

$h_{n+1} \leftarrow \mu(s(T, t_0, x_0, s_0, \psi_{h_n}))x(T, t_0, x_0, s_0, \psi_{h_n})$

if $|h_{n+1} - h_n| < \varepsilon$ then

return $h_{n+1}$

<figure><img src="image_4.png" /><figcaption>Figure 2: Optimal trajectories in state space, for the Monod growth function ($\mu_{\max} = 1.2$, $K = 7.1$) with $t_0 = 0$, $T = 2.5$, $D_{\max} = 0.7$ and $s_{\text{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)\}$.</figcaption></figure>

5 Numerical Simulations

In this section, we illustrate the fixed point algorithm with the growth functions of Monod and Haldane (2) with parameter values from [27].

In Figures 2 and 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 \in [0, s_{\text{in}}]} \mu(s)(s_{\text{in}} - s).$$

This is expected since $c = \mu(s(T))x(T)$ and using the expression for the singular curve (9) evaluated at final time we get $\mu'(s(T))(s_{\text{in}} - s(T)) = \mu(s(T))$ and we recognize this as a necessary condition for maximizing $s \mapsto \mu(s)(s_{\text{in}} - s)$.

Next, Figures 4, 5 and 6 each show optimal trajectories in state space for various 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 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 bang-bang-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 \mapsto c(t_0)$ appears to be continuous it
<figure><img src="image_5.png" /><figcaption>Figure 3: Optimal trajectories in state space, for the Haldane growth function ($\mu_{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)}.</figcaption></figure>

<table><caption>Table 1: Performance comparison with Bocop</caption><thead><tr><th>(x<sub>0</sub>, s<sub>0</sub>)</th><th>Biogas (ψ<sub>c</sub>)</th><th>Biogas (Bocop)</th><th>Relative Difference</th></tr></thead><tbody><tr><td>(3, 2)</td><td>3.2232</td><td>3.2235</td><td>9 · 10<sup>-5</sup></td></tr><tr><td>(3, 6)</td><td>5.3285</td><td>5.3290</td><td>9 · 10<sup>-5</sup></td></tr><tr><td>(1, 4.5)</td><td>1.8725</td><td>1.8729</td><td>2 · 10<sup>-4</sup></td></tr><tr><td>(7, 8)</td><td>4.7904</td><td>4.7933</td><td>6 · 10<sup>-4</sup></td></tr></tbody></table>

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 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 \mapsto 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 7 with graphs of error as function of iterations: $n \mapsto |h_{n+1}^k - h_n^k|$. 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 $\partial_{t_0}J(t_0, x_0, s_0, \psi_h)$ since this behaviour is not seen for Algorithm 2 which does not need the computation of $\partial_{t_0}J(t_0, x_0, s_0, \psi_h)$. Notice however, that Algorithm 2 requires more iterations but that the convergence accelerates at the end and reaches machine precision.

## 5.1 Comparison with Bocop

Finally, we compare our feedback to the control obtained with the open source toolbox for optimal control Bocop [16, 28]. This package implements a direct method that approximates the optimal control problem by a finite dimensional optimization problem using a time discretization.

Table 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 4 and 5.

Next, in Table 2 we show some computational times associated with our feedback (for various error tolerances $\epsilon$) 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
<table><caption>Table 2: Computation time (in seconds) comparison with Bocop</caption><thead><tr><th>(x<sub>0</sub>, s<sub>0</sub>)</th><th>CPU time (ψ<sub>c</sub>) ε = 10<sup>-4</sup></th><th>CPU time (ψ<sub>c</sub>) ε = 10<sup>-6</sup></th><th>CPU time (Bocop)</th></tr></thead><tbody><tr><td>(3, 2)</td><td>2.32</td><td>3.60</td><td>2.06</td></tr><tr><td>(3, 6)</td><td>1.37</td><td>1.74</td><td>1.74</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 4: Initial condition $(x_0, s_0) = (7, 8)$, on the left, optimal trajectories in state space and on the right, $t_0 \mapsto c(t_0)$ value of the Hamiltonian as a function of initial time. Monod growth function ($\mu_{max} = 1.2, K = 7.1$) with $T = 1, D_{max} = 0.7$ and $s_{in} = 10$.</figcaption></figure>

the control for a single initial time, where as our algorithm for a range of initial times.

# 6 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 perturbations 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 (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.

# Acknowledgements

A. Haddon was supported by a doctoral fellowship CONICYT - PFCHA / Doctorado Nacional / 2017-21170249. The first author was also supported by FONDECYT grant 1160567, and by Basal Program CMM-AFB 170001 from CONICYT- Chile. C. Hermosilla was supported by CONICYT- Chile through FONDECYT grant number 3170485.
<figure><img src="image_7.png" /><figcaption>Figure 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 \mapsto c(t_0)$ value of the Hamiltonian as a function of initial time. Monod growth function ($\mu_{\max} = 1.2, K = 7.1$) with $T = 2.5, D_{\max} = 0.7$ and $s_{\text{in}} = 10$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Initial condition $(x_0, s_0) = (200, 70)$, on the left, optimal trajectories in state space and on the right $t_0 \mapsto c(t_0)$ value of the Hamiltonian as a function of initial time. Haldane growth function ($\mu_{\max} = 0.74, K = 9.28, K_i = 256$) with $T = 2, D_{\max} = 1$ and $s_{\text{in}} = 100$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Error as a function of iterations: $n \mapsto |h_{n+1}^k - h_n^k|$ of Algorithm 1 on the left and of Algorithm 2 on the right. Monod growth function ($\mu_{\max} = 1.2, K = 7.1$) with $T = 1, D_{\max} = 0.7$ and $s_{\text{in}} = 10$. On the left initial condition ($x_0, s_0$) = (2, 2) and on the left ($x_0, s_0$) = (7, 8)</figcaption></figure>

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_{\text{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) \in \mathbb{R}: x + s = s_{\text{in}}\}$ is invariant for the dynamics (1) and thus, for initial conditions in $I$, the dynamics reduce to $\dot{s} = (D - \mu(s))(s_{\text{in}} - s)$. This special case was solved in [13], with the following assumptions

(H1) The function $s \mapsto \mu(s)(s_{\text{in}} - s)$ has a unique maximizer $s^*$ on $[0, s_{\text{in}}]$.

(H2) The upper bound on the controls is such that $D_{\max} > \mu(s^*)$.

The optimal control is then $D^*(s) = 0$ if $s > s^*$, $D^*(s) = D_{\max}$ if $s < s^*$ and $D^*(s) = \mu(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 (15) where the singular arc is reduced to $s = s^*$. In other words, it corresponds to the control $\psi_{h^*}$ where $h^*$ satisfies equation (9) for the singular arc with $s = s^*$, which in this case is $h^*\mu'(s^*) = \mu(s^*)^2$. Next, since $s^*$ is a maximizer we have $\mu'(s^*)(s_{\text{in}} - s^*) - \mu(s^*) = 0$ and therefore $h^* = \mu(s^*)(s_{\text{in}} - s^*)$.

To prove the optimality of $\psi_{h^*}$, we must now show that $h^*$ is a fixed point of the mapping $h \mapsto -\partial_{t_0}J(t_0, x_0, s_0, \psi_h)$. For this we first study the trajectories obtained with the feedback control $\psi_{h^*}$. We denote in the remainder of the section the right-hand side of the differential equation for $s(\cdot)$ with control $\psi_{h^*}$ as $f(s) := (\psi_{h^*}(s) - \mu(s))(s_{\text{in}} - s)$.

Notice that for $s > s^*$ we have $f(s) = -\mu(s)(s_{\text{in}} - s) < 0$ and for $s < s^*$ we have $f(s) = (D_{\max} - \mu(s))(s_{\text{in}} - s) > 0$ from assumption (H2). Thus, $s^*$ is reachable from any initial condition in $I$ with control $\psi_{h^*}$. We define the time $t^*$ when $s^*$ is reached, from a given initial condition $s_0 \in [0, s_{\text{in}}]$ and initial time $t_0$ with control $\psi_{h^*}$, that is, $t^* := \inf\{t \ge t_0 : s(t, t_0, s_0, \psi_{h^*}) = s^*\}$. Finally, note that with control $D = \mu(s^*)$ the point $s = s^*$ becomes a steady state. Therefore the trajectories with control $\psi_{h^*}$ are

$$s(t) = \begin{cases} s(t, t_0, s_0, \psi_{h^*}), & \text{for } t_0 \le t \le \min(t^*, T), \\ s^*, & \text{for } \min(t^*, T) \le t \le T. \end{cases}$$

We can now compute $\partial_{t_0}J(t_0, x_0, s_0, \psi_{h^*})$ and for this we need the following.
**Lemma 2.** For any initial condition $(x_0, s_0) \in I$, for the trajectories with control $\psi_{h^*}$ we have $\partial_{t_0}s(t) = -f(s(t))$ at time $t \in [t_0, t^*]$.

*Proof.* We can write the differential equation satisfied by $s(\cdot)$ as $s(t) = s_0 + \int_{t_0}^{t} f(s(\tau)) d\tau$ and differentiating we get $\partial_{t_0}s(t) = -f(s_0) + \int_{t_0}^{t} f'(s(\tau))\partial_{t_0}s(\tau) d\tau$. This is a linear differential equation and the solution is $\partial_{t_0}s(t) = -f(s_0) \exp\left(\int_{t_0}^{t} f'(s(\tau)) d\tau\right)$. Now, as $f(s(t))$ does not change sign for $t \in [t_0, t^*)$ and since $f(s(t))$ is the derivative of $s(t)$ we have $\int_{t_0}^{t} f'(s(\tau)) d\tau = \int_{s_0}^{s(t)} \frac{f'(s)}{f(s)} ds = \ln\left(\frac{f(s(t))}{f(s_0)}\right)$. $\square$

We are now in a position to prove the optimality of the feedback control proposed earlier.

**Proposition 1.** For any initial condition $(x_0, s_0) \in I$ and for any initial time $t_0$ such that $s^*$ is reachable, that is when $t^* \le T$, we have $\partial_{t_0}J(t_0, x_0, s_0, \psi_{h^*}) = -\mu(s^*)(s_{in} - s^*)$, so that $\psi_{h^*}$ is the optimal control.

*Proof.* We start by writing the cost as

$$J(t_0, x_0, s_0, \psi_{h^*}) = \int_{t_0}^{t^*} \mu(s(t))(s_{in} - s(t)) dt + (T - t^*)\mu(s^*)(s_{in} - s^*)$$

differentiating with respect to $t_0$ we get

$$\partial_{t_0}J(t_0, x_0, s_0, \psi_{h^*}) = -\mu(s_0)(s_{in} - s_0) + \int_{t_0}^{t^*} \partial_s(\mu(s(t))(s_{in} - s(t)))\partial_{t_0}s(t) dt.$$

Note that the terms with $\partial_{t_0}t^*$ cancel out because $s(t^*) = s^*$. Now, using Lemma 2 we get

$$\begin{align*}
\partial_{t_0} J(t_0, x_0, s_0, \psi_{h^*}) &= -\mu(s_0)(s_{in} - s_0) - \int_{t_0}^{t^*} \partial_s(\mu(s(t))(s_{in} - s(t)))\dot{s}(t) dt \\
&= -\mu(s_0)(s_{in} - s_0) - \int_{t_0}^{t^*} \frac{d}{dt}(\mu(s(t))(s_{in} - s(t))) dt \\
&= -\mu(s^*)(s_{in} - s^*).
\end{align*}$$

We conclude by recalling that $h^* = \mu(s^*)(s_{in} - s^*)$ and therefore $h^*$ is a fixed point of $h \mapsto -\partial_{t_0}J(t_0, x_0, s_0, \psi_h)$.

$\square$

## References

[1] Rehl, T., Müller, J.: Co2 abatement costs of greenhouse gas (ghg) mitigation by different biogas conversion pathways. Journal of environmental management 114, 13–25 (2013)

[2] Beddoes, J.C., Bracmort, K.S., Burns, R.T., Lazarus, W.F.: An analysis of energy production costs from anaerobic digestion systems on us livestock production facilities. USDA NRCS Technical Note (1) (2007)

[3] Guwy, A., Hawkes, F., Wilcox, S., Hawkes, D.: Neural network and on-off control of bi-carbonate alkalinity in a fluidised-bed anaerobic digester. Water Research 31(8), 2019–2025 (1997)
[4] Nguyen, D., Gadhamshetty, V., Nitayavardhana, S., Khanal, S.K.: Automatic process control in anaerobic digestion technology: A critical review. *Bioresource technology* **193**, 513–522 (2015)

[5] García-Diéguez, C., Molina, F., Roca, E.: Multi-objective cascade controller for an anaerobic digester. *Process Biochemistry* **46**(4), 900–909 (2011)

[6] Rodríguez, J., Ruiz, G., Molina, F., Roca, E., Lema, J.: A hydrogen-based variable-gain controller for anaerobic digestion processes. *Water Science and Technology* **54**(2), 57–62 (2006)

[7] Djatkov, D., Effenberger, M., Martinov, M.: Method for assessing and improving the efficiency of agricultural biogas plants based on fuzzy logic and expert systems. *Applied energy* **134**, 163–175 (2014)

[8] Dimitrova, N., Krastanov, M.: Nonlinear adaptive stabilizing control of an anaerobic digestion model with unknown kinetics. *International Journal of Robust and Nonlinear Control* **22**(15), 1743–1752 (2012)

[9] Sbarciog, M., Loccfuier, M., Wouwer, A.V.: An optimizing start-up strategy for a bio-methanator. *Bioprocess and biosystems engineering* **35**(4), 565–578 (2012)

[10] Bayen, T., Cots, O., Gajardo, P.: Analysis of an optimal control problem related to the anaerobic digestion process. *Journal of Optimization Theory and Applications* (2018)

[11] Bernard, O., Chachuat, B., Hélias, A., Rodriguez, J.: Can we assess the model complexity for a bioprocess: theory and example of the anaerobic digestion process. *Water science and technology* **53**(1), 85–92 (2006)

[12] Stamatelatou, K., Lyberatos, G., Tsiligiannis, C., Pavlou, S., Pullammanappallil, P., Svoronos, S.: Optimal and suboptimal control of anaerobic digesters. *Environmental Modeling & Assessment* **2**(4), 355–363 (1997)

[13] Ghouali, A., Sari, T., Harmand, J.: Maximizing biogas production from the anaerobic digestion. *Journal of Process Control* **36**, 79–88 (2015)

[14] Haddon, A., Harmand, J., Ramírez, H., Rapaport, A.: Guaranteed value strategy for the optimal control of biogas production in continuous bio-reactors. *IFAC-PapersOnLine* **50**(1), 8728–8733 (2017)

[15] Haddon, A., Ramírez, H., Rapaport, A.: First results of optimal control of average biogas production for the chemostat over an infinite horizon. *IFAC-PapersOnLine* **51**(2), 725–729 (2018)

[16] Team Commands, Inria Saclay: Bocop: an open source toolbox for optimal control (2017). URL http://www.bocop.org

[17] Gerdts, M.: Optimal Control and Parameter Identification with Differential-Algebraic Equations of index 1: User's Guide. Tech. rep., Institut für Mathematik und Rechneranwendung, Universität der Bundeswehr München (2011)

[18] Bokanowski, O., Desilles, A., Zidani, H.: ROC-HJ-Solver. a C++ Library for Solving HJ equations (2013). URL http://uma.ensta-paristech.fr/soft/ROC-HJ/

[19] Yane (2018). URL http://www.nonlinearmpc.com

[20] Lobry, C., Rapaport, A., Sari, T., et al.: The Chemostat: Mathematical Theory of Microorganism Cultures. John Wiley & Sons (2017)
[21] Bastin, G., Dochain, D.: On-line estimation and adaptive control of bioreactors. Elsevier (1991)

[22] Hermosilla, C.: Stratified discontinuous differential equations and sufficient conditions for robustness. Discrete and Continuous Dynamical Systems-Series A **35**(9), 23 (2015)

[23] Clarke, F.H., Ledyaev, Y.S., Stern, R.J., Wolenski, P.R.: Nonsmooth analysis and control theory, vol. 178. Springer Science & Business Media (2008)

[24] Clarke, F.: Functional analysis, calculus of variations and optimal control, vol. 264. Springer Science & Business Media (2013)

[25] Hermosilla, C., Zidani, H.: Infinite horizon problems on stratifiable state-constraints sets. Journal of Differential Equations **258**(4), 1430–1460 (2015)

[26] Bardi, M., Capuzzo-Dolcetta, I.: Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations. Springer Science & Business Media (2008)

[27] Bernard, O., Hadj-Sadok, Z., Dochain, D., Genovesi, A., Steyer, J.P.: Dynamical model development and parameter identification for an anaerobic wastewater treatment process. Biotechnology and bioengineering **75**(4), 424–438 (2001)

[28] Bonnans Frederic, J., Giorgi, D., Grelard, V., Heymann, B., Maindrault, S., Martinon, P., Tissot, O., Liu, J.: Bocop – A collection of examples. Tech. rep., INRIA (2017). URL http://www.bocop.org