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

Authors: Dorin Sendrescu
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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.


---

## Body

## 1. Introduction

A process carried out in a bioreactor can be defined as a set ofm biochemical reactions involving n components (with n≥m). The global dynamics can be represented by following the dynamical state-space model [1]:(1.1)dξdt=K⋅φ(ξ,t)-Dξ+F-Q,
where ξ∈ℜn×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 K·φ(ξ,t) is the rate of consumption and/or production of the components in the reactor, that is, the reaction kinetics and the term -Dξ+F-Q 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 Section2. 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)S+C⟶μsoX+G,(2.2)S⟶μsrX+E+G,(2.3)E+C⟶μeoX+G.In the reaction schemes (2.1)–(2.3), S is the glucose, X the biomass, C the dissolved oxygen, E the ethanol, and G 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. μso,μsr, and μeo 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]:(2.4)dXdt=(μso+μsr+μeo)X-DX,(2.5)dSdt=(-k1μso-k2μsr)X+D(Sin-S),(2.6)dEdt=(k3μsr-k4μeo)X-DE,(2.7)dCdt=(-k5μso-k6μeo)X-DC+OTR,(2.8)dGdt=(k7μso+k8μsr+k9μeo)X-DG+CTR,(2.9)dVdt=Fin=DV.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 V in the bioreactor. In (2.4)–(2.9), ki,i=1,9¯ are the positive yield coefficients, Sin is the glucose concentration on the feed, OTR is the oxygen rate, defined as OTR=KLa(C*-C) (where KLa is the mass transfer coefficient and C* the equilibrium concentration of the dissolved oxygen), CTR is the carbon dioxide transfer rate (defined as CTR=KVKLaG, with KV a transfer coefficient), Fin is the input feed rate, and D=Fin/V is the so-called dilution rate.The model (2.4)–(2.8) can be expressed in the matrix form:(2.10)[ẊṠĖĊĠ]=[111-k1-k200k3-k4-k50-k6k7k8k9]︸K[μsoμsrμeo]X-D[XSECG]+[0DSin0OTR0]-[0000CTR].The next notations will be used: the state vector isξ=[XSECG]T=[ξ1ξ2ξ3ξ4ξ5]T, the vector of reaction rates (the reaction kinetics) is Φ(ξ)=[μsoXμsrXμeoX]T; K is the matrix of yield coefficients, F=[0DSin0OTR]T is the vector of feeding rates, and Q=[0000CTR]T 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)Φ(ξ,t)=H(ξ)⋅μ(ξ,t),
where H(ξ) is a state-dependent diagonal matrix, whose elements correspond to the reactions’ reactants.The specific growth ratesμ(ξ,t) are considered as a vector of time-varying parameters. For the baker’s yeast bioprocess the form of matrix H(ξ) is very simple:(2.12)H(ξ)=[X000X000X].Then the model (2.10) becomes(2.13)dξdt=K⋅H(ξ)⋅μ(ξ,t)-D⋅ξ+F-Q.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(2.14)φ1(ξ)=μ1(ξ)⋅X,μ1(ξ)=min⁡(μs,μcmax⁡k5),
where the kinetic terms associated with the glucose consumption μs and with the respiratory capacity μcmax are modelled by Monod-type laws: μs=qsmaxS/(S+Ks), μcmax=qcmaxC/(C+Kc), with qsmax and qcmax the maximal values of the specific growth rates of glucose and oxygen, and Ks and Kc 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(2.15)φ2(ξ)=μ2(ξ)⋅X,μ2(ξ)=max⁡(0,μs-μcmax⁡k5).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(2.16)φ3(ξ)=μ3(ξ)⋅X,μ3(ξ)=max⁡(0,min⁡(μe,μcmax⁡-k5μsk6)),
where μe is the potential ethanol oxidative rate, modelled by Monod-type law: μe=qemaxE/(E+Ke), with qemax the maximal value of the specific growth rates of ethanol, and Ke 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 μs=μcmax/k5, 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 μs≤μcmax/k5 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 μs>μcmax/k5 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,E, C, and G 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, E, C, and G 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(2.17)[ξ̇1ξ̇2ξ̇3ξ̇4ξ̇5]=[11-k1-k20k3-k50k7k8][μ1μ2]ξ1-D[ξ1ξ2ξ3ξ4ξ5]+[0DSin0OTR0]-[0000CTR].The respirative partial model is stated as(2.18)[ξ̇1ξ̇2ξ̇3ξ̇4ξ̇5]=[11-k100-k4-k5-k6k7k9][μ1μ3]ξ1-D[ξ1ξ2ξ3ξ4ξ5]+[0DSin0OTR0]-[0000CTR].


