\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{subfigure,graphicx}
\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2019 (2019), No. 84, pp. 1--26.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2019 Texas State University.}
\vspace{8mm}}
\begin{document}
\title[\hfilneg EJDE-2019/84\hfil Shape optimization of bioreactors]
{Shape optimization of spatial chemostat models}
\author[M. Crespo, B. Ivorra, A. M. Ramos, A. Rapaport \hfil EJDE-2019/84 \hfilneg]
{Maria Crespo, Benjamin Ivorra, \\
Angel Manuel Ramos, Alain Rapaport}
\address{Maria Crespo \newline
Departamento de Matematica Aplicada a la Ingenier\'ia Industrial,
ETSII - Universidad Polit\'ecnica de Madrid,
C Jos\'e Guti\'errez Abascal 2, 28006 Madrid, Spain}
\email{maria.crespo@upm.es, Phone +34 913377370}
\address{Benjamin Ivorra \newline
Instituto de Matem\'atica Interdisciplinar \&
Departamento de An\'alisis Matem\'atico y Matem\'atica Aplicada,
Universidad Complutense de Madrid,
Plaza de Ciencias, 3, 28040 Madrid, Spain}
\email{ivorra@ucm.es}
\address{ Angel Manuel Ramos\newline
Instituto de Matem\'atica Interdisciplinar \&
Departamento de An\'alisis Matem\'atico y Matem\'atica Aplicada,
Universidad Complutense de Madrid,
Plaza de Ciencias, 3, 28040 Madrid, Spain}
\email{angel@mat.ucm.es}
\address{Alain Rapaport \newline
MISTEA, Universit\'e Montpellier,
INRA, Montpellier SupAgro. 2 Place P.Viala,
34060 Montpellier, France}
\email{rapaport@supagro.inra.fr}
\dedicatory{Communicated by Jesus Ildefonso Diaz}
\thanks{Submitted February 1, 2019. Published July 11, 2019.}
\subjclass[2010]{35Q93, 35K40, 65K10, 92E20}
\keywords{Shape optimization; optimal design; continuous bioreactor;
\hfill\break\indent spatial chemostat model;
advection-diffusion-reaction; hybrid genetic algorithm}
\begin{abstract}
In this work, we study the shape optimization of a continuous bioreactor in which a
substrate is degraded by a microbial ecosystem in a nonhomogeneous environment.
The bioreactor considered here is a three-dimensional vertically oriented cylindrical tank.
The behavior of reactants is described with a spatial chemostat model based on an
Advection-Diffusion-Reaction system while the fluid flow is modeled using incompressible
Navier-Stokes equations. We consider that the reaction rate between biomass and nutrient
shows either monotonic or non-monotonic behavior. We tackle an optimization problem
which aims to minimize the considered total reactor volume, with an output concentration
(at stationary state) maintained below a desired threshold, by choosing a suitable
bioreactor shape. We propose a methodology to create three different discrete
parametrizations of the bioreactor geometry and obtain the optimized shapes with the help
of a Hybrid Genetic Algorithm. We show that the optimized reactors exhibit height much
larger than width and their exterior wall is concavely curved (the concavity at the upper
part of the exterior wall being more pronounced for non-monotonic functions).
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks
\section{Introduction}\label{sec:Introduction}
Shape optimization has been extensively exploited in design engineering
\cite{Jameson1988,ElSayed2005,Isebe2008,ivorra2013,ivorra2016,pironneau1984,Poloni2000},
particularly in aeronautical \cite{Campana2006,Raino1999,nejati2003} and automotive areas \cite{Muyl2004,Zuo2013}. Traditionally, finding the optimal geometry of a particular device is
based on a trial and error approach, in which, a number of prospective configurations is
simulated and the results are compared. An alternative strategy relays in performing the mathematical modeling of the process, carrying out numerical simulations and solving the
desired optimization problem with an appropriate optimization algorithm.
Taking into consideration the exponential growth of the available computing power, this second approach provides a powerful computational tool able to simulate and analyze the efficiency of different geometry configurations.
In this work, we tackle the shape optimization of a continuous bioreactor. A bioreactor is a
vessel in which microorganisms (e.g., bacteria), called biomass, are used to degrade a
considered diluted substrate (e.g., nitrate). A reactor in which substrate is continually added
and product continually removed is called continuous bioreactor. The influence of the
bioreactor design into the process efficiency has received considerable attention in the literature
\cite{Gajardo2008,Ramirez2017,Rodriguez2014}. Most of the works are developed by
experimentalists (see, e.g., \cite{Bitog2011,Carvalho2006}) and focus on specific
biological reactions occurring in continuous flow systems. Among the different reactor geometry configurations reviewed in literature, the most popular are flat-plate reactors
\cite{Sierra2008,Tredici1998}, torus-shaped reactors \cite{Pramparo2008} and tubular
reactors \cite{Richmond1993,Tredici1998}. In 2008, Xia et al. \cite{Xia2008} showed that
flow conditions (regarding mass transfer, shear stress, mixing, etc.) are strongly influenced by
the reactor geometry, particularly at large scales. Nevertheless, computational fluid dynamics
has not been commonly used to its full capacity to optimize reactor performance.
Of particular interest are the works developed by Ansoni et al. \cite{Ansoni2016} and
Coenen et al. \cite{Coenen2013}. In \cite{Ansoni2016}, the authors consider
a tubular reactor and model its hydrodynamics with 3D Navier-Stokes equations.
They look for the optimal design (configuration of the inlet and the outlet pipes) of a given bioreactor, so that the dispersion of the residence time and the shear flow are minimized.
These two concepts, related to fluid dynamics, are linked to the reactor effectiveness.
In \cite{Coenen2013}, the authors consider a cylindrical photobioreactor and model the
dynamics of the organic compound with an advection-reaction equation. They look for
the optimal geometry (radius and height) so that the reactant concentration at the outlet
of the bioreactor is minimized.
We aim to solve the following design optimization problem: given the input reactant
concentrations and the flow rate to be processed, what is the minimal reactor volume
(and its shape) so that a desired output reactant concentration is attained?
This problem has been modeled using ordinary differential equations (see, e.g.,
\cite{dochain2001,Harmand2005,Harmand2003,Hill1989,Luyben1982}) by approximating
the behavior of a tubular device with a bioreactor composed by $N$ well-mixed tanks in series.
Then, the aim of the optimization becomes to find what are the volumes of the $N$ tanks
such that the total volume of the whole process is minimal. However, these studies suffer of
two important drawbacks:
\begin{itemize}
\item While the proposed results are valid for small and medium sized systems,
they do not describe the diffusion phenomena or the impact of fluid motions
that may occur in larger tanks.
\item The dimensioning parameters were not considered; only the total volume
of the systems was optimized. However, in a real case, design parameters
such as the diameter or the height of any biological or chemical system will
influence its performance.
\end{itemize}
To overcome these drawbacks, we propose to couple hydrodynamics with biological
phenomena occurring in a diffusive bioreactor. To do so, we use a particular spatial
modeling based on Navier-Stokes equation (describing the fluid dynamics) together with
an Advection-Diffusion-Reaction system (describing the behavior of the reactants in
the bioreactor). We give a methodology to create three different discrete
parametrizations of the bioreactor geometry and obtain the optimized shapes with
the help of a Genetic Multi-layer Algorithm \text{(GMA)}, a global optimization
method based on the
hybridization of a Multi-layer Secant method \text{(MSA)}
\cite{redondo1,redondo2,redondo3}, finding suitable initial conditions for a global
optimization algorithm, with a given Genetic Algorithm \text{(GA)}
\cite{deb2001,goldberg1989,sharapov2007,weise2009}. This kind of hybridization has
been already tested for solving different optimization problems
\cite{gomez2011,ivorra2007}
and, according to several numerical experiments, it seems that it achieves
better results
(in terms of computational time and precision) than the \text{GA} used alone.
The optimization problem is solved for monotonic and non-monotonic growth rate
functions, in order to analyze the influence of the reaction into the optimal
reactor configuration. In contrast to the previously cited experimental
studies \cite{Bitog2011,Carvalho2006,Richmond1993,Sierra2008}, here, we do not
specify beforehand a particular type of biological system but describe, in a
general way, a biological substrate-biomass reaction in a continuous reactor.
Compared to the works developed in
\cite{Ansoni2016,Coenen2013}, we couple the fluid flow with the biological
phenomena, while in\cite{Ansoni2016,Coenen2013} the authors only model one of
the two physics.
Furthermore, in this case, the reactor geometry is parametrized with five
variables (compared to the two-dimensional parametrization performed in
\cite{Ansoni2016,Coenen2013}) to be able to obtain a broader range of
possible bioreactor shapes.
This article is organized as follows:
Section \ref{sec:Modeling} we introduce a mathematical model describing the
dynamics of the bioreactor.
In Section \ref{sec:Optimization}, we state the optimization problem which aims to minimize the
considered reactor volume, with an outflow substrate concentration maintained to a
desired threshold, by choosing a suitable bioreactor shape.
In Section \ref{sec:Results}, we explain the numerical experiments carried out
for the optimization problem and show the results.
Section \ref{sec:Conclusions} draws the conclusions after the comparison
between the obtained optimized reactors.
\section{Mathematical Modeling}\label{sec:Modeling}
Let a vertical cylinder denoted by $\Omega^{\ast}$ be the domain of the bioreactor
in consideration. When the problem is initialized, there is a certain amount of
biomass inside $\Omega^{\ast}$ that reacts with the polluted water entering the
device through the inlet $\Gamma_{\rm in}^{\ast}$ (the upper boundary of the cylinder). Treated water leaves the reactor through the outlet $\Gamma_{\rm out}^{\ast}$ (the lower boundary of the cylinder). We denote $\Gamma_{\rm wall}^{\ast}= \delta \Omega^{\ast} \setminus (\Gamma_{\rm in}^{\ast} \cup \Gamma_{\rm out}^{\ast})$, where null flux is considered.
We present the following model to describe the dynamics in the reactor,
which includes advection-diffusion-reaction phenomena
(see \cite{Crespo2015,Crespo2017,Crespo2016,diaz2017}):
\begin{equation}\label{eq:ADV}
\begin{gathered}
\frac{{\rm d} S}{{\rm d} t}= \nabla \cdot (D_{\rm S} \nabla S
- \mathbf{u} S) - \mu(S) B \quad \text{in }\Omega^{\ast} \times (0,T),\\
\frac{{\rm d} B}{{\rm d} t}= \nabla \cdot (D_{\rm B} \nabla B - \mathbf{u} B)
+ \mu(S) B \quad \text{in }\in \Omega^{\ast} \times (0,T), \\
S(x,0)= S_{\rm 0}(x),\text{ }B(x,0)= B_{\rm 0}(x) \quad x \in \Omega^{\ast},\\
\mathbf{n} \cdot (- D_{\rm S}\nabla S + \mathbf{u} S)
= S_{\rm in} u_{3} \quad \text{in } \Gamma_{\rm in}^{\ast} \times (0,T),\\
\mathbf{n} \cdot (- D_{\rm B}\nabla B + \mathbf{u} B)=0 \quad \text{in }
\Gamma_{\rm in}^{\ast} \times (0,T),\\
\mathbf{n} \cdot (-D_{\rm S}\nabla S)
=\mathbf{n} \cdot (-D_{\rm B}\nabla B)=0 \quad \text{in }
\Gamma_{\rm out}^{\ast} \times (0,T),\\
\mathbf{n} \cdot (-D_{\rm S} \nabla S + \mathbf{u}S)
=\mathbf{n} \cdot (-D_{\rm B} \nabla B + \mathbf{u}B)=0 \quad
\text{in } \Gamma_{\rm wall}^{\ast} \times (0,T),
\end{gathered}
\end{equation}
where $T>0 $ (s) is the length of the time interval for which we want to model
the process, $S$ (kg/m$^{3}$), $B$ (kg/m$^{3}$) are the substrate and biomass
concentrations inside the reactor which diffuse throughout the water in the
vessel with diffusion coefficients $D_{\rm S}$ (m$^2$/s) and
$D_{\rm B}$ (m$^2$/s), respectively, $S_{0}$ (kg/m$^3$), $B_{0}$
(kg/m$^{3}$) are the concentrations of substrate and biomass inside
the bioreactor at the initial time, $S_{\rm in}$ (kg/m$^3$) is the
substrate concentration that enters the reactor and $\mathbf{n}$ is the
outward unit normal vector on the boundary of the domain $\Omega^{\ast}$.
Notice that besides the advection-diffusion terms, we also have a term
corresponding to the reaction of biomass and substrate, governed by the
growth rate function $\mu$ (s$^{-1}$). We work with the following growth rate
functions, which are extensively used in the literature --
the Monod function \cite{SmithWaltman,BookAlain} is defined in $[0,+\infty)$ by
\begin{equation}\label{eq:MONOD}
\mu(S) = \mu_{\rm max} \frac{S}{K_{\rm S}+S},
\end{equation}
where $\mu_{\rm max}$ ($s^{-1}$) is the maximum specific growth rate and
$K_{\rm S}$ (kg/m$^3$) is the half-saturation constant. The Haldane function
\cite{Andrews,BookAlain} is defined in $[0,+\infty)$ by
\begin{equation}\label{eq:HALDANE}
\mu(S)= \mu^{\ast}\frac{S}{K_{\rm S}+S+S^2/K_{\rm I}},
\end{equation}
where $\mu^{\ast}$ (s$^{-1}$) is the maximum specific growth rate in the
absence of inhibition and $K_{\rm I}$ (kg/m$^3$) is the inhibition constant.
Finally, vector $\mathbf{u}=(u_1,u_2,u_{3})$ (m/s) is the flow velocity,
which fulfills the following stationary Navier-Stokes equations for Newtonian
incompressible viscous fluids (see, e.g., \cite{glowinski2013})
\begin{equation}\label{eq:NS}
\begin{gathered}
-\eta \Delta \mathbf{u} + \rho (\mathbf{u}\cdot \nabla) \mathbf{u}
+ \nabla p = 0\quad \text{in } \Omega^{\ast},\\
\nabla \cdot \mathbf{u} =0 \quad \text{in } \Omega^{\ast},\\
\mathbf{u}=0 \quad \text{in } \Gamma_{\rm wall}^{\ast},\\
\mathbf{u}= - u_{\rm in}\text{ }E(x)\text{ }\mathbf{n} \quad
\forall x \in \Gamma_{\rm in}^{\ast},\\
\mathbf{n} \cdot (\eta \nabla \mathbf{u}) =0 \quad \text{in }
\Gamma_{\rm out}^{\ast},\\
p(x)=p_{\rm atm} \quad \forall x \in \Gamma_{\rm out}^{\ast},
\end{gathered}
\end{equation}
where $p$ is the pressure field (Pa); $p_{\rm atm}$ is the atmospheric pressure
(Pa); $\eta$ is the fluid dynamic viscosity (kg/m s); $u_{\rm in}$ (m/s) is
the maximum injection velocity; $E$ is the laminar flow inlet profile
(a paraboloid of revolution) equal to $0$ in the inlet border and unity in the
inlet center; $\rho$ is the fluid density (kg/ m$^3$), assumed constant
through the reactor (as done, e.g., in \cite{Ansoni2016,Coenen2013}).
\begin{remark} \rm
Note that the flow field $\mathbf{u}$ has been considered stationary in
order to reduce the computational complexity met when numerically solving
system \eqref{eq:ADV}-\eqref{eq:NS}. This assumption is supported by numerical
experiments, which seem to show that if we solve a time-dependent version
of \eqref{eq:NS} coupled with \eqref{eq:ADV}, variable $\mathbf{u}$ approximates
its stationary state much faster than variables $(S,B)$.
\end{remark}
\begin{remark} \rm
According to \cite{Crespo2015}, if $S_{\rm in}\in L^{\infty}(0, T)$,
$S_{\rm in} \geq 0$ in $(0,T)$, $S_0 \in L^{\infty}(\Omega^{\ast})$,
$S_0 \geq 0$ in $\Omega^{\ast}$, $B_0\in L^{\infty}(\Omega^{\ast})$ and
$B_0\geq 0$ in $\Omega^{\ast}$, there exists a unique solution
$(S,B) \in L^2(0,T; H^{1}(\Omega^{\ast}))^2
\cap \mathcal{C}(0,T,L^2(\Omega^{\ast}))^2
\cap L^{\infty}(\Omega^{\ast} \times (0,T))^2$ of system
\eqref{eq:ADV}-\eqref{eq:NS}.
\end{remark}
As we will see in Section \ref{sec:Results}, we aim to find stationary
solutions of system \eqref{eq:ADV}, which we denote by $(\hat{S},\hat{B})$.
A usual way to get them is to solve numerically \eqref{eq:ADV} and then
consider its solution for a time value $\hat{T}$ (large enough) as the
steady-state solution. This is usually computationally easier
(see, e.g., \cite{leveque2007}) and allows us to recover non-trivial
stationary solutions (different to $(\hat{S},\hat{B})=(S_{\rm in},0)$)
by choosing appropriate initial conditions.
\section{Optimization problem}\label{sec:Optimization}
Here, we optimize the main design parameters (reactor shape and total volume)
with respect to the output concentration. More precisely, in
Section \eqref{sec:GeneralProblem} we introduce the general formulation of
the considered continuous optimization problem.
Then, in Section \eqref{sec:NumericalProblem} we propose three particular
discrete implementations of this problem to be solved numerically in
Section \ref{sec:Results}.
\subsection{General problem}\label{sec:GeneralProblem}
Let us consider cylindrical bioreactors $\Omega^{\ast}$ which are empty solids
of revolution, and so, they can be described by using a 2D domain
$\Omega \subset \mathbb{R}^2$ (similar to the one depicted by
Figure \ref{fig:ShapeGeneral}) by using cylindrical coordinates $(r,z)$,
where $r$ is the distance to the cylinder axis. The simplified domain $\Omega$
is described as follows: $H$ (m) is the bioreactor height; $r$ (m) is the
radius of the inlet $\Gamma_{\rm in}$ and the outlet $\Gamma_{\rm out}$; $h$ (m)
is the height of the inlet and outlet pipes; $R_1$ (m) and $R_2$ (m) are
the radius of the bioreactor wall perpendicular to the inlet and outlet pipes,
respectively; the curve of the exterior wall corresponds to the graph of the
function $\psi: [h,H-h]\to [r,+\infty)$, which satisfies $\psi(h)=R_2$
and $\psi(H-h)=R_1$.
Since, in practice, the inlet and outlet pipes have standard dimensions
(depending on the desired industrial application), we assume that $r$ and $h$
have fixed values. Similarly, we take into account that the height and width
of the reactor cannot exceed certain values (for example, due to a limitation
of the physical space in an industrial factory).
\begin{figure}[ht]
\centering
\includegraphics[width=0.5\textwidth]{fig1}
\caption{Schematic representation of the bioreactor geometries used to
solve problem \eqref{ProblemGeneral}.
The exterior curve (depicted in blue), which corresponds to part of the
bioreactor exterior wall, is defined as $(z,\psi(z))$, where $z \in [h,H-h]$.}
\label{fig:ShapeGeneral}
\end{figure}
Given a prescribed output substrate concentration $S_{\rm lim}$ (kg/m$^3$),
we state the following optimization problem
\begin{equation} \label{ProblemGeneral}
\begin{gathered}
\text{Find }\phi^{\rm opt}\in \Phi, \text{ such that} \\
\operatorname{Vol}(\phi^{\rm opt})= \min_{\phi \in \Phi} \operatorname{Vol}(\phi),\\
\hat{S}_{\rm out}(\phi^{\rm opt})\leq S_{\rm lim},
\end{gathered}
\end{equation}
where $\phi=(H,R_1,R_2,\psi) \in \Phi$ defines a particular bioreactor
shape and $\Phi=\{ [H_{\rm min},H_{\rm max}]\times [r,R_{1,{\rm max}}]
\times[r,R_{2,{\rm max}}]\times \mathcal{C}([h,H-h],[r, R_{\rm max}])\}$
is the admissible space; $\operatorname{Vol}(\phi)$ (m$^3$) is the volume of the reactor,
computed as
\begin{equation}\label{eq:Vol}
\operatorname{Vol}(\phi)= \int_{\Omega^{\ast}(\phi)} {\rm d}x {\rm d}y {\rm d}z,
\end{equation}
with $\Omega\subset \mathbb{R}^2$ being the $(r,z)$-domain obtained with the
set $\phi$ and $\Omega^{\ast}(\phi)\subset \mathbb{R}^{3}$ is the corresponding
3D domain; and $\hat{S}_{\rm out}(\phi)$ (kg/m$^3$) denotes an average of the
concentration of substrate that leaves the bioreactor (at steady state),
computed as
\begin{equation}\label{eq:Sout}
\hat{S}_{\rm out}(\phi)
= \frac{\int_{\Gamma_{\rm out}^{\ast}} S_{\phi}(x,y,0,\hat{T})
|u_{3}(x,y,0)|{\rm d}x{\rm d} y}{\int_{\Gamma_{\rm out}^{\ast}}
|u_{3}(x,y,0)|{\rm d}x{\rm d} y},
\end{equation}
with $S_{\phi}(\cdot,\cdot,\cdot,\hat{T})$ the solution of system \eqref{eq:ADV},
obtained with the 3D domain $\Omega^{\ast}(\phi)$, at time $\hat{T}$ and $u_3$
the third component of the velocity vector in \eqref{eq:NS}.
\begin{remark} \rm
In practice, we aim to improve (in terms of reactor volume) a given bioreactor
which attains the prescribed output substrate concentration and we choose
the admissible space $\Phi$ so that it accounts for the geometry of the given vessel.
Therefore, if problem \eqref{ProblemGeneral} does not have a solution is because
the minimum does not exist (only an infimum value is ensured).
For the numerical experiments considered in Section \ref{sec:Results},
we look for one of the optimal solutions of the discretized
problems \eqref{ProblemP1}, \eqref{ProblemP2} and \eqref{ProblemP3},
whose existence is guaranteed because of the compactness of the selected
admissible spaces $\tilde{\Phi}^{i}$ ($i=1,2,3$).
\end{remark}
\subsection{Numerical problem}\label{sec:NumericalProblem}
Here, we present three discrete versions of the optimization problem
\eqref{ProblemGeneral}, related to three different discrete parametrizations
of the bioreactor geometry. The parametrization proposed in Section \ref{sec:P1}
allows us to model tubular shapes (typically used in the industry sector),
while the parametrizations in Sections \ref{sec:P2} and \ref{sec:P1} offer the
possibility to obtain a wider range of reactor geometries. We solve the discrete
optimization problems \eqref{ProblemP1}, \eqref{ProblemP2} and \eqref{ProblemP3}
by using the Hybrid Genetic Algorithm, and its parameters, presented in
Section \ref{sec:ga}.
For the sake of simplicity, the objective function and the restriction in
problem \eqref{ProblemGeneral} are combined into a new objective function
$J(\phi)$ (m$^3$) as
\begin{equation}\label{eq:J}
J(\phi)= \operatorname{Vol}(\phi)\Big( 1 + \beta
\frac{\max(\hat{S}_{\rm out}(\phi)- S_{\rm lim},0)}{S_{\rm lim}} \Big),
\end{equation}
where $\beta$ is a free parameter (usually large) and the term multiplied by
the coefficient $\beta$ is a barrier function used to penalize solutions
with $S_{\rm lim}$ smaller than an average of the substrate concentration
exiting the bioreactor. Introducing the penalty term in \eqref{eq:J} allows
us to solve problem \eqref{ProblemGeneral} with the unconstrained (low-cost)
optimization algorithm proposed in Section \ref{sec:ga}.
This approach has been previously used in the literature for solving
computationally expensive industrial optimization problems
(see, e.g. \cite{ivorra2009}). As it is typically done when performing
optimization with a penalty method, the value of the penalty parameter
$\beta$ is usually obtained by performing several tests (see
\cite[Chapter 15.1]{nocedal2006}).
\begin{figure}[ht]
\centering
\subfigure[Domain in problem \eqref{ProblemP1}]
{\includegraphics[width=0.42\textwidth]{fig2a}}
\subfigure[Domain in problem \eqref{ProblemP2}]
{ \includegraphics[width=0.56\textwidth]{fig2b}}
\\
\subfigure[Domain in problem \eqref{ProblemP3}]
{\includegraphics[width=0.58\textwidth]{fig2c}}
\caption{Schematic representation of the bioreactor geometries used
to solve the discrete problems \eqref{ProblemP1}, \eqref{ProblemP2}
and \eqref{ProblemP3}.}
\label{fig:reactorshape}
\end{figure}
\subsubsection{Parametrization 1}\label{sec:P1}
As a first approach, we consider bioreactor geometries as depicted in
Figure \ref{fig:reactorshape}-(a). The exterior wall corresponds to the segment
$[h, H-h]\times \{R\}$, where $R \in [r,R_{\rm max}]$ (m).
The variable $\phi$ in problem \eqref{ProblemGeneral} is taken as
$\phi=(H,R,R,\psi)$, where
$ \psi: [h,H-h] \to [r,R_{\rm max}]$ with
\[
\psi(z)= R.
\]
In this case, the bioreactor geometry only depends on parameters $H$ and $R$
and the optimization problem \eqref{ProblemGeneral} can be reformulated as
\begin{equation} \label{ProblemP1}
\begin{gathered}
\text{Find }\tilde{\phi}^{1,{\rm opt}}\in \tilde{\Phi}^{1}, \text{ such that} \\
J(\tilde{\phi}^{1,{\rm opt}})= \min_{\tilde{\phi}^{1} \in \tilde{\Phi}^{1}}
J(\tilde{\phi}^{1}),
\end{gathered}
\end{equation}
where $\tilde{\phi}^{1,{\rm opt}}= (H^{\rm opt},R^{\rm opt})$ and
$\tilde{\Phi}^1:= \{(H,R) \in [H_{\rm min}, H_{\rm max}]\times[r,R_{\rm max}]\}
\subset \mathbb{R}^2$ is the admissible space.
\subsubsection{Parametrization 2}\label{sec:P2}
As a second approach, we consider bioreactor geometries as depicted in
Figure \ref{fig:reactorshape}-(b). The exterior wall corresponds to a
semi-ellipse with center $(r,\frac{H}{2})$ and with lengths of the semi-axis
given by the pair $(R-r,\frac{H-2h}{2})$, where $R \in [r,R_{\rm max}]$ (m).
The variable $\phi$ in problem \eqref{ProblemGeneral} is taken as
$\phi=(H,r,r,\psi)$, where
$\psi: [h,H-h] \to [r,R_{\rm max}] $ with
\[
\psi(z)= r + (R-r)\sqrt{ 1 - \big(\frac{z-H/2}{h-H/2}\big)^2}.
\]
It is straightforward to see that, if $R \in [r,R_{\rm max}]$, then
$\psi \in \mathcal{C}([h,H-h],[r,R_{\rm max}])$.
As in problem \eqref{ProblemP1}, the bioreactor geometry only depends on
parameters $H$ and $R$ and the optimization problem \eqref{ProblemGeneral}
can be reformulated as
\begin{equation} \label{ProblemP2}
\begin{gathered}
\text{Find }\tilde{\phi}^{2,{\rm opt}}\in \tilde{\Phi}^2, \text{ such that} \\
J(\tilde{\phi}^{2,{\rm opt}})= \min_{\tilde{\phi}^2
\in \tilde{\Phi}^2} J(\tilde{\phi}^2),
\end{gathered}
\end{equation}
where $\tilde{\phi}^{2,{\rm opt}}= (H^{\rm opt},R^{\rm opt})$ and
$\tilde{\Phi}^2:= \{(H,R) \in [H_{\rm min}, H_{\rm max}]\times[r,R_{\rm max}]\}
\subset \mathbb{R}^2$ is the admissible space.
\subsubsection{Parametrization 3}\label{sec:P3}
As a third approach, we consider bioreactor geometries as depicted in
Figure \ref{fig:reactorshape}-(c). The shape of the exterior wall is a
quadratic \emph{B\'ezier curve} (see, for example, \cite{farin2014}),
associated with the control points $\mathbf{P}=(R_{\rm 1},H-h)$,
$\mathbf{Q}=(R_{\rm 2},h)$ and $\mathbf{E}=(E_{\rm 1},E_{\rm 2})$, where
$(E_1,E_2)\in [E_{1,{\rm min}},E_{1,{\rm max}}]\times [E_{2,{\rm min}},
E_{2,{\rm max}}]$), by the formula
\begin{equation}\label{beziercurve}
\mathbf{B}(\sigma)= (B_1(\sigma),B_2(\sigma))
= (1-\sigma)^2 \mathbf{P} + 2(1-\sigma)\sigma \mathbf{E}
+ \sigma^2\mathbf{Q}, \quad \sigma \in [0,1].
\end{equation}
The variable $\phi$ appearing in problem \eqref{ProblemGeneral} is taken as
$\phi=(H,R_1,R_2,\psi)$, where
$ \psi: [h,H-h] \to [r,R_{\rm max}] $ wiht
\[
\psi(z)= B_1(B_2^{-1}(z)).
\]
To assure that $\psi ([h,H-h])\subseteq [r,R_{\rm max}]$, we impose the
radius expansions $R_1$ and $R_2$ to lie in the segment $[r,R_{\rm max}]$.
Once $R_1$ and $R_2$ have been chosen we take $E_{1, {\rm min}}$ and
$E_{1,{\rm max}}$ as in Lemma \ref{lem:cuentasB1B2} below so that
$B_1 \circ B_2^{-1} \in \mathcal{C}([h,H-h],[r,R_{\rm max}])$.
Furthermore, we take $E_{2, {\rm min}}=h$ and $E_{2,{\rm max}}=H-h$.
Now, we define two new optimization parameters $\alpha_1,\alpha_2\in [0,1]$
such that
\begin{equation}\label{eq:alphaopt}
E_1= E_{1,{\rm min}} + \alpha_1 \cdot (E_{1,{\rm max}}-E_{1,{\rm min}})
\quad \text{and}\quad E_2= h + \alpha_2\cdot (H-2h).
\end{equation}
In that case, the bioreactor geometry only depends on parameters $H$, $R_1$, $R_2$,
$\alpha_1$ and $\alpha_2$. The solution of the optimization problem
\eqref{ProblemGeneral} is approximated by computing
\begin{equation} \label{ProblemP3}
\begin{gathered}
\text{Find }\tilde{\phi}^{3,{\rm opt}}\in \tilde{\Phi}^{3}, \text{ such that} \\
J(\tilde{\phi}^{3,{\rm opt}})
= \min_{\tilde{\phi}^{3} \in \tilde{\Phi}^{3}} J(\tilde{\phi}^{3}),
\end{gathered}
\end{equation}
where $\tilde{\phi}^{3,{\rm opt}}= (H^{\rm opt},R_1^{\rm opt},
R_2^{\rm opt},\alpha_1^{\rm opt},\alpha_2^{\rm opt})$ and
$$
\tilde{\Phi}^{3}:= \{(H,R_1,R_2,\alpha_1,\alpha_2)
\in [H_{\rm min},H_{\rm max}]\times [r,R_{1,{\rm max}}]
\times [r,R_{2,{\rm max}}]\times[0,1]^2\} \subset \mathbb{R}^{5}
$$
is the admissible space, with $R_{i,{\rm max}}\leq R_{\rm max}$ for $i=1,2$.
\begin{lemma}\label{lem:cuentasB1B2}
Let us denote $E_{1,{\rm min}}= r-\sqrt{(R_1-r)(R_2-r)}$ and
$E_{1,{\rm max}}= R_{\rm max}+\sqrt{(R_1-R_{\rm max})(R_2-R_{\rm max})}$.
If $(E_1,E_2)\in [E_{1,{\rm min}},E_{1,{\rm max}}]\times [h,H-h]$, then
$B_1 \circ B_2^{-1} \in \mathcal{C}([h,H-h],[r,R_{\rm max}])$.
\end{lemma}
\begin{proof}
We divide the proof in four steps:
\smallskip
\noindent\textbf{Step 1.} Let us prove that, if $E_2\in [h,H-h]$, then
$B_2([0,1])=[h,H-h]$.
To obtain the minimum and maximum values of $B_2(\sigma)$, $\sigma \in [0,1]$,
we compute the critical points $\sigma_2^{\ast}$ satisfying the equation
$ \frac{{\rm d} B_2}{{\rm d} \sigma}(\sigma_2^{\ast})=0$ on the interior of
$]0,1[$. When considering $E_2$ as a variable, one can see that
$\sigma_2^{\ast}$ depends on $E_2$ through the expression
$\sigma_2^{\ast}(E_2)= \frac{E_2-H+h}{2E_2-H}$, with corresponding value
$B_2(\sigma_2^{\ast}(E_2))= \frac{E_2^2+h^2-Hh}{2E_2-h}$.
Now, to find the lower and upper bounds for variable $E_2$
(assuring that $B_2(\sigma)\in [h,H-h]$ $\forall \sigma \in [0,1]$),
we respectively solve equations $B_2(\sigma_2^{\ast}(E_{2,{\rm m}}))= h$ and
$B_2(\sigma_2^{\ast}(E_{2,{\rm M}}))= H-h$.
It is easy to prove that the unique solutions of these equations are
$E_{2,{\rm m}}=h$ and $E_{2,{\rm M}}=H-h$. Finally, taking into account that
$ \frac{{\rm d} B^2_2}{{\rm d} \sigma^2}= 2H-4E_2$, it follows that
$ \frac{{\rm d}^2 B_2}{{\rm d} \sigma^2}\big|_{E_{\rm 2}=h}=2(H-2h)>0$ and
$ \frac{{\rm d}^2 B_2}{{\rm d} \sigma^2}\big|_{E_{\rm 2}=H-h}=2(2h-H)<0$,
and so, one can conclude that $E_{2,{\rm min}}=h$ and $E_{2,{\rm max}}=H-h$.
\smallskip
\noindent\textbf{Step 2.} Let us prove that the function $B_2:[0,1] \to [h,H-h]$
is injective.
Let $\sigma, \bar{\sigma} \in [0,1]$ satisfying $B_2(\sigma)=B_2(\bar{\sigma})$.
By definition, this implies that
$$
(1-\sigma)^2(H-h) + 2(1-\sigma)\sigma E_2 + \sigma^2 h
= (1-\bar{\sigma})^2(H-h) + 2(1-\bar{\sigma})\bar{\sigma} E_2 + \bar{\sigma}^2 h.
$$
Easy calculations lead to
$$
(H-h) \big(\sigma^2-\bar{\sigma}^2 - 2\sigma
+ 2 \bar{\sigma}\big) + 2E_2\big(\sigma-\bar{\sigma}-\sigma^2+\bar{\sigma}^2\big)
+ h\big(\sigma^2-\bar{\sigma}^2\big)=0.
$$
Denoting $x=\bar{\sigma}-\sigma$ and $y=\bar{\sigma}+\sigma$, the previous
equation can be rewritten as
\begin{gather*}
x(2-y)(H-h) + 2x(y-1) E_2 - xyh= 0 \\
\Leftrightarrow x \left( 2(H-h) -2E_2 + y(2E_2-H) \right)=0.
\end{gather*}
This implies that either $x=0$ or $y= \frac{2(H-h-E_2)}{H-2E_2}$.
In the second case, it is easy to see that $y= 1+ \frac{H-2h}{H-2E_2}$ and,
since we assume that $E_2>h$, it follows that $y>2$, but this enters in
a contradiction with the definition of $y$. Thus, we can conclude that $x=0$,
so $\sigma=\bar{\sigma}$ and the injectivity is proved.
\smallskip
\noindent\textbf{Step 3.} Let us prove that, if
$E_1\in [E_{1,{\rm min}},E_{1,{\rm max}}]$, then $B_1([0,1])=[r,R_{\rm max}]$.
Similarly to step 1, to obtain the minimum and maximum values of $B_1(\sigma)$,
$\sigma \in [0,1]$, we compute the critical points $\sigma_1^{\ast}$ satisfying
the equation $ \frac{{\rm d} B_1}{{\rm d} \sigma}(\sigma_1^{\ast})=0$ on the
interior of the domain. When considering $E_1$ as a variable, one can see that
$\sigma_1^{\ast}$ depends on $E_1$ through the expression
$\sigma_1^{\ast}(E_1)= \frac{R_1-E_1}{R_1+R_2-2E_1}$, with corresponding value
$B_1(\sigma_1^{\ast}(E_1))= \frac{R_1R_2 - E_1^2}{R_1+R_2-2E_1}$.
Now, to obtain lower and upper bounds for the variable $E_1$ (assuring that
$B_1(\sigma)\in [r,R_{\rm max}]$ $\forall \sigma \in [0,1]$), we respectively
solve equations $B_1(\sigma_1^{\ast}(E_{1,{\rm m}}))= r$ and
$B_1(\sigma_1^{\ast}(E_{1,{\rm M}}))= R_{\rm max}$.
Each of these equations has two solutions, given by
$E_{1,{\rm m}}^{\pm}= r \pm \sqrt{(R_1-r)(R_2-r)}$ and
$E_{1,{\rm M}}^{\pm}= R_{\rm max} \pm \sqrt{(R_1-R_{\rm max})(R_2-R_{\rm max})}$.
Taking into account that $ \frac{{\rm d}^2 B_1}{{\rm d} \sigma^2}= 2( R_1+R_2-2E_1)$,
$$
\frac{{\rm d}^2 B_1}{{\rm d} \sigma^2}\big|_{E_{\rm 1}
=E_{1,{\rm m}}^{\pm}} = 2(R_1+R_2-2r) \mp 4\sqrt{(R_1-r)(R_2-r)}
$$
and
$$
\frac{{\rm d}^2 B_1}{{\rm d} \sigma^2}\big|_{E_{\rm 1}=E_{1,{\rm M}}^{\pm}}
= 2(R_1+R_2-2R_{\rm max}) \mp 4\sqrt{(R_1-R_{\rm max})(R_2-R_{\rm max})},
$$
and so, one can conclude that $E_{1,{\rm min}}= E_{1,{\rm m}}^{-}$ and
$E_{1,{\rm max}}= E_{1,{\rm M}}^{+}$.
\smallskip
\noindent\textbf{Step 4.}
Let us prove that $B_1 \circ B_2^{-1} \in \mathcal{C}([h,H-h],[r,R_{\rm max}])$.
Since $B_1:[0,1]\to [r,R_{\rm max}]$ and $B_2:[0,1]\to[h,H-h]$ are continuous
functions and $B_2$ is injective, we conclude that $B_1 \circ B_2^{-1}$ is
well defined and continuous because it is the composition of continuous functions.
\end{proof}
\subsubsection{Optimization algorithm}\label{sec:ga}
In this section, we describe in detail the optimization algorithm and the parameters
used to solve numerical problems \eqref{ProblemP1}, \eqref{ProblemP2}
and \eqref{ProblemP3}. For the sake of simplicity, here, we consider a general
optimization problem
\begin{equation}
\min_{x \in \Theta} J(x) \label{problema}
\end{equation}
where $J:x \to \mathbb{R}$ is the fitness function; $x$ is the
optimization parameter, $\Theta \subset \mathbb{R}^{N}$, with $N\in \mathbb{N}$,
is the search space. Notice that in order to recover problems
\eqref{ProblemP1}, \eqref{ProblemP2} and \eqref{ProblemP3} one should replace
$x= \tilde{\phi}^{i}$ ($i=1,2,3$), $\Theta= \tilde{\Phi}^{i}$ ($i=1,2,3$) and
$N=2,2,5$ in the general problem formulation \eqref{problema}.
The proposed algorithm, called Genetic Multi-Layer Algorithm \text{(GMA)},
is a global optimization method based on a hybridization between a genetic
algorithm \text{(GA)} \cite{deb2001,goldberg1989,sharapov2007,weise2009}
(which performs a global search of the solution) and a multi-layer secant
method \text{(MSA)} \cite{redondo1,redondo2,redondo3} (which provides suitable
initial populations for the GA). A complete validation of these algorithms on
various industrial problems can be found in
\cite{gomez2011,Isebe2008,ivorra2006,ivorra20062,ivorra2007,ivorra2009,ivorra2013,
ivorra20132,ivorra2018} and references therein. Broadly speaking, GAs are search
techniques which try to solve problems similar to \eqref{problema} through a
stochastic process based on a natural selection process that mimics biological
evolution. The \text{GA}s have many advantages as, for example, they can solve
complex optimization problems (e.g., with high dimensional search space or
function with various with local minima). However, they exhibit lower accuracy
than other methods, such as gradient algorithms. Before explaining the
methodology used to enhance these inconveniences, we detail the \text{GA}
used in this work.
\subsection*{Genetic algorithm scheme:} \quad
\smallskip
\noindent\textbf{Step 1.} Inputs: The user must define six parameters:
$N_{\rm p} \in \mathbb{N}$, $N_{\rm g} \in \mathbb{N}$, $p_{\rm c} \in [0, 1]$,
$p_{\rm m}\in [0, 1]$, $\lambda \in \mathbb{R}$ and $\hat{g}\in \mathbb{N}$,
the meaning of which is clarified later in the following steps. In addition,
the user needs to provide a first set, called \emph{initial population} and
denoted by $X^0 = \{ x^0 _j \in \Phi , j=1,\dots,N_{\rm p}\}\in \Theta$,
\textit{of} $N_{\rm p}$ possible solutions of the optimization problem
\eqref{problema}. Each row $x_{j}^{0}$ in $X^{0}$ ($j=1, \ldots, N_{\rm p}$)
is called \emph{individual}, while each component $x_{j,k}^{0}$ of an individual
($k=1,\ldots,N$) is called \emph{gene}.
In our case $\Theta= \prod_{k=1}^{N} [l_{k},u_{k}]$, where $l_{k}$ and $u_{k}$
are respectively the lower and upper bounds \textit{of} the
gene $x_{j,k}^{i}$.
\smallskip
\noindent\textbf{Step 2.} Generation of new populations: Starting from the
initial population $X^{0}$, we recursively create $N_{\rm g} \in \mathbb{N}$
new populations, which we call \emph{generations}, by applying $4$ stochastic
steps, called \emph{selection}, \emph{crossover}, \emph{mutation} and \emph{elitism},
which are described in Steps 3.1, 3.2, 3.3 and 3.4, respectively.
More precisely, let $X^{i} = \{ x_{j}^{i}\in \Theta, j = 1,dots, N_{\rm p}\}$
with $i = 1, \dots , N_{\rm g}-1$, denotes the population at iteration $i$.
Thus, using the following $(N_{\rm p}, N)$-real valued matrix notation
\[
X^{i} = \begin{bmatrix}
x_1^{i}(1) & \cdots & x_1^{i}(N) \\
\vdots & \vdots & \vdots \\
x_{N_{\rm p}}^{i}(1) & \cdots & x_{N_{\rm p}}^{i}(N)
\end{bmatrix},
\]
$X^{i+1}$ is obtained by considering
\begin{equation}\label{ecgenetico}
X^{i+1}= (I_{N}-\mathcal{E}^{i}) (\mathcal{C}^{i}\mathcal{S}^{i}X^{i}
+ \mathcal{M}^{i}) + \mathcal{E}^{i}X^{i}
\end{equation}
where matrices $\mathcal{S}^{i}$, $\mathcal{C}^{i}$, $\mathcal{M}^{i}$,
$\mathcal{E}^{i}$ and $I_{N}$ are described below.
\smallskip
\noindent\textbf{Step 2.1} Selection: This operator is used to select
individuals according to their fitness value. There exist various selection
techniques (see, for instance, \cite{goldberg1989,sharapov2007,weise2009}),
among which we use the \emph{Roulette Wheel Selection} method. We randomly select
$N_{\rm p}$ individuals from $X^{i}$ with eventual repetitions.
Each individual $x_{j}^{i}\in X^{i}$, with $j=1, \ldots, N_{\rm p}$ has a
probability to be selected during this process which is given by
$J(x_{j}^{i})^{-1} / \sum_{k=1}^{N_{\rm p}} J(x_{k}^{i})^{-1}$.
This step can be summarized as
$$
X^{i+1,1}= \mathcal{S}^{i}X^{i},
$$
where $\mathcal{S}^{i}$ is a binary valued $(N_{\rm p},N_{\rm p})$-matrix
satisfying $S_{j,k}^{i}=1$ if the k-th individual in $X^{i}$ is the j-th
selected individual, and $S_{j,k}^{i}=0$ in other case.
\smallskip
\noindent\textbf{Step 2.2}
Crossover: This operator is used to create a new individual by combining the genes
of two existing individuals from the population $X^{i}$ (chosen during the
previous selection process). There are several methods for combining individuals
(see, for instance, \cite{deb2001,goldberg1989,sharapov2007}), among which we
use the \emph{Arithmetic Crossover} method.
For each pair of consecutive individuals (rows) $2j - 1$ and $2j$ in $X^{i+1,1}$,
with $1 \leq j \leq \operatorname{floor}(N_{\rm p}/2)$
(where $\operatorname{floor}(a)$
is the nearest integer lower than or equal to $a$), we determine, with a
probability $p_{\rm c}$, if those rows exchange data or if they are directly
copied into an intermediate population denoted by $X^{i+1,2}$. Thus, matrix
$\mathcal{C}^{i}$ is a real valued matrix of size $(N_{\rm p},N_{\rm p})$,
satisfying
$$
\mathcal{C}_{2j-1,2j-1}^{i}=\lambda_1\text{,}\quad
\mathcal{C}_{2j-1,2j}^{i}=1 - \lambda_1\text{,}\quad
\mathcal{C}_{2j,2j}^{i}=\lambda_2,\quad
\mathcal{C}_{2j,2j-1}^{i}=1 - \lambda_2,
$$
where $\lambda_1=\lambda_2=1$ with probability $1-p_{\rm c}$, or
$\lambda_1,\lambda_2$ are randomly chosen in $(0,1)$ considering a uniform
distribution, in other case. Other coefficients of $\mathcal{C}^{i}$ are set to
$0$. If $N_{\rm p}$ is odd, then we also set
$\mathcal{C}^{i}(N_{\rm p},N_{\rm p})=1$, and then the $N_{\rm p}$-th row of
$X^{i+1,1}$ is directly copied into $X^{i+1,2}$.
\smallskip
\noindent\textbf{Step 2.3} Mutation: This operator randomly modifies the value
of one or more genes of an individual from the population $X^{i+1,2}$
(obtained during the previous crossover process). It provides diversity in the
population and intends to avoid the premature convergence phenomenon
(i.e., population concentrated near a local minimum, see \cite{goldberg1989}).
Each individual can be mutated with a probability $p_{\rm m}$ given by the user.
There exist different techniques to randomly mutate individuals
(see, for instance, \cite{deb2001,sharapov2007}), among which we use the
\emph{Non-Uniform Mutation} method.
We decide, with a probability $p_{\rm m}$, if each row of $X^{i+1,2}$ is
randomly perturbed or not. This step is defined by
$$
X^{i+1,3}=X^{i+1,2} + \mathcal{M}^{i},
$$
where $\mathcal{M}^{i}$ is a real valued matrix with size $(N_{\rm p},N)$
and the $j$-th row satisfies
$$
\mathcal{M}^{i}_{j} =\begin{cases}
\vec{0} & \text{with probability } 1-p_{\rm m}\\
\Delta(g,x_{j}^{i}) & \text{in other case}
\end{cases}
$$
and the $k$-th component of the vector $\Delta(g,x_{j}^{i})$ is defined as
$$
\Delta(g,x_{j}^{i})
= \begin{cases}
(u_{k}-x^{i}_{j,k})(1 - \gamma^{(1 - \frac{g}{N_{\rm g}})^{\lambda}})
& \text{if } \tau=0 \\
(l_{k}-x_{j,k}^{i})(1 - \gamma^{(1 - \frac{g}{N_{\rm g}})^{\lambda}})
& \text{if } \tau=1 \end{cases}
$$
where $g$ is the current generation number, $\tau$ is a binary random number,
$\gamma$ is a uniform random number in $[0,1]$ and $\lambda$ is a parameter
given by the user, determining the degree of dependency on the iteration number.
This mutation method decreases the mutation rate as the generation number increases.
\smallskip
\noindent\textbf{Step 2.4} Elitism: This operator ensures that at least one
of the best individuals of the current generation is directly copied to the
next generation. The main advantage of elitism is that a decreasing convergence
is guaranteed. For more details about elitism methods see, for instance,
\cite{sharapov2007,weise2009}.
Let $x^{i}_{b}$, where $b \in {1,\dots,N_{\rm p}}$, be the individual in
$X^{i}$ with the lowest value of the fitness function (or, if there exist
various, one of those individuals selected randomly with a uniform distribution).
If $x^{i}_{b}$ has a lower fitness value than all the individuals in $X^{i+1,3}$,
it is directly copied at the $b$-th row of $X^{i+1}$. This step can be formalized as
$$
X^{i+1} = (I_{N}- \mathcal{E}^{i}) (X^{i+1,3})+ \mathcal{E}^{i} X^{i},
$$
where $I_N$ is the identity matrix of size $N$ and $\mathcal{E}^{i}$ is a
real-valued $(N_{\rm p},N_{\rm p})$-matrix such that $\mathcal{E}^{i}(b,b)=1$
if $x^{i}_{b}$ has a lower fitness value than all the individuals in
$X^{i+1,3}$ and 0 otherwise, $\mathcal{E}^{i}=0$ elsewhere.
The genetic search is terminated when $N_{\rm g}$ generations have been computed,
or after a number of generations specified by the user, $\hat{g}$, without
improvement of the fitness value (i.e., the fitness of the best element has not
decreased).
\smallskip
\noindent\textbf{Step 3.} Output: When \text{GA} stops, it returns, as
an output solution, the individual who has the lowest value for the objective
function $J$ among all the individuals in all the populations considered during
the whole evolving process, i.e.,
\begin{align*}
&GAO(X^{0};N_{\rm p};N_{\rm g}: p_{\rm m}; p_{\rm c};\lambda;\hat{g})\\
&= \operatorname{argmin} \big\{J(x_{i}^{j})|\text{ }x_{i}^{j}
\text{ is the -th row of } X^{i},\; i=1,\ldots,N_{\rm g},\;
j=1,\ldots,N_{\rm p} \big\}.
\end{align*}
As said at the beginning of this section, to accelerate the convergence and
improve the accuracy of the above-described \text{GA}, we combine it with
the \text{MSA} described below to build a hybrid algorithm, called \text{GMA}.
Its general scheme is as follows:
\medskip
\subsection*{Genetic multi-Layer algorithm scheme:}\quad
\smallskip
\noindent\textbf{Step 1.} Inputs: The user must define seven parameters:
$s_{\rm max} \in \mathbb{N}$, $N_{\rm p} \in \mathbb{N}$,
$N_{\rm g} \in \mathbb{N}$, $p_{\rm c} \in [0, 1]$, $p_{\rm m}\in [0, 1]$,
$\lambda \in \mathbb{R}$ and $\hat{g}\in \mathbb{N}$. $s_{\rm max}$
denotes the number of iterations of the \text{MSA}.
\smallskip
\noindent\textbf{Step 2.} Initial population: A first family of possible
solutions of the optimization problem \eqref{problema}, denoted
\textit{by} $X_1^0= \{x^0_{1,j} \in \Theta, j=1,\dots,N_{\rm p}\}$,
is randomly generated in the search space $\Theta$ considering a uniform
distribution.
\smallskip
\noindent\textbf{Step 3.} Main loop: For $s$ from 1 to $s_{\max}$:
\smallskip
\noindent\textbf{Step 3.1} We run the \text{GA} starting from the initial
population $X_{s}^{0}$ and obtain the optimal individual
$o_s=GAO (X_{s}^0 ,N_{\rm p},N_{\rm g},p_{\rm m},p_{\rm c},\lambda,\hat{g})$.
\smallskip
\noindent\textbf{Step 3.2} We build a new initial population for the \text{GA},
$X_{s+1}^0=\{ x_{s+1,j}^0 \in \Theta, j=1,\dots,N_{\rm p} \}$, by considering
a secant method between each element in $X_{s}^{0}$ and the optimal individual
$o_{s}$, i.e., for all $j\in \{ 1,\dots,N_{\rm p} \}$, if $J(o_s)=J(x_{s,j}^0)$
we set
$$
x_{s+1,j}^0=x^0_{s,j},
$$
else we set
$$
x_{s+1,j}^0= \operatorname{proj}_{\Theta}(x^0_{s,j}-J(o_s)
\frac{o_s-x_{s,j}^0}{J(o_s)-J(x_{s,j}^0)}),
$$
where $\operatorname{proj}_{\Theta}: \mathbb{R}^{N} \to \Theta$ is the
projection function for controlling that the new individuals fit into the
search space $\Theta$, defined as
$\operatorname{proj}_{\Theta}(x)(k)=\min(\max(x(k),l_{k}),$ $u_{k})$, with
$k=1,\dots,N$.
\smallskip
\noindent\textbf{Step 4.} Output: After $s_{\rm max}$ iterations, the
\text{GMA} returns the following output:
$$
GMAO(s_{\rm max},N_{\rm p},N_{\rm g},p_{\rm m},p_{\rm c},\lambda,\hat{g})
=\operatorname{argmin}\{ J(o_s)|\text{ }s= 1,\dots,s_{\rm max} \}.
$$
This algorithm tries to improve, individual by individual, the initial
population of the \text{GA}. More precisely, for each individual in
the initial population:
\begin{itemize}
\item If there is a significant evolution of the cost function between this
individual and $o_s$, the secant method generates a new individual close to
$o_s$ that performs a refined search near the current solution.
\item Otherwise, the secant method creates a new individual far from $o_s$,
to expand the exploration of the admissible space.
\end{itemize}
Moreover, when the \text{GMA} ends, its solution is improved by performing
$10$ iterations of the Steepest Descent \text{(SD)} algorithm, in which
the descent step size $\rho$ is determined using $10$ iterations of a dichotomy
method starting from $\rho_0 = 1$. This last layer of \text{SD} is carried
out in order to enhance the accuracy of the final solution.
This algorithm has been already tested for solving different computationally
expensive industrial optimization problems
(see, e.g., \cite{gomez2011,ivorra2013,ivorra2016,ivorra20182}).
Furthermore, it has been compared with other well-known metaheuristic
method and it exhibits better performance for a set of Benchmark problems
(see \cite{ivorra20132}). A Matlab version of the \text{GMA} presented
in this paper has been implemented in the free optimization package
``Global Optimization Platform'' (GOP), which can be downloaded
at http: //www.mat.ucm.es/momat/software.html.
\section{Numerical experiments}\label{sec:Results}
In this section, we first introduce the numerical solver used for computing
the solutions of system \eqref{eq:ADV}-\eqref{eq:NS}.
Then, in Section \ref{sec:NumCases} we describe the considered numerical
experiments based on the optimization problems \eqref{ProblemP1}, \eqref{ProblemP2}
and \eqref{ProblemP3}. Section \ref{sec:OptResults} presents the optimization
results, which are analyzed and compared in Section \ref{sec:ResultsComp}.
\subsection{Numerical implementation of the model}\label{sec:NumImplementation}
The cylindrical version of system \eqref{eq:ADV}-\eqref{eq:NS} was solved using
the software COMSOL Multiphysics 5.0 (http: //www.comsol.com),
based on the Finite Element method (see \cite{angel}), with Lagrange P2-P1
elements to stabilize the pressure and to satisfy the Ladyzhenskaya, Babouska,
and Brezzi stability condition. The 2nd-order Lagrange elements model the velocity
and concentration components, while linear elements represent the pressure.
At a first stage, we solve the stationary Navier-Stokes equations \eqref{eq:NS}
using Galerkin least square streamline and crosswind diffusion methods so as
to prevent numerical oscillations. At a second stage, the velocity field
(solution of \eqref{eq:NS}) is introduced as an input value in the transient
advection-diffusion-reaction system \eqref{eq:ADV}, which is then solved by
considering an upwind scheme. We use a direct damped Newton method to solve
the corresponding linear systems. A complete description of those techniques
can be found in \cite{glowinski2013}. The numerical experiments were carried
out in a 2.8Ghz Intel i7-930 64bits computer with 12Gb of RAM. We used a
triangular mesh with around 3000 elements, which produced significantly
accurate results with respect to finer meshes that turned out to be computationally
unreachable. We assumed that the solution of system \eqref{eq:ADV}
at finite time $\hat{T}= 10^{7}$ (s) could be considered as a reasonable
approximation the steady state $(\hat{S},\hat{B})$ of system \eqref{eq:ADV}.
Model variables \eqref{eq:Vol} and \eqref{eq:Sout} were estimated using the
functions \emph{Domain Integration} and \emph{Boundary Integration}
of COMSOL (based on a trapezoidal approximation of the integral),
respectively. Thus, the value of the cost function \eqref{eq:J} was an
output of the COMSOL model. After performing several numerical tests
(as described in Section \ref{sec:GeneralProblem},) we have chosen
$\beta=10^{9}$ to be the value of the penalty parameter in the cost function
as it has given good results for the considered numerical experiments.
In this work, the \text{GMA} has been applied with
$(s_{\rm max},N_{\rm g},N_{\rm p},p_{\rm c}, p_{\rm m}, \lambda, \hat{g})
= (100,10,10,0.4,0.2,1,25)$ for solving numerically problems
\eqref{ProblemP1}, \eqref{ProblemP2} and \eqref{ProblemP3}.
Those parameters had previously been successfully used for solving similar
optimization problems in \cite{MCrespoThesis,gomez2011,ivorra2013,ivorra2016}.
Depending on the considered case (detailed below), each function evaluation
in problems \eqref{ProblemP1}, \eqref{ProblemP2} and \eqref{ProblemP3} may
take from 15 up to 60 minutes. With a restriction of 3 months of computational
time to run the \text{GMA}, the number of function evaluations carried out
to solve problems \eqref{ProblemP1}, \eqref{ProblemP2} and \eqref{ProblemP3}
ranged between 2000 and 6000.
\subsection{Cases considered in this work}\label{sec:NumCases}
Model parameters were set as follows (see \cite{bellorivas,symon1971}):
$D_{\rm S}= 4.3\cdot 10^{-12}$(m$^2$/s), $D_{\rm B}= 5\cdot 10^{-10}$ (m$^2$/s),
$S_{\rm in}= 15$ (kg/m$^3$), $B_{0}= 1$ (kg/m$^3$), $S_{0}=15$ (kg/m$^3$),
$p_{\rm atm}= 10^{5}$ ({\rm Pa}), $\rho=10^{3}$ (kg/m$^3$),
$\eta=10^{-3}$ (kg/m s) and $u_{\rm in}=2.2\cdot 10^{-4}$ (m/s).
We consider four different reaction rate functions $\mu_1$, $\mu_2$, $\mu_3$
and $\mu_4$, which are described in Table \ref{table:reactions}
(see pages 132, 182 and 187 in \cite{Dochain2010}). In Figure \ref{fig:reactions},
we plot those four growth rate functions. We can observe that they have the
same order of magnitude but with different slopes.
\begin{figure}[ht]
\centering
\includegraphics[width=0.7\textwidth]{fig3}
\caption{Functions $\mu_1(S)$, $\mu_2(S)$, $\mu_{3}(S)$ and $\mu_{4}(S)$ (s$^{-1}$), detailed in Table \ref{table:reactions}, with $S \in [0,20]$ (kg/m$^3$).}
\label{fig:reactions}
\end{figure}
\begin{table}[ht]
\caption{Considered growth rate functions}
\label{table:reactions}
\centering
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{|c c |}
\hline
$\mu_1(\cdot)$ & $\mu_2(\cdot)$ \\ \hline
Monod function \eqref{eq:MONOD} & Monod function \eqref{eq:MONOD} \\
$\mu_{\rm max}= 9.17\cdot 10^{-5}$ s$^{-1}$
& $\mu_{\rm max}= 5.5\cdot 10^{-5}$ s$^{-1}$ \\
$K_{\rm S}=5$ kg/m$^3$ & $K_{\rm S}= 0.075$ kg/m$^3$ \\ \hline
\end{tabular}
\\
\begin{tabular}{|c c|}
\hline
$\mu_{3}(\cdot)$ & $\mu_{4}(\cdot)$\\ \hline
Haldane function \eqref{eq:HALDANE} & Haldane function \eqref{eq:HALDANE} \\
$\mu_{\rm max}=1.39\cdot 10^{-4}$ s$^{-1}$,
& $\mu_{\rm max}= 1.11\cdot 10^{-4}$ s$^{-1}$, \\
$K_{\rm S}=4$ kg/m$^3$, & $K_{\rm S}=0.5$ kg/m$^3$, \\
$K_{\rm I}= 3$ kg/m$^3$ & $K_{\rm I}=4$ kg/m$^3$ \\
\hline
\end{tabular}
\end{table}
When solving problems \eqref{ProblemP1} and \eqref{ProblemP2}, design
parameters $H_{\rm min} = 2$ (m), $H_{\rm max} = 10$ (m) and $R_{\rm max} = 5$ (m)
were taken to generate the admissible spaces $\tilde{\Phi}^{1}$ and
$\tilde{\Phi}^2$. On the other hand, when solving problem \eqref{ProblemP3},
the admissible space $\tilde{\Phi}^{3}$ was generated with design parameters
$H_{\rm min}=2$ (m) $H_{\rm max}=10$ (m) and $R_{1,{\rm max}} = R_{2,{\rm max}}= 3.5$
(m). In order to compute the values $E_{1,{\rm min}}$ and $E_{1,{\rm max}}$,
we chose $R_{\rm max}=5$ (m). In all cases we set $r=h=0.5$ (m) and
$S_{\rm lim}= 1$ (kg/m$^3$).
\subsection{Optimization results}\label{sec:OptResults}
In this section, we present the optimization results obtained when solving
problems \eqref{ProblemP1}, \eqref{ProblemP2} and \eqref{ProblemP3}.
Section \ref{sec:NumP1} puts together the numerical results obtained when
solving problems \eqref{ProblemP1} and \eqref{ProblemP2} because similar
conclusions were obtained. We point out that in all the optimized reactors
obtained in Sections \ref{sec:NumP1} and \ref{sec:NumP3},
the optimal solutions $\tilde{\phi}^{i,{\rm opt}}$ are such that the second
term in \eqref{eq:J} is zero, and therefore, the value $J(\tilde{\phi}^{i,opt})$
corresponds to the reactor volume
$\operatorname{Vol}(\tilde{\phi}^{i,{\rm opt}})$ $(i=1,2,3)$.
\begin{remark} \rm
We observe from Tables \ref{table:ResultsP1}-\ref{table:ResultsP3} that the value
$\hat{S}_{\rm out}(\tilde{\phi}^{i,{\rm opt}})$ $(i=1,2,3)$ is clearly smaller
than the prescribed value $S_{\rm lim}$ in all the considered cases, which
seems to indicate that smaller (and therefore better) domains could be obtained
with this value closer to $S_{\rm lim}$. This inaccuracy may be due to the lack
of numerical precision of the COMSOL model, which in turn is caused by the
restriction on the computational time (see Section \ref{sec:NumImplementation}
for more details). This fact highlights the difficulties tackled during the
numerical resolution of our optimization problem.
\end{remark}
\subsubsection{Parametrizations 1 and 2}\label{sec:NumP1}
Table \ref{table:ResultsP1} shows the optimal results when solving
problem \eqref{ProblemP1}, while the optimized shapes are depicted in
Figure \ref{fig:ResultsP1}. Similarly, Table \ref{table:ResultsP2}
shows the optimal results when solving problem \eqref{ProblemP2}, while
the optimized shapes are depicted in Figure \ref{fig:ResultsP2}.
\begin{table}[ht]
\caption{Value of the optimal parameters $(H^{\rm opt} ({\rm m})$ and
$R^{\rm opt} ({\rm m}))$ in $\tilde{\phi}^{1,{\rm opt}}$, solution of
\eqref{ProblemP1} with functions $\mu_{i}, i=1,\ldots,4$; outflow substrate
concentrations ($\hat{S}_{\rm out}(\tilde{\phi}^{1,{\rm opt}})$ (kg/ m$^3$);
and reactor volumes ($\operatorname{Vol}(\tilde{\phi}^{1,{\rm opt}})$
m$^3$).}
\label{table:ResultsP1}
\centering
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{| c c c c c |}
\hline
$\mu$ & $H^{\rm opt}$ & $R^{\rm opt}$
& $\hat{S}_{\rm out}(\tilde{\phi}^{1,{\rm opt}})$
& $\operatorname{Vol}(\tilde{\phi}^{1,{\rm opt}})$ \\ \hline
$\mu_1$ & $10$ & $0.61$ & $0.9528$ & $11.3063$\\ \hline
$\mu_2$ & $9.82$ & $0.66$ & $0.9773 $ & $12.8553$ \\ \hline
$\mu_{3}$ & $9.92$ & $1.98$ & $0.9718$ & $ 110.6468$ \\ \hline
$\mu_{4}$ & $9.9$ & $1.6$ & $0.9024$ & $72.3634$ \\ \hline
\end{tabular}
\end{table}
\begin{table}[ht]
\caption{Value of the optimal parameters $(H^{\rm opt} ({\rm m})$ and
$R^{\rm opt}$ (m) in $\tilde{\phi}^{2,{\rm opt}}$, solution of
\eqref{ProblemP2} with functions $\mu_{i}, i=1,\ldots,4$;
outflow substrate concentrations ($\hat{S}_{\rm out}(\tilde{\phi}^{2,{\rm opt}})$
(kg/ m$^3$); and reactor volumes ($\operatorname{Vol}(\tilde{\phi}^{2,{\rm opt}})$
(m$^3$).}
\label{table:ResultsP2}
\centering
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{| c c c c c |}
\hline
$\mu$ & $H^{\rm opt}$ & $R^{\rm opt}$
& $\hat{S}_{\rm out}(\tilde{\phi}^{2,{\rm opt}})$
& $\operatorname{Vol}(\tilde{\phi}^{2,{\rm opt}})$ \\ \hline
$\mu_1$ & $9.28$ & $0.728$ & $0.9893$ & $12.8479$\\ \hline
$\mu_2$ & $9.91$ & $0.73$ & $0.9558$ & $13.8268$ \\ \hline
$\mu_{3}$ & $9.98$ & $1.96$ & $0.9538$ & $80.2715$ \\ \hline
$\mu_{4}$ & $10$ & $1.36$ & $0.9216$ & $40.6650$ \\ \hline
\end{tabular}
\end{table}
\begin{figure}[ht]
\centering
\subfigure[$\mu_1$]{ \includegraphics[width=0.23\textwidth]{fig4a}}
\subfigure[$\mu_2$]{ \includegraphics[width=0.23\textwidth]{fig4b}}
\subfigure[$\mu_{3}$]{\includegraphics[width=0.23\textwidth]{fig4c}}
\subfigure[$\mu_{4}$]{ \includegraphics[width=0.23\textwidth]{fig4d}}
\caption{Shape of the optimized reactors, $\Omega(\tilde{\phi}^{1,{\rm opt}})$, where $\tilde{\phi}^{1,{\rm opt}}$ is the solution of problem \eqref{ProblemP1}.}
\label{fig:ResultsP1}
\end{figure}
\begin{figure}[ht]
\centering
\subfigure[$\mu_1$]{ \includegraphics[width=0.23\textwidth]{fig5a}}
\subfigure[$\mu_2$]{ \includegraphics[width=0.23\textwidth]{fig5b}}
\subfigure[$\mu_{3}$]{\includegraphics[width=0.23\textwidth]{fig5c}}
\subfigure[$\mu_{4}$]{\includegraphics[width=0.23\textwidth]{fig5d}}
\caption{Shape of the optimized reactors, $\Omega(\tilde{\phi}^{2,{\rm opt}})$,
where $\tilde{\phi}^{2,{\rm opt}}$ is the solution of problem \eqref{ProblemP2}.}
\label{fig:ResultsP2}
\end{figure}
From Figures \ref{fig:ResultsP1} and \ref{fig:ResultsP2} we observe that the
optimal reactors have height larger than width. This outcome is in line with
the results found in \cite{Richmond1993,Tredici1998}, where the authors
performed experimental studies to conclude that the most efficient reactor
was a tubular one with its height much greater than its radius.
Nevertheless, this strategy is not always applicable due to the practical
restriction on the reactor height.
When comparing the results obtained with the reaction functions,
we observe that the optimized reactors exhibit similar heights
and the main difference lies in the reactor radius. For instance,
the value of $\operatorname{Vol}(\tilde{\phi}^{1,{\rm opt}})$
(similarly, the value of $\operatorname{Vol}(\tilde{\phi}^{2,{\rm opt}})$)
is higher with $\mu_3$ than with $\mu_1$. This difference seems to be due to
the fact that function $\mu_3$ is qualitatively smaller than function
$\mu_1$ (see Figure \ref{fig:reactions}) and thus, the optimal volume must
be bigger to ensure that the prescribed value $S_{\rm lim}$ is reached.
The influence of the reactor width on the bioreactor dynamics is explained
in Remarks \ref{rmk:widthflow} and \ref{rmk:storage} later on.
\subsubsection{Parametrization 3}\label{sec:NumP3}
Table \ref{table:ResultsP3} shows the optimal results, while the optimized
shapes are depicted in Figure \ref{fig:ResultsP3}.
\begin{table}[ht]
\caption{Value of the optimal parameters
$(H^{\rm opt} ({\rm m}), R_1^{\rm opt} ({\rm m}), R_2^{\rm opt} ({\rm m}),
\alpha_1^{\rm opt}$ and $\alpha_2^{\rm opt}$) in $\tilde{\phi}^{3,{\rm opt}}$,
solution of \eqref{ProblemP3} with functions $\mu_{i}, i=1,\ldots,4$;
exterior control point coordinates ($E_1$ (m) and $E_2$ (m)), associated
with $\tilde{\phi}^{3,{\rm opt}}$ and computed using equation \eqref{eq:alphaopt}
and Lemma \ref{lem:cuentasB1B2}; outflow substrate concentrations
($\hat{\rm S}_{\rm out}(\tilde{\phi}^{3,{\rm opt}}) ({\rm kg}/{\rm m}^3))$;
and reactor volumes ($\operatorname{Vol}(\tilde{\phi}^{3,{\rm opt}}) ({\rm m}^3)$).}
\label{table:ResultsP3}
\centering
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{| c c c c c c |}\hline
$\mu$ & $H^{\rm opt}$ & $R_1^{\rm opt}$ & $R_2^{\rm opt}$
& $\alpha_1^{\rm opt}$ & $\alpha_2^{\rm opt}$ \\ \hline
$\mu_1$ & $9.6359$ & $0.5931$ & $0.9214$ & $1.6358\cdot 10^{-5}$ & $0.0059$ \\
\hline
$\mu_2$ & $9.0814$ & $0.5730$ & $1.1024$ & $0.0064$ & $0.0922$ \\
\hline
$\mu_{3}$ & $9.8879$ & $2.0432$ & $1.2402$ & $0.0093$ & $0.9940$ \\
\hline
$\mu_{4}$ & $9.6401$ & $1.0410$ & $2.0805$ & $1.9301\cdot 10^{-4}$ & $0.0273$ \\
\hline
\end{tabular}
\begin{tabular}{| c c c c c |}\hline
$\mu$ & $ E_1$ & $E_2$
& $\hat{\rm S}_{\rm out}(\tilde{\phi}^{3,{\rm opt}})$
& $\operatorname{Vol}(\tilde{\phi}^{3,{\rm opt}})$ \\
\hline
$\mu_1$ & $0.3022$ & $0.5512$ & $0.9540$ & $10.0790$ \\
\hline
$\mu_2$ & $0.4111$ & $1.2451$ & $0.8923$ & $12.0855$ \\
\hline
$\mu_{3}$ & $-0.3928$ & $8.8376$ & $0.9545$ & $27.1679$ \\
\hline
$\mu_{4}$ & $-0.421$ & $0.736$ & $0.9962$ & $22.3042$ \\
\hline
\end{tabular}
\end{table}
\begin{figure}[ht]
\centering
\subfigure[$\mu_1$]{\includegraphics[width=0.23\textwidth]{fig6a}}
\subfigure[$\mu_2$]{\includegraphics[width=0.23\textwidth]{fig6b}}
\subfigure[$\mu_{3}$]{\includegraphics[width=0.23\textwidth]{fig6c}}
\subfigure[$\mu_{4}$]{\includegraphics[width=0.23\textwidth]{fig6d}}
\caption{Shape of the optimized reactors, $\Omega(\tilde{\phi}^{3,{\rm opt}})$,
where $\tilde{\phi}^{3,{\rm opt}}$ is the solution of problem \eqref{ProblemP3}.}
\label{fig:ResultsP3}
\end{figure}
From Figure \ref{fig:ResultsP3} we observe that, as stated in Section
\ref{sec:NumP1}, the optimal reactors have height larger than width.
Moreover, the exterior wall of the optimized reactors is concave,
with a radius expansion observed at least in some limited part of the
reactor (as said previously, the influence of the reactor radius in the
bioreactor dynamics will be explained in Remarks \ref{rmk:widthflow} and
\ref{rmk:storage} below).
When comparing reaction functions, Figures \ref{fig:ResultsP3}-(a)
and \ref{fig:ResultsP3}-(c) seem to show that, for instance, the
radius expansion of the domain $\Omega(\tilde{\phi}^{3,{\rm opt}})$
is wider for growth rate function $\mu_3$ than for $\mu_1$, as observed
in Section \ref{sec:NumP1}. On the other hand, the main difference
between considering Monod ($\mu_1$ and $\mu_2$) or Haldane ($\mu_3$ and $\mu_4$)
reaction functions is observed in the concavity at the upper part of the
exterior wall (see Remark \ref{rmk:widthflow} for a physical interpretation).
\subsection{Comparison between the optimized reactors}\label{sec:ResultsComp}
Here, we compare the solutions obtained when solving the optimization
problems \eqref{ProblemP1}, \eqref{ProblemP2} and \eqref{ProblemP3}.
Figures \ref{fig:ResultsP1}, \ref{fig:ResultsP2} and \ref{fig:ResultsP3}
seem to show that the optimized reactors have height larger than width
(indeed, $H_{\rm max}$ set to 10 (m) limits the optimal shape height and the
optimal heights in all the considered cases approach this limit) and generally,
the optimal widths approach its lower bound (the minimum reactor radius
allowed was $r=0.5$ (m)). However, in some of the considered cases
(see Figures \ref{fig:ResultsP1}-(c), \ref{fig:ResultsP1}-(d),
\ref{fig:ResultsP2}-(c), \ref{fig:ResultsP2}-(d), \ref{fig:ResultsP3}-(c)
and \ref{fig:ResultsP3}-(d)), a radius expansion (at least in some
limited part of the reactor) is observed. We interpret that increasing
the reactor width favors the reaction due to two main reasons:
\begin{enumerate}
\item It helps that the vertical flow velocity decreases (in absolute value),
and so the time that the biomass and the substrate remain in contact for
reacting increases (see Remark \ref{rmk:widthflow} for a more detailed
analysis of the relation between the reactor width and the vertical flow).
\item It originates an area of biomass storage. For example, due to the
apparition of Dean vortices in this area (see, e.g., \cite{Dean1928})
the biomass located near the device exterior wall remains more time inside
the bioreactor (compared to the biomass located at the reactor center),
and so the amount of reaction between biomass and substrate increases
(see Remark \ref{rmk:storage} for an specific explanation about the
distribution of substances in the reactor).
\end{enumerate}
\begin{remark}\label{rmk:widthflow} \rm
To understand the influence of the bioreactor width on the vertical flow
velocity, we used four different domains, denoted by $\Omega_{i}$, $i=1,\ldots,4$.
The first reactor is cylindrical (depicted in Figure \ref{fig:velocity}-(a))
and the other three present a radius extension on the top, center and bottom
parts of the domain, (depicted in Figures \ref{fig:velocity}-(b)
to \ref{fig:velocity}-(d), respectively).
We solved system \eqref{eq:NS} with domains $\Omega_{i}$, $i=1,\ldots,4$
and denoted $u_{3,\Omega_{i}}$ (m/s) the vertical flow velocity
(third component of the velocity vector) obtained when solving system
\eqref{eq:NS} in the domain $\Omega_{i}$, evaluated at $r=0$ (i.e.,
symmetry streamline). Figure \ref{fig:velocity}-(e) represents $|u_{3,\Omega_{i}}|$,
$i\in \{1,\ldots,4\}$, which can be seen as functions of $z$.
We observe that, in regions where the reactor radius increases, the absolute
value of the vertical velocity decreases. This physical interpretation may
explain, for instance, the optimal domains $\Omega(\tilde{\phi}^{3,{\rm opt}})$
obtained with reactions $\mu_3$ and $\mu_4$ (see Figures \ref{fig:ResultsP3}-(c)
and \ref{fig:ResultsP3}-(d)), since the Haldane function shows inhibition for
large values of substrate (see Figure \ref{fig:reactions}) and the maximum
value of substrate appears at the reactor inlet.
\begin{figure}[ht]
\centering
\subfigure[$\Omega_1$]{\includegraphics[width=0.17\textwidth]{fig7a}}
\subfigure[$\Omega_2$]{\includegraphics[width=0.23\textwidth]{fig7b}}
\subfigure[$\Omega_{3}$]{\includegraphics[width=0.23\textwidth]{fig7c}}
\subfigure[$\Omega_{4}$]{\includegraphics[width=0.23\textwidth]{fig7d}}
\subfigure[Vertical velocity profile along the symmetry streamline obtained
for reactor domains $\Omega_{i}$, $i=1,\ldots,4$.]
{\includegraphics[width=0.6\textwidth]{fig7e}}
\caption{Influence of the reactor width into the vertical flow velocity.}
\label{fig:velocity}
\end{figure}
\end{remark}
\begin{remark}\label{rmk:storage} \rm
Figures \ref{fig:storageh2}-(a) and \ref{fig:storageh2}-(b) represent
the distributions of substrate and biomass at steady state, respectively,
computed with the optimal reactor $\Omega(\tilde{\phi}^{3,{\rm opt}})$,
obtained for the growth rate function $\mu_{3}$. One observes that
the substrate is mainly agglomerated in the area originated by the inlet
streamlines (see Figure \ref{fig:storageh2}-(d)). On the other hand,
the biomass becomes withdrawn from this central area and is mainly concentrated
around the reactor wall (where Dean vortices appear \cite{Dean1928}, see
Figure \ref{fig:storageh2}-(d)). Thus, a reaction front is created between
the central area and the outer part of the reactor (as shown in
Figure \ref{fig:storageh2}-(c)) favoring the reaction between the two species.
Although the optimization problems \eqref{ProblemP1}, \eqref{ProblemP2}
and \eqref{ProblemP3} have been solved for a singular pair of diffusion
coefficients $(D_{\rm S},D_{\rm B})$, numerical experiments seem to show that
the analysis of the distribution of substances in the reactor, performed above,
is suitable in the range of typical diffusion coefficients $D_{\rm S}$
(from $10^{-10}$ to $10^{-7}$ (m$^2$/s) \cite{Nogueira2006,Stewart2003,Venterea2000})
and $D_{\rm B}$ (from $10^{-13}$ to $10^{-7}$ (m$^2$/s)
\cite{Harding1986,Reiter1999,singh2013}).
\end{remark}
\begin{figure}[ht]
\centering
\subfigure[$\hat{S}(r,z)$ (kg/m$^3$)]{\includegraphics[width=0.23\textwidth]{fig8a}}
\subfigure[$\hat{B}(r,z)$ (kg/m$^3$)]{\includegraphics[width=0.23\textwidth]{fig8b}}
\subfigure[$\mu(\hat{S})\hat{B}$ (kg/m$^3$ h)]{\includegraphics[width=0.23\textwidth]{fig8c}}
\subfigure[Streamlines]{\includegraphics[width=0.19\textwidth]{fig8d}}
\caption{(a) substrate concentration (at steady state).
(b) biomass concentration (at steady state). (c) reaction (at steady state)
(d) streamlines. (a)-(d) associated with the optimal reactor
$\Omega(\phi^{3,{\rm opt}})$ obtained for the growth rate function $\mu_3$.}
\label{fig:storageh2}
\end{figure}
Now, it is of interest to compare the optimized reactors obtained when solving
the optimization problems \eqref{ProblemP1}, \eqref{ProblemP2} and \eqref{ProblemP3}.
In this direction, we denote by $d_{r}(\tilde{\phi}^{i},\tilde{\phi}^{j})$
($i\neq j$) the relative difference between the optimal reactor volumes
$\operatorname{Vol}(\tilde{\phi}^{i,{\rm opt}})$ and
$\operatorname{Vol}(\tilde{\phi}^{j,{\rm opt}})$, defined as
\begin{equation}\label{eq:dr}
d_{r}(\tilde{\phi}^{i},\tilde{\phi}^{j})= 100 \times
\frac{\operatorname{Vol}(\tilde{\phi}^{j,{\rm opt}})
-\operatorname{Vol}(\tilde{\phi}^{i,{\rm opt}})}{\operatorname{Vol}
(\tilde{\phi}^{i,{\rm opt}})}.
\end{equation}
Table \ref{table:comp} shows the comparison, in terms of reactor volume,
between the optimized reactors obtained when creating the domain with the
three proposed parametrizations. Additionally, values
$d_{r}(\tilde{\phi^{1}},\tilde{\phi^2})$ (resp.
$d_{r}(\tilde{\phi^{1}},\tilde{\phi^{3}})$) are included in Table \ref{table:comp}
in order to outline how much can be gained by using non-tubular reactors.
One observes that for Monod growth rate functions ($\mu_1$ and $\mu_2$),
the relative difference $d_{r}(\tilde{\phi}^{1},\tilde{\phi}^{j})$ $(j=2,3)$
is between $-10\%$ and $13\%$, so one can conclude that the variation
(in terms or reactor volume) between tubular and non-tubular reactors is
relatively low. On the other hand, for Haldane growth rate functions
($\mu_3$ and $\mu_4$), one observes that by using non-tubular reactors, one can
gain from $27\%$ up to $75\%$ in terms of reactor volume. Additionally,
one may notice that the optimized reactor volumes differ substantially depending
on the considered reaction function. Similar results have been obtained in
other works focusing on the analysis of continuous bioreactors modeled
through ordinary differential equations (see, e.g.\
\cite{dochain2001,Harmand2005,Harmand2003,Hill1989,Luyben1982})).
\begin{table}[ht]
\caption{Comparison, in terms of reactor volume $({\rm m}^3)$,
between optimized reactors obtained when solving problems \eqref{ProblemP1}
$(\operatorname{Vol}(\tilde{\phi}^{1,{\rm opt}}))$, \eqref{ProblemP2}
$(\operatorname{Vol}(\tilde{\phi}^{2,{\rm opt}}))$ and \eqref{ProblemP3}
$(\operatorname{Vol}(\tilde{\phi}^{3,{\rm opt}}))$; relative differences
$d_{r}(\tilde{\phi}^{1},\tilde{\phi}^{j})$ $(\%)$ ($j=2,3$) computed using
equation \eqref{eq:dr}.}
\label{table:comp}
\centering
\renewcommand{\arraystretch}{1.3}
\begin{tabular}{| c c c c |}
\hline
$\mu$ & $\operatorname{Vol}(\tilde{\phi}^{1,{\rm opt}})\text{ }({\rm m}^3)$
& $\operatorname{Vol}(\tilde{\phi}^{2,{\rm opt}})\text{ }({\rm m}^3)$
& $d_{r}(\tilde{\phi}^{1},\tilde{\phi}^2)\text{ }(\%)$ \\
\hline
$\mu_1$ & $11.3063$ & $12.8479$ & $+ 13.63$ \\
\hline
$\mu_2$ & $12.8553$ & $13.8268$ & $+ 7.56$ \\
\hline
$\mu_{3}$ &$ 110.6468$ & $80.2715$ & $- 27.45$ \\
\hline
$\mu_{4} $& $72.3334$ & $40.6650$ & $- 43.78$ \\
\hline
\end{tabular}
\\
\begin{tabular}{| c c c |}
\hline
$\mu$ & $\operatorname{Vol}(\tilde{\phi}^{3,{\rm opt}})\text{ }({\rm m}^3)$
& $d_{r}(\tilde{\phi}^{1},\tilde{\phi}^{3})\text{ }(\%)$ \\
\hline
$\mu_1$ & $ 10.0790 $ &$ - 10.85$\\
\hline
$\mu_2$ & $12.0855$ & $- 5.98$\\
\hline
$\mu_{3}$ & $27.1679$ & $-75.45$ \\
\hline
$\mu_{4} $& $22.3042$ & $-69.17$ \\
\hline
\end{tabular}
\end{table}
\section{Conclusions}\label{sec:Conclusions}
We have explored the shape design of a continuous biological reactor.
The main objective was to reduce the reactor volume, ensuring that a
prescribed output concentration value was reached at stationary state.
As a matter of generalization, we have not imposed a particular type of
biological dynamics in the reactor but proposed a general methodology
to be applied and adapted depending on the considered case. We have used a
mathematical model that couples hydrodynamics (described with the
incompressible Navier--Stokes equations in three dimensions) with biological
phenomena (described with an Advection-Diffusion-Reaction system).
Using the Finite Element Method, we have numerically computed the output
substrate concentration at steady state and the volume of a reactor associated
with a particular set of design parameters. Then, we have defined three discrete
optimization problems related to three different parametrizations of the device
geometry and solved them by using a Genetic Multi-layer Algorithm, a
self-implemented global optimization method based on the search of a suitable
initial population for a given genetic algorithm. The optimization problem
is solved for monotonic and non-monotonic growth rate functions, in order to
analyze the influence of the reaction into the optimal reactor configuration.
One of the proposed parametrizations allows us to model tubular shapes
(typically used in the industry sector), while the other two offer the possibility
to obtain a wider range of reactor geometries. We have taken into account that
the reaction between species may be modeled by either monotonic or non-monotonic
growth functions, and we have analyzed the influence of this factor on the
optimization results.
From a general point of view, the optimized reactors exhibit height much larger
than width and their exterior wall is concavely curved, with a radius expansion
observed at least in some limited part of the reactor. The magnitude of the
radius extension appears to be related to the reaction function.
The slower is the reaction the wider should be the device. The advantage
of the radius extensions in the reactor performance could be attributed to
two main factors:
\begin{itemize}
\item The width of the reactor helps to decrease the absolute value of the
vertical flow velocity, and consequently, increases the time of potential
reaction between substances.
\item The reactor corners may act as a biomass storage. The biomass located
near the reactor exterior wall is withdrawn slower from the device than the
biomass located near the device center, favoring the reaction between species.
\end{itemize}
When comparing the optimized reactors obtained for both monotonic and
non-monotonic growth rate functions, one observes that, if the reactor is
modeled with non-monotonic kinetics (e.g., Haldane function), the radius
expansion located at the top of the reactor is more pronounced. We interpret
that this difference relays on the fact that for large values of substrate
(i.e., at the inlet of the device) the Haldane reaction shows inhibition, and
so the radius expansion should be bigger to decrease the absolute value of
the vertical flow near the reactor inlet.
An interesting message is that when the degrees of freedom on the shape
parame\-trization is large enough, concavely curved reactors are systematically
the best ones, although convexly curved (generated with less degrees of freedom)
could already be better than perfect tubular shapes. The volume of the
optimized reactors using concavely curved geometries (instead of typical tubular
geometries) allow us to reduce the reactor volume between $25\%$ and $90\%$
depending on the considered case. This study shows to practitioners how much
could be gained compared to the best cylindrical shapes, and how this gain is
related to the reaction rate function (monotonic versus non-monotonic).
Of course, the economic cost relative to the production of non-conventional
should be taken into consideration for the best choice of the reactor geometry.
Including this new (or other) criteria into the design problem should be tackled
by using multi-objective approaches instead of single-objective approaches
(see \cite{Ferrandez2018}).
\subsection*{Acknowledgments}
This work was carried out thanks to the financial support of the Spanish
``Ministry of Economy and Competitiveness'' under project MTM2015-64865-P;
the research group MOMAT (Ref. 910480) supported by ``Universidad Complutense
de Madrid''; and the ``Junta de Andaluc\'ia'' and the European Regional
Development Fund through project P12-TIC301. The authors thank the French
LabEx Numev (convention ANR-10-LABX-20) for the postdoctoral grant of the
First Author at UMR MISTEA, Montpellier, France.
\begin{thebibliography}{10}
\bibitem{Andrews} J. F. Andrews;
\emph{A mathematical model for the continuous culture of
microorganisms utilizing inhibitory substrates}, Biotechnol. Bioeng.,
\textbf{10} (1968), no.~6, 707--723.
\bibitem{Ansoni2016} J. L. Ansoni, P.~Seleghim;
\emph{Optimal industrial reactor design:
development of a multiobjective optimization method based on a posteriori
performance parameters calculated from {CFD} flow solutions}, Adv. Eng.
Softw., \textbf{91} (2016), no.~C, 23--35.
\bibitem{Jameson1988} J.~Antony;
\emph{Aerodynamic design via control theory}, J. Sci. Comput.,
\textbf{3} (1988), no.~3, 233--260.
\bibitem{bellorivas} J. M. Bello, B.~Ivorra, A. M. Ramos, A.~Rapaport, J.~Harmand;
\emph{Bioreactor shape optimisation. modeling, simulation, and shape
optimization of a simple bioreactor for water treatment}, Les STIC pour
l{'}environnement 2011, Lavoisier, 2011, pp.~125--141.
\bibitem{Bitog2011} J. P. Bitog, I.-B. Lee, C.-G. Lee, K.-S. Kim,
H.-S. Hwang, S.-W. Hong, I.-H. Seo, K.-S. Kwon, E.~Mostafa;
\emph{Application of computational fluid
dynamics for modeling and designing photobioreactors for microalgae
production: A review}, Comput. Electron. Agr., \textbf{76} (2011), no.~2,
131--147.
\bibitem{Campana2006} E. F. Campana, D.~Peri, Y.~Tahara, F.~Stern;
\emph{Shape optimization in
ship hydrodynamics using computational fluid dynamics}, Comput. Methods Appl.
Mech. Eng., \textbf{196} (2006), no.~1--3, 634--651.
\bibitem{Carvalho2006} A. P. Carvalho, L. A. Meireles, F.X. Malcata;
\emph{Microalgal reactors: A
review of enclosed system designs and performances}, Biotechnol. Prog.,
\textbf{22} (2006), no.~6, 1490--1506.
\bibitem{Coenen2013} T.~Coenen, W. V. Moortel, F.~Logist, J.~Luyten,
J. F. M. van Impe, J.~Degr\`eve;
\emph{Modeling and geometry optimization of photochemical
reactors: Single- and multi-lamp reactors for {UV–H}$_2${O}$_2$ {AOP}
systems}, Chem. Eng. Sci., \textbf{96} (2013), 174--189.
\bibitem{MCrespoThesis} M.~Crespo;
\emph{Mathematical modeling and optimization of bioreactors and
liquid crystals}, Ph.D. thesis, Universidad Complutense de Madrid, 2017.
\bibitem{Crespo2015} M.~Crespo, B.~Ivorra, A. M. Ramos;
\emph{Existence and uniqueness of
solution of a continuous flow bioreactor model with two species}, RACSAM A:
Mat., \textbf{110} (2016), no.~2, 357--377.
\bibitem{Crespo2017} M.~Crespo, B.~Ivorra, A. M. Ramos;
\emph{Asymptotic stability of a coupled advection-diffusion-reaction system
arising in bioreactor processes}, Electr. J. Differ. Equ., \textbf{2017} (2017),
no.~194, 1--26.
\bibitem{Crespo2016} M.~Crespo, A. M. Ramos, B.~Ivorra, A.~Rapaport;
\emph{Modeling and
optimization of activated sludge bioreactors for wastewater treatment taking
into account spatial inhomogeneities}, J. Process Contr. ,\textbf{54} (2017),
118--128.
\bibitem{Dean1928} W. R. Dean;
\emph{The stream-line motion of fluid in a curved pipe (2$^{\rm nd}$
paper)}, Lond. Edinb. Dubl. Phil. Mag. Sci., \textbf{5} (1928), no.~30,
673--695.
\bibitem{deb2001} K.~Deb;
\emph{Multi-objective optimization using evolutionary algorithms},
Wiley Interscience Series in Systems and Optimization, Wiley, 2001.
\bibitem{diaz2017} J. I. D\'iaz;
\emph{Remarks on the preceding paper by {C}respo, {I}vorra and
{R}amos on the stability of bioreactor processes}, Electr. J. Differ. Equ.,
\textbf{2017} (2017), no.~195, 1--25.
\bibitem{Dochain2010} D.~Dochain;
\emph{Automatic control of bioprocesses}, ISTE, Wiley, 2010.
\bibitem{dochain2001} D.~Dochain, P.~Vanrolleghem;
\emph{Dynamical modelling and estimation in
wastewater treatment processes}, IWA Publishing, 2001.
\bibitem{ElSayed2005} M.~El-Sayed, T.~Sun, J.~Berry;
\emph{Shape optimization with computational
fluid dynamics}, Adv. Eng. Softw., \textbf{36} (2005), no.~9, 607--613.
\bibitem{farin2014} G.~Farin;
\emph{Curves and surfaces for computer-aided geometric design: A
practical guide}, Computer science and scientific computing, Elsevier
Science, 2014.
\bibitem{Ferrandez2018} M. R. Ferr{\'a}ndez, S.~Puertas-Mart\'in,
J. L. Redondo, B.~Ivorra, A. M. Ramos, P. M. Ortigosa;
\emph{High-performance computing for the optimization of
high-pressure thermal treatments in food industry}, J. Supercomput., (2018).
\bibitem{Gajardo2008} P.~Gajardo, H.~Ram\'irez, A.~Rapaport;
\emph{Minimal time sequential batch
reactors with bounded and impulse controls for one or more species}, SIAM J.
Control Optim., \textbf{47} (2008), no.~6, 2827--2856.
\bibitem{glowinski2013} R.~Glowinski;
\emph{Numerical methods for nonlinear variational problems},
Scientific Computation, Springer Berlin Heidelberg, 2013.
\bibitem{goldberg1989} D. E. Goldberg;
\emph{Genetic algorithms in search, optimization, and machine
learning}, Artificial Intelligence, Addison-Wesley, 1989.
\bibitem{gomez2011} S.~Gomez, B.~Ivorra, A. M. Ramos;
\emph{Optimization of a pumping ship
trajectory to clean oil contamination in the open sea.}, Math. Comput. Model,
\textbf{54} (2011), no.~1, 477--489.
\bibitem{Harding1986} S. E. Harding, P.~Johnson;
\emph{The brownian diffusion of dormant and
germinating spores of bacillus megaterium}, J. Appl. Bacteriol., \textbf{60}
(1986), no.~3, 227--232.
\bibitem{Harmand2005} J. Harmand, D.~Dochain;
\emph{The optimal design of two interconnected
(bio)chemical reactors revisited}, Comput. Chem. Eng., \textbf{30} (2005),
no.~1, 70--82.
\bibitem{BookAlain} J.~Harmand, C.~Lobry, A.~Rapaport, T.~Sari;
\emph{The chemostat:
Mathematical theory of microorganism cultures}, Wiley, 2017.
\bibitem{Harmand2003} J.~Harmand, A.~Rapaport, A.~Trofino;
\emph{Optimal design of interconnected
bioreactors: New results}, AIChE J., \textbf{49} (2003), no.~6, 1433--1450.
\bibitem{Hill1989} G. A. Hill, C. W. Robinson;
\emph{Minimum tank volumes for {CFST} bioreactors
in series}, Can. J. Chem. Eng., \textbf{67} (1989), no.~5, 818--824.
\bibitem{Isebe2008} D.~Isebe, P.~Azerad, F.~Bouchette, B.~Ivorra, B.~Mohammadi;
\emph{Shape optimization of geotextile tubes for sandy beach protection},
Int. J. Numer. Meth. Eng., \textbf{74} (2008), no.~8, 1262--1277.
\bibitem{ivorra2018} B.~Ivorra;
\emph{Application of the laminar navier--stokes equations for
solving 2d and 3d pathfinding problems with static and dynamic spatial
constraints: Implementation and validation in comsol multiphysics}, J. Sci.
Comput., \textbf{74} (2018), no.~2, 1163--1187.
\bibitem{ivorra20182} B.~Ivorra, M. R. Ferr\'andez, M.~Crespo, J. L. Redondo,
P. M. Ortigosa, J. G. Santiago, A.M. Ramos;
\emph{Modelling and optimization applied to the
design of fast hydrodynamic focusing microfluidic mixer for protein folding},
J. Math. Industry, \textbf{8} (2018), no.~4.
\bibitem{ivorra2016} B.~Ivorra, J.~L\'opez, A. M. Ramos, J.G. Santiago;
\emph{Design sensitivity
and mixing uniformity of a micro-fluidic mixer}, Phys. Fluids (1994-present)
\textbf{28} (2016), no.~1, 012005.
\bibitem{ivorra20062} B.~Ivorra, B.~Mohammadi, L.~Dumas, O.~Durand, P.~Redont;
\emph{Semi-deterministic vs. genetic algorithms for global optimization of
multichannel optical filters}, Int. J. Comput. Sci. Eng., \textbf{2} (2006),
no.~3, 170--178.
\bibitem{ivorra2009} B.~Ivorra, B.~Mohammadi, A. M. Ramos;
\emph{Optimization strategies in credit portfolio management},
J. Global Optim., \textbf{42} (2009), no.~2, 415--427.
\bibitem{ivorra20132} B.~Ivorra, B.~Mohammadi, A. M. Ramos;
\emph{A multi-layer line search method to improve the initialization
of optimization algorithms}, Eur. J. Oper. Res., \textbf{247} (2015), no.~3,
711--720.
\bibitem{ivorra2006} B.~Ivorra, B.~Mohammadi, J. G. Santiago, D. E. Hertzog;
\emph{Semi-deterministic and genetic algorithms for global optimization of
microfluidic protein folding devices}, Int. J. Numer. Meth. Eng., \textbf{66}
(2006), no.~2, 319--333.
\bibitem{ivorra2007} B.~Ivorra, A. M. Ramos, B.~Mohammadi;
\emph{Semideterministic global optimization method: Application to a
control problem of the burgers equation}, J. Optimiz. Theory Appl.,
\textbf{135} (2007), no.~3, 549--561.
\bibitem{ivorra2013} B.~Ivorra, J. L. Redondo, J. G. Santiago, P. M. Ortigosa,
A. M. Ramos;
\emph{Two- and three-dimensional modeling and optimization applied to the
design of a fast hydrodynamic focusing microfluidic mixer for protein
folding}, Phys. Fluids (1994-present), \textbf{25} (2013), no.~3, 032001.
\bibitem{leveque2007} R. J. LeVeque;
\emph{Finite difference methods for ordinary and partial
differential equations: Steady-state and time-dependent problems}, Society
for Industrial and Applied Mathematics, 2007.
\bibitem{Luyben1982} K. Ch. A. M. Luyben, J.~Tramper;
\emph{Optimal design for continuous stirred
tank reactors in series using {M}ichaelis {M}enten kinetics}, Biotechnol.
Bioeng., \textbf{24} (1982), no.~5, 1217--1220.
\bibitem{Raino1999} R. A. E. M\"akinen, J.~Periaux, J.~Toivanen;
\emph{Multidisciplinary shape
optimization in aerodynamics and electromagnetics using genetic algorithms},
Int. J. Numer. Methods Fluids, \textbf{30} (1999), no.~2, 149--159.
\bibitem{Muyl2004} F.~Muyl, L.~Dumas, V.~Herbert;
\emph{Hybrid method for aerodynamic shape
optimization in automotive industry}, Comput. Fluids, \textbf{33} (2004),
no.~5-6, 849--858.
\bibitem{nejati2003} V.~Nejati, K.~Matsuuchi;
\emph{Aerodynamics {D}esign and {G}enetic
{A}lgorithms for {O}ptimization of {A}irship {B}odies}, JSME Int. J. Ser. B
Fluids Therm. Eng., \textbf{46} (2003), no.~4, 610--617.
\bibitem{nocedal2006} J.~Nocedal, S.J. Wright;
\emph{Numerical optimization, 2$^{\rm nd}$ ed},
Springer New York, USA, 2006.
\bibitem{Nogueira2006} R.~Nogueira, L. F. Melo;
\emph{Competition between nitrospira spp. and
nitrobacter spp. in nitrite-oxidizing bioreactors}, Biotechnol. Bioeng.,
\textbf{95} (2006), no.~1, 169--175.
\bibitem{pironneau1984} O.~Pironneau;
\emph{Optimal shape design for elliptic systems}, Springer series
in computational physics, Springer-Verlag, 1984.
\bibitem{Poloni2000} C.~Poloni, A.~Giurgevich, L.~Onesti, V.~Pediroda;
\emph{Hybridization of a multi-objective genetic algorithm, a neural network
and a classical optimizer for a complex design problem in fluid dynamics},
Comput. Methods in Appl. Mech. Eng., \textbf{186} (2000), no.~2--4, 403--420.
\bibitem{Pramparo2008} L.~Pramparo, J.~Pruvost, F.~St\"uber, J.~Font,
A.~Fortuny, A.~Fabregat, P.~Legentilhomme, J.~Legrand, C.~Bengoa;
\emph{Mixing and hydrodynamics
investigation using {CFD} in a square-sectioned torus reactor in batch and
continuous regimes}, Chem. Eng. J., \textbf{137} (2008), no.~2, 386--395.
\bibitem{Ramirez2017} H.~Ram\'irez, A.~Rojas, D.~Jeison;
\emph{Productivity optimization of microalgae cultivation in a batch
photobioreactor process}, Math. Methods
Appl. Sci., (2017), 1--21.
\bibitem{angel} A. M. Ramos;
\emph{Introducci\'on al an\'alisis matem\'atico del m\'etodo de
elementos finitos}, Editorial Complutense, 2012.
\bibitem{redondo3} J. L. Redondo, J.~Fern{\'a}ndez, I.~Garc{\'\i}a, P.M. Ortigosa;
\emph{Parallel algorithms for continuous competitive location problems},
Optim. Method Soft., \textbf{23} (2008), no.~5, 779--791.
\bibitem{redondo1} J. L. Redondo, J.~Fern{\'a}ndez, I.~Garc{\'\i}a, P.M. Ortigosa;
\emph{A robust and efficient global optimization algorithm for planar
competitive location problems}, Ann. Oper. Res., \textbf{167} (2009), no.~1,
87--105.
\bibitem{redondo2} J. L. Redondo, J.~Fern{\'a}ndez, I.~Garc{\'\i}a, P.M. Ortigosa;
\emph{Solving the multiple competitive location and design problem on
the plane}, Evol. Comput., \textbf{17} (2009), no.~1, 21--53.
\bibitem{Reiter1999} B. R. Reiter;
\emph{Assessment of laboratory methods for quantifying aqueous
bacterial diffusion}, Master's thesis, University of Saskatchewan, 1999.
\bibitem{Richmond1993} A.~Richmond, S.~Boussiba, A.~Vonshak, R.~Kopel;
\emph{A new tubular reactor for mass production of microalgae outdoors},
J. Appl. Phycol., \textbf{5} (1993), no.~3, 327--332.
\bibitem{Rodriguez2014} J. C. Rodr\'iguez, H.~Ram\'irez, P.~Gajardo, A.~Rapaport;
\emph{Optimality of affine control system of several species in competition
on a sequential batch reactor},
Int. J. Control, \textbf{87} (2014), no.~9, 1877--1885.
\bibitem{sharapov2007} R. R. Sharapov;
\emph{Genetic algorithms: Basic ideas, variants and analysis},
InTech, 2007.
\bibitem{Sierra2008} E.~Sierra, F. G. Aci\'en, J.M. Fern\'andez,
J. L. Garc\'ia, C.~Gonz\'alez, E.~Molina;
\emph{Characterization of a flat plate photobioreactor for the
production of microalgae}, Chem. Eng. J., \textbf{138} (2008), no.~1--3,
136--147.
\bibitem{singh2013} V.~Singh;
\emph{Environmental hydrology}, Water Science and Technology Library,
Springer Netherlands, 2013.
\bibitem{SmithWaltman} H.~Smith, P.~Waltman;
\emph{The theory of the chemostat. in cambridge
studies in mathematical biology}, vol.~13, Cambridge: Cambridge University
Press, 1995.
\bibitem{Stewart2003} P. S. Stewart;
\emph{Diffusion in biofilms}, J. Bacteriol., \textbf{185} (2003),
no.~5, 1485--1491.
\bibitem{symon1971} K. R. Symon;
\emph{Mechanics}, Addison-Wesley World student series,
Addison-Wesley Publishing Company, 1971.
\bibitem{Tredici1998} M. R. Tredici, G. C. Zittelli;
\emph{Efficiency of sunlight utilization:
Tubular versus flat photobioreactors}, Biotechnol. Bioeng., \textbf{57}
(1998), no.~2, 187--197.
\bibitem{Venterea2000} R. T. Venterea, D. E. Rolston;
\emph{Mechanistic modeling of nitrite
accumulation and nitrogen oxide gas emissions during nitrification}, J.
Environ. Qual., \textbf{29} (2000), no.~6, 1741--1751.
\bibitem{weise2009} T.~Weise;
\emph{Global optimization algorithms - theory and application},
{S}econd ed., Self-Published, 2009.
\bibitem{Xia2008} J.-Y. Xia, S.-J. Wang, S.-L. Zhang, J.-J. Zhong;
\emph{Computational investigation of fluid dynamics in a recently developed
centrifugal impeller bioreactor}, Biochemical Engineering Journal,
\textbf{38} (2008), no.~3,
406--413.
\bibitem{Zuo2013} W.~Zuo;
\emph{An object-oriented graphics interface design and optimization
software for cross-sectional shape of automobile body}, Adv. Eng. Softw.,
\textbf{64} (2013), 1--10.
\end{thebibliography}
\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}
% end of file template.tex