Adaptive Control of a Fermentation Bioprocess for Lactic Acid Production

Abstract

This paper presents the design and the analysis of an indirect adaptive control strategy for a lactic acid production, which is carried out in continuous stirred tank bioreactors. Firstly, an indirect adaptive control structure based on the nonlinear process model is derived by combining a linearizing control law with a new parameter estimator. This estimator is used for on-line estimation of the bioprocess unknown kinetics, avoiding the introduction of a state observer. Secondly, a tuning procedure of estimator design parameters is achieved by stability analysis of the control scheme. The effectiveness and performance of estimation and control algorithms are illustrated by numerical simulations applied in the case of a lactic fermentation bioprocess for which kinetic dynamics are strongly nonlinear, time varying, and completely unknown, and not all the state variables are measurable.

1. Introduction

The control of biotechnological processes has been and remains an important problem attracting wide attention, the main engineering motivation being the improvement of the operational stability and production efficiency of such living processes.

It is well known that traditional control design involves a complicated mathematical analysis and has many difficulties especially in controlling highly nonlinear and time varying plants. A powerful tool for nonlinear controller design is the feedback linearization technique [1–3], but the use of it requires the complete knowledge of the process. In practice, there are many processes characterized by highly nonlinear dynamics; as a consequence, an accurate model for these processes is difficult to develop. Therefore, in recent years, it has been noticed a great progress in adaptive and robust adaptive control of nonlinear systems, due to their ability to compensate for parametric uncertainties. An important assumption in previous works on nonlinear adaptive control was the linear dependence on the unknown parameters [2, 4, 5].

In the modern industry, the development of advanced control strategies for bioprocesses is hampered by two major difficulties [6–9]. The first one is related to kinetics of these processes, which are strongly nonlinear and furthermore the kinetic and process parameters are often partially or completely unknown. Another difficulty lies in the absence, for many processes, of cheap and reliable instrumentation suited for real-time monitoring of the process variables. In order to overcome these difficulties, several strategies were developed, such as exactly linearizing control [6, 7, 10], adaptive approach proposed by Bastin and Dochain and applied in several works [6, 7, 10–12], optimal control, especially for fed-batch bioprocesses [6, 13], neural control [14], sliding mode control [15], robust and robust adaptive control [12, 16], and model predictive control [17]. Some problems which occur in implementation of numerical control of bioreactors are presented and analyzed in [18].

The difficulties encountered in the measurement of state variables impose the use of so-called “software sensors”—combinations between hardware sensors and software estimators [6, 7]. Note that these software sensors are used not only for the estimation of concentrations of some components but also for the estimation of kinetic parameters or even kinetic reactions. The interest for the development of software sensors for bioreactors is proved by the big number of publications and applications in this area [6–8, 19].

This paper, which is an extended version of Petre et al. [20], presents the design and the analysis of a nonlinear adaptive control strategy, capable to deal with the model uncertainties in an adaptive way, for a lactic fermentation bioprocess that is carried out in two continuous stirred tank bioreactors sequentially connected. In order to avoid the derivation of additional state observers, a new indirect adaptive control algorithm is presented. More exactly, the consumption rates of two limiting substrates are considered completely unknown and summarized in two unknown and time varying parameters, which are estimated by means of an appropriately observer-based estimator. This algorithm could be considered as an extension of the theory proposed by Bastin and Dochain [6] and improved by Ignatova et al. [11], since the substrate consumption rates are considered as completely unknown functions describing the whole kinetics.

The present work is focused in two directions. First, an adaptive control structure based on the nonlinear model of the process is designed as a combination of a linearizing control law and of a parameter estimator, used for the on-line estimation of bioprocess unknown kinetics. Second, by using the stability and convergence analysis of the proposed control scheme, a tuning procedure of the kinetic estimator design parameters is achieved.

The effectiveness and performance of both estimation and control algorithms are illustrated by simulations applied in the case of a lactic fermentation bioprocess for which kinetic dynamics are strongly nonlinear, time varying and completely unknown, and not all the state variables are measurable.

We note that the bioprocess model and the control goal are the same as in the Ben Youssef et al. [21], where an adaptive-multivariable predictive control law was proposed and analyzed. However, in the present work, an indirect multivariable adaptive linearizing controller is developed and analyzed. The control structure is derived by combining a linearizing control law with a new parameter estimator used for on-line estimation of the bioprocess unknown kinetics.

In the work by Ignatova et al. [11], by using the observer-based estimator proposed by Petre [22], an indirect adaptive scheme was derived and analyzed for a gluconic acid production bioprocess. In the present work, the multivariable control strategy presented in Petre [22], and further developed by Ignatova et al. [11], is improved and full stability and convergence proof are provided. Finally, a tuning procedure of kinetic estimator design parameters for a lactic acid production bioprocess was obtained.

The paper is organized as follows. Section 2 is devoted to a brief description and mathematical modeling of a lactic acid fermentation bioprocess. Some nonlinear and adaptive control strategies are proposed in Section 3. From the stability analysis of the adaptive control scheme, a tuning procedure of kinetic estimator design parameters is achieved. Simulations results presented in Section 4 illustrate the performance of the proposed control algorithms and, finally, Section 5 concludes the paper.

2. Process Modeling and Control Problem

