# 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." Automatica, 2019, 105, pp.80-88. 10.1016/j.automatica.2019.02.043. hal-02422853

HAL Id: hal-02422853

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

Submitted on 23 Dec 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.
# Gradient-based optimization of a rotating algal biofilm process

Pierre-Olivier Lamare<sup>a,b</sup> Nina Aguillon<sup>c</sup> Jacques Sainte-Marie<sup>b</sup> Jérôme Grenier<sup>d</sup>
Hubert Bonnefond<sup>d</sup> Olivier Bernard<sup>a</sup>

## 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;

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

trate through the very turbid medium [2]. 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 [3,4]. 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

Email addresses:

pierre-olivier.lamare@mines-paristech.fr
(Pierre-Olivier Lamare), aguillon@ljll.math.upmc.fr
(Nina Aguillon), jacques.sainte-marie@inria.fr
(Jacques Sainte-Marie), jerome.grenier@student.ecp.fr
(Jérôme Grenier),
hubert.bonnefond@inalve-conception.com
(Hubert Bonnefond), olivier.bernard@inria.fr
(Olivier Bernard).

<footnote><sup>a</sup>Université Nice Côte d'Azur, Inria BIOCORE, BP93, 06902 Sophia-Antipolis Cedex, France</footnote>

<footnote><sup>b</sup>Inria ANGE, 2 rue Simone Iff, CS 42112, 75589 Paris Cedex 12, France</footnote>

<footnote><sup>c</sup>Université Pierre et Marie Curie, Laboratoire Jacques-Louis Lions, 75005 Paris, France</footnote>

<footnote><sup>d</sup>Sorbonne Universités, UPMC Université Paris 06, INSU-CNRS, Laboratoire d'Océanographie de Villefranche, 181 Chemin du Lazaret, 06230 Villefranche-sur-Mer, France</footnote>

<footnote>* This paper was not presented at any IFAC meeting.</footnote>
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 [6].

<figure><img src="image_1.png" /><figcaption>Fig. 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 [7]. Therefore, various control problems are tackled with this method as ramp metering [8], air-traffic control [9], inverse pressure design and drag minimization problems [10]. Very recently, the method has been applied to the state and parameter estimation for traffic flow and overland flow system in [11]. Finally, there exist alternative approaches for gradient-based optimization of biological systems such as [12].

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 [9,13]. In [13] 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 [9], beyond the application, is the treatment of the control constraints. In [9], 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 [14,15]). 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} \quad (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) \quad (2) $$

$$ \frac{dC(t)}{dt} = -k_r C(t) + k_d \sigma I(t) B(t), \quad (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} \quad (4) $$

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

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

$$ (A(0, \cdot), C(0, \cdot)) = (A^0, C^0) \in E, \quad (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 [16] and focus on the dynamics of the photoinhibition rate C. To lead this analysis we use the singular perturbation theory (see e.g. [17]). We define

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

$$K_2 = \varepsilon^{-1}k_d, \quad (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]. \quad (9)$$

Then, with the change of time

$$\tilde{t} = \varepsilon t, \quad (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}}{d\tilde{t}} = -\left(\sigma \tilde{I} + \frac{1}{\tau}\right) \tilde{A} + \frac{1 - \tilde{C}}{\tau} \quad (11)$$

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

with the initial conditions

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

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

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

$$\tilde{A} = \frac{1 - \tilde{C}}{\tau \sigma \tilde{I} + 1}. \quad (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}}{d\tilde{t}} = 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]. \quad (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}(\tilde{t})+1}$. We have

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

Therefore, we have

$$\frac{d\tilde{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}} \left(\tau \sigma \tilde{I}(\tilde{t}) + 1\right)^{-2} + \frac{\varepsilon}{\tau \sigma \tilde{I}(\tilde{t}) + 1} \frac{d\tilde{C}(\tilde{t})}{d\tilde{t}}.$$

Using $\tilde{A}(\tilde{t}) = \tilde{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\tilde{A}(t)}{dt} = -\left(\sigma \tilde{I}(0) + \frac{1}{\tau}\right) \tilde{A}(t). \quad (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) - \tilde{C}(t) = O(\varepsilon) \quad (18)$$

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

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

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

where

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

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

with initial condition (14).

