Model-Free Extremum Seeking Control of Bioprocesses: A Review with a Worked Example

: Uncertainty is a common feature of biological systems, and model-free extremum-seeking control has proved a relevant approach to avoid the typical problems related to model-based optimization, e.g., time- and resource-consuming derivation and identiﬁcation of dynamic models, and lack of robustness of optimal control. In this article, a review of the past and current trends in model-free extremum seeking is proposed with an emphasis on ﬁnding optimal operating conditions of bioprocesses. This review is illustrated with a simple simulation case study which allows a comparative evaluation of a few selected methods. Finally, some experimental case studies are discussed. As usual, practice lags behind theory, but recent developments conﬁrm the applicability of the approach at the laboratory scale and are encouraging a transfer to industrial scale.


Extremum Seeking: A Real-Time Output Feedback Optimization Technique
Real-time optimization (RTO) of steady-state plants stems from computer-aided production engineering (CAPE) and aims at improving process performance through the optimization of a measurable criterion or objective function under economical, safety or quality constraints [1]. One of the most critical limitations of earlier RTO methods lies in model accuracy, and this has lead to adapted RTO techniques [2] which take possible plant-model mismatches into account. The resulting model adaptations may impact the optimization sequentially as in run-to-run methods developed by [3], or also impact the optimization criterion and the related constraint gradients which are modified to match with the true plant [4,5]. Direct input adaptation methods transforming the RTO into a feedback control problem constitute a third alternative exploiting the invariance properties of the steady-state plant. These latter methods include: • Self-optimizing control (SOC), which tracks specific invariants or set-points of the considered process to achieve an indirect optimization [6]. • Extremum seeking control (ESC), which aims at driving the system to optimal operating conditions corresponding to the extremum of a measurable convex objective function [7][8][9]. • Necessary conditions of optimality (NCO) tracking, which uses the NCO as invariants to enforce optimality conditions [10][11][12] and gradient-based optimization. They are therefore considered the run-to-run optimization parent branches of both SOC and ESC; • Dynamic real-time optimization (DRTO) [1,13,14], which originates from the integration of RTO strategies, achieving process operating condition updates, in a model predictive control (MPC) framework. The recent additions of economic costs and constraints has lead to the introduction of DRTO solutions in the form of nonlinear MPC strategies with economic objectives [14]. Figure 1 shows the extremum seeking control (ESC) scheme as a simple output feedback where the objective function y is assumed to be measurable and the controller drives the system to an optimum y * of the steady-state function y = h(u).

The Origins of ESC
Before delving into the specific subject of ESC of bioprocesses, this section gives a brief historical perspective of the development of ESC in a more general context. For more details and references, we direct the reader to the excellent review of Tan [8] which covers 90 years of developments starting in the 1920s, and the recent research monograph of Scheinker and Krstic [15].
The idea of using a perturbation signal at the system input and observing its effect at the output to estimate the slope of a nonlinear static map can be traced back to the work of the French engineer Leblanc in 1922 [16] in the context of power transfer from a transmission line to a tram car. However, it seems that it was only during World War II that the technique was investigated again by Russian scientists [17]. In 1951, Draper was the first to report research on the subject in the English literature [18]. The 50s and 60s were then flourishing with various works, such as [19,20], essentially contributing to the algorithm's performance without truly developing a systematic design framework. In this period, Blackman [21] published an interesting survey of those methods.
Basically, classical ESC methods can be subdivided into methods that involve an external perturbation (dither) signal and those which do not require one and resort to relay switches [18,22] or sliding mode techniques [23].
In 1971, a first Lyapunov stability analysis of a so-called peak-holding scheme, in which the gradient of the performance index is not estimated, but the variation of a parameter allows sensing the direction in which the performance index is evolving, was developed by Luxat and Lees [22]. Three decades were needed for a general stability analysis by Wang and Krstic [24] following the impulse of [25] in developing the framework of ES.
An important milestone in fostering the popularity of ES is the publication in 2003 of the research monograph [7] by Ariyur and Krstic. From an (real-life) application point of view, ES was applied successfully in various fields including maximum power point tracking in photovoltaic systems [26,27], drag reduction in formation flight [28], continuously variable transmission [29], ignition systems in combustion engines [30], vapor compression systems [31,32], wind energy conversion systems [33], cavity resonators [34,35], braking sytems in automotive industry [36].
The classical perturbation-based ES scheme suffers from the time-scale separation between the process which is considered as "fast" (as the process has to settle down quickly and be viewed almost as a static map by the controller), and the perturbation signal which must be slow with respect to the process dynamics. Successful applications mostly consider relatively fast electrical and mechanical systems where the slow convergence of ES is not detrimental. Applications with large time constants and long transient phases are more problematic, as is the case for bioprocesses which are the subject of the next sections. One way around the slow convergence is to include some prior knowledge about the model, leading to model-based ES, or some form of parameter estimation, which allows one to keep on model-free extremum seeking. This distinction is also addressed in the following.

