Feedback Control of a Continuous-Flow Stirred Tank Reactor with

Ind. Eng. Chem. Res. 2003, 42, 3765-3785

3765

Feedback Control of a Continuous-Flow Stirred Tank Reactor with Competing Autocatalators Woraritt Chaivorapoj, I4 nanc¸ Birol, Ali C ¸ ınar, and Fouad Teymour* Chemical and Environmental Engineering Department, Illinois Institute of Technology, Chicago, Illinois 60616

Two types of nonlinear feedback control schemes are introduced and analyzed for their capability of recovering the original state of an isothermal continuous-flow stirred tank reactor with one robust cubic autocatalytic species, perturbed by a temporary disturbance of an invading cubic autocatalytic species in the inflow. The control objectives are to eliminate the invading species from the system and to restore the original state of the host species. The extent of applicability of the control design to different nonrobust invading species is studied, when the controller is tuned for a specific invader. Moreover, a time-delay feature is suggested in one of the control schemes developed to achieve the control objectives in systems with poor detection of invading species. 1. Introduction An autocatalytic reaction

A + nB f (n + 1)B

(1)

is a chemical process where the presence of the product, B, affects its own production rate. It can represent both biological processes, such as reproduction of single cells by division1-5 or sexual reproduction in higher organisms,4,5 and chemical processes such as polymerization.6-8 Numerous theoretical studies have examined the steadystate and dynamic characteristics of autocatalytic reactions carried out in an isothermal continuous-flow stirred tank reactor (CSTR). The works of Gray and Scott on quadratic autocatalysis, A + B f 2B,9 and cubic autocatalysis, A + 2B f 3B,10 with autocatalyst decay, B f C, show exotic bifurcation diagrams, with isolas and mushrooms in plots of conversion versus residence time. Although their works lay the ground rules for the investigation of autocatalytic reactions through a dynamic systems perspective, their investigation was confined to systems comprising a single autocatalytic reaction. On the multiple-species front, the work of Fredrickson and Stephanopoulos11 studied mixed microbial cultures in homogeneous environments with time-invariant input. They consider a pure and simple competition of two species, in the sense that they interact only through a shared resource. They show that two microbial populations cannot coexist in a steady state under these conditions11,12 if neither of them is fed to the system. However, the coexistence of such competing microbial species is possible if the homogeneity of the environment condition is relaxed (e.g., by wall attachment13). In a recent study on the pure and simple competition of multiple autocatalysts, Birol and Teymour14 found that a stable steady-state coexistence is also not possible for this system if its input is timeinvariant and free of autocatalysts. Later, Birol et al.15 showed that, under these homogeneous conditions, dynamic coexistence was not possible either. * To whom correspondence should be addressed. Tel.: +1312-567-8947. Fax: +1-312-567-8874. E-mail: [email protected].

To ensure the long-term survival of a species in such an environment, we consider the case where small concentrations of that species are always present in the feed. Such a feed of species is conceivable if the investigated system is part of a reactor network and the fed species has a steady-state concentration in the feeding reactor. Similarly, pulses of invasion can result from bursts of temporary connectivity between the investigated reactor and a reactor where the invading species has a steady-state concentration. Thus, whenever the conditions in the environment turn favorable, the fed species will have a chance to thrive in the system. This host autocatalytic species is considered to be robust, given the fact that it can never be totally washed out of the system. Chaivorapoj et al.16 studied the effect of a temporary, nonsustained disturbance of another autocatalytic species in the feed on a CSTR system with one robust species. Their analysis shows that, after the disturbance is removed, the second nonrobust species could either disappear from the environment or coexist in it with the robust species, depending on the reproduction and death rates of these species and the magnitude of the perturbation. They found that, whenever the nonrobust invading species gets a foothold in the environment, the concentration level of the host robust species is decreased considerably. To combat this type of an invasion, they devised openloop control strategies that eliminate the invading species from the environment and restore the original high-concentration steady state of the host species. Stemming from the above-mentioned strategies, the purpose of this work is to design feedback control schemes to eliminate the invading species and to restore the original concentration level of the host species in an isothermal CSTR. 2. Problem Statement 2.1. Autocatalytic Model. Consider an isothermal CSTR in which a resource R is utilized by two distinct cubic autocatalytic species P and Q for their reproduction and in which both species decay to an inert (dead) form with different rates. This reaction scheme is given by

10.1021/ie020427+ CCC: $25.00 © 2003 American Chemical Society Published on Web 07/10/2003

3766

Ind. Eng. Chem. Res., Vol. 42, No. 16, 2003 krp

R + 2P 98 3P kdp

P 98 Dp krq

R + 2Q 98 3Q kdq

Q 98 Dq

(2) (3) (4) (5)

where Dp and Dq are the dead forms of species P and Q and krp (krq) and kdp (kdq) are the kinetic rate constants for the reproduction and death reactions of species P (Q), respectively. The only interaction between species P and Q is through the utilization of the same resource R, which is considered to be pure and simple competition, as defined by Fredrickson and Stephanopoulos.11 When the concentrations of resource and the two species in the CSTR are represented by capital letters and their concentrations in the feed by capital letters with subscript o, the mass balance equations can be written as

1 dR ) -krpRP2 - krqRQ2 + (Ro - R) dt′ θ

(6)

dP 1 ) krpRP2 - kdpP + (Po - P) dt′ θ