Lactic acid has traditionally been used in the food industry as an acidulating and/or preserving agent, and in the biochemical industry for cosmetic and textile applications [21, 23, 24]. Recently, lactic acid fermentation has received much more attention because of the increasing demand for new biomaterials such as biodegradable and biocompatible polylactic products. Two major factors limit its biosynthesis, affecting growth, and productivity: the nutrient limiting conditions and the inhibitory effect caused by lactic acid accumulation in the culture broth [21]. It was shown that, in addition to the inhibition effect of lactic acid, the richness of the fermentation medium may have a strong influence on growth dynamics and that the nutritional limitations during culture may interfere with inhibition by acid lactic [21]. A reliable model that explicitly integrates nutritional factor effects on both growth and lactic acid production in a batch fermentation process implementing Lb. casei was developed by Ben Youssef et al. [21] and it is described by the following basic differential equations:𝑋̇=𝜇𝑋−𝑘_𝑑𝑋,𝑃̇=𝜈_𝑝𝑋,𝑆̇=−𝑞_𝑠𝑋, (2.1) where 𝑋, 𝑆, and 𝑃 are, respectively, the concentrations of cells, substrate (glucose) and lactic acid. 𝜇, 𝜈_𝑝, and 𝑞_𝑠 correspond, respectively, to specific growth rate of the cells, specific rate of lactic acid production, and to specific rate of glucose consumption. 𝑘_𝑑 is the rate of cell death.

Since in lactic acid fermentation process the main cost of raw material comes from the substrate and nutrient requirements, in [21] some possible continuous-flow control strategies that satisfy the economic aspects of lactic acid production were investigated. The advantage of implementing continuous-flow process is that the main product, which is also an inhibitor, is continuously withdrawn from the system. Much more, according to microbial engineering theory, for a product-inhibited reaction like lactic acid or alcoholic fermentation [21, 25], a multistage system composed of many interconnected continuous stirred tank reactors, where in the different reactors some physiological states of microbial culture (substrate, metabolites) can be kept constant to some optimal values, may be a good idea.

Therefore, the model (2.1) can be extended to a continuous-flow process that is carried out in two continuous stirred tank reactors sequentially connected, as in the schematic view presented in Figure 1.

Figure 1 

A cascade of two reactors for the lactic acid production.

For this two-stage continuous flow bioreactor, the mathematical model is given by the following set of differential equations, each stage being of same constant volume V.

First Stage:
𝑋̇₁=𝜇₁−𝑘_𝑑𝑋₁−𝐷₁𝑋₁,𝑃̇₁=𝜈_𝑝1𝑋₁−𝐷₁𝑃₁𝑆̇₁,=−𝑞_𝑠1𝑋₁+𝐷₁₁𝑆₁^𝑖𝑛−𝐷₁𝑆₁,𝛼̇₁=𝐷₁₂𝛼₁^𝑖𝑛−𝐷₁𝛼₁, (2.2)(2.3)(2.4)(2.5) where 𝐷₁=𝐷₁₁+𝐷₁₂.

Second Stage:
𝑋̇₂=𝜇₂−𝑘_𝑑𝑋₂+𝐷₁𝑋₁−𝐷₁+𝐷₂𝑋₂,𝑃̇₂=𝜈_𝑝2𝑋₂+𝐷₁𝑃₁−𝐷₁+𝐷₂𝑃₂𝑆̇₂,=−𝑞_𝑠2𝑋₂+𝐷₁𝑆₁+𝐷₂𝑆₂^𝑖𝑛−𝐷₁+𝐷₂𝑆₂,𝛼̇₂=𝐷₁𝛼₁−𝐷₁+𝐷₂𝛼₂, (2.6)(2.7)(2.8)(2.9) where 𝑋_𝑖,𝑆_𝑖,𝑃_𝑖, and 𝛼_𝑖, (𝑖=1,2) are, respectively, the concentrations of biomass (cells), substrate, lactic acid, and enrichment factor in each bioreactor. 𝜇_𝑖,𝜈_𝑝𝑖, and 𝑞_𝑠𝑖(𝑖=1,2) correspond, respectively, to specific growth rate of the cells, specific rate of lactic acid production, and specific rate of glucose consumption in each bioreactor. 𝐷₁₁ is the first-stage dilution rate of a feeding solution at an influent glucose concentration 𝑆₁^𝑖𝑛. 𝐷₁₂ is the first-stage dilution rate of a feeding solution at an influent enrichment factor 𝛼₁^𝑖𝑛. 𝐷₂ is the influent dilution rate added at the second stage and 𝑆₂^𝑖𝑛 is the corresponding feeding glucose concentration. It can be seen that in the second stage no growth factor feeding tank is included since this was already feeding in the first reactor.

In the model (2.2)–(2.9), the mechanism of cell growth, the specific lactic acid production rate and the specific consumption rate in each bioreactor are given by [21]: 𝜇_𝑖=𝜇_max𝑖𝐾^gc_𝑃𝑖𝐾^gc_𝑃𝑖+𝑃_𝑖𝑆_𝑖𝐾^gc_𝑆+𝑆_𝑖𝑃1−^𝑖𝑃^gc_𝐶,𝜈_𝑝𝑖=𝜂𝜇_𝑖𝑆+𝛽^𝑖𝐾^rc_𝑆𝑖+𝑆_𝑖,𝑞_𝑠𝑖=𝜈_𝑝𝑖𝑌_𝑃𝑆,𝑖=1,2, (2.10) with 𝜇_max𝑖 the maximum specific growth rate, 𝐾^gc_𝑃𝑖 the lactic acid inhibition constant, 𝐾^gc_𝑆 the affinity constant of the growing cells for glucose, 𝑃^gc_𝐶 the critical lactic acid concentrations, 𝐾^rc_𝑆𝑖 is the affinity constant of the resting cells for glucose, and 𝑌_𝑃𝑆 is the constant substrate-to-product conversion yield. The superscript gc denotes the parameters related to growing cells and rc to that of the resting cells, and 𝜂 and 𝛽 are positive constants.

