Distribution-Based Identification of Yield Coefficients in a Baker’s Yeast Bioprocess

Abstract

A distribution-based identification procedure for estimation of yield coefficients in a baker’s yeast bioprocess is proposed. This procedure transforms a system of differential equations to a system of algebraic equations with respect to unknown parameters. The relation between the state variables is represented by functionals using techniques from distribution theory. A hierarchical structure of identification is used, which allows obtaining a linear algebraic system of equations in the unknown parameters. The coefficients of this algebraic system are functionals depending on the input and state variables evaluated through some test functions from distribution theory. First, only some state equations are evaluated throughout test functions to obtain a set of linear equations in parameters. The results of this first stage of identification are used to express other parameters by linear equations. The process is repeated until all parameters are identified. The performances of the method are analyzed by numerical simulations.

1. Introduction

A process carried out in a bioreactor can be defined as a set of 𝑚 biochemical reactions involving 𝑛 components (with 𝑛≥𝑚). The global dynamics can be represented by following the dynamical state-space model [1]:𝑑𝜉()𝑑𝑡=𝐾⋅𝜑𝜉,𝑡−𝐷𝜉+𝐹−𝑄, (1.1) where 𝜉∈ℜ^𝑛×1.

This model describes the behaviour of an entire class of biotechnological processes and is referred to as the general dynamical state-space model of this class of bioprocesses [1, 2]. In (1.1), the term 𝐾⋅𝜑(𝜉,𝑡) is the rate of consumption and/or production of the components in the reactor, that is, the reaction kinetics and the term −𝐷𝜉+𝐹−𝑄 represents the exchange with the environment, that is, the dynamics of the component transportation through the bioreactor. The strongly nonlinear character of the model (1.1) is given by the reaction kinetics. In many situations, the yield coefficients, the structure, and the parameters of the reaction rates are partially known or unknown. Many of the evolved control methods for these kind of systems—like model predictive control [3], robust or adaptive control [4, 5]—are based on good initial estimates of the yield parameters. This paper deals with the yield coefficients identification of baker’s yeast bioprocess.

In recent years, a progress has been made in the area of continuous-time system identification. Even if the most physical systems are naturally continuous, a much more attention has been paid to parameter estimation of discrete-time systems, mainly because they are better suited for numerical implementations. Continuous-time identification makes possible a more direct link to the physical properties and operation of the underlying systems [6, 7] and the direct estimation of physical parameters which have a clear significance [8]. One important direction in continuous-time system identification is to transform the system of differential equations to an algebraic system that reveals the unknown parameters. By using some measures, the direct computation of the input-output data derivatives can be completely avoided. For linear system identification, several methods are reported on this direction: identification based on the Laplace transformation and then use the Laguerre filter [9] or transforming the continuous-time system to the frequency domain [10].

A novel approach to identify the continuous-time system is the distribution-based approach, using deterministic distributions [11] or random distributions [12]. Through these approaches, derivatives are described in the sense of distribution theory and one constructs the input-output algebraic relationships using differential information produced in the distribution sense.

Identification of nonlinear continuous-time systems is far away more complicated. The traditional procedures are based on the Volterra functional series, expressed in time domain or frequency domain [7]. The parameter identification of deterministic nonlinear continuous-time systems (NCTSs), modelled by polynomial-type differential equations has been considered by numerous authors [13, 14]. They use the modulating functions as a sum of sinusoids, whose fundamental period has a fixed length, to convert differential equations into an algebraic form in the frequency domain. These methods can be applied only if the NCTS is exactly integrable. This means that all the terms can be written as pure derivatives of some computable function of the measured signals. The computational burden in applying the method of modulating functions is overtaken by using Fast Fourier Transform algorithms.

In this paper we propose an identification method for a class of NCTS considering that the unknown parameters can appear in rational relations with measured variables. Using techniques utilized in distribution approach the measurable functions and their derivatives are represented by functionals on a fundamental space of testing functions. Such systems are common in biotechnology [1, 15]. It is supposed that the system is described by state equations and all state variables are measurable (or they can be estimated using some state estimators [16]).

The most important approaches for the yield coefficients estimation make use of the state transformations based on the general structure [17–19]. The main idea of this paper is to use a hierarchical structure of identification to estimate the yield coefficients in the baker’s yeast bioprocess. First, some state equations are utilized to obtain a set of linear equations in some parameters. The results of this first stage of identification are used for expressing other parameters by linear equations. This process is repeated until all parameters are identified.

The paper is organized in the following way. The nonlinear dynamical model of a baker’s yeast process is given in Section 2. Section 3 presents the problem statement of continuous-time systems identification based on distribution approach and the hierarchical structure of yield coefficients identification. Section 4 is dedicated to the analysis of properties of the identification algorithm. Some experimental results are presented in Section 5, and conclusions in Section 6.

2. Nonlinear Dynamical Model of Baker’s Yeast Process

The yeast growth has already been the subject of numerous studies and papers. The reaction scheme of baker’s yeast bioprocess is [20]𝑆+𝐶^𝜇𝑜𝑠𝑆⟶𝑋+𝐺,_𝜇𝑟𝑠⟶𝑋+𝐸+𝐺,𝐸+𝐶⟶^𝜇𝑜𝑒𝑋+𝐺. (2.1)(2.2)(2.3)

