Multivariable Nonlinear Control of Biomass and Metabolite

Ind. Eng. Chem. Res. 2006, 45, 8985-8997

8985

Multivariable Nonlinear Control of Biomass and Metabolite Concentrations in a High-Cell-Density Perfusion Bioreactor Jean-Se´ bastien Descheˆ nes* and Andre´ Desbiens Laboratoire d’ObserVation et d’Optimisation des Proce´ de´ s (LOOP), De´ partement de ge´ nie e´ lectrique et de ge´ nie informatique, UniVersite´ LaVal, Que´ bec, Canada, G1K 7P4

Michel Perrier Department of Chemical Engineering, EÄ cole Polytechnique, Montre´ al, Canada, H3T 1J4

Amine Kamen NRC Biotechnology Research Institute (Bioprocess Sector), Montre´ al, Canada, H4P 2R2

This paper presents the development of a multivariable nonlinear adaptive controller for perfusion bioreactors, and its simulated behavior on a model that has been identified from experimental data. A contribution to the bioprocess model is also proposed, which is supported by experimental observations. The proposed control strategy is a multivariable approach to regulate the biomass and substrate concentrations using the fresh medium addition and direct bleeding streams as the manipulated variables. Level control would be ensured by a proportional integral (PI) control loop, using either the perfusion flow (draining flow that retains the cells in the reactor) or a nutrient-free phosphate-buffered saline (PBS) solution flow added to the reactor. The flow that is used for level control determines the operation mode of the reactor, being perfusion or chemostat. This paper presents the controller design, switching considerations between the operation modes, and parameter tuning guidelines. The controller tunings are essentially obtained by pole placement, based on linearization of the closed-loop dynamics. Simulation results prove the technique to be rather efficient, while the transitions between the operation modes are smooth and without risks. 1. Introduction In bioreactors, most implemented control strategies are monovariable, using either the dilution rate1-8 or feeding substrate concentration1-3 as the manipulated variable. However, there are reports of multivariable strategies that use both.1-3,9,10 Although monovariable strategies may only act to control a single variable, complex bioprocesses such as animal cell culturing would ideally require a far more controlled environment, because so many factors influence their yields. More global views on the simultaneous control of biomass and metabolite concentrations in a bioreactor were addressed,11,12 focusing on steady-state achievability12 or presenting a nonlinear control algorithm for a bioreactor operating in perfusion mode.11 However, both papers remarked that, to achieve arbitrary values for the variables, the control strategy should allow switching between the perfusion and chemostat modes. Reports in the literature about backstepping control algorithms for bioreactors are fairly limited: the only one found was written by Dochain and Perrier,13 who limited their study to the monovariable control of the biomass concentration in a continuous reactor through the use of the feeding substrate concentration. Gopaluni et al.14 presented a closely related report that focused on the multivariable control of a chemical continuously stirred tank reactor (CSTR). Control algorithms for bioprocesses usually include nonlinear feedback control1,4,5,7,9 and nonlinear model predictive control,6,10 and some studies were even made to apply linear controllers to bioreactors.2,3 The aim of this paper is to detail an algorithm for the multivariable nonlinear adaptive control of biomass and me* To whom correspondence should be addressed. Tel.: 418 656 2131, ext 12645. Fax: 418 656 3159. E-mail address: [email protected].

tabolite concentrations in a bioreactor. Basic backstepping design rules11 are used to build the controller. Because backstepping algorithms are typically sensitive to measurement noise, the objective will be to minimize this sensitivity. To achieve arbitrary steady-state values for both output variables simultaneously, the reactor will be allowed to switch between perfusion and chemostat modes as necessary.12 2. Bioreactor Model Consider the general bioreactor setting shown in Figure 1. Fresh media is added to the reactor by way of an input flow (Fin), which can be diluted by an addition flow (Fw) of a phosphate buffered saline (PBS) solution, to maintain adequate osmolality for the cell culture. The reactor is directly drained by a bleed flow (F1), while the perfusion flow (F2) also drains the reactor, although it is assumed to retain cells completely. The model considers the living biomass (XV), the dead biomass (Xd), and the substrate (S) concentrations. The biomass growth rate is denoted as µ, whereas k1 is the substrate consumption rate coefficient, kd is the death rate coefficient, and V is the reactor volume. A model can be derived from basic mass conservation equations, neglecting volume variations:

S˙ ) -k1µXV +

( ) ( ) () ()

Fin F1 + F2 Sin S V V

X˙ V ) (µ - kd)XV X˙ d ) kdXV -

10.1021/ie060582e CCC: $33.50 © 2006 American Chemical Society Published on Web 11/14/2006

F1 X V V

F1 X V d

(1a) (1b) (1c)

8986

Ind. Eng. Chem. Res., Vol. 45, No. 26, 2006

while the chemostat model (eq 3) becomes

Figure 1. Bioreactor setting.

The use of Fw and F2 were shown to be complementary.12 Indeed, to ensure control of the living biomass and substrate concentrations, as well as the reactor level in perfusion mode, F2 is manipulated by a proportional integral (PI) controller to maintain the level constant, whereas the concentrations XV and S are regulated by manipulating the flows Fin and F1. Because process dynamics are slow (order of days), and the reactor volume is ∼3 L, the constant volume operation hypothesis is easily justified. In this mode, because constant volume operation implies that the sum of F1 and F2 is equal to Fin, the model can be re-expressed as

