# Gradient-based optimization of a rotating algal biofilm process

Pierre-Olivier Lamare, Nina Aguillon, Jacques Sainte-Marie, Jérôme Grenier, Hubert Bonnefond, Olivier Bernard

► To cite this version:

Pierre-Olivier Lamare, Nina Aguillon, Jacques Sainte-Marie, Jérôme Grenier, Hubert Bonnefond, et al. "Gradient-based optimization of a rotating algal biofilm process." [Research Report] RR-9250, Inria - Sophia antipolis. 2019. hal-01990002

HAL Id: hal-01990002

https://inria.hal.science/hal-01990002

Submitted on 23 Jan 2019

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.
informatics mathematics
Inria

Gradient-based
optimization of a rotating
algal biofilm process

Pierre-Olivier Lamare, Nina Aguillon, Jacques Sainte-Marie, Jérôme Grenier, Hubert Bonnefond, Olivier Bernard

RESEARCH
REPORT

Nº 9250

January 2019

Project-Teams ANGE and
BIOCORE
<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# Gradient-based optimization of a rotating algal biofilm process

Pierre-Olivier Lamare, Nina Aguillon, Jacques Sainte-Marie, Jérôme Grenier, Hubert Bonnefond, Olivier Bernard

Project-Teams ANGE and BIOCORE

Research Report n° 9250 — January 2019 — 22 pages

**Abstract:** Microalgae are microorganisms which have been only recently used for biotechnological applications, especially in the perspective of biofuel production. Here we focus on the shape optimization and optimal control of an innovative process where the microalgae are fixed on a support. They are thus successively exposed to light and dark conditions. The resulting growth can be represented by a dynamical system describing the denaturation of key proteins due to an excess of light. A Partial Differential Equations (PDE) model of the Rotating Algal Biofilm (RAB) is then proposed, representing local microalgal growth submitted to the time varying light. An adjoint-based gradient method is proposed to identify the optimal (constant) process folding and the (time varying) velocity of the biofilm. When applied to a realistic case, the optimization points out a particular configuration which significantly increases the productivity compared to a base case where the biofilm is fixed.

**Key-words:** Adjoint-based optimization; PDE Control; microalgae; biofilm; biofuel; Rotating Algal Biofilm;

RESEARCH CENTRE
SOPHIA ANTIPOLIS – MÉDITERRANÉE

2004 route des Lucioles - BP 93
06902 Sophia Antipolis Cedex
Optimisation par méthode de gradient d'un procédé de
biofilm rotatif d'algues

**Résumé** : Les micro-algues sont des micro-organismes qui ont été que très récemment utilisés pour des applications bio-technologiques et plus précisément pour la production de bio-carburant. Dans ce rapport nous nous intéressons à l'optimisation de forme et au contrôle optimal d'un procédé innovant où les micro-algues sont fixées sur un support. Elles sont successivement exposées à la lumière et à l'obscurité. Le taux de croissance qui en résulte peut être représenté par un système dynamique décrivant la dénaturation de protéines clés dûe à un excès de lumière. Un modèle d'Équations aux Dérivées Partielles (EDP) pour le Biofilm Rotatif d'Algues (BRA) est proposé. Il représente la croissance locale des micro-algues soumis à une lumière variante en temps. Une méthode de gradient basée sur le calcul de l'adjoint du modèle est proposée pour identifier le repliement (constant) optimal du procédé et la vitesse (variante en temps) du bio-film. Une fois cette méthode utilisée dans un cas réaliste, l'optimisation donne une configuration qui améliore significativement la productivité comparativement au cas où le biofilm est fixe.

Mots-clés : Optimisation par méthode de gradient; contrôle d'EDP; micro-algues; biofilm; bio-carburant; Biofilm Rotatif d'Algues;
# 1 Introduction

Algae are microorganisms which were so far rarely used for biotechnological applications. In the last decade, their potential for production of bioproducts has been highlighted for addressing different markets such as pharmaceutics, food, feed, or even, at longer horizon, green chemistry and biofuels [24]. These organisms are generally cultivated under a planktonic form in open or closed photobioreactors. But the cell concentration remains low since for higher biomass light cannot penetrate through the very turbid medium [21]. Cell biomass thus only represent a tiny fraction of the liquid, generally around 0.1%. A high energy demand is then required by the harvesting and dewatering phases, and deeply limit the economic and environmental sustainability of these classical algal culturing process [11, 15]. In the recent years, surface attached microalgae biofilms have emerged [5]. They consist in growing microalgal cells fixed on a support and forming a biofilm. A biofilm is an assemblage of microbial cells that are irreversibly associated with a surface and enclosed in a matrix of extracellular polymeric substances. Biofilms are ubiquitous in nature, covering all kinds of surfaces, such as rocks, plants and sediments in seawater and freshwater environments. Algae in a biofilm is much more concentrated and straightforwardly harvested by a simple scrapping mechanism. Recently, innovative processes have been proposed where the algal biofilm is fixed on a rotating support so that cells are successively submitted to illuminated and dark periods (see Figure 1). The main advantage of this so called Rotating Algal Biofilm (RAB) is that it can better manage the way light is supplied to the microalgae and the subsequent dynamics of the photosystems harvesting light can be better exploited. When exposed during a too long period to high light intensity, some key proteins in the photosystems are denatured and must be resynthesized, which penalize the photosynthetic activity [16].

<figure><img src="image_3.png" /><figcaption>Figure 1: Schematic representation of the RAB. In this case the belt folding ratio is 1/3.</figcaption></figure>

The purpose of this work is to propose a first model of this new process, accounting for the dynamical character of the photosynthetic apparatus when submitted to varying light. The microalgae growth rate and the related productivity are given by the Han model representing the dynamics of the photosystems (Section 2). Roughly speaking productivity accounts for how much the biomass is growing. Then, a Partial Differential Equations (PDE) model based on the former derivations is designed for the RAB (Section 3). This nonlinear infinite dimensional model is further used to derive a gradient based strategy to estimate, depending on the Photon Flux Density (PFD) the optimal *belt folding ratio* (ratio of belt which is enlightened) and velocity (Section 4) with respect to a productivity criterion. This latter method is powerful for large system with a large number of control variables and sufficiently general to be applied to a wide variety of control optimization problem with PDEs constraints.

Originally, the adjoint method has been developed within the optimal control community with the seminal work of Lions [14]. Therefore, various control problems are tackled with this method
as ramp metering [22], air-traffic control [1], or inverse pressure design and drag minimization problems [17]. Moreover, this method has been widely spread through the data assimilation community (especially for state and parameter identification in weather forecast). One of the first work introducing the adjoint method for meteorological problem is [12]. Very recently, the method has been applied to the state and parameter estimation for traffic flow and overland flow system in [18]. Finally, there exist alternative approaches for gradient-based optimization of biological systems such as [19].

Here we use the classical methodology of gradient-based optimization but in the non-trivial framework of a very complex, nonlinear, realistic, and innovative model. The closest works in terms of the control problem are [1, 25]. In [25] the velocities vary in space and time, so that a relaxation technique can be used to transform the problem into a convex problem in the control variable. However, in the RAB the velocity cannot be space dependent due to the rigidity of the conveyor belt. Finally, the main difference of our work compared to [1], beyond the application, is the treatment of the control constraints. In [1], the cost function is augmented to penalize the violation of the constraints. Here, we favoured a projection step during the descent algorithm to satisfy the constraints.

## 2 Photoinhibition Dynamics

Photosystems are the functional units associated with the light part of photosynthesis. They are made of protein complexes and absorb light to generate a flux of energy and electrons. There are two types of photosystems in series PSI and PSII (see [6, 26]). PSII is generally assumed to play the main role in the photosynthetic dynamics. In the Han model, PSII is described by three different states: open and ready to harvest a photon (A), closed while processing the absorbed photon energy (B), or inhibited if two many photons have been absorbed simultaneously (C). Their respective evolutions satisfy the following dynamical system

$$ \frac{dA(t)}{dt} = -\sigma I(t)A(t) + \frac{B(t)}{\tau} \qquad (1) $$

$$ \frac{dB(t)}{dt} = \sigma I(t)A(t) - \frac{B(t)}{\tau} + k_r C(t) - k_d \sigma I(t)B(t) \qquad (2) $$