In the reaction schemes (2.1)–(2.3), 𝑆 is the glucose, 𝑋 the biomass, 𝐶 the dissolved oxygen, 𝐸 the ethanol, and 𝐺 the dissolved carbon dioxide. The first reaction scheme represents the respiratory growth on glucose; the second reaction scheme represents the fermentative growth on glucose, and finally the third reaction represents the respiratory growth on ethanol. 𝜇^𝑜_𝑠,𝜇^𝑟_𝑠, and 𝜇^𝑜_𝑒 are three specific growth rates that reflect the capacity of the baker’s yeast to exploit three catabolic pathways for energy and material sources.

The dynamical model of the fed-batch bioprocess with the reaction schemes (2.1)–(2.3) can be obtained by using the mass balance of the components, considering that the bioreactor is well mixed, the yield coefficients, are constant and the dynamic of the gas phase can be neglected [20]:𝑑𝑋=𝜇𝑑𝑡^𝑜𝑠+𝜇^𝑟_𝑠+𝜇^𝑜_𝑒𝑋−𝐷𝑋,𝑑𝑆=𝑑𝑡−𝑘₁𝜇^𝑜_𝑠−𝑘₂𝜇^𝑟_𝑠𝑋+𝐷𝑆_in−𝑆,𝑑𝐸=𝑘𝑑𝑡³𝜇^𝑟_𝑠−𝑘₄𝜇^𝑜_𝑒𝑋−𝐷𝐸,𝑑𝐶=𝑑𝑡−𝑘₅𝜇^𝑜_𝑠−𝑘₆𝜇^𝑜_𝑒𝑋−𝐷𝐶+OTR,𝑑𝐺=𝑘𝑑𝑡⁷𝜇^𝑜_𝑠+𝑘₈𝜇^𝑟_𝑠+𝑘₉𝜇^𝑜_𝑒𝑋−𝐷𝐺+CTR,𝑑𝑉𝑑𝑡=𝐹_in=𝐷𝑉. (2.4)(2.5)(2.6)(2.7)(2.8)(2.9)

The dynamical model is expressed as a set of differential equations, in terms of the concentrations of components, and with the additional (2.9), which represents the dynamics of the solution volume 𝑉 in the bioreactor. In (2.4)–(2.9), 𝑘_𝑖,𝑖=1,9 are the positive yield coefficients, 𝑆_in is the glucose concentration on the feed, OTR is the oxygen rate, defined as OTR=𝐾_𝐿𝑎(𝐶^∗−𝐶) (where 𝐾_𝐿𝑎 is the mass transfer coefficient and 𝐶^∗ the equilibrium concentration of the dissolved oxygen), CTR is the carbon dioxide transfer rate (defined as CTR=𝐾_𝑉𝐾_𝐿𝑎𝐺, with 𝐾_𝑉 a transfer coefficient), 𝐹_in is the input feed rate, and 𝐷=𝐹_in/𝑉 is the so-called dilution rate.

The model (2.4)–(2.8) can be expressed in the matrix form:𝑋̇𝑆̇𝐸̇𝐶̇𝐺̇=111−𝑘₁−𝑘₂00𝑘₃−𝑘₄−𝑘₅0−𝑘₆𝑘₇𝑘₈𝑘_9𝐾𝜇^𝑜_𝑠𝜇^𝑟_𝑠𝜇^𝑜_𝑒𝑋𝑆𝐸𝐶𝐺+0𝑋−𝐷𝐷𝑆_in00−0000.OTRCTR (2.10)

The next notations will be used: the state vector is 𝜉=[𝑋𝑆𝐸𝐶𝐺]^𝑇=[𝜉₁𝜉₂𝜉₃𝜉₄𝜉₅]^𝑇, the vector of reaction rates (the reaction kinetics) is Φ(𝜉)=[𝜇^𝑜_𝑠𝑋𝜇^𝑟_𝑠𝑋𝜇^𝑜_𝑒𝑋]^𝑇; 𝐾 is the matrix of yield coefficients, 𝐹=[0𝐷𝑆_in0OTR]^𝑇 is the vector of feeding rates, and 𝑄=[0000CTR]^𝑇 is the vector of rates of removal of the components in gaseous form.

The most difficult task for the construction of the dynamical model is the modelling of the reaction kinetics. The form of kinetics is complex, nonlinear, and in many cases unknown. The reaction rates can be expressed asΦ(𝜉,𝑡)=𝐻(𝜉)⋅𝜇(𝜉,𝑡), (2.11) where 𝐻(𝜉) is a state-dependent diagonal matrix, whose elements correspond to the reactions’ reactants.

The specific growth rates 𝜇(𝜉,𝑡) are considered as a vector of time-varying parameters. For the baker’s yeast bioprocess the form of matrix 𝐻(𝜉) is very simple:𝐻(𝜉)=𝑋000𝑋000𝑋. (2.12)

Then the model (2.10) becomes𝑑𝜉(𝜉)()𝑑𝑡=𝐾⋅𝐻⋅𝜇𝜉,𝑡−𝐷⋅𝜉+𝐹−𝑄. (2.13)

For the reaction rates of the baker’s yeast bioprocess there are many possible models [21, 22], and so forth. However, the process of yeast growth on glucose with ethanol production is generally described by the following three metabolic reactions, widely analyzed by Sonnleitner and Kappeli [20].

First, the reaction rate of the respiratory growth on glucose and the associated specific growth rate are𝜑₁(𝜉)=𝜇₁(𝜉)𝜇⋅𝑋,₁(𝜉)=min𝜇_𝑠,𝜇_𝑐max𝑘₅, (2.14) where the kinetic terms associated with the glucose consumption 𝜇_𝑠 and with the respiratory capacity 𝜇_𝑐max are modelled by Monod-type laws: 𝜇_𝑠=𝑞_𝑠max𝑆/(𝑆+𝐾_𝑠), 𝜇_𝑐max=𝑞_𝑐max𝐶/(𝐶+𝐾_𝑐), with 𝑞_𝑠max and 𝑞_𝑐max the maximal values of the specific growth rates of glucose and oxygen, and 𝐾_𝑠 and 𝐾_𝑐 the saturation parameters for glucose and oxygen uptake, respectively.

