Nonlinear Analysis of Substrate-Inhibited Continuous Cultures

Article pubs.acs.org/IECR

Nonlinear Analysis of Substrate-Inhibited Continuous Cultures Operated with Feedback Control on Dissolved Oxygen Giuseppe Olivieri,† Maria Elena Russo,‡ Pier Luca Maffettone,† Erasmo Mancusi,§,⊥ Antonio Marzocchella,*,† and Piero Salatino† †

Dipartimento di Ingegneria Chimica, dei Materiali e della Produzione Industriale, Università degli Studi di Napoli Federico II, P.le V. Tecchio n. 80, 80125 Napoli, Italy ‡ Istituto di Ricerche sulla Combustione C.N.R., P.le V. Tecchio n. 80, 80125 Napoli, Italy § Dipartimento di Ingegneria, Università degli Studi del Sannio, Piazza Roma, 21, I 82100 Benevento, Italy ⊥ Chemical and Food Engineering Department, Universidade Federal de Santa Catarina, 88040-970 Florianópolis, SC, Brazil ABSTRACT: The dynamic behavior of open- and closed-loop continuous bioreactors for cultures of substrate-inhibited aerobic microorganisms was studied. A proportional-integral feedback control scheme was adopted to stabilize the dissolved oxygen concentration by manipulating the dilution rate: the bioreactor operated as an oxystat. The effects of operating parameters (inlet substrate concentration and gas−liquid mass transfer rate) and of controller parameters (gain, reset time, and set-point) on solutions of mathematical models of both open- and closed-loop systems were characterized by means of parametric continuation technique and bifurcational analysis tools. The open-loop unstable steady state characterized by inhibiting levels of substrate in the bioreactor may be stabilized at any dissolved oxygen concentration by adopting a PI feedback controller. The bifurcation analysis provided useful guidelines to tune the control parameters. Stable periodic and multiperiodic regimes may appear through Hopf and Flip bifurcations, respectively, if gain and reset time are not properly set. Ranges of parameter values characterized by no solution were also assessed. Bifurcation maps of the gain vs the set-point value of oxygen concentration were proposed for different representative cases of gas−liquid mass transfer rate and inlet substrate concentration.

1. INTRODUCTION Rigorously controlled and reproducible conditions in bioreactors for laboratory operations are the prerequisite to collect reliable data related to both kinetics and stoichiometry of microorganism growth. Kinetics assessment based on data collected during batch cultures is affected by the continuous change in the environmental conditions (pH, biomass, substrate and metabolite concentrations, ...) inherent to this culturing method. Continuous steady-state operations provide a proper tool to characterize microorganism growth kinetics and stoichiometry;1−4 however, the operation of continuous bioreactors may be affected by the stability of the steady-state regimes.5−8 Cultures of microorganisms characterized by growth depending on just the concentration of the substrate and biomass (two-dimensional model) may be affected by both multiple steady states and periodic regimes. As a consequence, the bioconversion process frequently needs control strategies to ensure the prescribed environmental conditions. Unlike the chemostat, the auxostat (continuous cultures with classic feedback control) can operate at specific growth rates close to wash-out conditions to maximize biomass productivity.9 In particular, the control system becomes essential when the open-loop steady-state operation is unstable.10,11 Control strategies of bioreactors reported in the literature differ from one another in algorithms and manipulated and controlled variables.12−17 Notwithstanding some pioneering works highlighting the possibility to stabilize unstable steady states11,18−21 on the basis of theoretical analyses and numerical investigations of feedback control systems, only a few interesting experimental results regarding chemostats operated under unstable © 2013 American Chemical Society

steady-state conditions with different control strategies have been published so far. Some tests have actually been carried out in turbidostats,22,23 nutristats,24,25 and produstats.26 Online measurement devices of the feedback variable (substrate or biomass) are typically affected by several error sources, as bubble-swept effects of the probe, biofilm growth on the probe, clogging of the filter in the sampling line for computer-automated analysis. pH and DOT online measurements are quite reliable manipulated variables, thanks to the robustness of modern electrodes with respect to these error sources. Moreover, DOT is a good state variable, since oxygen is often a limiting substrate for aerobic processes. A comprehensive approach to describe the stability and nature of regime solutions is the systematic application of bifurcation analysis and continuation techniques.27−30 The bifurcation analysis has been extensively applied to characterize the behavior of nonlinear systems. Previous studies have been traditionally aimed at improving the understanding of the qualitative dynamics of the process units, for example.31 Over recent decades, bifurcation analysis has been proposed to address control synthesis problems in a rigorous framework. One of the best-known approaches in this area assumes that a control structure has been selected and it makes use of parametric continuation to identify regions of the controller parameter space that are characterized by unstable plant operation and undesired global dynamics, Received: Revised: Accepted: Published: 13422

March 10, 2013 August 9, 2013 August 11, 2013 August 12, 2013 dx.doi.org/10.1021/ie400782y | Ind. Eng. Chem. Res. 2013, 52, 13422−13431

Industrial & Engineering Chemistry Research

Article

