Dual-Input Slope Seeking Control of Continuous Micro-Algae Cultures with Experimental Validation

Featured Application: The production of algal biomass or of a product associated with biomass growth can be optimized by manipulating the dilution rate and the incident light intensity in indoor photo-bioreactors. However, model-based optimization is delicate in view of the time and experimental efforts required to develop sufﬁciently accurate dynamic models. Extremum seeking provides an alternative real-time optimization approach, where prior model knowledge is not required, and only a measurable (or estimable) performance index is needed online. Abstract: This paper investigates the application of adaptive slope-seeking strategies to dual-input single output dynamic processes. While the classical objective of extremum seeking control is to drive a process performance index to its optimum, this paper also considers slope seeking, which allows driving the performance index to a desired level (which is thus sub-optimal). Moreover, the consideration of more than one input signal allows minimizing the input energy thanks to the degrees of freedom offered by the additional inputs. The actual process is assumed to be locally approachable by a Hammerstein model, combining a nonlinear static map with a linear dynamic model. The proposed strategy is based on the interplay of three components: (i) a recursive estimation algorithm providing the model parameters and the performance index gradient, (ii) a slope generator using the static map parameter estimates to convert the performance index setpoint into slope setpoints, and (iii) an adaptive controller driving the process to the desired setpoint. The performance of the slope strategy is assessed in simulation in an application example related to lipid productivity optimization in continuous cultures of micro-algae by acting on both the incident light intensity and the dilution rate. It is also validated in experimental studies where biomass production in a continuous photo-bioreactor is targeted.


Introduction
Extremum Seeking Control (ESC) is a Real-Time Optimization (RTO) technique used to drive a measurable process performance index to its optimum following direct input adaptation [1]. In the classical scheme proposed in [2], a sinusoidal dither signal ensures the persistency of excitation of the input and filters are used to extract the gradient information from output measurements. This gradient is forced to zero on average so as to converge to a close neighborhood of the optimal operating point. An interesting feature of ESC is that limited a priori knowledge about the process under consideration is required in order to perform RTO. Although this approach has become popular following the stability proof provided in [3] and the intrinsic simplicity of the scheme [4][5][6], it suffers from the three time-scale separation (implying slow convergence) and its heuristic tuning procedure.
Various approaches have been proposed to improve the convergence rate. In [7], a proportional integral scheme is introduced to improve the transient behavior while [8] uses a recursive least square (RLS) algorithm to estimate online the parameters of a linear model approximating the process input/output static map. However, in [9], the authors point out that the latter approach generally leads to a biased estimation of the optimal operating point when dealing with dynamic systems. It is therefore necessary to account for those dynamics, and in order to preserve simplicity and avoid the need for extensive prior knowledge about the process, a black-box approach is preferable, for instance using Wiener-Hammerstein models [10,11].
Extension of ESC to multiple-inputs-single-output (MISO) process optimization has also provided interesting results both in simulation and experiments. In [12], the classical scheme is extended to achieve maximum power tracking of PV panels with an experimental validation. In [13], an extended Kalman filter is used to estimate the gradient and Hessian of a static quadratic dual-input-single-output (DISO) map to remove thermo-acoustic instabilities in an atmospheric combustor test rig. In applications such as PV panels or wind turbine power generation, it is sometimes required to operate at a given sub-optimal operating point so as to adjust the production to the actual load/demand. In this regard, slope seeking control (SSC) [14] provides the possibility of optimal and suboptimal setpoint tracking. In [15,16], the authors propose an adaptive SSC for a class of SISO dynamic processes, and a numerical validation is provided with biomass productivity optimization in microalgae cultures.
In this work, an extension of adaptive slope seeking a DISO dynamic processes is proposed. To this end, the process under consideration is assumed to be approachable at least locally with a quadratic block-oriented model in the form of a Hammerstein model. Adaptive slope seeking allows driving the process either at the extremum of the measurable performance index or an intermediate operating point specified via the slope of the performance index. The degrees of freedom offered by the additional input can be used in various ways, for instance, to minimize the input energy. A recursive extended least squares (RELS) algorithm is used to estimate online the model parameters. To alleviate the tedious and generally heuristic integrator gain tuning procedure of the classical filter-based ES/SOP control, a pole placement controller (PPC) is used, which was introduced by the authors in [15,16].
The performance of the slope seeking strategy is first assessed in a simulation example related to lipid productivity optimization in continuous cultures of micro-algae by acting on both the incident light intensity and the dilution rate. During the last decade, microalgae cultures have indeed received some attention from the real-time optimization research community (see for instance [17][18][19][20][21]) .
The proposed strategy is then further validated in a lab-scale experimental study of the cultures of the micro-algae Scenedesmus obliquus in a flat-panel photo-bioreactor (PBR).
The paper is organized as follows: Section 2 introduces the main assumptions regarding the dynamic process under consideration and describes the dual-input adaptive slope seeking strategy. In Section 3, the estimation algorithm used to infer the static map and dynamic model parameters is presented. Section 4 is devoted to slope reference generation and proposes two algorithms offering the possibility to minimize one or both inputs. Section 5 presents a pole placement controller applied to drive the process to the selected set-point. In Section 6, the performance of the slope-seeking strategy is assessed with a simulation example related to the optimization of the lipid production in continuous cultures of micro-algae in photobioreactors, while an experimental validation is proposed in Section 7. Finally, some conclusions and perspectives are drawn.