Second, the reaction rate of the fermentative growth on glucose and the associated specific growth rate are𝜑₂(𝜉)=𝜇₂(𝜉)𝜇⋅𝑋,₂(𝜉)=max0,𝜇_𝑠−𝜇_𝑐max𝑘₅. (2.15)

Finally, if the oxidation capacity is sufficiently high to oxidize both ethanol and glucose, then their consumption is possible. Then, the reaction rate of the respiratory growth on ethanol and the associated specific growth rate are𝜑₃(𝜉)=𝜇₃(𝜉)𝜇⋅𝑋,₃(𝜉)=max0,min𝜇_𝑒,𝜇_𝑐max−𝑘₅𝜇_𝑠𝑘₆, (2.16) where 𝜇_𝑒 is the potential ethanol oxidative rate, modelled by Monod-type law: 𝜇_𝑒=𝑞_𝑒max𝐸/(𝐸+𝐾_𝑒), with 𝑞_𝑒max the maximal value of the specific growth rates of ethanol, and 𝐾_𝑒 the saturation parameter for growth on ethanol.

The above kinetic model is based on the well-known bottleneck hypothesis developed by Sonnleitner and Kappeli [20]. This model assumes a limited capacity of yeast, leading to the production of ethanol under conditions of oxygen limitation and/or high glucose concentration. The glucose concentration that corresponds to the oxidative capacity is reached when 𝜇_𝑠=𝜇_𝑐max/𝑘₅, and it is called critical concentration. Depending on this critical concentration, two operating regimes can be observed. At glucose concentration lower than its critical value, the system is in the respirative regime. The glucose consumption rate is smaller than the maximal respiratory capacity 𝜇_𝑠≤𝜇_𝑐max/𝑘₅ and the rate of the oxidative glucose metabolism is given by the glucose consumption rate. Ethanol can be oxidized in parallel with glucose and the rate of oxidative ethanol metabolism is determined by the excess of respiratory capacity and the available ethanol. At glucose concentration higher than its critical value, the system is in the so-called respirofermentative regime. Then the glucose consumption rate is larger than the maximal respiratory capacity 𝜇_𝑠>𝜇_𝑐max/𝑘₅ and the respiratory capacity establishes the rate of the oxidative glucose metabolism. The excess of glucose is processed by the fermentative metabolism. Under lack of oxygen conditions, the fermentative pathway is predominating.

Since there are three unknown kinetics, the implementation of some state observers (such as asymptotic observers) require the online measurements of three state variables. Also, the online estimation algorithms for kinetic parameters (e.g., observer-based techniques) necessitate the knowledge of minimum three states. From the technological point of view, 𝐸, 𝐶, and 𝐺 are the online available state variables, so it is important to design the estimators by using these measurements. Unfortunately, according to the yield coefficient values, 𝐸, 𝐶, and 𝐺 are linearly dependent, with the corresponding yield matrix being ill conditioned [23]. This would result in sensitive algorithms to numerical errors, thus having performances degraded. To solve this problem, Pomerleau and Perrier [23] suggested a reformulation of model (2.10). This reformulation is based on the division of the complete process model (2.10) into two partial models, corresponding to the above-described mechanisms: the respirofermentative partial model for the ethanol production state of the process and the respirative partial model, corresponding to the ethanol consumption state of the process.

The fermentative partial model is stated as𝜉̇₁𝜉̇₂𝜉̇₃𝜉̇₄𝜉̇₅=11−𝑘₁−𝑘₂0𝑘₃−𝑘₅0𝑘₇𝑘₈𝜇₁𝜇₂𝜉₁𝜉−𝐷¹𝜉₂𝜉₃𝜉₄𝜉₅+0𝐷𝑆_in00−000OTR0CTR. (2.17)

The respirative partial model is stated as𝜉̇₁𝜉̇₂𝜉̇₃𝜉̇₄𝜉̇₅=11−𝑘₁00−𝑘₄−𝑘₅−𝑘₆𝑘₇𝑘₉𝜇₁𝜇₃𝜉₁𝜉−𝐷¹𝜉₂𝜉₃𝜉₄𝜉₅+0𝐷𝑆_in00−000OTR0CTR. (2.18)

3. Distribution-Based Identification of Yield Coefficients

3.1. Problem Statement

In this approach the set of nonlinear differential equations describing the state evolution is mapped into a set of linear algebraic equations with respect to the model parameters. Using techniques utilized in distribution approach, the measurable functions and their derivatives are represented by functionals on the fundamental space of testing functions [14]. The main advantages of this method are that a set of algebraic equations with real coefficients results and the formulations are free from boundary conditions. The idea of utilizing test functions in system identification was proposed by Pearson and Lee [13] in terms of modulating functions in order to ameliorate the noise handling for deterministic least-squares identification based on time-limited data.

Let us denote by Ω_𝑛 the fundamental space from the distribution theory [24] of the real functions 𝜔∶ℜ→ℜ,𝑡→𝜔(𝑡) with compact support 𝑇, having continuous derivatives at least up to the order 𝑛. Let 𝑞∶ℜ→ℜ,𝑡→𝑞(𝑡) be a function which admits a Riemann integral on any compact interval 𝑇 from ℜ. Using this function, a unique distribution (or generalized function)𝐹_𝑞∶Ω_𝑛⟶ℜ,𝜔⟶𝐹_𝑞(𝜔)∈ℜ (3.1) can be built by the relation:𝐹_𝑞(𝜔)=_ℜ𝑞(𝑡)𝜔(𝑡)𝑑𝑡,∀𝜔∈Ω_𝑛. (3.2)