The kinetic parameters of this model may be readjusted depending on the medium enrichment factor 𝛼_𝑖 as follows [21]:𝜇_max𝑖=𝜇_max𝛼_𝑖−𝛼₀𝐾_𝛼𝜇+𝛼_𝑖−𝛼₀,𝐾^gc_𝑃𝑖=𝐾^gc_𝑃max𝛼_𝑖−𝛼₀𝐾_𝛼𝑃+𝛼_𝑖−𝛼₀,𝐾^rc_𝑆𝑖=𝐾^rc_𝑆max𝛼_𝑖−𝛼₀𝐾_𝛼𝑆+𝛼_𝑖−𝛼₀, (2.11) where 𝛼₀ is the minimal nutritional factor necessary for growth, 𝐾_𝛼𝜇, 𝐾_𝛼𝑃 and 𝐾_𝛼𝑆 are saturation constants. 𝜇_max, 𝐾^gc_𝑃max and 𝐾^rc_𝑆max correspond to the limit value of each kinetic parameter. Note that the hyperbolic-type expressions in (2.11) quantify the nutritional limitations and are based on the growth model proposed by Bibal [26].

For this two-stage reactor configuration, the operating point of the continuous lactic acid fermentation process could be adjusted by acting on at least two control inputs, that is, the glucose feeding flow rates both in the first and in the second reactor. These considerations show the multivariable nature of the control problem.

3. Control Strategies

3.1. Control Objective

As it was formulated in the previous section, the control objective consists in adjusting the plant’s load in order to convert the glucose into lactic acid via fermentation, which is directly correlated to the economic aspects of lactic acid production. More exactly, considering that the process model (2.2)–(2.9) is incompletely known, its parameters are time varying and not all the states are available for measurements, the control goal is to maintain the process at some operating points, which correspond to a maximal lactic production rate and a minimal residual glucose concentration. After a detailed process steady-state analysis, in [21] it has been demonstrated that these desiderata can be satisfied if the operating point of lactic acid fermentation is kept around the points 𝑆^∗₁=3g/L and 𝑆^∗₂=5g/L. This choice is the best option in order to simultaneously satisfy the both objectives: maximal lactic production rate and minimal residual glucose concentration [21]. As control variables we select the glucose feeding flow rates both in the first and in the second reactor denoted by 𝐹₁ and 𝐹₂, respectively, where 𝐹₁=(𝐷₁−𝐷₁₂)𝑆₁^𝑖𝑛 and 𝐹₂=𝐷₂𝑆₂^𝑖𝑛. In this way we obtain a multivariable control problem with two inputs: 𝐹₁ and 𝐹₂ and two outputs: 𝑆₁ and 𝑆₂.

In the following, in Section 3.2 we will present the design of an exactly linearizing feedback controller. This controller will be used as a benchmark in order to compare its behavior with the behavior of the indirect multivariable adaptive controller developed in Section 3.3. Also, in this subsection, the parameter estimator used in the indirect multivariable adaptive controller will be designed. In Section 3.4 the tuning procedure of estimator design parameters via stability and convergence analysis of the control scheme will be derived.

3.2. Exactly Linearizing Feedback Controller

Firstly, we consider the ideal case where maximum prior knowledge concerning the process (kinetics, yield coefficients, and state variables) is available, and the relative degree of differential equations in process model is equal to 1. Assume that for the two interconnected reactors we wish to have the following first order linear stable closed loop (process + controller) behavior:𝑆^∗₁̇−𝑆₁̇𝑆^∗₂̇−𝑆₂̇+𝜆₁00𝜆₂𝑆^∗₁−𝑆₁𝑆^∗₂−𝑆₂=0,𝜆₁,𝜆₂>0, (3.1) where 𝑆^∗₁ and 𝑆^∗₂ are the desired values of 𝑆₁ and 𝑆₂.

Since (2.4) and (2.8) in the process model have the relative degree equal to 1, these can be considered as an input-output model and can be rewritten in the following form:𝑆₁̇𝑆₂̇=−𝑞_𝑠1𝑋₁−𝑞_𝑠2𝑋₂−𝐷₁0−𝐷₁𝐷₁+𝐷₂𝑆₁𝑆₂+𝐹₁𝐹₂. (3.2) Then from (3.1) and (3.2) one obtains the following exactly multivariable decoupling linearizing feedback control law:𝐹₁𝐹₂=𝑆^∗₁̇𝑆^∗₂̇+𝜆₁00𝜆₂𝑆^∗₁−𝑆₁𝑆^∗₂−𝑆₂+𝑞_𝑠1𝑋₁𝑞_𝑠2𝑋₂+𝐷₁0−𝐷₁𝐷₁+𝐷₂𝑆₁𝑆₂. (3.3)

The control law (3.3) applied to the process (3.2) will result in:𝑆₁̇𝑆₂̇=𝜆₁00𝜆₂𝑆^∗₁−𝑆₁𝑆^∗₂−𝑆₂+𝑆^∗₁̇𝑆^∗₂̇, (3.4) and leads to the following linear error models:𝑒₁̇=−𝜆₁𝑒₁,𝑒₂̇=−𝜆₂𝑒₂, (3.5) where 𝑒₁=𝑆^∗₁−𝑆₁ and 𝑒₂=𝑆^∗₂−𝑆₂ represent the tracking errors. It is clear that for 𝜆₁,𝜆₂>0, the error models (3.5) have an exponential stable point at 𝑒₁=0 and 𝑒₂=0.

3.3. An Indirect Multivariable Adaptive Controller

In the previous subsection, it was assumed that the functional forms of the nonlinearities as well as the process parameters are known. Since such prior knowledge is not realistic, in this subsection we consider that the kinetics functions are completely unknown. So, the substrate consumption rates −𝑞_𝑠1𝑋₁ and −𝑞_𝑠2𝑋₂ in (3.2) can be expressed as follows:−𝑞_𝑠1𝑋₁=𝜌₁,−𝑞_𝑠2𝑋₂=𝜌₂, (3.6) where 𝜌₁ and 𝜌₂ are two unknown time-varying parameters. In this case, the control law (3.3) becomes:𝐹₁𝐹₂=𝑆^∗₁̇𝑆^∗₂̇+𝜆₁00𝜆₂𝑆^∗₁−𝑆₁𝑆^∗₂−𝑆₂−𝜌₁𝜌₂+𝐷₁0−𝐷₁𝐷₁+𝐷₂𝑆₁𝑆₂, (3.7) where the unknown parameters 𝜌₁ and 𝜌₂ will be substituted by their on-line estimates 𝜌₁ and 𝜌₂ calculated by using an observer-based parameter estimator (OBE) [22] applied only for the dynamics of 𝑆₁ and 𝑆₂ from (2.4) and (2.8) rewritten as follows:𝑆₁̇=𝜌₁+𝐹₁−𝐷₁𝑆₁,𝑆₂̇=𝜌₂+𝐷₁𝑆₁+𝐹₂−𝐷₁+𝐷₂𝑆₂. (3.8)