Fin S˙ ) -k1µXV + (Sin - S) V X˙ V ) (µ - kd)XV X˙ d ) kdXV -

(2a)

F1 X V V

(2b)

F1 X V d

(2c)

In chemostat operation, Fw is manipulated instead of F2 (which is null) to achieve control of the reactor level. Because the amount of glucose brought to the reactor strictly originates from the fresh media addition flow (the glucose concentration in flow Fw is null), the following model is obtained:

S˙ ) -k1µXV + ) -k1µXV +

( ) ( ) ( ) ()

()

Fin Fw F1 Sin + ×0S V V V

Fin F1 S S V in V

X˙ V ) (µ - kd)XV X˙ d ) kdXV -

F1 X V V

()

F1 X V d

x˘ 1 ) θ1x1 - x1u1

(5a)

y 1 ) x1

(5b)

x˘ 2 ) θ2x1 + Sinu2 - x2u1

(5c)

y 2 ) x2

(5d)

Although the dead cell concentration is seemingly ignored in this control strategy, this does not mean that it should be completely disregarded. The algorithm presented here is, without any loss of generality, a low-level layer for eventual online optimization strategies. These higher-level algorithms should indeed make use of viability, among other offline measurements, for constraints handling and overall system diagnosis. Interestingly, the two models can be expressed as follows:

x˘ 1 ) f(θ1,x1) + g(x1)u1

(6a)

x˘ 2 ) f(θ1,θ2,x1,x2) + g(x1,x2)u1 + h(x1,x2)u2

(6b)

The process parameters θ1 and θ2 are assumed constant (or slowvarying),11 which results in x˘ 1 being independent of x2 and u2. Therefore, it is possible to design the multivariable controller sequentially, i.e., by first designing a control law to regulate y1 through u1, and then a second control law for y2 using u2. It might be pertinent to note that, at this point, the problem formulation remains general, because no specific structures are claimed for µ, k1, and kd. The controller design will remain based on these generalities. Because the model allows it, the backstepping approach,17 to which an integral action18 is added, will be used. The motivations for using this algorithm are that its design is based on stability considerations (Lyapunov stability), and its adaptive features provide mostly constant closed-loop performances. 3. Adaptive Control Strategy

(3a)

3.1. Controller Design. The first controller (for y1) is designed using the following Lyapunov function:

V1 )

(3b) (3c)

In practice, cell viability is the most difficult variable to estimate accurately. However, there are instruments available for the online measurement of XV15 and S.16 Thus, focusing on XV and S as the controlled variables, the previous dynamic models can be rearranged by redefining XV as y1, S as y2, F1/V as u1, Fin/V as u2, and (µ - kd) and -k1µ as parameters θ1 and θ2, respectively. The perfusion model (eq 2) becomes

x˘ 1 ) θ1x1 - x1u1

(4a)

y 1 ) x1

(4b)

x˘ 2 ) θ2x1 + (Sin - x2)u2

(4c)

y 2 ) x2

(4d)

() ( )

1 2 1 z1 + θ˜ 2 2 2γ1 1

(7)

where

z1 ) kp11 + ki1

∫ 1 dt

1 ) yr1 - y1 θ˜ 1 ) θ1 - θˆ 1

(8)

The reference trajectory (i.e., the filtered set point) is yr1; thus, 1 is the closed-loop error. Proportional and integral actions are included in the definition of the augmented error z1. The parameters kp1 and ki1 are the corresponding tuning knobs. The estimated value of θ1 is θˆ 1; thus, θ˜ 1 is the estimation error. The positive adaptation gain is γ1. Assuming that θ˙ 1 = 0 (i.e., that, in practice, the time variations of θ1 are slow), the Lyapunov function’s derivative is given by

Ind. Eng. Chem. Res., Vol. 45, No. 26, 2006 8987

V˙ 1 ) z1z˘ 1 +

θ˜ 1 θ˜ γ1 1

V2 ) z2z˘ 2 +

( )

) z1[kp1(y˘ r1 - θ1x1 + x1u1) + ki11] + θ˜ 1 -

θ˙ 1 γ1

(

( )

) z2{kp2[y˘ r2 - θ2x1 - (Sin - x2)u2] + ki22} + θ˜ 2 -

( ) )

) z1{kp1[y˘ r1 - (θ˜ 1 + θˆ 1)x1 + x1u1] + ki11} + θ˜ 1 ) z1[kp1(y˘ r1 - θˆ 1x1 + x1u1) + ki11] + θ˜ 1 -

θ˜ 2 θ˜ γ2 2

θ˙ 1 γ1

θ˙ 1 - kp1z1x1 γ1 (9)

To ensure stability, V˙ 1 is forced to be negative:

) z2{kp2[y˘ r2 - (θ˜ 2 + θˆ 2)x1 - (Sin - x2)u2] + ki22} +

( )

θ˜ 2 -