$$ \frac{dC(t)}{dt} = -k_r C(t) + k_d \sigma I(t) B(t), \qquad (3) $$

where A, B, and C are the relative frequencies of the three possible states and I is a continuous time-varying light signal. In the sequel we drop the t argument to alleviate the equations. Besides, the sum of A, B, and C is equal to 1 allowing to rewrite system (1)–(3) as

$$ \frac{dA}{dt} = - \left( \sigma I(t) + \frac{1}{\tau} \right) A + \frac{(1-C)}{\tau} \qquad (4) $$

$$ \frac{dC}{dt} = k_d \left( - \left( \frac{k_r}{k_d} + \sigma I(t) \right) C + \sigma I(t) (1-A) \right). \qquad (5) $$

We complete system (4) and (5) with initial conditions

$$ (A(0, \cdot), C(0, \cdot)) = (A^0, C^0) \in E, \qquad (6) $$

where $E = \{(A, C) \in \mathbb{R}_+^2 | A + C \in [0, 1]\}$. Photosynthetic production is described by the transition between open state and closed state. Excitation is assumed to occur at a rate $\sigma I$, with $\sigma$ [m$^2\mu E^{-1}$] being the effective cross section of the photosynthetic unit (PSU), whereas deexcitation
is assumed to occur at a rate of $\frac{1}{\tau}$, with $\tau$ [s] the turnover time of the electron transport chain. Photoinhibition occurring at high light irradiance corresponds to the transition from closed state to inhibited state. This process is assumed to occur at a rate of $k_d \sigma I$, with $k_d$ [–] a damage constant. The reverse transition from inhibited state to closed state accounts for the repair of damaged photosynthetic units by enzymatic processes in the cell, a mechanism that is assumed to occur at a constant rate $k_r$ [s$^{-1}$]. The photosynthesis rate is thus proportional to $\sigma I A$.

The dynamics of resting state A is fast compared to the dynamics of the photoinhibition C (due to the presence of the small multiplicative parameter $k_d$ in (5), see Table 1 below for an example of a possible value for this parameter). Since we focus on light variations at time scales larger than second, we can do a slow-fast approximation, and assume that A rapidly reaches its steady state [7] and focus on the dynamics of the photoinhibition rate C. To lead this analysis we use the singular perturbation theory (see e.g. [10]). We define

$$K_1 = \varepsilon^{-1}k_r \qquad (7)$$

$$K_2 = \varepsilon^{-1}k_d, \qquad (8)$$

where $\varepsilon$ is a small parameter of magnitude order of $k_r$ and $k_d$. Using parameters (7) and (8) equation (5) reads

$$\frac{dC}{dt} = \varepsilon K_2 \left[ - \left( \frac{K_1}{K_2} + \sigma I(t) \right) + \sigma I(t) (1 - A) \right]. \qquad (9)$$

Then, with the change of time

$$\tilde{t} = \varepsilon t, \qquad (10)$$

and by defining $\tilde{A}(\tilde{t}) = A(t)$, $\tilde{C}(\tilde{t}) = C(t)$, and $\tilde{I}(\tilde{t}) = I(t)$ we get the following system

$$\varepsilon \frac{d\tilde{A}}{dt} = -\left(\sigma\tilde{I} + \frac{1}{\tau}\right)\tilde{A} + \frac{1-\tilde{C}}{\tau} \qquad (11)$$

$$\frac{d\tilde{C}}{dt} = K_2 \left[ - \left( \frac{K_1}{K_2} + \sigma \tilde{I} \right) \tilde{C} + \sigma \tilde{I} (1 - \tilde{A}) \right], \qquad (12)$$

with the initial conditions

$$\tilde{A}(0) = A(0) = A^0 \qquad (13)$$

$$\tilde{C}(0) = C(0) = C^0. \qquad (14)$$

Setting $\varepsilon = 0$ in (11) we get

$$\tilde{A} = \frac{1 - \tilde{C}}{\tau \sigma \tilde{I} + 1}. \qquad (15)$$

Then, we substitute $\tilde{A}$ in (12) by its expression given in (15) we get the reduced dynamics for $\tilde{C}$

$$\frac{d\tilde{C}}{dt} = K_2 \left[ - \left( \frac{K_1}{K_2} + \frac{\tau (\sigma \tilde{I})^2}{\tau \sigma \tilde{I} + 1} \right) \tilde{C} + \frac{\tau (\sigma \tilde{I})^2}{\tau \sigma \tilde{I} + 1} \right]. \qquad (16)$$

Let us define the boundary layer correction as $\tilde{A}(t) = \tilde{A}(\tilde{t}) - \frac{1-\tilde{C}(\tilde{t})}{\tau\sigma\tilde{I}(t)+1}$. We have

$$\frac{d\tilde{A}(t)}{dt} = \frac{d\tilde{A}(\tilde{t})}{dt} - \frac{d}{dt} \left( \frac{1 - \tilde{C}(\tilde{t})}{\tau \sigma \tilde{I}(\tilde{t}) + 1} \right).$$
Therefore, we have

$$ \frac{d\hat{A}(t)}{dt} = -\left(\sigma\tilde{I}(\tilde{t}) + \frac{1}{\tau}\right)\tilde{A}(\tilde{t}) + \frac{1-\tilde{C}(\tilde{t})}{\tau} - \varepsilon\tau\sigma\frac{d\tilde{I}(\tilde{t})}{d\tilde{t}}(\tau\sigma\tilde{I}(\tilde{t})+1)^{-2} \\ + \frac{\varepsilon}{\tau\sigma\tilde{I}(\tilde{t})+1}\frac{d\tilde{C}(\tilde{t})}{d\tilde{t}}. $$

Using $\tilde{A}(\tilde{t}) = \hat{A}(t) + \frac{1-\tilde{C}(\tilde{t})}{\tau\sigma\tilde{I}(\tilde{t})+1}$ and setting $\varepsilon = 0$ which freezes $\tilde{t}$ at 0, we get the boundary layer model as

$$ \frac{d\hat{A}(t)}{dt} = -\left(\sigma\tilde{I}(0) + \frac{1}{\tau}\right)\hat{A}(t). \qquad (17) $$

We have the following result.

**Lemma 1.** Let A and C be the solution to (4)-(6), then for any fixed $T > 0$ there exists $\varepsilon_0$ such that $0 < \varepsilon < \varepsilon_0$ and

$$ C(t) - \check{C}(t) = O(\varepsilon) \qquad (18) $$

$$ A(t) - \frac{1 - \check{C}(t)}{\tau \sigma I(t) + 1} - \hat{A}(t/\varepsilon) = O(\varepsilon) \qquad (19) $$

hold uniformly for $t \in [0, T]$, where $\check{C}$ denotes the solution of equation

$$ \frac{d\check{C}}{dt} = -\alpha(I)\check{C} + \beta(I), \qquad (20) $$

where

$$ \alpha(I) = k_d \tau \frac{(\sigma I)^2}{\tau \sigma I + 1} + k_r \qquad (21) $$

$$ \beta(I) = k_d \tau \frac{(\sigma I)^2}{\tau \sigma I + 1}, \qquad (22) $$

with initial condition (14).

*Proof.* It is straightforward to see that the equilibrium point $\hat{A} = 0$ of the boundary layer model in (17) is asymptotically stable. By applying Thikhonov Theorem (see e.g. Theorem 11.1 in [10]), we get that for any fixed $T > 0$ there exists $\varepsilon_0$ such that $0 < \varepsilon < \varepsilon_0$ and system (11)-(14) has an unique solution $\tilde{A}(\tilde{t}, \varepsilon)$, $\tilde{C}(\tilde{t}, \varepsilon)$ on $[0, \varepsilon T]$, and

$$ \tilde{C}(\tilde{t}, \varepsilon) - \bar{C}(\tilde{t}) = O(\varepsilon) \qquad (23) $$

$$ \tilde{A}(\tilde{t}, \varepsilon) - \frac{1 - \bar{C}(\tilde{t})}{\tau \sigma \tilde{I}(\tilde{t}) + 1} - \hat{A}(\tilde{t}) = O(\varepsilon) \qquad (24) $$