(7)

dQ 1 ) krqRQ2 - kdqQ + (Qo - Q) dt′ θ

(8)

where θ is the reactor residence time. In general, the resource R, with concentration Ro, and a trace amount of species P, with concentration Po, are fed continuously to the reactor, whereas species Q is introduced to the reactor as a transient pulse disturbance in the feed. Thus, the concentration Qo, is set to zero in the following text, unless stated otherwise. The model equations can be simplified by introducing the dimensionless concentrations r ) R/Ro, p ) P/Ro, po ) Po/Ro, and q ) Q/Ro; dilated time t ) krpRo2t′; scaled parameters dp ) kdp/ krpRo2, dq ) kdq/krpRo2, and kq ) krq/krp; and dilated residence time T ) krpθRo2. We will assume, without loss of generality, a unit reproduction rate of species P, krp ) 1, because all of the following analyses and conclusions stay valid in the general setup. Therefore, when krp ) 1 and Qo ) 0 are selected, the dimensionless mass balance equations reduce to

1 dr ) -rp2 - kqrq2 + (1 - r) dt T

(9)

dp 1 ) rp2 - dpp + (po - p) dt T

(10)

dq 1 ) kqrq2 - dqq - q dt T

(11)

At steady state, where the derivative terms vanish, there will be two steady-state solutions for q in eq 11. The system will reduce to a single-species autocatalysis when q ) 0 and represents a two-species autocatalysis when q > 0. 2.1.1. Single-Species Autocatalysis. We will first consider the system with no invading species Q (q ) 0) in the reactor. Thus, the dimensionless mass balance equations reduce to the system investigated by Gray and Scott,10 namely,

dr 1 ) -rp2 + (1 - r) dt T

(12)


(13)

At steady state, it can be easily shown from eqs 12 and 13 that concentrations of r and p are related by

p)

1 - r + po 1 + dpT

(14)

After substitution of eq 14 into the steady-state version of eq 12, it can be shown that the concentration levels for r at steady state are given by the roots of the cubic polynomial

1 r3 - 2(1 + po)r2 + (1 + po)2 + (1 + dpT)2 r T 1 (1 + dpT)2 ) 0 (15) T

[

]

For any given values of po, T, and dp, there will be either one or three real solutions of r and p; thus, there will be one or three steady states. 2.1.2. Two-Species Autocatalysis. Next, let us consider the two-species autocatalysis (q > 0), where the second species is initially introduced to the system as a pulse disturbance in the feed. However, the effect of this temporary disturbance may not vanish after the pulse ends, owing to the ability of Q to replicate itself. The static structure of the system is thus altered by the presence of Q in the reactor, as demonstrated by bifurcation diagrams in the next section. If we let dq/dt ) 0 in eq 11, the nonzero steady-state concentrations q are related to r by

1 + dq T q) kq r

(16)

Substituting eq 16 into the steady-state version of eq 9 yields an expression for the concentration p as

p)

x

1 +d ) ( T q

2

kq r

2

+

(1 - r) rT

(17)

Note that the r terms in the denominators of eqs 16 and 17 will not impose any singularity, because the steadystate value of r cannot be physically equal to zero, as a constant supply of r is fed to the system. By substitution of eq 17 into the steady-state version of eq 10, the steady-state concentrations r are given by the roots of the fourth-order polynomial

T 2kq2r4 -2T 2kq2(po + 1)r3 + Tkq[kq(1 + dpT)2 + 2(1 + dqT)2 + Tkq(po + 1)2]r2 - Tkq[kq(1 + dpT)2 + 2(po + 1)(1 + dqT)2]r + (1 + dqT)2[kq(1 + dpT)2 + (1 + dqT)2] ) 0 (18)

Ind. Eng. Chem. Res., Vol. 42, No. 16, 2003 3767

Figure 1. Open-loop bifurcation diagrams for r, p, and q versus T planes for system parameters po ) 0.04, dp ) 0.03, dq ) 0.3840, and kq ) 64dq.

Now, for any given set of system parameters, kq, T, dp, dq, and po, there will be either zero, two, or four real roots of eq 18. Thus, this autocatalytic model can have a maximum of four steady states for the two-species case (q * 0) in addition to a maximum of three steady states resulting from the single-species case (q ) 0). The emergence of these additional steady states provides analytical evidence that the presence of the autocatalytic species Q can alter the static structure of the system. We will clarify this static change, by the use of bifurcation diagrams computed for a specific set of system parameters, in the following section. 2.2. Base Case. In this section, bifurcation diagrams and open-loop responses of the system are described briefly, to provide an insight into the steady state and dynamic structure of the uncontrolled system. A set of system parameters that provide complete static behavior of this system, i.e., seven steady states, is selected to be the base case for which feedback control strategies are devised. 2.2.1. Bifurcation Diagrams. As in the work of Gray and Scott,10 the residence time (T) is selected in this study as the main bifurcation parameter. Because of the difficulties of solving the fourth-order parametric polynomial of eq 18 analytically, CONTENT software17 is used to compute the bifurcation diagrams of the resource R, species P, and species Q concentrations versus T. For the base case, bifurcation diagrams of the system with single and two species are shown for the parameter set po ) 0.04, dp ) 0.03, dq ) 0.3840, and kq ) 24.576 (64dq) in Figure 1. Mushroom curves in (T, r) and (T, p) projections represent the single-species steady states, corresponding to a line at q ) 0 in the (T, q) projection. The solid curves in Figure 1 show the stable steady states and the dashed curves the unstable ones. The large isola in the (T, r) plane represents the Q-dominant steady states, which shows the regime in which the presence of the invading species Q suppresses the host species P in the reactor. The small isola, located inside the intersection set of the mushroom and the large isola in the (T, r) plane, represents coexistence steady states. The projection of this interaction isola on the (T, p) plane is located inside the mushroom, and on the (T, q) plane, it is located inside the big isola. Note that there are no stable steady states on this interaction

