In this paper an adaptive model-based algorithm is proposed for the on-line
estimation of reaction rates in stirred tank bioreactors. The main design condition imposes that the observation errors reflecting the mismatch between the estimated parameters and the 'true' values follow second-order dynamics of convergence. The gain matrices are shown to be functions of the state and of user-defined damping coefficients and natural periods of oscillation for second-order trajectories. The application of the algorithm is illustrated with a simple case-study involving the estimation of the specific reaction rate for a single substrate,
single product scheme.Junta Nacional de Investigação Científica e Tecnológica (JNICT)

Azevedo, S. Feyo de

Ferreira, E. C.

Oliveira, Rui Manuel Freitas

Ferreira, Eugénio C.

Universidade do Minho: RepositoriUM

  Copyright  IFAC Computer Applications in Biotechnology, Garmisch-Partenkirchen, Germany, 1995 259 MODEL-BASED ESTIMATION OF REACTION RATES  IN STIRRED TANK BIOREACTORS   R. Oliveira1, E.C. Ferreira2, S. Feyo de Azevedo1*   1Departamento de Engenharia Química, Faculdade de Engenharia da Universidade do Porto Rua dos Bragas, 4099 Porto Codex, PORTUGAL (e-mail: sfeyo@garfield.fe.up.pt) 2Departamento de Engenharia Biológica, Escola de Engenharia, Universidade do Minho Campus de Gualtar, 4709 Braga, PORTUGAL (e-mail: eugenio@deb.uminho.pt)     Abstract: In this paper an adaptive model-based algorithm is proposed for the on-line estimation of reaction rates in stirred tank bioreactors. The main design condition imposes that the observation errors reflecting the mismatch between the estimated parameters and the 'true' values follow second-order dynamics of convergence. The gain matrices are shown to be  functions of the state and of user-defined damping coefficients and natural periods of oscillation for second-order trajectories. The application of the algorithm is illustrated with a simple case-study involving the estimation of the specific reaction rate for a single substrate, single product scheme.  Keywords: Estimator; Adaptive; Convergence; Second-order systems; Kinetics    1. THE ESTIMATION PROBLEM  Model-based state observation and parameter estimation represent fundamental tools for monitoring and control of biotechnological processes. Dochain and Bastin (1990) established a general theoretical framework for the analysis of bioreactor dynamics. The concept of a general state space dynamical model and that of minimal kinetic modelling were proposed  ddtKH D Fξ ξ ρ ξ ξ= − + −( ) ( ) Q  (1)  where ξ is the state vector (the set of n component concentrations), K an (n×m) yield coefficients matrix, D the dilution rate, F the feed rate vector with dim(F)=n and Q the gaseous outflow rate vector with dim(Q)=n. The reaction rates are defined as ϕ(ξ)=Hξ)ρ(ξ) to take advantage of any possible knowledge of the kinetic model, being H(ξ) a m×r matrix of known functions of the state and ρ(ξ) a vector of r unknown functions of the state.   This general model constitutes the core element for the design of state observers and parameter estimators. A wide class of problems is covered, depending upon the degree of knowledge of the process model. The reaction rates ϕ are most often complex functions of the operating conditions and of the state of the process. Situations like 'unknown (some or all) yield coefficients', or 'unknown or partially known kinetic models' are extensively studied and discussed in the literature (Bastin and Dochain, 1990; Pomerleau and Perrier, 1990; Feyo de Azevedo et al., 1992). The construction of a suitable kinetic model may constitutes a difficult task, if not an impossible one. As such, there is a clear incentive to design estimators of the unknown functions ρ(ξ).  The algorithm presented below addresses a class of estimation problems which can be defined in the following four points -   i. The biotechnological process can be described by the state space dynamical model (1); ii. The yield coefficients (matrix K) are known and invariant; iii. The dilution rate D, the feed rates F and the  260 gaseous outflow rates Q are known on-line; iv. The state variables are known on-line either by measurement or by means of a state observer.  As such, the scope of the algorithm will be the on-line estimation of ρ(ξ) from the on-line knowledge of D, F, Q and ξ.  For this general problem observer-based estimators were proposed (Bastin and Dochain, 1990; Pomerleau and Perrier, 1990), involving trial and error procedures or employing pole placement arguments for the tuning of the regression gain matrices. The latter, albeit with a theoretical basis, tends to be a sensitive technique. In a previous paper (Oliveira et al., 1994) a tuning procedure has been discussed for the gain matrices of one such algorithm which aimed at imposing a second-order dynamics of convergence of the estimated reaction rates to their true values. It was shown that the characteristic parameters of the second-order responses obtained, i.e. the natural period of oscillations τi and the damping coefficients ζi were functions of the state variables, hence time-varying. In this work an adaptive scheme is proposed which imposes second-order convergence with constant user-defined damping coefficients and periods of oscillation, while obeying the stability requirements dictated by a Lyapunov's analysis (Bastin and Dochain, 1990; Narendra and Annaswamy, 1989).   2. THE OBSERVER-BASED ESTIMATOR  The proposed kinetics estimator is based on a partial-state observer which provides estimates of 'r' state space variables (ξ1) and on an additional equation for the updating of the 'r' parameter estimates. The observation error ( )ξ ξ1 1− !  which is supposed to reflect the mismatch between the estimated parameters and their true values, is used as the driven force in the updating law. The two relevant equations are stated as follows: ddtK H D F Q!( ! , ) ! ! ( ! )ξ ξ ξ ρ ξ ξ ξ1 1 1 2 1 1 1 1 1 1= − + − − −Ω (2a) ddt!( ! )ρ ξ ξ= −Ω2 1 1  (2b)  where !ξ1  denotes the on-line estimate of ξ1 and !ρ  the on-line estimate of ρ(ξ). Ω1 and Ω2 are square (rxr) tuning matrices for the control of stability and tracking properties of the algorithm. The defining equations for these matrices are the result of stability considerations and of the established design conditions (eqns. 17a and 17b below).   3. STABILITY ANALYSIS  The dynamics of the observation error ( )ξ ξ1 1− !  and of the tracking error ( )ρ ρ− !  are obtained by subtracting eqn. (2a) from eqn. (1) and by rearranging eqn. (2b), leading to the following non-linear error system: [ ]d dt K H H D( ! )( , ) ( ! , ) ! ( ! ) ( ! )ξ ξ ξ ξ ρ ξ ξ ρ ξ ξ ξ ξ1 1 1 1 2 1 2 1 1 1 1 1− = − − − − −Ω (3a) ddtddt( ! )( ! )ρ ρ ξ ξ ρ− = − − +Ω2 1 1  (3b) where "ρ  is considered as an external perturbation.  The point ξ ξ= ! 1  and ρ ρ= !  is an equilibrium point of the unperturbed system. Stating that Ω2 is independent of ρ, a linear approximation of the unperturbed system around such point gives: ( )ddtC K H( ! )( ! , , ! ) ( ! ) ( ! , )( ! )ξ ξ ξ ξ ρ ξ ξ ξ ξ ρ ρ1 1 1 2 1 1 1 1 1 2− = − − + −Ω (4a) ddt( ! ) ( ! )ρ ρ ξ ξ− = − −Ω2 1 1  (4b) being C( ! , , ! )ξ ξ ρ1 2  defined by [ ]C K H DI N(! , , ! ) ( , )!!ξ ξ ρ ∂ ξ ξ ρ∂ξ ξ ξ1 2 11 211 1=−=(4c)  The relevant coefficients matrix for the analysis of stability is: A C K H= −−(! , , ! ) (! , )ξ ξ ρ ξ ξ1 2 1 1 1 22 0ΩΩ (4d)  From the direct method of Lyapunov, it follows that the unperturbed system is exponentially stable if -  Condition 1 - the eingenvalues of matrix A have strictly negative real parts.  With the convenient and consistent definition of matrices Ω1 and Ω2 as:  Ω Ω1 = +C  (5a) ( )[ ]Ω Γ2 1 1 2 1= −K H ! ,ξ ξ  (5b)  matrix A can then be rewritten as: ( )( )[ ]AK HK H=−−−ΩΓ1 1 21 1 210! ,! ,ξ ξξ ξ  and the eingenvalues of A (λi, i=1,...,2r) are obtained  261 by solving in order to λ the polynomial equation    P2r(λ)=0, represented by: ( )( )[ ]− −− −=−ΩΓλ ξ ξξ ξ λI K HK H Irr1 1 21 1 21 0! ,! , (6)  A decoupled solution of (6) is possible by defining Ω and Γ as: Ω = ωIr  (7a) Γ = γIr  (7b)  then, eqn.(6) can be decoupled in "r" eqns.: ( )λ ω λ γ i  i+ + = 0       i=1,...,r (8) which accepts the solutions: λ ω ω γi =− ± −21242  (9) If the parameter ω is a real positive constant, then Condition 1 is obeyed, hence the unperturbed non-linear error system (3) is exponential stable  In addition, the perturbed system is globally stable (i.e. the output error is bounded for all t) if "ρ  is a continuously differentiable bounded function. The conditions under which this is verified were established by Dochain and Bastin (1990) -  Condition 2 - The dilution rate is bounded below: 0 < ≤ ∀D D t tmin ( )  (10)  Condition 3 - The feed rates are bounded: 0 < ≤ ∀ ∀F t F ti i t( ) ( )  (11)  Condition 4 - Each reaction involves at least one reactant that is neither a catalyst nor an autocatalyst.  Condition 5 - ρ(ξ) is a differentiable function of ξ.  It is thus clear that global stability is related to the prevalent experimental conditions.   4. TUNING AND DYNAMICS OF CONVERGENCE  The error system eqns. (3a) and (3b) is the basis for the analysis of convergence. Under the restriction of employing the linear approximation of this full error system, around ξ ξ= ! 1  and ρ = ρ^, it is possible to prove that second order convergence with constant parameters can be obtained.   Substituting eqns. (5a) and (5b) (for Ω1 and Ω2) in such linearized system gives: ddtK H( ! )( ! ) ( ! )ξ ξ ξ ξ ρ ρ1 1 1 1 1− = − − + −Ω (12a) [ ] ( )1111 ˆˆ ξξρ −Γ= −HKdtd   (12b)  Differentiating eqn. (12b) and combining with (12a) results in -   [ ] [ ] [ ] [ ]Γ ΩΓ− −+ + =1 12211 1 1K HddtK HddtK H K H! !!ρ ρ ρ ρ (13)  with definitions (7), eqn. (13) simplifies to 1 22γρ ωγρ ρ ρddtddt! ! !+ + =  (14)  Setting ω and γ as ωζτ=2 (15a) γτ=12  (15b)  then, eqn. (14) can be rewritten as: τρ ζτ ρ ρ ρ222 2ddtddt! ! !+ + =  (16)  where τ and ζ represent respectively the well known natural periods of oscillation and damping coefficients which characterise second order trajectories.   Finally, combining eqns. (5a), (5b), (7a), (7b), (15a) and (15b) leads to the setting of matrices Ω1 and Ω2 as: ( ) ( )Ω1 1 2 1 2 2! , , ! ! , , !ξ ξ ρ ξ ξ ρ ζτ= +C Ir  (17a) ( )[ ]Ω2 1 2 2 1 1 2 11(! , ) ! ,ξ ξ τ ξ ξ= −K H  (17b)   The adaptive observer-based kinetics estimator is thus constituted by eqns. (2a) and (2b), with the tuning matrices Ω1 and Ω2 defined by eqns. (17a) and (17b). Under the constraints of the validity of the linearized tangent observation and tracking error model, this scheme leads to a second-order dynamics of convergence of ρ^ to ρ(ξ). The only other condition is that matrix K H1 1 2(! , )ξ ξ  is full rank. This is in principle attained through an appropriate partition of the system state.        262 5. PROCEDURE IMPLEMENTATION  The general procedure implementation can be summarised in the following way:  (i) Derive appropriate process dynamical model and state observer; (ii) Choose suitable state partition ξ1 to satisfy the theoretical constraints; derive the observer-based kinetics estimator; choose a suitable sampling time T (problem dependent); (ii) At sampling time nT: read available information and compute the state either solely by direct measurement or by additionally applying a state observer procedure; (iii) Compute a linear interpolation function of the state between times (n-1)T and nT; (iv) Integrate eqns. (2), employing a finite difference 4th/5th order variable step integration scheme. At each internal integration step, update the regression gains as stated by eqns. (17), using where required the linearly interpolated values of step (iii)   The theoretical development and the numerical implementation are based in continuous-time equations. The switch to discrete time versions (Euler's approach) as employed by other authors (Bastin and Dochain, 1990, Pomerleau and Perrier, 1990) would raise specific stability questions in which the sampling period would play an important role. In particular, it can be shown that the sampling period would have an upper bound function of the tuning parameters ζ and τ. The requirement of too frequent measurements may pose practical difficulties. With today's availability of computational power there seems to be no reason to avoid higher order and more complex integration schemes. More detailed discussion is out of the scope of this work   6. A CASE-STUDY  5.1 The fermentation process  The use of algorithm (2) is illustrated through a simple application: the estimation of a microbial specific reaction rate in a simple biological culture which involves a single biomass (X) growing on a single substrate (S) and yielding a single product (P).   The reaction scheme is stated as follows: S X Pϕ → +  (18)  The process dynamics in a fed-batch fermenter is described as ddtXSPkkDXSPDSin= −−+1 0012ϕ  (19)  where D is the dilution rate (D=Fo/V), being Fo the input flow rate, V the broth volume in the fermenter) and Sin the substrate concentration in the feed. The reaction rate may be defined as:  ϕ = XSα (20) being α the unknown kinetic term.  Typical experimental conditions were: Sin= 10 g.l-1, V(0)= 3 l, X(0)= 0.1 g.l-1, S(0)= 0.1 g.l-1, P(0)=0.01 g.l-1, k1= 2.85 g gS X.−1 , k2= 1.55 g gP X.−1 . The input flow rate Fo was a forced square wave between 0.01 and 0.02 l.hr-1.   5.2 The observer-based estimator  In the present case the objective is the on-line estimation of the time-varying reaction rate term α from the on-line knowledge of X, S, P, Sin, V and Fo.  The application of algorithm (2) with ξ1=X leads to the following two eqns.: dXdtXS DX X X! ! ! ! ( ! )= − − −α ω1  (21a) ddtX X!( ! )α ω= −2  (21b)  From the defining eqns. (17a) and (17b), ω1 and ω2 are given by ω αζτ12= − +S D!  (22a) ωτ2 21= !XS (22b)   7. RESULTS AND CONCLUSIONS  To test the capability of the algorithm at imposing a second order dynamics of convergence, the 'true' specific rate has been assumed to be a square wave signal. The process has been implemented in a process simulator and data taken at sampling times of T=6 minutes.  Running the estimator with constant values of ω1 and ω2 (in eqns. 21), i.e. not employing the defining eqns. (22) for the adaptive algorithm, would lead to a variable second order convergence (Oliveira et al., 1994). Tests were performed, employing the adaptive scheme (eqns. 22), but with negative weights for  263 either ζ or τ. As expected, the estimator diverged, since such combinations lead to negative values of ω (eqn.15a), which violates stability condition 1 discussed above.  Finally, under the same operating conditions, the estimator has been tested for a range of positive values for the damping factor and for the natural period of oscillation. Illustrative results are presented in Figs. 1-6. The square wave line is the 'true' specific rate parameter, whereas the full line represents its on-line estimate. The accuracy of the estimates can be assessed from the ITAE error index (ITAE - integral of time-weighted absolute errors) given in the legend. The dynamics of convergence of α^ to α exhibits characteristics which are in full agreement with typical second order dynamic response to step inputs, viz.-  (a) increasing ζ, the response turns to be less oscillatory;  (b) decreasing τ the response becomes faster;  (c) ζ=1 defines the frontier between oscillatory and non-oscillatory behaviour;  (d) low values of ζ give rise to high peak responses. As expected, a choice of the damping factor around 0.8, together with a low natural period of oscillation, leads to good convergence.   Overall, results are good and coherent with the theory. The choice of tuning parameters has an intuitive and theoretically simple basis, since this type of second order response is widely observed in natural phenomena and well studied.   Further theoretical analysis is out of the scope of this paper. Work is in progress which aims in particular at testing the robustness and at establishing the domains of validity of the procedure proposed.   REFERENCES  Bastin, G., Dochain, D. (1990). On-Line Estimation and Adaptive Control of Bioreactors. Elsevier, Amsterdam. Feyo de Azevedo, S., Pimenta, P., Oliveira, F., Ferreira, E.C. (1992). Studies on On-Line State and Parameter Estimation through a Real-Time Process Simulator. In: Modelling and Control of Biotechnical Processes, (M.N. Karim, G. Stephanopoulos, Ed.), pp. 453-458, Pergamon Press, New York. Narendra, K.S., Annaswamy, A.M. (1989). Stable Adaptive Systems. Prentice-Hall Inc., London. Oliveira, R., Ferreira, E.C., Oliveira, F., Feyo de Azevedo, S. (1994). A Study on the Convergence of Observer-Based Kinetic Estimators in Fed-Batch Fermentations. In: The 5th Int. Symp. on Process Systems Engineering - PSE'94, (E.S. Yoon, Ed.), pp. 683-688, Kyongju, Korea. Pomerleau, Y. and Perrier, M. (1990). Estimation of Multiple Specific Growth Rates in Bioprocesses. AIChE J., 36, 207-215.   Acknowledgements - This work was partially supported by JNICT - Junta Nacional de Investigação Científica e Tecnológica, under contracts number 0064/B/91 (Instituto de Sistemas e Robótica) and BD/2501/93-RM.  264 -0.21.20 20t (hr)α(l/gh) Fig. 1. Results obtained with τ = 0.15  and ζ = 0.25 (ITAE = 14.0).  -0.21.20 20t (hr)α(l/gh) Fig. 3. Results obtained with τ = 0.15  and ζ = 1.0 (ITAE = 9.9).  -0.21.20 20t (hr)α(l/gh) Fig. 5. Results obtained with τ = 0.1  and ζ = 0.8 (ITAE = 4.6).   -0.21.20 20t (hr)α(l/gh) Fig. 2. Results obtained with τ=0.15  and ζ = 0.5 (ITAE = 8.1).  -0.21.20 20t (hr)α(l/gh) Fig. 4. Results obtained with τ = 0.15  and ζ = 1.5 (ITAE = 16.0).  -0.21.20 20t (hr)α(l/gh) Fig. 6. Results obtained with τ = 0.01  and ζ = 0.8 (ITAE = 0.73).  

Model-based estimation of reaction rates in stirred tank bioreactors

http://repositorium.sdum.uminho.pt/bitstream/1822/3760/1/CAB6%5b1%5d.pdf

Model-based estimation of reaction rates in stirred tank bioreactors

Abstract

Similar works

Full text

Available Versions

Universidade do Minho: RepositoriUM