hold uniformly for $\tilde{t} \in [0, \varepsilon T]$, where $\bar{C}$ is the solution to equation (16) with initial condition (14). Using the change of variable (10) we get that (18) and (19) hold for $t \in [0, T]$. This concludes the proof of Lemma 1. □
The net specific growth rate is given by the balance between photosynthesis and respiration $\mu(A, I) = k\sigma IA - R$ where the constant $R > 0$ denotes the respiration rate and $k > 0$ is a coefficient. This gives the following dynamics for the biomass $X$

$$ \frac{dX}{dt} = \mu(A, I)X. \qquad (25) $$

Now, as we are interested in the evolution of $C$ we expressed the function $\mu$ as a function of this latter state. Because $A = \frac{1-C}{\tau\sigma I+1}$, then the biomass dynamics writes

$$ \mu(A, I) = \gamma(I)C + \zeta(I) = \mu(C, I), \qquad (26) $$

where

$$ \gamma(I) = - \frac{k\sigma I}{\tau\sigma I + 1} \qquad (27) $$

$$ \zeta(I) = \frac{k\sigma I}{\tau\sigma I + 1} - R. \qquad (28) $$

In the rest of the paper we assume that the light signal $I$ is given by the function

$$ I(t) = -\frac{4}{T^2} (I_0 - 100) t^2 + \frac{4}{T} (I_0 - 100) t + 100, \qquad (29) $$

where $t \in [0, T]$ and $I_0$ is the maximal light intensity of the signal. The parameters used are given in Table 1.

<table><caption>Table 1: Parameter values used for the simulations.</caption><tr><td>I<sub>0</sub></td><td>1500</td><td>[µmol.m<sup>-2</sup>.s<sup>-1</sup>]</td></tr><tr><td>k<sub>r</sub></td><td>4.8e-4</td><td>[s<sup>-1</sup>]</td></tr><tr><td>k<sub>d</sub></td><td>2.99e-4</td><td>[-]</td></tr><tr><td>τ</td><td>1/0.1460</td><td>[s]</td></tr><tr><td>σ</td><td>0.0029</td><td>[m<sup>2</sup>.(µmol)<sup>-1</sup>]</td></tr><tr><td>k</td><td>3.6467e-4</td><td>[-]</td></tr><tr><td>R</td><td>0.05</td><td>[d<sup>-1</sup>]</td></tr></table>

## 3 PDE Modelling

The dynamics modelling of the RAB during the day can be represented by Partial Differential Equations (PDEs). Basically, it consists to decouple the system in two: an enlightened part and a dark part. For each part we consider two PDEs: a PDE for the photoinhibition dynamics and a PDE for the biomass dynamics. Therefore, we have a 4 × 4 system of transport-reaction PDEs linked at the boundary by a linear relationship. We have the following dynamics for the photoinhibition rate and the biomass, for the enlightened portion of the RAB,

$$ \partial_t C_1 + u(t) \partial_x C_1 = -\alpha(I(t)) C_1 + \beta(I(t)) \qquad (30) $$

$$ \partial_t X_1 + u(t) \partial_x X_1 = \gamma(I(t)) X_1 C_1 + \zeta(I(t)) X_1, \qquad (31) $$

respectively and $(t, x) \in \mathbb{R}^+ \times (0, \ell)$, where $\ell$ is the length of the belt exposed to light (see Figure 1). The functions $\alpha, \beta, \gamma$, and $\zeta$ are defined in (21), (22), (27), and (28) respectively. The light signal $I(t)$ depends on the time because in outdoor cultivation photons come from the sun.
In the above equations $u$ denotes the velocity of the RAB. We suppose that the velocity can vary during the process operation and that it can be adapted to the state of the system together with the light intensity. For the portion of the RAB in darkness we have

$$ \partial_t C_2 + u(t) \partial_x C_2 = -k_r C_2 \quad (32) $$

$$ \partial_t X_2 + u(t) \partial_x X_2 = -R X_2, \quad (33) $$

where $(t, x) \in \mathbb{R}^+ \times (0, (N-1)\ell)$. Here, $N\ell$ is the total length of the conveyor $(1/N$: ratio of belt which is enligthen, called the belt folding ratio). The space derivative corresponds to the transport part and the right-hand sides of (30)–(33) correspond to the reaction terms. Since we are working in bounded space domains $(0, \ell)$ and $(0, (N-1)\ell)$ we need a boundary condition for each equation to get a well-posed problem. It is given by

$$ C_1(t, 0) = C_2(t, (N-1)\ell), \quad C_2(t, 0) = C_1(t, \ell), \quad (34) $$

$$ X_1(t, 0) = X_2(t, (N-1)\ell), \quad X_2(t, 0) = X_1(t, \ell). \quad (35) $$

The above boundary conditions may be read as “what leaves the enlightened part enters the dark part and vice versa”. Finally, initial conditions shall be specified

$$ C_i(0, \cdot) = C_i^0 \in E \quad (36) $$

$$ X_i(0, \cdot) = X_i^0 \in L^\infty((0, 1); \mathbb{R}^+), \quad (37) $$

for $i = 1, 2$. We are interested to maximize the production of biomass between the intial time instant and $t = T$ and we get the productivity by ponderating it by the production duration (productivity per unit of time). Therefore, we state the average areal productivity of our process as

$$ \bar{\mathcal{P}}(u, N; T, \ell, I) = \frac{1}{T} \left[ \int_{0}^{(N-1)\ell} (X_2(T, x) - X_2^0(x)) dx + \int_{0}^{\ell} (X_1(T, x) - X_1^0(x)) dx \right]. \quad (38) $$

**Proposition 1.** The average areal productivity over the time interval $[0, T]$ can be rewritten as

$$ \bar{\mathcal{P}}(u, N; T, \ell, I) = \frac{1}{T\ell} \int_0^T \int_0^\ell (\gamma(I(s)) (X_1 C_1)(s, x) + \zeta(I(s)) X_1(s, x)) dx ds \\ - \frac{R}{T\ell} \int_0^T \int_0^{(N-1)\ell} X_2(s, x) dx ds. \quad (39) $$

*Proof.* Since $\int_0^T \partial_t X_i(t, x) dt = X_i(T, x) - X_i(0, x)$, using (31) and (33), the quantity $\bar{\mathcal{P}}$ in (38) can be rewritten as

$$ \bar{\mathcal{P}} = -\frac{1}{T} \left[ \int_0^T u(t) \int_0^{(N-1)\ell} \partial_x X_2(t, x) dx dt \int_0^T u(t) \int_0^\ell \partial_x X_1(t, x) dt dx \right] \\ + \frac{1}{T} \left[ -R \int_0^T \int_0^{(N-1)\ell} X_2(t, x) dx dt + \int_0^T \int_0^\ell \mu(C_1(t, x), I(t)) X_1(t, x) dx dt \right]. \quad (40) $$

Using boundary conditions (72) and the expression of $\mu$ given in (26) one gets (39). $\square$
# 4 Gradient-Based Optimization of the RAB

## 4.1 Statement of the Optimization Problem

In this section, the objective is to take benefit of the model to optimize process performances and especially productivity as given by the equation in (39).

We use a practical method for an efficient identification of the optimal conditions. More precisely, we will find these conditions with an adjoint-based optimization algorithm. We aim at optimizing both the velocity $u$ of the RAB and the shape parameter $N$ which relies on the belt folding ratio. The mixed parameter optimization and optimal control problem for the given time horizon $[0, T]$, given light signal $I$, and a given $\ell$, reads

$$ \max_{(u,N) \in \mathcal{U} \times D} \bar{P}(u, N; T, \ell, I), \quad (41) $$

where $D = [2, +\infty)$ and $\mathcal{U}$ is a precompact set of functions, which will be specified later. The main difficulty in optimizing the productivity $\bar{P}$ is the presence of the parameter $N$ in the bound of the second integral in (39). To get ride of this issue let us consider the following functional

$$ \begin{aligned} \mathcal{P}(u, N; T, \ell, I) = & \frac{1}{T} \int_0^T \int_0^1 (\gamma(I(s)) (\bar{X}_1 \bar{C}_1) (s, x) + \zeta(I(s)) \bar{X}_1 (s, x)) \, dxds \\ & - \frac{R(N-1)}{T} \int_0^T \int_0^1 \bar{X}_2(s, x) \, dxds, \end{aligned} \quad (42) $$