isola, which is similar to the variant of the problem with po ) 0.14 Because of the added consumption of the single resource R by species Q, the reproduction performance of species P in the coexistence case is decreased considerably, so that its isola in the (T, p) plane is trapped beneath the lower residence time branch of the mushroom and flattened near low concentration levels. As opposed to the robust nature of species P, species Q will not be able to survive in the reactor if the residence time of the system is not in the range for which its isola exists. 2.2.2. Open-Loop Dynamic Response. Next, we illustrate the dynamic response of the system to a disturbance, introduced as a rectangular concentration pulse of the invading species in the feed (qo), using Matlab18 simulations. Throughout this study, we fix the duration of the pulse in qo at 2 time units but allow its magnitude to change. The initial condition is selected at the upper stable steady state of the mushroom curve in p, for Ts ) 8.641 71, namely, rs ) 0.210 559, ps ) 0.658 678, and qs ) 0. This point is located in the multiple-steady-state region of the mushroom curve that coincides with the range of stable steady states of the isola; thus, the system has a total of seven steady states, of which three are stable. With qo ) 0.8 for t ∈ (50, 52), the open-loop response in Figure 2 shows that Q survives in the system with an asymptotic concentration of 0.208 812 while the concentration of P is reduced considerably. This simulation illustrates how the system migrates from a steady state on the mushroom of the host species to a steady state on the isola curve of the invading species at Ts ) 8.641 71. It should be noted that the diagrams for ro and T are included because these will reflect the control action in subsequent closed-loop plots. Also, the diagram for q is scaled to a maximum value of 1 × 10-6, which is assumed to be the detection limit capability of the system toward this species; thus, this diagram monitors the potential washout of q from the system. 2.3. Control Objectives. Having built an understanding for the system behavior, next we wish to alter it through feedback control, to achieve a desired process performance. Our control objective is to bring the system to an operating point with a high yield of the host species P. We take this operating point as an initial condition and disturb the system by a rectangular pulse of an invading species Q in the feed, qo. Figure 2 shows the open-loop response of the system for a pulse of amplitude qo ) 0.8 that lasts 2 time units. Therefore, we have to (1) eliminate the invading species Q from the reactor and (2) restore the initial desired operating point of the system. To accomplish this, we propose two types of nonlinear feedback control schemes and analyze them for their capability of achieving both control objectives. 3. Sequential Control 3.1. Control Strategy. In line with the two control objectives, our first approach is to achieve them sequentially through a primary controller that eliminates the invasion and a secondary controller, which is activated after the successful elimination of Q to restore the initial operating point. Because of their mathematical simplicity, we will consider both the primary and secondary controllers to be proportional feedback controllers. The first controller manipulates the residence time, T,

3768

Ind. Eng. Chem. Res., Vol. 42, No. 16, 2003

Figure 2. Open-loop response for Figure 1 with qo ) 0.8 for t ∈ (50, 52).

whenever the concentration of the invading species Q in the reactor is above a minimum detection limit (taken to be 1 × 10-6) to eliminate it. When the Q concentration is reduced below the minimum detection limit, the first controller assumes the elimination of the invasion and turns itself off. Figure 1 suggests that this elimination can be achieved by either increasing or decreasing the residence time, so that the system dynamics are shifted outside the stable region of the big isola curve in the bifurcation diagram; hence, Q will not be able to survive. The increasing (decreasing) T corresponds to the reducing (increasing) flow rate of the reactor feed, to starve Q to death (to wash it out). Therefore, to eliminate Q from the system, we can use either positive or negative feedback control for the first controller. Once Q is eliminated and the first controller is turned off, the second controller is activated to restore the initial high concentration of the host species because the system might be trapped at an undesired low-p stable steady state. The controlled and manipulated variables of the second controller are the concentration p of the host species P and the resource concentration ro in the feed, respectively. Unlike the first controller, the second controller can employ only a positive controller gain because we would like to increase the amount of resource available for P when its concentration is below the set point. In simulating the action of this control scheme, we consider physical constraints on the residence time and the resource concentration in the feed. Because a negative residence time is meaningless and a zero residence time corresponds to an infinite feed flow rate, we set the minimum achievable residence time to Tmin ) 0.01. Similarly, there is a maximum possible concentration of resource we can feed to the system, which we set to ro,max ) 2.5. The physical

implementation of this can be pictured as a reservoir of resource at ro,max being diluted to ro and fed to the reactor. 3.2. Control Laws. First of all, let us define a new scaled parameter: the dimensionless concentration of R in the feed ro ) Ro/Ros, where Ros is the initial steadystate feed concentration of R. Then, the differential equation for the resource, eq 9, is rearranged to

1 dr ) -rp2 - kqrq2 + (ro - r) dt T

(19)