Under these notations, the algorithm for on-line computation of 𝜌₁ and 𝜌₂ is particularized as follows:𝑆₁̇=𝜌₁+𝐹₁−𝐷₁𝑆₁+𝜔₁𝑆−𝑆₁,𝑆₂̇=𝜌₂+𝐷₁𝑆₁+𝐹₂−𝐷₁+𝐷₂𝑆₂+𝜔₂𝑆₂−𝑆₂,𝜌₁̇=𝛾₁𝑆−𝑆₁̇,𝜌₂=𝛾₂𝑆₂−𝑆₂, (3.9)(3.10)(3.11)(3.12) where 𝛾₁,𝛾₂>0 are the gains of the updating laws (3.11) and (3.12), and 𝜔₁,𝜔₂>0 are design parameters to control the stability and the tracking properties of the estimator.

Remark 3.1. The parameter error convergence for the estimator (3.9)–(3.12) is assured if some signals in system satisfy the well-known persistency of excitation condition [4, 5]: A piecewise continuous signal vector 𝜙∶ℜ⁺→ℜ^𝑛 is referred to as persistent excitation with a level of excitation 𝜅₀>0 if there exist constants 𝜅₁,𝑇₀>0 such that𝜅₁𝐼≥(1/𝑇₀)∫^𝑡+𝑇0_𝑡𝜙(𝜏)𝜙^𝑇(𝜏)𝑑𝜏≥𝜅₀𝐼, for all 𝑡≥0. Although the matrix 𝜙(𝜏)𝜙^𝑇(𝜏) is singular for each 𝜏, this condition requires that 𝜙(𝑡) varies in such a way with time that the integral of the matrix 𝜙(𝜏)𝜙^𝑇(𝜏) is uniformly positive definite over any time interval [𝑡, 𝑡+𝑇₀]. In our particular case it is necessary that, for example, 𝐷₁, 𝐷₂ or some kinetic parameters to be persistently exciting (which is generically fulfilled).

Finally, an indirect multivariable adaptive linearizing controller is obtained by combination of (3.9)–(3.12) and (3.7) rewritten as follows: 𝐹₁𝐹₂=𝑆^∗₁̇𝑆^∗₂̇+𝜆₁00𝜆₂𝑆^∗₁−𝑆₁𝑆^∗₂−𝑆₂−𝜌₁𝜌₂+𝐷₁0−𝐷₁𝐷₁+𝐷₂𝑆₁𝑆₂. (3.13)

A block diagram of the designed indirect multivariable adaptive control system is shown in Figure 2.

Figure 2 

The structure of the indirect adaptive control system.

Remark 3.2. Since the on-line measurement of the lactic acid concentrations is not possible, in order to provide on-line evaluation of the current lactic acid production, an asymptotic state observer is derived as follows.
Let us define the auxiliary variables 𝑧₁ and 𝑧₂ as: 𝑧₁=𝑃₁+𝑌_𝑃𝑆𝑆₁,𝑧₂=𝑃₂+𝑌_𝑃𝑆𝑆₂, (3.14) which use only the measurements of the outputs 𝑆₁ and 𝑆₂ that are available on-line. The dynamics of 𝑧_𝑖,𝑖=1,2, deduced from model (2.2)–(2.9), are expressed by the following linear stable equations: z₁̇̂=−𝐷₁z₁̂+𝑌_𝑃𝑆𝐹₁,z₂̇̂=𝐷₁z₁̂−𝐷₁+𝐷₂z₂̂+𝑌_𝑃𝑆𝐹₂, (3.15) that are independent of the process kinetics. Then, from (3.14) and (3.15), the on-line estimation of 𝑃₁ and 𝑃₂ are given by the following relations: 𝑃₁=𝑧₁−𝑌_𝑃𝑆𝑆₁,𝑃₂=𝑧₂−𝑌_𝑃𝑆𝑆₂. (3.16)

3.4. Tuning Procedure of Estimator Design Parameters via Stability and Convergence Analysis

Let us denote by 𝑦=[𝑆₁𝑆₂]^𝑇 the vector containing the concentrations of controlled feeding substrates and by 𝑢=[𝐹₁𝐹₂]^𝑇 the vector of the input control variables, that is, the vector of the mass feeds of substrates 𝑆₁ and 𝑆₂, respectively. Then (3.8) can be written as 𝑦̇=𝜌−𝐷𝑦+𝑢, (3.17) where 𝜌=[𝜌₁𝜌₂]^𝑇 is the unknown time-varying parameter vector and 𝐷 is the dilution rate matrix, whose structure is as follows:𝐷𝐷=¹0−𝐷₁𝐷₁+𝐷₂. (3.18)

The OBE (3.9)–(3.12) can be represented in a matrix form as follows: ̇,̇,𝑦=𝜌−𝐷𝑦+𝑢+Ω𝑦−𝑦𝜌=Γ𝑦−𝑦 (3.19) where 𝜌 is the estimation of unknown parameter 𝜌, and Ω,Γ∈ℜ^2×2 are positive definite diagonal matrices whose entries 𝜔₁,𝜔₂>0 and 𝛾₁,𝛾₂>0, respectively, represent the estimator design parameters.