where $\bar{X}_1, \bar{X}_2, \bar{C}_1$, and $\bar{C}_2$ satisfy the following system

$$ \partial_t \bar{C}_1 + \frac{u(t)}{\ell} \partial_x \bar{C}_1 = -\alpha(I(t)) \bar{C}_1 + \beta(I(t)) \quad (43) $$

$$ \partial_t \bar{X}_1 + \frac{u(t)}{\ell} \partial_x \bar{X}_1 = \gamma(I(t)) \bar{X}_1 \bar{C}_1 + \zeta(I(t)) \bar{X}_1 \quad (44) $$

$$ \partial_t \bar{C}_2 + \frac{u(t)}{(N-1)\ell} \partial_x \bar{C}_2 = -k_r \bar{C}_2 \quad (45) $$

$$ \partial_t \bar{X}_2 + \frac{u(t)}{(N-1)\ell} \partial_x \bar{X}_2 = -R \bar{X}_2, \quad (46) $$

in $(t, x) \in (0, \infty) \times (0, 1)$ subject to boundary conditions

$$ \bar{C}_1(t, 0) = \bar{C}_2(t, 1), \quad \bar{C}_2(t, 0) = \bar{C}_1(t, 1), \quad (47) $$

$$ \bar{X}_1(t, 0) = \bar{X}_2(t, 1), \quad \bar{X}_2(t, 0) = \bar{X}_1(t, 1), \quad (48) $$

and initial conditions

$$ \bar{C}_i(0, \cdot) = \bar{C}_i^0 \quad (49) $$

$$ \bar{X}_i(0, \cdot) = \bar{X}_i^0, \quad (50) $$

where

$$ \bar{C}_i^0(x) = C_i^0(\ell x), \quad (51) $$

$$ \bar{X}_i^0(x) = X_i^0((N-1)\ell x). \quad (52) $$
And let us consider the new objective

$$ \max_{(u,N) \in \mathcal{U} \times D} \mathcal{P}(u, N; T, \ell, I). \qquad (53) $$

We have the following lemma.

**Lemma 2.** The mixed parameter optimization and optimal control for the nonlinear infinite dimensional problems (41) and (53) are equivalent.

*Proof.* As equation (52) leave to believe we will use a change of variables, namely

$$ \tilde{x}_1 = \frac{x_1}{\ell}, \qquad \tilde{x}_2 = \frac{x_2}{(N-1)\ell}, \qquad (54) $$

$$ \bar{C}_i(t, \tilde{x}_i) = C_i(t, x_i), \qquad \bar{X}_i(t, \tilde{x}_i) = X_i(t, x_i), \qquad (55) $$

for $i = 1, 2$. One has

$$ \partial_t C_i(t, x_i) = \partial_t \bar{C}_i(t, \tilde{x}_i), \quad i = 1, 2, \qquad (56) $$

$$ \partial_{x_1} C_1(t, x_1) = \frac{1}{\ell} \partial_{\tilde{x}_1} \bar{C}_1(t, \tilde{x}_1), \qquad (57) $$

$$ \partial_{x_2} C_2(t, x_2) = \frac{1}{(N-1)\ell} \partial_{\tilde{x}_2} \bar{C}_2(t, \tilde{x}_2), \qquad (58) $$

with

$$ \left\{ \begin{array}{l} C_i(t, 0) = \bar{C}_i(t, 0), \quad i = 1, 2, \\ C_1(t, \ell) = \bar{C}_1(t, 1), \\ C_2(t, (N-1)\ell) = \bar{C}_2(t, 1). \end{array} \right. \qquad (59) $$

The same holds for $\bar{X}_1$ and $\bar{X}_2$. Therefore, system (30)–(33) becomes (43)–(46) through relationships (54)–(55). Then, it can be verified that the productivity in (39) is given by (42) for the dimensionless system (43)–(50). It concludes the proof of Lemma 2. $\square$

While the optimization problem (41) needs some elaborate analysis, the optimization problem (53) is rather simple to solve in comparison.

Let us note that system (43)–(50) is well-posed. More precisely, for the existence and uniqueness of solution to system (43)–(50) we will consider solution in a broad sense. The two following definitions give the precise meaning of this kind of solution.

**Definition 1.** The characteristic for the dynamics of $\bar{C}_1$ and $\bar{X}_1$ (resp. $\bar{C}_2$ and $\bar{X}_2$) is an absolutely continuous function $s \mapsto \zeta_1(s; t^*, x^*)$ (resp. $s \mapsto \zeta_2(s; t^*, x^*)$) which satisfies $\zeta_1(t^*; t^*, x^*) = x^*$ (resp. $\zeta_2(t^*; t^*, x^*) = x^*$) and the ordinary differential equation

$$ \frac{d}{ds} \zeta_1(s; t^*, x^*) = \frac{u(s)}{\ell}, \qquad (60) $$

(resp. $\frac{d}{ds} \zeta_2(s; t^*, x^*) = \frac{u(s)}{(N-1)\ell}$) almost everywhere on the domain where $\zeta_1(\cdot; t^*, x^*)$ (resp. $\zeta_2(\cdot; t^*, x^*)$) is defined.

**Definition 2.** Let $\bar{C}_1^0, \bar{X}_1^0, \bar{C}_2^0, \bar{X}_2^0 \in L^\infty((0, 1); \mathbb{R}^+)$ be given. A solution to the Cauchy problem associated with (43)–(50) is a vector function with components $\bar{C}_1, \bar{X}_1, \bar{C}_2, \bar{X}_2$ such that,
for every $(t^*, x^*) \in \mathbb{R}^+ \times (0,1)$, the former components are satisfying the ordinary differential equations

$$ \frac{d}{dt} \bar{C}_1(t, \zeta_1(t; t^*, x^*)) = -\alpha(I(t))\bar{C}_1(t, \zeta_1(t; t^*, x^*)) + \beta(I(t)) \quad (61) $$

$$ \frac{d}{dt} \bar{X}_1(t, \zeta_1(t; t^*, x^*)) = \zeta(I(t))\bar{X}_1(t, \zeta_1(t; t^*, x^*)) + \gamma(I(t))\bar{X}_1(t, \zeta_1(t; t^*, x^*)) \\ \times \bar{C}_1(t, \zeta_1(t; t^*, x^*)) \quad (62) $$

$$ \frac{d}{dt} \bar{C}_2(t, \zeta_2(t; t^*, x^*)) = -k_r \bar{C}_2(t, \zeta_2(t; t^*, x^*)) \quad (63) $$

$$ \frac{d}{dt} \bar{X}_2(t, \zeta_2(t; t^*, x^*)) = -R \bar{X}_2(t, \zeta_2(t; t^*, x^*)) , \quad (64) $$

for every $t \ge t^*$, such that $\zeta_1(t; t^*, x^*)$, $\zeta_2(t; t^*, x^*) \in (0,1)$. Moreover, the vector function components satisfy the initial conditions (50) for $t = 0$ together with the boundary conditions (47)-(48).

Then, the existence and uniqueness of solution accordingly to Definitions 1 and 2 to system (43)-(50) is guaranteed by the following proposition.

**Proposition 2.** Let $u \in C([0, \infty); \mathbb{R}^+ \setminus \{0\})$, $I \in C([0, \infty]; \mathbb{R}^+)$ a periodic function of period $T$, $\bar{C}_i^0 \in L^\infty((0,1); \mathbb{R}^+ \setminus \{0\})$ and $\bar{X}_i^0 \in L^\infty((0,1); \mathbb{R}^+ \setminus \{0\})$ for $i = 1, 2$. Moreover, we assume that $\bar{C}_i^0, \bar{X}_i^0$ are bounded. Then, there exists an unique bounded solution $(\bar{C}_1, \bar{C}_2, \bar{X}_1, \bar{X}_2)^\top \in C^0([0, \infty); L^\infty((0,1); \mathbb{R}^+))$ to system (43)-(50).

*Proof.* The important point to remark concerning system (43)-(50) is that it is a cascade of $2 \times 2$ linear hyperbolic systems. Indeed, the dynamics of $\bar{C}_i$, $i = 1, 2$, are completely independent of the dynamics of $\bar{X}_i$, $i = 1, 2$, while the dynamics of $\bar{X}_1$ depends on $\bar{C}_1$. Then, using Definitions 1 and 2 we can easily construct bounded positive solutions in the space of $L^\infty$ functions for $\bar{C}_i$, $i = 1, 2$, globally in time. Then, using $\bar{C}_1$ as an exogeneous input for equation (44) we can apply again the characteristics method and obtain that there exist global bounded positive solutions in time in $L^\infty$ for $\bar{X}_i$, $i = 1, 2$. The uniqueness of the solutions follows directly from Definition 2. This concludes the proof of Proposition 2. $\square$

From now on and for the sake of simplicity we will drop the bar in the notation of $\bar{C}_i$ and $\bar{X}_i$, $i = 1, 2$, and the dependence in the time variable whenever it is unnecessary.

## 4.2 Existence of an optimal control

The velocity is supposed to be time-varying. Thus, the optimization problem consists in finding an optimal point in an infinite-dimensional space. In our case this functional space is $\mathcal{U} := \{C^1(0, T) | u_{\min} \le u(t) \le u_{\max}\}$. This control space ensures that, for each fixed $N$, there exists an optimal control for our problem. More precisely, we have the following theorem.

**Theorem 1.** Let $N$ be fixed. There exists an optimal control $u^*$ in $\mathcal{U}$ maximizing the function $\mathcal{J}$.

*Proof.* We are considering the following control space

$$ \mathcal{U} := \{u \in C^1(0, T) : 0 < u_{\min} \le u(t) \le u_{\max}\} . $$

System satisfied by $X_1$ and $X_2$ can be written as

$$ \partial_t X_1 + \frac{u(t)}{\ell} \partial_x X_1 = (a(t) + b[u](t, x)) X_1 \quad (65) $$

$$ \partial_t X_2 + \frac{u(t)}{(N-1)\ell} \partial_x X_2 = -RX_2, \quad (66) $$
with $x \in (0, 1)$, where

$$a(t) = \zeta(I(t))) \qquad (67)$$