## 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Ωn the fundamental space from the distribution theory [24] of the real functions ω:ℜ→ℜ,t→ω(t) with compact support T, having continuous derivatives at least up to the order n. Let q:ℜ→ℜ,t→q(t) be a function which admits a Riemann integral on any compact interval T from ℜ. Using this function, a unique distribution (or generalized function)(3.1)Fq:Ωn⟶R,ω⟶Fq(ω)∈R
can be built by the relation:(3.2)Fq(ω)=∫Rq(t)ω(t)dt,∀ω∈Ωn.In distribution theory, the notion ofk-order derivative is introduced. If Fq∈Ωn, then its k-order derivative is a new distribution Fq(k)∈Ωn uniquely defined by the relations:(3.3)Fq(k)(ω)=(-1)kFq(ω(k)),∀ω∈Ωn,ω⟶Fq(k)(ω)=(-1)k∫Rq(t)ω(k)(t)dt∈R,
where(3.4)ω(k):R⟶R,t⟶ω(k)(t)=dkω(t)dtk
is the k-order time derivative of the test function.Whenq∈Ck(ℜ), then(3.5)Fq(k)(ω)=∫Rq(k)(t)ω(t)dt=(-1)k∫Rq(t)ω(k)(t)dt
that means the k-order derivative of a distribution generated by a function q∈Ck(ℜ) equals the distribution generated by the k-order time derivative of the function q. 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 supportT=(tatb),ta<tb, all of these accomplishing the condition:(3.6)ω(t)=0;∀t∈(-∞,ta]∪[tb,+∞).The nonzero restriction is one of the following three types.(1)
Exponential:(3.7)ω(t)=exp⁡(|tatb|(t-ta)(t-tb)).(2)
Sinusoidal:(3.8)ω(t)=sinp((t-tb)πtb-ta).(3)
Polynomial:(3.9)ω(t)=(t-ta)p(t-tb)p,
where p≥2 is an integer.Figures1 and 2 present the exponential type test function and its first-order derivative for T=[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 rateD can be externally modified, it will be considered the first component of the input vector and the other component of u is the concentration Sin:u=[u1u2]=[DSin].From the third and fourth equations of the fermentative partial model (2.17) one gets the following relations for the reaction rates:(3.10)μ1ξ1=1k5(-ξ̇4-Dξ4+OTR),μ2ξ1=1k3(ξ̇3+Dξ3).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 ofk3 and k5).
Replacing (3.10) in the first equation of the system (2.17) one get
(3.11)(ξ̇1+Dξ1)=θ3(ξ̇3+Dξ3)+θ5(-ξ̇4-Dξ4+OTR),
where θ3=1/k3 and θ5=1/k5.
One denotes(3.12)v(t)=ξ̇1(t)+Dξ1(t),w1(t)=ξ̇3(t)+Dξ3(t),w2(t)=-ξ̇4(t)-Dξ4(t)+OTR.
Multiplying these functions with a test functionωi(t) and integrating over ℜ one gets the following distributions:
(3.13)Fvi(ωi)=∫Rv(t)⋅ωi(t)dt=∫Rξ̇1(t)⋅ωi(t)dt+∫RDξ1(t)⋅ωi(t)dt=∫Rξ1(t)⋅ωi(1)(t)dt+∫RDξ1(t)⋅ωi(t)dt,(3.14)Fw1i(φi)=∫Rw1(t)⋅ωi(t)dt=∫Rξ̇3(t)⋅ωi(t)dt+∫RDξ3(t)⋅ωi(t)dt=-∫Rξ3(t)⋅ωi(1)(t)dt+∫RDξ3(t)⋅ωi(t)dt,Fw2i(ωi)=∫Rw2(t)⋅ωi(t)dt=∫R-ξ̇4(t)⋅ωi(t)dt+∫R(-Dξ4(t)+OTR)⋅ωi(t)dt=∫Rξ4(t)⋅ωi(1)(t)dt+∫R(-Dξ4(t)+OTR)⋅ωi(t)dt.
As one must identify two parameters it is enough to choose a finite numberN≥2 of test functions ωi,i=1:N and to build an algebraic system of equations,
(3.15)Fw⋅θ=Fv,
where, θ=[θ3θ5]T is the vector of unknown parameters, Fw is a N×2 matrix of real numbers:
(3.16)Fw=[Fw11(ω1),…,Fw1i(ωi),…,Fw1N(ωN)Fw21(ω1),…,Fw2i(ωi),…,Fw2N(ωN)]T
and Fv denotes an N-column real vector build from distributions calculated as in (3.13):
(3.17)Fv=[Fv1(ω1),…,Fvi(ωi),…,FvN(ωN)]T.
If rank(Fw)=2 then a unique solution is obtained:
(3.18)θ̂=(FwT⋅Fw)-1FwT⋅Fv.
One obtains k̂3=1/θ̂3 and k̂5=1/θ̂5.Step 2 (identification ofk1 and k2).
Considering knownk3=k̂3 and k5=k̂5 from Step1 and substituting (3.10) in the second equation of (2.17) one gets
(3.19)(ξ̇2+Dξ2-DSin)=θ1(1k̂5(-ξ̇4-Dξ4+OTR))+θ2(1k̂3(ξ̇3+Dξ3)),
where θ1=-k1 and θ2=-k2.The procedure of identification fromStep1 is repeated and one obtains k̂1=-θ̂1 and k̂2=-θ̂2.Step 3 (identification ofk7 and k8).
Considering knownk3=k̂3 and k5=k̂5 from Step1 and substituting (3.10) in the last equation of (2.17) one gets
(3.20)(ξ̇5+Dξ5+CTR)=θ7(1k̂3(ξ̇3+Dξ3))+θ8(1k̂5(-ξ̇4-Dξ4+OTR)),
where θ7=k7 and θ8=k8.The procedure of identification fromStep1 is repeated and one obtains k̂7=θ̂7 and k̂8=θ̂8.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 k4, k6, and k9.