In distribution theory, the notion of 𝑘-order derivative is introduced. If 𝐹_𝑞∈Ω_𝑛, then its 𝑘-order derivative is a new distribution 𝐹_𝑞^(𝑘)∈Ω_𝑛 uniquely defined by the relations:𝐹_𝑞^(𝑘)(𝜔)=(−1)^𝑘𝐹_𝑞𝜔^(𝑘),∀𝜔∈Ω_𝑛,𝜔⟶𝐹_𝑞^(𝑘)(𝜔)=(−1)^𝑘_ℜ𝑞(𝑡)𝜔^(𝑘)(𝑡)𝑑𝑡∈ℜ, (3.3) where𝜔^(𝑘)∶ℜ⟶ℜ,𝑡⟶𝜔^(𝑘)(𝑡)=𝑑^𝑘𝜔(𝑡)𝑑𝑡^𝑘 (3.4) is the 𝑘-order time derivative of the test function.

When 𝑞∈𝐶^𝑘(ℜ), then𝐹_𝑞^(𝑘)(𝜔)=_ℜ𝑞^(𝑘)(𝑡)𝜔(𝑡)𝑑𝑡=(−1)^𝑘_ℜ𝑞(𝑡)𝜔^(𝑘)(𝑡)𝑑𝑡 (3.5) that means the 𝑘-order derivative of a distribution generated by a function 𝑞∈𝐶^𝑘(ℜ) equals the distribution generated by the 𝑘-order time derivative of the function 𝑞. So, in place of the states and their time derivatives of a system one utilize the corresponding distributions and, in some particular cases, it is possible to obtain a system of equations linear in parameters. If the system is compatible the model parameters are structurally identifiable.

In our study, it were utilized three types of test functions characterized by a bounded support 𝑇=(𝑡_𝑎𝑡_𝑏),𝑡_𝑎<𝑡_𝑏, all of these accomplishing the condition:𝜔(𝑡)=0;∀𝑡∈−∞,𝑡_𝑎∪𝑡_𝑏,+∞. (3.6)

The nonzero restriction is one of the following three types.(1)Exponential: 𝜔(𝑡)𝑡=exp^𝑎𝑡_𝑏𝑡−𝑡_𝑎𝑡−𝑡_𝑏. (3.7)(2)Sinusoidal: 𝜔(𝑡)=sin^𝑝𝑡−𝑡_𝑏𝜋𝑡_𝑏−𝑡_𝑎. (3.8)(3)Polynomial: 𝜔(𝑡)=𝑡−𝑡_𝑎^𝑝𝑡−𝑡_𝑏^𝑝, (3.9) where 𝑝≥2 is an integer.

Figures 1 and 2 present the exponential type test function and its first-order derivative for 𝑇=[0.10.9].

Figure 1 

Exponential type test function.

Figure 2 

First-order derivative of exponential-type test function.

3.2. The Hierarchical Structure of the Identification Algorithm

Consider all state variables accessible for measurements. This assumption can be considered unrealistic but the unknown states can be estimated using different approaches. From the multitude of proposed solutions one can mention methods based on balance equations, observers (asymptotic observers, high-gain observers) [5], Kalman filters [25], neural networks [26], and fuzzy reasoning [27].

The dynamical systems (2.17) and (2.18) contain rational dependences between parameters and measured variables. To obtain linear equations in unknown parameters, the identification problem is split in several simpler problems [28]. Based on the specific structure of these systems, it is possible to group the state equations in such way to determine some interconnected relations. They are organized in a hierarchical structure. In the first stage some state equations are utilized to obtain a set of linear equations in some parameters. If these parameters are identifiable then they can be used as known parameters in the following stages. This process is repeated in the other stages until all the parameters are identified.

Because the dilution rate 𝐷 can be externally modified, it will be considered the first component of the input vector and the other component of 𝑢 is the concentration 𝑆_in∶𝑢=[𝑢₁𝑢₂]=[𝐷𝑆_in].

From the third and fourth equations of the fermentative partial model (2.17) one gets the following relations for the reaction rates:𝜇₁𝜉₁=1𝑘₅−𝜉̇₄−𝐷𝜉₄,𝜇+OTR₂𝜉₁=1𝑘₃𝜉̇₃+𝐷𝜉₃. (3.10)

These relations are replaced in the other three equations to obtain linear equations with respect to yield coefficients. So, one can split the identification problem in the following steps.

