- Email: [email protected]
Nonlinear analysis and control of a continuous fermentation process
Computers and Chemical Engineering 26 (2002) 659– 670 www.elsevier.com/locate/compchemeng
Nonlinear analysis and control of a continuous fermentation process G. Szederke´nyi a,*, N.R. Kristensen b, K.M. Hangos a, S. Bay Jørgensen b a
Systems and Control Research Laboratory, Computer and Automation Institute HAS, P.O. Box 63, H-1518 Budapest, Hungary b Department of Chemical Engineering, Technical Uni6ersity of Denmark, DK-2800 Lyngby, Denmark
Abstract Different types of nonlinear controllers are designed and compared for a simple continuous bioreactor operating near optimal productivity. This operating point is located close to a fold bifurcation point. Nonlinear analysis of stability, controllability and zero dynamics is used to investigate open-loop system properties, to explore the possible control difficulties and to select the system output to be used in the control structure. A wide range of controllers are tested including pole placement and LQ controllers, feedback and input–output linearization controllers and a nonlinear controller based on direct passivation. The comparison is based on time-domain performance and on investigating the stability region, robustness and tuning possibilities of the controllers. Controllers using partial state feedback of the substrate concentration and not directly depending on the reaction rate are recommended for the simple fermenter. Passivity based controllers have been found to be globally stable, not very sensitive to the uncertainties in the reaction rate and controller parameter but they require full nonlinear state feedback. © 2002 Published by Elsevier Science Ltd. Keywords: Nonlinear controllers; Fermenter; Bioreactors
1. Introduction The analysis and control of nonlinear process systems is a challenging and emerging interdisciplinary field of major practical importance. The most common way to control nonlinear process systems is to use either linear techniques on locally linearized versions of the nonlinear models or model-based predictive control (Seborg, 1999). A number of powerful and theoretically well grounded tools and techniques have become available for nonlinear control design in the field of systems and control theory (Isidori, 1995; van der Schaft, 1996). Such techniques have been applied successfully in other application areas. However, most often these techniques do require symbolic computation and may be non-feasible for real process systems. This may be one of the reasons why they are neither well known nor extensively applied on process systems. Bioreactors in particular exhibit strong nonlinear characteristics and their operation is known to be difficult to reproduce and to control. Therefore, a rela* Corresponding author.
tively simple continuous bioreactor is selected serve as a benchmark problem for advanced nonlinear analysis and control techniques. Many authors have examined various approaches of analyzing and controlling fedbatch (Kuhlmann, Bogle, & Chalabi, 1998; Boskovic & Narendra, 1995; Johnson, 1987; Vanimpe & Bastin, 1995) and continuous cultivation processes (Takamatsu, Hashimoto, Shioya, Mizuhara, Koike, & Ohno, 1975; Kuhlmann, Bogle, & Chalabi, 1997; Te Braake, van Can, Scherpen, & Verbruggen, 1998). The purpose of the present paper is to investigate nonlinear bioreactor dynamics for design of different nonlinear stabilizing control techniques of a relatively simple continuous cultivation process. For this purpose, linear methods are used as a reference for the nonlinear methods. The paper is structured as follows. In Section 2 the model analysis, including stability, controllability and zero dynamics is presented. After selection of possible control structures the different control designs are developed and their design properties and performance are compared in Section 3. It is concluded for the investigated model type that controllers using partial state feedback of substrate concentration can be recom-
0098-1354/02/$ - see front matter © 2002 Published by Elsevier Science Ltd. PII: S 0 0 9 8 - 1 3 5 4 ( 0 1 ) 0 0 7 9 3 - 1
G. Szederke´nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
660
mended from a combined performance point of view.
design,
tuning
and
2. Model analysis Model analysis constitutes a worthwhile preample for control design of a process system. The dynamic properties of an open-loop system provide information on possible control problems and difficulties and may yield guidelines for control structure design. Below follows presentation of the continuous cultivation process model and designed operating point. Thereafter, the linear and nonlinear stability and controllability analysis tools are presented and applied on the simple process benchmark. Finally, the zero dynamics are analyzed to provide information for selection of outputs for control structures.
2.1. Bioreactor state-space model
v(S) =vmax
S . K2S 2 + S+K1
(3)
The first equation originates from the biomass component mass balance, while the second is from the substrate component mass balance. They are coupled by the nonlinear growth rate function v(S)X which is the main source of nonlinearity and uncertainty in this simple model. The variables and parameters of the model together with their units and parameter values are given in Table 1. The parameter values are taken from Kuhlmann et al. (1998). The above model can easily be written in standard input-affine form with the centered state vector x= [X( S( ]T = [X−X0 S−S0]T consisting of the centered biomass and substrate concentrations. The centered input flowrate is chosen as manipulate input variable u= F( = F− F0 x; = f(x)+ g(x)u,
(4)
Æ (X( + X0)F0 Ç v(S( + S0) (X( + X0)− Ã Ã V f(x)= Ã Ã, Ã − v(S( + S0) (X( + X0) + (SF − (S( + S0))F0 Ã È Y V É
In order to be able to focus on the key issues in controller design and performance analysis for bioreactors, a relatively simple bioreactor is selected. The continuous reactor is assumed perfectly stirred and an unstructured biomass growth rate model with substrate inhibition kinetics is chosen. Despite its simplicity, this model exhibits some of the key properties which render bioreactors difficult to operate, and therefore control design around the reactor is desirable.
with (X0, S0, F0) being a steady-state operating point.
2.1.1. Nonlinear state-space model An isothermal nonlinear continuous fermenter is considered in this paper with constant volume V and constant physico-chemical properties. The dynamics of the process is given by the state-space model
2.1.2. Calculation of the optimal operating point The maximal biomass productivity XF is selected as the desired optimal operating point, i.e. the substrate cost is assumed negligible. This equilibrium point can be calculated from the nonlinear model:
dx XF =v(S)X− , dt V
(1)
S0 =
dS v(S)X (SF −S)F =− + , dt Y V
(2)
(5)
1 −2K1 + 2 K 21 + S 2FK1K2 + SFK1 , 2 SFK2 + 1
X0 = (SF − S0)Y,
(6) (7)
and the corresponding inlet feed flow rate is F0 = v(S0)V.
where
(8)
Substituting the parameter values from Table 1 gives
Table 1 Variables and parameters of the fermentation process model X S F V SF Y vmax K1 K2
Æ (X( + X0) Ç Ã − Ã V g(x)= Ã Ã, Ã(SF − (S( + S0)) Ã È V É
Biomass concentration Substrate concentration Feed flow rate Volume Substrate feed concentration Yield coefficient Maximal growth rate Saturation parameter Inhibition parameter
4 10 0.5 1 0.03 0.5
[g/l] [g/l] [l/h] [l] [g/l] – [l/h] [g/l] [l/g]
g S0 = 0.2187 , l
g X0 = 4.8907 , l
l F0 = 3.2089 . h
(9)
2.1.3. Linearized model In order to compare linear and nonlinear model analysis and control techniques, the linearized version around the above steady-state point (X0, S0, F0) of the nonlinear state Eqs. (4) and (5) is determined: x; = Ax + Bu, where
(10)
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
661
Fig. 1. Open loop behavior of the system.
A=
n (f (x
= x=0
Æ 0 Ã Ã vmaxS0 Ã− 2 È (K2S 0 +S0 + K1)Y
vmaxX0(K2S 20 −K1) Ç (K2S 20 +S0 +K1)2 Ã Ã, vmaxX0(K2S 20 −K1) F0 Ã − VÉ (K2S 20 +S0 +K1)2 −
(11) Æ X0 Ç Ã −V Ã B= g(0)= Ã Ã. ÃSF −S0 Ã È V É
0 0.4011 A= − 1.6045 1.2033
n
2.2.1. Stability A local stability analysis shows that the system is stable, within a neighborhood of the desired operating point, but because the point is very close the folded bifurcation point (X*, S*, F*)=(4.8775, 0.2449, 3.2128),
(12)
The system matrices in the optimal operating point are the following
complement the usual analysis based on locally linearized models.
−1.227 B= 2.4453
n
(13)
2.2. Stability and controllability analysis The stability and controllability properties of the system play key roles in control and control structure design not only at the desired operating point but also within the entire foreseeable operating region. Therefore, nonlinear analysis techniques are recommended to
this stability region is very small. This is illustrated in Fig. 1, which shows that the system moves to the undesired wash-out steady state when it is started from close neighborhood of the desired operating point (X(0)= 4.7907(g/l), S(0)= 0.2187(g/l)).
2.2.1.1. Stability analysis based on local linearization. Stability analysis based on local linearization around the operating point depends upon the eigenvalues of the linearized state matrix A in Eqs. (10) and (11). These eigenvalues are a complex conjugated pair in our case: u12 = − 0.60179 0.5306i.
(14)
We can see that the process is indeed stable around the operating point but the linear analysis does not provide any information on the extent of the stability region.
662
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
2.2.1.2. Nonlinear stability analysis. Nonlinear stability analysis is based on Lyapunov technique which aims at finding a positive definite scalar generalized energy function V(x) which has negative definite time derivative within the whole operating region. Most often a general quadratic Lyapunov function candidate is used in the form of V(x) =x TQx, with Q being a positive definite symmetric quadratic matrix, which usually is diagonal. This function is scalar-valued and positive definite everywhere. The stability region of an autonomous nonlinear system is then determined by the negative definiteness of its timederivative: dV (V (V = x; = f( (x), dt (x (x where f( (x)=f(x) in the open loop case (assuming zero input) and f( (x)=f(x) + g(x) C(x) in the closed loop case where C(x) is the static linear or nonlinear feedback law. The diagonal weighting matrix Q in the quadratic Lyapunov function is selected in a heuristic way: a state variable which does not produce overshoots during the simulation experiments gets a larger weight than another state variable with overshooting behavior. In the new norm defined by this weighting, a more accurate estimate of the stability region can be obtained. With this analysis, we cannot calculate the exact stability region but the results provide valuable information for selecting the controller type and tuning its parameters. The nonlinear stability analysis results in then time derivative of the quadratic Lyapunov function as a
function of the state variables, which is a two variable function in our case seen in Fig. 2. The stability region of the open-loop system is the region on the (x1, x2) plane over which the function is negative.
2.2.2. Controllability The controllability analysis is performed using both linear and nonlinear techniques. 2.2.2.1. Analysis based on local linearization. After calculating the Kalman-controllability matrix of the linearized model we find that the system is controllable (in the linear sense) within a narrow neighborhood of the required operating point, because the controllability matrix has rank 2. However, even in an important operating point with a maximum reaction rate very close to it S* = K1/K2, X*=(Sf − S*)Y the controllability matrix has only rank 1, i.e. the system is not locally controllable there. Thus, a nonlinear controllability analysis is desirable to provide an understanding of the loss of linear controllability. 2.2.2.2. Nonlinear controllability analysis. Nonlinear controllability analysis based on the generation of controllability distributions is used for identifying the singular points of the state-space around which control of the system is problematic or even impossible. The local controllability distribution is generated incrementally in a algorithmic way (Isidori, 1995) in two steps as follows. The initial distribution is Z0 = span{g}.
(15)
This is extended by Lie-brackets of f and g in the first step Z1 = Z0 + span{[ f, g]}=span{g, [ f, g]},
(16)
(g (f f(x)− g(x). (x (x
(17)
[ f, g](x)=
The second step then gives the following Z2 = Z1 + [ f, Z1]+ [g, Z1] =span{g, [ f, g], [ f, [ f, g]], [g, [ f, g]]}.
(18)
The Lie-product of f and g is: [ f, g](X( , S( ) =
Fig. 2. Time derivative of the Lyapunov function of the open loop system with Q= I.
Æ vmax(X( +X0) (K2S( 2 +2K2S( S0 +K2S 20 −K1) (− SF +S( +S0) Ç Ã Ã (K2S( 2 +2K2S( S0 +K2S 20 +S( +S0 +K1)2YV Ã Ã. 2 2 Ã − vmax(X( +X0) (K2S( +2K2S( S0 +K2S 0 −K1) (− SF +S( +S0) Ã (K2S( 2 +2K2S( S0 +K2S 20 +S( +S0 +K1)2V È É
(19)
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
2.2.2.3. Singular points. At the point [X( S( ]T = [− X0 SF −S0]T (X =0[g/l], S= SF ) all the elements of Z2 (and, of course, all the elements of Z0 and Z1) are equal to zero. This means, that the controllability distribution has rank 0 at this point. Moreover, this singular point is a steady state point in the state-space. Herefrom follows, that if the system reaches this (undesired) point, it is impossible to drive the process away by manipulating the input feed flow rate. If X( = −X0 (X =0 [g/l]) and S( "SF −S0 (S " SF ), Z2 has rank 1. From a practical point of view it means that if the biomass concentration decreases to 0 [g/l] then it cannot be increased by changing the input flow rate. Physically this also makes sense since it is not possible to produce biomass if all the microorganisms have been washed out. Both of the above singular points are such ‘wash-out’ states. Therefore, these are undesirable states. It is easy to calculate from Eq. (19) that in addition to the previous results Z1 also loses rank at the points S( =
1 − 2K2S0 9 2 K2K1 , 2 K2
but if we calculate the Lie-products [ f, [ f, g]] and [g, [ f, g]], we find that these singular points ‘disappear’ but the previous ones do not.
2.2.2.4. Nonsingular points. At any other point in the state-space including the desired operating point [X( S( ]T =[0 0]T the controllability distribution has rank 2. Thus, the system is controllable in a neighborhood of these points and we can apply state feedback controllers to stabilize the process. 2.3. Analyzing the zero dynamics of the fermenter for output selection Zero dynamics of nonlinear systems indicates the stability properties of the closed system when the output is forced to be zero. Thus, it also contains essentially the same information as zeros of linear time-invariant systems. In order to analyze zero dynamics we need to extend the original nonlinear state Eq. (4) by a nonlinear output equation y =h(x),
(20)
where y is the output variable and h is a nonlinear function given by the sensor choice. Then the zero dynamics of an input-affine nonlinear system can be analyzed using a suitable nonlinear coordinate transformation (Isidori, 1995) z=b(x):
n n
z1 y = , u(x) z2
(21)
663
where u(x) is a solution of the following partial different equation (PDE): Lgu(x)= 0,
(22)
where Lgu(x)=((u/(x) g(x) i.e. (u (u g1 + g2 = 0. (x1 (x2
(23)
In the case of the simple fermenter model in Eqs. (4) and (5), we can analytically solve the above equation to obtain: u(x)=F
V( −Sf + x2 + S0) x1 + X0
(24)
where F is an arbitrary continuously differentiable function. Then we can use the simplest possible coordinate transformation z= b(x) in the following form:
n
Æ Ç y à à z1 = ÃV(− Sf + x2 + S0) Ã. z2 Ã Ã È É x1 + X0
(25)
2.3.1. Selecting the substrate concentration as output If a linear function of the substrate concentration is chosen as output, i.e. z1 = y=ksx2,
(26)
where ks is an arbitrary positive constant then the inverse transformation x= b − 1(z) is given by
n
Æ − z2X0ks − Sf Vks + Vz1 + S0Vks Ç Ã Ã x1 z2ks Ã. =Ã z1 x2 Ã Ã É È ks
(27)
Thus the zero dynamics in the transformed coordinates can be computed as z; 2 = u: =Lfu(x)=
(u x; , (x
(28)
which gives z; 2 = Lfu(x)+ Lgu(x)u= Lfu(x)=Lfu(b − 1(z)),
(29)
since by construction Lgu(x)= 0 (see Eq. (22)). The above equation is constrained by y= ksx2 =z1 =0. Then the zero dynamics of the system is given by the differential equation z; 2 = Lfu(b − 1(0, z2))= −
(z2Y+ V)S0vmax , Y(K2 + S 20 + S0 + K1)
(30)
which is linear and globally stable. The equilibrium state of the zero dynamics is at z2 = − V/Y which together with z1 = 0 corresponds to the desired equilibrium state x1 = 0, x2 = 0 in the original coordinates. The above
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
664
analysis shows that if we manage to stabilize the substrate concentration either by a full state feedback or by on output feedback (partial state feedback) or e6en by a dynamic controller (which is outside the scope of this paper) then the o6erall system will be stable.
2.3.2. Selecting the biomass concentration as output The output in this case is a linear function of the biomass concentration: z1 = y=kxx1.
(31)
The zero dynamics of the system is given by: z; 2 = −
Vvmax(z 22YX0 + z2(YVSf + VX0) + Sf V 2) , (K2X z + z2(2K2X0Sf V + VX0) + V 2(K2S 2f + Sf + K1))Y 2 2 0 2
(32) which is only locally stable around the desired equilibrium state and the right hand side of Eq. (32) has singular points (where the denominator is 0). The stability region is independent of kx and can be determined using the parameters of the system.
2.3.3. Selecting a linear combination of biomass and substrate concentrations as output In this case, the output is a linear combination of the biomass and substrate concentrations: z1 = y=kxx1 + ksx2.
(33)
In this case, the zero dynamics is also locally stable around the desired equilibrium state and it again has singular points. Furthermore, a new undesired equilibrium state appears at z2 =Vkx /ks which can be inside the operating region depending on the values of kx and ks. In summary: Analysis of the zero dynamics shows that the best choice of output to be controlled is the substrate concentration and in6ol6ing the biomass concentration into the output generally introduces singular points into the zero dynamics and narrows the stability region. These theoretical issues will help in understanding the following simulation results.
3. Controller design and performance analysis The following controllers are designed and compared on the relatively simple fermentation process: pole placement controller, LQ controllers, feedback linearization based controllers and finally controllers based on direct nonlinear passivation.
3.1. The control problem statement and comparison 6iewpoints The above model analysis shows the following special properties of the control problem:
1. The desired setpoint is close to a bifurcation point which is a local singular point from the viewpoint of controllability. 2. The system is stable only within a narrow neighborhood of the desired operating point. Two physical limitations are also important: It is most often difficult to measure the biomass concentration and modelling uncertainty is present in the growth rate term v of the model equations. Therefore, we prefer controllers of the following type.The controller uses only the substrate concentration for feedback, i.e. it should be a partial state feedback controller.The closed-loop system should be robust with respect to uncertainties in the reaction kinetic function v. The performance analysis of the controllers is performed by extensive simulation studies and by investigating the overall behavior of the controllers. The controllers are then compared using the following performance criteria: Stability region To obtain a crude estimate of the extent of the stability region a simple Lyapunov function in the positive definite form of V(x)= x TQx is defined and the time derivative of this function is computed and analyzed as a function of the state variables. Time-domain performance Here the qualitative behavior of the responses, the presence of overshoots and the possible spurious steady states are investigated. Tuning possibilities Besides availability of guidelines for tuning the controllers, robustness with respect to the controller parameters is also investigated.
3.2. Pole placement controller
The purpose of this section is to provide a simple controller design approach for later comparison to examine the possibilities of stabilizing the system by partial feedback (preferably by feeding back the substrate concentration only).
3.2.1. Pole placement by full state feedback First, a full state feedback is designed such that the poles of the linearized model of the closed loop system are at [− 1 − 1.5]T. The necessary full state feedback gain is Kpp = [− 0.3747
0.3429],
(34)
A simulation run is shown in Fig. 3 starting from the initial state X(0)= 0.1(g/l) and S(0)= 0.5(g/l). It is clearly seen that the closed loop nonlinear system has an additional undesirable stable equilibrium point and the controller drives and stabilizes the system towards this point. This stable undesired operating point can be
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
665
substrate concentration is the output, let us consider the following static partial state feedback u= Kx where K=[0 k] i.e. we only use the substrate concentration for feedback. The stability region of the closedloop system can be investigated using the time-derivative of the quadratic Lyapunov function. Fig. 5 shows that the stability region of the closed loop system is quite wide. Furthermore, it can be easily shown that for e.g. k= 1 the only stable equilibrium point of the closed loop system (except for the wash-out steady state) is the desired operating point. The eigenvalues of the closed loop system with k= 1 are − 0.9741 and − 2.6746. Fig. 3. Centered state variables and input, full state feedback pole placement controller X(0)= 0.1(g/l), S(0)=0.5(g/l).
3.3. LQ control LQ-controllers are popular and widely used for process systems. They are known to stabilize any stabilizable linear time invariant system globally, that is over the entire state-space. This type of controller is designed for the locally linearized model of the process and minimizes the cost function J(x(t), u(t))=
&
(x T(t) Rxx(t)+u T(t) Ruu(t)) dt,
0
(35) where Rx and Ru (the design parameters) are positive definite weighting matrices of appropriate dimensions. The optimal input that minimizes the above functional is in the form of a linear full state feedback controller
Fig. 4. Time derivative of the Lyapunov function as a function of centered state variables q1 = 1, q2 = 0.1, pole placement controller, Kpp = [− 0.3747 0.3429].
easily calculated from the state equations and the parameters of the closed loop system: X = 3.2152(g/l), S= 3.5696(g/l). The time derivative of the Lyapunov function is shown in Fig. 4. The appearance of the undesired stable operating point warns us not to apply linear controllers based on locally linear models for a nonlinear system without a careful prior in6estigation.
3.2.2. Partial linear feedback Motivated by the fact that the zero dynamics of the fermenter is globally asymptotically stable when the
Fig. 5. Time derivative of the Lyapunov function as a function of centered state variables q1 =1, q2 =0.1, partial linear feedback, k= 1.
666
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
3.3.2. Expensi6e control The weighting matrices in this case were Rx =10·I2 × 2 and Ru = 1. There were no significant differences in terms of controller performance compared to the previous case. The full state feedback gain in this case was K= [− 1.5635 2.5571]. The stability region is again investigated using the time-derivative of a quadratic Lyapunov function. The time-derivative function as a function of the centered state variables for cheap and expensive control is seen in Figs. 6 and 7, respectively. Unlike the linear case where LQR always stabilizes the system, it is seen that the stability region does not cover the entire operating region. Indeed, a simulation run in Fig. 8 starting with an ‘unfortunate’ initial state exhibits unstable behavior for the nonlinear fermenter. Note that an LQ controllers is structurally the same as a linear pole placement controller (i.e. a static linear full state feedback). Therefore, undesired stable steadystates may also appear depending on the LQ-design. Fig. 6. Time derivative of the Lyapunov function as a function of centered state variables q1 = 1, q2 = 1, LQ controller, cheap control.
3.4. Local asymptotic linearization
stabilization
6ia
feedback
A nonlinear technique, feedback linearization is investigated in order to change the system dynamics into a linear behavior. Then, different linear controllers can be employed on the feedback linearized system.
3.4.1. Exact linearization 6ia state feedback In order to satisfy the conditions of exact linearization, first we have to find an artificial output function u(x) that is a solution of the PDE (Isidori, 1995). Lgu(x)= 0.
(36)
Fig. 7. Time derivative of the Lyapunov function as a function of centered state variables q1 = 1, q2 = 1, LQ controller, expensive control.
u = −Kx. The results for two different weighting matrix selections are investigated.
3.3.1. Cheap control In this case, the design parameters Rx and Ru are selected to be Rx =I2 × 2 and Ru =1. The resulting full state feedback gain is K =[− 0.6549 0.5899].
Fig. 8. Unstable simulation run, LQ controller, expensive control, X(0) =0.1(g/l), S(0) = 0.5(g/l).
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
667
Moreover, the second new coordinate z2 depends on v which indicates that the coordinate transformation is sensitive to uncertainties in the reaction rate expression. Simulations showed that it is hard to numerically evaluate the functions h and i and the partially closed loop system produced infeasible large input signals considering the constrains on substrate flow rate. The system can be exactly linearized theoretically, but the feedback is hard to compute in practice. Moreover, in an engineering sense it is not practically useful to linearize such an output function as i.
3.4.2. I/O linearization Here we are looking for more simple and practically useful forms of linearizing the input–output behavior of the system. The static nonlinear full state feedback for achieving this goal is calculated as Fig. 9. Time derivative of the Lyapunov function as a centered state variables q1 =1, q2 =1, linearizing the biomass concentration, k = 0.5.
Note that the above equation are exactly the same as the Eq. (22) used for determining the coordinate transformation for analyzing zero dynamics of the system. Let us choose the simplest possible output function again i.e. u(x) =
V(−SF +x2 +S0) . x1 +X0
(37)
Then the components of state feedback u = h(x)+ i(x)6 for linearizing the system are calculated as h(x) =
− L 2f u(x) , LgLfu(x)
1 , i(x) = LgLfu(x)
Lf h(x) 1 + 6, Lgh(x) Lgh(x)
(42)
provided that the system has relative degree 1 in the neighborhood of the operating point where 6 denotes the new reference input. As we will see the key point in designing such controllers is the selection of the output (h) function where the original nonlinear state Eq. (4) is extended by a nonlinear output equation y=h(x) where y is the output variable.
3.4.2.1. Controlling the biomass concentration. In this case h(x)= X( = x1s and h(x)= −
Lf h(x) Vvmax(x2 + S0) −F0, = Lgh(x) K2(x2 + S0)2 + x2 + S0 + K1 (43)
i(x)= −
V . x1 + X0
(38) (39)
(44)
The outer loop for stabilizing the system is the following
and the new coordinates are z1 = u(x)=
u= h(x)+ i(x)6 = −
6= − kh(x).
V(−SF + x2 +S0) , x1 + X0
z2 = Lfu(x) = vmaxV(S0X0 + S0x1 + X0x2 + x1x2 −YSFS0 −YSFx2 +YS 20 +2Yx2S0 + Yx 22) . Y(x1 + X0)(K2x 22 +2K2x2S0 +K2S 20 +x2 +S0 +K1)
(45)
It was found that the stabilizing region of this controller is quite wide but not global. The time-derivative of the Lyapunov function is shown in Fig. 9.
The state-space model of the system in the new coordinates is
3.4.2.2. Controlling the substrate concentration. In this case the chosen output is h(x)= S( = x2. The full state feedback is composed of the functions
z; 1 = z2
(40)
h(x)=
z; 2 = 6
(41)
−
which is linear and controllable. The exactly linearized model may seem simple put if we have a look at the new coordinates we can se that they are quite complicated functions of x depending on both state variables.
Vvmax(x2 + S0) (x1 + X0) +F0, Y(K2(x2 + S0)2 + x2 + S0 + K1) (SF − x2 − S0) (46)
i(x)=
V . Sf − x2 − S0
(47)
668
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
In the outer loop a negative feedback with the gain k= 0.5 was applied i.e. 6= − 0.5x2. The time derivative of the Lyapunov function of the closed loop system as a function of X( and S( is shown in Fig. 10. Note that for this case it was proven (see the zero dynamics analysis) that the closed loop system is globally stable except for the singular points where the biomass concentration is zero.
3.4.2.3. Controlling the linear combination of the biomass and the substrate concentrations. In this case the output of the system was chosen as h(x)= Kx (48) where the row vector K is calculated as the result of the previously described LQR cheap control design problem. Then the function h and i are given as h(x)=
− Kf(x) Kg(x)
(49)
i(x)=
1 Kg(x)
(50)
The value of K was [ −0.6549 0.5899]T and in the outer loop a negative feedback with the gain k=0.5 was applied. The time derivative of the Lyapunov function of the closed loop system is shown in Fig. 11.
3.5. Passi6ity based control of the system In order to design a controller which stabilizes the system over the entire operating region we turn back to the quadratic Lyapunov functions used for closed-loop stability region analysis before: V(x)=x TQx
(51)
with a positive define diagonal weighting matrix Fig. 10. Time derivative of the Lyapunov function as a centered state variables q1 = 1, q2 = 1, linearizing the substrate concentration, k = 0.5.
Q=
q1 0
0 q2
n
(52)
The control strategy is first to render the system lossless with an inner state feedback. Therefore, the input is decomposed in the following way: u = 6p(x)+ 6,
(53)
where 6p is the inner feedback and 6 is the new external input. The state equation of the partially closed loop system with 6 =0 is given by x; =f(x)+ g(x)6p(x).
(54)
We can exactly calculate 6p that makes the partially closed loop system lossless with respect to the storage function V in Eq. (51):
Fig. 11. Time derivative of the Lyapunov function as a centered state variables q1 = 1, q2 =1, linearizing the linear combination of the biomass and substrate concentrations, K =[ − 0.6549 0.5899]T, k = 0.5.
d V= 0, dt
(55)
(V (V x; = ( f(x)+ g(x)6p(x))= 0, (x (x
(56)
6p(x)=
((V/(x)f(x) q x f (x)+ q2x2 f2(x) =− 1 1 1 , ((V/(x)g(x) q1x1g1(x)+ q2x2g2(x)
6p(0)= 0.
(57)
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
669
0.5 in Fig. 12. Note that in this case q1 and q2 are design parameters.
3.6. E6aluation and comparison of the controllers The evaluation and comparison of the controllers is performed using the evaluation criteria introduced in the beginning of this section. The results are summarized as follows:
Fig. 12. Time derivative of the Lyapunov function as a function of centered state variables q1 = 1, q2 = 1, passivity based controller, k =0.5.
If we define the output as y= LgV(x)=
(V g(x), (x
(58)
then the partially closed loop system becomes passive (more precisely, lossless) with the storage function V with respect to the supply rate s(6, y)=6y. (59) Therefore, the system can be stabilized with the outer feedback 6= −ky,
k\0
(60)
Using the special quadratic form of V(x) in Eq. (51) we obtain 6= −k
(V g(x) =kx TQg(x), (x
k \0.
(61)
In the case of the simple fermenter model the following simple quadratic feedback is obtained: 6 = − kx TQg(x)
= −k − q1x1 k \ 0.
(x1 +X0) (SF −(x2 +S0)) +q2x2 , V V (62)
The above direct passivation based controller stabilizes the fermenter globally, which is seen from the time-derivative of the quadratic Lyapunov function corresponding to the parameter values q1 =q2 =1, k=
3.6.1. Stability region Stability region has been investigated using the quadratic Lyapunov function. It has been found that the stability region is Narrow, hard to determine for full state feedback pole placement controller and for LQR (expensive control) controllers. Wide, hard to determine for partial linear controller and for LQR (cheap control) controller. Wide, can be estimated using the zero-dynamics analysis for linearization of the biomass concentration controller and for Linearization of the linear combination of the biomass and substrate concentration. Global for linearization of the substrate concentration and for passivity based control. 3.6.2. Time-domain performance Time-domain performance can be characterized for the investigated controllers as follows Acceptable with the possibility of having undesired stable steady for full state feedback pole placement and for LQR controllers. Excellent for the input–output linearization type controllers of every kind: linearization of the biomass concentration, linearization of the substrate concentration and linearization of the linear combination of the biomass and substrate concentration. Good for all the other controllers. 3.6.3. Tuning possibilities and robustness with respect to the uncertain model parameters (K1 and K2) Tuning possibilities and robustness with respect to the uncertain model parameters (K1 and K2) have been found to fall into the following categories The linear controllers including pole placement, partial linear state feedback and LQR are difficult to tune but their performance is not very sensitive with respect to the uncertain parameters if these are in a limited region (9 10%). All kinds of the input–output linearization type controllers have excellent tuning possibilities in time domain but their robustness with respect to the model parameters is very poor. The Lyapunov function can be set arbitrarily using passivity based controllers and it is robust with respect to the variations in the model parameters.
670
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
4. Conclusion Different types of nonlinear controllers are designed and compared for a simple continuous bioreactor operating near optimal productivity. This operating point is close to a fold bifurcation point in this paper. A relatively simple cultivation model with constant volume and physico-chemical properties and with a substrate inhibited biomass growth rate expression is used to form an input-affine nonlinear state space model. Nonlinear analysis of stability, controllability and zero dynamics provides information not only in the vicinity of the desired operating point but also within the whole operating region. Analysis of the zero dynamics shows that the best choice of output to be controlled generally narrows the stability region for the investigated model. Nonlinear stability analysis based on the Lynapunov technique is used to obtain a crude estimate of the stability region for the controllers. A wide range of controllers are tested including pole placement and LQ controllers, feedback and input – output linearizing controllers and a direct nonlinear passivation controller. The comparison is based on time-domain performance and on investigating the stability region, robustness and tuning possibilities of the controllers. Controllers using partial state feedback of substrate concentration and not directly depending on the reaction rate are recommended for the relatively simple bioreactor investigated. Such controllers have a wide and predictable operating region and use the measurable state variable only. Passivity based controllers have been found to be globally stable, rather insensitive to uncertainties in
the reaction rate, however they require full nonlinear state feedback.
Acknowledgements This research is partially supported by the Hungarian Research Fund through grant No. T032 479 which is gratefully acknowledged.
References Boskovic, J. D., & Narendra, K. S. (1995). Comparison of linear, nonlinear and neural-network-based adaptive controllers for a class of fed-batch fermentation processes. Automatica, 31, 817– 840. Isidori, A. (1995). Nonlinear Control Systems. Berlin: Springer. Johnson, A. (1987). The control of fed-batch fermentation processes —a survey. Automatica, 23, 691 – 705. Kuhlmann, C., Bogle, D., & Chalabi, Z. (1997). On the controllability of continuous fermentation processes. Bioprocess Engineering, 19, 53 – 59. Kuhlmann, Ch., Bogle, I. D. L., & Chalabi, Z. S. (1998). Robust operation of fed batch fermenters. Bioprocesses Engineering, 19, 53 – 59. Seborg, D.E., Ad6ances In Control, Highlights of ECC’99, chapter A, Perspective on Advanced Strategies for Process Control (Revisted). Springer, 1999. Takamatsu, T., Hashimoto, I., Shioya, S., Mizuhara, K., Koike, T., & Ohno, H. (1975). Theory and Practice of optimal control in continuous fermentation processes. Automatica, 11, 141–148. te Braake, H. A. B., van Can, E. J. L., Scherpen, J. M. A., & Verbruggen, H. B. (1998). Control of nonlinear chemical processes using neural models and feedback linearization. Computers and Chemical Engineering, 22, 1113 – 1127. van der Schaft, Arjan (1996). L2-Gain and Passi6ity Techniques in Nonlinear Control. Berlin: Springer. Vanimpe, J. F., & Bastin, G. (1995). Optimal adaptive control of fed-batch fermentation processes. Control Engineering Practice, 3, 939 – 954.
Nonlinear analysis and control of a continuous fermentation process G. Szederke´nyi a,*, N.R. Kristensen b, K.M. Hangos a, S. Bay Jørgensen b a
Systems and Control Research Laboratory, Computer and Automation Institute HAS, P.O. Box 63, H-1518 Budapest, Hungary b Department of Chemical Engineering, Technical Uni6ersity of Denmark, DK-2800 Lyngby, Denmark
Abstract Different types of nonlinear controllers are designed and compared for a simple continuous bioreactor operating near optimal productivity. This operating point is located close to a fold bifurcation point. Nonlinear analysis of stability, controllability and zero dynamics is used to investigate open-loop system properties, to explore the possible control difficulties and to select the system output to be used in the control structure. A wide range of controllers are tested including pole placement and LQ controllers, feedback and input–output linearization controllers and a nonlinear controller based on direct passivation. The comparison is based on time-domain performance and on investigating the stability region, robustness and tuning possibilities of the controllers. Controllers using partial state feedback of the substrate concentration and not directly depending on the reaction rate are recommended for the simple fermenter. Passivity based controllers have been found to be globally stable, not very sensitive to the uncertainties in the reaction rate and controller parameter but they require full nonlinear state feedback. © 2002 Published by Elsevier Science Ltd. Keywords: Nonlinear controllers; Fermenter; Bioreactors
1. Introduction The analysis and control of nonlinear process systems is a challenging and emerging interdisciplinary field of major practical importance. The most common way to control nonlinear process systems is to use either linear techniques on locally linearized versions of the nonlinear models or model-based predictive control (Seborg, 1999). A number of powerful and theoretically well grounded tools and techniques have become available for nonlinear control design in the field of systems and control theory (Isidori, 1995; van der Schaft, 1996). Such techniques have been applied successfully in other application areas. However, most often these techniques do require symbolic computation and may be non-feasible for real process systems. This may be one of the reasons why they are neither well known nor extensively applied on process systems. Bioreactors in particular exhibit strong nonlinear characteristics and their operation is known to be difficult to reproduce and to control. Therefore, a rela* Corresponding author.
tively simple continuous bioreactor is selected serve as a benchmark problem for advanced nonlinear analysis and control techniques. Many authors have examined various approaches of analyzing and controlling fedbatch (Kuhlmann, Bogle, & Chalabi, 1998; Boskovic & Narendra, 1995; Johnson, 1987; Vanimpe & Bastin, 1995) and continuous cultivation processes (Takamatsu, Hashimoto, Shioya, Mizuhara, Koike, & Ohno, 1975; Kuhlmann, Bogle, & Chalabi, 1997; Te Braake, van Can, Scherpen, & Verbruggen, 1998). The purpose of the present paper is to investigate nonlinear bioreactor dynamics for design of different nonlinear stabilizing control techniques of a relatively simple continuous cultivation process. For this purpose, linear methods are used as a reference for the nonlinear methods. The paper is structured as follows. In Section 2 the model analysis, including stability, controllability and zero dynamics is presented. After selection of possible control structures the different control designs are developed and their design properties and performance are compared in Section 3. It is concluded for the investigated model type that controllers using partial state feedback of substrate concentration can be recom-
0098-1354/02/$ - see front matter © 2002 Published by Elsevier Science Ltd. PII: S 0 0 9 8 - 1 3 5 4 ( 0 1 ) 0 0 7 9 3 - 1
G. Szederke´nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
660
mended from a combined performance point of view.
design,
tuning
and
2. Model analysis Model analysis constitutes a worthwhile preample for control design of a process system. The dynamic properties of an open-loop system provide information on possible control problems and difficulties and may yield guidelines for control structure design. Below follows presentation of the continuous cultivation process model and designed operating point. Thereafter, the linear and nonlinear stability and controllability analysis tools are presented and applied on the simple process benchmark. Finally, the zero dynamics are analyzed to provide information for selection of outputs for control structures.
2.1. Bioreactor state-space model
v(S) =vmax
S . K2S 2 + S+K1
(3)
The first equation originates from the biomass component mass balance, while the second is from the substrate component mass balance. They are coupled by the nonlinear growth rate function v(S)X which is the main source of nonlinearity and uncertainty in this simple model. The variables and parameters of the model together with their units and parameter values are given in Table 1. The parameter values are taken from Kuhlmann et al. (1998). The above model can easily be written in standard input-affine form with the centered state vector x= [X( S( ]T = [X−X0 S−S0]T consisting of the centered biomass and substrate concentrations. The centered input flowrate is chosen as manipulate input variable u= F( = F− F0 x; = f(x)+ g(x)u,
(4)
Æ (X( + X0)F0 Ç v(S( + S0) (X( + X0)− Ã Ã V f(x)= Ã Ã, Ã − v(S( + S0) (X( + X0) + (SF − (S( + S0))F0 Ã È Y V É
In order to be able to focus on the key issues in controller design and performance analysis for bioreactors, a relatively simple bioreactor is selected. The continuous reactor is assumed perfectly stirred and an unstructured biomass growth rate model with substrate inhibition kinetics is chosen. Despite its simplicity, this model exhibits some of the key properties which render bioreactors difficult to operate, and therefore control design around the reactor is desirable.
with (X0, S0, F0) being a steady-state operating point.
2.1.1. Nonlinear state-space model An isothermal nonlinear continuous fermenter is considered in this paper with constant volume V and constant physico-chemical properties. The dynamics of the process is given by the state-space model
2.1.2. Calculation of the optimal operating point The maximal biomass productivity XF is selected as the desired optimal operating point, i.e. the substrate cost is assumed negligible. This equilibrium point can be calculated from the nonlinear model:
dx XF =v(S)X− , dt V
(1)
S0 =
dS v(S)X (SF −S)F =− + , dt Y V
(2)
(5)
1 −2K1 + 2 K 21 + S 2FK1K2 + SFK1 , 2 SFK2 + 1
X0 = (SF − S0)Y,
(6) (7)
and the corresponding inlet feed flow rate is F0 = v(S0)V.
where
(8)
Substituting the parameter values from Table 1 gives
Table 1 Variables and parameters of the fermentation process model X S F V SF Y vmax K1 K2
Æ (X( + X0) Ç Ã − Ã V g(x)= Ã Ã, Ã(SF − (S( + S0)) Ã È V É
Biomass concentration Substrate concentration Feed flow rate Volume Substrate feed concentration Yield coefficient Maximal growth rate Saturation parameter Inhibition parameter
4 10 0.5 1 0.03 0.5
[g/l] [g/l] [l/h] [l] [g/l] – [l/h] [g/l] [l/g]
g S0 = 0.2187 , l
g X0 = 4.8907 , l
l F0 = 3.2089 . h
(9)
2.1.3. Linearized model In order to compare linear and nonlinear model analysis and control techniques, the linearized version around the above steady-state point (X0, S0, F0) of the nonlinear state Eqs. (4) and (5) is determined: x; = Ax + Bu, where
(10)
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
661
Fig. 1. Open loop behavior of the system.
A=
n (f (x
= x=0
Æ 0 Ã Ã vmaxS0 Ã− 2 È (K2S 0 +S0 + K1)Y
vmaxX0(K2S 20 −K1) Ç (K2S 20 +S0 +K1)2 Ã Ã, vmaxX0(K2S 20 −K1) F0 Ã − VÉ (K2S 20 +S0 +K1)2 −
(11) Æ X0 Ç Ã −V Ã B= g(0)= Ã Ã. ÃSF −S0 Ã È V É
0 0.4011 A= − 1.6045 1.2033
n
2.2.1. Stability A local stability analysis shows that the system is stable, within a neighborhood of the desired operating point, but because the point is very close the folded bifurcation point (X*, S*, F*)=(4.8775, 0.2449, 3.2128),
(12)
The system matrices in the optimal operating point are the following
complement the usual analysis based on locally linearized models.
−1.227 B= 2.4453
n
(13)
2.2. Stability and controllability analysis The stability and controllability properties of the system play key roles in control and control structure design not only at the desired operating point but also within the entire foreseeable operating region. Therefore, nonlinear analysis techniques are recommended to
this stability region is very small. This is illustrated in Fig. 1, which shows that the system moves to the undesired wash-out steady state when it is started from close neighborhood of the desired operating point (X(0)= 4.7907(g/l), S(0)= 0.2187(g/l)).
2.2.1.1. Stability analysis based on local linearization. Stability analysis based on local linearization around the operating point depends upon the eigenvalues of the linearized state matrix A in Eqs. (10) and (11). These eigenvalues are a complex conjugated pair in our case: u12 = − 0.60179 0.5306i.
(14)
We can see that the process is indeed stable around the operating point but the linear analysis does not provide any information on the extent of the stability region.
662
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
2.2.1.2. Nonlinear stability analysis. Nonlinear stability analysis is based on Lyapunov technique which aims at finding a positive definite scalar generalized energy function V(x) which has negative definite time derivative within the whole operating region. Most often a general quadratic Lyapunov function candidate is used in the form of V(x) =x TQx, with Q being a positive definite symmetric quadratic matrix, which usually is diagonal. This function is scalar-valued and positive definite everywhere. The stability region of an autonomous nonlinear system is then determined by the negative definiteness of its timederivative: dV (V (V = x; = f( (x), dt (x (x where f( (x)=f(x) in the open loop case (assuming zero input) and f( (x)=f(x) + g(x) C(x) in the closed loop case where C(x) is the static linear or nonlinear feedback law. The diagonal weighting matrix Q in the quadratic Lyapunov function is selected in a heuristic way: a state variable which does not produce overshoots during the simulation experiments gets a larger weight than another state variable with overshooting behavior. In the new norm defined by this weighting, a more accurate estimate of the stability region can be obtained. With this analysis, we cannot calculate the exact stability region but the results provide valuable information for selecting the controller type and tuning its parameters. The nonlinear stability analysis results in then time derivative of the quadratic Lyapunov function as a
function of the state variables, which is a two variable function in our case seen in Fig. 2. The stability region of the open-loop system is the region on the (x1, x2) plane over which the function is negative.
2.2.2. Controllability The controllability analysis is performed using both linear and nonlinear techniques. 2.2.2.1. Analysis based on local linearization. After calculating the Kalman-controllability matrix of the linearized model we find that the system is controllable (in the linear sense) within a narrow neighborhood of the required operating point, because the controllability matrix has rank 2. However, even in an important operating point with a maximum reaction rate very close to it S* = K1/K2, X*=(Sf − S*)Y the controllability matrix has only rank 1, i.e. the system is not locally controllable there. Thus, a nonlinear controllability analysis is desirable to provide an understanding of the loss of linear controllability. 2.2.2.2. Nonlinear controllability analysis. Nonlinear controllability analysis based on the generation of controllability distributions is used for identifying the singular points of the state-space around which control of the system is problematic or even impossible. The local controllability distribution is generated incrementally in a algorithmic way (Isidori, 1995) in two steps as follows. The initial distribution is Z0 = span{g}.
(15)
This is extended by Lie-brackets of f and g in the first step Z1 = Z0 + span{[ f, g]}=span{g, [ f, g]},
(16)
(g (f f(x)− g(x). (x (x
(17)
[ f, g](x)=
The second step then gives the following Z2 = Z1 + [ f, Z1]+ [g, Z1] =span{g, [ f, g], [ f, [ f, g]], [g, [ f, g]]}.
(18)
The Lie-product of f and g is: [ f, g](X( , S( ) =
Fig. 2. Time derivative of the Lyapunov function of the open loop system with Q= I.
Æ vmax(X( +X0) (K2S( 2 +2K2S( S0 +K2S 20 −K1) (− SF +S( +S0) Ç Ã Ã (K2S( 2 +2K2S( S0 +K2S 20 +S( +S0 +K1)2YV Ã Ã. 2 2 Ã − vmax(X( +X0) (K2S( +2K2S( S0 +K2S 0 −K1) (− SF +S( +S0) Ã (K2S( 2 +2K2S( S0 +K2S 20 +S( +S0 +K1)2V È É
(19)
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
2.2.2.3. Singular points. At the point [X( S( ]T = [− X0 SF −S0]T (X =0[g/l], S= SF ) all the elements of Z2 (and, of course, all the elements of Z0 and Z1) are equal to zero. This means, that the controllability distribution has rank 0 at this point. Moreover, this singular point is a steady state point in the state-space. Herefrom follows, that if the system reaches this (undesired) point, it is impossible to drive the process away by manipulating the input feed flow rate. If X( = −X0 (X =0 [g/l]) and S( "SF −S0 (S " SF ), Z2 has rank 1. From a practical point of view it means that if the biomass concentration decreases to 0 [g/l] then it cannot be increased by changing the input flow rate. Physically this also makes sense since it is not possible to produce biomass if all the microorganisms have been washed out. Both of the above singular points are such ‘wash-out’ states. Therefore, these are undesirable states. It is easy to calculate from Eq. (19) that in addition to the previous results Z1 also loses rank at the points S( =
1 − 2K2S0 9 2 K2K1 , 2 K2
but if we calculate the Lie-products [ f, [ f, g]] and [g, [ f, g]], we find that these singular points ‘disappear’ but the previous ones do not.
2.2.2.4. Nonsingular points. At any other point in the state-space including the desired operating point [X( S( ]T =[0 0]T the controllability distribution has rank 2. Thus, the system is controllable in a neighborhood of these points and we can apply state feedback controllers to stabilize the process. 2.3. Analyzing the zero dynamics of the fermenter for output selection Zero dynamics of nonlinear systems indicates the stability properties of the closed system when the output is forced to be zero. Thus, it also contains essentially the same information as zeros of linear time-invariant systems. In order to analyze zero dynamics we need to extend the original nonlinear state Eq. (4) by a nonlinear output equation y =h(x),
(20)
where y is the output variable and h is a nonlinear function given by the sensor choice. Then the zero dynamics of an input-affine nonlinear system can be analyzed using a suitable nonlinear coordinate transformation (Isidori, 1995) z=b(x):
n n
z1 y = , u(x) z2
(21)
663
where u(x) is a solution of the following partial different equation (PDE): Lgu(x)= 0,
(22)
where Lgu(x)=((u/(x) g(x) i.e. (u (u g1 + g2 = 0. (x1 (x2
(23)
In the case of the simple fermenter model in Eqs. (4) and (5), we can analytically solve the above equation to obtain: u(x)=F
V( −Sf + x2 + S0) x1 + X0
(24)
where F is an arbitrary continuously differentiable function. Then we can use the simplest possible coordinate transformation z= b(x) in the following form:
n
Æ Ç y à à z1 = ÃV(− Sf + x2 + S0) Ã. z2 Ã Ã È É x1 + X0
(25)
2.3.1. Selecting the substrate concentration as output If a linear function of the substrate concentration is chosen as output, i.e. z1 = y=ksx2,
(26)
where ks is an arbitrary positive constant then the inverse transformation x= b − 1(z) is given by
n
Æ − z2X0ks − Sf Vks + Vz1 + S0Vks Ç Ã Ã x1 z2ks Ã. =Ã z1 x2 Ã Ã É È ks
(27)
Thus the zero dynamics in the transformed coordinates can be computed as z; 2 = u: =Lfu(x)=
(u x; , (x
(28)
which gives z; 2 = Lfu(x)+ Lgu(x)u= Lfu(x)=Lfu(b − 1(z)),
(29)
since by construction Lgu(x)= 0 (see Eq. (22)). The above equation is constrained by y= ksx2 =z1 =0. Then the zero dynamics of the system is given by the differential equation z; 2 = Lfu(b − 1(0, z2))= −
(z2Y+ V)S0vmax , Y(K2 + S 20 + S0 + K1)
(30)
which is linear and globally stable. The equilibrium state of the zero dynamics is at z2 = − V/Y which together with z1 = 0 corresponds to the desired equilibrium state x1 = 0, x2 = 0 in the original coordinates. The above
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
664
analysis shows that if we manage to stabilize the substrate concentration either by a full state feedback or by on output feedback (partial state feedback) or e6en by a dynamic controller (which is outside the scope of this paper) then the o6erall system will be stable.
2.3.2. Selecting the biomass concentration as output The output in this case is a linear function of the biomass concentration: z1 = y=kxx1.
(31)
The zero dynamics of the system is given by: z; 2 = −
Vvmax(z 22YX0 + z2(YVSf + VX0) + Sf V 2) , (K2X z + z2(2K2X0Sf V + VX0) + V 2(K2S 2f + Sf + K1))Y 2 2 0 2
(32) which is only locally stable around the desired equilibrium state and the right hand side of Eq. (32) has singular points (where the denominator is 0). The stability region is independent of kx and can be determined using the parameters of the system.
2.3.3. Selecting a linear combination of biomass and substrate concentrations as output In this case, the output is a linear combination of the biomass and substrate concentrations: z1 = y=kxx1 + ksx2.
(33)
In this case, the zero dynamics is also locally stable around the desired equilibrium state and it again has singular points. Furthermore, a new undesired equilibrium state appears at z2 =Vkx /ks which can be inside the operating region depending on the values of kx and ks. In summary: Analysis of the zero dynamics shows that the best choice of output to be controlled is the substrate concentration and in6ol6ing the biomass concentration into the output generally introduces singular points into the zero dynamics and narrows the stability region. These theoretical issues will help in understanding the following simulation results.
3. Controller design and performance analysis The following controllers are designed and compared on the relatively simple fermentation process: pole placement controller, LQ controllers, feedback linearization based controllers and finally controllers based on direct nonlinear passivation.
3.1. The control problem statement and comparison 6iewpoints The above model analysis shows the following special properties of the control problem:
1. The desired setpoint is close to a bifurcation point which is a local singular point from the viewpoint of controllability. 2. The system is stable only within a narrow neighborhood of the desired operating point. Two physical limitations are also important: It is most often difficult to measure the biomass concentration and modelling uncertainty is present in the growth rate term v of the model equations. Therefore, we prefer controllers of the following type.The controller uses only the substrate concentration for feedback, i.e. it should be a partial state feedback controller.The closed-loop system should be robust with respect to uncertainties in the reaction kinetic function v. The performance analysis of the controllers is performed by extensive simulation studies and by investigating the overall behavior of the controllers. The controllers are then compared using the following performance criteria: Stability region To obtain a crude estimate of the extent of the stability region a simple Lyapunov function in the positive definite form of V(x)= x TQx is defined and the time derivative of this function is computed and analyzed as a function of the state variables. Time-domain performance Here the qualitative behavior of the responses, the presence of overshoots and the possible spurious steady states are investigated. Tuning possibilities Besides availability of guidelines for tuning the controllers, robustness with respect to the controller parameters is also investigated.
3.2. Pole placement controller
The purpose of this section is to provide a simple controller design approach for later comparison to examine the possibilities of stabilizing the system by partial feedback (preferably by feeding back the substrate concentration only).
3.2.1. Pole placement by full state feedback First, a full state feedback is designed such that the poles of the linearized model of the closed loop system are at [− 1 − 1.5]T. The necessary full state feedback gain is Kpp = [− 0.3747
0.3429],
(34)
A simulation run is shown in Fig. 3 starting from the initial state X(0)= 0.1(g/l) and S(0)= 0.5(g/l). It is clearly seen that the closed loop nonlinear system has an additional undesirable stable equilibrium point and the controller drives and stabilizes the system towards this point. This stable undesired operating point can be
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
665
substrate concentration is the output, let us consider the following static partial state feedback u= Kx where K=[0 k] i.e. we only use the substrate concentration for feedback. The stability region of the closedloop system can be investigated using the time-derivative of the quadratic Lyapunov function. Fig. 5 shows that the stability region of the closed loop system is quite wide. Furthermore, it can be easily shown that for e.g. k= 1 the only stable equilibrium point of the closed loop system (except for the wash-out steady state) is the desired operating point. The eigenvalues of the closed loop system with k= 1 are − 0.9741 and − 2.6746. Fig. 3. Centered state variables and input, full state feedback pole placement controller X(0)= 0.1(g/l), S(0)=0.5(g/l).
3.3. LQ control LQ-controllers are popular and widely used for process systems. They are known to stabilize any stabilizable linear time invariant system globally, that is over the entire state-space. This type of controller is designed for the locally linearized model of the process and minimizes the cost function J(x(t), u(t))=
&
(x T(t) Rxx(t)+u T(t) Ruu(t)) dt,
0
(35) where Rx and Ru (the design parameters) are positive definite weighting matrices of appropriate dimensions. The optimal input that minimizes the above functional is in the form of a linear full state feedback controller
Fig. 4. Time derivative of the Lyapunov function as a function of centered state variables q1 = 1, q2 = 0.1, pole placement controller, Kpp = [− 0.3747 0.3429].
easily calculated from the state equations and the parameters of the closed loop system: X = 3.2152(g/l), S= 3.5696(g/l). The time derivative of the Lyapunov function is shown in Fig. 4. The appearance of the undesired stable operating point warns us not to apply linear controllers based on locally linear models for a nonlinear system without a careful prior in6estigation.
3.2.2. Partial linear feedback Motivated by the fact that the zero dynamics of the fermenter is globally asymptotically stable when the
Fig. 5. Time derivative of the Lyapunov function as a function of centered state variables q1 =1, q2 =0.1, partial linear feedback, k= 1.
666
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
3.3.2. Expensi6e control The weighting matrices in this case were Rx =10·I2 × 2 and Ru = 1. There were no significant differences in terms of controller performance compared to the previous case. The full state feedback gain in this case was K= [− 1.5635 2.5571]. The stability region is again investigated using the time-derivative of a quadratic Lyapunov function. The time-derivative function as a function of the centered state variables for cheap and expensive control is seen in Figs. 6 and 7, respectively. Unlike the linear case where LQR always stabilizes the system, it is seen that the stability region does not cover the entire operating region. Indeed, a simulation run in Fig. 8 starting with an ‘unfortunate’ initial state exhibits unstable behavior for the nonlinear fermenter. Note that an LQ controllers is structurally the same as a linear pole placement controller (i.e. a static linear full state feedback). Therefore, undesired stable steadystates may also appear depending on the LQ-design. Fig. 6. Time derivative of the Lyapunov function as a function of centered state variables q1 = 1, q2 = 1, LQ controller, cheap control.
3.4. Local asymptotic linearization
stabilization
6ia
feedback
A nonlinear technique, feedback linearization is investigated in order to change the system dynamics into a linear behavior. Then, different linear controllers can be employed on the feedback linearized system.
3.4.1. Exact linearization 6ia state feedback In order to satisfy the conditions of exact linearization, first we have to find an artificial output function u(x) that is a solution of the PDE (Isidori, 1995). Lgu(x)= 0.
(36)
Fig. 7. Time derivative of the Lyapunov function as a function of centered state variables q1 = 1, q2 = 1, LQ controller, expensive control.
u = −Kx. The results for two different weighting matrix selections are investigated.
3.3.1. Cheap control In this case, the design parameters Rx and Ru are selected to be Rx =I2 × 2 and Ru =1. The resulting full state feedback gain is K =[− 0.6549 0.5899].
Fig. 8. Unstable simulation run, LQ controller, expensive control, X(0) =0.1(g/l), S(0) = 0.5(g/l).
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
667
Moreover, the second new coordinate z2 depends on v which indicates that the coordinate transformation is sensitive to uncertainties in the reaction rate expression. Simulations showed that it is hard to numerically evaluate the functions h and i and the partially closed loop system produced infeasible large input signals considering the constrains on substrate flow rate. The system can be exactly linearized theoretically, but the feedback is hard to compute in practice. Moreover, in an engineering sense it is not practically useful to linearize such an output function as i.
3.4.2. I/O linearization Here we are looking for more simple and practically useful forms of linearizing the input–output behavior of the system. The static nonlinear full state feedback for achieving this goal is calculated as Fig. 9. Time derivative of the Lyapunov function as a centered state variables q1 =1, q2 =1, linearizing the biomass concentration, k = 0.5.
Note that the above equation are exactly the same as the Eq. (22) used for determining the coordinate transformation for analyzing zero dynamics of the system. Let us choose the simplest possible output function again i.e. u(x) =
V(−SF +x2 +S0) . x1 +X0
(37)
Then the components of state feedback u = h(x)+ i(x)6 for linearizing the system are calculated as h(x) =
− L 2f u(x) , LgLfu(x)
1 , i(x) = LgLfu(x)
Lf h(x) 1 + 6, Lgh(x) Lgh(x)
(42)
provided that the system has relative degree 1 in the neighborhood of the operating point where 6 denotes the new reference input. As we will see the key point in designing such controllers is the selection of the output (h) function where the original nonlinear state Eq. (4) is extended by a nonlinear output equation y=h(x) where y is the output variable.
3.4.2.1. Controlling the biomass concentration. In this case h(x)= X( = x1s and h(x)= −
Lf h(x) Vvmax(x2 + S0) −F0, = Lgh(x) K2(x2 + S0)2 + x2 + S0 + K1 (43)
i(x)= −
V . x1 + X0
(38) (39)
(44)
The outer loop for stabilizing the system is the following
and the new coordinates are z1 = u(x)=
u= h(x)+ i(x)6 = −
6= − kh(x).
V(−SF + x2 +S0) , x1 + X0
z2 = Lfu(x) = vmaxV(S0X0 + S0x1 + X0x2 + x1x2 −YSFS0 −YSFx2 +YS 20 +2Yx2S0 + Yx 22) . Y(x1 + X0)(K2x 22 +2K2x2S0 +K2S 20 +x2 +S0 +K1)
(45)
It was found that the stabilizing region of this controller is quite wide but not global. The time-derivative of the Lyapunov function is shown in Fig. 9.
The state-space model of the system in the new coordinates is
3.4.2.2. Controlling the substrate concentration. In this case the chosen output is h(x)= S( = x2. The full state feedback is composed of the functions
z; 1 = z2
(40)
h(x)=
z; 2 = 6
(41)
−
which is linear and controllable. The exactly linearized model may seem simple put if we have a look at the new coordinates we can se that they are quite complicated functions of x depending on both state variables.
Vvmax(x2 + S0) (x1 + X0) +F0, Y(K2(x2 + S0)2 + x2 + S0 + K1) (SF − x2 − S0) (46)
i(x)=
V . Sf − x2 − S0
(47)
668
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
In the outer loop a negative feedback with the gain k= 0.5 was applied i.e. 6= − 0.5x2. The time derivative of the Lyapunov function of the closed loop system as a function of X( and S( is shown in Fig. 10. Note that for this case it was proven (see the zero dynamics analysis) that the closed loop system is globally stable except for the singular points where the biomass concentration is zero.
3.4.2.3. Controlling the linear combination of the biomass and the substrate concentrations. In this case the output of the system was chosen as h(x)= Kx (48) where the row vector K is calculated as the result of the previously described LQR cheap control design problem. Then the function h and i are given as h(x)=
− Kf(x) Kg(x)
(49)
i(x)=
1 Kg(x)
(50)
The value of K was [ −0.6549 0.5899]T and in the outer loop a negative feedback with the gain k=0.5 was applied. The time derivative of the Lyapunov function of the closed loop system is shown in Fig. 11.
3.5. Passi6ity based control of the system In order to design a controller which stabilizes the system over the entire operating region we turn back to the quadratic Lyapunov functions used for closed-loop stability region analysis before: V(x)=x TQx
(51)
with a positive define diagonal weighting matrix Fig. 10. Time derivative of the Lyapunov function as a centered state variables q1 = 1, q2 = 1, linearizing the substrate concentration, k = 0.5.
Q=
q1 0
0 q2
n
(52)
The control strategy is first to render the system lossless with an inner state feedback. Therefore, the input is decomposed in the following way: u = 6p(x)+ 6,
(53)
where 6p is the inner feedback and 6 is the new external input. The state equation of the partially closed loop system with 6 =0 is given by x; =f(x)+ g(x)6p(x).
(54)
We can exactly calculate 6p that makes the partially closed loop system lossless with respect to the storage function V in Eq. (51):
Fig. 11. Time derivative of the Lyapunov function as a centered state variables q1 = 1, q2 =1, linearizing the linear combination of the biomass and substrate concentrations, K =[ − 0.6549 0.5899]T, k = 0.5.
d V= 0, dt
(55)
(V (V x; = ( f(x)+ g(x)6p(x))= 0, (x (x
(56)
6p(x)=
((V/(x)f(x) q x f (x)+ q2x2 f2(x) =− 1 1 1 , ((V/(x)g(x) q1x1g1(x)+ q2x2g2(x)
6p(0)= 0.
(57)
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
669
0.5 in Fig. 12. Note that in this case q1 and q2 are design parameters.
3.6. E6aluation and comparison of the controllers The evaluation and comparison of the controllers is performed using the evaluation criteria introduced in the beginning of this section. The results are summarized as follows:
Fig. 12. Time derivative of the Lyapunov function as a function of centered state variables q1 = 1, q2 = 1, passivity based controller, k =0.5.
If we define the output as y= LgV(x)=
(V g(x), (x
(58)
then the partially closed loop system becomes passive (more precisely, lossless) with the storage function V with respect to the supply rate s(6, y)=6y. (59) Therefore, the system can be stabilized with the outer feedback 6= −ky,
k\0
(60)
Using the special quadratic form of V(x) in Eq. (51) we obtain 6= −k
(V g(x) =kx TQg(x), (x
k \0.
(61)
In the case of the simple fermenter model the following simple quadratic feedback is obtained: 6 = − kx TQg(x)
= −k − q1x1 k \ 0.
(x1 +X0) (SF −(x2 +S0)) +q2x2 , V V (62)
The above direct passivation based controller stabilizes the fermenter globally, which is seen from the time-derivative of the quadratic Lyapunov function corresponding to the parameter values q1 =q2 =1, k=
3.6.1. Stability region Stability region has been investigated using the quadratic Lyapunov function. It has been found that the stability region is Narrow, hard to determine for full state feedback pole placement controller and for LQR (expensive control) controllers. Wide, hard to determine for partial linear controller and for LQR (cheap control) controller. Wide, can be estimated using the zero-dynamics analysis for linearization of the biomass concentration controller and for Linearization of the linear combination of the biomass and substrate concentration. Global for linearization of the substrate concentration and for passivity based control. 3.6.2. Time-domain performance Time-domain performance can be characterized for the investigated controllers as follows Acceptable with the possibility of having undesired stable steady for full state feedback pole placement and for LQR controllers. Excellent for the input–output linearization type controllers of every kind: linearization of the biomass concentration, linearization of the substrate concentration and linearization of the linear combination of the biomass and substrate concentration. Good for all the other controllers. 3.6.3. Tuning possibilities and robustness with respect to the uncertain model parameters (K1 and K2) Tuning possibilities and robustness with respect to the uncertain model parameters (K1 and K2) have been found to fall into the following categories The linear controllers including pole placement, partial linear state feedback and LQR are difficult to tune but their performance is not very sensitive with respect to the uncertain parameters if these are in a limited region (9 10%). All kinds of the input–output linearization type controllers have excellent tuning possibilities in time domain but their robustness with respect to the model parameters is very poor. The Lyapunov function can be set arbitrarily using passivity based controllers and it is robust with respect to the variations in the model parameters.
670
G. Szederke´ nyi et al. / Computers and Chemical Engineering 26 (2002) 659–670
4. Conclusion Different types of nonlinear controllers are designed and compared for a simple continuous bioreactor operating near optimal productivity. This operating point is close to a fold bifurcation point in this paper. A relatively simple cultivation model with constant volume and physico-chemical properties and with a substrate inhibited biomass growth rate expression is used to form an input-affine nonlinear state space model. Nonlinear analysis of stability, controllability and zero dynamics provides information not only in the vicinity of the desired operating point but also within the whole operating region. Analysis of the zero dynamics shows that the best choice of output to be controlled generally narrows the stability region for the investigated model. Nonlinear stability analysis based on the Lynapunov technique is used to obtain a crude estimate of the stability region for the controllers. A wide range of controllers are tested including pole placement and LQ controllers, feedback and input – output linearizing controllers and a direct nonlinear passivation controller. The comparison is based on time-domain performance and on investigating the stability region, robustness and tuning possibilities of the controllers. Controllers using partial state feedback of substrate concentration and not directly depending on the reaction rate are recommended for the relatively simple bioreactor investigated. Such controllers have a wide and predictable operating region and use the measurable state variable only. Passivity based controllers have been found to be globally stable, rather insensitive to uncertainties in
the reaction rate, however they require full nonlinear state feedback.
Acknowledgements This research is partially supported by the Hungarian Research Fund through grant No. T032 479 which is gratefully acknowledged.
References Boskovic, J. D., & Narendra, K. S. (1995). Comparison of linear, nonlinear and neural-network-based adaptive controllers for a class of fed-batch fermentation processes. Automatica, 31, 817– 840. Isidori, A. (1995). Nonlinear Control Systems. Berlin: Springer. Johnson, A. (1987). The control of fed-batch fermentation processes —a survey. Automatica, 23, 691 – 705. Kuhlmann, C., Bogle, D., & Chalabi, Z. (1997). On the controllability of continuous fermentation processes. Bioprocess Engineering, 19, 53 – 59. Kuhlmann, Ch., Bogle, I. D. L., & Chalabi, Z. S. (1998). Robust operation of fed batch fermenters. Bioprocesses Engineering, 19, 53 – 59. Seborg, D.E., Ad6ances In Control, Highlights of ECC’99, chapter A, Perspective on Advanced Strategies for Process Control (Revisted). Springer, 1999. Takamatsu, T., Hashimoto, I., Shioya, S., Mizuhara, K., Koike, T., & Ohno, H. (1975). Theory and Practice of optimal control in continuous fermentation processes. Automatica, 11, 141–148. te Braake, H. A. B., van Can, E. J. L., Scherpen, J. M. A., & Verbruggen, H. B. (1998). Control of nonlinear chemical processes using neural models and feedback linearization. Computers and Chemical Engineering, 22, 1113 – 1127. van der Schaft, Arjan (1996). L2-Gain and Passi6ity Techniques in Nonlinear Control. Berlin: Springer. Vanimpe, J. F., & Bastin, G. (1995). Optimal adaptive control of fed-batch fermentation processes. Control Engineering Practice, 3, 939 – 954.