V˙ 1 ) -c1z1

(

)

(17)

(10)

where c1 is a positive design parameter. Therefore, the combination of eqs 9 and 10 leads to

kp1(y˘ r1 - θˆ 1x1) + kp1x1u1 + ki1 ) -c1z1

Therefore, the control and adaptation laws are as follows:

u2 )

(11)

(

)

c2z2 ki22 1 + + y˘ r2 - θˆ 2y1 Sin - y2 kp2 kp2

(18a)

θˆ 2 ) -γ2kp2z2y1

and

θˆ 1 ) -γ1kp1z1x1

(12)

From eq 11, the control law is deduced:

u1 )

(

)

c1z1 ki11 1 - y˘ r1 + θˆ 1 x1 kp1 kp1

(13)

whereas eq 12 is the adaptation law. However, without an observer for x1, its best estimate available is xˆ 1 ) y1. Therefore, the realizable control and adaptation laws are

(

θˆ 2 γ2

) z2{kp2[y˘ r2 - θˆ 2x1 - (Sin - x2)u2] + ki22} + θˆ 2 θ˜ 2 - - kp2z2x1 γ2 ) -c2z22

2

θˆ 2 γ2

)

(18b)

where xˆ 1 ) y1 and xˆ ) y2. Again, the design parameter c2 and the adaptation gain γ2 must be positive. For the second operation mode, where Fw controls the reactor level, the same Lyapunov function V2 (eq 15) can be kept to design the controller. Its derivative is, instead, given by

V˙ 2 ) z2z˘ 2 +

θ˜ 2 θ˜ γ2 2

( )

) z2[kp2(y˘ r2 - θ2x1 + x2u1 - Sinu2) + ki22] + θ˜ 2 -

θˆ 2 γ2

c1z1 ki11 1 - y˘ r1 + θˆ 1 u1 ) y1 kp1 kp1

(14a)

) z2{kp2[y˘ r2 - (θ˜ 2 + θˆ 2)x1 + x2u1 - Sinu2] + ki22} +


(14b)

θ˜ 2 -

Because the dynamic model for y1 does not change, regardless of whether F2 or Fw is used to control the reactor volume, this same controller will be usable in both modes. For the second output variable, however, both scenarios must be anticipated. A first controller can be designed for the first model, i.e., while F2 is used for level control. Let the following Lyapunov function V2 be

1 1 2 θ˜ V2 ) z22 + 2 2γ2 2

(15)

(

) z2[kp2(y˘ r2 - θˆ 2x1 + x2u1 - Sinu2) + ki22] + θ˜ 2 -

∫ 2 dt

θˆ 2 γ2

θˆ 2 γ2

kp2z2x1 )

-c2z22

)

(19)

The resulting control and adaptation laws are

u2 )

where

( )

(

)

1 c2z2 ki22 + + y˘ r2 - θˆ 2y1 + y2u1 Sin kp2 kp2

(20a)

(16a)


2 ) yr2 - x2

(16b)

θ˜ 2 ) θ2 - θˆ 2

(16c)

where xˆ 1 ) y1, and the adaptation law is the same as that in eq 18. One can also notice the fairly slight differences between the control laws found in eqs 20 and 18: the error definition and z2 are the same, as are the reference trajectories. The fact that these terms remain the same for both operation modes

z2 ) kp22 + ki2

The Lyapunov function’s derivative is given by

(20b)

8988


allows their free use in either function that determines u2, thus simplifying the switching task. 3.2. Tunings and Closed-Loop Dynamics. Parameter tuning for adaptive backstepping controllers has always been awkward. No conventional techniques exist to calculate values, and the exact individual effects of most parameters are often vague: they are often chosen after trial-and-error simulations.13,14 In this paper, the closed-loop dynamics of the reactor under backstepping control will be approximated, and transfer functions will be derived to find guidelines for controller parameter tuning, which the intent to reduce sensitivity to measurement noise. To represent the measurement noise, the presence of output disturbances on the measured variables is considered. For the first output, y1, the dynamic model is given by

x˘ 1 ) θ1x1 - x1u1

(21a)

y1 ) x1 + p1

(21b)

This system always has a real pole at ki1/kp1, whereas the two other poles can be either real or complex. If s2 and s3 are real, they will be symmetrically apart from c1/2. Because these poles correspond to the disturbance rejection specification, it is best not to fix them too quick, because fast disturbance rejection also means output noise sensitivity. The most straightforward way to set the specifications is to force the relation s2 ) s3 ) c1/2. In summary, the following tuning technique is proposed: first, the desired specifications on the closed-loop dynamics are set. Two of the three poles will be set equal to c1/2: this gives the necessary value for c1. The third pole specification will give the ratio of ki1/kp1. For an arbitrary value of γ1, kp1 will be set according to the following rule, considering an initial steadystate value xj1:

c1

kp1 )

where p1 is the disturbance. Thus, the output error is

1 ) yr1 - y1 ) yr1 - x1 - p1

(22)

Substituting the adaptation and control laws from eq 14, it then follows that

˘ 1 ) y˘ r1 - x˘ 1 - p˘ 1