The set point of Q is zero, and the set point of P is the initial steady-state concentration of species P (ps). Then, we can write the control laws on T and ro as

T ) Ts + kCq(0 - q)

(20)

ro ) ros + kCr(ps - p)

(21)

where Ts is the initial steady-state residence time and ros ) 1. kCq and kCr are controller gains, which are turned on and off according to the control strategies described in section 3.1. As we have discussed above, kCq could be either positive to decrease T for washing Q out or negative to increase T for starving Q to death, but kCr is always positive. 3.3. Tuning the Primary Controller. Next, we will analyze the bifurcation diagram of the closed-loop system with the primary controller turned on and the secondary turned off. On the basis of this analysis, we will tune the controller and justify the tuning through some case studies. 3.3.1. Closed-Loop Bifurcation Diagrams. Turning on the primary controller by setting kCq * 0 and turning off the secondary controller by setting kCr ) 0,


Figure 3. Closed-loop bifurcation diagram on r, p, and q versus kCq planes for system parameters po ) 0.04, dp ) 0.03, dq ) 0.3840, kq ) 64dq, and T ) 8.641 71 with nonzero q.

we can rewrite the model equations for the closed-loop system as

ros - r dr ) -rp2 - kqrq2 + dt Ts + kCq(-q)

(22)

po - p dp ) rp2 - dpp + dt Ts + kCq(-q)

(23)

q dq ) kqrq2 - dqq dt Ts + kCq(-q)

(24)

Note that, in the open-loop case, the residence time of the system T was used as the bifurcation parameter, which is now governed by the control law of eq 20. In that equation, Ts is the residence time at the set point and q is a state variable. Thus, the only free parameter is the controller gain kCq, which is used as the bifurcation parameter. Continuation of the steady states of the system with the host and the invading species of the previous example results in the bifurcation diagram of Figure 3. As before, the dashed and solid lines represent unstable and stable steady states, respectively. For kCq ) 0, the steady-state solutions of the closed-loop and open-loop systems are identical. Furthermore, we can use the control law in eq 20 to find the steady-state value of the residence time and convert the plots of r, p, and q versus kCq of Figure 3, into plots of r, p, and q versus T as shown in Figure 4. These plots are similar in appearance to those of Figure 1, but one should note that they differ markedly in their stability structure. This difference results from the additional relation in eq 20 that defines the closed-loop system. To make use of the fact that -kCq is the slope of the linear relation between q and T given in eq 20, let us replot the q versus T projection of the bifurcation diagram on a linear scale (Figure 5). The five straight

lines, corresponding to kCq ) -869.5, -46.3, 11.3, 12.8, and 179.6, connect the desired steady-state point, where q ) 0 and T ) Ts, to the bifurcation points where the steady-state structure of the closed-loop system changes. Using the controller gain kCq, we can classify the static structure of the system in four groups as in Table 1. Therefore, if we select a control gain kCq > 179.6 or kCq < -46.3, corresponding to the first two categories of Table 1, we are guaranteed to eliminate Q from the system. Note that, because the set point for q is zero, the positive and negative control gains can also be interpreted as negative and positive feedback, respectively. Also note that, multiple stable steady states can be found only for negative feedback control with 11.3 < kCq < 12.8. Moreover, if the system has a stable offset for a selected kCq, it can be estimated by the use of the bifurcation diagram of Figure 5. 3.3.2. Closed-Loop Dynamic Response. Now, let us consider some examples of the closed-loop response of the system when the primary controller is turned on to illustrate the dynamics of elimination of the invading species Q. In the following examples, the timing of the pulse in qo is fixed at t ∈ (50, 52) and its magnitude is varied. Various kCq values will be selected to illustrate different types of dynamic behavior, corresponding to different regions of the bifurcation diagram of Figure 5. For kCq ) 12, the static structure of the system is in the region of multiple stable steady states. If we introduce a pulse of magnitude qo ) 0.8, the closed-loop response of the system is as shown in Figure 6, where we observe that there is still some amount of Q left in the system, i.e., q ) 0.332 353 as t f ∞, as suggested by the bifurcation diagram. Note that this offset in q is higher than that of the open-loop response because its intersection point in Figure 5 is higher than that for kCq ) 0. Because kCq is positive, T is decreased to wash

3770


Figure 4. Closed-loop bifurcation diagrams for r, p, and q versus T in a logarithmic scale.

Figure 5. Closed-loop bifurcation diagram of q versus T in a linear scale for nonzero q. The straight lines represent the slope of -kCq with the origin at T ) 8.641 71 and q ) 0.

Q out, when its concentration is higher than the minimum detection limit at t = 51. T settles down to a constant value for t > 80 because the system is stabilized at a new steady state, yielding a constant offset in q. Because the secondary controller remains turned off, with kCr ) 0, there is no adjustment of ro; hence, it stays constant at ro ) 1.

For kCq ) -20 and the same pulse of the previous example (Figure 7), Q is removed by increasing T, which corresponds to starving Q to death. Even though, for kCq ) -20, the bifurcation diagram of Figure 5 shows a stable nonzero q, the controller successfully perturbs the system dynamics to such an extent that it gives the host species a survival edge. Another possible interpretation

Ind. Eng. Chem. Res., Vol. 42, No. 16, 2003 3771 Table 1. Classification of the Static Structure of the System with Nonzero q kCq range(s)

static structure for q * 0

