Control of Microalgae Growth in Artificially Lighted Photobioreactors Using Metaheuristic-Based Predictions

A metaheuristic algorithm can be a realistic solution when optimal control problems require a significant computational effort. The problem stated in this work concerns the optimal control of microalgae growth in an artificially lighted photobioreactor working in batch mode. The process and the dynamic model are very well known and have been validated in previous papers. The control solution is a closed-loop structure whose controller generates predicted control sequences. An efficient way to make optimal predictions is to use a metaheuristic algorithm, the particle swarm optimization algorithm. Even if this metaheuristic is efficient in treating predictions with a very large prediction horizon, the main objective of this paper is to find a tool to reduce the controller’s computational complexity. We propose a soft sensor that gives information used to reduce the interval where the control input’s values are placed in each sampling period. The sensor is based on measurement of the biomass concentration and numerical integration of the process model. The returned information concerns the specific growth rate of microalgae and the biomass yield on light energy. Algorithms, which can be used in real-time implementation, are proposed for all modules involved in the simulation series. Details concerning the implementation of the closed loop, controller, and soft sensor are presented. The simulation results prove that the soft sensor leads to a significant decrease in computational complexity.


Introduction
Microalgae biotechnology has attracted increased interest in recent years due to its potential to be used in a wide range of applications. Microalgae can be used in human and animal nutrition and also in the production of a wide range of added-value compounds [1,2]. They can also be involved in environmental applications such as the treatment of wastewater or the production of biofuels. In addition, their capacity to biomitigate carbon dioxide makes microalgae cultivation one of the most eco-friendly biotechnologies.
Light energy is required to perform water photolysis, one of the most energydemanding reactions in nature. Photosynthetic growth is very attractive due to its simple and cheap requirements. Microalgae may be cultivated in photobioreactors that have at least one transparent side through which the light radiates the culture [3].
This work is focused on controlling microalgae growth in a specific artificially lighted photobioreactor (PBR). Therefore, dynamic models for photobioreactors and control approaches have been of interest to us. Many references concerning microalgae cultures present reliable mathematical models validated with experimental data and propose control strategies [4][5][6].
Light is the main factor that restricts the photosynthetic growth process. It creates a heterogenous light field inside the photobioreactor. The microalgae growth rate decreases • A specific artificially lighted PBR for microalgae growth has to be controlled over a given control horizon. This batch PBR is described in [10,11], and sufficient details are given in Section 2 and Appendix A.

•
The dynamic process model (PM) concerning the microalgae growth is presented and validated in [10,11]. Section 2 recalls the basic elements that led to this PM used in our work. • A bi-criteria optimal control problem concerning the PBR at hand is stated in Section 3.1.
The optimum criterion requires minimizing the amount of light (the energy) while a specified quantity of biomass is produced.

•
The optimal control problem must be solved through a closed-loop control structure. Only in this way will the controller's output depend on the real state of the PBR (e.g., the biomass concentration) in each sampling period.
The closed-loop control structure proposed in this work is the Receding Horizon Control (RHC) structure [16,17], a general and efficient control structure involving a prediction technique. A module of the Controller, called the Predictor, generates sequences of control output values, called predicted control sequences (pcs), that will determine the process's optimal trajectories. This solution entails two other options:

•
Generating an optimal pcs requires a significant computational effort, as in our control problem. For this reason, the Predictor will resort to a metaheuristic algorithm [18,19] that can cope with high computational complexity. An adaptive version of the hybrid topology particle swarm optimization [20,21] algorithm will be integrated into the Predictor. • Despite the ability of the particle swarm optimization algorithm to efficiently generate optimal predictions, we propose a soft sensor to further decrease the computational complexity, especially for very large prediction horizons.
where α is the linear scattering modulus. 3.4 is devoted to the proposed soft sensor and presents its necessity, structure, functioning, and implementation.
To keep the presentation self-contained, Section 4 briefly describes PSO and its advances, namely, hybrid topology and adaptation of particle speed. To conduct the simulation study presented in Section 5, we implemented the proposed algorithms using the MATLAB language and system. Many algorithms were implemented realistically, like in a real-time application. Section 6 presents and discusses the simulation results. These results prove that the control solutions and algorithms work well and that the soft sensor makes a real contribution.

Photosynthetic Growth of Microalgae: The Process Dynamic Model
The objective of this section is to specify the PBR's dynamic model, which will be used in the following sections to state and solve a specific optimal control problem. The basic aspects of the considered process model and the PBR's parameters have already been presented and validated (see [10,11] and Appendix A).
The microalgae are cultivated in photobioreactors of various geometries with transparent regions allowing light to penetrate the culture. The incident light intensity 0 q is measured at the point where the light touches the microalgae culture. It is attenuated inside the PBR, with the available light for any depth, , of the culture being named irradiance ( ). Models that describe the attenuation of light inside microalgal cultures are named radiative models. This work considers a continuously stirred flat-plate PBR lighted on one side, already presented in [10,11]. The radiative model considered here was presented in [7]. The irradiance is related to the biomass concentration (denoted as X for the moment), the depth of the culture (z, see Figure 1), and few coefficients: the mass absorption ( a E ), the mass scattering ( S E ), and the backward scattering fraction (b).

( )
where α is the linear scattering modulus. The attenuation of light creates a heterogenous light field; the growth decreases along with the increase in the culture's depth. Coupling between radiative models and growth kinetics models is present in many papers and has been investigated by numerous authors. It has been agreed that light is the most important factor that governs growth and must be considered a substrate. Thus, the specific growth rate of microalgae, , is a function of the light available inside the culture: ( ( )). Because there is a different for any , the coupling between the radiative and kinetic models can be done in two ways: The attenuation of light creates a heterogenous light field; the growth decreases along with the increase in the culture's depth. Coupling between radiative models and growth kinetics models is present in many papers and has been investigated by numerous authors. It has been agreed that light is the most important factor that governs growth and must be considered a substrate. Thus, the specific growth rate of microalgae, µ, is a function of the light available inside the culture: µ(G(z)). Because there is a different µ for any z, the coupling between the radiative and kinetic models can be done in two ways:

•
An average irradiance, G avg , can be computed and used to calculate the specific growth rate, µ G avg ; • Various local photosynthetic responses can be calculated µ(G(z)) for various depths z and used to obtain the average photosynthetic response µ.
In this paper, the second approach is considered, even though both approaches give similar results. An inhibition model gives the average specific growth rate: where µ max is 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 ∈ [0 L]). The specific growth rate represents the ratio between the newly formed biomass and the existing one, thus being given in h −1 . To express the biomass concentration, we considered the volumetric growth rate given below: r x is the volumetric growth rate that balances growth (i.e., r xp = the photosynthetic growth rate) and decay (i.e., r xd = the volumetric decay rate). µ d is the specific decay rate, which describes processes such as cell death or respiration (opposite to photosynthesis). µ d is considered constant here.
Batches are done regularly at a constant incident light intensity. Still, many experimental setups allow light variation over a wide range of intensities. In these conditions, by expressing the specific growth rate as a function of light, the biomass is an explicit function of the input, q 0 .
In this paper, we have adopted the hypothesis that the incident light intensity is not constant over the entire control horizon and consequently is denoted by q(t). We emphasize that q(t), which will replace q 0 , is the control input for the BPR considered in the closed-loop control structure proposed in this paper.
The amount of light consumed in the process is calculated by simply integrating the incident light intensity over the interval [t 0 , t]: where A is the lighted surface of the reactor.
In this work, we have considered the PBR dynamic model to be the following two ordinary differential equations: where x 1 (t) is the biomass concentration (in g/L − ) and x 2 (t) the amount of light consumed up to moment t (in µmol/m 2 /s). The nonlinear character of the first differential equation will be clearly expressed in Section 3.1.
The efficiency of a batch can be evaluated by the variable "biomass yield on light energy", denoted by Y. It is defined by the ratio between the newly produced biomass and the amount of light used in the considered interval: where x 1 (t 0 ) is the biomass concentration at the beginning of the batch and V is the working volume of the PBR. Other PBR parameters, such as the specific growth rate and existing biomass, will complete the dynamic model.

Optimal Control Problem
The elements defining an optimal control problem (OCP) and their integration into real applications are presented in [22,23]. In the following, we describe the OCP-called the photobioreactor optimal control (PBROC)-that was solved in this work. The three parts that define the PBROC problem are given hereafter.
(1) The dynamic model of the optimized process As stated in Section 2, the PBR depth is divided into equal k L segments (in our implementation, k L = 100 for L = 4 cm). The corresponding depth values are z i , i = 1, . . . , k L . This discretization involves the calculation of the following constants (see Equation (1)) used for the calculation of the state variables: The PBR's dynamic model consists of two state equations as follows: . . where Equations (8) replace the continuous Equations (1)- (3). The constant C = 3600 × 10 −6 is present in Equation (7) because the light intensity (µmol/m 2 /s) will be finally expressed in mol photons/m 2 /h. We recall that q(t) is the control input for this process. With a consecrated notation, it holds u(t) = q(t) The state variables have the following significance: x 1 (t) = the biomass concentration; x 2 (t) = the amount of light consumed up to moment t. The second state variable x 2 (t) will be considered within the process model only when the PBR productivity must be evaluated.
In our work, we are interested in two PBR parameters that can be considered as output variables: the specific growth rate and the biomass existing in the PBR. These parameters are calculated as follows: SGR indirectly depends on t through the intermediary of x 1 (t) and q(t).
(2) Constraints For the PBROC, the following constraints have to be met: initial conditions: bound constraints: The constants t 0 and t f are the initial and final times of the control horizon; q m and q M are, respectively, the minimum and maximum values of the admissible light intensity. Generally speaking, a dynamic system's evolution might be subjected to other types of constraints: path constraints, terminal equality constraints, algebraic equality, etc.
(3) Optimum criterion In our context, the optimal control seeks the control input u(·) minimizing an objective function I. We recall hereafter its continuous general form: The first part is a Lagrange-type term that measures the quality along the trajectory of the dynamic system; the second part is a Mayer-type term that measures the quality of the trajectory in its final state.
The proposed OCP involves finding out the control variable u(·) that meets all constraints and minimizes the objective function I. In this paper, a value I(u, t 0 , t f ) is associated with each evolution of the system (Equations (6) and (7)): The constant m 0 is an imposed value for the biomass m(t) at the final moment t f . The value x 2 (t f ) is the total amount of light that has radiated the microalgae culture. Our OCP has the optimum criterion J defined below: Equations (14) and (15) correspond to a bilocal optimization problem with a fixed final time. The minimization of the first term corresponds to a small amount of light, while a small value for the terminal penalty means achieving the goal m 0 . Taking into account the second state equation from Equations (7) and (15), we can render the optimum criterion in the form The constants w 1 and w 2 are the new scale factors, whose setting is an important issue. These constants must sufficiently penalize the non-fulfilment of the bilocal constraint but avoid the minimization of the light term falling into eclipse.
For PBR users, it is important to have minimal final biomass, that is This constraint regards the state variable x 1 because it can be written as follows: Remark 1. The control algorithm (containing the Predictor module) will generate a set of system trajectories. The constraint in Equation (17) can be seen as a path constraint defining the admissible trajectories.
Adding the inequality in Equation (17) to the set of constraints in Equations (11)-(13), we render the optimum criterion in its final form: Sensors 2021, 21, 8065 7 of 28 Therefore, we have a weighted combination of the two criteria that transforms our OCP into a standard optimization problem with a single optimal criterion. To summarize, the PBROC can be seen as a procedure with: • Input data: Equations (5)- (13), (17) and (19), all the PBR constructive parameters and constants; • Output data: the optimal control variable q(t), t 0 ≤ t ≤ t f .
NB: This work aims to give a closed-loop solution to the PBROC problem, which is more complex than finding the optimal control input sequence applied from the known initial state.