A DISO Adaptive Slope-Seeking Strategy
The control strategy is built upon the assumption that the process can be described by a block-oriented model, which is made of the interconnection of nonlinear static building blocks and linear time-invariant dynamic building blocks. A great variety of configurations are possible, including series and parallel structures. In this study, we restrict ourselves to the classical Hammerstein model, where a nonlinear static block followed by a linear dynamical model represent process and sensor dynamics. Such block-oriented models are often used to model nonlinear systems [22], and in the following, we will consider a class of systems that can be modeled locally using quadratic static building blocks and second-order linear time invariant dynamic blocks: Assumption 1. The process input/output static map can be approximated locally by a quadratic function of the form: where u represents the input vector and x the output of the nonlinear static block, A ∈ R m×m , b ∈ R m×1 , c ∈ R and u ∈ R m×1 , with m being the number of inputs.
This first assumption implies that (i) driving the system to a stable steady-state manifold is feasible and that (ii) the system in the steady state is likely to be a convex/concave function of the input. In bioprocess applications, these statements find their justifications in the bifurcation analysis developed in [23], which can be considered as a basis for further investigations of bioprocesses [24]. Assumption 2. The actuator, process and sensor dynamics can be approximated locally by linear time-invariant dynamic blocks. In particular, they can be described by proper and stable second-order transfer functions with unitary steady-state gains: with K = 1+β 1 +β 2 1+γ 1 +γ 2 and z −1 the back-shift operator such that z −1 x(k) = x(k − 1).
The relatively slow dynamics of (bio)processes allows considering the latter in steady state and to consequently isolate the fast actuator dynamics in a separate block. The loworder approximation of the process nonlinearity and dynamics is justified by the fact that an adaptive control approach will be developed, where the parameters of the block-oriented models are estimated online. The estimation of these parameters requires persistency of excitation of the input signals, which corresponds to the following assumption: Assumption 3. Input signals are persistently exciting and allow closed-loop parameter estimation of the block-oriented models.
In practice, this is achieved using appropriate dither signals in the gradient estimation algorithm, such as pseudo-random binary sequences [25,26]. This point will be further elaborated in the subsequent sections.
Using the notation of Assumptions 1 and 2, the quadratic DISO Hammerstein model can be described by the following state-space model: where J is the measured performance index. A last assumption is now required regarding the existence of a unique extremum (maximum) point of the measurable performance index: Assumption 4. The performance index is concave, or in other words, there exists a unique set of inputs u opt ∈ R m×1 such that: ensuring that there exists a stabilizing input vector u opt driving the process performance index to its maximum J max .
As discussed in [27], this restriction can be relaxed in the multi-modal case provided an adequate dither signal design.
In the adaptive SSC depicted in Figure 1, a RELS algorithm is used to estimate online the static map matrices [A, b, c], which feed a slope reference generator [15] whose role is to convert a desired reference J re f into a corresponding gradient reference ∇J re f . This allows driving the process to the optimal performance index (where ∇J = 0) or to an intermediate, sub-optimal index (where ∇J = 0). The regulation is effected by comparing the gradient estimate ∇J provided by the recursive estimator to the slope reference. The set-point deviation = ∇J − ∇J re f has to be driven to zero by a controller such that: Note that the optimal value of the performance index does not need to be known a priori. The controller will achieve optimal performance by enforcing ∇J re f = 0. When a suboptimal operating point is selected, a non-zero reference is specified for the slope of the performance index. It is in principle possible to operate the plant to the left of the optimum (positive slope) or to the right (negative slope), depending on the process dynamics and stability behavior. As none of these values are a priori known, this procedure will require some form of trial and error by the process operator. It should also be noticed that (5) is only valid if Assumption 4 is locally verified, i.e., in a concave neighborhood. However, information about the achievable optimum and the sensitivity in the vicinity of the optimum is quickly collected in the course of the first experiment. In addition, when a suboptimal operation point is selected, the energy of one or both input signals can be reduced. This is interesting from an economical point of view if the considered process consumes expensive resources. This latter point will be elaborated in the future.
A detailed stability analysis of this scheme is outside the scope of this paper, but stability should be guaranteed under Assumptions 1-4. More details on stability analysis of discrete-time ESC can be found in [2,28], and some results for multivariable ESC can be found in [29,30].