(-∞, -869.5) (179.6, ∞)

no steady states

(-869.5, -46.3)

no stable steady states

(-46.3, 11.3) (12.8, 179.6)

single stable steady state

(11.3, 12.8)

multiple stable steady states

of this is that the size of the disturbance is not large enough to allow Q to have a foothold in the system. Note that the operating point of the system is located on the top part of the mushroom curve for p of Figure 1. Because increasing T moves the operating point of the system forward, after the invading species becomes extinct, the system steady state will still be on the highconcentration part of the mushroom, which will return to its initial value as T is decreased. In this scenario, the system will not be trapped at the lower steady state in p. Consequently, when the first controller is shut down (when q < 1 × 10-6), the initial condition is already restored without necessary action from the secondary controller. However, if the magnitude of the pulse in qo is increased from 0.8 to 18 (which creates severe conditions for P), the system is attracted by the stable nonzero q on the big isola; hence, the first controller with kCq ) -20 cannot get Q starved to death (Figure 8). Notice that, while this pulse magnitude may be argued to be too high to be realistic in most of the cases, the same effect can be achieved by a longer pulse duration of lower magnitude. Again, a constant offset in Q leads to a settling of T for t > 80. To enable the controller to eliminate Q, kCq is further decreased to -50, at which, according to Table 1, no stable coexistence steady states are found, and the closed-loop response in Figure 9 follows suit in restoring the initial state. For kCq ) 120 and a pulse of magnitude qo ) 0.8, the closed-loop response in Figure 10 shows an offset in q of about 0.064 085. Note that kCq ) 120 is located in the third category of Table 1, with a unique nonzero stable steady state for Q. In this simulation, T is decreased until it hits the minimum constraint at 0.01 (maximum flow rate) for t ∈ (50, 52). This offset can be eliminated if kCq is further increased to 181, where there exist no nonzero stable steady states for Q. The resulting closed-loop response is shown in Figure 11, where Q is successfully washed out of the reactor. Contrary to the case where negative controller gains are used, the operating point is shifted across the extinction point of the mushroom curve in p, driving P to the lower steadystate region. As a result, the system will be trapped below the ignition point at lower T of the mushroom curve in p; thus, increasing T back to its initial value will not restore the concentration of P to its operating point. This makes a case for the necessity of the secondary controller, which will be activated to resolve the problem after the primary controller successfully removes the invasion and is shut down, leaving the concentration p away from its operating point. 3.4. Tuning the Secondary Controller. Now, let us consider the second controller. Recall that the first controller will be shut down (kCq ) 0) when the dimensionless concentration of Q is less than 1 × 10-6, which is assumed to be close enough to zero. Then, the second controller will kick in only when the initial

condition of the system is not restored as a result of steady-state multiplicity. 3.4.1. Closed-Loop Bifurcation Diagrams. The possible steady states of the closed-loop system now fall on the mushroom curve of single-species steady states in Figure 1. Similar to the primary controller, we will use the bifurcation diagrams of the closed-loop system to tune the secondary controller gain, kCr, that ensures the restoration of the initial condition. The main free parameter of this closed-loop system is kCr, with T fixed at its initial value and ro determined by eq 21. By substitution of eq 21 into the model equations (11)-(19) with q ) 0, the closed-loop system is described by

dr 1 ) -rp + [ros + kCr(ps - p) - r] dt T

(25)


(26)

where ps is the set point, or the initial value of p, and ros ) 1. For kCr ) 0, steady states of p are those of the open-loop system given as 0.658 678, 0.121 144, and 0.046 065 for Ts ) 8.641 71. Recall that the upper steady-state solution in p is the set point (ps ) 0.658 678). The closed-loop bifurcation diagram with kCr as the main free parameter can be obtained as in Figure 12. Solid and dashed lines represent stable and unstable regions, respectively. The solid horizontal line in the (kCr, p) plane represents the behavior of the upper steady-state solution as expected for the set point. The middle dashed and lower solid curves in the (kCr, p) plane represent the behavior of the middle and lower steady-state solutions, as they merge at a limit point at kCr = 0.6. Thus, the set point becomes a unique stable solution for the closed-loop system when kCr > 0.6, allowing the secondary controller to restore the set point of the system, regardless of initial conditions. 3.4.2. Closed-Loop Dynamic Response. The closedloop response of Figure 11 showed that, as the invading species Q was eliminated from the system at t = 120, the final steady state was trapped at the lower steady state of p. Therefore, in the two-controller strategy, the primary controller is switched off after Q is eliminated, and the secondary controller is activated simultaneously in order to restore the set point. The following closedloop response examples show the dynamics of restoration attempts of the set point in t ∈ (120, 250), using values of kCr selected from different regions of the (kCr, p) plane of Figure 12. At kCr ) 0.1, where both of the upper and lower steady states in p are stable, the closedloop response in Figure 13 shows a considerable offset of 0.610 026 in p from the desired value ps ) 0.658 678 and an increase in ro from 1 to 1.06. However, when kCr is further increased to 0.8 beyond the limit point in Figure 12, the closed-loop response in Figure 14 shows that the secondary controller can restore P to its set point value by first increasing ro to 1.5 and then gradually decreasing it to its original value. 3.5. Controller Performance. From the above analyses, the implementation of this control scheme, for a system with the studied parameter set, satisfies both control objectives. Accomplishment of both control objectives, for this specific invading species with dq ) 0.3840 and kq ) 24.576 (64dq), is guaranteed by the implementation of the primary controller with either kCq > 179.6 or kCq < -46.3 and the secondary controller with kCr > 0.6.

