Model Based Optimal Control of the Photosynthetic Growth of Microalgae in a Batch Photobioreactor

: The present paper investigates the optimal control of the photosynthetic growth process in an artiﬁcial light photobioreactor operated in batch mode, the objective being to ﬁnd an optimal incident light intensity for which the consumption of light energy, for any amount of newly formed biomass, is minimal. By using a simple and reliable model for the photosynthetic growth of microalgae of microalgae, predictions can be made on the quantity of produced biomass and on the amount of light consumed, whose ratio gives the biomass yield on light energy. This variable is unimodal on the allowed range of incident light intensities and has been used as objective function. An improved objective function is proposed by using the speciﬁc growth rate and a weighing factor that allows obtaining the desired amount of biomass while the light energy consumption is optimal. A closed-loop control structure has been designed based on the developed optimization algorithm. The optimal controller has been validated in simulation, comparing different lengths of the optimization horizon and the sampling period. It was found that a bigger sampling period, for the cases where there is no online information on the biomass concentration, does not signiﬁcantly affect the productivity. The optimization algorithm can be used either online or ofﬂine, being useful for various experimental setups.


Introduction
Microalgae are microscopic organisms found in fresh or sea water, but also in nonaquatic habitats.The photosynthetic organisms, eukaryotic or prokaryotic, express varied coloration (e.g., green, blue, brown, red) due to the wide range of intracellular pigments [1,2].Microalgae biotechnology is continuously growing based on the capacity of these organisms to be used in industrial applications, apart from being used for human and animal nutrition [3,4].Added value compounds can be produced by microalgae, such as pigments [2,5], long chain polysaturated fatty acids [6,7], proteins [3,8,9], lipids [10,11], carbohydrates [12,13], etc.Another important role of microalgae is their involvement in environmental applications such as wastewater treatment [14].They can be also involved in the production of different types of biofuels (e.g., biodiesel, biohydrogen, bioethanol, or methane), within the context of an integrated biorefinery, to achieve economic viability [15][16][17].Furthermore, the microalgae's ability to bio-mitigate carbon dioxide [18] is an additional reason to count microalgae biotechnology as very promising from a scientific, technological, and economic point of view.
Most species of microalgae are photoautotrophic, which means that they are using light (solar or artificial) as a source of energy and some inorganic chemicals such as CO 2 , salts of Energies 2022, 15, 6535 2 of 15 nitrogen, salts of phosphorous, etc. Water photolysis is one of the most energy-demanding reactions in nature.The water is split, the hydrogen is used in the photosynthesis transport chain and the O 2 is released as residue.All of these reactions are named light-dependent reactions.CO 2 is required to form primary metabolites (e.g., carbohydrates, lipids, etc.), reactions that do not require light but require the protons and electrons produced by the photodissociation of water.Due to the simple metabolic requirements needed for the growth of microalgae, the photosynthetic metabolism is very attractive.
Microalgae are cultivated in photobioreactors (PBRs) which have at least one transparent side for the light to enter in contact with the culture.They have various shapes (e.g., flat-pate, tubular, columns, etc.), are open or closed, and serve laboratory or even industrial purposes [19,20].PBRs can work in discontinuous (batch), semi-continuous (fed-batch), or continuous modes, each having advantages and disadvantages.This paper focuses on the modeling and optimization of batch processes, these being often used in laboratory studies and even industrial applications for developing and validating mathematical models and control strategies.
Regardless of their low metabolic requirements, their robustness, and the versatile reactors in which they are cultivated, there are several obstacles to overcome for making the microalgal biotechnology economically feasible.In its support, reliable mathematical models validated on experimental data and control strategies must be developed [21][22][23].Targeting high concentrations of microalgal biomass (or high productivities) often requires optimal control.Various studies have been made on the optimal control of continuous photobioreactors, naturally [24][25][26] or artificially lighted [27], where the control variable was the feeding flow rate (in continuous and semi-continuous operating modes).Other key factors of the microalgae cultures, such as the pH control, have also been the subject of optimization techniques; in this case, the input variable being the CO 2 flow rate [28].This paper focuses on the batch mode operation of artificially illuminated PBRs, the control variable being the incident light intensity on the PBR surface.
The light available inside the culture of microalgae is one of the main factors that constrain the growth process because it creates a heterogenous light field.Light is attenuated along the depth of the culture, resulting in decreasing growth rates.The models that can describe the light attenuation inside the microalgae cultures are named radiative models and are fundamental in photobioreactor studies [23,29,30].The radiative models must be coupled with the biological process through kinetic growth models which result in models that express the specific growth rate as a function of irradiance (the light available inside the microalgae culture) [23,31,32].Often, the photobioreactors are operated discontinuously due to certain advantages such as obtaining a high concentration of biomass, lower risk of contamination, etc.They are artificially lighted with efficient light sources that are also dimmable.LED panels are more popular due to their good efficiency and robustness.There are also commercial light growth panels specifically developed for photosynthesis, with powerful COB LEDs (chip on board).By controlling the incident light intensity, a new category of controllers named lumostats can be developed.Lumostatic operation of the PBRs involve the control of variables such as light-to-microalgae ratio, biomass yield on light energy, etc. [32][33][34][35].The productivity of a lumostatic batch can be better when compared with batch cultures operated at constant incident light intensity.Nevertheless, the lumostatic operation of an artificially lighted PBR requires optimization to minimize the amount of light consumed [34][35][36].
The present paper aims to find the optimal set of commands for which the consumption of light energy, for any amount of newly formed biomass, is minimal.A closed-loop optimal controller has been proposed and validated in simulation.The paper has four sections, and after a brief introduction, the mathematical model used for prediction and optimization is presented, followed by the formulation of the optimization problem.The last section presents the optimal controller along with its validation in various scenarios.