( )(

) y˘ r1 - θ1x1 - x1

With these tunings, the resulting estimation and regulation dynamics become independent of the value of γ1 (the demonstration is straightforward; eq 28 is used to replace kp1 in eqs 26 and 12). For the second output, there are two possible cases:

x˘ 2 ) θ2x1 + (Sin - x2)u2

)

1 c1z1 ki11 + + y˘ r1 + x1θˆ 1 - p˘ 1 y1 kp1 kp1 (23)

u2 )

)

(

(29)

)

Sin - x2 c2z2 ki22 + + y˘ r2 - θˆ 2y1 Sin - y2 kp2 kp2

and

Y1(s) ) Yr1(s) +

x˘ 2 ) θ2x1 + Sinu2 - x2u12

s3P1(s) s3 + [(ki1/kp1) + c1]s2 + [(c1kil/kp1) + γ1k2p1xj12]s + γ1kp1ki1xj12 (24) As expected, the transfer function between Yr1(s) and Y1(s) is unity. For disturbance rejection, the characteristic equation is

(

(

c2z2 ki22 1 + + y˘ r2 - θˆ 2y1 Sin - y2 kp2 kp2

∴ x˘ 2 ) θ2x1 +

Linearization of the closed-loop system (the details of which are presented in the Appendix) converges toward the following transfer function:

s3 +

(28)

2xj1xγ1

) (

)

ki1 c1ki1 + c1 s2 + + γ1k2p1xj12 s + γ1kp1ki1xj12 ) 0 kp1 kp1 (25)

u2 )

∴ x˘ 2 ) θ2x1 +

) ( )

γ1k2p1xj12

˘ 2 ) y˘ r2 - y˘ 2 ) y˘ r2 - x˘ 2 - p˘ 2 ) y˘ r2 - θ2x1 -

s + γ1kp1ki1xj1 (26) 2

ki1 kp1

s2 )

c1 x + 2

s3 )

c1 xc21 - 4γ1k2p1xj12 2 2

-

4γ1k2p1xj12 2

)

and

˘ 2 ) y˘ r2 - θ2x1 -

c21

(

(Sin - x2) c2z2 ki22 + + y˘ r2 kp2 (Sin - y2) kp2 θˆ 2y1 - p˘ 2 (31)

and leads to

s1 )

(30)

c2z2 ki22 + + y˘ r2 - θˆ 2y1 + y2u1 - x2u1 kp2 kp2

ki1 c1ki1 + c1 s 2 ) + (s + s1)(s + s2)(s + s3) ) s + kp1 kp1 3

)

The two following error dynamics are then determined:

Setting the three poles to s1, s2, and s3 forces the relation

(

(

1 c2z2 ki22 + + y˘ r2 - θˆ 2y1 + y2u1 Sin kp2 kp2

c2z2 ki22 - y˘ r2 + θˆ 2y1 - y2u1 + x2u1 kp2 kp2 p˘ 2

) -(θ˜ 2 + θˆ 2)x1 (27) ) -θ˜ 2x1 -

c2z2 ki22 + θˆ 2(x1 + p1) - (x2 + kp2 kp2 p2)u1 + x2u1 - p˘ 2

c2z2 ki22 + θˆ 2p1 - p2u1 - p˘ 2 kp2 kp2

(32)


The following transfer functions are determined after linearization of the closed-loop systems:

this leads to Fin ) F1, or, stated otherwise, u1 ) u2. From eq 1, the following relationship then holds at steady state:

Y2(s) ) Yr2(s) + -θˆ 2s ki2 c2ki2 s3 + + c2 s 2 + + γ2k2p2xj12 s + γ2kp2ki1xj12 kp2 kp2 P1(s) + ... 2

( ( ‚‚‚ +

) (

)

(

θˆ 2xj1

s3 -

)

)

s2 (Sin - xj2) P1(s) ki2 c2ki2 3 2 2 2 2 s + + c2 s + + γ2kp2xj1 s + γ2kp2ki1xj1 kp2 kp2 (33)

( (

) (

)

)

and

Y2(s) ) Yr2(s) + -θˆ 2s2 ki2 c2ki2 s3 + + c2 s 2 + + γ2k2p2xj12 s + γ2kp2ki1xj12 kp2 kp2 P1(s) + ...

( (

) (

)

)

‚‚‚ + (s3 + θˆ 1s2) P2(s) ki2 c2ki2 s3 + + c2 s 2 + + γ2k2p2xj12 s + γ2kp2ki2xj12 kp2 kp2 (34)

( (

) (

)

)

Interestingly, for both operation modes, the closed-loop poles are identical. Setting these poles in the same way as that done previously leads to

ki2 kp2

(35a)

s2 )

c2 xc22 - 4γ2k2p2xj12 + 2 2

(35b)

s3 )

c2 xc22 - 4γ2k2p2xj12 2 2

(35c)

s1 )

Thus, the preceding tuning method can be used for both operation modes. 3.3. Switching. The study of the continuity in the behavior of the system at the switching point between perfusion and chemostat modes is reliant on the transfer functions at that point. The transfer functions, relative to output y1 (eq 24), have their parameters remaining constant or varying continuously with the operating point xj1. As for the parameters of the transfer functions that are related to output y2, these also vary continuously, as long as the following holds at the switching point:

θ h1 )

-θ2x1 (Sin - xj2)

(36)

The equivalence can be shown starting from the model that has been presented as eq 1. At the switching point, both flows F2 and Fw are null. With the reactor volume being kept constant,

0)θ h 2xj1 + uj2Sin - uj1xj2

(37a)

0)θ h 1xj1 - uj1xj1

(37b)

Combining and manipulating these equations directly leads to eq 36. The controller that manipulates F2 and Fw also has a continuous behavior at the switching point. Indeed, because these flows are never used together, the negative output of the controller can be treated as a positive value of Fw, while its positive output is treated as F2. 4. Simulations 4.1. Model Justifications. Simulation parameters are chosen to match, as closely as possible, the known properties of the HEK 293-SF cells grown in NSFM13 medium, which is a mammalian cell line and medium composition used at the NRC Biotechnology Research Institute (BRI) for adenovirus production. To best represent the behavior of a high-cell density culture, some aspects were considered. First, the growth rate function was chosen as a Contois model,1,19 because cell concentration effects are expected to become more significant at higher densities. Second, as often confirmed by experimental data, the biomass concentration (total and/or living) has a tendency to overshoot at the beginning of runs: to the knowledge of the authors, the phenomenon itself was never reported in the literature (hence, not modeled), although it has been observed by many scientists at the BRI over the 15 past years, with different cell lines. The following equation relates the Contois model structure, where S and X are the substrate and biomass concentrations (respectively), and µmax and KC are constants:

µ(S,X) )

µmaxS KCX + S

(38)

Figure 2 shows a perfusion run over 46 days, operated under three preparation batches of the same media. Before the switch to batch 2, a peak is briefly observed. Viability at this moment is ∼80%, and, as Figure 3 shows, this somewhat low viability is not due to low supplies of glucose, nor is it due to high concentrations of lactate. Moreover, the culture was not under oxygen limitations, nor was ammonia in sufficient concentration to inhibit growth (results not shown). For this run, however, amino acids were not investigated. Figure 4 shows a second perfusion run operated entirely under the same batch of media preparation. Once again, the same overshoot phenomenon can be shortly observed on the living biomass concentration, just before the recirculation pump is activated. Figure 5 shows the glucose and lactate concentrations during the culture, once again to confirm that there is no limitation or inhibition from these factors. Here again, the culture was not under oxygen limitations, nor was ammonia present in sufficient concentration to hamper growth. The amino acids were also investigated, and none were shown to be within limitations ranges (results not shown). Other perfusion runs have shown this same overshoot phenomenon, not only with 293-SF cells, but with hybridomas and other cells as well (data not shown). However, no models were found in the literature to account for this phenomenon. Figure 6 shows a perfusion simulation with arbitrary parameter values, under a Contois growth model and a constant death rate parameter: the culture that is shown operates under constant

8990


Figure 2. Biomass concentrations in first perfusion run.

Figure 3. Viability and glucose/lactate concentrations in the first perfusion run.

bleed and feed rates. As expected, no overshoot is observed on the biomass concentration. Figure 7 shows a perfusion simulation under similar conditions, with a death rate that now varies in proportionality to the total cell concentration:

kd ) k′d(XV + Xd)

(39)

Such a modification of the model allows the overshoot to be observed, thus representing practical observations. Hopefully, it would also prove to better fit the experimental data. At the start, the two cultures from Figures 2-5 were intended to be chemostat runs. The setup was actually a chemostat, although, eventually, the cells started to coat the bleed pipes, thus retaining the cells more and more in the reactor. It took until some point during run 2 to actually notice the problem, which was then solved by starting the recirculation pump at a

relatively high flow rate to ensure the cells to be drawn out more adequately. From this moment, the reactor was 100% sure to function in chemostat mode. The data were still fitted to the model, whereas the pipe-jamming phenomenon was taken into account and was attributed to first-order dynamics. Figure 8 shows a good fit between the experimental and simulation data. The estimated retention rate during the culture is shown in Figure 9. Table 1 gives the model parameters used in simulations: µmax is the maximum growth rate, KC is the saturation constant in the Contois growth rate model, parameter k′d is the cell death proportionality constant, and k1 is the substrate consumption rate coefficient. Sin is the glucose concentration in the fresh media composition, and V is the reactor operation volume.


Figure 4. Biomass concentrations in the second perfusion run.

Figure 5. Viability and glucose/lactate concentrations in the second perfusion run.

Figure 6. Culture simulation with a constant death rate.

8992


Figure 7. Culture with the death rate proportional to the total cell concentration.

Figure 8. Fitted experimental data from perfusion run 2. Table 1. Process Parameters Used in Simulations parameter µmax KC k′d k1 Sin V

value 0.0504 h-1 1.87 mM/L/106 cells 0.000758 h-1/106 cells 5.2 30 mM/L 2.65 L

Table 2. Parameter Tunings for the Adaptive Controller

Figure 9. Estimated cell retention ratio during experiment 2.