Receding Horizon Control Structure
For the PBROC problem stated before, we need a closed-loop solution, i.e., an optimal controller that sends a value q(t), t 0 ≤ t ≤ t f toward the BPR in each sampling period. The current value q(t) will be computed based on the current state of the BPR such that the process will meet all the constraints (Equations (5)- (13) and (17)) and the optimum criterion (Equation (19)).
The control structure adopted in this work to achieve the closed-loop solution is RHC (see [16,17]) as presented in Figure 2. RHC involves the calculation of the optimum criterion over the interval [t k t f ], where [t k t k+1 ] is the current sampling period. Therefore, the controller has to predict the future evolution of the process (see [16]). That is why the RHC structure includes a process model, which allows predictions to be made using the current real state x(k) acquired from the process. The controller selects its control output u(k) (for the process, it is the control input) as the prediction's first element and sends it toward the process. The process feedback arrives after a sampling period when the new state x(k + 1) is achieved.
Adding the inequality in Equation (17) to the set of constraints in Equations (11)-(13), we render the optimum criterion in its final form: Therefore, we have a weighted combination of the two criteria that transforms our OCP into a standard optimization problem with a single optimal criterion. To summarize, the PBROC can be seen as a procedure with: • Input data: Equations (5)- (13), (17) and (19), all the PBR constructive parameters and constants; • Output data: the optimal control variable 0 ( ), f q t t t t ≤ ≤ .
NB: This work aims to give a closed-loop solution to the PBROC problem, which is more complex than finding the optimal control input sequence applied from the known initial state.

Receding Horizon Control Structure
For the PBROC problem stated before, we need a closed-loop solution, i.e., an optimal controller that sends a value q(t), t0 ≤ t ≤ tf toward the BPR in each sampling period. The current value q(t) will be computed based on the current state of the BPR such that the process will meet all the constraints (Equations (5)- (13) and (17)) and the optimum criterion (Equation (19)).
The control structure adopted in this work to achieve the closed-loop solution is RHC (see [16,17]) as presented in Figure 2. RHC involves the calculation of the optimum crite- t t + is the current sampling period. Therefore, the controller has to predict the future evolution of the process (see [16]). That is why the RHC structure includes a process model, which allows predictions to be made using the current real state ( ) x k acquired from the process. The controller selects its control output ( ) u k (for the process, it is the control input) as the prediction's first element and sends it toward the process. The process feedback arrives after a sampling period when the new state ( 1) x k + is achieved.

Controller's Predictions
In the following, we use the discrete form of the PBROC problem (Equations (5)-(13), (17) and (19)) in a way that is adapted to the presentation of the next sections. The control horizon is finite [0, H × T], where H is a positive integer and T is the sampling period. The discrete moments t k = k × T will be specified simply by k = 0, 1, . . . , H.
In this work, we consider a realistic case of the control input q(t), which is constant within a sampling period T: Hence, the control input is a step function. The continuous process model (Equations (5)-(9)) can always be converted into a discrete one as follows: where T , and f is a two-dimensional vector function. Any overall control sequence Q 0 meets all constraints and determines the system evolution over the control horizon. It is defined as follows: The optimum criterion (Equation (19)) has its discrete form presented hereafter.
The min operator is applied for the set of sequences that leads to final biomass greater or equal to m 0 .
At each moment k, the controller makes predictions over the predicted horizon [k, H]. A predicted control sequence has the following form: where q(k + i|k) , i = 0, . . . , H − k − 1 is the predicted value for the control input q(k + i) based on our knowledge up to moment k.
The controller uses an algorithm to generate such control sequences and to determine the optimal one. In this work, the optimal predictions will be made by the APSO algorithm.

Remark 3.
Note that the state x(k) is not estimated because it is acquired from the process at every sampling moment. This state, the process model (Equation (21)), and a generated pcs(k) allow estimation of all the intermediary states and computation of the objective function.
The control input q(k + i|k) is kept constant within the sampling period [k + i, k + i + 1], such that a pcs represents a step function.

Remark 4.
It should also be noted that q(k + i|k) = q(k + i); the value q(k + i|k) is a future control input predicted at the present moment, whereas the future real control input q(k + i) is unknown at the prediction moment. We can assert the same thing for the state variables.
Using the current state x(k) = x(k|k), the generated pcs(k) and the process model (Equation (21)), the controller can calculate the corresponding predicted state sequence, pss(k), as defined below: Equation (25) refers to the prediction horizon [k, H] and the state trajectory that starts with x(k). The sequence pss(k) has a length greater by one unit than pcs(k).
Within the RHC loop, the controller has to make an optimal prediction, that is, to calculate the optimal control sequence (ocs) that minimizes the objective function, as below: The optimal state sequence oss(k) is the pss(k) corresponding to the optimal control sequence (ocs). The controller can now establish the current control output using ocs(k). This is the first element of the sequence ocs(k): Table 1 presents an outline of the controller's actions for every sampling period. Its structure is a general one that uses a generic function ("Predictor").  Obtain the current value of the state vector, x(k).
Send u(k) towards the dynamic system 5 Shift the prediction horizon and wait for the next sampling period To determine the optimal prediction in line #2 is not a simple job because of the computational complexity the Predictor faces. Except for the first value of the sequence ocs(k), the remaining part is not used. This fact has a simple explanation: ocs also proves that an admissible final state is accessible from the current state x(k). The same thing must be proved for the next moment (k + 1), starting from a new state variable acquired from the process (Remark 3).