3772


Figure 6. Closed-loop response for kCq ) 12 with qo ) 0.8 for t ∈ (50, 52).

Figure 7. Closed-loop response for kCq ) -20 with qo ) 0.8 for t ∈ (50, 52).

3.6. Continuation of Special Points. In the closedloop bifurcation diagram of Figure 5, the intersections of the main isola and the lines representing kCq ) 179.6, 12.8, 11.3, and -46.3 are points where the closed-loop system changes its stability structure. We have made use of this information to tune the controller. Despite the fact that this set of controller gains is specific to the host and the invading species investigated, we can exploit this information to study the extent of applicability of this control strategy when other types of

invading species are introduced to the system. To that end, we perform a two-parameter continuation on the Hopf bifurcation (HB) and limit (LP) points used previously for controller tuning. The two free parameters used are the reproduction and death rates of the invading species. Let us first consider the continuation of the LP at kCq ) 179.6 (Figure 15). The solid curve in Figure 15 represents the locus of kq and dq values for which a stability exchange occurs at kCq ) 179.6 in a T-q


Figure 8. Closed-loop response for kCq ) -20 with qo ) 18 for t ∈ (50, 52).

Figure 9. Closed-loop response for kCq ) -50 with qo ) 18 for t ∈ (50, 52).

bifurcation diagram of the closed-loop system. This curve divides the parameter space of Figure 15 into two regions. Any invading species with a kq and dq combination selected from the lower region will be successfully eliminated by this control strategy when kCq > 179.6, while successful elimination is not guaranteed for kq and dq combinations selected from the upper region of the diagram. The region below the solid curve represents a shift to the right and/or shrinkage of the big isola of

Figure 5;14 thus, a controller gain kCq > 179.6 would ensure successful elimination. Another implication of this two-parameter continuation diagram is that a primary controller with kCq > 179.6 will be successful in eliminating any invading species with kq/dq < 44.8 regardless of dq, which is defined by the region under the dashed tangent line below the LP continuation curve. Moreover, this plot also implies that the controller we considered can even successfully eliminate a

3774




certain family of infinitely stable invading species (dq ) 0) with kq < 10.6 because the solid curve intersects the kq axis around 10.6. Next, let us consider the continuation of the HB at kCq ) -46.3 (Figure 16), which again divides the parameter space into two regions. Similar to the previous case, the upper and lower regions define kq and dq values for unremovable and removable invading species, respectively. Note that, in this case, there is no intersec-

tion of the curve with the kq axis, which reasonably implies that infinitely stable invading species (dq ) 0) cannot be starved to death by increasing T. Also contrary to the previous case, the continuation plot does not suggest a threshold value for a reproduction to death rate ratio, for which successful elimination is guaranteed. However, previous analysis of the open-loop system shows that no species with kq/dq < 16 can have a steady-state presence in the environment.14


Figure 12. Closed-loop bifurcation diagrams of r and p versus kCr for system parameters po ) 0.04, dp ) 0.03, dq ) 0.3840, kq ) 64dq, and T ) 8.641 71.

Figure 13. Closed-loop response for kCq ) 181 and kCr ) 0.1 with qo ) 0.8 for t ∈ (50, 52).

4. Effect of the Detection Limit on the Controller Performance Throughout this study, the minimum detection limit in our simulations was fixed at 1 × 10-6, but in a realistic system such as a living organism, this limit may not be as low. This uncertainty in the detection

limit provides the motivation to investigate its effect on the performance of our control design. To that end, we fix the controller parameters kCq and kCr at 181 and 0.8, respectively. Recall that this set of controller gains, developed for the base case, achieves both control objectives. Now, we will attempt to develop an analytical

3776


Figure 14. Closed-loop response for kCq ) 181 and kCr ) 0.8 with qo ) 0.8 for t ∈ (50, 52).

Figure 15. Solid curve as the locus of LP at kCq ) 179.6 on a kq-dq continuation plot, which partitions the parameter space in terms of the success of the tuned controller. Dashed lines represent species with constant reproduction-to-death ratios indicated.

method to calculate the maximum possible detection limit that allows this controller to successfully remove a given invading species. Also, we will propose a remedy for a system with poor detection limit. 4.1. Maximum Detection Limit. From the model equations of the system, we know that, unless given as an initial condition, the concentration of the resource in the environment cannot exceed the feed level. Therefore, with r < 1 in eq 11, we can write the inequality

dq 1 < kqq2 - dqq - q dt T

(27)

Note that the right-hand side of the inequality (27) is negative for

q
kp/ dp). In terms of static structure of the system, while po is increased to eliminate Q, the mushroom curve of single-species steady states is translated to a single solution curve,10 resulting in the disappearance of steady-state multiplicity from the system. This guar-

antees the restoration of the initial steady state, after the elimination of Q. A possible way of manipulating po can be done by recovering P from the outflow and recycling it to the inflow. The only error signal for this controller is the difference between the dimensionless concentration of species P in the reactor (p) and its initial value or set point (ps). This difference is an indicator for both the presence of the second species in

3780