Step 1 (identification of 𝑘₃ and 𝑘₅). Replacing (3.10) in the first equation of the system (2.17) one get 𝜉̇₁+𝐷𝜉₁=𝜃₃𝜉̇₃+𝐷𝜉₃+𝜃₅−𝜉̇₄−𝐷𝜉₄+OTR, (3.11) where 𝜃₃=1/𝑘₃ and 𝜃₅=1/𝑘₅.
One denotes 𝑣(𝑡)=𝜉̇₁(𝑡)+𝐷𝜉₁(𝑡),𝑤₁(𝑡)=𝜉̇₃(𝑡)+𝐷𝜉₃(𝑡),𝑤₂(𝑡)=−𝜉̇₄(𝑡)−𝐷𝜉₄(𝑡)+OTR. (3.12)
Multiplying these functions with a test function 𝜔_𝑖(𝑡) and integrating over ℜ one gets the following distributions: 𝐹^𝑖_𝑣𝜔_𝑖=_ℜ𝑣(𝑡)⋅𝜔_𝑖(𝑡)=𝑑𝑡_ℜ𝜉̇₁(𝑡)⋅𝜔_𝑖(𝑡)𝑑𝑡+_ℜ𝐷𝜉₁(𝑡)⋅𝜔_𝑖(𝑡)𝑑𝑡=_ℜ𝜉₁(𝑡)⋅𝜔_𝑖⁽¹⁾(𝑡)𝑑𝑡+_ℜ𝐷𝜉₁(𝑡)⋅𝜔_𝑖(𝑡)𝐹𝑑𝑡,_𝑖𝑤1𝜑_𝑖=_ℜ𝑤₁(𝑡)⋅𝜔_𝑖(𝑡)=𝑑𝑡_ℜ𝜉̇₃(𝑡)⋅𝜔_𝑖(𝑡)𝑑𝑡+_ℜ𝐷𝜉₃(𝑡)⋅𝜔_𝑖(𝑡)𝑑𝑡=−_ℜ𝜉₃(𝑡)⋅𝜔_𝑖⁽¹⁾(𝑡)𝑑𝑡+_ℜ𝐷𝜉₃(𝑡)⋅𝜔_𝑖(𝑡)𝐹𝑑𝑡,_𝑖𝑤2𝜔_𝑖=_ℜ𝑤₂(𝑡)⋅𝜔_𝑖(𝑡)𝑑𝑡=_ℜ−𝜉̇₄(𝑡)⋅𝜔_𝑖(𝑡)𝑑𝑡+_ℜ−𝐷𝜉₄(𝑡)+OTR⋅𝜔_𝑖(𝑡)=𝑑𝑡_ℜ𝜉₄(𝑡)⋅𝜔_𝑖⁽¹⁾(𝑡)𝑑𝑡+_ℜ−𝐷𝜉₄(𝑡)+OTR⋅𝜔_𝑖(𝑡)𝑑𝑡. (3.13)(3.14)
As one must identify two parameters it is enough to choose a finite number 𝑁≥2 of test functions 𝜔_𝑖,𝑖=1∶𝑁 and to build an algebraic system of equations, 𝐹_𝑤⋅𝜃=𝐹_𝑣, (3.15) where, 𝜃=[𝜃₃𝜃₅]^𝑇 is the vector of unknown parameters, 𝐹_𝑤 is a 𝑁×2 matrix of real numbers: 𝐹_𝑤=𝐹¹_𝑤1(𝜔₁),…,𝐹^𝑖_𝑤1(𝜔_𝑖),…,𝐹^𝑁_𝑤1(𝜔_𝑁)𝐹¹_𝑤2(𝜔₁),…,𝐹^𝑖_𝑤2(𝜔_𝑖),…,𝐹^𝑁_𝑤2(𝜔_𝑁)^𝑇 (3.16) and 𝐹_𝑣 denotes an 𝑁-column real vector build from distributions calculated as in (3.13): 𝐹_𝑣=𝐹¹_𝑣(𝜔₁),…,𝐹^𝑖_𝑣(𝜔_𝑖),…,𝐹^𝑁_𝑣(𝜔_𝑁)^𝑇. (3.17) If rank(𝐹_𝑤)=2 then a unique solution is obtained: 𝜃=𝐹^𝑇_𝑤⋅𝐹_𝑤⁻¹𝐹^𝑇_𝑤⋅𝐹_𝑣. (3.18) One obtains 𝑘₃=1/𝜃₃ and 𝑘₅=1/𝜃₅.

Step 2 (identification of 𝑘₁ and 𝑘₂). Considering known 𝑘₃=𝑘₃ and 𝑘₅=𝑘₅ from Step 1 and substituting (3.10) in the second equation of (2.17) one gets 𝜉̇₂+𝐷𝜉₂−𝐷𝑆_in=𝜃₁1𝑘₅−𝜉̇₄−𝐷𝜉₄+OTR+𝜃₂1𝑘₃𝜉̇₃+𝐷𝜉₃, (3.19) where 𝜃₁=−𝑘₁ and 𝜃₂=−𝑘₂.

The procedure of identification from Step 1 is repeated and one obtains 𝑘₁=−𝜃₁ and 𝑘₂=−𝜃₂.

Step 3 (identification of 𝑘₇ and 𝑘₈). Considering known 𝑘₃=𝑘₃ and 𝑘₅=𝑘₅ from Step 1 and substituting (3.10) in the last equation of (2.17) one gets 𝜉̇₅+𝐷𝜉₅+CTR=𝜃₇1𝑘₃𝜉̇₃+𝐷𝜉₃+𝜃₈1𝑘₅−𝜉̇₄−𝐷𝜉₄+OTR, (3.20) where 𝜃₇=𝑘₇ and 𝜃₈=𝑘₈.

The procedure of identification from Step 1 is repeated and one obtains 𝑘₇=𝜃₇ and 𝑘₈=𝜃₈.

Since the respirative partial model (2.18) has a similar structure with the fermentative partial model (2.17), in a similar manner one obtains the parameters 𝑘₄, 𝑘₆, and 𝑘₉.

4. Analysis of the Algorithm Properties

In the following one presents some consistency and numerical aspects related to the presented algorithm and a short overview of the classical identification procedure widely used for this kind of biotechnological systems.

4.1. Identifiability