The controller described in section 3.1 will be tested in simulation with the aforementioned model. The controller parameter values that have been used are presented in Table 2, being set for an operation concentration of 6 × 106 cells/mL.

a

parameter

value

Tref kp1 ki1 kp2 ki2 γ1 γ2 c1 c2

10 h 0.167a 0.0167a 0.167a 0.0167a 0.01 0.01 0.2 0.2

Calculated for a cell concentration of ∼6 × 106 cells/mL.

Consider two distinct steady states, for perfusion and chemostat operation, respectively. In both cases, the substrate concentration is regulated at 16 mM/L. Using these tunings, the


Figure 10. Comparison of output perturbation responses in perfusion operation mode.

Figure 11. Comparison of output perturbation responses in chemostat operation mode.

linearized closed-loop transfer functions for perfusion operation are

Y1(s)

)

Yr1(s) Y1(s)

)1

(40a)

s3 s3 + 0.3s2 + 0.03s + 0.001

(40b)

)

0.138s3 s + 0.3s + 0.03s + 0.001

(40c)

)

s3 + 0.0592s3 s3 + 0.3s2 + 0.03s + 0.001

(40d)

P1(s) Y2(s)

Yr2(s)

)

P1(s) Y2(s)

Y2(s)

P2(s)

3

2

Keeping the same tunings for chemostat operation (2 × 106 cells/mL) would lead to

Y1(s) Yr1(s) Y1(s)

Y2(s) Yr2(s)

)1

(41a)

)

s3 s3 + 0.3s2 + (19/900)s + (1/9000)

(41b)

)

0.211s2 s3 + 0.3s2 + (19/900)s + (1/9000)

(41c)

P1(s) Y2(s)

)

P1(s) Y2(s) P2(s)

)

s3 + 0.0389s2 19 s + (1/9000) s3 + 0.3s2 + 900

(41d)

Figure 10 shows the step responses of the linearized models, compared to those of the actual nonlinear system to output perturbations in perfusion operation. Figure 11 shows the responses in chemostat operation.

For both cases, the transfer functions obtained are good approximations of the actual system responses. Thus, the pole placement method can be considered to be reliable. Figure 12 shows the system under white measurement noises on both measured variables, with a variance of 0.01. The objective of this test is to compare the measurement noise sensitivity using the proposed tunings with those used by Descheˆnes et al.11 It is expected that the sensitivity to measurement noise can be reduced by adequately placing the regulation poles. Indeed, for a measurement noise of variance 0.01 on both measured variables, the variances observed on u1 and u2 were, respectively, 0.32 and 0.94,11 while they are now 0.016 and 0.10. Because this last paper showed the sensitivity to noise on the parameters as not being significant, this aspect was not included for this specific test (see Figure 12). For the simulation behavior to closely match reality, further tests shall include additive noise on the parameters θ1 and θ2. Because the controller is still rather noise-sensitive, both measurements will be filtered by a 2-h-time-constant first-order filter: this completes the simulation setup and control strategy tunings.11 A main concern about this strategy is what might occur at the switching point. Figures 13 and 14 show the response to a set-point change toward this state, where, at equilibrium, the feed and bleed flows are equal. Figure 13 shows a satisfying behavior of the system at that point. Figure 14 shows the level controller output (note that our simulator treats its positive value as Fw, and its negative value as a positive F2) and the outputs of the two controllers for the feed rate (u2). Because the latter are actually equal at this equilibrium state, there can be no risk associated with the switching between the two possible outputs. As seen, switching does not cause any stability problem, nor does it create competition between the two controllers. An important test to run on this algorithm is to operate far from the initial set point for which the tunings were actually

8994


Figure 12. System under measurement of white noise.

Figure 13. Simulation at a chemostat limit.

calculated. Figure 15 shows the simulation of a set-point change in biomass from 6 × 106 cells/mLto 2 × 106 cells/mL, followed by a set-point change in substrate concentration from 16 mM/L to 20 mM/L. It is observed that the control itself is not affected by operating away from the initial set-point, although parameter estimations show a severe lag. However, because the actual goal of the algorithm is the control of the measured variables, estimation of the parameters is not significantly important. 6. Summary and Conclusion This paper has presented the implementation of a multivariable adaptive backstepping control strategy on a simulated

bioreactor, with parameters actually fitted on experimental data. A contribution to bioreactor modeling was proposed, which explains a certain experimental cell culture behavior. Guidelines for controller parameter tuning were developed through linear approximation of the closed-loop dynamics. Results showed a reduced sensitivity to measurement noise, and a certain robustness of the control algorithm itself. Even though the parameter estimates were shown to lag severely the further one operated from the original set points, this did not affect the control performances. The algorithm is highly flexible over operation ranges for both biomass and substrate concentrations, and switching between operation modes proved smooth and without risks.


Figure 14. Control input for level control and outputs of controllers for u2.

Figure 15. Effect of operating far from the initial steady-state conditions.

7. Appendix

It then follows that

This section details the linearization developments for the closed-loop transfer functions on the first output (y1). To represent the measurement noise, the presence of an output disturbance is introduced. The following dynamic model represents the system:

y˘ r1 - y˘ 1 ) y˘ r1 - θ1x1 -

