Bioprocess Engineering 5 (1990) 225-230

Bi0pr0cessEngineering 9 Springer-Verlag 1990

Control in bioreactors showing gradients S. R. Weijers, G. H o n d e r d a n d K. Ch. A. M . Luyben, D e l f t

Abstract. In large-scale bioreactors gradients often occur as a result of non-ideal mixing. This phenomenon complicates design and control of large-scale bioreactors. Gradients in the oxygen concentration can be modeled with a two-compartment model of the liquid phase. Application of this model had been suggested for the control of the dissolved oxygen concentration with a batch gluconic acid fermentation process as the model system. The control system consists of a controller, an observer and a parameter estimator. In this work, the controller design is reconsidered and, in simulation experiments, the performance of the control system has been investigated, When the parameter values are known, the controller in combination with the observer works adequately. The parameter estimator, however, yields incorrect parameters, which are caused by a coupling between two parameters. This causes a deviation of the estimated states from the process states. The simulation results suggest that a priori knowledge of the parameters is required for application of the model for control and state estimation.

List of symbols a I

aco, in C1 C2 Ce Cg Cg, i,, C~ H k' kc

m o l m -a mol m - 3 mol m - 3 mol m - 3 mol m a C - e q m -3

ke

s- 1

kid OTR OUR Pl P2 Po V1 V2 Vr Vo Yo~

s-t molm as-1 m o l m -3 s -1

Z Zc

m3 m3 m3 m3 C - e q mol l

coefficient for computation of set point for C1 coefficient for computation of feedforward DO concentration in aerated part DO concentration in non-aerated part electrode response oxygen concentration in exhaust gas oxygen concentration inlet gas biomass concentration Henry coefficient root locus gain controller gain electrode constant mass transfer coefficient oxygen transfer rate oxygen uptake rate smallest process pole larger process pole "gas" pole volume aerated part volume non-aerated part total reactor volume volume gas phase (holdup + head space) yield of biomass on oxygen zero of transfer function C1/Cg, i. zero controller

Greek letters s- ~ #max (Pc m 3 s- ~ (P0 m 3 s-

maximal specific growth rate circulation flow gas flow

Subscripts 1 2 c 9

in aerated part in non-aerated part controller gas phase

1 Introduction In a e r o b i c bioprocesses it is often necessary to c o n t r o l the dissolved o x y g e n (DO) c o n c e n t r a t i o n . In industrial f e r m e n t a tions, large-scale reactors are used for reasons of e c o n o m y . O n a large scale, h o w e v e r , gradients in the DO c o n c e n t r a t i o n m a y o c c u r due to imperfect mixing~ As a result, it is n o t possible to find a sensor l o c a t i o n representative which can serve the w h o l e reactor. This m a k e s c o n t r o l on a large scale m o r e difficult t h a n c o n t r o l in small scale applications. C o n t r o l of b i o r e a c t o r s is further c o m p l i c a t e d by the lack of sensors, e.g., for the active biomass. U n m e a s u r a b l e q u a n tities m u s t be r e c o n s t r u c t e d f r o m o t h e r m e a s u r e m e n t s via mass balances a n d / o r state estimation. T h e fact that process p a r a m e t e r s m a y v a r y d u r i n g the course of the process is a n o t h e r aspect of b i o t e c h n o l o g i c a l processes. Because the process p a r a m e t e r s m u s t be k n o w n for c o n t r o l l e r design, and possibly for state e s t i m a t i o n and on-line o p t i m i z a t i o n , it is desirable to estimate these p a r a m e t e r s on-line. A m e c h a n i s t i c m a t h e m a t i c a l m o d e l m i g h t be used to solve the a b o v e m e n t i o n e d problems. It has been s h o w n [1] t h a t gradients can be described by a t w o - c o m p a r t m e n t model. In this model, the a q u e o u s phase of the reactor is c o n c e p t u a l l y divided into an aerated a n d a n o n - a e r a t e d part, b o t h well mixed. T h e aerated p a r t represents the section in the r e a c t o r n e a r the stirrer, where the sensor is located. T h e c o m p a r t m e n t s e x c h a n g e o x y g e n via a circulation flow (q~c, Fig. 1).