## 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Ωn the fundamental space from the distribution theory [24] of the real functions ω:ℜ→ℜ,t→ω(t) with compact support T, having continuous derivatives at least up to the order n. Let q:ℜ→ℜ,t→q(t) be a function which admits a Riemann integral on any compact interval T from ℜ. Using this function, a unique distribution (or generalized function)(3.1)Fq:Ωn⟶R,ω⟶Fq(ω)∈R
can be built by the relation:(3.2)Fq(ω)=∫Rq(t)ω(t)dt,∀ω∈Ωn.In distribution theory, the notion ofk-order derivative is introduced. If Fq∈Ωn, then its k-order derivative is a new distribution Fq(k)∈Ωn uniquely defined by the relations:(3.3)Fq(k)(ω)=(-1)kFq(ω(k)),∀ω∈Ωn,ω⟶Fq(k)(ω)=(-1)k∫Rq(t)ω(k)(t)dt∈R,
where(3.4)ω(k):R⟶R,t⟶ω(k)(t)=dkω(t)dtk
is the k-order time derivative of the test function.Whenq∈Ck(ℜ), then(3.5)Fq(k)(ω)=∫Rq(k)(t)ω(t)dt=(-1)k∫Rq(t)ω(k)(t)dt
that means the k-order derivative of a distribution generated by a function q∈Ck(ℜ) equals the distribution generated by the k-order time derivative of the function q. 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 supportT=(tatb),ta<tb, all of these accomplishing the condition:(3.6)ω(t)=0;∀t∈(-∞,ta]∪[tb,+∞).The nonzero restriction is one of the following three types.(1)
Exponential:(3.7)ω(t)=exp⁡(|tatb|(t-ta)(t-tb)).(2)
Sinusoidal:(3.8)ω(t)=sinp((t-tb)πtb-ta).(3)
Polynomial:(3.9)ω(t)=(t-ta)p(t-tb)p,
where p≥2 is an integer.Figures1 and 2 present the exponential type test function and its first-order derivative for T=[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 rateD can be externally modified, it will be considered the first component of the input vector and the other component of u is the concentration Sin:u=[u1u2]=[DSin].From the third and fourth equations of the fermentative partial model (2.17) one gets the following relations for the reaction rates:(3.10)μ1ξ1=1k5(-ξ̇4-Dξ4+OTR),μ2ξ1=1k3(ξ̇3+Dξ3).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 ofk3 and k5).
Replacing (3.10) in the first equation of the system (2.17) one get
(3.11)(ξ̇1+Dξ1)=θ3(ξ̇3+Dξ3)+θ5(-ξ̇4-Dξ4+OTR),
where θ3=1/k3 and θ5=1/k5.
One denotes(3.12)v(t)=ξ̇1(t)+Dξ1(t),w1(t)=ξ̇3(t)+Dξ3(t),w2(t)=-ξ̇4(t)-Dξ4(t)+OTR.
Multiplying these functions with a test functionωi(t) and integrating over ℜ one gets the following distributions:
(3.13)Fvi(ωi)=∫Rv(t)⋅ωi(t)dt=∫Rξ̇1(t)⋅ωi(t)dt+∫RDξ1(t)⋅ωi(t)dt=∫Rξ1(t)⋅ωi(1)(t)dt+∫RDξ1(t)⋅ωi(t)dt,(3.14)Fw1i(φi)=∫Rw1(t)⋅ωi(t)dt=∫Rξ̇3(t)⋅ωi(t)dt+∫RDξ3(t)⋅ωi(t)dt=-∫Rξ3(t)⋅ωi(1)(t)dt+∫RDξ3(t)⋅ωi(t)dt,Fw2i(ωi)=∫Rw2(t)⋅ωi(t)dt=∫R-ξ̇4(t)⋅ωi(t)dt+∫R(-Dξ4(t)+OTR)⋅ωi(t)dt=∫Rξ4(t)⋅ωi(1)(t)dt+∫R(-Dξ4(t)+OTR)⋅ωi(t)dt.
As one must identify two parameters it is enough to choose a finite numberN≥2 of test functions ωi,i=1:N and to build an algebraic system of equations,
(3.15)Fw⋅θ=Fv,
where, θ=[θ3θ5]T is the vector of unknown parameters, Fw is a N×2 matrix of real numbers:
(3.16)Fw=[Fw11(ω1),…,Fw1i(ωi),…,Fw1N(ωN)Fw21(ω1),…,Fw2i(ωi),…,Fw2N(ωN)]T
and Fv denotes an N-column real vector build from distributions calculated as in (3.13):
(3.17)Fv=[Fv1(ω1),…,Fvi(ωi),…,FvN(ωN)]T.
If rank(Fw)=2 then a unique solution is obtained:
(3.18)θ̂=(FwT⋅Fw)-1FwT⋅Fv.
One obtains k̂3=1/θ̂3 and k̂5=1/θ̂5.Step 2 (identification ofk1 and k2).
Considering knownk3=k̂3 and k5=k̂5 from Step1 and substituting (3.10) in the second equation of (2.17) one gets
(3.19)(ξ̇2+Dξ2-DSin)=θ1(1k̂5(-ξ̇4-Dξ4+OTR))+θ2(1k̂3(ξ̇3+Dξ3)),
where θ1=-k1 and θ2=-k2.The procedure of identification fromStep1 is repeated and one obtains k̂1=-θ̂1 and k̂2=-θ̂2.Step 3 (identification ofk7 and k8).
Considering knownk3=k̂3 and k5=k̂5 from Step1 and substituting (3.10) in the last equation of (2.17) one gets
(3.20)(ξ̇5+Dξ5+CTR)=θ7(1k̂3(ξ̇3+Dξ3))+θ8(1k̂5(-ξ̇4-Dξ4+OTR)),
where θ7=k7 and θ8=k8.The procedure of identification fromStep1 is repeated and one obtains k̂7=θ̂7 and k̂8=θ̂8.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 k4, k6, and k9.