Necessity
Generally speaking, the control output value is subjected to the constraint in Equation (13), where the bounds q m and q M are mainly technological limits. The light intensity q(k) can be supplied between these bounds. For example, in this work, q m = 50 and q M = 2000 µmol·m −2 ·s −1 , which is a very large interval. Therefore the computational complexity of the algorithm making the predictions is very important.
We define the admissibility domain (D k ) of the control input as the interval containing all the values that lead to a favorable process evolution (e.g., meeting certain dynamic constraints). For example, we can consider only the input values that entail positive values for the growth rate of microalgae. We emphasize that the admissibility domain depends on the current state of the process and the "favorable" dynamic constraint at hand.
If the current state of the process entails, for certain physical reasons, a narrower interval D k ⊂ [q m , q M ] such that q(k) ∈ D k , this situation could be used to diminish the computational complexity. The process control input will be looked for in a smaller interval. Hence, the constraint in Equation (13) could be replaced by NB: The process control input may be considered, at the same time, to be the control output q(k) issued out from the controller because it is the connection between the controller and the process.
Using narrower intervals, we expect to improve the Predictor's computational complexity, which in this simulation study is based on the APSOA. This improvement will be proved in the next sections.
This section proposes a soft sensor to reduce the admissibility domain of the process's control input according to its current state. The sensor returns information about the domain D k meeting Equation (30). This is based on: • measurement of the biomass concentration (x 1 (k)) and • the numerical integration of the PM over the next sampling period ([k, k + 1]) for a certain number of light intensity values. These values are considered successively as the control output for the current sampling period: For example, in our simulations, we have adopted nl = 40. The sensor's estimations are based on a real measure x 1 (k) and refer only to the current sampling period to be realistic. In the next sampling period, the sensor will be based on the new real measure x 1 (k + 1).
The sensor makes the integrations and generates its output OUT(k), which the Predictor can use in two different ways to define a smaller interval D k (as described later in this section). Figure 3 shows how the SENSOR is included in the controller. It also shows how data circulate among modules and between the controller and the process in sampling period k. To understand this figure, let us suppose that the biomass concentration x 1 (k) is acquired from the process at the beginning of the current sampling period. All the modules, including the SENSOR, have the data necessary for their calculations and estimations. Finally, the controller computes the output u(k) = q*(k), which is sent toward the process (red lines are used for the information exchanged between the controller and the process).

Output Parameters of the Soft Sensor
To understand what the soft sensor evaluates and what information is returned, we shall analyze the PBR evolution over a period of time, for example, [0, 50] hours. Figure 4 depicts the evolution of the state variables during the interval [0, 50] hours. The initial values of the state variables and all the system parameters are presented in Appendix A. At the moment k = 50, the following values are obtained through numerical integration of the PM: The latter evolves due to the control output u(k) and will produce a new biomass concentration returned to the controller after a sampling period (x 1 (k + 1)).

Output Parameters of the Soft Sensor
To understand what the soft sensor evaluates and what information is returned, we shall analyze the PBR evolution over a period of time, for example, [0, 50] hours. Figure 4 depicts the evolution of the state variables during the interval [0, 50] hours. The initial values of the state variables and all the system parameters are presented in Appendix A. At the moment k = 50, the following values are obtained through numerical integration of the PM: One may wonder if there is a narrower interval, D(50), within which the controller would search for the next control output u(50). To respond to this question, a certain number of values for light intensity u(50) = q(50) = q 0 can be tested. Considering the initial state (Equation (32)), we carried out the numeric integration of the PM over the next sampling period (one hour) for all these values (see Equation (31), with nl = 40). We noted the final values for the biomass concentration, SGR, and biomass yield on light energy.  In the first mode (mode = 1), the light intensity, q(k), avoids the values placed at the On the left, the values for which the specific growth is negative are avoided because the biomass will decrease. On the other side, very large light intensity values lead to growth saturation inside the PBR and a certain inefficiency when considering biomass yield on light energy. In this mode, the light intensity belongs to the following interval:  The soft sensor returns information that can define narrower intervals D k ⊂ [q m , q M ] for the current sampling period. We can call this interval the admissibility domain for the control output at moment k. There are two modes of using this information.
In the first mode (mode = 1), the light intensity, q(k), avoids the values placed at the extremities of the interval [q m , q M ]. On the left, the values for which the specific growth is negative are avoided because the biomass will decrease. On the other side, very large light intensity values lead to growth saturation inside the PBR and a certain inefficiency when considering biomass yield on light energy. In this mode, the light intensity belongs to the following interval:  The second mode of using the sensor (mode = 2) is related to the biomass yield on light energy (Equation (4)). Figure 5 shows Y as a function of light intensity. In our case, the value ( ) 451 opt q k = µmol·m −2 ·s −1 determines a maximum value for Y. This information can be useful in control applications that aim to maximize Y over larger intervals. In the PBROC problem, apparently, this information would not be useful because we have to minimize the amount of light while a specified quantity of biomass is produced. Nevertheless, the strategy to control the PBR using light intensity values around the optimal value ( opt q ) would have a beneficial influence. In this sensor mode, the light intensity belongs to the following interval: In this work, after a few tests, we chose the value 0.2 r = , which offers a sufficient range for the control output and keeps a certain optimal behavior. The simulations proved that this mode produces very good results. Let us remark that if the value of Y is near optimal in the current sampling period, the amount of light is near minimal for the newly produced biomass. This is why the global optimum criterion can be determined using the interval in Equation (34).
The soft sensor's third parameter is the value ( ) opt q k . It belongs to the calling application to compute the extremities of the interval in Equation (34) using an appropriate value r. Finally, the list of parameters returned by the soft sensor is The value q 1 M is a fraction of the light intensity that produces the maximum SGR. In this work, we have considered p = 0.8. Therefore, in our example, a reasonable choice is to consider D(50) = [150 , 1101] µmol·m −2 ·s −1 , like in Figure 5 where the red lines define the two values (q 1 m and q 1 M ). Aiming to use this mode, the first two output parameters of the soft sensor are q 1 m and q 1 M . The second mode of using the sensor (mode = 2) is related to the biomass yield on light energy (Equation (4)). Figure 5 shows Y as a function of light intensity. In our case, the value q opt (k) = 451 µmol·m −2 ·s −1 determines a maximum value for Y. This information can be useful in control applications that aim to maximize Y over larger intervals. In the PBROC problem, apparently, this information would not be useful because we have to minimize the amount of light while a specified quantity of biomass is produced. Nevertheless, the strategy to control the PBR using light intensity values around the optimal value (q opt ) would have a beneficial influence. In this sensor mode, the light intensity belongs to the following interval: In this work, after a few tests, we chose the value r = 0.2, which offers a sufficient range for the control output and keeps a certain optimal behavior. The simulations proved that this mode produces very good results. Let us remark that if the value of Y is near optimal in the current sampling period, the amount of light is near minimal for the newly produced biomass. This is why the global optimum criterion can be determined using the interval in Equation (34).
The soft sensor's third parameter is the value q opt (k). It belongs to the calling application to compute the extremities of the interval in Equation (34) using an appropriate value r. Finally, the list of parameters returned by the soft sensor is According to one of the two modes described above, it is up to the calling application to use this information.