as multistability.32−35 Regions of state multiplicity or instability are identified, thus enabling the choice of the controller parameter values and preventing transitions to undesired solution regimes or to stabilize a desired unstable regime. A further approach, also known as bifurcation control,36 incorporates information about the bifurcation characteristics of the uncontrolled plant into the synthesis of the control structure. In this context, control strategies capable of tailoring the bifurcation characteristics of nonlinear processes are formulated.37−41 This contribution is part of a research program on the characterization of microorganism kinetics by means of continuous cultures. Olivieri et al.42 successfully applied the closed-loop operation of a continuous bioreactor under stabilized steadystate conditions, adopting the oxystat configuration,43 with oxygen adopted as the controlled variable. They were able to assess the specific growth rate of an aerobic Pseudomonas strain as a function of phenol, the inhibiting substrate. Controller parameters (gain and reset time) and operating parameters (inlet substrate concentration and gas−liquid mass transfer coefficient) were accurately selected to stabilize the continuous culture at the desired set point by checking the stability of the state ex post. This task was made difficult because kinetic parameters were unknown. The success of the experimental campaign and the awareness that the experimental method can be applied to numerous aerobic cultures have fueled interest in a more rigorous theoretical investigation of controlled bioreactors. The aim of this work is the bifurcational analysis of a bioreactor model operated according to both a conventional open-loop chemostat and a novel closed-loop oxystat. Growth kinetics and stoichiometry were assumed to be consistent with Haldane44 and Pirt45 models, respectively. The potential of the closed-loop oxystat coupled with the simple classic unstructured model was pointed out by Olivieri et al.42 as regarding the characterization of substrate-inhibited growth kinetics of microorganisms. In this article, the analysis is based on parametric continuation of regime solutions to characterize the bifurcational scenario of both openand closed-loop systems. The effects of control and operating parameters are described by means of solution and bifurcation diagrams. A picture of the steady and dynamic regimes is proposed to support experimental investigations. A thorough description of the steady-state and dynamic regimes is reported for the operating and key control parameters typically adopted in experimental investigations.

V. The specific growth rate of the microorganism is described by a model relationship characterized by Haldane-type substrate inhibition and Monod-type oxygen limitation: μ(S , O2 ) = μ M

S KS + S +

S2 KI

O2 K O2 + O2

(1)

2.2. Model Equations. The mass balances for biomass, substrate, and dissolved oxygen extended to the liquid volume of the continuous stirred tank reactor (CSTR) are dX = −DX + μ(S , O2 ) ·X dt

(2a)

μ(S , O2 ) dS = D(S IN − S) − ·X dt YX /S

(2b)

dO2 μ(S, O2 ) = (D + KLaL)(O2Eq − O2 ) − ·X dt YX /O2

(2c)

where D is the dilution rate, μ the specific growth rate, KLaL the volumetric gas−liquid oxygen transfer coefficient referred to the liquid volume unit, and YX/S and YX/O2 are the substrate and oxygen to biomass yields, respectively. Substrate and oxygen uptake for cell maintenance is included in the stoichiometry of the process.45 Accordingly, the biomass yields are described according to the relationships: m 1 1 = M + S YX /S μ YX /S (3a) 1 YX /O2

=

YM X/S

1 YXM/O2

+

mO2 μ

(3b)

YM X/O2

where and are the maximum values of oxygen to biomass and substrate to biomass yields of YX/S and YX/O2 and mS and mO2 are the substrate maintenance coefficient and the oxygen maintenance coefficient. The set of eqs 1, 2, and 3 is an ordinary differential equation model, and the phase space is R3. The open-loop dynamics of the chemostat may be simulated by the resolution of this set of equations. The model of the closed-loop dynamics of the chemostat includes the equation set 1−3 and a function describing the PI feedback control based on the state variable O2. In previous works, the authors had investigated the oxystat configuration by manipulating the inlet substrate concentration.46 In this study, the control of the dilution rate eq 4 is adopted:

2. THEORETICAL FRAMEWORK The model is based on the mass balance for (carbon source) substrate, oxygen, and cells coupled with the constitutive equations. Two systems are proposed: an open-loop bioreactor and a closedloop bioreactor. The latter is the open-loop bioreactor equipped with a feedback controller. 2.1. Model Assumptions. The model relies on the following assumptions: I. The reactor is well-mixed with respect to the liquid phase; II. The reactor is continuously fed with a sterile liquid stream bearing (carbon source) substrate at constant concentration (SIN); III. The oxygen concentration in the feeding stream is in equilibrium with air (OEq 2 ); IV. An air stream is continuously fed into the reactor at a constant flow rate;

⎡ ⎛ OSET − O 1 D(t ) = DS⎢1 + K C⎜ 2 SET 2 + ⎢⎣ TR ⎝ O2

∫0

t

OSET − O2 2 OSET 2

⎞⎤ d t ⎟⎥ ⎠⎥⎦ (4)

where OSET 2 is the set-point of the O2, KC is the controller gain, TR is the controller reset-time, and DS is the steady-state value of D at O2 = OSET 2 . Model equations are rewritten in dimensionless form, introducing the following dimensionless variables:

13423

S=

O2 S X ; ; X = IN M ; O2 = Eq IN S S YX /S O2

μ=

μ ; t = t ·μ M /(2 KS/KI + 1) μ M /(2 KS/KI + 1)

(5)

dx.doi.org/10.1021/ie400782y | Ind. Eng. Chem. Res. 2013, 52, 13422−13431


Article

Table 1. Dimensionless Mathematical Model

Dynamic simulations of the closed-loop system were carried out to investigate the transient response of the bioreactor in a few selected cases. The initial conditions are:

The specific growth rate and the time are made dimensionless with reference to the maximum specific growth rate, μ = μM/ (2(KS/KI)1/2 + 1), provided that O2 ≫ KO2 A dimensionless state variable, θ, is also introduced to assess the dynamics of the closed-loop system. It is defined as θ=

∫0

t


⎧ S = S0 ≥ 0 ⎪ ⎪ t = t 0 → ⎨ O2 = O20 ≥ 0 ⎪ ⎪ ⎩ X = X0 ≥ 0