$$b[u](t, x) = \gamma(I(t))C_1(t, x), \qquad (68)$$

and the boundary conditions

$$\begin{pmatrix} X_1(0,t) \\ X_2(0,t) \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} X_1(1,t) \\ X_2(1,t) \end{pmatrix}. \qquad (69)$$

Without loss of generality we consider the following system

$$\partial_t y_1 + u(t) \partial_x y_1 = (a(t) + b[u](t, x)) y_1 \qquad (70)$$

$$\partial_t y_2 + u(t) \partial_x y_2 = cy_2, \qquad (71)$$

with $x \in (0, 1)$ and the boundary conditions

$$\begin{pmatrix} y_1(0,t) \\ y_2(0,t) \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} y_1(1,t) \\ y_2(1,t) \end{pmatrix}. \qquad (72)$$

where $0 < a, c < 0$, and $-N \le b[\cdot](t, x) \le 0$ for all $(t, x) \in \mathbb{R}^+ \times (0, 1)$. Moreover, $b[u] \in C^1([0, T]; X)$ for all $u \in \mathcal{U}$. Let us consider the objective function given by

$$\mathcal{J} = \int_{0}^{T} \int_{0}^{1} L(y_{1}(t, x)) \, dxdt. \qquad (73)$$

Let us denote $L^2((0, 1); \mathbb{R}^2)$ by $X$ equipped with the usual scalar product

$$\langle (\phi_1, \phi_2)^\top, (\psi_1, \psi_2)^\top \rangle = \int_0^1 (\phi_1 \psi_1 + \phi_2 \psi_2) \, dx. \qquad (74)$$

For every $t$, let us consider the following linear operator $A(t) : D(A(t)) \to X$

$$D(A(t)) = \left\{ \phi \in H^1(0,1); \begin{pmatrix} \phi_1(0) \\ \phi_2(0) \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} \phi_1(1) \\ \phi_2(1) \end{pmatrix} \right\} \\ (A(t)\phi)(x) = u(t)\phi'(x) - \begin{pmatrix} b[u](t,x) - a(t) & 0 \\ 0 & -c \end{pmatrix} \phi(x). \qquad (75)$$

and $D(A(t)) = D(A(0))$. Let us denote $\mathcal{C}$ the set of operators $A(t)$. It is clear that the mapping $t \mapsto A(t)$ is continuous in the uniform operator topology. Then, using Theorem 5.1 in [20] system (70)–(72) has a unique solution in $C^1([0, T]; H^1((0, 1); \mathbb{R}^2))$. Let $(u_\beta)$, $0 \le \beta \le 1$ be a sequence in $\mathcal{U}$. Since, $\mathcal{U}$ is compact (by the Arzelà-Ascoli Theorem) there exists a sequence $(u_\beta')$ of $(u_\beta)$ so that $u_{\beta'}(t) \to u^*(t) \in \mathcal{U}$ uniformly for $0 \le t \le T$ as $\beta' \to 0$. As wrote before, there exists a classical solution to the initial value problem

$$\begin{cases} \dot{y}(t) + A_{\beta}(t)y(t) = 0, \\ y(0) = y^0. \end{cases} \qquad (76)$$

By Theorem 3.1 in [2], the solution reads $y_\beta(t) = U_\beta(t, 0)y^0$ where $U_\beta(t, 0)$ is the solution operator associated with (76), and $U_\beta(t, 0)y^0$ is continuously differentiable in $t$. Thus, $U(t, 0) = \lim_{\beta \to 0} U_\beta(t, 0)$ is the strong limit of $C^1$ evolution operators. Therefore, the set $\mathcal{C}$ of generators is strongly compact in the sense of [23]. Since $A(t)$ are single valued and closed and $D(A(t)) = D(A(0))$, then by applying Theorem 3.1 and Theorem 5.1 of [23] we get that there exists an optimal control $u^*$ in $\mathcal{U}$ maximizing the function $\mathcal{J}$. This concludes the proof of Theorem 1. $\square$
Let us note that our system of RAB satisfies the assumption of Theorem 1. Indeed, I is continuous. And by standard argument, one has that $C_1 \in C^1([0, T]; X)$.

The mixed optimization problem is generally numerically solved with a gradient descent algorithm. The derivation of the gradient of the cost function with respect to the velocity and the folding of the RAB cannot be derived as for a vector space function. The so-called “adjoint-method” may be used [12]. This method is based on an adjoint-system which has to be solved backward in time. Notice the adjoint variables play the role of Lagrange multiplier associated with the optimization constraints. This method is strongly linked to the calculus of variations theory. The interest reader is reported to the introductory monograph [13]. Finally, the continuous gradient is expressed as a function of the adjoint-state and nominal state solutions.

From now on we rewrite the velocity $u$ as the sum of a positive constant term and an unknown function:

$$u(t) = u_0 + \delta(t), \qquad (77)$$

where $u_0 \in \mathbb{R}^+$ and $\delta \in C^1((0, T); \mathbb{R})$ such that $u \in C^1((0, T); \mathbb{R}^+)$. The term $u_0$ turns out to be a key control which is more efficiently managed when explicit in (77) (see Remark 2 below). Besides, $u_0$ imposes the rotation direction.

**Remark 1.** The velocity is assumed to be continuously differentiable to ensure that the optimal control problem as a solution. But, in practice it would be continuously differentiable. Indeed, the velocity is induced by an engine and the inertia of the system leads to an equation of the form $\dot{m} = F$, where the right-hand side term $F$ represents the balance between engine couple and friction terms, both terms being continuous.

## 4.3 Adjoint-Based Gradient Derivation

In order to consider a minimization problem we will consider the cost function in (42) times minus one i.e.

$$\mathcal{J}(u_0, \delta, N) = -\mathcal{P}. \qquad (78)$$

We have the following result.

**Proposition 3.** The gradient of the cost function, stated in (78), with respect to $U = (u_0, \delta, N)^\top$ is given by

$$\nabla_{u_0} \mathcal{J}(U) = -\frac{1}{T\ell} \int_0^T \int_0^1 \left( p_1 \partial_x C_1 + q_1 \partial_x X_1 + \frac{p_2 \partial_x C_2}{(N-1)} + \frac{q_2 \partial_x X_2}{(N-1)} \right) dxds \qquad (79)$$