The Radiative Model
The photoautotrophic microalgae rely on light as a source of energy and can grow on inorganic nutrients such as dissolved carbon dioxide and salts of nitrogen, phosphorous, magnesium, etc.The mathematical modeling of their growth is, thus, closely related to the availability of light inside the culture.
The light required for photosynthesis has a spectral region that corresponds to the visible light (400 to 700 nm), and it is named Photosynthetically Active Radiation (PAR).Photons with other wavelengths either have not enough or they have too much energy to undertake photosynthesis and can damage the cell.Since the photons are the ones that interact with the cell photosystems, the incident PAR, I 0 , is measured in µmole photons•m −2 •s −1 (to simplify µmol•m −2 •s −1 ) and hereafter named incident light intensity.The incident light intensity that touches the culture of microalgae at z = 0 attenuates inside the PBR due to absorption phenomena and it is named irradiance and is denoted by I. z is the depth of the culture, z ∈ [0, L] and L is the depth of the PBR.The propagation of the light inside the PBR depends on a series of factors such as medium properties, reactor geometry, the shape of cells, the concentration of pigments [21], etc.The attenuation of light inside of a microalgae culture (i.e., the irradiance) is generally described by the Lambert-Beer law, which relates the irradiance to the depth of the culture, z, and the properties of the liquid environment: where the constant ε is an extinction coefficient specific to the aquatic environment.The measuring unit of ε is m −1 so that the exponential factor is dimensionless.The mathematical modeling of the irradiance has been approached by various authors [21,29] aiming to describe more accurately the behavior of light in very dense media, such as microalgae cultures, where it is absorbed by the pigments but also scattered by the cells.In addition, the geometry of the PBR and the design of the lighting source are heavily influencing the attenuation of the light inside the microalgae cultures.These types of models are named radiative models.
The model adopted in this study relates the attenuation of light with the concentration of biomass, X, as follows [29]: where α = E a /(E a + 2bE s ) is the linear scattering modulus.For the exponential to be dimensionless, the measuring unit of the coefficient ε must be m 2 •kg −1 (same as E a and E s as reported in Table 1).The extinction coefficient, ε, is calculated by considering the mass absorption and mass scattering coefficients (i.e., E a and E s ) and the backward scattering fraction, b. Figure 1 shows the attenuation of the light with the depth of the culture for various concentrations of biomass.To better illustrate the attenuation of light, the normalized irradiance, I(z)/I 0 has been calculated.
The radiative model, in Equation (2), has been developed for flat-plate PBRs, and lighted on one side.To illustrate the normalized irradiance, in Figure 1, a depth of 0.054 m has been considered.This corresponds to an airlift rectangular PBR, used in previous studies to obtain experimental data and to identify the mathematical models [37].Figure 2 illustrates the airlift PBR along with its main components.The values of the parameters in the radiative model are given subsequently in Table 1.The radiative model, in Equation ( 2), has been developed for flat-plate PBRs, and lighted on one side.To illustrate the normalized irradiance, in Figure 1, a depth of 0.054 m has been considered.This corresponds to an airlift rectangular PBR, used in previous studies to obtain experimental data and to identify the mathematical models [37].Figure 2 illustrates the airlift PBR along with its main components.The values of the parameters in the radiative model are given subsequently in Table 1.

The Kinetic Growth Model
The heterogeneous field created by the light inside the PBR results in growth rates that are decreasing along the depth of the culture.The coupling between the radiative models and the growth kinetics is an active research topic in PBR studies [23,31,32].Because light energy is the most important factor for photosynthetic growth, it is considered a substrate.The specific growth rate related to photosynthesis,  , must thus be expressed as a function of the irradiance.However, the energy absorbed by microalgae is not all used for growth, a part of it being used for maintenance and, thus, a negative term is defined as the specific decay rate,  .The specific growth rate, , has, therefore, two terms, one associated with growth and one to decay:  The radiative model, in Equation ( 2), has been developed for flat-plate PBRs, and lighted on one side.To illustrate the normalized irradiance, in Figure 1, a depth of 0.054 m has been considered.This corresponds to an airlift rectangular PBR, used in previous studies to obtain experimental data and to identify the mathematical models [37].Figure 2 illustrates the airlift PBR along with its main components.The values of the parameters in the radiative model are given subsequently in Table 1.

Figure 2. (A)
. Front view of an airlift PBR with a LED growth light panel and their main components [37].(B).The lateral view of the PRB showing the attenuation of the light inside the microalgae culture.