Implementation of the Soft Sensor
The SENSOR function is described by the pseudocode presented in Table 2. To estimate its output parameters, the SENSOR has the biomass concentration as the input variable. Other important information is available through global variables: the current moment, k, sampling period, T, and the PBR's dynamic model SPM (simplified to PM). The latter includes only Equations (6) and (9) concerning the biomass concentration and SGR.
Initialize the light intensity vector (q 0 (l), l = 1, . . . , nl) according to Equation (31). 2 The biomass concentration can be determined online due to a linear correlation between the dry matter and the turbidity of the culture. The dry matter can be determined daily, and the turbidity measured online (see [10]).
After the numerical integration of the SPM, to estimate its output parameters, the SENSOR uses very simple and fast methods whose precision is satisfactory. The greater the value of nl, the more precisely the three output parameters are determined. A very high precision for these parameters is not mandatory because the two intervals D k will be used in a stochastic environment.
Line #1 can be implemented, for example, by a loop that generates the following vector: The numerical integration for the nl light intensities is carried out within lines #3-#6. The function EvalState is called at line #4 and achieves the SPM numerical integration over a sampling period ([0 T]), starting from the initial state X 0 and applying the light intensity q 0 (j). The integration results are the biomass concentration (Bc(j)), SGR(j), and Y(j) at the final moment, T.
Lines #7-#13 calculate q 1 m (k) as being the first light intensity that involves a positive value of SGR. The value q 1 M (k) is determined in lines #14-#22 as being q 0 (j 2 ), to which the SGR that is the closest to Sg80 (defined in line #14) corresponds. The third output parameter is calculated in lines #23-#24. The value q opt (k) ensures the maximum value of Y.
A real implementation of the SENSOR-using the MATLAB language and system-is presented within the files "SENSOR.m" and "XXEvalState.m" inside the folder "PSO_pred_senz".

Brief Description of PSO
PSO is a well-known metaheuristic used in many applications [24][25][26]. To keep the presentation self-contained, we present hereafter a few elements defining this metaheuristic for the readers that are newcomers in the field of PSO. Its detailed structure is given in Appendix B. A particle swarm system aims to optimize an objective function with n decision variables through a simulated movement of the particles within the search space where the objective function is defined. The latter models the "environment" shape, i.e., the altitude of "valleys" and "hills", which characterizes the optimization problem (OP). Each particle "flies" over new regions and updates its information. In the initial version of PSO, the particles communicate among themselves through the intermediary of a global variable, memorizing the swarm's best position "encountered".
As the search process proceeds, the exploration decreases, and the exploitation is intensified. Finally, the particles converge to the global best solution of the OP. The main elements concerning the analysis of PSO algorithms are presented in [24]. The convergence of PSO algorithms is addressed in [23,24]. Results guide the choice of algorithm parameters.
The swarm is composed of N particles, coded through a three-component vector, usually denoted by X i , V i , P best i . Each component is an n-dimensional vector representing the position, speed, and best personal position reached in the search process (see [20]). It holds as: An iterative process simulates the evolution of the particle swarm. P gbest denotes the particles' best position up to the current step, called the "global best" position. When the algorithm ends, that is when convergence is achieved; P gbest is the OP solution.
The PSO algorithm's main action is updating the speed and position of particles at each iteration of the search process. This update is carried out using the following equations: w is inertia weight, C 1 and C 2 are acceleration coefficients, rand 1 and rand 2 are random numbers in the interval [0, 1], and t is the step number.

Adaptive PSO Algorithm
The PSO algorithm used in this paper uses two advances: hybrid topology and adaptation of particle speed. Hybrid topology particle swarm optimization (HTPSO) [19,24] improves the PSO metaheuristic with better communication abilities among particles. It is also enhanced with adaptive updating of the particle speed by a continuous change of some algorithm parameters.
The proposed algorithm will be referred to as APSOA (Adaptive PSO Algorithm), implicitly including hybrid topology.
APSOA also uses a swarm's local topology, regarded as a communication network. The local topology means the existence, for any particle #i, of a "social neighborhood", i.e., a set of 3-5 particles that inform particle #i about their best personal experience. These neighborhoods are settled deterministically or randomly (see [19,24]). Each particle will decide the local best position P lbest i , which is the best experience of the particles belonging to the "social neighborhood" (including the particle itself). It holds that: . . , p d lbest i . . . , p n lbest i ; d = 1, . . . , n; A new term, containing C 3 and rand 3 , appears in the speed equation (see [24]). The updating of the speed and position is performed using the following equations: An efficient technique can enhance the APSOA: adaptation of particle speed. This technique modifies the coefficients C 1 , C 2 , C 3 , and w during the iterative process (see [19]). These coefficients are adapted to the phase of the search process and prepared to help the algorithm's convergence. A linear increase for the coefficients C 1 , C 2 , and C 3 between their minimum and maximum values will adapt the particles' speed. At the same time, the parameter w decreases as follows: T is the estimated or maximum number of steps until convergence. Because the structure of APSOA is well known and has already been presented in many papers, including [19,24], this structure is described in Table A1 in Appendix B. However, some specific characteristics are mentioned hereafter: • APSOA is organized like a function of five input parameters; • The predicted sequence length (h) is an input parameter because RHC involves predicted sequences with different lengths; • k and x 0 define the moment and initial state of the process, for which the optimal sequence is determined.