Figure 21. Bifurcation diagram on r, p, and q versus kCp projections of the closed-loop system with parameters po ) 0.04, dp ) 0.03, dq ) 0.3840, kq ) 64dq, and T ) 8.641 71 with q * 0.

the system and the deviation from the set point of p. Unlike the primary controller of the former scheme, this single proportional controller can only be implemented with a positive controller gain. The control law governing the feed level of po can be written as

po ) pos + kCp(ps - p)

(32)

where the positive quantities pos and kCp are the initial feed concentration of P and the controller gain, respectively. 5.2. Controller Tuning. We analyze the performance of the controller for different controller gains using bifurcation diagrams, with kCp being the free parameter. Then, we illustrate the closed-loop dynamic behavior of the system for various kCp values and compare the results with the two-stage controller developed in section 3. 5.2.1. Closed-Loop Bifurcation Diagrams. By substituting the control law of eq 32 into the process dynamics (eqs 11-19), we obtain the closed-loop system equations as

1 dr ) -rp2 - kqrq2 + (ro - r) dt T

(33)

dp 1 ) rp2 - dpp + [pos + kCp(ps - p) - p] (34) dt T dq 1 ) kqrq2 - dqq - q dt T

(35)

where pos ) 0.04 is the trace amount of P in the feed when the system is operating at steady state as in the previous section. Recall that the set point of P, ps, is the upper steady-state solution on the mushroom curve

in p for Ts ) 8.641 71, i.e., ps ) 0.658 678. For this residence time, the open-loop system has seven steady states, four of which are two-species steady states (q * 0) and three of which are single-species steady states (q ) 0). Continuation of these seven steady-state solutions results in the bifurcation diagrams of Figures 21 and 22 for two-species and single-species steady states, respectively. Note that the curves in the (kCp, r) and (kCp, p) planes of single-species steady states in Figure 22 correspond to a line at q ) 0 in the (kCp, q) plane. It is evident from the plot of q versus kCp of the two-species bifurcation diagram that the invading species can have a nonzero steady state in the system for a finite range of the controller gain. Similarly, investigating the plot of p versus kCp of a single-species bifurcation diagram, we observe that, beyond a certain controller gain, the set point level of the host species concentration is the only steady state, which is also always stable. Thus, when a controller gain that is greater than the higher of the two limiting controller gains, i.e., kCp = 1.3 for this example, is selected, the closed-loop system has a unique (stable) steady state, which is also the set point of the system. 5.2.2. Closed-Loop Dynamic Response. The following closed-loop responses illustrate the dynamic behavior of implementation of the single controller with kCp for different regions of the bifurcation diagrams of Figures 21 and 22. With the type of pulse we investigated in the open-loop response of Figure 2, i.e., qo ) 0.8 for t ∈ (50, 52), selecting kCp ) 0.5 gives the closedloop response of Figure 24, which shows an offset in q of about 0.25; hence, the manipulated variable, po, is not adjusted back to its initial value. If kCp is increased to 1, with the same type of the pulse qo, Q is eliminated from the reactor (Figure 25). Even though for kCp ) 1 the invading species has a stable nonzero steady state,


Figure 22. Bifurcation diagrams for r and p versus kCp projections of the closed-loop system with parameters po ) 0.04, dp ) 0.03, dq ) 0.3840, kq ) 64dq, and T ) 8.641 71 with q ) 0.

Figure 23. Bifurcation diagrams of (kCp, q) and (kCp, p).

the perturbation size is not large enough for Q to have a foothold in the system and is quickly outnumbered by P. Also, the set point is restored by the controller after the system reaches the steady state because the only stable steady state with q ) 0 is the set point for

kCp > 0.0145 (Figure 22). However, by increasing the perturbation size extremely to 18, we observe an offset of 0.17 in q in the closed-loop response of Figure 26 for kCp ) 1. This offset in q can be eliminated by increasing kCp beyond the limiting gain of 1.3, say, to 1.4 (as in

3782


Figure 24. Closed-loop response for kCp ) 0.5 with qo ) 0.8 for t ∈ (50, 52).

Figure 25. Closed-loop response for kCp ) 1 with qo ) 0.8 for t ∈ (50, 52).

Figure 27), for which the only stable steady state of the system is the set point. 5.3. Continuation of the LP. In the single feedback control scheme, the minimum successfully applicable controller gain is found at a LP in the bifurcation diagram of two-species steady states in Figure 23. This LP in Figure 23 occurs at kCp ) 1.3, beyond which only

the stable steady state of the set point exists. Continuation of the locus of this LP in the kq-dq plane will illustrate the family of removable and unremovable invading species for this control design, as shown in Figure 28. This curve divides the parameter space into two regions. Based on the closed-loop simulations with the same type of pulse in qo as above, the upper and lower


Figure 26. Closed-loop response for kCp ) 1 with qo ) 18 for t ∈ (50, 52).

Figure 27. Closed-loop response for kCp ) 1.4 with qo ) 18 for t ∈ (50, 52).

regions define kq and dq of unremovable and removable invading species for this controller with kCp > 1.3, respectively. Unlike previous cases, the line of constant kq/dq ) 64 has two intersections with the curve in Figure 28, which interestingly indicates that, in this family of invaders, only species with dq above or below certain values can be successfully eliminated. Notice that, however, this curve intercepts the kq axis, meaning that