Under the previous notations, the exactly multivariable decoupling linearizing feedback control law (3.3), respectively, (3.7) can be written as:𝑢=𝑦^∗̇+Λ𝑦^∗−𝑦−𝜌+𝐷𝑦, (3.20) where 𝑦^∗=[𝑆^∗₁𝑆^∗₂]^𝑇, and the indirect adaptive linearizing control law (3.13) can be written as: 𝑢=𝑦^∗̇+Λ𝑦^∗−𝑦−𝜌+𝐷𝑦, (3.21) where Λ=^𝜆10_0𝜆2>0.

Case 1. By defining the errors 𝑦=𝑦−𝑦 and 𝜌=𝜌−𝜌, in the case when the unknown parameter 𝜌 is constant, the following error system can be derived from (3.17) and (3.19): 𝑦̇=−Ω𝑦+𝜌,𝜌̇=−Γ𝑦. (3.22)
It can be seen that the error system in (3.22) is a linear time-invariant system that can be rewritten as: ̇̇=𝑦𝜌−Ω𝐼₂−Γ0𝑦𝜌, (3.23) where 𝐼₂ is the identity matrix of order two. The stability of the equilibrium (𝑦=0, 𝜌=0) of (3.22) can be proved by using the well-known Hurwitz criterion. So, by using the coefficients of the characteristic polynomial associated to (3.23), given by 𝑃(𝜆)=𝜆⁴+𝜔₁+𝜔₂𝜆³+𝜔₁𝜔₂+𝛾₁+𝛾₂𝜆²+𝛾₂𝜔₁+𝛾₁𝜔₂𝜆+𝛾₁𝛾₂=𝑎₄𝜆⁴+𝑎₃𝜆³+𝑎₂𝜆²+𝑎₁𝜆+𝑎₀, (3.24)after straightforward calculation we obtain that all Hurwitz determinants are positive for any 𝜔₁,𝜔₂>0 and 𝛾₁,𝛾₂>0. Indeed, by using the coefficients of the above characteristic polynomial we can construct the following Hurwitz determinant: Δ₄=𝑎₃𝑎₁𝑎00₄𝑎₂𝑎₀00𝑎₃𝑎₁00𝑎₄𝑎₂𝑎₀=𝜔₁+𝜔₂𝛾₂𝜔₁+𝛾₁𝜔₂001𝜔₁𝜔₂+𝛾₁+𝛾₂𝛾₁𝛾₂00𝜔₁+𝜔₂𝛾₂𝜔₁+𝛾₁𝜔₂001𝜔₁𝜔₂+𝛾₁+𝛾₂𝛾₁𝛾₂. (3.25) Now we construct and calculate all Hurwitz determinants as follows: Δ₁=𝑎₃=𝜔₁+𝜔₂Δ>0,₂=𝑎₃𝑎₁𝑎₄𝑎₂=𝜔₁+𝜔₂𝛾₂𝜔₁+𝛾₁𝜔₂1𝜔₁𝜔₂+𝛾₁+𝛾₂=𝜔₁𝜔₂𝜔₁+𝜔₂+𝛾₂𝜔₁+𝛾₁𝜔₂Δ>0,₃=𝑎₃𝑎₁0𝑎₄𝑎₂𝑎₀0𝑎₃𝑎₁=𝜔₁+𝜔₂𝛾₂𝜔₁+𝛾₁𝜔₂01𝜔₁𝜔₂+𝛾₁+𝛾₂𝛾₁𝛾₂0𝜔₁+𝜔₃𝛾₂𝜔₁+𝛾₁𝜔₂=𝜔₁𝜔₂𝜔₁+𝜔₂𝛾₂𝜔₁+𝛾₁𝜔₂+𝛾₁−𝛾₂²Δ>0,₄=𝛾₁𝛾₂⋅Δ₃>0. (3.26)
According to Hurwitz criterion [27, 28], the system (3.23) is stable since all Hurwitz determinants are positive if 𝜔₁,𝜔₂>0 and 𝛾₁,𝛾₂>0. Then, the matrix associated to (3.23) is a Hurwitz matrix and the equilibrium (𝑦=0, 𝜌=0) of (3.22) is uniformly asymptotically stable.

