It is well known that the biotechnology is one of the fields that over the last three decades has a very high and quick development. Therefore, due to their advantages, the control of industrial bioprocesses has been an important practical problem attracting wide attention. The main motivation in applying control methods to such living systems is to improve their operational stability and production efficiency. The operation in Stirred Tank Reactors (STR) has been and it is still a widely used technology in fermentation bioprocesses. But, other new technologies such as fixed bed, fluidized bed or air lift reactors, are considered for bioprocesses operation. These reactors present several advantages over the “classical” STRs. For instance, the fixed bed and fluidized bed reactors are characterized by higher production performance, i.e. larger production capacity and higher productivity (Bastin & Dochain, 1990; Bastin, 1991; Bouaziz & Dochain, 1993).
From mathematical point of view, the dynamics of these processes are characterized by partial differential equations and therefore are classified as distributed parameter systems (Bastin & Dochain, 1990; Bouaziz & Dochain, 1993; Christofides, 2001; Dochain et al., 1992). For instance, the concentrations of the reactants and products are not anymore homogeneous in the whole reactor, like in STRs, but are characterized by a spatial profile along the reactor. It is clear that the distributed parameter feature of these systems makes the control problem even more difficult (Bouaziz & Dochain, 1993; Christofides, 2001; Dochain et al., 1992; Petre & Selişteanu, 2007a; Slotine & Li, 1991). Therefore, the modelling and simulation of them, which are the objectives of this chapter, are associated with the formulation of the process model using partial differential equations (PDEs), in general nonlinear, and with the computation of the solution. Since in the most of cases for this kind of systems, with the exception of a few simple cases, there are no analytical methods for finding the solution of the involved equations, the following alternatives are commonly employed (Christofides, 2001; Petre & Selişteanu, 2007a; Slotine & Li, 1991; Vilas, 2008):
To assume that these processes behave like lumped parameter systems (the states are only time dependent).
To use classical numerical methods like finite differences, finite elements or finite volumes.
These methods are based on discretization techniques which allow us to approximate the infinite set of numbers that represent a continuous function by means of a finite set of parameters (Bouaziz & Dochain, 1993; Christofides, 2001; Dochain et al., 1992; Petre & Selişteanu, 2007a; Slotine & Li, 1991; Vilas, 2008).
The first option is only valid when the spatial distribution is negligible as compared with the time evolution, for instance in reactors where the homogenization of the medium is achieved by means of stirring devices (Dochain & Vanrolleghem, 2001; Vilas, 2008). Nevertheless, in the remaining cases it is necessary to use the second alternative. Its main inconvenience is that the numerical solution is computationally involved (especially in 2D or 3D spatial domains) making the approach unsuitable for real time tasks like control or online optimization (Bouaziz & Dochain, 1993; Christofides, 2001; Dochain et al., 1992; Petre & Selişteanu, 2007a; Slotine & Li, 1991; Vilas, 2008).
An alternative to these classical numerical methods is the development of some techniques for the projection of the PDEs onto a low dimensional subspace. In accordance to these techniques, the original PDEs are transformed into a set of ordinary differential equations (ODEs) known as reduced order model (Aksikas et al., 2007; Americano da Costa Filho et al., 2009; Bouaziz & Dochain, 1993; Christofides, 2001; Dochain et al., 1992; Hoo & Zheng, 2001; Petre et al., 2007; Shvarstman et al., 2000).
As a result, the first objective of this chapter is to provide the mathematical tools, which are used for most of numerical methods, for solving PDEs and, on this basis, to give a brief outline of the most commonly employed techniques. Among the different alternatives, some of them based on Garlerkin scheme will be described and used in this chapter. In particular, the finite element method will be chosen on the basis of its flexibility and reduced order models since they are the most efficient (Vilas, 2008). These reduced order models obtained in this way can be used either for the process simulation or computation of their solution.
As we mentioned above, over recent years, a considerable research effort was concentrated on the design of control policies for distributed process systems (Christofides, 2001). Standard approaches to control this kind of systems are based on the spatial discretization of the original set of PDEs to obtain a set of ODEs. This allows us to employ standard finite-dimensional methods just described above to construct the controller (Christofides, 2001; Dochain et al., 1992; Hoo & Zheng, 2001; Shvarstman et al., 2000). However, this approach can result in a set of ODEs of high dimensionality which could make the approach unsuitable for real time applications. Also, the controllability and observability properties would depend on the number of discretization points as well as its location and may lead to a poor control quality (Christofides, 2001). Due to these disadvantages, new methods based on spectral decomposition techniques, which take into account the spatially distributed nature of these systems, have developed (Aksikas et al., 2007; Shi et al., 2006). This approach uses the Galerkin method so as to approximate the system by a low-dimensional set of ODEs to design the controller (Aksikas et al., 2007; Shvarstman et al., 2000). In (Americano da Costa Filho et al., 2009) this approach is used in combination with the Lyapunov's direct method to derive stabilizing controllers applied in the case of chemical processes that are carried out in tubular reactors.
This chapter is an extended work of the research achieved in some works of the authors: (Petre, 2003; Petre & Selişteanu, 2005, 2007a, 2007b; Petre et al., 2007, 2008; Petre, 2008; Selişteanu & Petre, 2004), and deals with the approximation and simulations of the dynamical model for a class of nonlinear propagation bioprocesses.
First, the dynamics of a class of propagation bioprocesses involving n components and m reactions that are carried out in fixed bed reactors without dispersion is analyzed. Since the dynamics of these bioprocesses are described by partial differential equations, either for simulation but especially for their controlling, one method consists of approximation of these infinitely order models by finite order models. These approximate models are in fact a set of ordinary differential equations obtained here by orthogonal collocation method. More exactly, infinitely dimension of the initial parameter distributed model will be reduced by approximating the partial derivative equation of each reaction component by a finite number, equal to p+1, of ordinary differential equations at p+1 discrete spatial positions along the bioreactor. These points are chosen as zeros of some orthogonal polynomials. Since it is difficult to know the connections between the original distributed parameter model and its approximate version (Christofides, 2001), our objective is to analyze the behaviour of both models to observe their intrinsic dynamical properties. This is realized by simulations conducted in the case of a fixed bed bioreactor without diffusion (Dochain et al., 1992; Petre et al., 2007; Petre & Selişteanu, 2007a).
In the following the control problem of these classes of propagation bioprocesses is analyzed. Since the biotechnological processes have a nonlinear nature, to control these processes some nonlinear control techniques will be used. These techniques not only improve the linear control methods and allow the analysis of strong nonlinearities but it also allow us to deal with model uncertainties and even the controller design may result simpler than in its linear counterpart (Slotine & Li, 1991). A widely extended nonlinear control technique is the feedback linearization (Isidori, 1995; Khalil, 2002; Slotine & Li, 1991), which makes use of algebraic transformations to obtain a closed loop linear system in which the conventional control techniques can be applied. The main inconveniences of this technique are two: firstly, the tracking control problem may lead to complex transformations and secondly, model uncertainty may affect the control performance. Therefore, it is necessary to apply adaptive and robust-adaptive control techniques able to drive the system to the desired reference despite the presence of uncertainties.
Consequently, by using the obtained results in (Petre & Selişteanu, 2005, 2007a; Petre et al., 2007, 2008; Petre, 2008; Selişteanu & Petre, 2006), to control the mentioned propagation bioprocesses, in this chapter a class of nonlinear adaptive controllers are designed based on their finite order models. The nonlinear controller design is based on the input-output linearizing technique. The information required about the process is the measurements of the state variables and its relative degree. It must be noted that if for the analyzed process there are no accessible state variables, these will be estimated by using an appropriate state observer.
Numerical simulations conducted in the case of a fixed bed reactor are included to illustrate the performances of the presented adaptive control strategies.
All simulations are achieved by using the development, programming and simulation environment MATLAB (registered trademark of The MathWorks, Inc., USA).
The chapter is organized as follows. Section 2 introduces the distributed parameter dynamical model for the class of fixed bed reactors. Its reduction to an ordinary differential equation system by orthogonal collocation method is presented in Section 3. A detailed analysis of obtained results by application of this method in the case of a fixed bed reactor without diffusion is presented in Section 4. The adaptive control strategies of propagation bioreactors are developed in Section 5, the performances of the designed adaptive controllers being presented in Section 6. Finally, concluding remarks and further research directions are presented.
2. Dynamical model of fixed bed bioreactors
A fixed bed bioreactor is a reactor where the biomass is immobilized on fixed carriers such as polymers, porous glass or ceramics. Consider a fixed bed bioreactor without dispersion operating in plug flow conditions as shown in Fig. 1, in which takes place a single autocatalytic growth reaction , where is the specific growth rate, with one limiting substrate S and one biomass population X.
To achieve the model of this bioreactor consider a section of the reactor with length dz located at a distance z (0 ≤ z ≤ L) from the bioreactor input. Assuming that along the length L of the bioreactor the cross section is constant and equal to A, then, the volume of this section is dV = Adz.
The mass balance of substrate concentration S around this section is given by:
where the term is the time variation of the amount of substrate in the elementary volume dV, FS and are the influent and the effluent substrate into, respectively from the volume dV, where F is the hydraulic flow rate, and the term represents the amount of substrate consumed by the biomass in the volume dV, where k1 is the yield coefficient.
Similarly, the mass balance of the biomass concentration X in the volume dV is given by:
where the term is the time variation of the amount of biomass in the elementary volume dV and term represents the amount of biomass produced in the volume dV. Since the biomass is immobilized, the hydraulic terms FX and have to disappear in this equation.
The equations (3) constitute the distributed parameter dynamical model of the analyzed fixed bed bioreactor. For completeness, we must to define the limit and initial conditions as:
where Sin(t) is the influent substrate concentration and X0(z) is the initial immobilized biomass concentration.
Assume now that two reactions take place in bioreactor: (i) an autocatalytic growth reaction with one limiting substrate S and one biomass population X with a reaction rate , as in the previous case; (ii) a death reaction of microorganisms X → Xd, where Xd is the non-active biomass. If we assume that the non-active biomass leaves the bioreactor, the distributed parameter dynamical model of this fixed bed bioprocess will be described as:
where kd is the death coefficient. The limit and initial conditions are defined as:
Let us define the state vector with the following partitions:
If we denote by the reaction rate vector, the model (5) can be rewritten as:
with the limit and the initial conditions:
From the two above examples, one can deduce that in the case of a fixed bed bioreactor in which m biochemical reactions with n reactants take place, among which n1 are microorganisms fixed on some supports and which remain within the reactor, and n2 other components flow through the reactor, the distributed parameter dynamical model will be described as:
with the following limit conditions:
Usually in (11) and (12), is the fixed biomass concentration vector, is the other concentration vector, is the influent concentration vector, is the reaction rate vector, and are the yield coefficient matrices having appropriately structures and dimensions.
3. Approximation of the dynamical model via orthogonal collocation
Since the model (11) is infinitely dimensional, in this section we will reduce the model order by approximating it by a set of ordinary differential equations. From (11) one can see that the state variables and are functions of time and space, that is and . We shall reduce the dimension of the model (11) by approximating the partial derivative equation of each component of and by a finite number, equal to p+1, of ordinary differential equations at p+1 discrete spatial positions along the bioreactor. To do this, we expand each variable as a finite sum of products of some time functions and space functions as:
where , ; are the values of at some discrete spatial positions along the bioreactor, called collocation points, and the basis functions are chosen as orthogonal functions (e.g. Lagrange polynomials) such as:
The integer p in (13) corresponds to the number of interior collocation points determined by collocation method. The points z = z0 and z = zp+1 correspond to the input (z = 0) and the output (z = L) of the reactor.
It must be noted that the collocation method offers two important advantages: first, its implementation is easier and second, the nature and physical dimension of the state variables remain unchanged after the reduction procedure. Moreover, the orthogonal methods preserve mass balances (Christofides, 2001; Dochain et al., 1992).
According to collocation method, the partial derivative of with respect to z appearing in (11) can be written as:
By introducing (13)-(15) into (11) and (12), each partial derivative equation is transformed into p+1 differential equations at the p interior collocation points and at the output of the reactor. Thus it is obtained the following n(p+1) order system of one order ordinary differential equations:
4. Application to a fixed bed bioreactor
In this section, the above presented collocation method will be applied in the case of the distributed parameter dynamical model of the fixed bed bioreactor described by (5)-(6). We will consider four interior collocation points, i.e. p = 4. Firstly, we define the values of the concentration of X, S and Xd and of the specific growth rate at each interior collocation point and at the output of the reactor, zi, as: Xi = X(z = zi), Si = S(z = zi), Xdi = Xd(z = zi), .
For p = 4, the relations in (13) are particularized as follows:
Let us define the state vectors at collocation points as:
Numerous simulation experiments have been performed on the example presented above to illustrate the dynamical behaviour of the two classes of models (exactly and reduced). The values of bioreactor dimensions and of process parameters used in simulation are (Petre & Selişteanu, 2007a; Petre et al., 2008): L = 1 m, A = 0.02 m2, k1 = 0.4, kd = 0.05 h-1. For the specific growth rate we have chosen a Contois model:
with = 0.35 h-1 and KC = 0.4.
In fact, in simulations we compared the reduced order model obtained by orthogonal collocation method with another approximation of the original model (5) obtained by a finite difference method, where the derivatives of variable with respect to time and space are approximated by finite differences as:
where Δt and Δz are discretization intervals. The choice of the discretization intervals is a very important problem (Bouaziz & Dochain, 1993). In order to have a satisfactory accuracy, a high number of discretization points may be necessary, but this choice requires excessive computer time.
The graphics in Fig. 2 show the response of the bioprocess to a step of the influent substrate concentration Sin from 7.5 to 10 g/l at time t = 5 s, for three values of space discretization points (1: M = 100, 2: M = 200, 3: M = 400).
To choose the best approximation, for the simulation of the reduced model (16) different solutions have been tested. So, as orthogonal basis functions have been chosen Lagrange polynomials (14), which are assumed to be reliable and easy to compute. The collocation points have been chosen as the zeros of the Jacobi polynomials and of the Legendre polynomials. The Jacobi polynomials can be computed by using the following recursive expression (Corduneanu, 1981; Dochain et al., 1992):
The Legendre polynomials can be computed by using the following recursive expression (Corduneanu, 1981):
Different values of the parameters and of the Jacobi polynomials and different numbers of interior collocation points have been considered. We found out that the best choice of and is and . We also found out that a number p = 4 of interior collocation points are sufficient for correctly simulating the process.
The simulation results performed in the same conditions as in the first experiment are presented in Fig. 3. These graphics show the behaviour of reduced order process model for four different sets of values of and of the Jacobi polynomials of order 4 as follows:; ; ; .
Thus, we obtain a set of four Jacobi polynomials, whose zeros are given by the following four sets of values:
The initial simulation conditions have been chosen so as to correspond to a steady state obtained from:
Note that in all the simulations, the influent flow rate is constant, F = 2 l/s.
From Fig. 2 and 3 one can observe that the reduced order model performed by orthogonal collocation method using only p = 4 interior collocation positions obtained as solutions of the Jacobi polynomial with and constitutes the best approximate of initial model given by our simulations. This can be better observed from graphics in Fig. 5 where the comparative behaviour of the two classes of models is presented.
5. Adaptive control of the propagation bioprocesses
In this section, the control problem of a class of propagation bioprocesses that are carried out in fixed bed reactors without dispersion is presented. The nonlinear adaptive controllers are designed based on the finite order model (16) obtained from exactly model (11) by using the orthogonal collocation method. It can be see that the model (16) may be rewritten as (Bastin & Dochain, 1990; Petre, 2008):
where is the state vector, is the yield coefficient matrix, is the reaction rate vector, is the influent flow rate vector, and is the dilution matrix.
5.1. Problem statement
For the bioreactors described by the model (16) the control objective is to regulate the concentration of a single component at the bioreactor output, under the following conditions:
The control input is the influent flow rate .
The controlled variable is measured not only at the bioreactor output, but also at every interior collocation point and at the reactor input (only in the case of external substrate).
The yield coefficients are positive constants (some of them being unknown).
For simplicity, we will denote by y the concentration of the controlled component, by the value of y at each interior collocation point , , i.e. , and by the value of the controlled component at the output of bioreactor . Using these notations, may by expressed as a linear combination of state variables and as:
where is a vector with appropriately dimension, used to select the controlled variable.
Using (16), the dynamics of in (27) is given by:
Consider that for the bioprocess described by the model (20), the controlled variable is the substrate concentration at the output of the bioreactor, that is . Since the state vector is now given by
then the entries of the vector c in (27) will be: with , or .
The dynamics of the concentration is given by:
It is easy to verify that the term in (30) is a linear combination only of variables at the interior collocation points The term contains the influent concentrations at the input of the bioreactor. With the condition (iv), the last term in (30) can be rewritten as (Dochain et al., 1992):
where and contain the unknown and known reaction rates respectively, and and are given by:
As a conclusion, contains all the unknown parameters and contains the known reaction rates. Then, the dynamics of output takes the form:
5.2. Exactly linearizing controller
As it was mentioned above, the control objective is to regulate the concentration of variable at the output of the bioreactor at a desired value by acting on the feeding substrate flow rate F.
Controller design is achieved by using the input-output linearizing technique. Remember that the input-output linearizing principle (Isidori, 1995) consists in the calculus of a nonlinear control law such that the behaviour of closed loop system (controller + process) is the same as the behaviour of a linear stable system. Assume that for the closed loop system we wish to have the following first-order linear stable dynamics:
Firstly, we consider the ideal case, where maximum prior knowledge concerning the process is available. In particular we suppose that the parameters in (34) are known and all the state variables are available for on-line measurements. It can be seen that equation (34) has the relative degree equal to 1 (Isidori, 1995). Then, from (34) and (35), the above closed-loop dynamics will be achieved by implementing the following exactly linearizing nonlinear control law:
The control law (36) leads to the following linear error model:
with , which for has an asymptotic stable point at . But the use of control law (36) requires the complete knowledge of the process. It is well known that because of the reaction rates, the bioprocesses are characterized by highly nonlinear dynamics and furthermore the kinetic and process parameters are often partially or completely unknown; as a consequence an accurate model for these processes is difficult to develop. Therefore, in recent years, it has been noticed a great progress in adaptive and robust adaptive control of bioprocesses, due to their ability to compensate for parametric uncertainties (Bastin & Dochain, 1990; Petre & Selişteanu, 2005; Petre, 2008; Selişteanu & Petre, 2006). Consequently, in the following section, for the class of the presented fixed bed bioreactors, we will develop an adaptive control algorithm considering some real realistic conditions.
5.3. Adaptive control of fixed bed bioreactors
If the parameters in (36) are assumed unknown (see the conditions (iii) and (iv)), these will be replaced by their on-line estimates . Then the control law (36) becomes an adaptive control law given by:
The estimates can be on-line calculated by using, for example, a linear regressive parameter estimator (Bastin & Dochain, 1990; Petre & Selişteanu, 2005; Petre, 2008), described here by the following equations:
where stands for regressor matrix, is a positive and symmetric gain matrix, and , named forgetting coefficient, and are design parameters to control the stability and convergence properties of the estimator (Petre & Selişteanu, 2005; Selişteanu & Petre, 2006).
6. Simulation results
The performances of the designed nonlinear adaptive controllers were verified by several simulation experiments performed upon the fixed bed bioreactor described by the model (5). The values of bioreactor dimensions and process parameters used in simulation are the same as in Section 4 (Petre & Selişteanu, 2007a; Petre et al., 2008). Also, for the specific growth rate we have chosen a Contois model (21), with 0.35 h-1 and 0.4.
The interior collocation points of the reduced model (20) have been chosen as zeros of the Jacobi polynomials given in Section 4. For , and , the abscises of the four interior collocation points are: z1 = 0.3121, z2 = 0.5789, z3 = 0.8130, z4 = 0.9627. Of course, these values will determine the values of the entries in the matrices and .
The exactly linearizing control law (36) takes the form:
The behaviour of the closed loop system in the ideal case, when all the parameters are completely known, is presented in Fig. 6.
The initial simulation conditions correspond to a process steady state regime. So, for the interior collocation points, the used values are: X1(0) = 44.1051 mg/l, X2(0) = 19.8101 mg/l, X3(0) = 9.8169 mg/l, X4(0) = 6.2634 mg/l; X5(0) = 5.6010 mg/l; S1(0) = 2.9403 g/l, S2(0) = 1.3207 g/l, S3(0) = 0.6545 g/l, S4(0) = 0.4176 g/l, S5(0) = 0.3734 g/l, and X0(0) = 0 mg/l, S0(0) = Sin(0) = 7.5 g/l.
To verify the regulation properties of the controller, for the reference variable a piece-wise constant variation was considered as:
The value of the gain parameter in (41) is . The system evolves in open loop from the time to time , after which the system is closed by using the control law (41). The influent substrate concentration acts as a perturbation given by with g/l for s and g/l for s.
Assume now that the death parameter is known, and the specific growth rate is unknown. Assume also that in (41) can be rewritten as:
where is considered as an unknown positive parameter. It is clear that if should be known, then is a function of bioreactor state given by:
Assume also that at the output of the bioreactor the only measured variable is the substrate concentration S5. It can be seen that the practical implementation of the control law (41) requires the knowledge of the state , and of the specific reaction rate .
Since in (41) the variable is not directly measurable, this will be substituted by its estimate . For the estimation of , independent of the unknown specific reaction rate , we use an asymptotic state observer (Petre & Selişteanu, 2005), which can be derived as follows. Let us define the auxiliary state z as:
The dynamic of deduced from model (20) is expressed by the following linear stable equation:
Then, the adaptive version of the control law (41) is given by:
The adaptive algorithm given by (48), (39) and (49) was implemented under the same conditions as in the first case. The values of the controller design parameters used in simulations are: , , , and the initial conditions are: = 0.12, , = 2.75 g/l. Much more, in order to test the behaviour of the adaptive controlled system in more realistic circumstances, we considered that the measurements of controlled variable in all interior collocation points and also at the input and the output of the bioreactor are corrupted with an additive white noise with zero average (5% from their nominal values).
The simulation results are shown in Fig. 7. As in the first case, the system evolves in open loop starting from to time 10 s, after that the system is closed by using the above adaptive algorithm. The perturbation has the same evolution as in the ideal case.
From the graphics in Fig. 7 one can deduce that even if the initialization of and are different from their ideal values (given by and , the adaptive controller is efficiently both in regulation of controlled variable and in rejection of the perturbation despite the very high load variations of . The behaviour of controlled variables and of control inputs is comparable with the results obtained in the free noise simulation. One can observe also a good behaviour both of the proposed state observer (49), (45) and parameter estimator (39).
Moreover, it was proven that the adaptive algorithm given by (48), (39) and (49) is robust, that is even though the process model (5) has uncertainty parameters, the behaviour of closed loop system is good. It was verified that if the death coefficient suffers variations by comparison to its nominal value (e.g. h-1) the obtained results are still good.
The approximation of the infinitely order dynamical model for a class of nonlinear propagation bioprocesses described by partial differential equations was examined. These approximate models consist of a set of ordinary differential equations obtained by orthogonal collocation method. The results obtained by application of this method in the case of a fixed bed reactor showed that by an appropriately choosing of the collocation points along the reactor, the behaviour of the reduced order model is very close to the behaviour of original infinitely order model.
After that, the obtained reduced order model was used to design some control algorithms for these types of reactors. The controller design is based on the input-output linearization technique. The obtained algorithm was tested in the controlling problem of substrate concentration for a propagation bioprocess that is carried out in a fixed bed reactor.
The simulation obtained results demonstrated that the designed adaptive algorithms used in control of propagation bioreactors yield good results closely comparable to those obtained in the case when the process parameters are completely known and/or time invariable.
Moreover, these algorithms prove to be robust as well yielding good results even though the measurable variables are affected by noises and/or the model parameters suffer variations between wide limits. It must be also noted that these algorithms can relatively easily be extended to other types of distributed parameters bioreactors: fluidized bed, air lift reactors.
This work was supported by CNCSIS–UEFISCDI, Romania, project number PNII–IDEI 548/2008.