Optimal Parametric Sensitivity Control for a Fed Batch Reactor. J.D. Stigter∗ and K.J. Keesman∗ Systems and Control Group Wageningen University and Research Centre Bomenweg 4 6703 HD Wageningen The Netherlands Fax: +31 317 484819 e-mail:
[email protected]
Keywords: Identification of Non-Linear Systems, Biotechnical Processes, Control and Optimization,
well known Fisher information matrix (FIM) Ztf Λ(tf ) =
Abstract The paper presents a method to derive an optimal parametric sensitivity controller for optimal estimation of a set of parameters in an experiment. The method is demonstrated for a fed batch bioreactor case study for optimal estimation of the saturation constant KS and, albeit intuitively, the parameter combination µmax Y .
1
Introduction
A model structure f may essentially be viewed as a mapping of the input signals (u), via the states (x) and parameters (θ), to the system outputs (y). In parametric models the output sensitivity with respect to a (set of) parameter(s) θ determines whether a parameter can be estimated from the input/output data. If the sensitivity of y with respect to θ is small or even zero, then the intstrumentation may not be well chosen or the input sequence u(t) may not excite the parametric sensitivities sufficiently. Since our aim is to reconstruct the true parameter values from input/output observations it is intuitively appealing to focus on what information of the parameter vector θ is observable in the output signal through ‘off-line’ or ‘on-line’ otsch and van den Hof calculation of the sensitivities dx dθ . D¨ (1996) have shown that, for a linear system model, identifiability can be interpreted in terms of controllability and observability of the linear system model Σ augmented with the system of associated sensitivities driven by the state x(t) and (possible) input u(t). It is also known that the gradient dy dθ provides access to a measure for the information content of the specific dataset. More specifically, the gradient appears in the
T
∂y(τ ) ∂y(τ ) Q−1 dτ ∂θ ∂θ
(1)
0
where Q is a weighting matrix, usually in diagonal form with each diagonal element inversely proportional to the variance of the associated measurement noise of sensor i (i = 1, . . . , m). The FIM provides a measure (after some criterion has been chosen ‘a priori’) of the information content of the data for the specific experiment conducted. A natural question is then, of course, how the input sequence u(t) can be chosen in such a way that the parameters can be optimally estimated. This is the well known problem of experimental design which is a classical problem in the identification literature, (Norton, 1986; Ljung, 1987).
1.1
Definitions
Let the model structure be defined by the dynamical equation x(t) ˙ = f (x(t), u(t), θ)
(2)
where x(t) is the n-dimensional state vector, u(t) the rdimensional input vector, θ the q-dimensional vector of model parameters, and f the, possibly non-linear, model structure that maps the inputs, via the states and model parameters, to the outputs of the system model. The question of identifiability of the parameter set θ can be stated more explicitely by the question ‘Given a model structure f : Is it possible to distinguish uniquely the values of the parameter vector θ on the basis of a carefully designed experiment? While this question may be answered by existing methods, e.g. the Taylor series expansion, (Pohjanpalo, 1978), our concern here is to find a feedback law that optimally controls the associated system of parametric sensitivities which can be derived straightforwardly from the model (2), leading to x˙ θ (t) ≈
∂f ∂f xθ (t) + ∂x ∂θ
(3)
Surprisingly enough the idea of finding such a feedback law on the basis of the model equation (2) augmented with the senitivity functions (3) has hardly been pursued in the identification literature.
2 2.1
Optimal Sensitivity Control A Fed-Batch Experiment
In the following we will derive a singular controller that maximally excites the sensitivity of the state with respect to the parameter θ (assuming a directly observable state x). In other words, we will maximally excite the sensitivity dy dθ to a specific parameter in the model structure f so that the gradient dy dθ is optimally present in the data set generated by the designed experimental setup. The thought experiment taken here as an example is a fed batch bioreactor experiment which can be used, for example, to determine the respiration rate of a population of bacteria feeding on a supplied substrate, (Vanrolleghem et al., 1995; Dochain et al., 1995). Observation of the substrate concentration can then be used to determine certain combinations of characteristic biokinetic parameters in terms of process yield Y in grams of biomass per gram of substrate, maximum specific growth rate µmax in grams of biomass per minute, and the saturation constant KS in grams per liter. Optimal input profiles (in terms of a Fisher design criterium) for the biokinetic parameters have been obtained for this setup in an ‘ad hoc’ manner Versyck (2000). To avoid a measurement problem of substrate (encountered in wastewater treatment) the study focusses on directly observable substrates such as sugar and ammonium. While observation of other substrates is difficult to achieve in practice this does not form a serious constraint for the method proposed here since the method can be extended to include, for example, oxygen and/or oxygen uptake rate (OUR) observations instead of substrate. Assume a bioreactor with biomass (X) that grows on the substrate, continuously fed into the bioractor dynamically with a dynamical feed rate u(t) in grams per liter. Further assume that the growth process does not contribute substantially to the biomass X over the timespan considered (which is in the order of magnitude of several minutes) so that the biomass X may be assumed constant over the interval [t0 , tf ] where tf marks the end of the experiment. The growth dynamics include Monod kinetics so that the dynamical model reads x˙ 1 (t) = −
x1 (t) µmax + u(t) Y Ks + x1 (t)
(4)
where x1 (t) is the substrate concentration in the bioreactor. The reactor is assumed completely mixed (a so called continuously stirred tank reactor or CSTR) and, as said, the consumption of the substrate x1 (t) is assumed to be
directly observable. The parameters considered in the foland the parameter KS lowing are the combination µmax Y which is sufficient. It is known that the parameter KS is most difficult to estimate from a set of observations (Holmberg and Ranta, 1982; Vanrolleghem et al., 1995) since this parameter is strongly correlated with the parameter µmax . dx1 ) is most profoundly observable Also the sensitivity ( dK S in the non-linear regime of the model, i.e. when the substrate concentration is of the same order of magnitude as 1 (t) the saturation constant KS . Define x2 (t) , dx dKS and 1 (t) x3 (t) , d(µdx . The sensitvities of the parameter KS max /Y ) are derived as and the parameter combination µmax Y
µmax (x1 (t) − KS x2 (t)) (5) Y (KS + x1 (t))2 µmax KS x1 (t) x˙ 3 (t) = x3 (t) + (6) Y (KS + x1 (t))2 KS + x1 (t)
x˙ 2 (t) =
In the sequel we will focus on finding a singular controller dx1 for optimal excitation of the sensitivity dK and we will S µmax assume that the combination Y is known for reasons of simplicity. It could be added that this assumption is in tune with the practical conditions since an obvious strategy to estimate the combination µmax is to saturate the reY actor completely with substrate so that the growth model is in the linear regime.
2.2
A Non-Linear Singular Control Problem
dx1 In order to optimally excite the sensitivity x2 (t) = dK S 2 we maximize x2 (t) and define the following Hamiltonian x˙ 1 (t) 2 H(x(t), λ(t)) , −x2 (t) + ( λ1 (t) λ2 (t) ) · (7) x˙ 2 (t)
where λ1 (t) and λ2 (t) are the co-states defined as ∂H µmax KS λ1 (t) λ˙ 1 (t) = − = + 2 ∂x1 Y (KS + x1 (t)) [x1 (t) − KS (1 + 2x2 (t))]λ2 (t) (KS + x1 (t))3 ∂H µmax KS λ2 (t) λ˙ 2 (t) = − = + 2x2 (t) ∂x2 Y (KS + x1 (t))2
(8) (9)
Since H does not explicitely depend on time a first integral of the problem is H = constant. Also, since the final time tf is assumed unknown and no terminal conditions are specified (determining the value of the co-states at tf ) this constant can be assumed equal to zero. Since the problem is linear in the control variable u(t) a singular control law that minimizes the Hamiltonian H over all possible input sequences u(t) can be derived by setting (A. E. Bryson (Jr.), 1999) ∀i ∈ N :
di dH =0 dti du
(10)
In order to determine u(t) explicitely only two differentiations are needed. For i = 0 we get λ1 (t) = 0 so that λ1 (t) is the switching function for this problem. Consistency d dH between the condition H = 0 and dt du = 0 also eliminates λ2 (t), so that from the case i = 1, the singular arc condition (or interior boundary condition) can be derived as x1 (t) = KS (1 + 2x2 (t))
(11)
Finally, the case i = 2 determines the optimal input u? (t) as u? (t) =
µmax Y
3
Results and Discussion
Since the parameter KS is not ‘a priori’ known it is important to investigate thesensitivity of (13)–(14) to the KS ? initial condition x (0) = . It was found through 0 numerical simulation that an increase or decrease of x1 (0) does not substantially change the trajectory of the sensitivity x2 (t), (see figure 1). The figure demonstrates that, indeed, the initial condition x1 (0) = KS does not introduce a major difference for the evolution of the sensitivity dx1 1 dKS for initial conditions x1 (0) = 2KS and x1 (0) = 2 KS . In a second numerical simulation exercise the system (4)–
(12)
This surprisingly simple ‘feedback law’, together with the singular condition (11), determines a trajectory in state space on which KS can be optimally identified. The condition (11) determines when to switch from a ‘bang’ input to the singular control law (12). The optimal control law u? (t) allows a reduction of the original set of equations (on a singular arc) giving KS µmax Y KS + x?1 (t) µmax (x?1 (t) − KS x?2 (t)) x˙ ?2 (t) = Y (KS + x?1 (t))2 x˙ ?1 (t) =
(13) (14)
where x?i denotes the optimal trajectory in state space on a singular arc, i.e. given dH du = 0. It is immediately clear that the optimal controlled system does not have a point of equilibrium and, therefore, is never in steady state. Since the initial condition x2 (0) , 0 it is easily deduced that x1 (0) = KS immediately puts the system in singular mode so that the optimal control law u? can be applied at the very beginning of the experiment. One can solve equation (13), given an initial condition x?1 (0) = KS , analytically giving r µmax ? x1 (t) = −KS + 2KS t + 4KS2 (15) Y ? x1 (t) ‘lives’ on a singular arc for Since the solution x?2 (t) which condition (11) holds an analytical solution for x?2 (t) follows immediately as r µmax ? x2 (t) = −1 + t+1 (16) 2KS Y Of course, the singular arc condition (11) depends on KS meaning that in order to determine the locus of the singular arc satisfying this condition, knowledge of the parameter KS is needed. Violation of constraint (11) will therefore be investigated in more detail for sensitivity to errors in an estimate of KS at time t0 . This, together with some other simulation results, will be presented in the next section.
Figure 1: Substrate and sensitivity KS = 2.0, and Y = 0.75.
dx1 (t) dKS
for µmax = 0.5,
(5) was simulated for three constant input values, namely u(t) = 21 u? , u(t) = u? , and u(t) = 32 u? . The resulting sensitivities are plotted in figure 2 from which it can clearly be observed that, indeed, the optimal input u? = µmax Y excites the information content for estimation of the parameter KS best. It could finally be noted that the optimal solution for optimal identification of KS may well be used in an identification experiment where both µmax and KS are to be Y estimated. Intuitively (but also through simulation) one can show that a initial injection of substrate (first ‘bang’) will excite the sensitivity x3 (t) substantially. The experiment could therefore be organized as follows: (i) Apply a ‘bang’, i.e. u(t) = umax , at the beginning of the experiment in order to identify the parameter combination µmax Y . (iii) Apply a second ‘bang’, i.e. u(t) = 0, after a short period of time ts and observe whether the singularity condition (11) holds. (ii) Once the singularity condition (11) is satisfied, switch to the control u? = µmax and estimate KS . Y The above algorithm was simulated for a switching time ts (from the first ‘bang’ to the second ‘bang’) of 15
parameter values with maximal sensitivity or minimum uncertainty. For the case of a more complicated sensitivity controller for which the analytical solution is not easy to derive one should consider a numerical scheme that optimizes the sensitivities on basis of a gradient dH du (A. E. Bryson (Jr.), 1999). Finally, it could be mentioned that the presented methodology may be developed in a linear model setting so that optimal estimation trajectories and/or controllers can be calculated leading to truly ‘on-line’ optimal experimental designs.
References Figure 2: Sensitivity x2 (t) for three different constant in3 µmax ? puts, namely u(t) = 12 µmax Y , u(t) = 2 Y , and u (t) = µmax Y .
A. E. Bryson (Jr.) (1999). Dynamic Optimization. Addison Wesley. Dochain, D., Vanrolleghem, P. A., and Daele, M. V. (1995). Structural Identifiability of Biokinetic Models of Activated Sludge Respiration. Water Research, 29(11):2571–2578. D¨otsch, H. G. M. and van den Hof, P. M. J. (1996). Test for Local Structural Identifiability of High-order Nonlinearly Parametrized State Space Models. Automatica, 32(6):875–883. Frieden, B. R. (1998). Physics from Fisher Information: A Unification. Cambridge University Press. Holmberg, A. and Ranta, J. (1982). Procedures for Parameter and State Estimation of Microbial Growth Process Models. Automatica, 18(2):181–193. Ljung, L. (1987). System Identification – Theory for the User. Prentice Hall.
Figure 3: A simulated identification experiment.
Norton, J. P. (1986). An Introduction to Identification. Academic Press.
minutes for which the results are presented in figure 3. First, a bang-input umax = 5 µmax was applied and afY ter ts = 15 min this controller was shut off. The resulting trajectory forms an excellent starting point for an optimal experiment in which the parameter combination µmax and Y KS are to be estimated.
Pohjanpalo, H. (1978). System Identifiability Based on the Power Series Expansion of the Solution. Mathematical Biosciences, 41:21–33.
4
Conclusions
A first attempt to include parametric sensitivities in the control loop has been made for a fed-batch reactor. The study has lead to a satisfactory controller for estimation of KS , assuming Monod kinetics for the substrate consumption. The simple control law u? (t) = µmax gives Y satisfactory results in the sense that it optimally excites the sensitivity of substrate with respect to the parameter KS . The simplicity of the controller allows simple implementation in an ‘on-line’ estimation schedule in which the estimator controls its own input in order to find the
Vanrolleghem, P. A., Daele, M. V., and Dochain, D. (1995). Practical identifiability of a biokinetic model of activated sludge respiration. Water Research, 29(11):2561–2570. Versyck, K. J. (2000). Dynamic Input Design for Optimal Estimation of Kinetic Parameters in Bioprocess Models. PhD thesis, Katholieke Universiteit Leuven.