ESC of Bioprocesses
The classical bioprocess optimization problem [37,38] considers a culture of microorganisms or cells in a continuous bioreactor where the following generic reaction takes place where X, S, and P denote, respectively, the biomass (microorganism/cell population), the limiting substrate, and the product of interest. k 1 and k 2 are stoichiometric coefficients. The reaction kinetics can be described by various laws, among which the most popular are Monod law (2a) describing substrate activation; Contois law (2b) describing substrate activation and biomass inihibition; and Haldane law (2c), describing substrate activation at low concentrations, and inhibition at higher concentrations: Mass balances lead to a set of ordinary differential equations: The measurable (via direct measurements or state estimation) performance index can be linked to the productivity of the biomass, the product of interest, or the reaction rate: Optimization of bioprocesses was first addressed in 1999 by [37] in a numerical study using a classical ES scheme with a sinusoidal dither signal and a combination of high-pass and low-pass filters (further developed in [7]). In the same period, the experimental work of Akesson [39,40] applied a simple and practical probing method in the spirit of the peak-holding algorithm of Luxat and Lees [22], the idea being to send pulses of substrate and sense whether it is consumed by the micro-organisms. This work can be considered as the seed of a more general ES theory using nonlinear programming and developed in [41].

Extremum Seeking Approaches
Depending on the available knowledge about the model structure, different extremum seeking approaches have been developed. Logically, the better the system is known, the more efficient the controller may be. Roughly, the ES controllers can be classified into two main streams: • Model-based strategies, wherein a model-based controller takes advantage of the knowledge about the model structure and either estimates the unknown parameters from online information, or exploits some robustness properties. • Model-free strategies, wherein online information is directly exploited to estimate the gradient of the measurable cost function, and to drive the process to the optimum using feedback control.
No prior knowledge about the kinetic model is required, and the optimum is simply reached by driving the process wherein the estimated cost function gradient vanishes.
Of course, the two streams are not strictly disjunct, and for instance, one can imagine taking advantage of any prior knowledge to improve a model-free strategy. While this paper mostly focuses on the second category of methods, the first stream has brought a number of important results that are briefly reviewed in the next section.

Model-Based Strategies: Uncertain Trajectory Tracking
Model-based ES strategies rely on a hypothetical kinetic structure and aim at driving the system to an optimal operating point or along an optimal operating trajectory. Since kinetic parameters are unknown or at least uncertain, the control law requires their estimation by adaptive laws or the consideration of model uncertainties and robust design.

Adaptive Output Feedback Control
The concept of adaptive ESC of nonlinear systems with uncertain parameters has been developed by Guay and Zhang [42], who introduced an "inverse optimal" approach to solve the optimization problem. The application of this method to continuous bioreactors has been proposed in the wake of this latter publication in [43] starting from the simple microbial growth model (3a), (3b), and (4c). These equations can be rewritten in the following compact form, replacing the kinetic term by its measurement:ẋ Define the dilution rate as the system input (u = D); then an extremum seeking control law is proposed to drive and stabilize the system in the neighborhood of y * = f (µ * ). Analyze the behavior of the output time derivativeẏ = h(x, s, u, θ) where θ contains the stoichiometric parameters k 1 and k 2 but also the kinetic parameters involved in the hypothetical structure of µ. Since these parameters are unknown, parametric adaptive laws are inferred from the overall loop, including the extremum search combined with the adaptive controller. Stability and convergence studies resort to Lyapunov stability theory coupled to LaSalle-Yoshizawa's theorem, LaSalle's invariance principle, or Barbalat's Lemma [44][45][46]. Further application results can be found in [38,47], which propose model-based extremum seeking controllers optimizing either the reaction rates or the biomass production of continuous and fed-batch processes, characterized by substrate activation and saturation following Monod kinetics (2a). Further, Refs [48][49][50][51] have extended these studies to Haldane kinetics (2c), which represents the additional effect of substrate inhibition at high concentrations. More complex kinetic structures have been considered in [52] where biomass productivity is optimized under the assumption of microbial overflow metabolism [53]. In this case, growth inhibition results from the accumulation of by-products reducing the cell respiratory capacity, and the kinetics involve switching between different models. As discussed in the survey [47], adaptive extremum seeking can be generalized to the situations where no prior knowledge is assumed, using universal approximators such as neural networks [54].

Sliding-Mode Control
Another approach consists of optimizing the cost function by sliding-mode control along the cost function manifold and choosing its gradient as switching function σ = ∂y ∂u . Based on this idea, Lara-Cisneros et al. [55] proposed the optimization of fed-batch bioreactors modeled by (3) using a sliding-mode extremum seeking control combined with a high-gain observer estimating the uncertain kinetics. This study was extended to a two-stage anaerobic digestion model in [56]. In the same spirit, Vargas et al. [57] tackle the problem of biomass productivity optimization in bacteria cultures, presenting overflow metabolism using a two-level controller with super-twisting observers as process output estimators. Further, an input-output linearizing sliding-mode control is derived for well-defined relative degree systems with minimum phase in [58,59].