Case 2. In the case when the unknown parameter 𝜌 is time-varying, the stability analysis of the proposed indirect multivariable adaptive linearizing controller is carried out under the following realistic assumptions.(A1) The substrate concentrations 𝑆₁ and 𝑆₂, the components of the dilution rate matrix 𝐷, and the substrate feed rates are on-line measured. (A2) The substrate consumption rates 𝜌₁ and 𝜌₂ are fully unknown.(A3) The reference signal vector 𝑦^∗ and its time derivative are piecewise continuous bounded functions of time, and the measured substrate feed rates are continuous bounded functions of time. (A4) The measurements of substrate vector, denoted 𝑦_𝑚(𝑡), are corrupted by an additive noise vector 𝑤(𝑡)=[𝑤₁𝑤₂]^𝑇: 𝑦_𝑚(𝑡)=𝑦(𝑡)+𝑤(𝑡), for all 𝑡.(A5) The components of time varying parameter vector 𝜌(𝑡) are bounded as: 0>𝜌_𝑖(𝑡)≥𝜌^∗_𝑖, 𝜌^∗_𝑖<0, 𝑖=1,2, for all 𝑡.(A6) The time derivatives of the components of 𝜌(𝑡) are bounded as: |𝜌̇_𝑖(𝑡)|≤𝑚_1𝑖, 𝑚_1𝑖>0, 𝑖=1,2, for all 𝑡.(A7) The measurement noise vector is bounded as: |𝑤(𝑡)|≤𝑚_2𝑖, 𝑚_2𝑖>0, 𝑖=1,2, for all 𝑡.(A8) The components of dilution rates matrix 𝐷 are known and bounded: 0≤𝐷_𝑖(𝑡)≤𝐷_𝑖^max, 𝑖=1,2, 0≤𝐷₁₂(𝑡)≤𝐷^max₁₂, for all 𝑡.
Under these assumptions, the parameter estimator (3.19) and the adaptive control law (3.21) are rewritten as follows: 𝑦̇=𝜌−𝐷𝑦_𝑚+𝑢+Ω𝑦_𝑚,̇𝑦−𝑦𝜌=Γ_𝑚−𝑦,𝑢=𝑦^∗̇+Λ𝑦^∗−𝑦_𝑚−𝜌+𝐷𝑦_𝑚. (3.27)(3.28)(3.29)
From (3.17), (3.27), and (3.28) under the condition (A4), the following error system can be derived as follows: 𝑦̇=𝑦̇−𝑦̇=𝜌−𝐷𝑦+𝑢−𝜌+𝐷𝑦_𝑚−𝑢−Ω𝑦_𝑚()̇̇−𝑦=𝜌+𝐷−𝑦+𝑦+𝑤−Ω𝑦+𝑤−𝑦=−Ω𝑦+𝜌+𝐷−Ω𝑤,𝜌=𝜌̇−𝜌=𝜌̇−Γ𝑦+𝑤−𝑦=−Γ𝑦−Γ𝑤+𝜌̇, (3.30) which can be rewritten in the following form: 𝑦̇₁=−𝜔₁𝑦₁+𝜌₁+𝐷₁−𝜔₁𝑤₁,𝜌̇₁=−𝛾₁𝑦₁−𝛾₁𝑤₁+𝜌̇₁,𝑦̇₂=−𝜔₂𝑦₂+𝜌₂+𝐷₁+𝐷₂−𝜔₂𝑤₂−𝐷₁𝑤₁,𝜌̇₂=−𝛾₂𝑦₂−𝛾₂𝑤₂+𝜌̇₂. (3.31)

Defining the state vector 𝑥 as 𝑥=[𝑥^𝑇₁𝑥^𝑇₂]^𝑇 with 𝑥₁=[𝑦₁𝜌₁]^𝑇 and 𝑥₂=[𝑦₂𝜌₂]^𝑇, then (3.31) can be rewritten as: 𝑥̇=𝐴𝑥+𝑏, (3.32) where 𝐴=diag{𝐴_𝑖}, 𝑖=1,2 and:𝐴_𝑖=−𝜔_𝑖1−𝛾_𝑖0𝑏,𝑏=¹𝑏₂,𝑏₁=𝐷₁−𝜔₁𝑤₁−𝛾₁𝑤₁+𝜌̇₁,𝑏₂=𝐷₁+𝐷₂−𝜔₂𝑤₂−𝐷₁𝑤₁−𝛾₂𝑤₂+𝜌̇₂. (3.33)

The two eigenvalues of 𝐴_𝑖 denoted by 𝑣_𝑖1 and 𝑣_𝑖2 are related to 𝜔_𝑖 and 𝛾_𝑖 as follows: 𝑣_𝑖1+𝑣_𝑖2=−𝜔_𝑖<0, 𝑣_𝑖1𝑣_𝑖2=𝛾_𝑖>0, 𝑖=1,2. It can be seen that 𝑣_1𝑖,𝑣_2𝑖<0. The next assumption is considered for the estimator design parameters.(A9) The design parameters 𝜔_𝑖 and 𝛾_𝑖 are chosen such that 𝐴_𝑖 has real distinct eigenvalues [6, 11], for example, 𝑣_𝑖2<𝑣_𝑖1<0, where 𝑣_𝑖1,2=(−𝜔_𝑖±𝜔²_𝑖−4𝛾_𝑖)/2, that implies 0<𝛾_𝑖<𝜔²_𝑖/4.

The system (3.32) describes the dynamics of the estimation errors of the observer-based estimator including the observations of the measured substrates and the estimations of their unknown consumption rates 𝜌_𝑖.

We can consider that the estimator (3.27)-(3.28) is a set of quasi-independent second order linear systems because one controlled substrate concentration is associated to one unknown parameter of the vector 𝜌 only. Since the dynamics of estimation errors (3.32) is described also by a linear differential equation with a block diagonal matrix with two 2×2 stable blocks, then, we can have the following stability result, which is an extension of the result presented in [6, 11]:

Theorem 3.3. Under assumptions (A1)–(A9), the estimation errors 𝑦_𝑖 and 𝜌_𝑖 are asymptotically bounded for all 𝑡 as follows:(i)lim_𝑡→∞sup𝑦_𝑖(𝑡)≤2𝛿𝐵_𝑖1𝜔²_𝑖−4𝛾_𝑖+𝐵_𝑖2𝛾_𝑖,𝑖=1,2, (3.34)(ii)lim_𝑡→∞sup𝜌_𝑖(𝑡)≤𝐵_𝑖1+𝜔_𝑖𝐵_𝑖2𝛾_𝑖,𝑖=1,2 (3.35)and then the tracking error 𝑦^∗−𝑦 is bounded as follows:(iii)lim_𝑡→∞𝑦^∗_𝑖−𝑦_𝑖≤𝑀_𝑖𝑖𝜆_𝑖,𝑖=1,2, (3.36)where 𝐵₁₁=𝑚₂₁(𝐷₁^max+𝜔₁), 𝐵₂₁=(𝐷₁^max+𝐷₂^max+𝜔₂)𝑚₂₂+𝐷₁^max𝑚₂₁, 𝐵_𝑖2=𝑚_2𝑖𝛾_𝑖+𝑚_1𝑖, 𝛿=(𝑣_1𝑖/𝑣_2𝑖)^{(𝑣1𝑖/(𝑣1𝑖−𝑣2𝑖))}−(𝑣_1𝑖/𝑣_2𝑖)^{(𝑣2𝑖/(𝑣1𝑖−𝑣2𝑖))}, and 𝑀_𝑖𝑖>0, 𝑖=1,2 are constant bounds which will be defined.