Predictor Structure
The structure of the Predictor is presented in Table 3 as the pseudocode of the function Predictor_SZ, which has the following input parameters: • k is the current moment when the controller calls the Predictor to make an optimal prediction; • x 0 (k) is the initial state (biomass concentration) acquired from the process; • x 1 m and x 1 M are the two extremities of the interval D k determined by the SENSOR (whatever the utilization mode). Table 3. Algorithm of the Predictor function.
Initialization; /*space reservation for each particle*/ The predicted sequence will have a variable dimension, h, according to line #2, where n corresponds to the total control horizon.
The Predictor has two main tasks: • To set the two vectors' values, the control output limits for the predicted sequences (lines #3-#6). The two vectors (x m and x M ) have the same length as the predicted sequences. The values from the SENSOR are assigned to the first element of these vectors, corresponding to the first control output u * (k). The other limits are q m and q M .

•
To call the function APSOA with the appropriate arguments. The best prediction (ocs) is returned to the controller.
To describe the Predictor's structure simply, APSOA was presented before as a function (a separate program unit). In our real implementation, the APSOA code is included in the Predictor code. Details are given in Appendix B.

Simulation Study
The organization of similar simulations is presented in [27][28][29]. All of them concern applications of computational intelligence [30].

Study Objectives and Preliminaries
All the control tools developed in previous sections are available now to tackle the PBROC problem stated in Section 3.1. This simulation study has a few objectives, as listed below:

1.
To implement and simulate a closed loop based on Receding Horizon Control devoted to solving the PBROC problem.

2.
To study the feasibility of the prediction technique using a metaheuristic algorithm. In our case, we have chosen the Adaptive PSO Algorithm. When metaheuristic algorithms are involved in optimal control, the main impediment often encountered is the computational complexity, which affects the controller's computing time for the control output.

3.
To validate the hypothesis that the computational complexity of the prediction will diminish by introducing the admissibility domain D k for the control output.

4.
To validate the hypothesis that the proposed soft sensor can realistically provide the admissibility domains using the biomass concentration measured from the process.
Objectives 1 and 2 will be analyzed together as well as objectives 3 and 4 because they are intrinsically connected.
The simulation study will be performed with an application that emulates the closedloop functioning. Some modules will be implemented realistically, like in a real-time application, with only a few adjustments devoted to the simulation.

•
The Controller, which by definition includes the PM, is also connected to the real process (PBR). The soft sensor is connected to the process as well. Both modules are realistically implemented but should have connections with the process (PBR). Only these connections are simulated. • In our application, the real process is also simulated using the PM (see Remark 2). The red lines in Figure 3 also show the simulated connections that create the closed loop.

•
The sequence of sampling moments is simulated.

Closed-Loop Implementation
The algorithm RHC_Closed_Loop, which simulates the closed-loop solution of the PBROC problem and will allow us to analyze the proposed tools, is described in Table 4.   18 (1); /* optimal control output q * (k)*/ /*Send u(k) towards the process */ 20 x next ← RealProcessStep(u(k), x 0 , k) ; /*Wait for the next sampling moment and obtain the next state of the process*/ 21 After the appropriate initializations for the PM, APSOA, and simulation, the algorithm is mainly structured through a "while-loop" that follows the sampling moments k = 0, . . . , t f − 1.
The variable mode, set in line #2, specifies how the SENSOR is used. If mode = 0, the SENSOR is not called, and the admissibility domain D k is the widest possible (#11). The SENSOR is called within line #12, as though mode = 1. If it is actually mode = 2, lines #13-#15 make the corrections.
The function "RealProcessStep(u(k), x 0 , k)" simulates the state evolution of the real process. It integrates the PM over the next sampling period, starting from the initial state x 0 and having the step function u(k) as a control input.
Lines #19 and #20 contain comments corresponding to the red lines from Figure 3. The vector "state" has (n + 1) elements and stores the process states for k = 0, 1, . . . , n. Finally, it represents the global oss (see Section 3.3) obtained with RHC.
Details concerning the implementation of the RHC_Closed_Loop algorithm are given in Appendix C.

Simulation of the Tandem SENSOR-Predictor
To better understand how the SENSOR and Predictor modules work together, we present in this section a simulation of their standalone execution (which is not included in the closed-loop control). We calculate the optimal control sequence at the moment k = 0, that is, an ocs with n = 120 control output, considering different values for m 0 . The SENSOR (mode = 2) is called before the Predictor call as below: This relates to lines #12-#18 of the RHC_Closed_Loop algorithm (Table 4). Considering the initial biomass concentration x 1 (0) = 0.36 g/L, the PM is integrated using ocs. Table 5 presents the results.
For each execution, five values are displayed: the target m 0 (g), newly produced biomass (g), value of the optimum criterion (J), amount of Light x 2 (H) (mol photons/m 2 /h), and the number of calls (Ncalls) of the objective function. We can conclude that the target m 0 is precisely achieved in a small number of objective function calls.

Closed-Loop Simulation without SENSOR
In this section, the simulation results of the RHC structure for solving the PBROC problem are presented. We give hereafter the input data for the problem we have solved.
The PM presented in Section 3.1 is characterized by the physical and constructive parameters given in Appendix A. In addition, we have the following initial data:

•
Control horizon: 120 h; t f = 120 h; n = H = 120; • Sampling period: 1 h; T = 1; • Light intensity bounds: q m = 50; q M = 2000 (µmol·m −2 ·s −1 ); • Initial biomass concentration: x 1 (0) = 0.36 g/L. Desiderata: • To fulfil the optimum criterion in Equation (27), where w 1 = 1; w 2 = 10; • To ensure the newly produced biomass: The APSOA has the parameters presented in Appendix B and tuned for this application. Their values result clearly from the first lines of the script RHC_PSO_PBRJ13_SZ.m. Very important for the algorithm's computational complexity is the number of adopted particles in the swarm: Because the prediction is based on a stochastic algorithm, to analyze the simulation results, the RHC_Closed_Loop program was ran 30 times for each instance of the PBROC problem. Practical details about this operation are given in Appendix C.1.
A realistic measure of this application's computational complexity is the number of calls of the objective function during the control horizon. This fact is true, especially when the objective function involves numerical integrations, which have significant complexity and are time-consuming. This measure is also adequate for comparison between different versions of applications. That is why the RHC_Closed_Loop totalizes the number of calls for each sampling period.
The results of this simulation series without SENSOR (controller in mode = 0) are given in Table 6.
For each execution, four values are displayed: the value of the optimum criterion (J), the amount of light (Light) over the entire control horizon (mol photons/m 2 /h), newly produced biomass (g), and the average number of calls (Ncalls). The latter equals the total number of calls divided by the number of sampling periods (n = 120). The sampling periods have very different numbers of calls: the greater the value of k, the smaller the number of calls. However, it is easier to perceive and compare the total number of calls divided by 120 (Ncalls).