The Recursive Estimation Algorithm
This section develops a Recursive Extended Least Squares (RELS) estimator of the parameters of the quadratic Hammertein model.

Basic Concepts
Consider the output-error auto-regressive (OEAR) system [31]: where y s = ψ T s θ s is the noise-free output of the process (in the previous developments, it represents the noise-free performance index J) and is the output of an Auto-Regressive (AR) process with ν, a white noise with zero mean and variance σ 2 , and C(z −1 ) = 1 + c 1 z −1 + c 2 z −2 + · · · + c n c z −n c is a stable polynomial of order n c . Equation (6) can be rewritten as: where w(k) = ψ T n θ n is the noise regression model, where In the following, it is assumed that n c = 2. Given structure (8), the unknown parameters θ can be estimated online using a RELS algorithm [31] in the form: where e(k) is the a priori estimation error, y meas (k) the measured output, and P(k) the covariance matrix of the parameter estimation error θ − θ(k) at sample k. As a result, the RELS computes online estimates of both the noise-free and noise model regression parameters contained, respectively, in θ s and θ n . The forgetting factor λ determines how fast past data are disregarded, and typically 0.85 ≤ λ ≤ 0.99. An estimation of η(k) is provided by the a posteriori estimation error The RELS algorithm is generally coupled with a constant trace algorithm, ensuring a sufficient correction gain to track slowly time varying parameters. A regularization of the covariance matrix P(k) is obtained as: In addition, a projection algorithm [32] can be added to enforce estimator stability by introducing prior knowledge about some parameter bounds: where : The domain D is defined for i ∈ 1, . . . , n θ , n θ being the number of parameters such that:

Application to the Quadratic Hammerstein Model
Using the quadratic Hammerstein model introduced in (3) with m = 2 and without coupling terms (we consider a 12 = 0 to simplify the developments and because we observed that these terms have little impact on the performance of the controllers), , the development of (3a) gives: This expression can be inserted into (3b), which yields: This equation can be cast in the form of a product between a regressor of informative signals ψ s and the parameter vector θ s : The model parameters can be extracted from the vector θ s in (17c) as: The second-order transfer function representing the system dynamics can be considered with either two, one, or no zeros (which implies, respectively, that γ 1 = 0 and γ 2 = 0, or γ 1 = 0 and γ 2 = 0, or γ 1 = γ 2 = 0). Let us denote RM n p m z the regression model with n p poles and m z zeros. Table 1 lists some possible simplifications.