$$\nabla_{\delta} \mathcal{J}(U) = -\frac{1}{T\ell} \int_{0}^{1} \left( p_{1} \partial_{x} C_{1} + q_{1} \partial_{x} X_{1} + \frac{p_{2} \partial_{x} C_{2}}{(N-1)} + \frac{q_{2} \partial_{x} X_{2}}{(N-1)} \right) dx \qquad (80)$$

$$\nabla_N \mathcal{J}(U) = \frac{u}{(N-1)^2 T \ell} \int_0^T \int_0^1 (p_2 \partial_x C_2 + q_2 \partial_x X_2) dxds + \frac{R}{T} \int_0^T \int_0^1 X_2 dxds, \qquad (81)$$

where $(p_1, p_2, q_1, q_2)^\top$ satisfy the following system

$$-\partial_t p_1 - \frac{u}{\ell} \partial_x p_1 = -\alpha(I) p_1 + \gamma(I) (q_1 - 1) X_1 \qquad (82)$$

$$-\partial_t q_1 - \frac{u}{\ell} \partial_x q_1 = (\gamma(I)C_1 + \zeta(I))(q_1 - 1) \qquad (83)$$

$$-\partial_t p_2 - \frac{u}{(N-1)\ell} \partial_x p_2 = -k_r p_2 \qquad (84)$$

$$-\partial_t q_2 - \frac{u}{(N-1)\ell} \partial_x q_2 = -Rq_2 + R(N-1), \qquad (85)$$
with the boundary conditions

$$p_1(t, 1) = \frac{p_2(t, 0)}{N-1}, \quad p_2(t, 1) = (N-1)p_1(t, 0), \qquad (86)$$

$$q_1(t, 1) = \frac{q_2(t, 0)}{N-1}, \quad q_2(t, 1) = (N-1)q_1(t, 0), \qquad (87)$$

and the initial conditions

$$p_1(T, \cdot) = p_2(T, \cdot) = q_1(T, \cdot) = q_2(T, \cdot) = 0. \qquad (88)$$

*Proof.* For $U = (u_0, \delta, N)^\top \in \mathbb{R}^+ \times C^1((0, T); \mathbb{R}^+) \times \mathbb{R}^+$, we denote by $C_i$ and $X_i$, $i = 1, 2$, the corresponding solution to system (43)–(48). For a perturbation $\tilde{U} = (\tilde{u}_0, \tilde{\delta}, \tilde{N}) \in \mathbb{R}^+ \times C^1((0, T); \mathbb{R}^+) \times \mathbb{R}^+$ we denote by $\tilde{C}_i$ and $\tilde{X}_i$, $i = 1, 2$, the solution to system (43)–(48) for the variables $U + \varepsilon \tilde{U}$. Denoting

$$\begin{aligned} \hat{C}_i &= \lim_{\varepsilon \to 0} \frac{\tilde{C}_i(t, x) - C_i(t, x)}{\varepsilon} \\ \hat{X}_i &= \lim_{\varepsilon \to 0} \frac{\tilde{X}_i(t, x) - X_i(t, x)}{\varepsilon} \\ \hat{\mathcal{J}} &= \lim_{\varepsilon \to 0} \frac{\mathcal{J}(U + \varepsilon \tilde{U}) - \mathcal{J}(U)}{\varepsilon}, \end{aligned}$$

we have

$$\hat{\mathcal{J}} = -\frac{1}{T} \int_{0}^{T} \int_{0}^{1} \left( \gamma(I) \hat{X}_{1}C_{1} + \gamma(I) X_{1}\hat{C}_{1} + \zeta(I) \hat{X}_{1} \right) dxds + \frac{R(N-1)}{T} \int_{0}^{T} \int_{0}^{1} \hat{X}_{2}dxds \\ + \frac{R\tilde{N}}{T} \int_{0}^{T} \int_{0}^{1} X_{2}dxds.$$

The variables $\hat{C}_i$ and $\hat{X}_i$ satisfy the so-called tangent model,

$$\begin{aligned} \partial_t \hat{C}_1 + \frac{u_0 + \delta}{\ell} \partial_x \hat{C}_1 &= -\frac{\tilde{u}_0 + \tilde{\delta}}{\ell} \partial_x C_1 - \alpha(I) \hat{C}_1 \\ \partial_t \hat{X}_1 + \frac{u_0 + \delta}{\ell} \partial_x \hat{X}_1 &= -\frac{\tilde{u}_0 + \tilde{\delta}}{\ell} \partial_x X_1 + \gamma(I) \hat{X}_1 C_1 + \gamma(I) X_1 \hat{C}_1 + \zeta(I) \hat{X}_1 \\ \partial_t \hat{C}_2 + \frac{u_0 + \delta}{(N-1)\ell} \partial_x \hat{C}_2 &= -\frac{\tilde{u}_0 + \tilde{\delta}}{(N-1)\ell} \partial_x X C_2 + \frac{\tilde{N}u}{(N-1)^2\ell} \partial_x C_2 - k_r \hat{C}_2 \\ \partial_t \hat{X}_2 + \frac{u_0 + \delta}{(N-1)\ell} \partial_x \hat{X}_2 &= -\frac{\tilde{u}_0 + \tilde{\delta}}{(N-1)\ell} \partial_x X_2 + \frac{\tilde{N}u_0 + \delta}{(N-1)^2\ell} \partial_x X_2 - R \hat{X}_2, \end{aligned}$$

where $\tilde{u}(t) = \tilde{u}_0 + \tilde{\delta}(t)$, with the associated boundary conditions

$$\begin{aligned} \hat{C}_1(t, 0) &= \hat{C}_2(t, 1), \quad \hat{X}_1(t, 0) = \hat{X}_2(t, 1), \\ \hat{C}_1(t, 1) &= \hat{C}_2(t, 0), \quad \hat{X}_1(t, 1) = \hat{X}_2(t, 0), \end{aligned}$$

and the initial conditions

$$\hat{C}_1(0, \cdot) = \hat{C}_2(0, \cdot) = \hat{X}_1(0, \cdot) = \hat{X}_2(0, \cdot) = 0.$$
Next we multiply the equations for $\hat{C}_1$, $\hat{C}_2$, $\hat{X}_1$, and $\hat{X}_2$ by some functions $p_1$, $p_2$, $q_1$, and $q_2$ respectively, representing the adjoint variables (see, for instance, [3] page 601, for the precise meaning of the adjoint system of an infinite dimensional system), then we integrate by part over $[0, T] \times (0, 1)$ and we obtain that $p_1$, $p_2$, $q_1$, and $q_2$ satisfy system (82)-(85) with the boundary conditions (86)-(87) and the initial conditions (88) (for $t = T$ since the system (82)-(85) has to be considered backward in time). Now, using that

$$ \lim_{\epsilon \to 0} \frac{\mathcal{J}(U + \epsilon \tilde{U}) - \mathcal{J}(U)}{\epsilon} = \nabla \mathcal{J}(U) \cdot \tilde{U}, $$

we get the formal expression for $\nabla \mathcal{J}(U)$ given by (79), (80), and (81). This concludes the proof of Proposition 3. □

**Remark 2.** Let us remark that for $t = T$, the initial conditions (88) for the adjoint system yields

$$ \int_{0}^{1} \left( \frac{p_{2}(T, x) \partial_{x} C_{2}(T, x)}{(N-1)} + \frac{q_{2}(T, x) \partial_{x} X_{2}(T, x)}{(N-1)} \right) dx \\ + \int_{0}^{1} (p_{1}(T, x) \partial_{x} C_{1}(T, x) + q_{1}(T, x) \partial_{x} X_{1}(T, x)) dx = 0. $$

Hence, without the constant term $u_0$ in (77) we would not be able to control the final velocity of the rotating biofilm. Besides, with constant initial conditions of the PDEs the same observation holds for the starting velocity.

## 4.4 Optimization Algorithm

In practice, we applied the following strategy to initialize and stop the gradient-based minimal search. In the numerical implementation of the descent algorithm we implement the gradient of $\mathcal{J}$ for the full-discretized system (discretization in time and space of the equations). There are two reasons guiding the full-discretization choice:

* It is a well-known fact that a space semi-discretization is not valid for a transport-reaction equation. Indeed, the derived equation get properties that the non-discretized equation does not have, such as the observability in any time while the non-discretized equation has only the property to be observable in a given minimum time $T > 0$. Therefore, the only possibility is a time semi-discretization. Hence, the resulting equation is a discrete-time system and not an ODE. Applying the maximum principle of Pontryagin (see [13]) can be difficult with the derived discrete time system.

* But most importantly, at the end we have to discretize the equation for a practical implementation and it is well-known that it is more difficult to get stability results on the full-discretization than in the semi-discretization framework.

The strategy to derive the gradient in the continuous case is essentially the same as for the discrete case. In realistic operation $u$ cannot reach any constant maximum value. Therefore, we impose an upper bound $u_{max} > 0$ on the velocity. Besides, we also exclude the velocity $u \equiv 0$ which corresponds to the uncontrolled case. Hence, we impose a lower bound $u_{min} > 0$. At each iteration, once the velocity is updated we proceed with a projection step of the computed velocity on the interval $[u_{min}, u_{max}]$, with the method proposed in [8]. The optimization program has been implemented in a C++ code. The Lax-Wendroff scheme (see e.g. [9]) is used to solve the
<figure><img src="image_4.png" /><figcaption>Figure 2: Growth rate (26) for the steady states to (20).</figcaption></figure>

transport-reaction PDEs (43)–(48). This scheme has been favored because it avoids a numerical diffusion in the solution and its accuracy is of second order in space and time. The number of points for the spatial discretization is 50 corresponding to $dx = 0.02$. The time step $dt$ is governed by a classical Courant-Friedrichs-Lewy (CFL) condition ensuring the stability of the discrete scheme. This condition dictates a direct relation between the spatial mesh and the time mesh. Hence, by refining the spatial discretization, the time discretization will be thinner, giving a higher computation time.

# 5 Numerical Results

## 5.1 Process Exploitation Hypotheses

Figure 2 shows that, with the parameters from Table 1, photoinhibition is present at steady states computed from the equation in (20), namely $C^* = \frac{\beta(I)}{\alpha(I)}$. This means that, for an algae permanently exposed to a PFD larger than 168 µmol.m⁻².s⁻¹, the growth rate µ decreases when the value of irradiance I increases.

This situation will be our reference case, and the gain in productivity with the RAB process will be computed in reference to it. The belt length exposed to light is fixed to $l = 1$ m. Let us recall that the light signal considered in this paper will be the one given in (29). Therefore, N will give the total length of the belt in meters (we recall that it is given by $Nl$).

The final time, to assess the productivity at the end of the day, is 28800 seconds (8 hours). The upper bound for the velocity is 0.3 m.s⁻¹, to reduce energy need for circulating the belt, and the lower bound to be 0.005 m.s⁻¹ to avoid belt immobilization, and thus ensure nutrient supply to the microalgae. The tolerance $\epsilon$ of the gradient algorithm is fixed to $10^{-4}$.

## 5.2 Optimization of the RAB

The simulations were run on a 2012 commercial laptop with 4 Gb of RAM and a 2.5 GHz Intel Core i5 processor. First, we optimize the process with homogeneous initial conditions, namely,

$$X_{1,2}^0 = 10 \ [\text{g/m}], \qquad (89)$$

$$C_{1,2}^0 = \kappa, \qquad (90)$$
where $\kappa \in [0, 1]$. The starting point for the optimization is $(u, \delta, N)^{\top} = (0.01, 0, 3)^{\top}$. The convergence results obtained when starting from different constant initial conditions $\kappa$ are displayed on Figure 5. It turns out that less than 15 iterations are necessary to converge towards an optimal solution. And it took around one hour and forty minutes CPU time for a run.

As expected the productivity is better with a lower initial photoinhibition rate, see Figures 5. The optimization process converges towards a $N$ close to 7.5 and this for each initial conditions of the PDEs tested, see Figures 5. Therefore, the results of the optimization procedure seem to be independent of the initial photoinhibition rate. The average light on the RAB can be assessed by dividing the averaged impinging light (1030 $\mu$mol.m$^{-2}$.s$^{-1}$) by this optimal $N_{opt}$, yielding 137 $\mu$mol.m$^{-2}$.s$^{-1}$. The growth rate at equilibrium for this averaged light is approximatively 2.79 d$^{-1}$ (Figure 2). It is worth remarking that this value deviates only from 1.06 % from the optimum growth rate at steady state ($\mu_{opt} \approx 2.82$ d$^{-1}$). In other words, the optimal rate $N_{opt}$ provides an average irradiance received by the cells, close to the optimal one obtained at equilibrium in a fixed biofilm permanently exposed to light ($N = 1$). Thus, the RAB provides an optimal way to dilute light along time.

For every initial conditions tested, the optimal velocity computed is constantly equal to $u_{max}$ (0.3 m.s$^{-1}$). We have carried out other simulations with a higher $u_{max}$ and we observe the same behavior. Besides, this observation is in accordance with the flashing effect already reported in the literature [7], demonstrating that the faster light varies, the higher the growth rate.

**Remark 3.** From an optimal control point of view the saturation of the velocity on its upper constraints looks like a bang-bang control. Unfortunately, there does not exist any result similar to the Pontryagin's Maximum Principle (PMP) for nonlinear systems in infinite dimension [4].

To better understand the role of the velocity we run a simulation where the only parameter to optimize is $N$, the velocity $u$ is fixed to 0.001 m.s$^{-1}$, the initial condition for the photoinhibition rate is 0.5. As shown by Figure 6 the optimization process stops after 4 iterations, and $N$ stays small compared to the case where the velocity is higher (cf Figure 5). This result should be understood as follows: for a low velocity, having a large $N$ means staying a long time in the dark and thus a large loss of biomass by respiration. Therefore, the velocity should be above a certain threshold. Once the velocity is large enough, the most impacting parameter is $N$ which mainly drives the optimal solution. Finally, the adaptation of the velocity does not provide a significant gain, and it does therefore not seem mandatory.

Moreover, in order to minimize the energy input in the RAB process, we computed the (time-varying) lowest velocity for which the deviations with respect to the optimal productivity is no more than 5%. The result is depicted on Figure 3. As we can see, the average value of this “energy saving velocity” is lower than the constant velocity computed from optimization process. It means that we may find some velocity controls for which the productivity remains close to the optimal productivity and which is less energy consuming. Adding a penalty term for the velocity in the cost function (42) would also result in decreasing the optimal velocity.

So far, we have supposed that the light signal $I$ is time dependent and exogenous to the system, since the photons come from the sun. Now, we led some numerical experiments with a spatial dependent light signal:

$$I(x) = I_0 + A \sin(\omega x). \qquad (91)$$

We have observed that the influence of the three parameters $I_0$, $A$, and $\omega$ is of second order compared to $N$ and $u$. For instance, for fixed $N$ and $u$ the RAB is less efficient with $I$ given by (91) than with $I \equiv I_0$.

Then, we have considered the situation of non-homogeneous initial conditions for the biomass
<figure><img src="image_5.png" /><figcaption>Figure 3: Velocity found when allowing a deviation of 5 % with respect to the optimal productivity for an initial condition $C_i^0 = 0.3$.</figcaption></figure>

density and the photoinhibition rate:

$$X_1^0(x) = 0.24 \cos (2x\pi) + 6.6 \qquad (92)$$

$$X_2^0(x) = 0.24 \cos (2(x + 1)\pi) + 6.6 \qquad (93)$$

$$C_1^0(x) = 0.25 \cos (2x\pi) + 0.25 \qquad (94)$$

$$C_2^0(x) = 0.25 \cos (2(x + 1)\pi) + 0.25. \qquad (95)$$

Here, the starting point for the optimization is $(u, \delta, N)^T = (0.1, 0, 3)^T$. As for the non-homogeneous initial conditions case the optimal velocity computed is constantly equal to the upper bound $u_{max}$. The evolution of the cost function and of the ratio $N$ are depicted on Figure 7, respectively. The ratio is smaller than in the homogeneous initial conditions case. This is due to the fact there is a proportion of the biofilm with a small density of biomass.