dt (6)

In particular, θ is the integral term of the control action in eq 4 in dimensionless form: Da(t) =

Under steady-state conditions, the left side of the differential eqs T.1.1−T.1.6 is equal to 0, and each systemopen-loop system and closed-loop systemis reduced to a set of nonlinear integral−algebraic equations. Steady-state solutions of the equation set are denoted as XS, SS and O2S. 2.3. Computational Method. Bifurcational analysis was performed using a continuation algorithm provided by the toolbox MATCONT 47 associated with MATLAB 2010 implementing both the open-loop system (T.1.1−T.1.3) and the closed-loop system (T.I.4−6). Solution diagrams of the state variables XS, SS, and O2S for 1D parametric continuation and bifurcation diagrams for 2D parametric continuation were adopted to describe the simulated dynamics. The solution diagrams report the regime solution values for the selected state variable versus the bifurcation parameter for static regimes. The maximum and minimum values recorded during oscillations are used to characterize periodic regimes. The bifurcation study is carried out adopting the Damkhöler group, the control gain, and the oxygen set point as bifurcation parameters. The convention adopted for the solution and bifurcation diagrams is the following: solid lines represent stable stationary solutions; dashed lines, unstable stationary solutions; solid circles, stable limit cycles; open circles, unstable limit cycles. The bifurcation diagram (also known as a bifurcation map) illustrates the loci of bifurcation points (saddle-node, transcritical, Hopf, etc.) in the parameter space.48 Some simulation results are also reported to illustrate interesting situations.

DaS

(

KC


+

1 θ ;R

)+1

(7)

Table 1 reports the dimensionless form of the model of both the open-loop system (T.1.1−T.1.3) and the closed-loop system T.1.4−T.1.6. Equation T.1.7 is the dimensionless form of the kinetic growth rate eq 1. The dimensionless groups are listed in Table 2 along with the values used in the computations. In particular, some data come from a previous experimental investigation by Olivieri et al.29,42 An auxiliary equation is also coupled to the set T.1.4−T.1.6 OSET − O dθ = 2 SET 2 dt O2

(9)

(8)

As a consequence, the dynamics of the closed-loop system is described by the equation set T.1.4−T.1.6 and eq 8, an ODE problem in the phase space R4. The range of the dimensionless groups Da, 2L , α, KC, ;R and OSET 2 are also reported in Table 2. The group α is set at a value >1 to focus exclusively on the growth-inhibiting effect of the substrate. The dimensionless variables X, S, and O2 vary in the interval [0, 1] because both are normalized against the maximum value attainable. 13424



Article

Table 2. Values of Input Parameters of the Computations model parameter

dimensionless parameter

μM KS KI

0.26 h−1 a 0.005 g/La 0.20 g/La

K O2

10−4 g/Lb

YM X/S YM X/O2

β=

KI KS

6.32

K O2

0.0128

2O2 =

O2Eq

1.27 g/ga 0.17 g/gb

γS/O = 2

OEq 2 mS

m O2

D

30.3

YXM/O2O2Eq

−3

7.8 × 10 g/L 0.017 g/(g h)a

mSYXM/S

4S =

μ /(2 KS/KI + 1)

4 O2 = −1

0.198−1980 h

SIN

0.032−3.2 g/L

KLaL

1.98−1980 h−1

KC TR

0.1−10 (1.98 × 10−4)−19.8

OSET 2

0−7.8 × 10−3 g/L

0.109

M

0.040 g/(gh)c

mO2YXM/O2

0.348

M

μ /(2 KS/KI + 1) M

Da =

a

YXM/S KSKI

μ /(2 KSKI + 1)

100−104

D 100−102

IN

α=

S KSKI

2L =

KLaL μ /(2 KSKI + 1)

KC

;R = TR μ M (2 KSKI + 1)

OSET = 2

101−104

M

OSET 2

0.1−10 10−1−101 0−1

O2Eq

Olivieri et al.42 bOlivieri et al.29 cmO2 = mSYOStoich , where YOStoich = 2.38 g/g. 2/S 2/S

3. RESULTS AND DISCUSSION 3.1. Open-loop system − chemostat. 3.1.2. Dynamics. The dynamics of the system T.1.1−T.1.3 in terms of the state variable X(t), S(t), and O2(t) depends on the value of the three dimensionless parameters Da, 2L , and α. Figure 1 reports the solution diagram of the state variables as a function of Da at 2L = 200 and α = 5. It is possible to recognize the typical behavior of continuous cultures of cells characterized by the Haldane-type specific growth kinetics.10,13,29,42,49 In particular, three solution branches exist, denoted as 1, 2, and 3, respectively. At low Da, a single solution (3) characterized by stable steady state is observed: it is the trivial wash-out solution (XS = 0, SS = 1, and O2S = 1). For 1 < Da < 1.4, the three solutions (1, 2, and 3) are observed. The wash-out regime (solution 3) coexists with two regimes characterized by SS < 1 as a consequence of both the transcritical (TR2−3) and the saddle-node (SN1−2) bifurcations. At Da > 1.4, two solutions are physically consistent: solutions 1 and 2. The wash-out solution becomes unstable at the transcritical bifurcation characterized by Da = 1.4, and the new stable solution 1 appears. According to the number of solutions, three intervals of Da are identified: (I) solution 3 (stable) is present; (II) solution 1 (stable) and solution 3 (unstable) are present; (III) solution 1 (stable), solution 2 (unstable), and solution 3 (stable) are present. Under the operating conditions investigated, O2S is always higher than KO2 (Figure 1C). It appears that the cell growth rate does not significantly depend on oxygen concentration; therefore, the SS-vs-Da plot may look like the plot of μ vs S, typical of the