**PROOF.** It is straightforward to see that the equilibrium point $\tilde{A} = 0$ of the boundary layer model in (17) is asymptotically stable. By applying Thikhonov Theorem (see e.g. Theorem 11.1 in [17]), 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) \quad (23) $$

$$ \tilde{A}(\tilde{t}, \varepsilon) - \frac{1 - \bar{C}(\tilde{t})}{\tau \sigma \tilde{I}(\tilde{t}) + 1} - \hat{A}(t) = O(\varepsilon) \quad (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 I A - 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. \quad (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}{\sigma\tau I + 1}$, then the biomass dynamics writes

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

where

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

$$ \zeta(I) = \frac{k\sigma I}{\tau\sigma I + 1} - R. \quad (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}{I} (I_0 - 100) t + 100, \quad (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). It consists in decoupling 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. Specifically, one has 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)) \quad (30) $$

$$ \partial_t X_1 + u(t) \partial_x X_1 = \gamma(I(t)) X_1 C_1 + \zeta(I(t)) X_1, \quad (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 enligthened, 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 X_i(0, \cdot) = X_i^0 \in L^\infty((0, 1); \mathbb{R}^+), \quad (36) $$

for $i = 1, 2$. We are interested to maximize the production of biomass between the initial 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{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 (37) $$
**Proposition 2** The average areal productivity over the time interval $[0, T]$ can be rewritten as

$$ \bar{P}(u, N; T, \ell, I) = -\frac{R}{T\ell} \int_0^T \int_0^{(N-1)\ell} X_2(s, x)dxds \\ + \frac{1}{T\ell} \int_0^T \int_0^\ell (\gamma(I(s))(X_1C_1)(s, x) \\ + \zeta(I(s))X_1(s, x))dxds. \quad (38) $$

**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{P}$ in (37) rewrites as

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

Using boundary conditions (35) and the expression of $\mu$ given in (26) one gets (38).

# 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 (38).

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 (40) $$

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 (38). To get ride of this issue we consider the functional

$$ \mathcal{P}(u, N; T, \ell, I) = -\frac{R(N-1)}{T} \int_0^T \int_0^1 \bar{X}_2(s, x) dxds $$

$$ + \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, \quad (41) $$

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 (42) $$

$$ \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{C}_1 \quad (43) $$

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

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

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 (46) $$

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

and initial conditions

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

where

$$ \bar{C}_i^0(x) = C_i^0(\ell x), \quad \bar{X}_i^0(x) = X_i^0((N-1)\ell x). \quad (49) $$

And let us consider the new objective

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

We have the following lemma.

**Lemma 3** The mixed parameter optimization and optimal control for the nonlinear infinite dimensional problems (40) and (50) are equivalent.

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

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

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

for $i = 1, 2$. By computing the partial derivatives w.r.t. $t$ and $x$ of $\bar{C}_i$ and $\bar{X}_i$, one has that system (30)–(33) becomes (42)–(45) through relationships (51)–(52). Then, it can be verified that the productivity in (38) is given by (41) for the dimensionless system (42)–(48). It concludes the proof of Lemma 3.

While the optimization problem (40) needs some elaborate analysis, the optimization problem (50) is rather
simple to solve in comparison. Let us mention that system (42)-(48) is well-posed. The proof is detailed in [18] and consists in using the characteristics method and noting that system (42)-(48) has a cascade structure: $C_1$ influences $X_1$, $C_2$ influences $X_2$. 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 4** Let $N$ be fixed. There exists an optimal control $u^*$ in $\mathcal{U}$ maximizing the function $\mathcal{J}$.

The proof of Theorem 4 is detailed in [18] and consists in writing the solution of our system with a solution operator. This solution operator depends on the velocity $u$. Using the compactness of $\mathcal{U}$ and by the Arzelà-Ascoli Theorem any sequence $(u_\beta)$ in $\mathcal{U}$ has a convergent sub-sequence $(u_{\beta'})$. By considering the associated sequence of solution operator $U_{\beta'}$, it may be shown that the set of generators $C$ is strongly compact in the sense of [19]. Then, by applying Theorem 3.1 and Theorem 5.1 of [19] we get that there exists an optimal control $u^*$ in $\mathcal{U}$ maximizing the function $\mathcal{J}$.

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 [20]. 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 [21]. 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), \quad (53)$$

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 (53) (see Remark 7 below). Besides, $u_0$ imposes the rotation direction.

**Remark 5** 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 (41) times minus one i.e.

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

We have the following result.

**Proposition 6** The gradient of the cost function, stated in (54), 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}{M} + \frac{q_2 \partial_x X_2}{M} \right) dx ds \quad (55)$$