Remark 6.
It can be seen that newly produced biomass equals almost exactly 3 g for all executions. This fact is the consequence of how the constraint in Equation (41) is checked. Equation (27) defining the optimum criterion over the prediction horizon [k, H] specifies that only the predicted control sequences that lead to m(H) ≥ m 0 are considered. In our implementation of the objective function (like a programming unit), the constraint in Equation (41) is treated like a trajectory constraint (Remark 1). If the constraint in Equation (41) is not met, the value returned by the objective function is infinity, and the trajectory is not admissible. Otherwise, it equals the sum of the two terms (see file RHC_EvalFitnessJ13.m, which implements the objective function). On the other side, the optimum criterion in Equation (27) looks for a minimum for both terms and their sum. The smaller the amount of light, the smaller the value m(H) − m 0 (for physical reasons). Therefore, naturally, APSOA will find a minimum for which we have m(H) − m 0 ≈ 0.
Because the constraint in Equation (41) is met for all the executions, the produced biomass cannot discriminate among the quality of simulations. Therefore, we shall determine the typical one among the 30 simulations using, of course, the value of the optimum criterion. We shall consider the simulation that produces the closest value to the average optimum criterion as the typical execution. The minimum, average, and maximum values and standard deviation for the optimum criterion (J) are given in Table 7. In our simulation series, there is an execution that yields J typical = 9.226. The other simulation results (line #29 in Table 6) are presented hereafter: J typical = 9.226; Light = 9.226; ∆m = 3.0000; NCalls = 798.25 For each execution, the 120 values of the quasi-optimal control output are recorded for the final simulation of the closed loop. The light intensity for the typical execution is depicted in Figure 6a.
In our simulation series, there is an execution that yields Jtypical = 9.226. The other simulation results (line #29 in Table 6) are presented hereafter: Jtypical = 9.226; Light = 9.226; Δm = 3.0000; NCalls = 798.25 For each execution, the 120 values of the quasi-optimal control output are recorded for the final simulation of the closed loop. The light intensity for the typical execution is depicted in Figure 6a. This control action involves the typical evolution of the state variables in Equations (6) and (7) and the newly produced biomass depicted in Figure 6b.
The simulations have proved that the APSOA converges for all executions. This fact is a characteristic of this metaheuristic with the considered advances. Despite the computational complexity of the first sampling period, when the pcs's length is maximal, the controller calculates the control variable in a few decades of seconds (depending on processor speed). This amount of time is very satisfactory for a sampling period of 1 h. Hence, the control structure could be a solution even for real-time control. This control action involves the typical evolution of the state variables in Equations (6) and (7) and the newly produced biomass depicted in Figure 6b.
The simulations have proved that the APSOA converges for all executions. This fact is a characteristic of this metaheuristic with the considered advances. Despite the computational complexity of the first sampling period, when the pcs's length is maximal, the controller calculates the control variable in a few decades of seconds (depending on processor speed). This amount of time is very satisfactory for a sampling period of 1 h. Hence, the control structure could be a solution even for real-time control.

Closed-Loop Simulation Using the Soft Sensor (mode = 1)
The RHC_Closed_Loop program with SENSOR (mode = 1) was ran 30 times for the same instance of the PBROC problem as in the previous section. Practical details about this simulation are given in Appendix C.2. The results of this simulation series are given in Table 8. Table 8 has two particularities. The first one concerns column ∆m, which always equals exactly 3 g. The second one is the consequence of the first one; columns J and Light are identical because the second term of J is practically null. The constant w 2 determines these particularities because it privileges the minimization of J's second term.

Remark 7.
A smaller value of w 2 entails a relaxation for the minimization of the second term; finally, the value of J will be greater, and the values of J and Light will be different. This relaxation will lead to a relaxation of the searching process and, consequently, a smaller number of objective function calls. Hence, we have a spare strategy to diminish the number of calls, choosing a smaller value of w 2 . The price to pay is greater light and a greater newly produced biomass (to a small extent). Table 9 presents the new statistics for the optimum criterion (J). The average value is practically the same as in Table 6, but what is noteworthy is the standard deviation, which is smaller by 43%. The new values of J are positioned in a smaller interval. The difference between J max and J min equals 0.292, smaller than 0.484 as in Table 6. This is the consequence of the controller's commands, which belong to smaller intervals D k . The typical simulation produces the following results (line #23 in Table 8): J typical = 9.271; Light = 9.271; ∆m = 3.0000; NCalls = 661.375.
The average value of Ncalls from Table 8 is 734. The typical simulation of the closed loop with SENSOR (controller in mode = 1) is illustrated in Figure 7.   Remark 8. The average number of calls of the objective function equals 734, which is smaller than 791 by 7.2%. This decrease is expected, and although not spectacular, it proves the improvement in computational complexity due to the soft sensor. The efficiency of the soft sensor is confined to a certain extent because it has a limited action. At moment k, the pcs covers the interval [k, H] and has the form given by Equation (24). The soft sensor improves the bounds only for the first component pcs(k); the other H-k-1 elements will have the initial bounds given by Equation (13). The SENSOR's estimations are based on a real measure x 1 (k) and could not be extended to the other pcs components because the next real states are unknown. The estimations refer to the current sampling period to keep them realistic.

Remark 9.
The increase in computation complexity is negligible when using the SENSOR because the PM's integration covers only the current sampling period. It should be noted that the controller already acquires the biomass concentration to implement the closed loop.

Closed-Loop Simulation Using the Soft Sensor (mode = 2)
This time, the SENSOR calculates the best biomass yield on light energy for the current sampling period. The RHC_Closed_Loop program with SENSOR (mode = 2) was ran 30 times for the same instance of the PBROC problem as in the previous sections. Practical details are given in Appendix C.3 about how the simulation results were generated. The results are presented in Table 10.  Table 10 has only two columns because the light values equal the J values and the newly produced biomass is always 3 g, as in Table 7. The minimum, average, maximum and typical values and standard deviation for the optimum criterion (J) are given in Table 11. The typical simulation produces the following results (line #6 in Table 10): J typical = 8.429; Light = 8.429; ∆m = 3.0; NCalls = 660.
The average value of Ncalls is 686.8, smaller than 791 from Table 6 (without SENSOR) by 13.1%. This decrease in computational complexity is even more important when we consider the argumentation of Remark 8.
The typical simulation of the closed loop with SENSOR (controller in mode = 2) is illustrated in Figure 8.
To explain the SENSOR's good result, we recall that we want to minimize the amount of light on the entire control horizon while the constraint in Equation (41) is fulfilled. This desideratum is equivalent to maximizing Y over the entire control horizon while meeting the constraint in Equation (41).

Remark 10.
Maximizing Y in each sampling period does not imply the maximum value of Y over the interval [0, H]. It is a "greedy" strategy that, in general, does not ensure the global optimum. Nevertheless, if the value of Y is near optimal in the current sampling period, the amount of light is near minimal in relation to the newly produced biomass. At the same time, the predicted sequence guarantees that the constraint in Equation (41) is met because it generates an admissible trajectory.
After a few exploratory tests of Equation (34), we chose the value r = 0.2, which offers a sufficient range for the control output to face the minimization of J.
Remark 11. The appearance of Figure 8, which entails a certain monotony of the control output, is due to the monotony of the value q opt as a function of initial biomass concentration, which is considered by the sensor within its estimations.

Conclusions
For the PBR and its PM, which is well known and has been validated in previous papers, we have stated the PBROC problem. The implementation of its solution requires metaheuristic-based predictions. The Adaptive PSO Algorithm was chosen for its convergence speed. Despite the ability of APSOA to efficiently generate optimal predictions (in a reasonable time), we have proposed a soft sensor to cope with a very large prediction horizon (e.g., 120).
The main contribution of our work is the soft sensor, which was conceived to diminish the Predictor's computational complexity. The SENSOR is based on measurement of the biomass concentration and numerical integration of the process model over the current sampling period. The returned information concerns the specific growth rate of microalgae and the biomass yield on light energy. Our proposal generated two modes of using the soft sensor to reduce the admissibility intervals (D(k)) for the process control input (u(k)).
The simulation study proved that our main contribution, the SENSOR, reduced the computational complexity significantly, expressed by the number of calls of the objective function.

Conclusions
For the PBR and its PM, which is well known and has been validated in previous papers, we have stated the PBROC problem. The implementation of its solution requires metaheuristic-based predictions. The Adaptive PSO Algorithm was chosen for its convergence speed. Despite the ability of APSOA to efficiently generate optimal predictions (in a reasonable time), we have proposed a soft sensor to cope with a very large prediction horizon (e.g., 120).
The main contribution of our work is the soft sensor, which was conceived to diminish the Predictor's computational complexity. The SENSOR is based on measurement of the biomass concentration and numerical integration of the process model over the current sampling period. The returned information concerns the specific growth rate of microalgae and the biomass yield on light energy. Our proposal generated two modes of using the soft sensor to reduce the admissibility intervals (D(k)) for the process control input (u(k)).
The simulation study proved that our main contribution, the SENSOR, reduced the computational complexity significantly, expressed by the number of calls of the objective function.
To conduct the simulation study presented in this paper, we implemented the proposed algorithms using the MATLAB language and system. The Supplementary Materials (attached to this article) include scripts and files that show how the closed loop, Controller, SENSOR, and Predictor (including the APSOA) can be implemented with good results. Many algorithms and their implementation could be used in a future real-time application, at least as a starting point.
For future work, we could investigate the conception of sensors for other processes as a general idea, aiming for the same objective of reducing the admissibility intervals. The difficulty would be in disclosing the biophysical aspects that justify the limitation of the admissibility domains according to the current state of the process. The initial biomass concentration x 1 (0) = 0.36 g/L.
The simulation in Section 3.4.3 was carried out using q 0 = 502.3 µmol/m 2 /s (the light intensity that radiates the PBR).
The mathematical model was validated (see [10]) on microalgae strains such as Chlamydomonas reinhardtii and Scenedesmus quadricauda (see [9]) cultivated in laboratory flat-plate photobioreactors in similar conditions. The biomass concentration, measured offline as dry mass, can be correlated with turbidity, measured online by a sensor. The PBRs are lighted on one side by dimmable LED panels, which the process computer can control. The calibration of the LED panel can be done with a portable light meter; the incident light intensity is measured inside the reactor at z = 0.

Appendix B
Description of APSOA and Implementation of Predictor Table A1. Description of adaptive particle swarm optimization algorithm. Update the particles' speed using Equation (39) 19 Speed_limitation v d i (t + 1) ) 20 Position_ limitation x d i (t + 1) /*using x m , x M */  The APSOA code is included in the Predictor code in our implementation, which is the script "RHC_Predictor_SZ.m". The objective function is implemented within the file "RHC_EvalFitnessJ13.m". Both files are inside the folder "PSO_pred_senz".