## 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(FwT·Fw)-1 or, equivalently, Fw 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 integralIf=∫abf(x)dx where f(x) is a given function and [a,b] 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 If is a weighted sum of function values: If=∑i=0N-1ci·f(xi), where coefficients {ci} are derived from interpolating polynomial fitted for points {xi,fi}. A Newton-Cotes formula of any degree n can be constructed. One of the simplest integration methods is the trapezoidal rule. The trapezoidal rule is based on linear interpolation of f(x) at x1=a and x2=b, that is,  f(x) is approximated by(4.1)p(x)=f(a)+(x-a)f(b)-f(a)b-a.
The integral of p(x) equals the of trapezoid with base (b-a) times the average height(4.2)12(f(a)+f(b))
hence(4.3)∫abf(x)dx=b-a2(f(a)+f(b)).
To increase the accuracy one subdivides the interval [a,b] and assumes that fi=f(xi) is known on a grid of equidistant points x0=a,xi=x0+ih,xn=b, where h=(b-a)/n is the step length. The trapezoidal approximation for the ith subinterval is(4.4)∫xixi+1f(x)dx=h2(fi+fi+1)+Ri,
where(4.5)Ri=-h312f′′(ςi),ςi∈[xi,xi+1].
Summing the contributions for each subinterval [xi,xi+1],i=0:n, gives(4.6)∫abf(x)dx=h2(f0+fn)+h∑i=2i=n-1fi+RT.
The global truncation error is(4.7)RT=-h312∑i=0i=n-1f′′(ςi)=-h212(b-a)f′′(ς),ς∈[a,b].
This shows that by choosing h small enough we can make the truncation error arbitrarily small. In other words, we have asymptotic convergence when h→0.In the composite Simpson’s rule one divides the interval[a,b] into an even number n=2 m steps of length h and uses the formula(4.8)∫abf(x)dx=h3(f0+4S1+2S2+fn)+RT,
where(4.9)S1=f1+f3+⋯+fn-1,S2=f2+f4+⋯+fn-2
are sums over odd and even indices, respectively. The remainder is(4.10)RT=h4180(b-a)fIV(ς),ς∈[a,b].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(4.11)H(z)=h2⋅1+z-11-z-1,H(z)=h3⋅1+4z-1+z-21-z-1,respectively.These transfer functions have the amplitude-frequency Bode characteristics from Figure3.Amplitude-frequency bode characteristics: (a) trapezoidal rule; (b) Simpson’s rule.
(a)(b)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 Figure3 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 φi(ξ) that are (nonlinear) functions of the state variables. These functions enter the model in the form Kφi(ξ) (where K is a constant matrix), which is a set of linear combinations of the same nonlinear functions φi(ξ). 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)(KaKb)=TK,
with Ka∈Rp×m of full row rank matrix (i.e., p=rank(K)),Kb∈R(n-p)×m and T a permutation matrix. The induced partitions of the vectors ξ and u are(4.13)(ξaξb)=Tξ(uaub)=Tu.Model (2.1) is then partitioned into two submodels:(4.14)dξadt=K⋅φ(ξa,t)-Dξa+ua,dξbdt=K⋅φ(ξb,t)-Dξb+ub.
Then with the state transformation(4.15)ξa=ξa,z=Cξa+ξb,
one transform the initial model into(4.16)ξ̇a=Kaφ(ξa,z-Cξa)-Dξa+ua,ż=-Dz+Cua+ub,
where the (n-p)×p matrix C is the unique solution of(4.17)CKa+Kb=0
that is,(4.18)C=-KbKa*,
where Ka* is a generalized or pseudoinverse of Ka. This subsystem (4.16) can be augmented with an equation derived from (2.7) as follows:(4.19)ż=-Dz+Cua+ub,ξb=z-Cξa.It can be considered as a linear time-varying (ifd varies in the course of time) model with state z, input (ξa,ua,ub), and output ξb. It is nonlinearly parametrized by the yield coefficients but linearly reparametrized by the nonzero entries of C. When data of the signals ξa, ua, ub, and ξb 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 ki recovered afterwards from (4.17). The identification makes sense only if the yield coefficients ki are structurally uniquely identifiable with the auxiliary model (4.19).