$$\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}}{M} + \frac{q_{2} \partial_{x} X_{2}}{M} \right) dx \quad (56)$$

$$\nabla_N \mathcal{J}(U) = \frac{u}{M^2 T \ell} \int_0^T \int_0^1 (p_2 \partial_x C_2 + q_2 \partial_x X_2) dx ds + \frac{R}{T} \int_0^T \int_0^1 X_2 dx ds, \quad (57)$$

where $M = N - 1$ and $(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 \quad (58)$$

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

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

$$-\partial_t q_2 - \frac{u}{(N-1)\ell} \partial_x q_2 = -Rq_2 + R(N-1), \quad (61)$$

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), \quad (62)$$

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

and the initial conditions

$$p_1(T, \cdot) = p_2(T, \cdot) = q_1(T, \cdot) = q_2(T, \cdot) = 0. \quad (64)$$
**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 (42)–(47). 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 (42)–(47) 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

$$ 
\begin{aligned}
\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.
\end{aligned}
$$

The variables $\hat{C}_i$ and $\hat{X}_i$ satisfy the so-called tangent model. Due to the space limitation, we only give the expression of the system satisfied by $\hat{C}_1$ and $\hat{C}_2$

$$ 
\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{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 \\
&\quad + \frac{\tilde{N}u}{(N-1)^2\ell} \partial_x C_2 - k_r \hat{C}_2,
\end{aligned}
$$

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

$$ \hat{C}_1(t, 0) = \hat{C}_2(t, 1), \quad \hat{C}_1(t, 1) = \hat{C}_2(t, 0), \quad (65) $$

and the initial conditions

$$ \hat{C}_1(0, \cdot) = \hat{C}_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, [22] 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 (58)–(61) with the boundary conditions (62)–(63) and the initial conditions (64) (for $t = T$

since the system (58)–(61) has to be considered backward in time). Now, using that

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

we get the formal expression for $\nabla \mathcal{J}(U)$ given by (55), (56), and (57). This concludes the proof of Proposition 6.

**Remark 7** Let us remark that for $t = T$, the initial conditions (64) 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 (53) 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 [21]) 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 [23].

The optimization program has been implemented in a C++ code. The Lax-Wendroff scheme (see e.g. [24]) is used to solve the transport-reaction PDEs (42)–(47). 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. By refining the spatial discretization, the time discretization will be thinner, giving a higher computation time.

# 5 Numerical Results

## 5.1 Process Exploitation Hypotheses

The set of parameters given in Table 1 are such 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 \quad [\text{g/m}], \qquad (66)$$

$$C_{1,2}^{0} = \kappa, \qquad (67)$$

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 A.1. 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 A.1. The optimization process converges towards a $N$ close to 7.5 and this for each initial conditions of the PDEs tested, see Figures A.1. 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 µmol.m⁻².s⁻¹) by this optimal $N_{opt}$, yielding 137 µmol.m⁻².s⁻¹. The growth rate at equilibrium for this averaged light is approximately 2.79 d⁻¹. 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⁻¹). 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⁻¹). 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 [16], demonstrating that the faster light varies, the higher the growth rate.

**Remark 8** *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 [25].*

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⁻¹, the initial condition for the photoinhibition rate is 0.5. As shown by Figure A.2 the optimization process stops after 4 iterations, and $N$ stays small compared to the case where the velocity is higher (cf Figure A.1). 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 op-
timal productivity is no more than 5%. The result is depicted on Figure 2. 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 (41) 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). \quad (68)$$

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 (67) than with $I \equiv I_0$.

Then, we have considered the situation of non-homogeneous initial conditions for the biomass density and the photoinhibition rate:

$$X_1^0(x) = 0.24 \cos(2x\pi) + 6.6 \quad (69)$$

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

$$C_1^0(x) = 0.25 \cos(2x\pi) + 0.25 \quad (71)$$

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