The Kinetic Growth Model
The heterogeneous field created by the light inside the PBR results in growth rates that are decreasing along the depth of the culture.The coupling between the radiative models and the growth kinetics is an active research topic in PBR studies [23,31,32].Because light energy is the most important factor for photosynthetic growth, it is considered a substrate.The specific growth rate related to photosynthesis,  , must thus be expressed as a function of the irradiance.However, the energy absorbed by microalgae is not all used for growth, a part of it being used for maintenance and, thus, a negative term is defined as the specific decay rate,  .The specific growth rate, , has, therefore, two terms, one associated with growth and one to decay:

The Kinetic Growth Model
The heterogeneous field created by the light inside the PBR results in growth rates that are decreasing along the depth of the culture.The coupling between the radiative models and the growth kinetics is an active research topic in PBR studies [23,31,32].Because light energy is the most important factor for photosynthetic growth, it is considered a substrate.The specific growth rate related to photosynthesis, µ p , must thus be expressed as a function of the irradiance.However, the energy absorbed by microalgae is not all used for growth, a part of it being used for maintenance and, thus, a negative term is defined as the specific decay rate, µ d .The specific growth rate, µ, has, therefore, two terms, one associated with growth and one to decay: The PBR is considered a continuous stirred tank reactor (CSTR), which means that the concentration of any component is uniform throughout the reactor volume and can be expressed through ordinary differential equations.Nevertheless, the heterogenous light field brings a second derivation variable, the depth z of the culture, which requires partial differential equations.To avoid a complex mathematical formulation, an averaged growth rate could be defined in two manners: the averaged specific growth rate, µ I avg , can be calculated with an averaged irradiance, I avg .An analytical expression for the average of the irradiance (Equation ( 2)) along the depth z of the culture can be easily defined [30,38]; -local photosynthetic responses, µ(I(z)), can be calculated for any depth z of the culture.These local photosynthetic responses are simply averaged into an average photosynthetic response, µ (or average specific growth rate) [31].The denote an averaged value.
The present study is based on the second approach, an averaged photosynthetic response being calculated from 100 local photosynthetic responses.The specific growth rate related to photosynthesis is expressed through a classical inhibition model: where µ o is a kinetic parameter related to the maximum specific growth rate, K S is the saturation constant and K I is the inhibition constant.L is the depth of the PBR (z The specific decay rate, µ d , is considerably smaller than µ o and, even though it varies with illumination conditions, it is considered constant here. Figure 3 shows the specific growth rate decrease along the depth z of the culture for a series of biomass concentrations.The averaged specific growth rates have also been displayed in the legend.When the culture is dense (e.g., 2.0 g•L −1 ), the light penetrates less than a quarter of the culture, leaving the rest of it in dark.In the dark zone of the culture the specific growth rate is negative, µ p < µ d .The average specific growth rate is very low (i.e., 0.002 h −1 ) for this concentration of biomass, and an incident light intensity lower than 500 µmol•m −2 •s −1 leads to the decrease of the biomass concentration.The average of the specific growth rates is justifiable because, in a CSTR, a cell is migrating continuously between the light and the dark zones.

𝜇 𝜇 𝜇
The PBR is considered a continuous stirred tank reactor (CSTR), which means that the concentration of any component is uniform throughout the reactor volume and can be expressed through ordinary differential equations.Nevertheless, the heterogenous light field brings a second derivation variable, the depth  of the culture, which requires partial differential equations.To avoid a complex mathematical formulation, an averaged growth rate could be defined in two manners: the averaged specific growth rate,   , can be calculated with an averaged irradiance,  .An analytical expression for the average of the irradiance (Equation ( 2)) along the depth  of the culture can be easily defined [30,38]; -local photosynthetic responses,    , can be calculated for any depth  of the culture.These local photosynthetic responses are simply averaged into an average photosynthetic response, 〈〉 (or average specific growth rate) [31].The 〈 〉 denote an averaged value.
The present study is based on the second approach, an averaged photosynthetic response being calculated from 100 local photosynthetic responses.The specific growth rate related to photosynthesis is expressed through a classical inhibition model: where  is a kinetic parameter related to the maximum specific growth rate,  is the saturation constant and  is the inhibition constant. is the depth of the PBR ( ∈ 0  ).The specific decay rate,  , is considerably smaller than  and, even though it varies with illumination conditions, it is considered constant here.
Figure 3 shows the specific growth rate decrease along the depth  of the culture for a series of biomass concentrations.The averaged specific growth rates have also been displayed in the legend.When the culture is dense (e.g., 2.0 g•L −1 ), the light penetrates less than a quarter of the culture, leaving the rest of it in dark.In the dark zone of the culture the specific growth rate is negative,   .The average specific growth rate is very low (i.e., 0.002 h −1 ) for this concentration of biomass, and an incident light intensity lower than 500 μmol•m −2 •s −1 leads to the decrease of the biomass concentration.The average of the specific growth rates is justifiable because, in a CSTR, a cell is migrating continuously between the light and the dark zones.To express the biomass as concentration, the volumetric growth rate, r x , is defined:

The Mass Balance Model
Integrating the volumetric growth rate, the concentration of biomass, X, in g•L −1 can be obtained: Energies 2022, 15, 6535 6 of 15 The differential equation for biomass, Equation ( 6), is given for a discontinuous process (batch).A batch process cannot be controlled because the reactor is filled with culture medium at the beginning, and the process takes place for a specific amount of time without any intervention.The artificially lighted PBRs, in contrast with the solar PBRs, have the incident light intensity, I 0 , as input variable and can be controlled to obtain optimal results.These reactors are used in laboratories for data acquisition, modeling, developing and validating control strategies, production of metabolites, etc.In most cases, a batch process in an artificially lighted PBR takes place at a constant incident light intensity, but many light growth panels are now dimmable and can be controlled to obtain optimal productivities or optimal consumption of light energy.
The specific growth rate can be a function of other variables besides the light, such as substrates, pH, temperature, etc., and more equations can be added to the model.However, in this study, we consider that the substrates are in concentrations that are not limiting or inhibiting the growth and the other key factors are controlled.The biomass equation can, thus, be decoupled from other equations of substrates, for example, and be manipulated individually.
By using the incident light intensity as the control variable, the PBRs can be operated in lumostatic mode, which is controlling the I 0 for maintaining constant certain variables such as light-to-microalgae ratio [32,34], etc.The batch cultures with controlled incident light intensities are termed lumostatic batches.The efficiency of a lumostatic batch can be evaluated through the biomass yield on light energy, Y, a performance index that has the following expression [32]: where X 0 is the concentration of biomass, in g•L −1 , at the beginning of the batch (t = 0), and V is the working volume of the PBR.The total light consumed in the batch process, Q, is calculated by integrating the incident light intensity over the batch time t 0 t f : where A is the lighted surface of the PBR.Here, I 0 (t) is the incident light intensity that is variable over the control horizon, while I 0 is used for a constant incident light intensity.Y has, thus, the meaning of grams of biomass produced per moles of photons consumed (a conversion constant equal to 3600 × 10 −6 is needed to express Q in moles, where 3600 converts seconds to hours and 10 −6 converts micromoles to moles, i.e., Q = A•q•3600 × 10 −6 ).
The mathematical model used in this study consists of four ordinary differential equations (Equation ( 9)) and a series of algebraic equations (Equations ( 2)-( 5) and ( 7)): where x 1 is the concentration of biomass in g•L −1 , x 2 is the newly produced biomass in grams, x 3 is the amount of light consumed in moles of photons and x 4 is the biomass yield on light energy in grams per moles of photons.The fourth equation is in fact algebraic but is expressed as a differential equation.The integration of the mass balance model must be done with an integrator able to solve DAEs (differential-algebraic equations) such as ode15s or ode23t in Matlab [39].
The mathematical model for photosynthetic biomass growth has been validated on Chlamydomonas reinhardtii [23] and Scenedesmus quadricauda [37], strains that have been cultivated in rectangular laboratory PBRs lighted on one side.The parameters of the model used in this study are presented in Table 1.

Parameter
Value Unit

The Optimal Control Problem
The objective of this paper is to find an optimal incident light intensity for which the consumption of light energy, for any amount of newly formed biomass, is minimal.To this end, reliable measurements for biomass must be available.There are direct (i.e., microscopic cell count) and indirect (e.g., dry mass, optical density, etc.) measurements of biomass concentration.While the direct measurement of biomass concentration requires microbiology skills, the dry mass is the most used measurement due to its simplicity, which requires filtering known volumes of the culture and weighing the filters before and after drying at 105 • C. The biggest disadvantage of this indirect measurement is the long time required for drying, in most cases only one value per day being available.The need for online biomass concentration measuring can be satisfied by hardware (e.g., turbidity sensors, submersible spectrophotometers, etc.) or software sensors.However, both offline and online measurements must be correlated with the dry mass.
The limited availability of online biomass measuring solutions is because the signals of the above-mentioned sensors are influenced by many factors, such as the geometry of the reactor, the composition of the culture medium, the shape of the cells, etc.A turbidity sensor, for example, is an optical sensor able to measure the suspended solids, but its signal is disturbed by the gas bubbles in the culture.Provided that it can be installed in a zone where its signal is not disturbed, the turbidity is linearly correlated with the dry mass [22], requiring the identification of two parameters, the slope, and the offset.Excellent correlation can be obtained during the exponential growth phase, while the lag, stationary, or death phases require a different set of parameters for the correlation.In a process, each new dry mass value can be used to adjust the correlation, resulting in better estimations of biomass concentration.
When no hardware sensor for biomass is available in the process, which is the case for most bioprocesses, a software sensor can be used [36].The software sensors can be based on state estimation, regression analysis, artificial neural networks, statistical machine learning [40], etc.However, the process mechanism can be also used to estimate the biomass concentration.In this spirit, the mathematical model presented in the previous section can be used.Two of its parameters could be reidentified daily on the dry mass measurements to obtain a very good fitting, namely µ 0 and K S (Equation ( 4)).To have an estimation, the mathematical model requires the initial conditions and the current process time.
We assume from here on that a value of the biomass concentration is available.Knowing the biomass concentration and the volume of the reactor, V, the newly produced biomass, X T , can be easily calculated (Equation ( 9)).To achieve the goal of finding the optimal incident light intensity for which the biomass yield on light energy is maximal, Energies 2022, 15, 6535 8 of 15 the light source must be dimmable and controlled by the process computer.The light sources can be calibrated with light meters, correlating the input voltage with the incident light intensity in µmol•m −2 •s −1 .To calculate the amount of light consumed, Q, the model requires the incident light intensity, I 0 , the lighted surface of the PBR, and the current process time.By dividing the newly produced biomass by the amount of light consumed, the biomass yield on light energy can be easily calculated (Y = X T /Q).
Integrating the mathematical model over a wide range of incident light intensities, 1 to 2000 µmol•m −2 •s −1 , it can be found that the biomass yield on light energy, Y, is a unimodal unidimensional function whose maximum gives the optimal incident light intensity, I 0(opt) (Figure 4).The optimization of Y can be done with unidimensional optimization methods [41] such as the golden section search, the parabolic interpolation, etc.The constraints to which the input variable I 0 is subjected have been selected from a technological point of view, 1 µmol•m −2 •s −1 to avoid division by zero (Equation ( 7)), and 2000 µmol•m −2 •s −1 is the maximal sun radiation on a summer day.cess time.
We assume from here on that a value of the biomass concentration is available.Knowing the biomass concentration and the volume of the reactor, , the newly produced biomass,  , can be easily calculated (Equation ( 9)).To achieve the goal of finding the optimal incident light intensity for which the biomass yield on light energy is maximal, the light source must be dimmable and controlled by the process computer.The light sources can be calibrated with light meters, correlating the input voltage with the incident light intensity in μmol•m −2 •s −1 .To calculate the amount of light consumed, , the model requires the incident light intensity,  , the lighted surface of the PBR, and the current process time.By dividing the newly produced biomass by the amount of light consumed, the biomass yield on light energy can be easily calculated (   ⁄ ).
Integrating the mathematical model over a wide range of incident light intensities, 1 to 2000 μmol•m −2 •s −1 , it can be found that the biomass yield on light energy, , is a unimodal unidimensional function whose maximum gives the optimal incident light intensity,  (Figure 4).The optimization of  can be done with unidimensional optimization methods [41] such as the golden section search, the parabolic interpolation, etc.The constraints to which the input variable  is subjected have been selected from a technological point of view, 1 μmol•m −2 •s −1 to avoid division by zero (Equation ( 7)), and 2000 μmol•m −2 •s −1 is the maximal sun radiation on a summer day.
Figure 4 displays the response of the system at various incident light intensities which have been obtained on a one-hour horizon.The simulations have been made on a horizon of one hour, mentioning that simulations on wider horizons are similar in terms of  , as can be seen in Figure 5.The initial  The simulations have been made on a horizon of one hour, mentioning that simulations on wider horizons are similar in terms of I 0(opt) , as can be seen in Figure 5.The initial conditions used for simulations are X(0) = 1.1386 g•L −1 , X T (0) = 4.4690 g (not displayed in Figure 4), Q(0) = 9.18 g•mol photon −1 , which means that Y(0) = X T (0)/Q(0) = 0.4868.
The mathematical model for the photosynthetic growth of microalgae can give even more information on the process.The specific growth rate, µ, provides valuable bounds for the control variable I 0 : a lower bound, I 0(min) , under which the biomass would decrease.Below I 0(min) the specific growth rate is negative (Figure 4), and an upper bound, I 0(max) , which is set in Figure 4 at 80% of the maximum growth rate.
The biomass yield on light energy, Y, can thus be used in the optimization process.The optimal control problem starts from the photosynthetic growth model which has been presented in the second section.It consists of a system of four differential equations, therefore the state vector is: Energies 2022, 15, 6535 9 of 15 while  0 1.1386 g•L requires ≈ 331 μmol•m •s .This points out the need for an in-crease in the incident light intensity as the biomass grows.Working with constant light throughout the photosynthetic growth process could require more light energy to obtain the same amount of biomass compared to an increasing light profile.This problem has been addressed in [32] by controlling the light-to-microalgae ratio, but the present paper takes it a step further by investigating the optimization of a lumostatic batch (i.e., a batch with variable incident light intensity).The disadvantage of using the biomass yield on light energy, , as objective function is that it lacks a scaling factor to push the concentration of biomass higher. assures that the consumption of light energy is optimal, but this could result in low concentrations of biomass.In this regard, the information of the specific growth rate, , which is already used to determine the lower bound  , can be used as follows: where  is a positive scaling factor.  has the effect of shifting the optimum to the right and, along with the scaling factor, any increasing   profile can be obtained.Thus, higher concentrations of biomass can be obtained, while the consumed amount of light energy is minimal.
The desired increasing light profile can be obtained by iteratively optimizing the process for any given current time and state.In this regard, a closed-loop control system has been proposed in the following section.

The Control Structure
The objective of this work is to provide a closed-loop solution for the optimal control problem presented in the previous section; in other words, an optimal controller that provides an incident light intensity value to the process, for each sampling period.Figure 6 presents the closed-loop control structure proposed in this work.The optimal controller calculates the optimal incident light intensity,  , over the optimization horizon   , with  0, … ,  where  and  are the beginning and the end of the photosynthetic growth process. is a positive constant, expressed in hours, which is the length The input variable for the system is the incident light intensity: The optimal control problem is subjected to the following constraints: the control horizon: t 0 ≤ t ≤ t f , where t 0 and t f are the initial and the final time of the control horizon (the batch period), -the initial conditions: x(t 0 ) = x 0 , -the set of lower and upper bounds: I 0(min) ≤ I 0 (t) ≤ I 0(max) , with t 0 ≤ t ≤ t f .The lower bound, I 0(min) , is critical because the biomass decreases under this value, while the upper bound is not, and can remain constant (e.g., 2000 µmol•m −2 •s −1 ).Even though I 0(max) can inhibit the microalgae growth, it is attenuated inside the culture.The lower and the upper bounds are required by most of the unidimensional optimization methods [41], e.g., fminbnd function in Matlab.
The general form of the optimal control problem which searches the control u(•) that maximizes the objective function J is: where the term L(x, u, t) is the intermediate cost function.The optimal control aims to find the state and control trajectories x and u so that J(x, u) is maximized: Finding the optimal incident light intensity for which the consumption of light energy, for any amount of newly formed biomass, is minimal, requires maximizing the biomass yield on light energy, thus the optimal control input has the following form: The optimization and control horizons are usually different, while the optimization can be performed on the entire batch time t 0 t f , the control horizon is usually lower, resulting that the u(t) is the vector of input variables.If t k is the current time, t 0 ≤ t k ≤ t f , the control horizon can be t k t k+p .Figure 5  As has been the case previously (Figure 4), the results are the final values of the states, after the integration of the model for one hour with each Optimization on three horizons is presented, 1, 84, and 168 h, taking a seven-day cultivation process as a reference.The initial condition for biomass is X(0) = 0.3 g•L −1 , while the other states have been set to zero (Q(0) should be set to a very small positive value to avoid division by zero).It can be observed that the optimal incident light intensity is similar, regardless of the length of the optimization horizon.Comparing Figure 4 with 5, it can be observed that higher concentrations of biomass require higher concentrations of light to reach an optimal yield, X(0 This points out the need for an increase in the incident light intensity as the biomass grows.Working with constant light throughout the photosynthetic growth process could require more light energy to obtain the same amount of biomass compared to an increasing light profile.This problem has been addressed in [32] by controlling the light-to-microalgae ratio, but the present paper takes it a step further by investigating the optimization of a lumostatic batch (i.e., a batch with variable incident light intensity).The disadvantage of using the biomass yield on light energy, Y, as objective function is that it lacks a scaling factor to push the concentration of biomass higher.Y assures that the consumption of light energy is optimal, but this could result in low concentrations of biomass.In this regard, the information of the specific growth rate, µ, which is already used to determine the lower bound I 0(min) , can be used as follows: where w is a positive scaling factor.µ(t) has the effect of shifting the optimum to the right and, along with the scaling factor, any increasing I 0 (t) profile can be obtained.Thus, higher concentrations of biomass can be obtained, while the consumed amount of light energy is minimal.
The desired increasing light profile can be obtained by iteratively optimizing the process for any given current time and state.In this regard, a closed-loop control system has been proposed in the following section.

The Control Structure
The objective of this work is to provide a closed-loop solution for the optimal control problem presented in the previous section; in other words, an optimal controller that provides an incident light intensity value to the process, for each sampling period.Figure 6 presents the closed-loop control structure proposed in this work.The optimal controller calculates the optimal incident light intensity, I 0(opt) , over the optimization horizon [t k t k+n ], with k = 0, . . ., f where t 0 and t f are the beginning and the end of the photosynthetic growth process.n is a positive constant, expressed in hours, which is the length of the optimization horizon.To calculate the current control variable, I k 0(opt) , the optimal controller requires the current state of the process, x k , and current time, t k .Thus, the optimal controller must predict the future evolution of the process, and calculate the optimal incident light intensity based on these predictions.This requires a process model to be included in the optimal controller.The lower bound, I 0(min) , can be found through a simple numerical root finding method (e.g., the bisection method) because µ intersects the Ox axis only once.There is no specific method to find the upper bound, I 0(max) , in Figure 4  incident light intensity based on these predictions.This requires a process model to be included in the optimal controller.The lower bound,  , can be found through a simple numerical root finding method (e.g., the bisection method) because  intersects the Ox axis only once.There is no specific method to find the upper bound,  , in Figure 4 being set at 80% of the specific growth rate obtained for 2000 μmol•m −2 •s −1 .The current lower and upper bounds,    , are calculated with the process model before the optimization procedure.The sampling period of the photosynthetic growth process is , where  is a positive constant, expressed in hours, which is the length of the sampling period (usually  ).The process generates a new set of values for the state variables,  , which are used by the optimal controller to calculate a new input, along with the current time and current bounds.

Simulation Results
A seven-day (168 h) photosynthetic batch process has been used as reference to validate the optimal controller in simulation.The initial conditions used for all simulations in this section are  0 0.3 g•L −1 , while the other three states have been set to zero (same as for Figure 5).
Figure 7 displays the results obtained with the closed loop optimal controller on a 1h sampling period,  1,   .The behavior observed in Figure 5 is also observed in closed loop, and the length of the optimization horizon is not significantly influencing the process.As the optimization horizon is bigger, slightly lower biomass values are obtained.For example, the 20-h optimization horizon results in a biomass concentration of only 0.96% lower in comparison with the 1-h optimization horizon.Also, lower  values give similar  , because  decreases throughout the batch culture.A higher optimization horizon increases computational complexity of the algorithm and thus a 1-h horizon is recommended.
In addition to the optimization horizon, the sampling period must be determined.Given that the biological process is slow, a 1-h period has been considered.This period can be shortened by the capacity of the process computer to measure the states and calculate a new control variable.Nevertheless, longer sampling periods have been investigated for cases where an online biomass sensor is not available, and the optimization procedure cannot be performed online.This would require the daily reidentification of the mathematical model based on dry mass measurements and the optimization of the process based on the model predictions, the operator being able to insert a new command every 12 h.The sampling period of the photosynthetic growth process is t k t k+p , where p is a positive constant, expressed in hours, which is the length of the sampling period (usually n ≥ p).The process generates a new set of values for the state variables, x k+1 , which are used by the optimal controller to calculate a new input, along with the current time and current bounds.

Simulation Results
A seven-day (168 h) photosynthetic batch process has been used as reference to validate the optimal controller in simulation.The initial conditions used for all simulations in this section are X(0) = 0.3 g•L −1 , while the other three states have been set to zero (same as for Figure 5).
Figure 7 displays the results obtained with the closed loop optimal controller on a 1-h sampling period, p = 1, [t k t k+1 ].The behavior observed in Figure 5 is also observed in closed loop, and the length of the optimization horizon is not significantly influencing the process.As the optimization horizon is bigger, slightly lower biomass values are obtained.For example, the 20-h optimization horizon results in a biomass concentration of only 0.96% lower in comparison with the 1-h optimization horizon.Also, lower Y values give similar I 0(opt) , because Y decreases throughout the batch culture.Figure 8 shows the evolution of the process with the closed loop optimal controller on 1-h and 12-h sampling periods.Similar results are obtained in terms of biomass growth, which proves that longer sampling periods are feasible for systems that have certain limitations.A higher optimization horizon increases computational complexity of the algorithm and thus a 1-h horizon is recommended.
In addition to the optimization horizon, the sampling period must be determined.Given that the biological process is slow, a 1-h period has been considered.This period can be shortened by the capacity of the process computer to measure the states and calculate a new control variable.Nevertheless, longer sampling periods have been investigated for cases where an online biomass sensor is not available, and the optimization procedure cannot be performed online.This would require the daily reidentification of the mathematical model based on dry mass measurements and the optimization of the process based on the model predictions, the operator being able to insert a new command every 12 h.
Figure 8 shows the evolution of the process with the closed loop optimal controller on 1h and 12-h sampling periods.Similar results are obtained in terms of biomass growth, which proves that longer sampling periods are feasible for systems that have certain limitations.Figure 8 shows the evolution of the process with the closed loop optimal controller on 1-h and 12-h sampling periods.Similar results are obtained in terms of biomass growth, which proves that longer sampling periods are feasible for systems that have certain limitations.
8. Simulation results on a 7-day batch with different sampling periods (1-h sampling period-blue and 12-h sampling period-red).
The results presented in Figures 7 and 8 are obtained using the biomass yield on light energy as objective function.Regardless of the optimization horizon and the sampling period, the concentration of biomass is slightly lower than 1.5 g•L −1 .Further increasing the sampling period results in significantly lower concentrations of biomass, as can be seen in Table 2.A constant incident light intensity of 133.17 μmol•m −2 •s −1 throughout the The results presented in Figures 7 and 8 are obtained using the biomass yield on light energy as objective function.Regardless of the optimization horizon and the sampling period, the concentration of biomass is slightly lower than 1.5 g•L −1 .Further increasing the sampling period results in significantly lower concentrations of biomass, as can be seen in Table 2.A constant incident light intensity of 133.17 µmol•m −2 •s −1 throughout the entire batch (see also Figure 5) results in a biomass concentration 35.48% lower compared with the 1-h sampling period.Higher biomass concentrations can be obtained with higher incident light intensities, and this can be achieved by using the objective function in Equation (15).
Figure 9 shows the simulation results on a seven-day batch for different weighing factors.A higher w results in a higher concentration of biomass, providing that the amount of light energy consumed is minimal.Because the second term of the objective function w * µ(t), is considerably smaller than Y(t), both objective functions have generated similar results (3rd and 4th graphs of Figure 9).The scaling factor thus provides the means to adapt the optimization algorithm to the technological limitations of the systems (i.e., the maximal output of the light source) and allows for the obtaining of a desired amount of biomass.
amount of light energy consumed is minimal.Because the second term of the objective function  *   , is considerably smaller than   , both objective functions have generated similar results (3rd and 4th graphs of Figure 9).The scaling factor thus provides the means to adapt the optimization algorithm to the technological limitations of the systems (i.e., the maximal output of the light source) and allows for the obtaining of a desired amount of biomass.

Conclusions
The paper proposes a closed-loop control system for the process of photosynthetic growth of microalgae.A simple mathematical model with four state variables has been proposed for predicting the evolution of biomass concentration, the total quantity of newly formed biomass, the total amount of consumed light, and the biomass yield on light energy.It has been found that the biomass yield on light energy, , is a unimodal function on the admissible range of incident light intensities,  , which can be used to optimize the lumostatic batch culture.Maximizing  results in finding the optimal incident light intensity for which the consumption of light energy, for any amount of newly formed biomass, is minimal.A second objective function has been proposed by adding a second term

Conclusions
The paper proposes a closed-loop control system for the process of photosynthetic growth of microalgae.A simple mathematical model with four state variables has been proposed for predicting the evolution of biomass concentration, the total quantity of newly formed biomass, the total amount of consumed light, and the biomass yield on light energy.It has been found that the biomass yield on light energy, Y, is a unimodal function on the admissible range of incident light intensities, I 0 , which can be used to optimize the lumostatic batch culture.Maximizing Y results in finding the optimal incident light intensity for which the consumption of light energy, for any amount of newly formed biomass, is minimal.A second objective function has been proposed by adding a second term to Y, which consists of using the specific growth rate along with a weighing factor.This new term shifts the optimum to the right and allows for the obtaining of higher concentrations of biomass.While the first objective function ensures an optimal consumption of light energy, the second objective function also targets a higher amount of newly formed biomass.The developed optimization algorithms have been integrated into a closed-loop structure.The optimal controller is generating an increasing light profile which is more efficient, in terms of consumed energy, compared with a constant light batch.It has been proven that a 1-h optimization horizon is sufficient for obtaining good results and reducing the complexity of the computation.The sampling period has also been investigated, providing a solution for the processes with certain limitations.Finally, it has been shown that tuning the weighing factor has led to the achievement of the desired increasing light profile.The results have been validated in simulation; a reliable solution has been proposed for the optimal control of lumostatic batches.

Figure 1 .
Figure 1.The normalized irradiance attenuation with the depth of the culture, at an incident light intensity of 500 μmol•m −2 •s −1 .

Figure 2 .
Figure 2. (A).Front view of an airlift PBR with a LED growth light panel and their main components [37].(B).The lateral view of the PRB showing the attenuation of the light inside the microalgae culture.

Figure 1 .
Figure 1.The normalized irradiance attenuation with the depth of the culture, at an incident light intensity of 500 µmol•m −2 •s −1 .

Figure 1 .
Figure 1.The normalized irradiance attenuation with the depth of the culture, at an incident light intensity of 500 μmol•m −2 •s −1 .

Figure 2 .
Figure 2. (A).Front view of an airlift PBR with a LED growth light panel and their main components [37].(B).The lateral view of the PRB showing the attenuation of the light inside the microalgae culture.

Figure 3 .
Figure 3.The decrease of the specific growth rate with the depth of the culture, at an incident light intensity of 500 μmol•m −2 •s −1 .

Figure 3 .
Figure 3.The decrease of the specific growth rate with the depth of the culture, at an incident light intensity of 500 µmol•m −2 •s −1 .

Figure 4
Figure 4 displays the response of the system at various incident light intensities which have been obtained on a one-hour horizon.The simulations have been made on a horizon of one hour, mentioning that simulations on wider horizons are similar in terms of I 0(opt) , as can be seen in Figure5.The initial conditions used for simulations are X(0) = 1.1386 g•L −1 , X T (0) = 4.4690 g (not displayed in Figure4), Q(0) = 9.18 g•mol photon −1 , which means that Y(0) = X T (0)/Q(0) = 0.4868.The mathematical model for the photosynthetic growth of microalgae can give even more information on the process.The specific growth rate, µ, provides valuable bounds for the control variable I 0 :
being set at 80% of the specific growth rate obtained for 2000 µmol•m −2 •s −1 .The current lower and upper bounds, I k 0(min) ≤ I k 0(opt) ≤ I k 0(max) , are calculated with the process model before the optimization procedure.

Figure 6 .
Figure 6.The closed-loop control scheme of the photosynthetic growth process.

Figure 6 .
Figure 6.The closed-loop control scheme of the photosynthetic growth process.

Figure 8 .
Figure 8. Simulation results on a 7-day batch with different sampling periods (1-h sampling periodblue and 12-h sampling period-red).

Figure 9 .
Figure 9. Simulation results on a seven-day batch with different weighing factors (see Equation (15)).With blue  1, with red  0.7 and with yellow  0.4.

Figure 9 .
Figure 9. Simulation results on a seven-day batch with different weighing factors (see Equation (15)).With blue w = 1, with red w = 0.7 and with yellow w = 0.4.

Table 2 .
Influence of the sampling period on the concentration of biomass.