Haldane kinetic model. Under these conditions, the saddle-node bifurcation is characterized by S = 1/α, as expected. Effects of both 2L and α on the bifurcation scenarios were investigated. Solution diagrams assessed at 2L and α ranging in the intervals reported in Table 2 confirm the dynamic behavior described in Figure 1, except for the Da at which bifurcations occurred. Table 3 reports the Da at the saddle-node and transcritical bifurcations for different values of 2L and α. The Da of the saddle-node bifurcation increases with α and decreases with 2L . The Da of the transcritical bifurcation increases with α increases and does not depend on 2L . Altogether, the Da range of existence of the three steady-state solutions widens as both α and 2L increase. Figure 2 reports SS vs O2S for the nontrivial steady states 1 and 2 at different values of 2L and α. The effects of aeration on the chemostat sensitivity to both substrates (carbon source and oxygen) may be analyzed. The plot in Figure 2 for 2L = 1000 and α = 5 shows that the nontrivial solutions are characterized by a narrow interval of O2S, close to the saturation level. The range of values of O2S increases when α increases and 2L decreases. In particular, for 2L = 200 and α = 20 (Figure 2), it clearly appears that the culture is operated under oxygen limitation conditions because of the high carbon source conversion. The results discussed in the previous paragraphs are particularly relevant to the definition of the strategy of closedloop operations. The higher the O2 setting, the higher the capacity of the controller to keep the O2 at the set-point. In other words, the large interval of O2S reflects a higher sensitivity of 13425



Article

Figure 2. Open-loop system: O2-vs-S phase portrait of steady states 1 and 2 as a function of 2L and α.

Figure 3. Open-loop system: Da-vs-2L bifurcation diagram at α = 20. Figure 1. Open-loop system: solutions diagram as a function of Da, 2L = 200, and α = 5.

Region II : the parameter plane region over the TR2−3 curve characterized by solution 1 (stable) and solution 3 (unstable); Region III : the region bounded by both SN1−2 and TR2−3 curves and embodying the parameter sets consistent with solutions 1 (stable), 2 (unstable), and 3 (stable). The extent of region III increases as 2L and α increase. 3.2. Closed-Loop System: Oxystat. 3.2.1. Dynamics. The analysis of the dynamic features of the closed-loop system is reported in this section. 2L and α are the operating parameters, and O2SET, KC, and ;R are the controller parameters. It is worth observing that according to eq 8, the Damkhöler number depends on state variables X, S, O2, and θ for the closed-loop system and not on input parameters, as it is for open-loop systems. Figure 4 reports the solutions diagram of X, S, and O2 (left column) as a function of O2SET at 2L = 200, α = 5, KC = 1, and ;R = 1. Results of DaS are also reported (right column) as X vs DaS, S vs DaS, and O2 vs DaS phase portraits for a direct comparison with the solution diagrams reported in Figure 1 for the open-loop system. In line with the results reported in Figure 2, the steady-state solutions 1 and 2 are still present for the closed-loop system as O2SET > 0.42 (region II). No solution is found for O2SET < 0.42 (region 0). The comparison of the plots X vs Da, S vs Da, and O2 vs Da represented for open-loop and closed-loop systems shows that

Table 3. Open-Loop System: Da at the Saddle-Node and Transcritical Bifurcations at Different Values of 2L and α α=5

2L = 50 2L = 200 2L = 1000

α = 20

α = 50

SN

TR

SN

TR

SN

TR

1.252 1.028 1.014

1.402 1.402 1.402

2.805 1.874 1.027

3.210 3.210 3.210

5.856 3.806 1.382

6.858 6.858 6.858

the controlled variable with respect to the manipulated variable (Da). 3.1.2. Bifurcation diagram. The bifurcation diagram is a synoptic view of the dynamics and steady-state solutions reported in the previous section. Figure 3 reports the loci of both TR2−3 and SN1−2 bifurcations in the parameter plane (2L , Da) for α = 20. As 2L decreases, the saddle-node bifurcation curve and the transcritical curve merge into a cusp point (CP1−2−3). Three regions bounded by SN1−2 and TR2−3 curves may be identified in the parameters plane. They are listed as follows, according to their Da intervals introduced in Figure 1: Region I : the region includes a set of parameters below SN1−2 and TR2−3 curves and consistent with solution 3 (stable) only; 13426



Article