Identifiability is a necessary prerequisite for model identification; it concerns uniqueness of the model parameters determined from the input-output data, under ideal conditions of noise-free observations and error-free model structure. A remarkable feature of distribution-based identification procedure is that it provides a linear reparameterization of the input-output relation of the nonlinear system. This reparametrization of the system involves a very simple identifiability condition to be accomplished, that is, the existence of the matrix (𝐹^𝑇_𝑤⋅𝐹_𝑤)⁻¹ or, equivalently, 𝐹_𝑤 is of full rank.

4.2. Consistency

Obviously, consistency of the estimates is directly influenced by the precision of numerical integration used to compute the value of the distributions. There are available numerous integration methods (often called numerical quadrature) with various degrees of precision. One presents shortly the techniques used in the simulations to the approximate calculation of a definite integral 𝐼_𝑓=∫^𝑏_𝑎𝑓(𝑥)𝑑𝑥 where 𝑓(𝑥) is a given function and [𝑎,𝑏] a finite interval. Interpolatory quadrature formulas, where the nodes are constrained to be equally spaced, are called Newton-Cotes formulas. These are especially suited for integrating a tabulated function (such is our case). Newton-Cotes numerical integration rule 𝐼_𝑓 is a weighted sum of function values: 𝐼_𝑓=∑^𝑁−1_𝑖=0𝑐_𝑖⋅𝑓(𝑥_𝑖), where coefficients {𝑐_𝑖} are derived from interpolating polynomial fitted for points {𝑥_𝑖,𝑓_𝑖}. A Newton-Cotes formula of any degree 𝑛 can be constructed. One of the simplest integration methods is the trapezoidal rule. The trapezoidal rule is based on linear interpolation of 𝑓(𝑥) at 𝑥₁=𝑎 and 𝑥₂=𝑏, that is, 𝑓(𝑥) is approximated by𝑝(𝑥)(𝑎)+()𝑓(𝑏)(𝑎)=𝑓𝑥−𝑎−𝑓.𝑏−𝑎 (4.1) The integral of 𝑝(𝑥) equals the of trapezoid with base (𝑏−𝑎) times the average height12𝑓(𝑎)+𝑓(𝑏) (4.2) hence^𝑏_𝑎𝑓(𝑥)𝑑𝑥=𝑏−𝑎2𝑓(𝑎)+𝑓(𝑏). (4.3) To increase the accuracy one subdivides the interval [𝑎,𝑏] and assumes that 𝑓_𝑖=𝑓(𝑥_𝑖) is known on a grid of equidistant points 𝑥₀=𝑎,𝑥_𝑖=𝑥₀+𝑖ℎ,𝑥_𝑛=𝑏, where ℎ=(𝑏−𝑎)/𝑛 is the step length. The trapezoidal approximation for the 𝑖th subinterval is^𝑥𝑖+1_𝑥𝑖𝑓(𝑥)ℎ𝑑𝑥=2𝑓_𝑖+𝑓_𝑖+1+𝑅_𝑖, (4.4) where𝑅_𝑖ℎ=−³𝑓12^′′𝜍_𝑖,𝜍_𝑖∈𝑥_𝑖,𝑥_𝑖+1. (4.5) Summing the contributions for each subinterval [𝑥_𝑖,𝑥_𝑖+1],𝑖=0∶𝑛, gives^𝑏_𝑎𝑓(𝑥)ℎ𝑑𝑥=2𝑓₀+𝑓_𝑛+ℎ^{𝑖=𝑛−1}_𝑖=2𝑓_𝑖+𝑅_𝑇. (4.6) The global truncation error is𝑅_𝑇ℎ=−³12^{𝑖=𝑛−1𝑖=0}𝑓^′′𝜍_𝑖ℎ=−²()𝑓12𝑏−𝑎^′′(𝜍),𝜍∈[𝑎,𝑏]. (4.7) This shows that by choosing ℎ small enough we can make the truncation error arbitrarily small. In other words, we have asymptotic convergence when ℎ→0.

In the composite Simpson’s rule one divides the interval [𝑎,𝑏] into an even number 𝑛=2 m steps of length ℎ and uses the formula^𝑏_𝑎𝑓(𝑥)ℎ𝑑𝑥=3𝑓₀+4𝑆₁+2𝑆₂+𝑓_𝑛+𝑅_𝑇, (4.8) where𝑆₁=𝑓₁+𝑓₃+⋯+𝑓_𝑛−1,𝑆₂=𝑓₂+𝑓₄+⋯+𝑓_𝑛−2 (4.9) are sums over odd and even indices, respectively. The remainder is𝑅_𝑇=ℎ⁴()𝑓180𝑏−𝑎^IV(𝜍),𝜍∈[𝑎,𝑏]. (4.10)

This shows that one have gained two orders of accuracy compared to the trapezoidal rule, without using more function evaluations.

Let us study properties of Newton-Cotes formulas in frequency domain (or spectral properties). Newton-Cotes rules are symmetric (hence linear phase) digital filters with finite impulse response. One of the most important characteristic of digital filter is magnitude/frequency response-function which shows how much filter damps or amplifies magnitude of particular frequency contained in input data.

So, for trapezoidal rule and Simpson’s rule one obtains𝐻(𝑧)=ℎ2⋅1+𝑧⁻¹1−𝑧⁻¹(𝑧)=ℎ,𝐻3⋅1+4𝑧⁻¹+𝑧⁻²1−𝑧⁻¹,respectively. (4.11)

These transfer functions have the amplitude-frequency Bode characteristics from Figure 3.

(a)

(b)

Figure 3 

Amplitude-frequency bode characteristics: (a) trapezoidal rule; (b) Simpson’s rule.