Proof. From assumptions (A6)–(A8), it follows that the input vector 𝑏=[𝑏^𝑇₁,𝑏^𝑇₂]^𝑇 of the system (3.32) is bounded. Therefore, (3.32) behaves as a linear time-invariant differential system with a bounded input, whose time response is: 𝑥_𝑖(𝑡)=Φ_𝑖(𝑡)𝑥_𝑖(0)+^𝑡₀Φ_𝑖(𝑡−𝜏)𝑏_𝑖(𝜏)𝑑𝜏,𝑖=1,2, (3.37) where Φ_𝑖(𝑡) is the state transition matrix associated to 𝐴_𝑖 defined as: Φ_𝑖(𝑡)=𝑒^𝐴𝑖𝑡=𝜑^𝑖₁₁(𝑡)𝜑^𝑖₁₂(𝑡)𝜑^𝑖₂₁(𝑡)𝜑^𝑖₂₂(𝑡), (3.38)with 𝜑^𝑖₁₁(𝑡)=(1/(𝜈_𝑖1−𝜈_𝑖2))(𝜈_𝑖1𝑒^{𝜈𝑖1𝑡}−𝜈_𝑖2𝑒^{𝜈𝑖2𝑡}), 𝜑^𝑖₁₂(𝑡)=(1/(𝜈_𝑖1−𝜈_𝑖2))(𝑒^{𝜈𝑖1𝑡}−𝑒^{𝜈𝑖2𝑡}), 𝜑^𝑖₂₁(𝑡)=−𝛾_𝑖𝜑^𝑖₁₂(𝑡), 𝜑^𝑖₂₂(𝑡)=𝜔_𝑖𝜑^𝑖₁₂(𝑡)+𝜑^𝑖₁₁(𝑡). Since 𝐴_𝑖 is a stable matrix, then the state 𝑥_𝑖 is bounded (see [5, Lemma 3.3.5] and [6, Theorem A2.6]). A bound of each component of 𝑥_𝑖(𝑡) in (3.37) is given by: 𝑥_𝑖1(𝑡)≤𝜑^𝑖₁₁(𝑡)𝑥_𝑖1(0)+𝜑^𝑖₁₂(𝑡)𝑥_𝑖2(0)+^𝑡₀𝜑^𝑖₁₁(𝑡−𝜏)𝑏_𝑖1(𝑡)𝑑𝜏+^𝑡₀𝜑^𝑖₁₂(𝑡−𝜏)𝑏_𝑖2(𝑡𝑥)𝑑𝜏,_𝑖2(𝑡)≤𝜑^𝑖₂₁(𝑡)𝑥_𝑖1(0)+𝜑^𝑖₂₂(𝑡)𝑥_𝑖2(0)+^𝑡₀𝜑^𝑖₂₁(𝑡−𝜏)𝑏_𝑖1(𝑡)𝑑𝜏+^𝑡₀𝜑^𝑖₂₂(𝑡−𝜏)𝑏_𝑖2(𝑡)𝑑𝜏. (3.39) Then, taking the limit of the above expressions for 𝑡→∞ and using Technical Lemma A2.2 from [6], we obtain the following inequalities: lim_𝑡→∞sup𝑥_𝑖1(𝑡)≤2𝐵_𝑖1𝛿𝜔²_𝑖−4𝛾_𝑖+𝐵_𝑖2𝛾_𝑖,lim_𝑡→∞sup𝑥_𝑖2(𝑡)≤𝐵_𝑖1+𝜔_𝑖𝛾_𝑖𝐵_𝑖2, (3.40)where 𝐵_𝑖1>0 and 𝐵_𝑖2>0 are the bounds of 𝑏_𝑖1 and 𝑏_𝑖2, that is, |𝑏_𝑖1|≤𝐵_𝑖1 and |𝑏_𝑖2|≤𝐵_𝑖2, 𝑖=1,2, with 𝐵_𝑖1 and 𝐵_𝑖2 defined in Theorem 3.3. Thus, the boundedness of 𝑦_𝑖 and 𝜌_𝑖 given in (i) and (ii), respectively, is proven.
To prove (iii) we proceed as follows. We can see that the adaptive control law (3.29) with parameter estimator (3.27) and (3.28) applied to system (3.17) leads to the following dynamics of the closed-loop system: 𝑦̇−𝑦^∗̇+Λ𝑦−𝑦^∗+Λ𝑤=𝑦̇+Ω𝑦+𝑤. (3.41) Since Ω is a diagonal stable matrix, then from (3.41) it can be seen that the tracking error 𝑦^∗−𝑦 is the output of a linear stable filter driven by the observer error 𝑦. From the first equation in (3.30) we can write: 𝑦̇+Ω(𝑦+𝑤)=𝜌+𝐷𝑤. Now, by using (3.35), and assumptions (A7) and (A8), the following inequalities can be written: 𝑦̇₁+𝜔₁𝑦₁+𝑤₁=𝜌₁+𝐷₁𝑤₁≤𝜌₁+𝐷₁𝑤₁≤𝐵₁₁+𝑚₂₁𝐷₁^max+𝜔₁𝐵₁₂𝛾₁=𝑀₁,𝑦̇₂+𝜔₂𝑦₂+𝑤₂=𝜌₂−𝐷₁𝑤₁+𝐷₁+𝐷₂𝑤₂≤𝜌₂+𝐷₁𝑤₁+𝐷₁+𝐷₂𝑤₂≤𝐵₂₁+𝜔₂𝐵₂₂𝛾₂+𝑚₂₂𝐷₁^max+𝐷₂^max+𝑚₂₁𝐷₁^max=𝑀₂. (3.42) From (3.41) and (3.42) one finds that: 𝑦_𝑖̇−𝑦^∗_𝑖̇+𝜆_𝑖𝑦_𝑖−𝑦^∗_𝑖≤−𝜆_𝑖𝑤_𝑖+𝑀_𝑖=𝜆_𝑖𝑚_2𝑖+𝑀_𝑖=𝑀_𝑖𝑖. (3.43)
From (3.43), after some straightforward calculation it can be obtained that lim_𝑡→∞|𝑦^∗_𝑖−𝑦_𝑖|≤𝑀_𝑖𝑖/𝜆_𝑖 and the proof is complete.