On the other hand, solution 3 is always unstable. Dynamic simulations carried out by setting initial conditions close to wash-out (S0 ∼ 1, O20 ∼ 1, X0 ∼ 0) show that batchwise conditions are reached in about t = 10: S and O2 approach 0 and 1, respectively, whatever the initial conditions (Figure 5A); the asymptotic value of X (Figure 5A and B) depends on the initial conditions (Figure 5B). In particular, the transient regime reported is characteristic of a large number of simulations with slightly different initial conditions. This result can be considered quite satisfactory, since the dimension of the closed-loop system (four) does not allow drawing the basin of attraction of each stable solution. This behavior is due to the integral term of the control action, which is responsible for the nonautonomous character of the closed-loop system. Figure 6 reports solutions 1 and 2 as a function of KC, 2L and α at O2SET = 0.9 and ;R = 1. The analysis of Figure 6A shows that regime 2 is stable for KC > 0.674 (region II). At KC = 0.674, a supercritical Hopf bifurcation appears. For KC < 0.674 (region III), regime 2 becomes unstable, and a stable periodic regime, named 4, emerges. The period of the oscillations (T4) is reported in Figure 6A. T4 increases continuously from 25 to 125 when KC decreases from 0.674 to 0.1. The effects of 2L on the dynamic features of closed-loop systems may be deduced by analyzing the solution diagram reported in Figure 6B. Stable oscillating regimes are found for 2L belonging to the interval 301−1076 (region III). The onset of stable oscillating regimes is characterized by supercritical Hopf bifurcations. The oscillation period T4 does not change significantly with 2L and it is ∼20. Under the operating conditions investigated, no periodic solutions are found for α ranging between 1 and 50 (Figure 6C). Figure 7 reports an interesting scenario found through simulations carried out for 2L = 1000, α = 20, KC = 1, and ;R = 1. It should be observed that the pair of values 2L = 1000, α = 20 is characterized by a solution diagram, SS vs O2S, for open-loop systems (Figure 2) that does not strongly differ from that drawn for 2L = 200, α = 5. However, the solution scenarios shown in Figures 4 and 7 are very different. As expected, no solution is found for O2SET < 0.48 (region 0). Steady-state solutions 1 (unstable) and 2 (stable) are found for O2SET > 0.48. At O2SET = 0.631, the stable solution 2 becomes unstable, and a periodical regime appears (4) due to a supercritical Hopf bifurcation (H2−4). As O2SET approaches 0.682, the period of solution 4 dramatically increases. Solution 4 vanishes at O2SET = 0.682 under the effect of a global homoclinic bifurcation (Hom4).

Figure 4. Closed-loop system assessed at 2L = 200, α = 5, KC = 1, and ;R = 1. (A) Solution diagram and (B) phase portrait of selected variables.

the stability of steady-state solutions 1 and 2 is reversed. Indeed, the stable (unstable) solution branch 1 (2) in the open-loop system (Figure 1) becomes unstable (stable) under controlled conditions. The trivial solution found for the open-loop system (branch 3 in Figure 1) deserves further discussion. Reactor wash-out is prevented by the control action, even though it is theoretically possible with the closed-loop system. On one hand, the integral term of the control action (eq 7) yields a continuous increase of Da at S = 1, O2 = 1, and X = 0: Da → +∞ for O2 = 1. Therefore, a time discontinuity may be observed because the reactor starts to be operated in batchwise mode.

Figure 5. Closed-loop system: state variable X, S, O2 vs t. 2L = 200, α = 5, O2SET=0.6, KC = 1 and ;R = 1. (A) X0 = 10−4, S0 = 1, O20 = 1 and (B) X0 = 10−4, O20 = 1. 13427



Article

Figure 6. Closed-loop system. Solution diagrams assessed at O2SET = 0.9 and ;R = 1. Periodic solutions: the maximum and minimum values of oscillations are represented by dot couples (A). ;R = 200, α = 5, (B) α = 5, KC = 1, and (C) 2L = 200, KC = 1.

Four regions may be identified by the SNI−II, HII−IV, and HomIV bifurcation curves, each of them characterized in terms of number of solutions and their stability. Region 0 :The closed-loop operation is always unsuccessful, since there is no solution. According to Figure 4B, the set-point O2SET is too low with respect to the potential of the reactor. Under these operating conditions, there is no DaS, which can result in the desired O2 = O2SET (see Figure 2 related to the open-loop operation). If the O2SET is selected in this region and an arbitrary/incorrect value of DaS is chosen, the system will approach a steady state characterized by the lowest admissible values of O2 and Da valid for the open-loop operation. Under these conditions, the control will fail because O2 ≠ O2SET and Da ≠ DaS. Region II :The region is characterized by the two steady-state solutions, 1 (unstable) and 2 (stable). Region III :The region is characterized by the three solutions: steady-state solutions 1 and 2 and periodic solution 4. The closed-loop operation is not successful because both the steady-state solutions 1 and 2 are unstable. The

A detailed analysis of the dynamics of the closed-loop system is reported in a zoom-in of the solution diagram (Figure 8). At O2SET = 0.6803, a homoclinic bifurcation (Hom5) introduces a 2-periodic oscillating solution (5). The stability features of solution 4 are complicated by the subcritical period doubling bifurcation detected at O2SET = 0.6808 (PDIV−V) before the Hom4 bifurcation point. An unstable 2-periodic regime (5) is present for O2SET < 0.6808, and solution 4 becomes unstable before vanishing at the Hom4 point. 3.2.2. Bifurcation diagram. Figure 9 reports the bifurcation diagrams in the O2SET vs KC space for ;R = 1. The attention is focused on the dynamics-assessed setting 2L and α at values within the practical interval.50 According to the typical values reported in the literature, the pairs of KLaL and SIN are (200, 5), (200, 20), (1000, 5), and (1000, 20). Stability regions are identified with reference to SN1−2, H2−4, and Hom4 bifurcation curves. The bifurcation curves Flip4−5 and Hom5 are not reported for two reasons: the limited range of existence and the unstable nature of the reactor when operated within these curves. 13428



Article

Figure 9. Closed-loop system: KC vs O2SET bifurcation diagram for ;R = 1.