x˘ 1 ) θ1x1 - x1u1

(A1a)

y1 ) x1 + p1

(A1b)

for which the adaptation and control laws were given in eq 14. Thus, the output error is

1 ) yr1 - y1 ) yr1 - x1 - p1

(A2)

and its derivative is given by

˘ 1 ) y˘ r1 - x˘ 1 - p˘ 1

(

)

1 c1z1 ki11 ) y˘ r1 - θ1x1 - x1 + + yr1 + x1θˆ 1 - p˘ 1 y1 kp1 kp1

(A3)

x1

)

ki1 c1 z1 + 1 + y˘ r1 + k k (x1 + p1) p1 p1 x1θˆ 1 - p˘ 1

) y˘ r1 - (θ˜ 1 + θˆ 1)x1 -

) y˘ r1 - θ˜ 1x1 -

(

c1 x1z1 ki1 x11 kp1 x1 + p1 kp1 x1 + p1 x1y˘ r1 + x1θˆ 1 - p˘ 1 x1 + p1

ki1 x11 c1 x1z1 kp1 x1 + p1 kp1 x1 + p1 x1y˘ r1 - p˘ 1 x1 + p1

(A4)

At this step, a linear approximation must be made. Using a firstorder Taylor series around steady-state values (which are represented by a bar over the variable), one finds

8996


θ˜ 1x1 = θ˜ 1x1 + θ˜ 1(x1 - xj1) + x1(θ˜ 1 - θ˜ 1)

( ) ( )

(A5a)

xj1jz1 x1z1 x1z1 ∂ = + (x - xj1) + ... x1 + p1 xj1 + pj1 ∂x1 x1 + p1 1

( )

... +

[

[

]

-xj1jz1

(xj1 + pj1)

1

2

(p1 - pj1) +

[

]

1

xj1 (xj1 + pj1)

]

(z1 - jz1)

Under initial conditions such that the system is at steady state and no disturbances are present, although the model is exact, this leads to jz1 ) 0, pj1 ) 0, and θ h 1 ) 0:

θ˜ 1x1 = xj1θ˜ 1

(A6a)

x1z1 = z1 x1 + p1

(A6b)

The same assumptions lead to 1 ) 0 and y˘ r1 ) 0. Thus, the following nonlinear functions can be approximated in the same way:

x11 = 1 x1 + p1 x1y˘ r1 = y˘ r1 x1 + p1

(A7a) (A7b)

It then follows that

() () () ()

c1 ki1 z1 - y˘ r1 - p˘ 1 y˘ r1 - y1 ) y˘ r1 - xj1θ˜ 1 kp1 kp1 1 (A8) c1 ki1 ) -xj1θ˜ 1 z1 1 - p˘ 1 kp1 kp1 Deriving a second time allows the estimation dynamics to appear:

c1 ki1 z˘ 1 - ˘ 1 - p¨1 kp1 kp1

y¨r1 - y¨1 ) -xj1θ˜ 1 ) xj1θˆ 1 -

c1 ki1 z˘ 1 - ˘ 1 - p¨1 kp1 kp1

) xj1(-γ1kp1z1y1) -

c1 ki1 z˘ ˘ - p¨1 kp1 1 kp1 1

which is an expression that should be linearized:

= yj1z1

(A10)

The initial condition on y1 is

yj1 ) xj1 + pj1 ) xj1 Equation A9 now becomes

c1 ki1 z˘ 1 ˘ - p¨1 kp1 kp1 1

(A12)

[

s2Yr1(s) - s2Y1(s) ) -γkp1xj12 kp1(Yr1(s) -Y1(s)) +

]

ki1 (Y (s) - Y1(s)) + ... s r1 (A13) c1 ... - [kp1s(Yr1(s) - Y1(s)) + ki1(Yr1(s) - Y1(s))] kp1 ki1 s(Y (s) - Y1(s)) - s2P1(s) kp1 r1 Then, by regrouping,

()

( )

ki1 c1ki1 sYr1(s) + c1sYr1(s) + Y (s) + kp1 kp1 r1 ki1 Yr1(s) γ1k2p1xj12Yr1(s) + γ1kp1ki1xj12 ) s2Y1(s) + sYr1(s) + s kp1 c1ki1 Y (s) + γ1k2p1xj12Y1(s) + c1sY1(s) + kp1 1 Y1(s) γ1kp1ki1xj12 - s2P1(s) s

s2Yr1(s) +

( )

( )

( )

∴ Y1(s) ) Yr1(s) ) s3P1(s) s3 + [(ki1/kp1) + c1]s2 + [(c1ki1/kp1) + γ1k2p1xj12]s + γ1kp1ki1xj12 (A14) The same procedure is applied for Y2(s). However, the details will not be provided in this paper, because they are very similar to those just provided here. Acknowledgment The authors would like to acknowledge the “Fonds Que´be´cois de Recherche sur la Nature et les Technologies” (FQRNT) for funding, the Biotechnology Research Institute (BRI) of the National Research Council Canada (NRC) for participation in the project, and the technical staff at BRI for practical issues on animal cell culture. Literature Cited

