- Email: [email protected]
Procedure for asymptotic state and parameter estimation of nonlinear distributed parameter bioreactors
Applied Mathematical Modelling 31 (2007) 1293–1307 www.elsevier.com/locate/apm
Procedure for asymptotic state and parameter estimation of nonlinear distributed parameter bioreactors T. Damak Unite´ de Commande Automatique, Ecole Nationale d’Inge´nieurs de Sfax, Sfax BP W 3038, Tunisia Received 1 November 2004; received in revised form 1 November 2005; accepted 20 February 2006 Available online 5 July 2006
Abstract Due to measurement problem in biological variables, such as the biomass concentration and its specific growth rate, the on-line measurement of these variables is not ensured and this presents a great difficulty, especially when we consider the control problem. This paper presents and analyses an approach for estimating these biological variables through an example of a nonlinear distributed parameter bioreactor. An orthogonal collocation method is applied in the distributed system to obtain a system of ordinary differential equations. A simulation study shows the feasibility and the robustness of the estimator used for this nonlinear process. 2006 Elsevier Inc. All rights reserved. Keywords: State and parameter estimation; Nonlinear systems; Distributed parameter bioreactor; Orthogonal collocation; Asymptotic stability
1. Introduction When studying a control problem of biotechnological processes, many difficulties occur. In fact, models of these processes are in general nonlinear in character and not very accurate. Furthermore, unavailability of measurements specially the biomass concentration and its specific growth rate causes problems. To overcome these difficulties, several estimation algorithms of different type have been proposed, e.g. extended adaptive Luenberger observer [1], gradient based filter [2], modulated gain estimator [3], extended Kalman filter [4], state and parameter estimation based on sliding mode [5,6] and nonlinear observers [7,8]. One solution is proposed in this paper. The system to be studied is a fixed bed bioreactor. It is described by a set of nonlinear distributed parameter equations introduced in [5,9]. The solution used consists in: – applying a functional approximation method in order to get an ordinary differential system, in an optimal or quasi-optimal way: the orthogonal collocation method [5,10] has been used to approximate a mathematical model of hyperbolic type; E-mail address: [email protected] 0307-904X/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2006.02.014
1294
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
– deriving an estimator which gives simultaneously the estimates of variable parameters and concentrations of the bioreactor reactants; this estimator has been built from local asymptotic stability properties. The paper is organized as follows: Section 2 describes the distributed parameter process with its mathematical model. The reduction of this model to ordinary differential equations is presented in Section 3. Section 4 describes the estimation problem. The estimation procedure is developed and discussed in Section 5. The performance of the estimator is illustrated by simulation experiments in Section 6. 2. Distributed parameter system model Consider a fixed bed bioreactor with two reactions described as follows [6,11,12]: – an autocatalytic growth reaction with one limiting substrate and one biomass whose concentrations are respectively s and x: s → x;
– a death reaction of the micro-organisms: x ! xd ; where xd is the non-active biomass concentration. The micro-organisms are fixed on some support and remain within the reactor tube (Fig. 1). The substrate liquid flows through the reactor; the micro-organisms grow by consuming the substrate whose flow rate can be the control variable [13]. Assuming that the non-active micro-organisms leave the reactor and neglecting diffusion phenomenon, the above reactions can be described by the following distributed parameter system of hyperbolic type: oxðz; tÞ ¼ ½lðx; sÞ k d xðz; tÞ; ot osðz; tÞ osðz; tÞ ¼ u k 1 lðx; sÞxðz; tÞ ot oz
ð1Þ
over the space domain 0 < z 6 1 with the boundary condition s(0, t) = sin(t). In system (1), the parameters are defined as follows: l kd k1 u
specific growth rate death coefficient yield coefficient F/A (F is the feeding flow rate and A is the cross section of the bioreactor)
Though the above model is of academic type, its structure is truly representative of biotechnological processes, e.g. in the fields of biological wastewater treatment [14] or of industrial fermentations. packing
biomass
F(t)
F(t)
sin(t)
sout(t) z=0
z=1
Fig. 1. Schematic view of a fixed bed reactor.
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
1295
3. Lumped parameter system model By using the orthogonal collocation method [5,10,15] system (1) is transformed into a set of nonlinear ordinary differential equations. This functional approximation method is particularly suited for digital computation when compared with other approximation methods. It consists in expanding the model variables in terms of orthogonal polynomials Lj(z) (e.g. orthogonal Lagrange polynomials) [16–18]: xðz; tÞ ¼
N þ1 X
Lj ðzÞxj ðtÞ;
j¼0
sðz; tÞ ¼
N þ1 X
ð2Þ Lj ðzÞsj ðtÞ;
j¼0
xj(t) and sj(t) are concentrations evaluated at point z = zj. Polynomials Lj(z) have the following property: Lj ðzi Þ ¼ 1
if j ¼ i;
Lj ðzi Þ ¼ 0
if j 6¼ i:
ð3Þ
The collocation points (zj, j = 0, 1, . . . , N + 1) are chosen, as follows: 0 ¼ z0 < z1 < zN < zN þ1 ¼ 1:
ð4Þ ðp;qÞ PN
The N internal collocation points are obtained by calculating the zeros of orthogonal Jacobi polynomials [10] defined by Z 1 ðp;qÞ p zq ð1 zÞ zk P N ðzÞdz ¼ 0 k ¼ 0; 1; . . . ; N 1; ð5Þ 0
z0 is the inlet of the reactor; zN+1 is the outlet of the reactor. The determination of the number and the location of collocation points has been developed in [5,19]. In system (1), the space derivative of s(z, t) is then expressed as follows: N þ1 osðz; tÞ X dLj ðzÞ ¼ sj ðtÞ: oz dz j¼0
ð6Þ
Setting lij ,
dLj ðz ¼ zi Þ dz
ð7Þ
and applying the collocation method, the (N + 1) dimensional differential system is derived: dxi ¼ ðli k d Þxi ; dt N þ1 X dsi ¼ u lij sj uli0 s0 k 1 li xi dt j¼1
ð8Þ
ði ¼ 1; 2; . . . ; N þ 1Þ; where li , l(xi, si); s0 = sin is the boundary condition. 4. Estimation problem In what follows, we will assume that parameters kd and k1 are known and constants. On the contrary, the specific growth rate l is known to be a complex function of concentrations x and s, and of some biological parameters; a lot of mathematical models have been proposed [9] considering that l ¼ lðh; s; xÞ;
1296
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
where h is the vector of kinetic parameters; we will assume that the derivative of l with respect to h is bounded. Generally speaking, l can be considered as a time-varying parameter l(t) to be estimated in real time. The substrate concentration s(z, t), which is always less than the feeding concentration sin, and the biomass concentration are positives. The velocity u(t) is bounded (0 < u(t) 6 umax). Moreover, we will assume that we can measure only substrate concentrations and that sensors are located at collocation points. So, the output equation is written as y ¼ ½s1 ; s2 ; . . . ; sN þ1 :
ð9Þ
In general, sensors are located at different points of the reactor but not necessarily at the collocation points. In that case, the output variable y can be estimated with a good precision by using the interpolation function (2), which is needed when applying the collocation procedure. One method of estimating simultaneously the specific growth rate l and the biomass concentration x(z, t), which is the main objective of this study, is presented in the next section. 5. Estimation procedure The state estimation model is given by the following equations: d^xi ¼ ð^ li k d Þ^xi aTi es ; dt N þ1 X d^si ^i^xi uli0 sin bTi es ¼ u lij^sj k 1 l dt j¼1 ði ¼ 1; 2; . . . ; N þ 1Þ;
ð10Þ
^i are respectively the estimates of xi, si and li; es is the error on substrate concentration at where ^xi ; ^si and l each collocation point zi: eTs ¼ ½ es1 es2 . . . esN þ1 with: esi ¼ ^si si ; ai and bi are the corresponding adaptive gains. The specific growth rate in the estimation model (10) has the same structure as in the original system model: ^ ¼ lð^ l h; ^x; sÞ: ð11Þ Then, the estimation problem consists of determining gains ai and bi (i = 1, 2, . . . , N + 1) in such a way that exi ¼ ^xi xi ! 0; ^i li ! 0; eli ¼ l esi ¼ ^si si ! 0 when t ! 1:
ð12Þ
Without specifying the structure of the l function, an error system can be derived for estimating biomass concentrations and specific growth rates at collocation points, by using a first order Taylor series expansion of l [20]: ^i li ¼ HTi ehi þ K i ; exi ; eli ¼ l
ð13Þ
where eh i ¼ ^ hi h 2 Rm ; ol 2 Rm ; Hi ¼ oh h;xi ;si ol Ki ¼ 2 R; ox h;xi ;si m being the number of parameters in h.
ð14Þ ð15Þ ð16Þ
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
1297
Subtracting (8) from (10) and using (13), we get the error system: dexi ¼ ðli þ K i^xi k d Þexi aTi es þ ^xi HTi ehi ; dt N þ1 X desi ¼ k 1 ðli þ K i^xi Þexi u lij esj bTi es k 1^xi HTi ehi dt j¼1
ð17Þ
ði ¼ 1; 2; . . . ; N þ 1Þ: 5.1. Concentration estimation Defining: ai ,li þ K i^xi k d ; bi , k 1 ðli þ K i^xi Þ ¼ k 1 ðai þ k d Þ
ð18Þ
ði ¼ 1; 2; . . . ; N þ 1Þ; system (17) is written as follows: N þ1 X dexi ¼ ai exi aij esj þ ^xi HTi ehi ; dt j¼1 N þ1 X desi ¼ bi exi ðulij þ bij Þesj k 1^xi HTi ehi dt j¼1
ð19Þ
ði ¼ 1; 2; . . . ; N þ 1Þ and finally, with a specific distribution of the error vector eT ¼ ½ ex1 es1 ex2 es2 . . . exN þ1 esN þ1 ;
ð20Þ
we can write (19) in a condensed manner: de ¼ Ae þ Bxeh ; dt es ¼ Ce:
ð21Þ
• The matrix A is structured into (N + 1)2 matrices Aij of the form: ai aii Aii ¼ if i ¼ j ¼ 1; 2; . . . ; N þ 1; bi ðulii þ bii Þ " # 0 aij if i 6¼ jði; j ¼ 1; 2; . . . ; N þ 1Þ: Aij ¼ 0 ðulij þ bij Þ • B is a matrix of dimension (2N + 2) Æ (N + 1): 3 2 1 0 0 6 k 1 0 0 7 7 6 7 6 1 0 7 6 0 7 6 6 k 1 0 7 B¼6 0 7: 6 . . 7 7 6 . . . 7 6 . 7 6 4 0 0 1 5 0
0
k 1
1298
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
• x is a matrix of dimension (N + 1) Æ [m(N + 1)]: 3 2 ^x1 HT1 O O 7 6 ^x2 HT2 O 7 6 O 7 6 7 6 . . 7: 6 . . w¼6 . . 7 7 6 . .. 7 6 . . 5 4 . O x is bounded. • C is a matrix of 2 0 1 60 0 6 6. . 6 C ¼ 6 .. .. 6 6 .. .. 4. . 0 •
eTh
¼
½ eTh1
^xN þ1 HTN þ1
O
0 eTh2
dimension (N + 1) Æ (2N + 2): 3 0 0 0 0 0 1 0 07 7 .. .. 7 .. .. 7 . . 7: . . 7 .. .. 7 .. .. . .5 . . 0
0
...
eThNþ1
0 1
.
The estimation error gains ai and bi must ensure the local asymptotic stability of the dynamical error matrix A. Proposition. With an appropriate choice of aij and bij (i 5 j) the matrix A will have a block-diagonal structure; we propose: aij ¼ 0; ulij þ bij ¼ 0
ð22Þ
ði; j ¼ 1; 2; . . . ; N þ 1; i 6¼ jÞ: Choosing aij = 0 (i 5 j) is compatible with the fact that the biomass concentration at point zi depends only on the substrate concentration at this point. Choosing ulij + bij = 0 (i 5 j) corresponds to uncorrelation between the different substrate concentration errors esi ði ¼ 1; 2; . . . ; N þ 1Þ. So, the block-diagonal matrix A is structured as follows: A ¼ diag½A1 A2 Ai . . . AN þ1 with
ai Ai ¼ bi
aii ulii bii
i ¼ 1; 2; . . . ; N þ 1:
ð23Þ
If each sub-matrix Ai is asymptotical stable, the whole error system (21) will also be asymptotical stable. The characteristic polynomial of Ai is given by P i ðkÞ ¼ k2 kðai ulii bii Þ þ bi aii ai ðulii þ bii Þ:
ð24Þ
The gains aii and bii must be bounded; assuming that their variations are also bounded, then considering the i assumptions given in Section 4, the sub-matrix Ai and its time derivate dA are bounded. The asymptotic stadt bility is obtained if the following conditions hold [21]: ri ¼ ai ulii bii < 0 ðsum of eigenvaluesÞ; pi ¼ bi aii ai ðulii þ bii Þ > 0 ðproduct of eigenvaluesÞ:
ð25Þ
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
1299
Choosing two eigenvalues k1i and k2i with negative real part, ri and pi are then fixed. So, from (25) we get: bi aii ¼ pi þ ai ðai ri Þ: Since bi = k1(ai + kd) we obtain the gains aii and bii: aii ¼
pi ai ri þ a2i k 1 ðai þ k d Þ
bii ¼ ai ulii ri
i ¼ 1; 2; . . . ; N þ 1;
i ¼ 1; 2; . . . ; N þ 1:
ð26Þ ð27Þ
The expression (26) can also be written as follows: aii ¼
1 ðai k1i Þðai k2i Þ : k1 ai þ k d
ð28Þ
Taking into account the proposition (22), the formulation of the state estimation model (10) is simplified: d^xi ¼ ð^ li k d Þ^xi aii ð^si si Þ; dt N þ1 X d^si ^i^xi uli0 sin þ ¼ ðulii þ bii Þ^si k 1 l bij sj dt j¼1
ð29Þ
ði ¼ 1; 2; . . . ; N þ 1Þ: Gains aii are given by (26) or (28); gains bij are determined by (22) and (27). Moreover, the error system (19) becomes: dexi ¼ ai exi aii esi þ ^xi HTi ehi ; dt desi ¼ bi exi ðai ri Þesi k 1^xi HTi ehi : dt
ð30Þ
5.2. Parameter estimation The resolution of the state estimation model (29) needs the knowledge of the estimate of the specific growth rate l depending on the parameter vector h at each point zi. Defining xi ¼ ^xi Hi , the system (30) can be written as follows: dexi ¼ ai exi aii esi þ xTi ehi ; dt desi ¼ bi exi ðai ri Þesi k 1 xTi ehi : dt
ð31Þ
The parameter estimation should be derived from the knowledge of the error on the substrate concentration estimate: d^ hi ¼ ci ðtÞesi dt
i ¼ 1; 2; . . . ; N þ 1;
ð32Þ
where the gain vector ci has to be determined. One solution consists in calculating ci in such a way that the augmented error vector ½ exi esi eli tends to ^ i li . zero, when time tends to infinity (asymptotic stability), with eli ¼ l In the following, we are considering two cases according as l is a function of one parameter or more. If we assume that h is constant or slowly variable, then: _ _ e_ hi ¼ ^ hi h_ ffi ^ hi :
1300
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
5.2.1. l is a function of one parameter h The whole error differential system is written as 2 3 2 32 3 e_ xi exi ai aii xi 6 7 6 76 7 4 e_ si 5 ¼ 4 bi ðai ri Þ k 1 xi 54 esi 5; 0 ci 0 ehi e_ hi
ð33Þ
whose characteristic polynomial is given by P i ðkÞ ¼ k3 ri k2 þ ðk 1 ci xi þ pi Þk þ k 1 k d ci xi :
ð34Þ
The asymptotic stability is ensured if the following conditions hold (ri being <0): ð35Þ ð36Þ
k 1 k d ci xi > 0; ri ðk 1 ci xi þ pi Þ k 1 k d ci xi > 0: In (35), xi can be positive or negative; if we choose ci ðtÞ ¼ Ci xi ðtÞ ¼ Ci^xi Hi
ð37Þ
with Ci positive, then (35) is always fulfilled. Choosing an adaptive gain Ci given by Ci ¼ x12 , yields to a constant coefficient characteristic polynomial: i
P i ðkÞ ¼ k3 ri k2 þ ðk 1 þ pi Þk þ k 1 k d :
ð38Þ
Then the condition (36) becomes: ri ðk 1 þ pi Þ k 1 k d > 0;
ð39Þ
which holds with an appropriate choice of ri and pi. One can notice that (39) is fulfilled by choosing ri < kd, which corresponds to a realistic condition, since the dynamics of the concentration estimator must be sufficiently rapid with respect to the system dynamics. 5.2.2. l is a function of several parameters The whole error differential system is now of the form: 32 3 2 3 2 exi e_ xi ai aii xTi 76 7 6 7 6 T 4 e_ si 5 ¼ 4 bi ðai ri Þ k 1 xi 54 esi 5; ehi e_ hi 0 ci 0
ð40Þ
whose characteristic polynomial is given by P i ðkÞ ¼ km1 ½k3 ri k2 þ ðk 1 xTi ci þ pi Þk þ k 1 k d xTi ci ;
ð41Þ
m being the dimension of h. The asymptotic stability of the whole system is no more ensured because there exist null eigenvalues (weak stability). Conditions (35) and (36) are now generalized: k 1 k d xTi ci > 0;
ri ðk 1 xTi ci
þ pi Þ
ð42Þ k 1 k d xTi ci
> 0:
ð43Þ
As before, the first condition (42) is fulfilled if: ci ðtÞ ¼ Ci xi ðtÞ;
ð44Þ
where Ci is a diagonal positive definite matrix. To obtain a characteristic polynomial with constant coefficients, we can use an adaptive gain Ci defined by 1 1 1 Ci ¼ diag 2 2 2 : ð45Þ xi1 xi2 xim
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
1301
The characteristic polynomial is then given by P i ðkÞ ¼ km1 ½k3 ri k2 þ ðmk 1 þ pi Þk þ mk 1 k d :
ð46Þ
The condition (43) becomes: ri ðmk 1 þ pi Þ mk 1 k d > 0;
ð47Þ
which also holds with an appropriate choice of ri and pi. It can be shown that the error on concentration estimates tends to zero; the error on the parameter ^hi does ^i tends to not diverge, but ^ hi tends in general towards wrong values, due to the null eigenvalues. Nevertheless l the true value of li, because, when exi and esi tend asymptotically to zero, we can deduce from (30) that ^xi HTi ehi ! 0 ^i li ¼ HTi ehi þ K i exi ! 0, when t ! 1. then, eli ¼ l In conclusion, the parameter estimator is given by the following differential system: d^ hi ¼ Ci^xi Hi esi ; dt
ð48Þ
where Ci is a diagonal positive definite matrix defined by (45). With this choice, the only tuning parameters of the resulting adaptive state and parameter estimator are ri and pi. The above results depend on the true values of li(xi, si, hi), Ki(xi, si, hi) and Hi(xi, si, hi). Since these true values are supposed to be unknown, they should be substituted by the corresponding estimated values b i ð^xi ; si ; ^ b i ð^xi ; si ; ^ ^i ð^xi ; si ; ^ l hi Þ; K hi Þ and H hi Þ. This substitution can be justified by using the L2-robustness theorem of Zeng et al. [22]. 6. Simulation results Let us consider for example, the specific growth rate given by the Contois model: lðh; x; sÞ ¼
lm s ; kcx þ s
ð49Þ
T
where h ¼ ½ lm k c is a constant parameter vector; lm denotes the maximum specific growth rate and kc is the saturation constant. In this case, we have: " #T ^ s s x i i i bi ¼ ^ mi H l ; 2 ^k c ^xi þ si ð^k ci ^xi þ si Þ i ð50Þ ^k c si i b i ¼ ^ K lm i : ð^k c ^xi þ si Þ2 i
6.1. Numerical data For simulation needs, parameters were set to the following numerical values [11]: A ¼ 0:02 m2 ;
k d ¼ 0:05 h1 ;
k 1 ¼ 0:4
lm ¼ 0:35 h1 ;
k c ¼ 0:4:
For the process simulation, and then generating data, we have used a model different from the lumped model given by the orthogonal collocation. This model is based on a high-dimensional finite difference method. The scheme used (central difference) is the following: oX X iþ1 X i1 : ð51Þ ¼ oz z¼zi 2Dz
1302
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
To have a good accuracy, an appropriate space step size Dz, equal to 0.002, is used corresponding to 500 space intervals. For the lumped model used in the estimator, four internal collocation points are defined. They were optimally located (p = 0; q = 4) [5] at ½ 0:312
0:579
0:813
0:963 :
The initial conditions for the process are chosen so as to correspond to a steady-state with F = 2 l/h, and sin(0) = s0(0) = 5 g COD/l. At the collocation points we get the following initial values: sð0Þ ¼ ½ 1:96
0:88
xð0Þ ¼ ½ 29:4 13:2
0:44
0:28
6:55 4:18
0:25 T in g COD=l; T
3:73 in g VSS=l:
The generating data, i.e., substrate measurements, are made with a sampling period of Te = 0.1 h. To obtain more realistic results, a zero mean-valued Gaussian noise with a standard deviation 0.005 is added into the measurement. The numerical solution method used to solve the set of ordinary differential equations is a fourth order Runge–Kutta method with a fixed time step equal to 0.01 h. 6.2. Simulation runs The practical adaptive state and parameter estimator given above is studied via simulations. The tuning parameters of the estimator are setting to the values ri = 5.05 and pi = 0.25. For testing the estimator convergence and robustness properties, simulation experiments, in which the input variable u and the parameters lm and kc are modified (step functions), are considered. Figs. 2–5 represent respectively the evolution, at the outlet of the reactor (point z5 = 1.0), of: – – – –
the the the the
substrate concentration s and its estimate, biomass concentration x and its estimate, specific growth rate l and its estimate, parameters lm and kc and their estimates,
Fig. 2. The evolution of substrate concentration at z = 1.0 (kc and lm are unknown).
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
1303
Fig. 3. The evolution of biomass concentration at z = 1.0 (kc and lm are unknown).
Fig. 4. The evolution of specific growth rate at z = 1.0 (kc and lm are unknown).
corresponding to the following changes: – – – – – – –
at at at at at at at
time time time time time time time
t = 25 h: flow rate step from 2 to 2.5 l/h, t = 50 h: flow rate step from 2.5 to 2 l/h, t = 75 h: flow rate step from 2 to 2.5 l/h, t = 100 h: flow rate step from 2.5 to 2 l/h; lm changes from 0.35 to 0.42 h1, t = 125 h: flow rate step from 2 to 2.5 l/h, t = 150 h: flow rate step from 2.5 to 2 l/h, t = 175 h: flow rate step from 2 to 2.5 l/h; kc changes from 0.4 to 0.48.
As can be seen from the figures, we can verify that the concentration estimates, i.e., the state variables of the system (s, x), tend in a satisfactory way to the actual values in spite of perturbation steps and noisy measurements.
1304
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
Fig. 5. The evolution of biological parameters at z = 1.0 (kc and lm are unknown).
Fig. 6. The evolution of substrate concentration at z = 0.312 (kc is known, lm unknown).
^ tends also to the actual value although the biological parameters l ^m and ^k c do The specific growth rate estimate l not converge to the true values as shown before in the case where l is a function of several parameters. This unidentifiability of the system parameters can be explained in another way [22] by considering the structure of the specific growth rate of Contois’s model lm s lðh; x; sÞ ¼ : ð52Þ kcx þ s For constant substrate and biomass concentrations there are many values lm and kc that give the same value l for the specific growth rate. Consequently, depending on the initial values of the parameters, estimates converge to some feasible values that give the true value to the specific growth rate.
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
1305
Fig. 7. The evolution of biomass concentration at z = 0.312 (kc is known, lm unknown).
Fig. 8. The evolution of specific growth rate at z = 0.312 (kc is known, lm unknown).
A second simulation experiment is done in the case where one of the two system parameter estimates is assumed to be known (e.g. kc). The results of the corresponding run are shown in Figs. 6–9 at point z1 = 0.312. The realized changes are identical to the ones given above except at time t = 175 h where the step change of kc is not considered. ^m ) converges to its correct value. This We can notice, in this case, that the parameter to be estimated (i.e. l can be seen in Fig. 9. It is possible to reconstruct the evolution of concentrations x(z, t) and s(z, t) at any point of the reactor by using expressions (2).
1306
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
Fig. 9. The evolution of biological parameters at z = 0.312 (kc is known, lm unknown).
7. Conclusion In this paper, it has been shown how it is possible to design an asymptotic estimator of state and time-varying parameters in the case of a nonlinear distributed parameter bioreactor. The objective was to have good estimates of the biomass concentration and its specific growth rate from measurements of the substrate concentration. The structure of the estimator is based upon an approximated model of the bioreactor behaviour. The excellent properties of the estimator, especially robustness to parameter variations, model mismatch and measurement errors, were tested through simulations. The results obtained show that the concentration estimates converge in a good way to the true values. The specific growth rate (time-varying parameter) is fairly well estimated, even if it is a function of several unknown parameters, which should be identified if needed. References [1] G. Bastin, D. Dochain, Adaptive estimation of microbial growth rate, in: H.A. Barker, P.C. Young (Eds.), Preprints 7th IFAC/ IFORS Symposium on Identification and System Parameter Estimator, York, United Kingdom, July 3–7, vol. 2, Pergamon Press, 1988, pp. 1161–1166. [2] G. Chamilothois, Y. Sevely, Adaptive control of biomass and substrate concentration in a continuous-flow fermentation process, RAIRO, Automatique, Productique, Informatique Industrielle 22 (1988) 159–175. [3] B. Dahhou, G. Roux, I. Queinnec, Adaptive control of a continuous fermentation process, in: IMACS Symposium, March 7–10, 1991, Proceedings 2, Lille, France, pp. 738–743. [4] J. Flaus, Estimation de l’e´tat de Bioproce´de´s a` Partir de Mesures Indirectes; The`se de Docteur de l’Institut National de Grenoble, France, 1990. [5] T. Damak, Mode´lisation, Estimation et Commande de Proce´de´s Biotechnologiques de Type Hyperbolique, The`se de Doctorat de l’universite´ Paul Sabatier, Toulouse, France, 1994. [6] C.T. Chen, C.S. Dai, Robust controller design for a class of non linear uncertain chemical processes, J. Process Contr. 11 (2001) 469– 482. [7] M. Farza, M. Nadri, H. Hammouri, approche non line´aire pour l’estimation d’e´tats et de parame`tres dans les proce´de´s chimiques et biotechnologiques; Syste`me d’information Mode´lsation, in: Optimisation Commande en Ge´nie des Proce´de´s, SIMO, Octobre 24–25, 2002, Tolouse, France. [8] W. Wu, Adaptive-like control methodologies for CSTR system with dynamic actuator constraints, J. Process Contr. 13 (2003) 525– 537. [9] G. Bastin, D. Dochain, On Line Estimation and Adaptive Control of Bioreactors, Elsevier, Amsterdam, 1990. [10] J. Villadsen, M.L. Michelsen, Solution of Differential Equation Models by Polynomial Approximation, Prentice Hall, London, 1978. [11] D. Dochain, J.P. Babary, N. Tali-Maamar, Modelling and adaptive control of nonlinear distributed parameter bioreactors via orthogonal collocation, Automatica 28 (5) (1992) 873–883.
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
1307
[12] T. Damak, J.P. Babary, M.T. Nihtila¨, Observer design and sensor location in distributed parameter bioreactors. in: Preprints, IFAC Symposium on Dynamics and Control of Chemical Reactors, Distillation Columns and Batch Processes, DYCORD+92, Maryland, USA, April 26–29, 1992, pp. 315–320. [13] J.P. Babary, S. Julien, M.T. Nihtila¨, et al., New boundary conditions and adaptive control of fixed-bed bioreactors, Chem. Eng. Process. 38 (1999) 35–44. [14] J.P. Babary, S. Bourrel, Sliding mode control of a denitrifying biofilter, Appl. Math. Modell. 23 (1999) 609–620. [15] T.T. Lee, F.Y. Wang, R.B. Newell, Robust model-order reduction of complex biological processes, J. Process Contr. 12 (2002) 807– 821. [16] M.T. Nihtila¨, T. Damak, J.P. Babary, On-line estimation of time delay via orthogonal collocation, Simul. Pract. Theor. 5 (1997) 101– 120. [17] H. Ben Bacha, M. Bouzguenda, T. Damak, et al., A study of a water desalination station using the SMCEC technique: production optimisation, Renew. Energy 21 (2000) 523–536. [18] H. Ben Bacha, T. Damak, M. Bouzguenda, et al., A methodology to design and predict operation of a solar collector for a solarpowered desalination unit using the SMCEC principle, Desalination 156 (2003) 305–313. [19] N. Tali-Maamar, T. Damak, J.P. Babary, et al., Application of a collocation method for simulation of distributed parameter bioreactors, Math. Comput. Simul. 35 (1993) 303–319. [20] Y. Zeng, B. Dahhou, M.T. Nihtila¨, et al., Microbial specific growth rate control via MRAC method, Int. J. Syst. Sci. 24 (11) (1993) 1973–1985. [21] J.L. Willems, Stability Theory of Dynamical Systems, Nelson, London, 1970. [22] Y. Zeng, B. Dahhou, M.T. Nihtila¨, Adaptive control of a nonlinear fermentation process via MRAC technique, Appl. Math. Modell. 17 (1993) 58–69.
Procedure for asymptotic state and parameter estimation of nonlinear distributed parameter bioreactors T. Damak Unite´ de Commande Automatique, Ecole Nationale d’Inge´nieurs de Sfax, Sfax BP W 3038, Tunisia Received 1 November 2004; received in revised form 1 November 2005; accepted 20 February 2006 Available online 5 July 2006
Abstract Due to measurement problem in biological variables, such as the biomass concentration and its specific growth rate, the on-line measurement of these variables is not ensured and this presents a great difficulty, especially when we consider the control problem. This paper presents and analyses an approach for estimating these biological variables through an example of a nonlinear distributed parameter bioreactor. An orthogonal collocation method is applied in the distributed system to obtain a system of ordinary differential equations. A simulation study shows the feasibility and the robustness of the estimator used for this nonlinear process. 2006 Elsevier Inc. All rights reserved. Keywords: State and parameter estimation; Nonlinear systems; Distributed parameter bioreactor; Orthogonal collocation; Asymptotic stability
1. Introduction When studying a control problem of biotechnological processes, many difficulties occur. In fact, models of these processes are in general nonlinear in character and not very accurate. Furthermore, unavailability of measurements specially the biomass concentration and its specific growth rate causes problems. To overcome these difficulties, several estimation algorithms of different type have been proposed, e.g. extended adaptive Luenberger observer [1], gradient based filter [2], modulated gain estimator [3], extended Kalman filter [4], state and parameter estimation based on sliding mode [5,6] and nonlinear observers [7,8]. One solution is proposed in this paper. The system to be studied is a fixed bed bioreactor. It is described by a set of nonlinear distributed parameter equations introduced in [5,9]. The solution used consists in: – applying a functional approximation method in order to get an ordinary differential system, in an optimal or quasi-optimal way: the orthogonal collocation method [5,10] has been used to approximate a mathematical model of hyperbolic type; E-mail address: [email protected] 0307-904X/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2006.02.014
1294
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
– deriving an estimator which gives simultaneously the estimates of variable parameters and concentrations of the bioreactor reactants; this estimator has been built from local asymptotic stability properties. The paper is organized as follows: Section 2 describes the distributed parameter process with its mathematical model. The reduction of this model to ordinary differential equations is presented in Section 3. Section 4 describes the estimation problem. The estimation procedure is developed and discussed in Section 5. The performance of the estimator is illustrated by simulation experiments in Section 6. 2. Distributed parameter system model Consider a fixed bed bioreactor with two reactions described as follows [6,11,12]: – an autocatalytic growth reaction with one limiting substrate and one biomass whose concentrations are respectively s and x: s → x;
– a death reaction of the micro-organisms: x ! xd ; where xd is the non-active biomass concentration. The micro-organisms are fixed on some support and remain within the reactor tube (Fig. 1). The substrate liquid flows through the reactor; the micro-organisms grow by consuming the substrate whose flow rate can be the control variable [13]. Assuming that the non-active micro-organisms leave the reactor and neglecting diffusion phenomenon, the above reactions can be described by the following distributed parameter system of hyperbolic type: oxðz; tÞ ¼ ½lðx; sÞ k d xðz; tÞ; ot osðz; tÞ osðz; tÞ ¼ u k 1 lðx; sÞxðz; tÞ ot oz
ð1Þ
over the space domain 0 < z 6 1 with the boundary condition s(0, t) = sin(t). In system (1), the parameters are defined as follows: l kd k1 u
specific growth rate death coefficient yield coefficient F/A (F is the feeding flow rate and A is the cross section of the bioreactor)
Though the above model is of academic type, its structure is truly representative of biotechnological processes, e.g. in the fields of biological wastewater treatment [14] or of industrial fermentations. packing
biomass
F(t)
F(t)
sin(t)
sout(t) z=0
z=1
Fig. 1. Schematic view of a fixed bed reactor.
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
1295
3. Lumped parameter system model By using the orthogonal collocation method [5,10,15] system (1) is transformed into a set of nonlinear ordinary differential equations. This functional approximation method is particularly suited for digital computation when compared with other approximation methods. It consists in expanding the model variables in terms of orthogonal polynomials Lj(z) (e.g. orthogonal Lagrange polynomials) [16–18]: xðz; tÞ ¼
N þ1 X
Lj ðzÞxj ðtÞ;
j¼0
sðz; tÞ ¼
N þ1 X
ð2Þ Lj ðzÞsj ðtÞ;
j¼0
xj(t) and sj(t) are concentrations evaluated at point z = zj. Polynomials Lj(z) have the following property: Lj ðzi Þ ¼ 1
if j ¼ i;
Lj ðzi Þ ¼ 0
if j 6¼ i:
ð3Þ
The collocation points (zj, j = 0, 1, . . . , N + 1) are chosen, as follows: 0 ¼ z0 < z1 < zN < zN þ1 ¼ 1:
ð4Þ ðp;qÞ PN
The N internal collocation points are obtained by calculating the zeros of orthogonal Jacobi polynomials [10] defined by Z 1 ðp;qÞ p zq ð1 zÞ zk P N ðzÞdz ¼ 0 k ¼ 0; 1; . . . ; N 1; ð5Þ 0
z0 is the inlet of the reactor; zN+1 is the outlet of the reactor. The determination of the number and the location of collocation points has been developed in [5,19]. In system (1), the space derivative of s(z, t) is then expressed as follows: N þ1 osðz; tÞ X dLj ðzÞ ¼ sj ðtÞ: oz dz j¼0
ð6Þ
Setting lij ,
dLj ðz ¼ zi Þ dz
ð7Þ
and applying the collocation method, the (N + 1) dimensional differential system is derived: dxi ¼ ðli k d Þxi ; dt N þ1 X dsi ¼ u lij sj uli0 s0 k 1 li xi dt j¼1
ð8Þ
ði ¼ 1; 2; . . . ; N þ 1Þ; where li , l(xi, si); s0 = sin is the boundary condition. 4. Estimation problem In what follows, we will assume that parameters kd and k1 are known and constants. On the contrary, the specific growth rate l is known to be a complex function of concentrations x and s, and of some biological parameters; a lot of mathematical models have been proposed [9] considering that l ¼ lðh; s; xÞ;
1296
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
where h is the vector of kinetic parameters; we will assume that the derivative of l with respect to h is bounded. Generally speaking, l can be considered as a time-varying parameter l(t) to be estimated in real time. The substrate concentration s(z, t), which is always less than the feeding concentration sin, and the biomass concentration are positives. The velocity u(t) is bounded (0 < u(t) 6 umax). Moreover, we will assume that we can measure only substrate concentrations and that sensors are located at collocation points. So, the output equation is written as y ¼ ½s1 ; s2 ; . . . ; sN þ1 :
ð9Þ
In general, sensors are located at different points of the reactor but not necessarily at the collocation points. In that case, the output variable y can be estimated with a good precision by using the interpolation function (2), which is needed when applying the collocation procedure. One method of estimating simultaneously the specific growth rate l and the biomass concentration x(z, t), which is the main objective of this study, is presented in the next section. 5. Estimation procedure The state estimation model is given by the following equations: d^xi ¼ ð^ li k d Þ^xi aTi es ; dt N þ1 X d^si ^i^xi uli0 sin bTi es ¼ u lij^sj k 1 l dt j¼1 ði ¼ 1; 2; . . . ; N þ 1Þ;
ð10Þ
^i are respectively the estimates of xi, si and li; es is the error on substrate concentration at where ^xi ; ^si and l each collocation point zi: eTs ¼ ½ es1 es2 . . . esN þ1 with: esi ¼ ^si si ; ai and bi are the corresponding adaptive gains. The specific growth rate in the estimation model (10) has the same structure as in the original system model: ^ ¼ lð^ l h; ^x; sÞ: ð11Þ Then, the estimation problem consists of determining gains ai and bi (i = 1, 2, . . . , N + 1) in such a way that exi ¼ ^xi xi ! 0; ^i li ! 0; eli ¼ l esi ¼ ^si si ! 0 when t ! 1:
ð12Þ
Without specifying the structure of the l function, an error system can be derived for estimating biomass concentrations and specific growth rates at collocation points, by using a first order Taylor series expansion of l [20]: ^i li ¼ HTi ehi þ K i ; exi ; eli ¼ l
ð13Þ
where eh i ¼ ^ hi h 2 Rm ; ol 2 Rm ; Hi ¼ oh h;xi ;si ol Ki ¼ 2 R; ox h;xi ;si m being the number of parameters in h.
ð14Þ ð15Þ ð16Þ
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
1297
Subtracting (8) from (10) and using (13), we get the error system: dexi ¼ ðli þ K i^xi k d Þexi aTi es þ ^xi HTi ehi ; dt N þ1 X desi ¼ k 1 ðli þ K i^xi Þexi u lij esj bTi es k 1^xi HTi ehi dt j¼1
ð17Þ
ði ¼ 1; 2; . . . ; N þ 1Þ: 5.1. Concentration estimation Defining: ai ,li þ K i^xi k d ; bi , k 1 ðli þ K i^xi Þ ¼ k 1 ðai þ k d Þ
ð18Þ
ði ¼ 1; 2; . . . ; N þ 1Þ; system (17) is written as follows: N þ1 X dexi ¼ ai exi aij esj þ ^xi HTi ehi ; dt j¼1 N þ1 X desi ¼ bi exi ðulij þ bij Þesj k 1^xi HTi ehi dt j¼1
ð19Þ
ði ¼ 1; 2; . . . ; N þ 1Þ and finally, with a specific distribution of the error vector eT ¼ ½ ex1 es1 ex2 es2 . . . exN þ1 esN þ1 ;
ð20Þ
we can write (19) in a condensed manner: de ¼ Ae þ Bxeh ; dt es ¼ Ce:
ð21Þ
• The matrix A is structured into (N + 1)2 matrices Aij of the form: ai aii Aii ¼ if i ¼ j ¼ 1; 2; . . . ; N þ 1; bi ðulii þ bii Þ " # 0 aij if i 6¼ jði; j ¼ 1; 2; . . . ; N þ 1Þ: Aij ¼ 0 ðulij þ bij Þ • B is a matrix of dimension (2N + 2) Æ (N + 1): 3 2 1 0 0 6 k 1 0 0 7 7 6 7 6 1 0 7 6 0 7 6 6 k 1 0 7 B¼6 0 7: 6 . . 7 7 6 . . . 7 6 . 7 6 4 0 0 1 5 0
0
k 1
1298
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
• x is a matrix of dimension (N + 1) Æ [m(N + 1)]: 3 2 ^x1 HT1 O O 7 6 ^x2 HT2 O 7 6 O 7 6 7 6 . . 7: 6 . . w¼6 . . 7 7 6 . .. 7 6 . . 5 4 . O x is bounded. • C is a matrix of 2 0 1 60 0 6 6. . 6 C ¼ 6 .. .. 6 6 .. .. 4. . 0 •
eTh
¼
½ eTh1
^xN þ1 HTN þ1
O
0 eTh2
dimension (N + 1) Æ (2N + 2): 3 0 0 0 0 0 1 0 07 7 .. .. 7 .. .. 7 . . 7: . . 7 .. .. 7 .. .. . .5 . . 0
0
...
eThNþ1
0 1
.
The estimation error gains ai and bi must ensure the local asymptotic stability of the dynamical error matrix A. Proposition. With an appropriate choice of aij and bij (i 5 j) the matrix A will have a block-diagonal structure; we propose: aij ¼ 0; ulij þ bij ¼ 0
ð22Þ
ði; j ¼ 1; 2; . . . ; N þ 1; i 6¼ jÞ: Choosing aij = 0 (i 5 j) is compatible with the fact that the biomass concentration at point zi depends only on the substrate concentration at this point. Choosing ulij + bij = 0 (i 5 j) corresponds to uncorrelation between the different substrate concentration errors esi ði ¼ 1; 2; . . . ; N þ 1Þ. So, the block-diagonal matrix A is structured as follows: A ¼ diag½A1 A2 Ai . . . AN þ1 with
ai Ai ¼ bi
aii ulii bii
i ¼ 1; 2; . . . ; N þ 1:
ð23Þ
If each sub-matrix Ai is asymptotical stable, the whole error system (21) will also be asymptotical stable. The characteristic polynomial of Ai is given by P i ðkÞ ¼ k2 kðai ulii bii Þ þ bi aii ai ðulii þ bii Þ:
ð24Þ
The gains aii and bii must be bounded; assuming that their variations are also bounded, then considering the i assumptions given in Section 4, the sub-matrix Ai and its time derivate dA are bounded. The asymptotic stadt bility is obtained if the following conditions hold [21]: ri ¼ ai ulii bii < 0 ðsum of eigenvaluesÞ; pi ¼ bi aii ai ðulii þ bii Þ > 0 ðproduct of eigenvaluesÞ:
ð25Þ
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
1299
Choosing two eigenvalues k1i and k2i with negative real part, ri and pi are then fixed. So, from (25) we get: bi aii ¼ pi þ ai ðai ri Þ: Since bi = k1(ai + kd) we obtain the gains aii and bii: aii ¼
pi ai ri þ a2i k 1 ðai þ k d Þ
bii ¼ ai ulii ri
i ¼ 1; 2; . . . ; N þ 1;
i ¼ 1; 2; . . . ; N þ 1:
ð26Þ ð27Þ
The expression (26) can also be written as follows: aii ¼
1 ðai k1i Þðai k2i Þ : k1 ai þ k d
ð28Þ
Taking into account the proposition (22), the formulation of the state estimation model (10) is simplified: d^xi ¼ ð^ li k d Þ^xi aii ð^si si Þ; dt N þ1 X d^si ^i^xi uli0 sin þ ¼ ðulii þ bii Þ^si k 1 l bij sj dt j¼1
ð29Þ
ði ¼ 1; 2; . . . ; N þ 1Þ: Gains aii are given by (26) or (28); gains bij are determined by (22) and (27). Moreover, the error system (19) becomes: dexi ¼ ai exi aii esi þ ^xi HTi ehi ; dt desi ¼ bi exi ðai ri Þesi k 1^xi HTi ehi : dt
ð30Þ
5.2. Parameter estimation The resolution of the state estimation model (29) needs the knowledge of the estimate of the specific growth rate l depending on the parameter vector h at each point zi. Defining xi ¼ ^xi Hi , the system (30) can be written as follows: dexi ¼ ai exi aii esi þ xTi ehi ; dt desi ¼ bi exi ðai ri Þesi k 1 xTi ehi : dt
ð31Þ
The parameter estimation should be derived from the knowledge of the error on the substrate concentration estimate: d^ hi ¼ ci ðtÞesi dt
i ¼ 1; 2; . . . ; N þ 1;
ð32Þ
where the gain vector ci has to be determined. One solution consists in calculating ci in such a way that the augmented error vector ½ exi esi eli tends to ^ i li . zero, when time tends to infinity (asymptotic stability), with eli ¼ l In the following, we are considering two cases according as l is a function of one parameter or more. If we assume that h is constant or slowly variable, then: _ _ e_ hi ¼ ^ hi h_ ffi ^ hi :
1300
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
5.2.1. l is a function of one parameter h The whole error differential system is written as 2 3 2 32 3 e_ xi exi ai aii xi 6 7 6 76 7 4 e_ si 5 ¼ 4 bi ðai ri Þ k 1 xi 54 esi 5; 0 ci 0 ehi e_ hi
ð33Þ
whose characteristic polynomial is given by P i ðkÞ ¼ k3 ri k2 þ ðk 1 ci xi þ pi Þk þ k 1 k d ci xi :
ð34Þ
The asymptotic stability is ensured if the following conditions hold (ri being <0): ð35Þ ð36Þ
k 1 k d ci xi > 0; ri ðk 1 ci xi þ pi Þ k 1 k d ci xi > 0: In (35), xi can be positive or negative; if we choose ci ðtÞ ¼ Ci xi ðtÞ ¼ Ci^xi Hi
ð37Þ
with Ci positive, then (35) is always fulfilled. Choosing an adaptive gain Ci given by Ci ¼ x12 , yields to a constant coefficient characteristic polynomial: i
P i ðkÞ ¼ k3 ri k2 þ ðk 1 þ pi Þk þ k 1 k d :
ð38Þ
Then the condition (36) becomes: ri ðk 1 þ pi Þ k 1 k d > 0;
ð39Þ
which holds with an appropriate choice of ri and pi. One can notice that (39) is fulfilled by choosing ri < kd, which corresponds to a realistic condition, since the dynamics of the concentration estimator must be sufficiently rapid with respect to the system dynamics. 5.2.2. l is a function of several parameters The whole error differential system is now of the form: 32 3 2 3 2 exi e_ xi ai aii xTi 76 7 6 7 6 T 4 e_ si 5 ¼ 4 bi ðai ri Þ k 1 xi 54 esi 5; ehi e_ hi 0 ci 0
ð40Þ
whose characteristic polynomial is given by P i ðkÞ ¼ km1 ½k3 ri k2 þ ðk 1 xTi ci þ pi Þk þ k 1 k d xTi ci ;
ð41Þ
m being the dimension of h. The asymptotic stability of the whole system is no more ensured because there exist null eigenvalues (weak stability). Conditions (35) and (36) are now generalized: k 1 k d xTi ci > 0;
ri ðk 1 xTi ci
þ pi Þ
ð42Þ k 1 k d xTi ci
> 0:
ð43Þ
As before, the first condition (42) is fulfilled if: ci ðtÞ ¼ Ci xi ðtÞ;
ð44Þ
where Ci is a diagonal positive definite matrix. To obtain a characteristic polynomial with constant coefficients, we can use an adaptive gain Ci defined by 1 1 1 Ci ¼ diag 2 2 2 : ð45Þ xi1 xi2 xim
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
1301
The characteristic polynomial is then given by P i ðkÞ ¼ km1 ½k3 ri k2 þ ðmk 1 þ pi Þk þ mk 1 k d :
ð46Þ
The condition (43) becomes: ri ðmk 1 þ pi Þ mk 1 k d > 0;
ð47Þ
which also holds with an appropriate choice of ri and pi. It can be shown that the error on concentration estimates tends to zero; the error on the parameter ^hi does ^i tends to not diverge, but ^ hi tends in general towards wrong values, due to the null eigenvalues. Nevertheless l the true value of li, because, when exi and esi tend asymptotically to zero, we can deduce from (30) that ^xi HTi ehi ! 0 ^i li ¼ HTi ehi þ K i exi ! 0, when t ! 1. then, eli ¼ l In conclusion, the parameter estimator is given by the following differential system: d^ hi ¼ Ci^xi Hi esi ; dt
ð48Þ
where Ci is a diagonal positive definite matrix defined by (45). With this choice, the only tuning parameters of the resulting adaptive state and parameter estimator are ri and pi. The above results depend on the true values of li(xi, si, hi), Ki(xi, si, hi) and Hi(xi, si, hi). Since these true values are supposed to be unknown, they should be substituted by the corresponding estimated values b i ð^xi ; si ; ^ b i ð^xi ; si ; ^ ^i ð^xi ; si ; ^ l hi Þ; K hi Þ and H hi Þ. This substitution can be justified by using the L2-robustness theorem of Zeng et al. [22]. 6. Simulation results Let us consider for example, the specific growth rate given by the Contois model: lðh; x; sÞ ¼
lm s ; kcx þ s
ð49Þ
T
where h ¼ ½ lm k c is a constant parameter vector; lm denotes the maximum specific growth rate and kc is the saturation constant. In this case, we have: " #T ^ s s x i i i bi ¼ ^ mi H l ; 2 ^k c ^xi þ si ð^k ci ^xi þ si Þ i ð50Þ ^k c si i b i ¼ ^ K lm i : ð^k c ^xi þ si Þ2 i
6.1. Numerical data For simulation needs, parameters were set to the following numerical values [11]: A ¼ 0:02 m2 ;
k d ¼ 0:05 h1 ;
k 1 ¼ 0:4
lm ¼ 0:35 h1 ;
k c ¼ 0:4:
For the process simulation, and then generating data, we have used a model different from the lumped model given by the orthogonal collocation. This model is based on a high-dimensional finite difference method. The scheme used (central difference) is the following: oX X iþ1 X i1 : ð51Þ ¼ oz z¼zi 2Dz
1302
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
To have a good accuracy, an appropriate space step size Dz, equal to 0.002, is used corresponding to 500 space intervals. For the lumped model used in the estimator, four internal collocation points are defined. They were optimally located (p = 0; q = 4) [5] at ½ 0:312
0:579
0:813
0:963 :
The initial conditions for the process are chosen so as to correspond to a steady-state with F = 2 l/h, and sin(0) = s0(0) = 5 g COD/l. At the collocation points we get the following initial values: sð0Þ ¼ ½ 1:96
0:88
xð0Þ ¼ ½ 29:4 13:2
0:44
0:28
6:55 4:18
0:25 T in g COD=l; T
3:73 in g VSS=l:
The generating data, i.e., substrate measurements, are made with a sampling period of Te = 0.1 h. To obtain more realistic results, a zero mean-valued Gaussian noise with a standard deviation 0.005 is added into the measurement. The numerical solution method used to solve the set of ordinary differential equations is a fourth order Runge–Kutta method with a fixed time step equal to 0.01 h. 6.2. Simulation runs The practical adaptive state and parameter estimator given above is studied via simulations. The tuning parameters of the estimator are setting to the values ri = 5.05 and pi = 0.25. For testing the estimator convergence and robustness properties, simulation experiments, in which the input variable u and the parameters lm and kc are modified (step functions), are considered. Figs. 2–5 represent respectively the evolution, at the outlet of the reactor (point z5 = 1.0), of: – – – –
the the the the
substrate concentration s and its estimate, biomass concentration x and its estimate, specific growth rate l and its estimate, parameters lm and kc and their estimates,
Fig. 2. The evolution of substrate concentration at z = 1.0 (kc and lm are unknown).
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
1303
Fig. 3. The evolution of biomass concentration at z = 1.0 (kc and lm are unknown).
Fig. 4. The evolution of specific growth rate at z = 1.0 (kc and lm are unknown).
corresponding to the following changes: – – – – – – –
at at at at at at at
time time time time time time time
t = 25 h: flow rate step from 2 to 2.5 l/h, t = 50 h: flow rate step from 2.5 to 2 l/h, t = 75 h: flow rate step from 2 to 2.5 l/h, t = 100 h: flow rate step from 2.5 to 2 l/h; lm changes from 0.35 to 0.42 h1, t = 125 h: flow rate step from 2 to 2.5 l/h, t = 150 h: flow rate step from 2.5 to 2 l/h, t = 175 h: flow rate step from 2 to 2.5 l/h; kc changes from 0.4 to 0.48.
As can be seen from the figures, we can verify that the concentration estimates, i.e., the state variables of the system (s, x), tend in a satisfactory way to the actual values in spite of perturbation steps and noisy measurements.
1304
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
Fig. 5. The evolution of biological parameters at z = 1.0 (kc and lm are unknown).
Fig. 6. The evolution of substrate concentration at z = 0.312 (kc is known, lm unknown).
^ tends also to the actual value although the biological parameters l ^m and ^k c do The specific growth rate estimate l not converge to the true values as shown before in the case where l is a function of several parameters. This unidentifiability of the system parameters can be explained in another way [22] by considering the structure of the specific growth rate of Contois’s model lm s lðh; x; sÞ ¼ : ð52Þ kcx þ s For constant substrate and biomass concentrations there are many values lm and kc that give the same value l for the specific growth rate. Consequently, depending on the initial values of the parameters, estimates converge to some feasible values that give the true value to the specific growth rate.
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
1305
Fig. 7. The evolution of biomass concentration at z = 0.312 (kc is known, lm unknown).
Fig. 8. The evolution of specific growth rate at z = 0.312 (kc is known, lm unknown).
A second simulation experiment is done in the case where one of the two system parameter estimates is assumed to be known (e.g. kc). The results of the corresponding run are shown in Figs. 6–9 at point z1 = 0.312. The realized changes are identical to the ones given above except at time t = 175 h where the step change of kc is not considered. ^m ) converges to its correct value. This We can notice, in this case, that the parameter to be estimated (i.e. l can be seen in Fig. 9. It is possible to reconstruct the evolution of concentrations x(z, t) and s(z, t) at any point of the reactor by using expressions (2).
1306
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
Fig. 9. The evolution of biological parameters at z = 0.312 (kc is known, lm unknown).
7. Conclusion In this paper, it has been shown how it is possible to design an asymptotic estimator of state and time-varying parameters in the case of a nonlinear distributed parameter bioreactor. The objective was to have good estimates of the biomass concentration and its specific growth rate from measurements of the substrate concentration. The structure of the estimator is based upon an approximated model of the bioreactor behaviour. The excellent properties of the estimator, especially robustness to parameter variations, model mismatch and measurement errors, were tested through simulations. The results obtained show that the concentration estimates converge in a good way to the true values. The specific growth rate (time-varying parameter) is fairly well estimated, even if it is a function of several unknown parameters, which should be identified if needed. References [1] G. Bastin, D. Dochain, Adaptive estimation of microbial growth rate, in: H.A. Barker, P.C. Young (Eds.), Preprints 7th IFAC/ IFORS Symposium on Identification and System Parameter Estimator, York, United Kingdom, July 3–7, vol. 2, Pergamon Press, 1988, pp. 1161–1166. [2] G. Chamilothois, Y. Sevely, Adaptive control of biomass and substrate concentration in a continuous-flow fermentation process, RAIRO, Automatique, Productique, Informatique Industrielle 22 (1988) 159–175. [3] B. Dahhou, G. Roux, I. Queinnec, Adaptive control of a continuous fermentation process, in: IMACS Symposium, March 7–10, 1991, Proceedings 2, Lille, France, pp. 738–743. [4] J. Flaus, Estimation de l’e´tat de Bioproce´de´s a` Partir de Mesures Indirectes; The`se de Docteur de l’Institut National de Grenoble, France, 1990. [5] T. Damak, Mode´lisation, Estimation et Commande de Proce´de´s Biotechnologiques de Type Hyperbolique, The`se de Doctorat de l’universite´ Paul Sabatier, Toulouse, France, 1994. [6] C.T. Chen, C.S. Dai, Robust controller design for a class of non linear uncertain chemical processes, J. Process Contr. 11 (2001) 469– 482. [7] M. Farza, M. Nadri, H. Hammouri, approche non line´aire pour l’estimation d’e´tats et de parame`tres dans les proce´de´s chimiques et biotechnologiques; Syste`me d’information Mode´lsation, in: Optimisation Commande en Ge´nie des Proce´de´s, SIMO, Octobre 24–25, 2002, Tolouse, France. [8] W. Wu, Adaptive-like control methodologies for CSTR system with dynamic actuator constraints, J. Process Contr. 13 (2003) 525– 537. [9] G. Bastin, D. Dochain, On Line Estimation and Adaptive Control of Bioreactors, Elsevier, Amsterdam, 1990. [10] J. Villadsen, M.L. Michelsen, Solution of Differential Equation Models by Polynomial Approximation, Prentice Hall, London, 1978. [11] D. Dochain, J.P. Babary, N. Tali-Maamar, Modelling and adaptive control of nonlinear distributed parameter bioreactors via orthogonal collocation, Automatica 28 (5) (1992) 873–883.
T. Damak / Applied Mathematical Modelling 31 (2007) 1293–1307
1307
[12] T. Damak, J.P. Babary, M.T. Nihtila¨, Observer design and sensor location in distributed parameter bioreactors. in: Preprints, IFAC Symposium on Dynamics and Control of Chemical Reactors, Distillation Columns and Batch Processes, DYCORD+92, Maryland, USA, April 26–29, 1992, pp. 315–320. [13] J.P. Babary, S. Julien, M.T. Nihtila¨, et al., New boundary conditions and adaptive control of fixed-bed bioreactors, Chem. Eng. Process. 38 (1999) 35–44. [14] J.P. Babary, S. Bourrel, Sliding mode control of a denitrifying biofilter, Appl. Math. Modell. 23 (1999) 609–620. [15] T.T. Lee, F.Y. Wang, R.B. Newell, Robust model-order reduction of complex biological processes, J. Process Contr. 12 (2002) 807– 821. [16] M.T. Nihtila¨, T. Damak, J.P. Babary, On-line estimation of time delay via orthogonal collocation, Simul. Pract. Theor. 5 (1997) 101– 120. [17] H. Ben Bacha, M. Bouzguenda, T. Damak, et al., A study of a water desalination station using the SMCEC technique: production optimisation, Renew. Energy 21 (2000) 523–536. [18] H. Ben Bacha, T. Damak, M. Bouzguenda, et al., A methodology to design and predict operation of a solar collector for a solarpowered desalination unit using the SMCEC principle, Desalination 156 (2003) 305–313. [19] N. Tali-Maamar, T. Damak, J.P. Babary, et al., Application of a collocation method for simulation of distributed parameter bioreactors, Math. Comput. Simul. 35 (1993) 303–319. [20] Y. Zeng, B. Dahhou, M.T. Nihtila¨, et al., Microbial specific growth rate control via MRAC method, Int. J. Syst. Sci. 24 (11) (1993) 1973–1985. [21] J.L. Willems, Stability Theory of Dynamical Systems, Nelson, London, 1970. [22] Y. Zeng, B. Dahhou, M.T. Nihtila¨, Adaptive control of a nonlinear fermentation process via MRAC technique, Appl. Math. Modell. 17 (1993) 58–69.