The simulations carried out by changing ;R in the range 0.1− 10 show that the bifurcation map does not change appreciably. The ;R investigated interval was not extended further because it reflected the behavior of the controller: the integral and the proportional behavior dominate at ;R = 0.1 and 10, respectively. It is worth observing that in line with eq 4 at ;R > 10, the integral action in the controller is negligible, and at ;R < 0.1, the integral action is predominant. 3.2.3. Controller Operation Guidelines. Region B increases as α increases and 2L decreases. The first parameter is related to the maximum oxygen demand; the second, to the oxygen pumping into the reaction volume. The increase in α and the decrease in 2L reduce the minimum (lower boundary of the interval) of the oxygen concentration in the bioreactor, as reported in Figure 2 for open-loop systems. The minimum of O2, that is, the SNI−II bifurcation, decreases as α increases and 2L decreases. These observations would suggest that the best selection of operating conditions is high α and low 2L to operate the reactor under the steady solution 2. Moreover, the lowest possible value of the O2SET should be selected to reduce the risk of oscillations (solution 4). The results of the dynamic analysis of the closed-loop system suggest that particular attention should be paid to the selection of operating conditions, including controller parameters, whenever aerobic bacteria characterized by substrate-inhibited growth kinetics are continuously cultivated in a bioreactor. The bioreactor can be successfully operated under steady-state conditions if the operating conditions fall in region II, where solution 2 is stable. An example of successful application of the closed-loop system is the experimental work by Olivieri et al.42 The main suggestions for the selection of the operating conditions are reported in the following: (a) A high value of inlet substrate concentration and low value of gas−liquid mass transfer rate should be adopted to extend the range of admissible values for the oxygen concentration set-point. (b) The set-point of the oxygen concentration should be as low as possible but larger than the minimum value associated with oxygen limitation. This selection is required when the

Figure 7. Closed-loop system: solutions diagram for 2L = 1000, α = 20, KC = 1, and ;R = 1.

Figure 8. Closed-loop system: inset of O2 vs O2SET solution diagram of Figure 7.

periodic solution 4 is stabilized by the control action, but the oscillation amplitude is quite large, more than five times the characteristic time-scale of the system. Region IV : No stable solution is present. Steady-state solutions 1 and 2 are unstable. The period of the periodic solution 4 increases so much that the periodic solution vanishes through the homoclinic bifurcation Hom4 with the unstable solution 1. 13429



Article

KI = kinetic constant; g/L KLaL = gas−liquid mass transfer coefficient; h−1 2L = gas−liquid mass transfer coefficient; dimensionless KS = substrate kinetic constant; g/L KO2 = oxygen kinetic constant; g/L KO2 = oxygen kinetic constant; dimensionless μ = specific growth rate; h−1 μM = kinetic constant; h−1 mO2 = oxygen maintenance coefficient; h−1 4 O2 = mO2YXM/O2(2 KS/KI )/μ M = oxygen maintenance coefficient; dimensionless mS = substrate maintenance coefficient; h−1 4S = mSYX /S M(2 KS/KI + 1)/μ M = substrate maintenance coefficient; dimensionless O2 = oxygen concentration; g/L O2 = O2/OEq 2 = oxygen concentration; dimensionless OEq = oxygen concentration in equilibrium with air; g/L 2 OSET = set-point value of oxygen concentration; g/L 2 Eq OSET = OSET 2 2 /O2 = set-point value of oxygen concentration; dimensionless θ = integral error S = substrate concentration; g/L SIN = inlet substrate concentration; g/L S = S/SIN = substrate concentration; dimensionless t = time; h t = time; dimensionless T = period of oscillating solutions; dimensionless TR = controller reset time; h ;R = TR μ M /(2 KS/KI + 1) = controller reset time; dimensionless X = biomass concentration; g/L X = X/SINYM X/S = biomass concentration; dimensionless YX/O2 = oxygen-to-biomass true growth yield YM X/O2 = oxygen-to-biomass maximum growth yield YX/S = substrate-to-biomass true growth yield YM X/S = substrate-to-biomass maximum growth yield

parameters of the substrate-inhibition kinetics have to be assessed. (c) The controller gain should be >1 to avoid oscillating conditions. (d) The reset time should be ∼1 to balance the proportional and the integral actions of the controller. It should be pointed out that the basin of attraction of stable solutions does not coincide with the phase space. Indeed, whenever the closed-loop bioreactor is operated at initial conditions close to the trivial solution 3, it switches to batchwise operation. The biomass grows until complete depletion of the substrate, and the final value of the biomass concentration depends on the initial conditions.

4. CONCLUSIONS The dynamic analysis of a continuous aerobic bioreactor was carried out. The study concerned both a traditional bioreactor and a bioreactor equipped with a controller of the feeding stream rate, known as open- and closed-loop systems, respectively. The controller modulated the feeding stream rate as a function of the dissolved oxygen concentration. The effects of substrate concentration in the bioreactor feeding, gas−liquid mass transfer rate, controller gain, reset time, and set point of the controlled variable were characterized. Solution diagrams of state variables as a function of the continuation parameter selected were assessed. Bifurcation diagrams were produced to select the best set of values for the operating parameters. The results show that the steady-state conditions, unstable when operated in the open-loop system, can be stabilized in closed-loop system. The bifurcation analysis suggests that the success of the control asks for an optimal selection of both the gain and the reset time of the controller. In particular, these parameters should be selected to avoid oscillating solutions as well as the absence of solutions for the continuous operation. The investigation confirms the effectiveness of the control strategy adopted when continuous cultures are utilized to assess the kinetics and stoichiometry of microorganism growth. Experimental investigations on nonlinear growth kinetics (e.g., substrate inhibition27) can be extensively adopted for different species/substrates to fulfill the constant demand of reliable kinetic data to support the design of novel bioprocesses. It is evident that nonlinear analysis tools offer a useful support to define proper conditions of the experimental campaign, reducing experiments’ failure and false responses due to hidden periodic regimes.

■

Superscript

S steady-state value

■

REFERENCES