Integration of ES to Model Predictive Control
The optimization problems introduced in Sections 2.1.1 and 2.1.2 do not consider constraints, which can be unrealistic in the context of real-life applications. Indeed, industrial processes are run under specific conditions implying media compositions, mininum and maximum authorized production rates, and energy consumption costs. Application of ES control to constrained optimization problems of control-affine nonlinear systems such as (3) was first proposed in [60] using a flatness-based approach. The optimal trajectory is generated by adaptive ES, and a state and input-constrained dynamic nonlinear programming problem is solved using an interior-point technique. This approach is applied to batch processes in [61].
The study of [13] generalizes the approach by integrating ES to model predictive control (MPC) considering a time-varying economic cost function based on a model presenting parametric uncertainty sets and state constraints. Uncertainty issues are solved by considering min-max optimization to determine the worst/best solutions with respect to the parameter uncertainty sets and a Lipschitz-based adaptive approach imposing Lipschitz bounds. This method is extended to economic model predictive control (EMPC) in [62], which does not track a specific trajectory but directly solves the time-varying economic optimization problem.

Model-Free Extremum Seeking Control
Model-free ES strategies rely on the existence of a measurable and convex objective function but do not require prior knowledge about the process model. Their principle is to estimate the gradient of the objective function and to force this estimate to zero (a vanishing gradient being the first necessary condition for an extremum of a differentiable function).