## 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(FwT·Fw)-1 or, equivalently, Fw 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 integralIf=∫abf(x)dx where f(x) is a given function and [a,b] 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 If is a weighted sum of function values: If=∑i=0N-1ci·f(xi), where coefficients {ci} are derived from interpolating polynomial fitted for points {xi,fi}. A Newton-Cotes formula of any degree n can be constructed. One of the simplest integration methods is the trapezoidal rule. The trapezoidal rule is based on linear interpolation of f(x) at x1=a and x2=b, that is,  f(x) is approximated by(4.1)p(x)=f(a)+(x-a)f(b)-f(a)b-a.
The integral of p(x) equals the of trapezoid with base (b-a) times the average height(4.2)12(f(a)+f(b))
hence(4.3)∫abf(x)dx=b-a2(f(a)+f(b)).
To increase the accuracy one subdivides the interval [a,b] and assumes that fi=f(xi) is known on a grid of equidistant points x0=a,xi=x0+ih,xn=b, where h=(b-a)/n is the step length. The trapezoidal approximation for the ith subinterval is(4.4)∫xixi+1f(x)dx=h2(fi+fi+1)+Ri,
where(4.5)Ri=-h312f′′(ςi),ςi∈[xi,xi+1].
Summing the contributions for each subinterval [xi,xi+1],i=0:n, gives(4.6)∫abf(x)dx=h2(f0+fn)+h∑i=2i=n-1fi+RT.
The global truncation error is(4.7)RT=-h312∑i=0i=n-1f′′(ςi)=-h212(b-a)f′′(ς),ς∈[a,b].
This shows that by choosing h small enough we can make the truncation error arbitrarily small. In other words, we have asymptotic convergence when h→0.In the composite Simpson’s rule one divides the interval[a,b] into an even number n=2 m steps of length h and uses the formula(4.8)∫abf(x)dx=h3(f0+4S1+2S2+fn)+RT,
where(4.9)S1=f1+f3+⋯+fn-1,S2=f2+f4+⋯+fn-2
are sums over odd and even indices, respectively. The remainder is(4.10)RT=h4180(b-a)fIV(ς),ς∈[a,b].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(4.11)H(z)=h2⋅1+z-11-z-1,H(z)=h3⋅1+4z-1+z-21-z-1,respectively.These transfer functions have the amplitude-frequency Bode characteristics from Figure3.Amplitude-frequency bode characteristics: (a) trapezoidal rule; (b) Simpson’s rule.
(a)(b)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 Figure3 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 φi(ξ) that are (nonlinear) functions of the state variables. These functions enter the model in the form Kφi(ξ) (where K is a constant matrix), which is a set of linear combinations of the same nonlinear functions φi(ξ). 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)(KaKb)=TK,
with Ka∈Rp×m of full row rank matrix (i.e., p=rank(K)),Kb∈R(n-p)×m and T a permutation matrix. The induced partitions of the vectors ξ and u are(4.13)(ξaξb)=Tξ(uaub)=Tu.Model (2.1) is then partitioned into two submodels:(4.14)dξadt=K⋅φ(ξa,t)-Dξa+ua,dξbdt=K⋅φ(ξb,t)-Dξb+ub.
Then with the state transformation(4.15)ξa=ξa,z=Cξa+ξb,
one transform the initial model into(4.16)ξ̇a=Kaφ(ξa,z-Cξa)-Dξa+ua,ż=-Dz+Cua+ub,
where the (n-p)×p matrix C is the unique solution of(4.17)CKa+Kb=0
that is,(4.18)C=-KbKa*,
where Ka* is a generalized or pseudoinverse of Ka. This subsystem (4.16) can be augmented with an equation derived from (2.7) as follows:(4.19)ż=-Dz+Cua+ub,ξb=z-Cξa.It can be considered as a linear time-varying (ifd varies in the course of time) model with state z, input (ξa,ua,ub), and output ξb. It is nonlinearly parametrized by the yield coefficients but linearly reparametrized by the nonzero entries of C. When data of the signals ξa, ua, ub, and ξb 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 ki recovered afterwards from (4.17). The identification makes sense only if the yield coefficients ki 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 as(5.1)NMSE=1N∑j=1NNMSE(θ̂j),withNMSE(θ̂j)=(θ̂j-θj*θj*)2,
where N is the number of estimated parameters, θ̂j is the jth 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 periodTs, 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:Ts∈{1min,10min,20min},(b)
input signals:u∈{0.1+0.05*sin(t),0.1+0.005*t,0.15*exp(9*10-3t)},(c)
initial conditions:
IC1: [1 0.2 0.15 0.0066 0.008]IC2: [0.1 0.3 0 0.002 0.01 2]IC3: [0.5 0.5 3 0.01 0.001 1],(d)
test functions:
Type1: exponentialType2: sinusoidalType3: 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 Figures4–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.Test functions (a) and their derivatives (b) used for identification.
(a)(b)Figure 6
Noisy measurements of state variableE: SNR = 20 dB.The simulation results presented in Tables1–5 leads to the following remarks.(1)
Sensitivity to the sampling period: results given in Table1 show that the performance of the algorithm deteriorates with the increasing sampling period.(2)
Sensitivity to the inputs type: results given in Table2 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 Tables3 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: Tables5 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.
Tsk1k2k3k4k5k6k7k8k9NMSE1 min0.49010.05000.10000.71991.20010.64010.81010.11001.11011.6046e-00810 min0.50220.04970.10000.70791.20380.65080.83570.10941.12872.8414e-00420 min1.21370.03810.09940.01011.49801.27102.32050.08402.20440.9722“Real” value0.490.050.10.721.20.640.810.111.11Table 2
Influence of the input type.
Input typek1k2k3k4k5k6k7k8k9NMSE0.1 + 0.05 * sin(t) (noise free)0.49010.05000.10000.71991.20010.64010.81010.11001.11011.6046e-0080.1 + 0.05 * sin(t)  (SNR = 20 dB)0.45310.05170.10000.69261.24590.66440.79690.11121.15230.00150.1 + 0.005 *t (noise free)0.49020.05000.10000.71981.20000.64020.80960.11001.11036.9366e-0080.1 + 0.005 *t (SNR = 20 dB)0.38420.05310.10020.63191.07300.71830.67540.11281.24590.01500.15 * exp(9 * 10−3t)  (noise free)0.49020.05000.10000.71981.20000.64020.80970.11001.11036.3148e-0080.15 * exp(9 * 10−3t)  (SNR = 20 dB)0.48820.04880.10040.66601.07380.68800.69000.11111.19320.0056“Real” value0.490.050.10.721.20.640.810.111.11Table 3
Influence of the initial conditions.
Initial conditionsk1k2k3k4k5k6k7k8k9NMSEIC10.49010.05000.10000.71991.20010.64010.81010.11001.11011.6046e-008IC20.49010.05000.10000.71971.19990.64010.80980.11001.11043.1454e-008IC30.49080.05000.10000.71781.20020.64190.81580.10991.11349.1629e-006“Real” value0.490.050.10.721.20.640.810.111.11Table 4
Influence of the test functions type.
Test functionk1k2k3k4k5k6k7k8k9NMSEtype 10.49020.04990.10000.71981.19990.64010.81020.11001.11023.2506e-008type 20.49010.05000.10000.71991.20000.64010.81020.11001.11011.5545e-008type 30.49080.05000.10000.71951.19930.64050.80960.11001.11085.4343e-007“Real” value0.490.050.10.721.20.640.810.111.11Table 5
Influence of the noise (distribution-based method).
SNRk1k2k3k4k5k6k7k8k9NMSEnoise free0.49010.05000.10000.71991.20010.64010.81010.11001.11011.6046e-00850 dB0.48240.05020.10000.71471.19650.64470.80370.11021.11825.4504e-00540 dB0.43930.05090.10020.68031.13580.67530.75680.11021.17120.003030 dB0.47070.04940.10070.68611.09530.67010.85400.10651.16220.002220 dB0.45310.05170.10000.69261.24590.66440.79690.11121.15230.001510 dB0.55670.05290.09930.55871.77580.78340.98870.11191.35870.0502“Real” value0.490.050.10.721.20.640.810.111.11Table 6
Influence of the noise (“classical” method).
SNRk1k2k3k4k5k6k7k8k9NMSEnoise free0.49030.05000.10000.71951.20010.64040.81090.11001.11053.1151e-0750 dB0.48760.04890.10030.70371.17480.65330.81990.10721.12543.2795e-0440 dB0.45010.05740.09860.64401.24320.70240.84740.11951.18210.007230 dB0.39230.06590.09660.05273.22901.18750.83930.13381.74310.551720 dB0.01330.07110.09830.06833.27321.28680.53920.13031.85780.769910 dB0.31150.07760.09770.98076.11952.03550.76250.12542.72353.5841“Real” value0.490.050.10.721.20.640.810.111.11


## 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.



---

*Source: 789156-2012-03-01.xml*