(1) Owens, J. D.; Legan, J. D. Determination of the Monod substrate saturation constant for microbial growth. FEMS Microbiol. Rev. 1987, 46, 419−432. (2) Fraleigh, S. P.; Bungay, H. R.; Clesceri, L. S. Continuous culture, feedback control and auxostats. Trends Biotechnol. 1989, 7, 159−164. (3) Napoli, F.; Olivieri, G.; Russo, M. E.; Marzocchella, A.; Salatino, P. Continuous lactose fermentation by Clostridium acetobutylicum − Assessment of acidogenesis kinetics. Bioresour. Technol. 2011, 102, 1608−1614. (4) Napoli, F.; Olivieri, G.; Russo, M. E.; Marzocchella, A.; Salatino, P. Continuous lactose fermentation by Clostridium acetobutylicum − Assessment of energetics and product yields of the acidogenesis. Enzyme Microb. Technol. 2012, 50, 165−172. (5) Yano, T.; Koga, S. Dynamic behavior of the chemostat subject to substrate inhibition. Biotechnol. Bioeng. 1969, 11, 139−153. (6) Crooke, P. S.; Wei, C. J.; Tanner, R. D. The effect of specific growth rate and yield expression on the existence of oscillatory behavior of a continuous fermentation model. Chem. Eng. Commun. 1980, 6, 333− 342. (7) Agrawal, P.; Lee, C.; Lim, H. C.; Ramkrishna, D. Theoretical investigations of dynamic behaviour of isothermal continuous stirred tank biological reactor. Chem. Eng. Sci. 1982, 37 (3), 453−462.

AUTHOR INFORMATION

Corresponding Author

*Tel: +39 081 7682541. Fax: +39 081 5936936. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

NOMENCLATURE α = SIN/(KSKI)1/2 = inlet substrate concentration; dimensionless β(KIKS)1/2 = substrate kinetic constant; dimensionless D = dilution rate; h−1 Da = μM/D(2(KSKI)1/2 + 1) = Damkhöler number 1/2 M Eq γS/O2 = (YM X/S(KSKI) )/(YX/O2O2L ) = substrate-to-oxygen yield; dimensionless KC = controller gain 13430



Article

(8) Ajbar, A. On the existence of oscillatory behaviour in unstructured models of bioreactors. Chem. Eng. Sci. 2001, 56, 1991−1997. (9) Gostomoski, P.; Mühlemann, M.; Lin, Y. H.; Mormino, R.; Bungay, H. Auxostats for continuous culture research. J. Biotechnol. 1994, 37, 167−177. (10) Edwards, V. H. The influence of high substrate concentrations on microbial kinetics. Biotechnol. Bioeng. 1970, 12, 679−712. (11) Edwards, V. H.; Ko, C. R.; Balogh, A. Dynamics and control of continuous microbial propagators to subject substrate inhibition. Biotechnol. Bioeng. 1972, 14, 939−974. (12) Agrawal, T. H.; Lim, H. C. Analyses of various control schemes for continuous bioreactors. Adv. Biochem. Eng. Biot. 1984, 30, 61−90. (13) Suzuki, S.; Shimizu, K.; Matsubara, M. On the parameter spaceclassification of the dynamic behaviour of a continuous microbial flow reactor. Chem. Eng. Commun. 1985, 33, 325−335. (14) Tsao, J. H.; Wu, W. T. Global control of continuous stirred tank bioreactor. Chem. Eng. J. 1994, 56, B59−B64. (15) Zhao, Y.; Skogestad, S. Comparison of various control configurations for continuous bioreactors. Ind. Eng. Chem. Res. 1997, 36, 697−705. (16) Lim, H. C., Lee, K. S. Control of bioreactor system. In Biotechnology; 2nd Ed.; Rhem, H. J., Reed, G., Eds.; Wiley: New York, 2001. (17) Tian, Y.; Sun, K.; Chen, L. S.; Kasperski, A. Studies on the dynamics of continuous bioprocess with impulsive state feedback control. Chem. Eng. J. 2010, 157, 558−567. (18) Di Biasio, D. H.; Lim, H. C.; Weigand, W. A.; Tsao, G. T. Phaseplane analysis of feedback control of unstable steady states in a biological reactor. AIChE J. 1978, 24, 686−693. (19) Chen, L. H.; Chang, H. S. 1984. Global stabilization of a biological reactor by linear feedback control. Chem. Eng. Commun. 1984, 27, 231− 254. (20) Chang, H. C.; Chen, L. H. Bifurcation characteristic of nonlinear system under conventional PID control. Chem. Eng. Sci. 1984, 7 (8), 1127−1142. (21) Shimizu, K.; Matsubara, M. Conditions for the phase-plane analysis of feedback control of chemostat. Biotechnol. Bioeng. 1985, 27, 519−524. (22) Di Basio, D. H.; Lim, H. C.; Weigand, W. A. An experimental investigation of stability and multiplicity of steady states in biological reactor. AIChE J. 1981, 27, 284−290. (23) Jayakumar, S.; Lim, H. C. Multiple steady states of Methylomonas mucosa for continuous production of polysaccharides. J. Biotechnol. 1989, 12, 21−36. (24) Rutgers, M.; Bogte, J. J.; Breure, A. M.; Van Andel, J. G. Growth and enrichment of pentachlorophenol-degrading microorganisms in the nutristat, a substrate concentration-controlled continuous culture. Appl. Environ. Microbiol. 1993, 59, 3373−3377. (25) Rutgers, M.; Gooch, D. D.; Breure, A. M.; Van Andel, J. G. Assessment of inhibition kinetics of the growth of strain P5 on pentachlorophenol under steady state conditions in a nutristat. Arch. Microbiol. 1996, 165, 194−200. (26) Schröder, M.; Müller, C.; Posten, C.; Deckwer, W. D.; Hecht, V. Inhibition kinetics of phenol degradation from unstable steady state data. Biotechnol. Bioeng. 1997, 54, 567−576. (27) Zhang, Y. C.; Henson, M. A. Bifurcational analysis of continuous bioreactor models. Biotechnol. Prog. 2001, 17, 647−670. (28) Russo, M. E.; Maffettone, P. L.; Marzocchella, A.; Salatino, P. Bifurcational and dynamical analysis of a continuous biofilm reactor. J. Biotechnol. 2008, 135, 295−303. (29) Olivieri, G.; Russo, M. E.; Marzocchella, A.; Salatino, P. Modelling of an aerobic biofilm reactor with double-limiting substrate kinetics: bifurcational and dynamical analysis. Biotechnol. Prog. 2011, 27, 1599− 1613. (30) Ma, Y. F.; Xiu, Z. L.; Sun, L. H.; Feng, E. M. Hopf bifurcation and chaos analysis of a microbial continuous culture model with time delay. Int. J. Nonlinear Sci. Numer. Simul. 2011, 7 (3), 305−308.