226

Bioprocess Engineering 5 (1990)

l

Cg

Cg

Cx

C2

C•

OTR

Fig. 1. Two-compartment model for gradients on a large scale

The reactor is aerated by a constant gas flow ((Po), with an oxygen concentration Co, i,. The exhaust gas concentration C o and the DO concentration in the aerated part, C1, can be measured. The model has been used for the control of the DO concentration in a batch gluconic acid fermentation process [2]. The DO concentration in the non-aerated compartment, C z , and the biomass concentration, Cx, were estimated with an observer using the model, because they cannot be measured. The volume of the aerated part V1, the mass transfer coefficient k, a and the circulation flow (Pc may vary during the course of the process, so these parameters were estimated in a separate parameter estimator. The estimated C z was controlled on set point, but the estimate did not always converge on the concentration in the process. As a result, the DO concentration was not always controlled correctly. The subject of this paper is the applicability of the twocompartment model for control, state and parameter estimation. Because proper functioning of the state estimation is very important, attention is paid to the state estimation, through studying the interaction between state and parameter estimation. Controller and observer design have been reconsidered.

2 Gluconic acid process and model Growth of Gluconobacter oxydans and its production of gluconic acid proceed well under the following environmental conditions [3]: pH=5; temperature = 2 4 - 26 ~ glucose concentration >>5 #M. The oxygen concentration in the medium influences the growth and product formation rate. When the oxygen concentration is greater than or equal to 50 #M and the growth rate is almost maximal, it is desired to keep the DO concentration at least at 50 #M. For this purpose a controller is used. If it is assumed that the controller will function correctly, then the growth rate is constant and maximal. The other model hypotheses and assumptions for the formulation of the mathematical model are:

the biomass concentration Cx is homogeneous, because the growth rate is much smaller than the mixing rate; C 1 can be measured with an electrode; the response of the electrode for C 1 is a first-order response; the gas concentration C o can be measured; the gas phase and both compartments are assumed to be completely stirred; the growth is maximal, maintenance can be neglected, yield is constant; the product formation is maximal under the afore-mentioned conditions and is not included in the model; the oxygen consumption rate is proportional to the growth rate. Combination of mass balances and kinetic equations results in the following model equations:

dx=#max Cx u

(1) ~+kxa

C2 -- -- #max C

Yo.

-C1

"~

(C2--C1)

(Pc ~+V2 (C1-C2)

d~ = Vl k 1

+ ~(P. ( C o , i,-Co)

d e =ke(C 1-Ce)

(2)

(3) (4)

(5)

3 Controller setup and structure The DO concentration must not be lower than 50 #M. During the course of the batch, the biomass concentration increases and so does the oxygen uptake rate (OUR). During the first fermentation period, from t = 0 - t 1, air is supplied to the fermentor. The increasing OUR causes the DO concentration to decrease (Fig. 2). At t = t a , the concentration in the nonaerated part is 50 #M. Without further action, the concentration would quickly decrease further and the fermentation would have to be stopped. However, because at this moment the productivity of the fermentor is high, it is desired to continue the production. For this reason, from t = tl, the OTR is increased by mixing pure oxygen with the inlet gas to keep the concentration C 2 on set point. Now it is possible to continue until t = te, when the oxygen fraction in the inlet gas is one (corresponding to a concentration of 40 mM). At that moment, the fermentation is stopped. From t = t~, the concentration C 2 is constant, so dC2/ dt = 0. From Eq. (3) then follows:

#max Cx(t). (6) Ct(t)=C2 set+ V2 OUR(t), with OUR(t)= YT]o~ '

(Pc

From the time derivative dC1/dt and Eq. (2), the variation of the gas concentration Co can be computed. Continuing in this way, the trajectory for Co,~,, required to keep C2

Honderd et al.: Control in bioreactors showing gradients '

~

227 For the state estimation an observer is used. This is a copy (model) of the process, in which the process course is computed. The estimation error is fed back to the observer input to correct the estimates. For the parameter estimation a direct grid search is used [4]. This is an iterative method, which does not use gradient information. The method is only to be used when the number of parameters is small, say 2 or 3, because otherwise the method becomes too slow.

20

0

0

I

l

I

I

[

t\

tl

te

0 0.5 0 0.25

0

Fig. 2. Course of the concentrations during a batch, t = t~ : start of the control

4 S y s t e m analysis

The computation of the PI controller settings has been based upon the transfer functions C 1(s)/Co, i, (s) and C 2 (s)/Cg, i, (s). These transfer functions are derived here. Unless stated otherwise, the parameter values shown in the symbol list are used in this paper. In Fig. 4 a flow diagram of the model is given. The OUR has been omitted from this diagram, because it is not important in the computation of the transfer functions. The electrode response has also been omitted. For the determinant of the system can be written:

Cg,in

A=I

Vl k l a )

[\Vo

q~

VoH ]

q~c

~+

s-X

+kla

Fig. 3. Controller structure

-V~IV9 rj 1 ~g/Vg 1 / ~ a C9,in ~ 1

1

C9

1

-1/V~ ~

(\Vg

+

VoH J ~

(~

V2]

} k,a

1/V2 1/s ~

C~

Vlkla'~'(q~

C2

-1

Because the term V1 k~ a/(Vo H) is negligible compared to % / V o, the determinant can be simplified to Eq. (8):

Fig. 4. Flow diagram for computation of the transfer functions

A = ( s + q~_go)(s+pl)(s+p j s-3,

(8)

on set point, can be calculated:

Co, i, (t) = H C2, set + aco,i, 0 UR (t),

with Pl and P2 given in Eq. (9): (7)

the coefficient aco,~. is being a function of the process parameters. In Fig. 3 the controller structure is shown. The computed trajectory for C~ is used as the set point for a control loop, together with a feedback control loop for C 2. Because the OUR increases exponentially, C1 and C0,~, have to increase exponentially, too. If only the two PI loops are used, compensation of the exponential disturbation is not possible. For this reason, feedforward is used: the controller output Co, ~, is calculated with Eq. (7) from the estimated parameters and 0 UR. On this open-loop control the two feedback loops are superimposed. The processes in- and output are fed to a state and a parameter estimator. The state estimator uses the estimated parameters and vice versa. The estimated parameters are used to calculate the controller settings in the PI control loops, the feedforward control Co, ~, and the set point for C~.

)

+_

+~+kla

-4*k I

(9)

The poles, computed with Eqs. (8) and (9), approximate the exact values within 1 percent. The exact values have been computed with TRIP, a TRansformation and Identification Package [5]. The transfer functions now directly follow from the flow diagram. Multiplication of the numerator and denominator by s 3 gives Eqs. (10) and (11):

c~ (s) kl a . (G/(H"Vo) (s + qor G,i,

(s+p,)(s+p2)(s+%/vp

C=(s) k t a. q~c" qoo/(H" V2" Vo) co,,. = (~ + p~) (s + p J (~+ ~o/v~)"

(1o)

di)

228

Bioprocess Engineering 5 (/990)

~0.1

The location of the zero and the poles of the transfer function C1/Co,~, are shown in Fig. 5. The transfer function C2/Cg,~, has the same poles, but no zero.

5 Controller design

]

X

-0.2 P2

5.1 Controller settings

For the control in the two feedback loops, PI controllers have been chosen. The controller settings are based upon the (SISO) transfer functions C 1/Co, i, and C 2/Co, in" Because the process parameters may vary during the course of the process, formulas for the settings as a function of the parameters have been derived. For the controllers can be written: G~(s) = kp s + k 1 _ k~(s + z~) ,

s

X9

i

-0.1

X

pgz Pl 0

Fig. 5. Poles and zero for the transfer function C1/Co, I,; pl = -0.0074, z= -0.025, po= -0.033, P2 -0.168 =

sl

0.1

(12)

s

with z~= p t. With the controller zero, the most important pole, Pl is compensated. The resulting root loci for both control loops are shown in Fig. 6. In both cases, the controller gain is chosen such that the dominant poles of the system have a relative damping of 0.7. For the loop on C 1, these poles are located approximately

I -0.2

)" >-

I~ .............. "'....~0"1

on: sl=re(sl)(l+_j) ,

with

re(s1)= p z + p ~ 2

(13)

and for the control loop on C2: sz=re(sz)(l_+j),

with

re(s2) =po

(14)

Now the root locus gains can be computed. For the loop on C 1 the root locus gain is approximated by Eq. (15): k' 1,c=o.7 =

Isl-OI [Sl-pol [Sl-p2[

Isl-zl

~ Ist-011sl-p2l,

(15)

and for the loop on C2 by Eq. (16):

k~,(=o. 7 = [$2--0 [ [Sz--pe [ [S2--P2 ] ~ Is2--O[ 2 [$2--P2 [ . (16) The controller gains kcl and kc2 can be calculated from this root locus gain and the root locus gain of the transfer functions: kcl

k'1,~=O.7k,~ ' with '

kc2 kzd=~

'

with

9I

X

b-0.2

a a~ 'q%. k ] = kH

(17)

k 1 a (Po (P~ k ~ - H Vg V2

(18)

The controller has been implemented in a digital computer, with a sampling time of 5 s. The time constant of the fastest control loop is approximately 30 s. Discretization effects can be neglected if the sample time is less than or equal to 5 s.

,, / ]0

Fig. 6a and b. Root loci of the PI-control loops when Pt is compensated. The controller zero and Pl have been omitted

5.2 Observer

The observer is implemented in a digital computer, thus a discrete observer must be used. First, the process output is predicted with the controller output and the model of the system, by integration of the model Eq. (1-5): 2 ( k + 1 I - ) = {b(k) ~(k[ + ) + F ( k ) u(k)

prediction

(19)

Then the prediction is corrected by feedback of the prediction error: 2(k+ 1 ] + ) = 2 ( k + 11 - ) + K [ y ( k + l ) - C

2(k+11 - ) ] correction

(20)

Each of the model equations is corrected by the estimation error of C o and Ce, thus the matrix K is a 5 x 2 matrix. For the choice of the matrix K, the steady state of the Kalman matrix has been used [6]. Assuming that the variance of the system noise of the state equations is 10-5, and the measurement noise for C o and Ce is 10 -3 and 10 -5, respectively, the matrix shown in Table 1 was computed, in which the factors for C x have been multiplied by a factor 2.

Honderd et al. : Control in bioreactors showing gradients

229

Table 1. Observer matrix K

Feedback on:

Equation for:

Co - Cg

Ce__~

Cx

C I

C2

C0

Ce

- 0.006 --0.5

0.0047 0.4

0.0054 0.46

0.071 0.03

0.0003 0.62

manner, enabling study of the interaction of the control system components. First, the controller was tested for the case that all the states and parameters are known. Then the combination of controller and observer was tested and finally the complete control system. The process was simulated by the model Eqs. (1 5). 6.1 States and parameters known

Table 2. Performance of different controller configurations

Measure

Feedforward Feedback C 1 Feedback C 2

Feedback C, Feedback C 2

Feedback C 2

IAE Error at t~ [%]

0.02 0.00

13.12 3.50

71.25 19.02

0.25

v______

0 0.005

0.1

Feedforward control with feedback on C~ and C 2 , has been compared with control without feedforward, for the cases that feedback on C z is used, and feedback on C 1 and C z . In Table 2 the results are summarized, with the relative error in C2 at the end of the batch and the IAE criterion [8] as performance measure. When feedforward is applied, C2 can be controlled very well; the error on t e reduces to zero (within the numerical precision). This means an improvement compared to the old controller, where only feedback was applied. For applying feedforward, the parameters and the O U R must be well known. If deviations occur, it can be expected that (exponentially increasing) deviations will result.

/ ~c

.--------

r-,

~c

/ 0

I

0

0 0.001

10000 Time (s)

201)00

Fig. 7. Controller, observer and parameters estimation; (Pcestimated by pseudo steady state assumption

5.3 Parameter estimation

The model parameters V1 and kl a are estimated as independent fit parameters. The model error, defined by (21), is minimized: k

J=

Z

(Ce-Ce) z,

with:

j = time window.

(21)

i=k-j

Every 5 rain, the model error is computed over a time window of 10 rain with the estimated Cx, C 2 and C~ and the measured C o and C e as starting values. The circulation flow q~c may vary as well. The sensitivity of the criterion for q~c is too low, however, to use this parameter as an independent fit parameter. Instead, it has been suggested to compute (Pc from a pseudo steady state assumption for C2. If d C z / d t = 0, from Eq. (3, 22) can be derived: ~'2 ~o2 ~c = (C, - C2)"

(22)

This relation, of course, only holds when C2 is constant. 6 Simulation results

For the simulations, the program M U S I C [7] has been used. The control system has been implemented in a stepwise

6.2 States estimated, parameters known

The combination of controller and observer has been tested with exactly known parameter values. For a (worst case) begin estimate of x(o)= 0, the estimates have converged to the process values after approximately one-third of the batch period, all within I per mill. At the end of the batch, the estimated C 2 equals the process C a (within numerical accuracy), so C2 is controlled adequately. It is concluded that the combination of controller and observer functions well if the parameters are known. In results reported earlier [2], an offset in the estimate of C2 was observed. 6.3 States and parameters estimated

Finally, the controller with observer and parameter estimation was tested. In Fig. 7, the result is shown when (Pc, 1/1 and k 1 a are estimated, with correct parameter and state estimates at t = 0. Of the state variables, only C2 is shown. The other state estimates (Cx, C o, C1 and Ce) converged to the process states and have been omitted. The estimate of ~oc tends away from the real (Pc. In the estimates of k~ a and I/1 a coupling is observed, as reported earlier. As a result, the estimate of C 2 diverges. For the estimation of %, application of a pseudo steady state assumption had been suggested. Because at the beginning of the batch period C z is not constant, the assumption does not hold, and ~oc is not estimated accurately. This might be solved by taking the time derivative of C2 into account. There is, however, a more serious problem in the estimation of (Pc:for proper estimation of cpc, C 2 from the process would be needed. Only an estimate of C 2 is available, which has been computed with an a priori estimate of (Pc. In effect, the

230

Bioprocess Engineering 5 (1990)

a priori estimated qoc is corrected by the change in the estimate of V2. It is concluded that (Pc cannot be estimated by the pseudo steady state assumption. In the estimation of k 1 a and V1 a coupling is observed. A deviation in the estimate of k 1 a is accompanied with a deviation of the estimate for 1/1 in the opposite direction. This is also the case when (Pc is fixed and correct. These and earlier results give the impression that not 1/1 and k 1 a, but their product is estimated correctly. F r o m further simulations it appeared that this is not a rule: the product is often estimated incorrectly. In order to further elucidate the interaction of parameter and state estimation a series of simulations revealed the parameter values with which the state estimation is correct. When the value of the product V1. k 1 a and the value of the quotient V2/CPcwere both correct, the state estimation of C 2 was correct (within 3 per mill. at t~), even when the parameter deviation of k l a was a factor 2 ( k l a = 0 . 1 , V1 =1.25 9 10 - 3 , ~oc=2.8125 910-~). This means that coupling of the parameters V1 and k~ a will not result in a significant error in the estimate of C 2, if the values for V1 9 k I a and V2/q~~ can be estimated accurately. If one of the parameters is known and correct, the other two parameters can be directly calculated from these two quantities. Thus, if one of the parameters would be known, the other two should be estimated correctly. F r o m simulations with ~oc fixed and correct, however, this possibility seems unlikely, as it was observed that in that case V1 and k~ a were not (as a rule) estimated properly, resulting in divergence of the estimation of C2. This means that, even if q)c would be known a priori, the coupling of Vt an kl a would cause an incorrect estimation of C2. The results give the impression that measurement of C~ does not yield enough information to estimate the parameters and C2. Finally, simulations have been done with an error in one or more of the state estimates at t = 0. F r o m these simulations it became clear that the product of V1 and kl a is not always estimated correctly. Furthermore, it was observed that errors in the state estimates at t = 0 result in errors in the parameter estimates, and that the state and parameter estimation interact. The reason for this may be the difference in dynamic behavior of process and observer if the estimation error 5 0 . The parameter estimator fits the observer output to the process output. The observer, with feedback of the prediction error, has dynamic behaviour a different form than the process (if the error r 0). This effect might be circumvented by applying recursive techniques like extended Kalman filtering.

process states and parameters are known, the DO concentration can be controlled better by application of feedforward in combination with two feedback PI control loops than without feedforward. When the process parameters are known, the controller in combination with an observer also given very good resuits. N o offset in the estimate of C2, reported earlier, was observed. Because some of the model parameters may vary, it is desirable to estimate these parameters on-line. In the estimation of the parameters V1 and k1 a, however, a coupling exists and these parameters are not estimated correctly. A steady state assumption, suggested for the computation of ~0c from the other state and parameter estimates cannot be applied. Finally, interaction between observer and parameter estimator was observed. These effects together cause the state estimate of C a to diverge from the process value. The other state estimates (C x, C o, C~ and Ce), however, converge to the process states. These results suggest that measurement of C 1 along does not yield enough information for estimation of both parameters and the oxygen concentration in the non-aerated compartment and more a priori knowledge of the parameters is needed than suggested earlier.

7 Conclusions

K. Ch. A. M. Luyben Dept. of Biochemical Engineering Technical University Delft POB 5045 Julianalaan 67, NL-2628 BC Delft The Netherlands

For control of the DO concentration in reactors showing gradients, use of a two-compartment model for control, state and parameter estimation has been suggested. When all the

References 1. Oosterhuis, N. M. G.: Scale up of bioreactors, a scale down approach. Ph.D. Thesis Delft University of Technology, Huisdrukkerij Suiker Unie (1984) 2. van Breugel, J.: Optimization and control of biotechn, processes using mechanistic mathematical models. Ph.D. Thesis, Delft University of Technology (1986) 3. Olijve, W.: Glucose metabolism in Gluconobacter oxydans. Ph.D. Thesis, Delft University of Technology (1978) 4. Beech, G.: FORTRAN-IV in Chemistry, an introduction to computer assisted methods. N.Y.: J. Wiley & Sons 1975 5. van den Bosch, P. P. J.: TRIP version 4.0. Delft University of Technolgy 1985 6. Gelb, A.: Applied Optimal Estimation. Cambridge: MIT Press 1974 7. Soeterboek, A. R. M.: MUSIC versie 2.4. Delft University of Technology 1987 8. Dorf, R.: Modern Control Systems, 4 tla ed. Reading, Mass.: Addison-Wesley Publ. Comp. 1987 Received August 24, 1989 G. Honderd (corresponding author) S. R. Weijers Laboratory for Control Engineering Faculty of Electrical Engineering Technical University Delft POB 5031 2600 GA Delft