Here, the starting point for the optimization is $(u, \delta, N)^\top = (0.1, 0, 3)^\top$. 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 A.3, 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 3 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\text{mol.m}^{-2}\text{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

<figure><img src="image_2.png" /><figcaption>Fig. 2. 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>

<figure><img src="image_3.png" /><figcaption>Fig. 3. 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_0$ in (29).</figcaption></figure>

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

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
<figure><img src="image_4.png" /><figcaption>Fig. A.1. Evolution of the cost function $J$ in (54) (solid lines) and of the parameter $N$ (dashed lines) with respect to the iteration for different homogeneous initial conditions.</figcaption></figure>

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.

## References

[1] R. H. Wijffels and M. J. Barbosa, "An outlook on microalgal biofuels," *Science*, vol. 329, no. 5993, pp. 796-799, 2010.

[2] C. Posten and S. F. Chen, Eds., *Microalgae Biotechnology*. Springer, 2016.

[3] L. Lardon, A. Helias, B. Sialve, J.-P. Steyer, and O. Bernard, "Life-cycle assessment of biodiesel production from microalgae," *Environmental Science & Technology*, vol. 43, no. 17, pp. 6475-6481, 2009.

[4] 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*, vol. 13, no. 5, pp. 655-671, 2011.

<figure><img src="image_5.png" /><figcaption>Fig. A.2. Evolution of the cost function $J$ in (54) (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⁻¹.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. A.3. Evolution of the cost function $J$ in (54) (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 (68) and (69) and the photoinhibition rate given in (70) and (71).</figcaption></figure>

[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*, vol. 150, pp. 195-201, 2013.

[6] 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*, vol. 45, no. 1, pp. 59-69, 1996.

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

[8] 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*, vol. 167, no. 2, pp. 733-760, 2015.

[9] A. M. Bayen, R. L. Raffard, and C. J. Tomlin, "Adjoint-based control of a new Eulerian network model of air traffic flow," *IEEE Transactions on Control Systems Technology*, vol. 14, no. 5, pp. 804-818, 2006.

[10] 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*, 2000, p. 667.

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

[12] K. M. Passino, "Biomimicry of bacterial foraging for distributed optimization and control," *IEEE Control Systems*, vol. 22, no. 3, pp. 52–67, 2002.

[13] D. B. Work and A. M. Bayen, "Convex formulations of air traffic flow optimization problems," *Proceedings of the IEEE*, vol. 96, no. 12, pp. 2096–2112, 2008.

[14] B.-P. Han, "Photosynthesis-irradiance response at physiological level: A mechanistic model," *Journal of Theoretical Biology*, vol. 213, no. 2, pp. 121–127, 2001.

[15] X. Wu and J. C. Merchuk, "A model integrating fluid dynamics in photosynthesis and photoinhibition processes," *Chemical Engineering Science*, vol. 56, no. 11, pp. 3527–3538, 2001.

[16] P. Hartmann, Q. Béchet, and O. Bernard, "The effect of photosynthesis time scales on microalgae productivity," *Bioprocess and Biosystems Engineering*, vol. 37, no. 1, pp. 17–25, 2014.

[17] H. K. Khalil, *Nonlinear Systems*. New Jersey: Prentice-Hall, 1996, vol. 2, no. 5.

[18] P.-O. Lamare, N. Aguillon, J. Sainte-Marie, J. Grenier, H. Bonnefond, and O. Bernard, "Gradient-based optimization of a rotating algal biofilm process," 2019, technical Report Available at: https://.

[19] 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*, no. 3–4, pp. 229–246, 1974.

[20] F.-X. Le Dimet and O. Talagrand, "Variational algorithms for analysis and assimilation of meteorological observations: Theoretical aspects," *Tellus A: Dynamic Meteorology and Oceanography*, vol. 38, no. 2, pp. 97–110, 1986.

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

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

[23] R. Herzog and K. Kunisch, "Algorithms for PDE-constrained optimization," *GAMM-Mitteilungen*, vol. 33, no. 2, pp. 163–176, 2010.

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

[25] Y. V. Egorov, "Necessary conditions for optimal control in Banach spaces," *Matematicheskii Sbornik*, vol. 64, no. 1, pp. 79–101, 1964, (in Russian).