(A9)

() ()

z1y1 = jz1yj1 + jz1(y1 - yj1) + yj1(z1 - jz1)

() ()

which can be converted to Laplace transforms:

x1z1 ∂ ∂ x1z1 (p1 - pj1) + (z - jz1) ∂p1 x1 + p1 ∂z1 x1 + p1 1 (A5b) -xj1jz1 xj1jz1 jz1 + + (x - xj1) + ... = xj1 + pj1 xj1 + pj1 (xj + pj )2 1

... +

y¨r1 - y¨1 ) -γkp1z1xj12 -

(A11)

(1) Bastin, G.; Dochain, D. On-line Estimation and AdaptiVe Control of Bioreactors; Elsevier: Amsterdam, 1990. (2) Zhao, Y.; Skogestad, S.; Comparisons, of various control configurations for continuous bioreactors. Ind. Eng. Chem. Res. 1997, 36, 697-705. (3) Wang, Z. Q.; Skogestad, S.; Zhao, Y. Exact linearization control of continuous bioreactors: a comparison of various control structures. In Proceedings of the Second IEEE Conference on Control Applications; Vancouver, British Columbia, Canada, 1993; pp 107-112. (4) Mailleret, L.; Bernard, O.; Steyer, J. P. Nonlinear adaptive control for bioreactors with unknown kinetics. Automatica 2004, 40, 1379-1385. (5) Szederkenyi, G.; Kristensen, N. R.; Hangos, K. M.; Jorgensen, S. B. Nonlinear analysis and control of a continuous fermentation process. Comput. Chem. Eng. 2002, 26, 659-670. (6) Parker, R. S.; Doyle, F. J. Optimal control of a continuous bioreactor using an empirical nonlinear model. Ind. Eng. Chem. Res. 2001, 40, 19391951. (7) Konstantinov, K. B.; Tsai, Y. S.; Moles, D.; Matanguihan, R. Control of long-term perfusion chinese hamster ovary cell culture by glucose auxostat. Biotechnol. Prog. 1996, 12, 100-109.

Ind. Eng. Chem. Res., Vol. 45, No. 26, 2006 8997 (8) Gostomski, P.; Muhlemann, M.; Lin, Y. H.; Mormino, R.; Bungay, H. Auxostats for continuous culture research. J. Biotechnol. 1994, 37, 167177. (9) Zhang, Y.; Zamamiri, A. M.; Henson, M. A.; Hjortso, M. A. Cell population models for bifurcation analysis and nonlinear control of continuous yeast bioreactors. J. Process Control 2002, 12, 721-734. (10) Zhu, G. Y.; Zamamiri, A.; Henson, M. A.; Hjortso, M. A. Model predictive control of continuous yeast bioreactors using cell population balance models. Chem. Eng. Sci. 2000, 55, 6155-6167. (11) Descheˆnes, J. S.; Desbiens, A.; Perrier, M.; Kamen, A. Adaptive backstepping control of biomass and metabolite concentrations in a perfusion bioreactor, Asia-Pacific J. Chem. Eng. 2006, in press. (12) Descheˆnes, J. S.; Desbiens, A.; Perrier, M.; Kamen, A. On Simultaneous Control of Biomass and Metabolite Concentration in Perfusion Bioreactors. In Proceedings of the DCDIS 4th International Conference on Engineering Applications and Computational Algorithms, Guelph, Ontario, Canada, 2005; pp 663-667. (13) Dochain, D.; Perrier, M. Adaptive backstepping nonlinear control of bioprocesses. In Proceedings of ADCHEM, Hong Kong, 2003; pp 7782. (14) Gopaluni, R. B.; Mizumoto, I.; Shah, S. L. A robust nonlinear adaptive backstepping controller for a CSTR. Ind. Eng. Chem. Res. 2003, 42, 4628-4644.

(15) Zeiser, A.; Elias, C. B.; Voyer, R.; Jardin, B.; Kamen, A. On-line monitoring of physiological parameters of insect cell cultures during the growth and infection process, Biotechnol. Prog. 2000, 16, 803-808. (16) Ozturk, S. S.; Thrift, J. C.; Blackie, J. D.; Naveh, D. Real-time monitoring and control of glucose and lactate concentrations in a mammalian cell perfusion reactor. Biotechnol. Bioeng. 1997, 53, 372-378. (17) Kstic´, M.; Kanellakopoulos, I.; Kokotovic´, P. Nonlinear and AdaptiVe Control Design, Wiley-Interscience: New York, 1995. (18) Aarset M. F.; Strand, J. P.; Fossen, T. L. Nonlinear vectorial observer backstepping with integral action and wave filtering for ships. In IFAC Control Applications in Marine Systems, Fukuoka, Japan, 1998; pp 77-82. (19) Contois, D. E. Kinetics of bacterial growth: relationship between population density and specific growth rate of continuous cultures. J. Gen. Microbiol. 1959, 21, 40-50.

ReceiVed for reView May 10, 2006 ReVised manuscript receiVed October 4, 2006 Accepted October 6, 2006 IE060582E

Multivariable Nonlinear Control of Biomass and Metabolite

Recommend Documents