As one can see classical Newton-Cotes formulas suppress high frequencies (noise) in the input data. In that sense they can be considered low-pass filters. From Figure 3 one observes that trapezoidal rule offers a better suppression of noise than the Simpson’s rule, so in the algorithm implementation one used the trapezoidal rule for numerical integration.

4.3. “Classical” Identification Procedure

The “classical” identification approach of yield coefficients is a two-step procedure under the assumption that full state measurements are available [1]. This method is based on a state transformation that allows reformulating the dynamical model into separate submodels. The first submodel depends only on the reaction structure and is independent of the kinetics. It can be linearly reparametrized and used for the identification of the yield coefficients by means of linear regressions, provided suitable identifiability conditions are satisfied. We present briefly this method for estimating the yield coefficients and we use it for comparison in the next section.

The general dynamical model given in (2.1) represents a particular class of nonlinear state-space models. The nonlinearity lies in the reaction rates 𝜑_𝑖(𝜉) that are (nonlinear) functions of the state variables. These functions enter the model in the form 𝐾𝜑_𝑖(𝜉) (where 𝐾 is a constant matrix), which is a set of linear combinations of the same nonlinear functions 𝜑_𝑖(𝜉). This particular feature can be exploited to separate the nonlinear part from the linear part of the model by an adequate linear state transformation. More precisely, one chooses a nonsingular partition𝐾_𝑎𝐾_𝑏=𝑇𝐾, (4.12) with 𝐾_𝑎∈𝑅^𝑝×𝑚 of full row rank matrix (i.e., 𝑝=rank(𝐾)),𝐾_𝑏∈𝑅^{(𝑛−𝑝)×𝑚} and 𝑇 a permutation matrix. The induced partitions of the vectors 𝜉 and 𝑢 are𝜉_𝑎𝜉_𝑏𝑢=𝑇𝜉^𝑎𝑢_𝑏=𝑇𝑢. (4.13)Model (2.1) is then partitioned into two submodels:𝑑𝜉_𝑎𝜉𝑑𝑡=𝐾⋅𝜑_𝑎,𝑡−𝐷𝜉_𝑎+𝑢_𝑎,𝑑𝜉_𝑏𝜉𝑑𝑡=𝐾⋅𝜑_𝑏,𝑡−𝐷𝜉_𝑏+𝑢_𝑏. (4.14) Then with the state transformation𝜉_𝑎=𝜉_𝑎,𝑧=𝐶𝜉_𝑎+𝜉_𝑏, (4.15) one transform the initial model into𝜉̇_𝑎=𝐾_𝑎𝜑𝜉_𝑎,𝑧−𝐶𝜉_𝑎−𝐷𝜉_𝑎+𝑢_𝑎,𝑧̇=−𝐷𝑧+𝐶𝑢_𝑎+𝑢_𝑏, (4.16) where the (𝑛−𝑝)×𝑝 matrix 𝐶 is the unique solution of𝐶𝐾_𝑎+𝐾_𝑏=0 (4.17) that is,𝐶=−𝐾_𝑏𝐾^∗_𝑎, (4.18) where 𝐾^∗_𝑎 is a generalized or pseudoinverse of 𝐾_𝑎. This subsystem (4.16) can be augmented with an equation derived from (2.7) as follows:𝑧̇=−𝐷𝑧+𝐶𝑢_𝑎+𝑢_𝑏,𝜉_𝑏=𝑧−𝐶𝜉_𝑎. (4.19)

It can be considered as a linear time-varying (if 𝑑 varies in the course of time) model with state z, input (𝜉_𝑎,𝑢_𝑎,𝑢_𝑏), and output 𝜉_𝑏. It is nonlinearly parametrized by the yield coefficients but linearly reparametrized by the nonzero entries of 𝐶. When data of the signals 𝜉_𝑎, 𝑢_𝑎, 𝑢_𝑏, and 𝜉_𝑏 are available, the auxiliary model (4.19) can be used to identify the yield coefficients independently of the knowledge of the reaction rates. The model (4.19) can be used to perform the identification of the nonzero entries of C by a linear regression technique, with the yield coefficients 𝑘_𝑖 recovered afterwards from (4.17). The identification makes sense only if the yield coefficients 𝑘_𝑖 are structurally uniquely identifiable with the auxiliary model (4.19).

5. Simulation Results

The efficacy of our distribution-based approach is shown by simulations, comparing it with this “classical” approach reviewed in the previous section. Also, the influence of sampling period, test functions type, initial conditions, and input type is analysed. To compare statistical performances of the different approaches and simulation conditions the empirical normalized mean square error (NMSE) was used, which is defined as1NMSE=𝑁^𝑁𝑗=1NMSE𝜃_𝑗,withNMSE𝜃_𝑗=𝜃_𝑗−𝜃^∗_𝑗𝜃^∗_𝑗², (5.1) where 𝑁 is the number of estimated parameters, 𝜃_𝑗 is the 𝑗th element of the estimated parameter vector, while the “*” superscript denotes the true value of the parameter.

The model given by (2.17) and (2.18) was integrated using a fourth-order Runge-Kutta routine with a sampling period of 1 minute.

In order to study the sensitivity of the estimation method to the sampling period 𝑇_𝑠, to the initial conditions, to the test functions type, to the input type, and to the noise, the following parameters were used:(a)sampling period: 𝑇_𝑠∈{1min,10min,20min},(b)input signals: 𝑢∈{0.1+0.05∗sin(𝑡),0.1+0.005∗𝑡,0.15∗exp(9∗10⁻³𝑡)},(c)initial conditions: IC₁: [1 0.2 0.15 0.0066 0.008] IC₂: [0.1 0.3 0 0.002 0.01 2] IC₃: [0.5 0.5 3 0.01 0.001 1],(d)test functions: Type1: exponential Type2: sinusoidal Type3: polynomial,(e)noise: {noise-free, zero-mean white (SNR = 50 dB, 40 dB, 30 dB, 20 dB, 10 dB)}.