4. Simulation Results and Comments

The performance of the designed multivariable adaptive controller (3.9)–(3.13), by comparison to the exactly linearizing controller (3.3) (which yields the best response and can be used as benchmark), has been tested by performing extensive simulation experiments.

For a proper comparison, the simulations were carried out by using the process model (2.2)–(2.11) under identical and realistic conditions. The simulations were designed so that several set point changes on the controlled variables 𝑆₁ and 𝑆₂ occurred, nearby the operational point 𝑆^∗₁=3g/L, 𝑆^∗₂=5g/L.

The kinetic parameter values used in simulations are presented in Table 1 [21], and the initial values of variables considered in simulations in Table 2.

Table 1 

Parameter values.

Table 2 

Initial conditions.

The system’s behavior was analyzed assuming that the dilution rates of the two feeding substrates and the influent enrichment factor act as perturbations of the form:𝐷₁(𝑡)=𝐷₁₀⋅1−0.15sin𝜋𝑡,𝐷25₂(𝑡)=𝐷₂₀⋅1+0.15cos𝜋𝑡,𝛼50_1𝑖𝑛(𝑡)=𝛼^𝑖𝑛₁₀⋅1+0.25sin𝜋𝑡.20 (4.1)

Also, one of the kinetic parameters—the maximum specific growth rate—is considered time varying:𝜇_max(𝑡)=𝜇_max0⋅1−0.1sin𝜋𝑡40. (4.2)

The values of the estimator design parameters, chosen so that to satisfy the tuning conditions of assumption (A9), are: 𝜆₁=𝜆₂=0.45, 𝜔₁=1.75, 𝜔₂=0.75, 𝛾₁=0.25, 𝛾₂=0.1.

The behavior of the indirect adaptive controlled system by comparison to the behavior of closed loop system with exactly linearizing controller is shown in Figure 3, the time evolution of the two controlled variables 𝑆₁ and 𝑆₂, respectively, and in Figure 4, the control inputs 𝐹₁ and 𝐹₂, respectively.

Figure 3 

Time evolution of outputs—adaptive control (2) by comparison to exactly linearizing control (1).

Figure 4 

Profiles of control inputs—adaptive control (2) by comparison to exactly linearizing control (1).

In order to test the behavior of the indirect adaptive controlled system in more realistic circumstances, we considered that the measurements of both controlled variables (𝑆₁ and 𝑆₂) are corrupted with an additive white noise with zero average (5% from their nominal values). The simulation results in this case, conducted in the same conditions as in previous simulation, are presented in Figures 5 and 6, respectively. The simulation shows that the behavior of controlled variables and of control inputs is comparable with the results obtained in the free noise simulation.

Figure 5 

Evolution of outputs—adaptive control (noisy data) (2) by comparison to linearizing control (1).

Figure 6 

Control inputs—adaptive control (noisy data) (2) by comparison to linearizing control (1).

The time evolution of the estimates of the unknown functions (3.6) provided by the observer-based estimator is presented in Figures 7 and 8, in both simulation cases.

Figure 7 

Profile of estimates of unknown function ρ 1: (1) without and (2) with noisy measurements.

Figure 8 

Profile of estimates of unknown function ρ 2: (1) without and (2) with noisy measurements.

From Figures 7 and 8, it can be noticed that the time evolution of estimates for noisy measurements of 𝑆₁ and 𝑆₂ is similar with the time profiles in free noise case.

The production of lactic acid can be evaluated by using the data provided by the asymptotic observer (3.15) and (3.16). The time evolution of the lactic acid concentration estimates in the two reactors is depicted in Figure 9.

Figure 9 

Profile of estimates of lactic acid concentration: (1) without and (2) with noisy measurements.

From graphics in Figures 3–6 it can be seen that the behavior of overall system with indirect adaptive controller, even if this controller uses much less a priori information, is good, being very close to the behavior of closed loop system in the ideal case when the process model is completely known. Note also the regulation properties and ability of the controller to maintain the controlled outputs 𝑆₁ and 𝑆₂ very close to their desired values, despite the very high variations for 𝐷₁, 𝐷₂, 𝛼₁^𝑖𝑛, and 𝜇_max, time variation of process parameters and the influence of noisy measurements.

In order to compare our proposed strategy with some classical control technique, several simulations were performed. As exemplification, Figures 10 and 11 present the outputs and control inputs for the case when a classical PI (proportional-integrative) control is used, with the best tuning of control parameters. As can be noticed, the behavior of the classical control is strongly influenced by the noisy measurements and by the variation of kinetic and process parameters.

Figure 10 

Evolution of outputs—PI control (noisy data) (2) by comparison to PI (free noise) (1).

Figure 11 

Control inputs—PI control (noisy data) by comparison to PI (free noise).

5. Conclusion

This paper is concerned with the design of an indirect multivariable adaptive control for a lactic acid fermentation process that is carried out in a cascade of two stirred tank fermentation reactors.

The controller was achieved by assumption that the consumption rates of two limiting substrates are completely unknown and are summarized in two fully unknown time varying parameters which are estimated by means of appropriately observer-based estimators. A simple tuning method of kinetic estimator design parameters was achieved by stability and convergence analysis of the control scheme.

The simulation results showed that the performances of the proposed multivariable adaptive controller are very good. The derived control method can be applied also for a large class of square nonlinear bioprocesses.

Acknowledgment

This work was supported by CNCSIS-UEFISCSU, Romania, Project no. PNII-IDEI 548/2008.