The Classical Pertubation-Based ES and the Use of a Filter Bank
One of the earlier forms of model-free ES is the perturbation-based ES of [7] shown in Figure 2 and described by the following equations : whereû is the estimated input; u =û + d stands for the input signal applied to the process, including a sinusoidal dither signal d = A sin(ωt); y is the measurable objective function; y − η is the high-pass filtered output, that is, without a continuous component; δ is the demodulated signal resulting from the multiplication of y − η by d; andξ ≈ 1 k I ∂û ∂t is the gradient estimation, obtained as the output of an optional low-pass filter. The use of a bank of filters leads us to designate this scheme as "BOF" (a wink to the French speaking reader, but there is no prejudice in this denomination) in the sequel of this review (see Figure 2). The design parameters of this ESC have to obey some rules: the cut-off pulsations of the low-pass filter ω l and high-pass filter ω h should be smaller than the dither signal pulsation ω; i.e., they should be located in the interval ]0 ω[. On the other hand, the dither signal pulsation ω should be small with respect to the process dynamics-for the process to be in periodic steady state-and its magnitude sufficiently large to overcome the information attenuation by the filters.

Static Map
These rules reflect the time-scale separation [7] already highlighted in the introduction to this review, where: • The process has the fastest dynamics, which, in the limit, can be seen as a static map; • The dither signal has intermediate dynamics;

•
The gradient estimator has the slowest dynamics.
It is important to note that the issue of the time-scale separation worsens in cases of systems with multiple inputs [7]. In the particular case of multi-input single-output systems with r inputs, it is necessary to design r+2 2 (rounded to the next lower integer) dither frequencies [63], thereby requiring more time separations and slowing down the gradient estimator.
A suitable choice of the design parameters ensures an exponential convergence of y in a O(ω + A) neighborhood of the optimum y * = f (u * ) [64], where u * is the corresponding optimum input.
Optimizing control of bioreactors using BOF-ES was first reported in [37], where a bifurcation analysis of model (3) with Monod and Haldane kinetics was achieved. This study highlights the operational limits and suggests the use of wash-out filters to guarantee the system's stability. Subsequently, Nguang and Chen [65] proposed an illustrative application of the same strategy to continuous fermentation processes and [66] to anaerobic digestion.
However, perturbation-based techniques and the BOF-ESC are only successful if the system can be seen as almost static, that is, slowly perturbed by the dither signal. The slow convergence resulting from this time scale separation can be a serious drawback in biological applications where the dominant time constants can be very large. Moreover, the characterization of the effects of the dither signal frequency and magnitude on the convergence, as presented in [64], is a significant additional issue for general nonlinear systems. This has motivated further studies, such as [67], which concluded that the convergence accuracy of BOF-ES is only concerned with the magnitude of the dither signal for a specific class of nonlinear systems represented by Wiener and Hammerstein block-oriented models. Furthermore, Tan et al. [68] considered several dither signals and discussed their impacts on the domains of attraction, convergence, speed, and accuracy.
To alleviate the undesirable time scale separation effect, Krstic [69] suggested introducing controllers to speed-up the closed-loop dynamics or to compensate for the phase-shift introduced by the system dynamics at the perturbation frequency. This latter idea has been proven achievable by [70], depending on the system relative degree, and thus, for a restricted class of nonlinear systems.
In the last decade, several practical design guidelines have been proposed to improve the convergence rate, without true analysis, based on the selection of the operating conditions [71], or the use of several dither signals to enrich the frequency content [72,73]. In [74], an application of BOF-ESC to yeast fed-batch cultures with overflow metabolism is presented, where the convergence speed was enhanced thanks to an appropriate manipulation of the substrates presenting fast dynamics.
The first experimental validation of BOF-ES in the field of bioprocess optimization can be attributed to [75] optimizing microalgae growth at pre-industrial scale by manipulating the pH. Overall, BOF-ES can definitely be considered as a useful, practical, and easy-to-implement solution to optimize bioprocesses as long as operating conditions, dithers and controlled variables are carefully selected with respect to system dynamics.

Recursive Parameter Estimators
A simple solution to improve the convergence speed of perturbation-based ES algorithms is to replace the bank of filters, that is, one part of the time scale separation, in a recursive estimator as shown in Figure 3, where: This ES structure convergence is now only dependent on the system dynamics and the dither signal parameter selection. An illustration of this recursive strategy, using a continuous recursive least-squares (RLS) algorithm with forgetting factor [25,76], can be found in [74], which discusses the application of BOF-ES and RLS-ES to biomass productivity optimization in fed-batch cultures. The main conclusion of this study is that the recursive estimator provides a more robust behavior and a faster convergence. The continuous ES-RLS scheme is described as follows: where K is the positive adaptation gain; R is the inverse covariance matrix acting as a directional adaptation gain; λ is the forgetting factor used in order to avoid covariance wind-up issues due to the absence of bounds in R growth (if λ = 0,Ṙ ≥ 0 [76]). Later, Chioua et al. [77] characterized the convergence rate and proved the stability of the ES loop with RLS as the gradient estimator for nonlinear systems, independently of their possible decompositions in block-oriented models (Wiener/Hammerstein).
In the same time period, Guay and Dochain [78] have reformulated the ES description as a time-varying control problem for general nonlinear systems, using the combination of a continuous parameter estimator based on a generic prediction model of the output y (not requiring a priori knowledge on the process), and an original integral control law. The main improvement of the method lies in the use of the dither signal, which, following a favorable contribution from the natural time-variation of the unknown parameter (gradient), becomes no longer necessary. Guay and Dochain [78] therefore show that the suppression of the intermediate dither dynamics is also possible, limiting the convergence speed to the process dynamics and the estimator performance.

Empirical ESC: Probing Control
The model-free methods presented in Sections 2.2.1 and 2.2.2 rely on the existence of a convex and measurable objective function. When the available process instrumentation does not allow this latter condition to be verified, one may resort to more practical or even empirical software solutions. For instance, optimizing biomass productivity of cells exhibiting overflow metabolism requires one to measure an objective function related to the cell respiratory capacity which should be maximal in order to optimize substrate consumption [53]. This critical measurement can, however, only be achieved indirectly [74] and the convex objective function is reconstructed using measurement signals and some prior knowledge about the reaction stoichiometry. The performance of the ES technique therefore depends on the stoichiometric coefficient estimation accuracy, which may become critical. In the early 2000s, Akesson et al. [39] proposed an interesting empirical probing control method detecting the possible saturation of the respiratory capacity of Escherichia coli, and in turn, inhibitory by-product (acetate) accumulation [39], by measuring the oxygen uptake following a sudden change in the substrate (glucose) uptake caused by feed pulses. Indeed, when the cell's respiratory capacity is saturated, the dissolved oxygen concentration (measured variable) does not vary any longer, since cells are not able to assimilate more oxygen. In analogy with classical ESC, probing control uses these pulse sequences as dither whose magnitude and frequency are designed in accordance with the sensitivity level of dissolved oxygen variation detections, and the combined respiratory capacity and substrate uptake maximization define the extremum of an unknown objective function. Probing control became popular during the first decade of the 21st century, promoted by several experimental studies in various fields such as human health [79] and agriculture [80].
The characterization of overflow metabolism for yeasts (S. cerevisiae, P. pastoris) and bacteria (E. coli) was also provided by [81] using probing control and a scaling-up approach of the method was proposed by [82].

Convergence Issues and Acceleration
The ESC problem, as defined in the previous sections, assumes basic practical constraints such as fast process dynamics (with respect to the dither and the optional filters) and a single measurable objective function. However, several issues may arise when the process natural dynamics are detrimental to objective-related time restrictions or when several competing objectives must be handled. In the following, some practical ways to manage and accelerate the ESC loop convergence are addressed, starting with a simple feedback form.

Extremum Seeking Control: A Simple Output Feedback Form
A minimal ESC loop description, as illustrated in Figure 4, considers two blocks: the gradient estimator and the integral controller pushing this gradient to a reference value (in the case of extremum seeking, zero) ξ re f such that the error = ξ re f −ξ is zero. A simple pole-placement procedure can be used to define the ES closed-loop dynamics, i.e., to design the integral constant k I . As discussed in [83], a controllable state-space representation of the gradient evolution can be inferred as: where w is a function of the input reference u re f and the parameter θ estimated by the recursive estimator from Figure 3, both assumed to provide the desired gradient ξ re f . By comparing Figure 4 with Figure 3, we can easily infer that A = 0, B = 1, and C = 1. By introducing an integral action, the system reads: When applying classical pole-placement procedures, the closed-loop dynamics can be specified in the time domain, yielding the integral (extremum seeking) control gain. However, it must be noticed that some unexpected transient dynamics may not be taken into account and u re f must therefore be chosen in accordance with the domain of attraction and stability analysis, often restricting the input operating range to a specific interval [u min u max ]. Consequently, if the closed-loop response presents transients due to neglected system dynamics, such that the input is out of this permissible range, the overall system may fall into unstable regions. To minimize the risks of loop destabilization, a sufficient level of robustness should be imposed, with large gain margins and limited overshoots. Figure 4. Extremum seeking as a simple output feedback.

Proportional-Integral Extremum Seeking
In an alternative attempt to accelerate the convergence of the ESC, Guay and Dochain [84] first proposed to combine the time-varying estimation method from [78] with a proportional-integral controller, as shown in Figure 5. This method generalizes standard ESC using the integral action as "extremum seeker" and proportional action as the fast optimizer of the objective function, amplifying the estimated gradient descent or ascent effect. A general formalism of proportional-integral extremum seeking control scheme (PIESC) has been provided by the same authors in [85]. The closed-loop PIES controller reads: where k p and k i are respectively the proportional and integral gains. A second PIESC approach has also been proposed by [86], using the classical perturbation-based BOF scheme. In this study, the parameter (gradient) estimator block in Figure 5 is replaced by an output derivativeẏ estimator modeled by a high-gain high-pass filter block with transfer function ω h s s+ω h and ω h >>. The filter output may then be read as: The tuning rules proposed by [86] may be summarized as follows: • The PIESC with BOF converges to a 1 ω neighborhood of the optimum and ω should therefore be chosen sufficiently high; • ν(0) should be chosen equal to y ω to avoid sudden initial jumps arising fromν; • ω h should be of a greater order of magnitude than ω to favor the derivative effect.
An application of the PIESC to the optimization problem of compressor discharge temperature setpoint selection for a vapor compressor system, minimizing power consumption, has recently been achieved by [87], who showed the performance of the method in realistic operating conditions (noise, input/sensor quantization), through experimental validation on a room air conditioner. As a result, the PIESC is able to achieve the optimizing control in a time-scale in the order of the process dynamics. Since no PIESC design for bioprocesses has been made available in the literature, a first attempt is provided in the current review in Section 4.1 and numerical results are in Section 4.4.

Dynamic Modeling with Block-Oriented Representations
As stated in Sections 2.2.1 and 2.2.2, the convergence accuracy of BOF-ES and RLS-ES applied to a specific class of nonlinear systems represented by Wiener/Hammerstein (W/H) block-oriented models, is only impacted by the dither signal magnitude. Semi-global stability results of the application of observer-based extremum seeking to W/H plants is presented as an extension of the work of [88] in [89,90], showing that the ES scheme requires the knowledge of the relative orders of both input and output dynamics of the plant. From the latter results, a discrete version of the method for W/H plants with knowledge limited to the linear dynamics is also proposed, suggesting guidelines for the choice of the ES parameters, first setting the dither signal magnitude and then adapting the ES integral gain. The latter also presents an experimental validation of the method on a spark-ignition engine. Eventually, the same authors proposed directly adding some model knowledge in the adaptation scheme of the fast extremum seeking method in [91]. Interestingly, this allows improving convergence not by completely removing the time-scale separations of, in decreasing speed order, the optimizer and the estimator, but by reversing them, that is, by accelerating the estimation and subsequently the optimization. Another fast algorithm using proportional-integral phasor extremum seeking was proposed by [92] for W/H systems; its semi-global and practically asymptotical stability was shown. The principle is to drive the system to its optimum by considering the measured output phasor, from which a proportional approximation of the static cost function gradient is obtained, under the assumption of a priori known input and output dynamic block relative orders. It should be noticed that considering the phasor allows taking both the gradient and the phase shift, respectively related to the static map and the plant dynamics, into account. This estimator is then combined to a PI-ES scheme, in the spirit of [86] with an original amplitude update mechanism reducing the perturbation amplitude as the algorithm approaches the optimum vicinity, in the spirit of [93], which will be discussed in the next subsection.
In the field of biosystems, Feudjio et al. [83] proposed an advanced auxiliary model-recursive prediction error method (AM-RPEM)-based ES algorithm for systems represented by block-oriented models. In this algorithm's structure, represented in Figure 6, a slope generator is coupled to the controller to reach any steady-state belonging to the static map, making ES a particular case of the slope-seeking strategy. The considered class of bioprocesses is characterized by the following assumptions: • The bioprocess is operated in continuous mode with a constant volume (the inflow is equal to the outflow); • The measurable performance index is the biomass production y = DX, i.e., the product of the dilution rate D imposed by the inflow and outflow pumps by the biomass concentration X; • The bioprocess can be approximated by a Hammerstein model [94] presenting a quadratic static nonlinearity of the form: v = M 1 u 2 + M 2 u + M 3 (13) where u = D and v represent the static map input and output signals, respectively. This representation assumes that the performance index can be approximated, at least locally, by a quadratic form. For a maximum, it is required that M 1 < 0. The Hammerstein model also includes dynamics, which can be described by a transfer function. • A discrete second-order transfer function may suffice to represent these dynamics: where z −1 is the backward shift operator such that z −1 y(k) = y(k − 1) and K is a gain which can be chosen so as to ensure unitary steady-state gain, i.e., G(1) = 1, which implies K = 1+β 1 +β 2 1−α . The two poles defined by the parameters β 1 and β 2 represent the actuator and process dynamics.
Under these assumptions, the unknown model parameters can be estimated using a simple recursive least-squares (RLS) method, as originally proposed in [74], or in order to improve the convergence and precision of the estimates, the recursive prediction error method (RPEM) from [95,96], where, unlike classical RLS algorithms which provide as estimates the linear combinations of the unknown parameters, the Hammerstein model parameters are directly inferred. However, since analytically computing parametric sensitivities can be tedious, an auxiliary model is used to simplify the evaluation of the sensitivities. Another option is also provided, considering an output-error auto-regressive model [97] which is mostly useful in the presence of relatively large process/measurement noise.

Multivalued Cost Functions and Competing Objectives
An intrinsic property of the previously discussed systems is the existence of a sole extremum. For complex (bio)processes, it is legitimate to consider possible multiple local extrema of the steady-state objective function. In [93], this problem is tackled using an adaptive dither magnitude update law, specifying the conditions of convergence to the global extremum. The proposed solution consists of designing and initializing the dither signal magnitude law to allow the bifurcation towards the submap of the objective function where the global extremum is located. Since ES convergence is affected by the dither signal magnitude, the latter is adaptively decreased with time, as the algorithm approaches the global extremum. Starting from a simplified ES scheme drawn in Figure 7, the closed-loop representation can be written, denoting σ = ωt: where f stands for the nonlinear state space function of the system, x is the state vector, δ is a positive design parameter, A(t) is the varying magnitude of the dither signal with A(0) > 0, and g(A) is a locally Lipschitz function that is zero at zero and positive otherwise. Fast and slow dynamics are inferred by considering x in equilibrium, that is, x eq = l(û + Asin(σ)), which is assumed to be asymptotically stable, leading to the following reduced system: Considering the corresponding average system over one 2π period and introducing τ = δσ, (16) becomes: Tan et al. [93] have shown that there exists a continuous function p(u, A) such that Q av = A p(u, A). The standard singular perturbation technique is unfortunately not applicable (for more details about this statement, see [93]) and it is therefore recommended to analyze the convergence properties of the ES scheme for a specific initialization A(0) = A 0 , assuming that there exists an isolated and unique real solution u = (A) to Q av = 0 such that is continuous and ∂p ∂ < 0 for all A ≥ 0 and u * = (0) is the global optimum. Following this latter assumption, it may be demonstrated that the convergence of the extremum seeking scheme of Figure 7 to u * is practically semi-global. The global extremum of a multimodal objective function can therefore be reached by the perturbation-based ES following judicious choices of the initial dither signal magnitude and the function g(A). A practical interpretation of the latter results in optimal operating mode seeking for bioprocesses is discussed in [98], where a multivalued cost function originating from the combination of both yield and production objective functions is derived from model (3) where only the biomass and the substrate are considered. The kinetics are modeled by a simple Monod law and an exponential factor representing growth inhibition by the substrate (which, physically, can be compared to the Haldane kinetics (2c)). Consequently, the resulting cost functions present two distinct quasi-steady-state maps between which the ES algorithm must be able to operate a switch when initialized on the map which does not contain the global extremum. Considering only the stable branch, the set-valued cost function Q(u, A) is, in this particular case, defined as a set of two continuous single-valued functions such that, omitting u and A for the sake of clarity, Q = {Q 1 , Q 2 } with the conditions: For each value of u 1 and u 2 , there are asymptotically stable equilibria x = 1 (u 1 ) and x = 2 (u 2 ); To enable some switch conditions from one single-valued cost to the other, in addition to the application of the results from [93], it is concluded that restrictive conditions involving the initial value A 0 and g(u, A) should be used in order to make the strategy successful, also implying a loss of performance.
The existence of several extrema may also result from a dynamic separation interpretation. Even if a process presents a single and a priori global optimum on the steady-state map, the isolated fast dynamic analysis may reveal the existence of a transient extremum. An illustration based on the CANON process, used to remove ammonium from concentrated wastewater streams, is provided by [99] where the fast dynamics, called the boundary layer model, are represented by the substrate and flow dynamics, while biomass growth is assumed to be the slow dynamics. In this example, even when the biomass is not at steady-state, there exists a transient optimum on the boundary layer, and more interestingly, it turns out that locally optimizing the transient objective leads to the global optimization of the steady-state objective function, which comes to optimize an instantaneous objective without taking the future trajectory into account. A phasor-based extremum seeking controller inspired by [100] is therefore applied and the resulting algorithm is called "greedy" ES in analogy to greedy algorithms in computer science [101] that select the best instantaneous options without regard to future decisions.

An Illustrative Example: Optimization of Biomass Productivity
The aim of this section is to provide a worked illustrative example of the respective performances of the model-free ES algorithms introduced in Sections 2.2 and 3 where biomass production (4a) is optimized using model (3) and considering three classical kinetic structures: Monod (2a), Haldane (2c), and Contois (2b) laws. All the parameter values used in the following simulations are included in Table 1. Regarding the bifurcation analysis of system (3) with Monod and Haldane kinetics, the reader may refer to [37] for detailed results which can be summarized as follows (the wash-out equilibrium (x e , s e ) = (0, S in ) is not considered): under the assumption of Monod kinetics, wherein the "e" index denotes equilibrium values and "*" index values at optimum, and under the assumption of Haldane kinetics for which there exist two equilibrium trajectories. A similar analysis using the Contois kinetic law (2b) can be achieved and provides the following original results. The new equilibrium reads: and the equilibrium definition condition is D e < µ max . The steady-state extremum y * e can be calculated considering ∂y e ∂D e = 0, which leads to the following optimal dilution rate: The stability analysis of the equilibrium requires one to determine the Jacobian of system (3) while considering the Contois kinetic law, which yields: and in steady-state, we obtain: The eigenvalues of (23) are: and the stability conditions on D e are then:

Bank of Filters with PI Control
A PIESC with BOF tuning is suggested, based on [86], as follows: where d = Asin(ωt) is the dither signal and the output derivativeẏ is computed by an additional high-pass filter (12), of the form:η where η is an intermediate variable, the image of the output derivative, and ω h is the filter cut-off frequency which should be chosen to be greater than ω, the dither signal angular frequency. Even if some tuning rules are proposed in [86] to design ω and ω h , there is unfortunately no indication concerning the tuning of k p , k i , and the dither signal magnitude A which must be designed by trial and error. In order to illustrate the advantages of the PI method with respect to the classical BOF-ES from [7], two parametrizations are considered in the following and are summarized in Table 2. The first case corresponds to the simple BOF-ES proportional control (k p = 0) strategy, for which the parameters are optimally selected by trial and error to provide the best transient behavior (fastest convergence).
For an accurate analysis of the choice of the dither signal magnitude and frequency, the reader may refer to [88] for semi-global stability results and [67] for the particular cases of Wiener/Hammerstein approximations. The second case corresponds to the BOF-ES PI control application with a specific parametrization respecting the tuning rules of [78].

Recursive Least-Square Estimation with PI Control
A simple recursive least-square extremum seeking algorithm (RLS-ES) combined with the PI control structure introduced in Section 4.1 is proposed in the same spirit as [85]. The PIRLS-ES controller slightly differs from (26) which is rewritten: whereγ replacesẏ and denotes the estimator of the cost function gradient y u . Indeed, taking the partial derivative of y = ux (u = D) with respect to u yields: Multiplying (29) by u, we obtain: where ϕ = [u 1], θ = [θ 1 θ 2 ] = γ − u 2 ∂x ∂u , and γ = ∂y ∂u . It must be noticed that when reaching the extremum steady state y * , γ = 0, and from (29), y * = θ * 2 = −u * 2 ∂x ∂u u=u * ,x=x * . This time, the cost function gradient is assumed to be estimated by a continuous RLS with forgetting factor [76] which reads: The chosen parametrization, which is also set by trial and error, is reported in Table 3. Table 3. PIRLS-ES controller parameters.

Hammerstein Model and Pole-Placement
The recursive least-square method of Section 4.2 is now considered under the assumption of a block-oriented model representation as in [83]. The microbial growth model (3) is approximated by a Hammerstein model with a static map x(t) = F(u(t)) followed by a first-order (strictly proper and stable) transfer function describing the system/sensor dynamics G(s) = 1 1+τs . Since measurements are collected at discrete times, a discrete form of G(s) is considered using the matched pole-zero method, leading to G(z) = K 1 z−α with K 1 = 1 − α and α = e − Ts τ , T s being the sampling period. The corresponding scheme will be denoted BOM-ES (block oriented model extremum seeking) in the following. The static map is approximated by the linear form : where x is an intermediate state variable, u = D the input, ξ the gradient, and c a constant. The first-order dynamics reads: Combining (31) and (32) results in: The regression (33) and associated RLS problem provide the gradient estimateξ = θ 3 1−θ 2 of ξ wherê ξ ≈ ∂y ∂u , the gradient estimation, is driven to zero in average by the extremum seeking control loop. A dither signal d = Asin(ωt) is used to ensure a persistency of excitation [25]. As a general rule, a minimum of n 2 distinct sinusoids is necessary for the identification of n parameters [63]. As suggested in Section 3.1, a simple pole-placement procedure is applied. The BOM-ES parameters are summarized in Table 4. A good compromise between convergence speed and robustness is to impose a rise time of 20 h to the closed-loop system, resulting in an integral gain k I = 5 × 10 −4 . Since 3 parameters of (33) are estimated, the dither signal is chosen as the sum of two sinusoids (which allows identifying up to 4 parameters) d = A 1 sin(ω 1 t) + A 2 sin(ω 2 t) and the measurement sampling time is set to T s = 0.1 d.

Numerical Results
The several ES algorithms designed in Sections 4.1-4.3 are now compared when optimizing the biomass production, y = Dx in the presence of different kinetic structures, such as Monod, Haldane, and Contois. First, Figures 8 and 9 show the results of the application of the BOF, PI-BOF, PI-RLS, and BOM extremum seeking methods to the Monod case.
Not surprisingly, BOM and PI-BOF algorithms achieved the fastest convergences of, respectively, 50 and 25 h; PI-RLS and BOF converged within, respectively, 300 and more than 1000 h. However, Figures 8 and 9 also show that PI-BOF reaches the optimum neighborhood with a significant offset (see the zoom box of Figure 9) while the convergence of BOM, even if slightly slower, is more accurate. These first results are confirmed in the Haldane case in Figures 10 and 11, even if the BOM algorithm performs better this time in terms of convergence time and accuracy, and PI-RLS and BOF respectively converged within 200 and 600 h.
Nevertheless, the algorithm behaviors are interestingly emphasized when using the Contois kinetic law, as shown in Figure 12. While the performances of the BOM and BOF strategies are identical to the previous cases, PI-BOF and PI-RLS show larger offsets with respect to optimal state and dilution rate equilibria. This can be attributed to the gradient values which are, however, still very close to 0. In the particular case of Contois kinetics, the convergence diagram of Figure 13 shows that the equilibrium slopes are much smoother compared to the previous Monod and Haldane cases. A small shift in a close neighborhood of the optimum, keeping y ≈ y * , may then lead to important offsets on the dilution rate and therefore on the states.
Observations based on these several numerical illustrations can be summarized as follows: • The BOM strategy is by far the most accurate since it allows taking unknown output dynamics into account. • The PI-BOF strategy allows solving the extremum seeking problem almost instantaneously and therefore provides the fastest convergence, but as a drawback, has a technically inherent and significant offset (see [78]). • The PI-RLS strategy presents a smaller offset and is overall easy to implement while achieving reasonable convergence times. • The BOF strategy, even if in the original it was the most intuitive and easiest to implement method, presents several drawbacks, such as long convergence time, tedious parameter tuning, and possible offsets.

Real-Life Applications: Microalgae Production Optimization
Most of the publications discussed in this review propose simulation results in order to provide some practical assessment and recommendations in view of future experimental implementation. Among the very few experimental validations reported in the literature, two investigations aiming at maximizing microalgae culture productivity by extremum seeking have been proposed in [75,102].

Microalgae Growth Rate Optimization
In [75], the method consists of modulating the pH to optimize the measured Nannochloropsis Oculata microalgae biomass growth rate, using the CO 2 "flow as actuator" in an internal loop regulating the pH. In this regard, a perturbation-based filter bank ES was first designed in a simulation and some tuning guidelines were suggested as follows: • Quantification of the system dynamics and uncertainties; • Selection of the dither frequency, and to some extent, its structure; • Trial and error optimization routine defining the high-pass filter cut-off frequency; • If the level of noise requires it, determination of the low-pass filter cut-off frequency.
The use of a low-pass filter should, however, be avoided as much as possible since it is a source of convergence slow-down (the reader may refer to the time-scale separations explained in Section 2.2.1). Following the initial simulation results, an experimental validation with a 27-liter tubular bubble column photobioreactor reproducing a down-scaled industrial process was achieved. The main results show that the ES scheme managed to bring the system to an acceptable pH range, as predicted by other works from the literature.

Microalgae Productivity Optimization
In a recent experimental investigation achieved by the authors [102], a RLS-ES controller was designed, assuming a Hammerstein representation between the biomass X of the strain Dunaliella tertiolecta, as measured output, and the dilution rate D, assumed to be the sole input. The biomass productivity, defined as P = DX is then maximized and its gradient estimate is provided by: assuming the static mapping X(k) = m D(k) + b. It must be noticed that a filtered biomass variable X f may also be considered leading to similar performance. The estimated optimal productivity and dilution rate can be inferred as: The experimental set-up was a 13-liter flat photobioreactor shown in Figure 14. For more details regarding the monitoring devices, the reader may refer to [102]. A ten-day experiment was performed, leading to the results presented in Figures 15-17, wherein the gradient was indeed pushed to zero on average within 4 days leading to the maximum productivity level. Figure 17 compares the predicted static map with the convergence results, proving the method's accuracy.    Following these encouraging first results, further experiments were scheduled, using the BOM-ES method from Section 4.3, and are currently in validation.

Conclusions
During the last two decades, extremum seeking (ES) has become a very popular research subject, and has lead to an increasing number of real-life applications to electro-mechanical systems. The application of ES to biosystems is, however, quite tedious due to the uncertain nature of the "living" in a broad sense, and the slow dynamics of the bioprocesses. In view of the model uncertainties, model-free ES appears as an appealing approach to achieve "plug-and-play" process optimization. Even if the seminal perturbation-based techniques based on the use of high and low pass filtering provide guaranteed results if a correct parameterization is selected, the three time-scale separations entail very lengthy convergence times. Most of the recent ES approaches therefore focus on ES convergence issues, acting on the controller structure and design, on parameter estimation, or on the dither signal definition. This review paper discussed some of these techniques using an illustrative simulation example, where the biomass production of a simple microbial growth model was maximized considering several kinetic structures with different complexity levels and corresponding static maps.
Only a few real-life applications of ES to bioprocesses are described in the open literature, namely, those to continuous microalgae cultures in photo-bioreactors. Hopefully, these experimental studies will pave the way for future studies at pilot and industrial scales.

Conflicts of Interest:
The authors declare no conflict of interest.