Some of the signals used in the identification algorithm are presented in Figures 4–6. Figure 4 presents the noise-free system time response for the fermentative partial model (2.17). The measured states are represented in percent with respect to their maximal values. Three test function of sinusoidal type are presented in Figure 5. In Figure 6 some noisy measurements are presented.

Figure 4 

Time response for noise-free system.

(a)

(b)

Figure 5 

Test functions (a) and their derivatives (b) used for identification.

Figure 6 

Noisy measurements of state variable 𝐸: SNR = 20 dB.

The simulation results presented in Tables 1–5 leads to the following remarks.(1)Sensitivity to the sampling period: results given in Table 1 show that the performance of the algorithm deteriorates with the increasing sampling period.(2)Sensitivity to the inputs type: results given in Table 2 conclude that the algorithm is not too sensitive to the type of input. However, as was expected, the best performances are obtained for a sinusoidal input type.(3)Sensitivity to the initial conditions and test function type: as shown in Tables 3 and 4, initial conditions and type of the test functions have very little influence on the accuracy of the estimates.(4)Sensitivity to the noise: Tables 5 and 6 present the simulation results for distribution-based method and for “classical” method, respectively. The results are still very good in the distribution-based case even for the high level of noise (SNR = 10 dB), while the “classical” approach is less robust to increasing noise level.

Table 1 

Influence of the sampling period.

Table 2 

Influence of the input type.

Table 3 

Influence of the initial conditions.

Table 4 

Influence of the test functions type.

Table 5 

Influence of the noise (distribution-based method).

Table 6 

Influence of the noise (“classical” method).

6. Conclusions

The paper presents a distribution-based identification procedure for offline estimation of yield coefficients in a baker’s yeast bioprocess. This procedure is a functional type method, which transforms a differential system of equations to an algebraic system in unknown parameters. The relation between the state variables of the system is represented by functionals using techniques from distribution theory based on test functions from a finite-dimensional fundamental space. The identification algorithm has a hierarchical structure, which allows obtaining a linear algebraic system of equations in the unknown parameters. The coefficients of this algebraic system are functionals depending on the input and state variables and are evaluated through some testing functions from distribution theory.

According to the proposed procedure, first, only some state equations are evaluated throughout testing functions to obtain a set of linear equations in some parameters. The results of this first stage of identification are utilized for expressing other parameters by linear equations. This process is repeated until all parameters are identified.

The influence of the sampling period, initial conditions, test functions type, input type, and noise on the parameters estimates was empirically analyzed.

The proposed method presents two main advantages: there is necessary no information about the parameters of the reaction rates to identify the yield coefficients of the baker’s yeast bioprocess and the algorithm provides very good results even the measurements are noise-contaminated because the evaluation of states derivatives is completely avoided.

Nomenclature

𝐶:	Dissolved oxygen concentration (g/L)
𝐶∗:	Equilibrium concentration of the dissolved oxygen (g/L)
CTR:	Carbon dioxide transfer rate (h−1g/L)
𝐷:	Dilution rate (h−1)
𝐸:	Ethanol concentration (g/L)
𝐹:	Vector of feeding rates
𝐹in:	Input feed rate (Lh−1)
𝐺:	Dissolved carbon dioxide concentration (g/L)
𝐾:	Matrix of yield coefficients
𝑘𝑖,𝑖=1,9:	Yield coefficients
𝐾𝐿𝑎:	Mass transfer coefficient (h−1)
𝐾𝑒:	Saturation parameter for growth on ethanol (g/L)
𝐾𝑠, 𝐾𝑐:	Saturation parameters for glucose/oxygen uptake, respectively (g/L)
𝐾𝑉:	Transfer coefficient
OTR:	Oxygen transfer rate (h−1g/L)
𝑄:	Vector of rates of removal of the components in gaseous form
𝑞𝑐max:	Maximal oxygen specific growth rate (gO2/(gbiomass⋅h))
𝑞𝑒max:	Maximal ethanol specific growth rate (gethanol/(gbiomass⋅h))
𝑞𝑠max:	Maximal glucose specific growth rate (gglucose/(gbiomass⋅h))
𝑆:	Glucose concentration (g/L)
𝑆in:	Glucose concentration on the feed (g/L)
𝑉:	Culture volume in the reactor (L)
𝑋:	Biomass concentration (g/L)
Φ:	Vector of reaction rates (reaction kinetics)
𝜇𝑜𝑠,𝜇𝑟𝑠, 𝜇𝑜𝑒:	Specific growth rates of baker’s yeast bioprocess (h−1)
𝜇:	Vector of specific growth rates
𝜇𝑐max:	Kinetic term associated with the respiratory capacity
𝜇𝑒:	Potential ethanol oxidative rate
𝜇𝑠:	Kinetic term associated with the glucose consumption
𝜉:	State vector
Ω𝑛:	The fundamental space from the distribution theory
ℜ:	The set of real numbers
𝜔(𝑡):	Test function
𝐹𝑞:	Distribution (generalized function)
𝜃:	The vector of unknown parameters
𝜃𝑖,𝑖=1,9:	The unknown parameters.

Acknowledgment

This work was supported by CNCSIS-UEFISCSU, Project number PN II-RU TE 106, 9/04.08.2010.