(31) Mancusi, E.; Merola, G.; Crescitelli, S.; Maffettone, P. L. Multistability and hysteresis in an industrial ammonia reactor. AIChE J. 2000, 46 (4), 824−828. (32) Alhumaizi, K.; Elnashaie, S. E. Effect of control loop configuration on the bifurcation behaviour and gasoline yield of industrial fluid catalytic cracking (FCC) units. Math. Comput. Mod. 1997, 25, 37−56. (33) di Bernardo, M. Bifurcation analysis for control system applications. In Bifurcation Control: Theory and Applications; SpringerVerlag: Berlin, 2003. (34) Hahn, J.; Monnigmann, M.; Marquardt, W. A method for robustness analysis of controlled nonlinear systems. Chem. Eng. Sci. 2004, 59, 4325−4338. (35) Zavala-Tejeda, V.; Flores-Tlacuahuac, A.; Vivaldo-Limac, E. The bifurcation behavior of a polyurethane continuous stirred tank reactor. Chem. Eng. Sci. 2006, 61, 7368−7385. (36) Chen, G.; Moiola, J. L.; Wang, H. O. Bifurcation control: theories, methods and appolications. Int. J. Bifurcation Chaos Appl. Sci. Eng. 2000, 10 (3), 511−548. (37) Abed, E. H.; Fu, J. H. Local Feedback stabilization and bifurcation control I. Hopf bifurcation. Syst. Control Lett. 1986, 7, 11−17. (38) Wang, H. O.; Abed, E. H. Bifurcation control of a chaotic system. Automatica 1995, 31, 1213−1226. (39) Wang, X.; di Bernardo, M.; Stoten, D. P.; Lowenberg, M. H.; Charles, G. Bifurcation tailoring via Newton flow-aided adaptive control. Int. J. Bifurcation Chaos Appl. Sci. Eng. 2003, 13, 677−684. (40) Altimari, P.; Russo, L.; Mancusi, E.; Di Bernardo, M.; Crescitelli, S. Optimal reference trajectory shaping and robust gain-scheduling for transition control of nonlinear processes. Ind. Eng. Chem. Res. 2009, 48 (20), 9128−9140. (41) Altimari, P.; Mancusi, E.; Di Bernardo, M.; Russo, L.; Crescitelli, S. Tailoring the bifurcation diagram of nonlinear dynamical systems: An optimization based approach. Int. J. Bifurcation Chaos Appl. Sci. Eng. 2010, 20 (4), 1027−1040. (42) Olivieri, G.; Russo, M. E.; Di Donato, A.; Marzocchella, A.; Salatino, P. Unstable steady state operations of substrate inhibited cultures by dissolved oxygen control. J. Biotechnol. 2010, 156, 302−308. (43) Hospodka, J. Oxygen-absorption rate-controlled feeding of substrate into aerobic microbial cultures. Biotechnol. Bioeng. 1966, 8, 117−134. (44) Haldane, J. B. S. Enzymes; Longmans, Green and Co.: London, 1930. (45) Pirt, S. J. The maintenance energy of bacteria in growing cultures. Proc. R. Soc. B 1965, 163, 224−231. (46) Olivieri, G.; Russo, M. E.; Mancusi, E.; Marzocchella, A.; Maffettone, P. L.; Salatino, P. Non-linear analysis of feedback controlled aerobic cultures. Chem. Eng. Trans. 2013, 32, 811−816. (47) Dhooge, A.; Govaerts, W.; Kuznetsov, Yu. MATCONT: A MATLAB package for numerical bifurcation analysis of ODEs. ACM T Math Software 2003, 29, 141−164. (48) Kubíček, M.; Marek, M. Computational Methods in Bifurcation Theory and Dissipative Structure; Springer Verlag: New York, 1983. (49) Viggiani, A.; Olivieri, G.; Siani, L.; Di Donato, A.; Marzocchella, A.; Salatino, P.; Barbieri, P.; Galli, E. An airlift biofilm reactor for the biodegradation of phenol by Pseudomonas stutzeri OX1. J. Biotechnol. 2006, 123 (4), 464−477. (50) Chisti, Y Mass Transfer. In Encyclopedia of Bioprocess Technology; Flickinger, M. C., Drew, S. W., Eds.; Wiley: New York, 1999; pp. 1607− 1640.

13431


Nonlinear Analysis of Substrate-Inhibited Continuous Cultures

Recommend Documents