invading species with kq < 1.32 and no decay (dq ) 0) can also be eliminated by this controller with kCp > 1.3. 6. Conclusions A CSTR populated with a host autocatalytic species P is disturbed with pulses of an invading species Q that utilizes the same resource R and consequently alters

3784


Figure 28. Solid curve as the locus of LP at kCp ) 1.3 for kq-dq continuation, and the dashed line represents species with kq/dq ) 64.

the steady state of the system. Control strategies designed to successfully reject this disturbance are proposed and tested. Two control objectives are defined requiring the elimination of the invading species and ensuring the return of the host species to its set point. One strategy involves a primary controller that controls the level of q by manipulating the reactor residence time and a secondary controller that restores the set point of p by manipulating the feed resource concentration. The two-controller approach is necessitated by the multiplicity of steady states present in the open loop, which might lead to the impossibility of achieving both control objectives with a single controller using the selected manipulated variables. Bifurcation analysis is utilized to tune the controllers and shows that, for proportional control, both positive and negative gains can be used and that large enough gains can totally eliminate offsets. Because of this characteristic, we did not consider the addition of an integral mode to the controller; however, if excessive signal noise is present, the use of proportional-integral control will be necessary. The positive and negative gains correspond respectively to increasing and decreasing of the residence time, and both are intended to provide enough destabilization of the steady states of existence of Q. Another strategy, using a single controller that controls p by manipulating the feed concentration po, was considered. Again bifurcation analysis showed that beyond a minimum gain all offset is eliminated. This strategy relies on the competitive exclusion principle,14 which in this case means that if P thrives in the system, it will consume resources at a level that would not sustain the survival of Q. This strategy has the advantage that it does not require a detection mechanism for Q but suffers from the requirement of feeding P, which is not always realistic. Both strategies rely on controlling p to restore the original state of the system, which under most circumstances could also be achieved by controlling the other state variable r, even when multiple steady states exist. In this case, however, such a strategy cannot ensure success because of the fact that new steady states emerge in the parameter space when Q is introduced.

This leads to the possibility that the same level of r can be achieved for a number of combinations of p and q and possibly other invading species. Variants of these control strategies, including proportional-integral and proportional-integral-derivative, will be tested in a networked CSTR environment in future studies. Literature Cited (1) Varetto, L. Studying Artificial Life with a Molecular Automaton. J. Theor. Biol. 1998, 193, 257. (2) Kauffman, S. Self-Replication: Even Peptides do it. Nature 1996, 382, 496. (3) Orgel, L. E. Molecular replication. Nature 1992, 358, 203. (4) Kauffman, S. At Home in the Universe; Oxford University Press: New York, 1995. (5) Sharov, A. A. Signs and Values. Proc. IEEE ISIC/CIRA/ ISAS Joint Conf. 1998, 677. (6) Wills, P. R.; Kauffman, S.; Stadler, B. M. R.; Stadler, P. F. Selection Dynamics in Autocatalytic Systems: Templates Replicating Through Binary Ligation. Bull. Math. Biol. 1998, 60, 1073. (7) da Silva, L.; Mundim, K. C.; Tsallis, C. Effects of cross-links on the autocatalytic polymerization of RNA-like chains. Physica A 1998, 259, 415. (8) Schuster, P.; Fontana, W. Chance and necessity in evolution: lessons from RNA. Physica D 1999, 133, 427. (9) Gray, P.; Scott, S. K. Autocatalytic Reactions in the Isothermal CSTR. Chem. Eng. Sci. 1983, 38, 29. (10) Gray, P.; Scott, S. K. Autocatalytic Reactions in the Isothermal CSTR. Chem. Eng. Sci. 1984, 39, 1087. (11) Fredrickson, A. G.; Stephanopoulos, G. Microbial Competition. Science 1981, 213, 972. (12) Hansen, S. R.; Hubbell, S. P. Single Nutrient Microbial Competition: Qualitatively Agreement Between Experimental and Theoretically Forecast Outcomes. Science 1980, 207, 1491. (13) Baltzis, B. C.; Fredrickson, A. G. Competition of Two Microbial Population for a Single Resource in a Chemostat When One of Them Exhibits Wall Attachment. Biotechnol. Bioeng. 1983, 25, 2419. (14) Birol, I˙ .; Teymour, F. Static and Dynamics of Multiple Cubic Autocatalytic Reactions. Physica D 2000, 144, 279. (15) Birol, I˙ .; Parulekar, S. J.; Teymour, F. Effect of Environment Partitioning on the Survival and Coexistence of Autocatalytic Replicators. Phys. Rev. E 2002, 66, 51916.

Ind. Eng. Chem. Res., Vol. 42, No. 16, 2003 3785 (16) Chaivorapoj, W.; Birol, I˙ .; Teymour, F. Competition Between Robust and Non-Robust Autocatalytic Replicators. Ind. Eng. Chem. Res. 2002, 41, 3630. (17) Kuznetsov, Yu. A.; Levitin, V. V. CONTENT: A Multiplatform Environment for Analyzing Dynamical Systems, Dynamical Systems Lab.; Centrum voor Wiskunde en Informatica: Amsterdam, The Netherlands, 1996.

(18) Matlab, version 5.3; The MathWorks Inc.: USA, 1999.

Received for review June 10, 2002 Revised manuscript received April 7, 2003 Accepted May 23, 2003 IE020427+

Feedback Control of a Continuous-Flow Stirred Tank Reactor with

Recommend Documents