Finally, we have displayed the gain between the considered rotating microalgal biofilm process and a biofilm permanently exposed to the light ($N = 1$) on Figure 4 in function of different $I_0$ in (29). More precisely, we have computed the productivity for different initial photoinhibition rate conditions at each $I_0$. Then, we have considered the average value of these productivity. As we can see the gain is always positive meaning that the investigated technology gives better productivity than a biofilm permanently exposed to light. Besides, the gain increases as the maximum light intensity becomes large. For instance, for $I_0 = 2500 \mu mol.m^{-2}.s^{-1}$ the gain is greater than 100 %. Moreover, on the same figure, the average $N$ for the different initial photoinhibition rate conditions at each $I_0$ is depicted. As expected, greater is the maximum light intensity greater is the optimized $N$. Indeed, a high maximum light intensity implies much more inhibition for photosystem PSII than at a low maximum light intensity. Thus, much more time in the dark is needed by the damaged photosystems to recover, meaning a larger $N$ should be consider for the belt.

# 6 Conclusion and Perspectives

We have explored a new process to enhance microalgae productivity in the objective of reducing costs and environmental impact for the production of biofuel. We have demonstrated that the process acts as a light diluter in time, in order to reach an average light (per cycle) close to the optimal constant light for photosynthesis.
<figure><img src="image_6.png" /><figcaption>Figure 4: Gain [%] between a fixed biofilm permanently exposed to light and the RAB (solid line) and the average *N* (dashed line) in function of different maximum light intensity *I*₀ in (29).</figcaption></figure>

This study also highlights the power of an adjoint-based gradient method to identify the optimal working mode of a complex process described by PDE. This approach can be extended to a more general case, where algae are not fixed on a conveyor, but are advected in a classical raceway cultivation process. The full-discretization of the equations is suitable for the 3D shape optimization also because it becomes straightforward to account for the problem constraints (for instance, possitivity or mass conservation). Then, cost-functions depending on the energy consumption should be considered in the future. In order to treat this perspective a deeper analysis should be led to derive a relevant criterion to penalize the energy demand in the biological production.

In conclusion, this class of problems has many industrial perspectives with a variety of applications for renewable resources. Theoretical applications developed for these systems are still in their infancy and we are convinced that solving this problem is an important step for the improvement of microalgae based processes.

## Acknowledgments

The authors are grateful to the Inria Project Lab *Algae In Silico* and the Ademe project PHYTORECOLTE for their support.
# A Appendix: Numerical Results of the Optimization Process

In this section, we collect figures obtained from the optimization procedure for the examples of Section 5.

<figure><img src="image_7.png" /><figcaption>Figure 5: Evolution of the cost function $J$ in (78) (solid lines) and of the parameter $N$ (dashed lines) with respect to the iteration for different homogeneous initial conditions.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Evolution of the cost function $J$ in (78) (solid line) and of the parameter $N$ (dashed line) with respect to the iteration for homogeneous initial conditions when optimizing only $N$ for a velocity equal to 0.001 m.s$^{-1}$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Evolution of the cost function $J$ in (78) (solid line) and of the parameter $N$ (dashed line) with respect to the iteration for non-homogeneous initial conditions for the biomass density given in (92) and (93) and the photoinhibition rate given in (94) and (95).</figcaption></figure>

References

[1] A. M. Bayen, R. L. Raffard, and C. J. Tomlin. Adjoint-based control of a new Eule-rian network model of air traffic flow. *IEEE Transactions on Control Systems Technology*, 14(5):804–818, 2006.

[2] M. G. Crandall and A. Pazy. Nonlinear evolution equations in Banach spaces. *Israel Journal of Mathematics*, 11(1):57–94, 1972.

[3] R. F. Curtain and H. Zwart. *An Introduction to Infinite-Dimensional Linear Systems Theory*, volume 21. Springer Science & Business Media, 2012.

[4] Y. V. Egorov. Necessary conditions for optimal control in Banach spaces. *Matematicheskii Sbornik*, 64(1):79–101, 1964. (in Russian).

[5] M. Gross, W. Henry, C. Michael, and Z. Wen. Development of a rotating algal biofilm growth system for attached microalgae growth with in situ biomass harvest. *Bioresource Technology*, 150:195–201, 2013.

[6] B.-P. Han. Photosynthesis-irradiance response at physiological level: A mechanistic model. *Journal of Theoretical Biology*, 213(2):121–127, 2001.

[7] P. Hartmann, Q. Béchet, and O. Bernard. The effect of photosynthesis time scales on microalgae productivity. *Bioprocess and Biosystems Engineering*, 37(1):17–25, 2014.

[8] R. Herzog and K. Kunisch. Algorithms for PDE-constrained optimization. *GAMM-Mitteilungen*, 33(2):163–176, 2010.

[9] C. Hirsch. *Numerical Computation of Internal and External Flows: The Fundamentals of Computational Fluid Dynamics*. Butterworth-Heinemann, 2007.

[10] H. K. Khalil. *Nonlinear Systems*, volume 2. Prentice-Hall, New Jersey, 1996.

[11] L. Lardon, A. Helias, B. Sialve, J.-P. Steyer, and O. Bernard. Life-cycle assessment of biodiesel production from microalgae. *Environmental Science & Technology*, 43(17):6475–6481, 2009.
[12] F.-X. Le Dimet and O. Talagrand. Variational algorithms for analysis and assimilation of meteorological observations: Theoretical aspects. *Tellus A: Dynamic Meteorology and Oceanography*, 38(2):97-110, 1986.

[13] D. Liberzon. *Calculus of Variations and Optimal Control Theory: A Concise Introduction*. Princeton University Press, 2011.

[14] J.-L. Lions. *Optimal Control of Systems Governed by Partial Differential Equations*. Springer-Verlag, Berlin, 1971.

[15] T. M. Mata, A. A. Martins, S. K. Sikdar, and C. A. V. Costa. Sustainability considerations of biodiesel based on supply chain analysis. *Clean Technologies and Environmental Policy*, 13(5):655-671, 2011.

[16] E. Molina Grima, J. M. Fernández Sevilla, J. A. Sánchez Pérez, and F. García Camacho. A study on simultaneous photolimitation and photoinhibition in dense microalgal cultures taking into account incident and averaged irradiances. *Journal of Biotechnology*, 45(1):59-69, 1996.

[17] S. Nadarajah and A. Jameson. A comparison of the continuous and discrete adjoint approach to automatic aerodynamic optimization. In *38th Aerospace Sciences Meeting and Exhibit*, page 667, 2000.

[18] V. T. Nguyen, D. Georges, and G. Besançon. State and parameter estimation in 1-D hyperbolic PDEs based on an adjoint method. *Automatica*, 67:185-191, 2016.

[19] K. M. Passino. Biomimicry of bacterial foraging for distributed optimization and control. *IEEE Control Systems*, 22(3):52-67, 2002.

[20] Amnon Pazy. *Semigroups of Linear Operators and Applications to Partial Differential Equations*, volume 44. Springer Science & Business Media, 2012.

[21] C. Posten and S. F. Chen, editors. *Microalgae Biotechnology*. Springer, 2016.

[22] J. Reilly, S. Samaranayake, M. L. Delle Monache, W. Krichene, P. Goatin, and A. M. Bayen. Adjoint-based optimization on a network of discretized scalar conservation laws with applications to coordinated ramp metering. *Journal of Optimization Theory and Applications*, 167(2):733-760, 2015.

[23] M. Slemrod. Existence of optimal controls for control systems governed by nonlinear partial differential equations. *Annali della Scuola Normale Superiore di Pisa, Classe di Scienze*, (3-4):229-246, 1974.

[24] R. H. Wijffels and M. J. Barbosa. An outlook on microalgal biofuels. *Science*, 329(5993):796-799, 2010.

[25] D. B. Work and A. M. Bayen. Convex formulations of air traffic flow optimization problems. *Proceedings of the IEEE*, 96(12):2096-2112, 2008.

[26] X. Wu and J. C. Merchuk. A model integrating fluid dynamics in photosynthesis and photoinhibition processes. *Chemical Engineering Science*, 56(11):3527-3538, 2001.
informatics mathematics
Inria

RESEARCH CENTRE
SOPHIA ANTIPOLIS – MÉDITERRANÉE

2004 route des Lucioles - BP 93
06902 Sophia Antipolis Cedex

Publisher
Inria
Domaine de Voluceau - Rocquencourt
BP 105 - 78153 Le Chesnay Cedex
inria.fr

ISSN 0249-6399