Persistency of Excitation
To avoid practical parameter identifiability issues using the RLS algorithm, it is recommended to use an input signal with a sufficiently rich frequency bandwidth as discussed in [25,26]. While using a pseudo-random binary signal (PRBS) with limited amplitude and sufficient signal to noise ratio is generally advised for linear systems, amplitude-modulated pseudo-random binary signals and multi-levels pseudo-random signals (MPRS) have been shown to be more appropriate for nonlinear system identification [33,34].
A PRBS is a periodic and deterministic signal with white-noise-like properties [35]. This signal is generated by adding a sequence of shift register specific outputs in a XOR gate in closed-loop. The sequence length is computed as L = 2 N − 1, where N is the number of shift registers in the sequence. The sequence duration is obtained as L T prbs with T prbs = p T s , the PRBS sampling period expressed as a multiple p ∈ N + of the process sampling time T s . Further details about PRBS generation can be found in [31,35].
In this work, it is proposed to generate MPRS by adding three PRBS with different bandwith (sampling time) and shift register sequences inducing a longer period and both amplitude and holding time modulations. Furthermore, this approach allows maximal amplitude and minimal holding time specifications. For DISO processes identification, using such design procedure may also help reducing spectral overlap between inputs when required.

Slope Reference Generator
As introduced in Section 2, the purpose of the slope reference generator (see Figure 1) is to convert the desired performance index set-point J re f ≤ J max into corresponding Consider the quadratic static map, where the dependency on time k is omitted to simplify the notation, and more importantly, the actuator/sensor dynamics is neglected (as we are interested in generating the slope reference signals). where The cost function gradient is given by: A set U 2 (u 1 ) of control inputs u 2 driving the process performance index J to a desired value J re f can be computed as the set of solutions of the second order polynomial Equation (19): From an operational and economical point of view, it can be interesting to drive the process output to a desired set-point while minimizing one of the inputs or both.
Considering the positiveness of u 2 , the gradient of U 2 with respect to u 1 reads: The value of u 1 minimizing u 2 such that J = J re f is computed by setting which yields: (25) and the corresponding minimal value of u 2 is computed by replacing u 1 by u s 1 in (22). It should be noticed that the positiveness of the inputs u 1 and u 2 , combined with (25) and the concavity of J (a 11 and a 22 are strictly negative parameters), imposes that b 1 and b 2 are positive. Therefore, the derivative (23) is negative when u 1 < u s 1 and positive when u 1 > u s 1 , and (22) presents a minimum in u s 1 .
These results are compiled in Algorithm 1, which summarizes the design procedure for slope reference generation either for extremum seeking (ESC) or suboptimal setpoint (SOP) with minimization/maximization of one input (u 2 ). Algorithm 1 Adaptive slope reference generation algorithm for DISO processes ESC and SOP with minimization/maximization of one input if Objective == 'ESC' then 3: Select a desired Performance index set-point: J re f (k) .

6:
Compute u s 1 k , the value of u 1 minimizing (maximizing) u 2 k such that J = J re f (k) given parameter estimates at iteration k: and substitute in (20) and (22) to obtain: Deduce the gradient reference: ∇J re f (k) = 2 a 11 u s 1 k + b 1 2 a 22 u min 2 k + b 2 9: endif 10: return ∇J re f (k) Algorithm 2 considers a suboptimal setpoint tracking with minimization of both inputs. Starting from arbitrary initial conditions u 1 0 and u 2 0 , the algorithm iterates to estimate the minimal values u min Select a desired Performance index set-point: J re f (k) .

3:
Compute u min 1 k and u min 2 k , respectively, based on u min 2 k−1 and u min 1 k−1 as: Deduce the gradient reference: return ∇J re f (k)

Controller Design
The model under consideration requires two parallel SISO controllers: one for the regulation of the input/output pair u 1 / ∂J ∂u 1 and another for u 2 / ∂J ∂u 2 . To this end and to replace the classical ESC integral controller of Figure 1, a pole placement controller (PPC) introduced by the authors in [15,16] is investigated.
At time k, the gradient of the performance index with respect to the input signal u i , i ∈ {1, 2} is given by: This expression can be cast into the form of a state-space model with the following correspondence: To ensure zero steady-state error under step disturbances, an integral action is introduced and pole-placement allows imposing a first-order closed-loop dynamic behavior. The dominant pole imposes the settling time t s and the non-dominant pole is selected to be much faster, so as to allow for a fast transient. Using the certainty equivalence principle [31,36], a ii is replaced by its estimate a ii and the controller gain is updated at each innovation.
Although the control algorithm depends on the estimation provided by the recursive algorithm, the PPC design is based on only one parameter a ii , which allows one to maintain some robustness level. Further, placing the closed-loop poles allows taking the a priori neglected dynamics into account. There is indeed a three time-scale defined by the estimator (the fastest), the linear dynamics, and the controller (the slowest), which should be respected during the pole placement procedure. The following simulation study will highlight these elements.

A Simulation Study: Optimization of Micro-Algae Lipid Productivity
Microalgae growth in a bioreactor can be described by the Droop model [37], which assumes that extracellular nutrient S is stored in a so-called internal quota Q, before being used for biomass X growth. The Droop Model is usually the cornerstone of more elaborated models accounting for light-related phenomena [38].
For the model under consideration, nutrient uptake ρ is described by a Michaelis-Menten type kinetics with internal quota limitation: whereas the growth rate is modelled by Droop kinetics and photo-inhibition is incorporated through Haldane kinetics: In these expressions, K S is the half-saturation constant with respect to substrate S, Q min is the quota threshold below which no growth is possible, and Q max is the maximum nitrogen quota. ρ max , µ max , K sI and K iI are, respectively, the maximum uptake and growth rates, the half-saturation and the inhibition constants with respect to incident light intensity.
To keep the model simple, photo-acclimation is described by a first-order dynamics: Mass balances in the chemostat lead to the following set of ordinary differential equations: where S in represents the concentration of inorganic substrate (nitrogen) in the inlet flow and D the dilution rate. The model parameters are partly taken from [39] and listed for convenience in Table 2. A classical performance index for continuous bioreactor is biomass productivity: However, this index does not reflect the lipid productivity. As lipid and carbohydrate production increases under nitrogen stress, and thus under internal quota limitation, the following performance index is introduced: which combines biomass productivity and internal quota depletion. The objective of ESC is to find in real-time and without prior knowledge about the model, with the inputs (I, D) maximizing the performance index P l . Figure 2 and Table 3 provide the results of a bifurcation analysis. The concave steadystate manifold validates all the assumptions listed in Section 2. A Hammerstein model is therefore a valid representation even though other block-oriented structures could be considered as well (see for instance the discussion in [16] for the SISO case). These results also show that the optimum of P l corresponds to a more concentrated outflow (higher biomass concentration for a smaller dilution rate) than the optimum of P and a smaller internal nitrogen quota (which can be linked to the production of lipid). For lipid extraction, it is desirable to work with a concentrated solution to improve the operations [40]. Furthermore, Figure 2 shows that the extremum of P (32) corresponds to operating conditions close to the wash-out of the reactor. Conversely, the extremum of P l (33) offers a comfortable security margin.

Asymptotic Observer Design
In contrast to P, P l is not directly measurable as Q is an intracellular variable. However, X and S can be measured online using NIR and UV probes combined with appropriated spectrometers [41], and a software sensor can be developed for Q. To not lose the advantage of ESC and require prior model knowledge, a robust estimator has to be designed notwithstanding the lack of knowledge regarding the kinetic structures and model parameter values. Such a robust estimator of the internal quota has been developed by [42] using super-twisting observers. Here, for the sake of simplicity, an asymptotic observer (AO) is proposed as follows:˙Q =ρ −μQ (34) whose rate of convergence is linked to the growth rate. Indeed, defining the estimation errorQ = Q −Q:˙Q It is therefore necessary to develop estimators of the reaction rates ρ and µ. This can be achieved in a simple way using forward discretization (with sampling period T s ) of the biomass Equation (31a): with ψ ρ = T s X k−1 , θ µ = µ. In the same way, a forward discretization of the substrate Equation (31b) yields: ψ ρ = −T s X k−1 , θ ρ = ρ. RELS with noise model and constant trace algorithm (Equations (6)-(11)) can be used to estimate online the reaction rates considered as slowly varying parameters. This pragmatic approach is less rigorous than the solution proposed in [42] but works quite well in practice and is preferred here for its great simplicity. Moreover, even if a new online measurement is required, it can be achieved using an available (at reasonable costs) substrate probe. This approach is preferred to a dynamic model, which, in essence, would be uncertain and require the same types of adaptation schemes.

Numerical Results
Initial conditions are set as X(0) =1 gX·L −1 , S(0) = 0.1 gN·L −1 , Q(0) = 0.05 gN·gX −1 , D(0) = 0.4 d −1 , I(0) = 600 µmol/(m 2 ·s), and S in = 0.1 gN·L −1 . The initial estimates for the AO and RELS estimating ρ and µ are, respectively, Q(0) = 0.07 gN·gX −1 , ρ(0) = 0.1 gN·gX −1 ·d −1 and µ(0 The AO is first run together with the two RELS for ρ and µ during the first 10 days in open loop. After this initialization period, the adaptive slope-seeking strategy is activated. SOP is first performed following a step sequence [t switch , P l ] = {(10, 8)(15, 6)(20, 7)(25, 8)}. Then, ESC is performed during the rest of the simulation starting from t = 30 d. The controller settling time is set to t s = 2 d, and the second eigenvalue of each pole placement controllers λ c = −100 d −1 . The constant trace algorithm is defined with k 0 = 1000 for the RELS of the adaptive slope seeking algorithm and the two RELS estimating ρ and µ. Measurement noises on X and S are defined as white noises with σ 2 X = 0.007 2 and σ 2 S = 0.001 2 (this corresponds to about 1 % of the variable ranges, which is realistic). The excitation signals are provided by MPRS obtained by adding 3 PRBS as suggested in Section 3.3. For the incident light intensity, T prbs i = p i T s , with p 1 = 2 (amplitude = 10 µmol/(m 2 ·s)), p 2 = 5 (amplitude = 25 µmol/(m 2 ·s)), p 3 = 11 (amplitude = 15 µmol/(m 2 ·s)), whereas for the dilution rate p 1 = 3 (amplitude = 0.02 d −1 ), p 2 = 8 (amplitude = 0.01 d −1 ), p 3 = 13 (amplitude = 0.02 d −1 ). Figure 3 shows that the AO provides a fast and accurate estimation of the internal quota despite noisy estimates computed by the 2 RELS. It should also be highlighted that the convergence of the AO is ensured by the good cell growth rate maintained all along the culture.  Figure 4 shows that the productivity reference trajectory (i.e., a succession of step changes) is well followed by the controller (see the evolution after 10 days of open-loop operation). The incident light intensity is set at the optimum, and its value quickly settles, whereas the dilution rate is varied so as to achieve the different productivity levels. In this figure, the blue dashed curves correspond to the suboptimal open-loop signals determined by the slope generator, while the black continuous lines represent the actual signals generated by the controller, which also include the PBRS dither signals. Figure 5 shows the phase plane plot evolution, while Figure 6 shows the gradient evolution.     Figure 8 shows the lab-scale photobioreactor where the experiments are carried out to assess the performance of SSC. The PBR is equipped with several plunging probes for the online monitoring of the micro-algae concentration as well as the control of the temperature and pH of the medium. To optimize the homogeneity of the micro-algae culture, air-lift flows are achieved on both PBR sides. The PBR vessel is immersed in a water bath whose temperature is regulated by a chiller/heater and a pump that circulates the water to the chiller and back to the water bath. pH is regulated by injecting carbon dioxide. The incident light is generated by a horticultural light fixture with LED lights providing appropriate intensity for the micro-algae growth. Quantum sensors are used to measure the light intensity at the PBR surface, and a PI controller is designed to regulate the incident light intensity. To ensure a continuous mode operation, two peristaltic pumps control the inlet and outlet flows so as to achieve some required dilution rate. The control algorithm is implemented on a LabView platform for real-time computation, and Labjack cards are used for data acquisition.

An Experimental Validation: Optimization of Micro-Algae Biomass Productivity
To validate the slope-seeking control (SSC) in a practical setting, the dilution rate and light intensity are used as manipulated inputs, while the measurable productivity is considered as the performance index. It is therefore a simpler situation than the one considered in the previous simulation experiments. Here, there is no need to measure or estimate the internal quota Q, as it is not involved in the cost function. A sample experiment is described in what follows, which is conducted over 8.6 days, with a sampling period of T s = 0.1 d = 2.4 h. After inoculation, the culture must first settle until biomass concentration reaches a steady-state corresponding to a constant dilution rate D(0) = 0.1 d −1 and a constant average incident light I av = 500 µE/(m 2 ·s). The inlet substrate concentration is S in = 0.1 gN·L −1 . Temperature and pH are kept constant, respectively, at 26 • C and 7.3. After the steady state has occurred (in this case for X(0) = 0.88 gX·L −1 ), the slope seeking control is activated with a productivity setpoint of P = 0.2 gX/(L·d). The dilution rate and incident light are then manipulated by the controller so as to achieve the desired level of productivity. After 7.6 days, the setpoint is changed to P = 0.1 gX/(L·d) to challenge the tracking performance of the SSC. The time evolution of the several signals is shown in Figures 9 and 10. During the first phase, i.e., when the setpoint P = 0.2 gX/(L·d) is imposed, the system converges within 1 day to the desired productivity level, but the steady-state is only truly established after 3 days, as shown by the biomass concentration evolution in Figure 10. The steady-state values of biomass, dilution rate, and incident light are, respectively, approximately X = 0.5 gX/L, D = 0.35 d −1 and I = 380 µE/(m 2 ·s). When imposing the reference P = 0.1 gX/(L·d), the new steady-state is achieved faster, i.e., within 1 day, and the corresponding values are X = 0.39 gX/L, D = 0.2 d −1 and I = 150 µE/(m 2 ·s). This behavior can be explained by the fact that the RELS estimates need some time to converge when the first setpoint is imposed. The second setpoint can be achieved faster as the RELS estimation is probably close to the correct values and is able to track well the changes. We can only speculate about this observation as the values of the parameters are not known a priori.  Figure 11 shows the evolution of the gradient components with respect to the dilution rate and light intensity. The latter is to be minimized, that is, pushed to zero, as shown in Figure 11, where this objective is achieved within 1 day and maintained despite some variations due to the slower convergence of the dilution rate component to 0.2 (comparable to the biomass convergence rate, that is, 3 days). Following the step change in productivity reference, the gradient experiences a sudden change, but the signal converges fast back to zero while the dilution component of the gradient reaches the new value of 0.4. Regarding the application of the method in an industrial environment, a short practical guideline could be established, listing the steps to be achieved in the following order: • Run first the extremum seeking (tracking a zero gradient) to determine the maximum achievable productivity level; • Use the resulting productivity level as a reference to set achievable sub-optimal levels corresponding to non-zero slopes/gradients calculated by the controller; • Sequentially reset the controller, setting new sub-optimal productivity levels until the desired operating conditions are met.

Conclusions
In this paper, an adaptive slope seeking strategy that allows extremum seeking control and suboptimal setpoint tracking for dual-input single-output processes is proposed. In this setup, the process is represented by Hammerstein models combining a quadratic static block with a second-order time-invariant dynamic block. A recursive extended least square algorithm is used to estimate the static map parameters in the presence of noisy measurements. Based on those estimates, 2 algorithms are provided for the generation of slope references, with the possibility of minimizing one or both inputs. The slope-seeking strategy is applied successfully in simulation to the optimization of the lipid productivity in continuous cultures of micro-algae in photobioreactors. Further, it was validated at a laboratory scale in experimental tests aiming at biomass productivity optimization. Further work might consider the extension of this scheme to situations with more than two inputs, paying particular attention to the selection of appropriate excitation signals.