Computation of Optimal Identification Experiments for Nonlinear

Numerical solutions can be obtained using direct methods, which transform the original problem into a nonlinear programming (NLP) problem via .... COM...
1 downloads 4 Views 63KB Size
Ind. Eng. Chem. Res. 2002, 41, 2425-2430

2425

PROCESS DESIGN AND CONTROL Computation of Optimal Identification Experiments for Nonlinear Dynamic Process Models: a Stochastic Global Optimization Approach Julio R. Banga,*,† Karina J. Versyck,‡ and Jan F. Van Impe§ Process Engineering Group, Instituto de Investigaciones Marinas (CSIC), C/Eduardo Cabello 6, 36208 Vigo, Spain, Chemical Production Engineering, Janssen Pharmaceutica, Janssen Pharmaceuticalaan 3 (Building 210), B-2440 Geel, Belgium, and BioTeC (Bioprocess Technology and Control), Chemical Engineering Department, Katholieke Universiteit Leuven, Kasteelpark Arenberg 22, B-3001 Heverlee, Belgium

The problem of optimal experimental design (OED) for parameter estimation of nonlinear dynamic systems is considered. It is shown how this problem can be formulated as a dynamic optimization (optimal control) problem where the performance index is usually a scalar function of the Fisher information matrix. Numerical solutions can be obtained using direct methods, which transform the original problem into a nonlinear programming (NLP) problem via parametrizations. However, because of the frequent nonsmoothness of the cost functions, the use of gradient-based methods to solve this NLP might lead to local solutions. Stochastic methods of global optimization are suggested as robust alternatives. A case study considering the OED for parameter estimation in a fed-batch bioreactor is used to illustrate the performance and advantages of two selected stochastic algorithms. 1. Introduction Sound mathematical models are the fundamental elements of modern methods in process systems engineering. The development of a model is usually an iterative process which involves two basic steps: designing the model structure and estimating the model parameters for a given structure from experimental data. Here we will consider this latter task of parameter estimation considering dynamic models of chemical and biochemical processes. These processes are generally described by systems of differential and algebraic equations (DAEs), with the particularity that these models are frequently highly nonlinear and stiff. The parameter estimation problem can be formulated as finding the parameter vector p to minimize some functional Jwls of the weighted squared errors between a vector of measured outputs ym(t) and the model predictions y(p,t):

Jwls(p) )

∫0t (ym(t) - y(p,t))TQ(t) (ym(t) - y(p,t)) dt f

(1)

with Q(t) a square weighting matrix which is often defined as the inverse of the measurement error covariance matrix. Performing experiments to obtain a rich enough set of ym(t) is a costly and time-consuming activity, very especially when considering real industrial units. The * Corresponding author. E-mail: [email protected]. Phone: +34-986-214473. Fax: +34-986-292762. † CSIC. ‡ Janssen Pharmaceutica. § Katholieke Universiteit Leuven.

purpose of optimal experimental design (OED) is to devise the necessary dynamic experiments in such a way that the parameters are estimated from the resulting experimental data with the best possible statistical quality. The latter is usually a measure of the accuracy and/or decorrelation of the estimated parameters. Although the OED applied to linear steady-state models is a well-known subject, in the more challenging case of nonlinear dynamic models, such as those arising from (bio)chemical processes, fewer references are available.1-5 2. Mathematical Statement of the OED Problem Mathematically, the OED problem can be formulated as a dynamic optimization (optimal control) problem. The objective is to find a set of time-varying input variables (controls) for the experiments optimizing the quality of the estimated parameters in some statistical sense. This optimal control problem is subject to a set of equality constraints (the DAEs describing the system dynamics) and a set of inequality path and/or point constraints on the state variables. These inequality constraints are usually associated with restrictions concerning issues such as practical implementation, safety, and/or model validity. Consider the so-called Fisher information matrix F:

F)

∂y ∂y (t)) Q(t) ( (t)) dt ∫0t (∂p ∂p f

T

(2)

with between the brackets the sensitivity matrix of the model outputs (vector y) with respect to the model parameters (vector p) and Q the inverse of the measurement error covariance matrix. For systems linear

10.1021/ie010183d CCC: $22.00 © 2002 American Chemical Society Published on Web 04/17/2002

2426

Ind. Eng. Chem. Res., Vol. 41, No. 10, 2002

in the parameters and under certain noise assumptions, the inverse of this matrix F evaluated in the true process parameters p* is the error covariance matrix of the best linear unbiased estimator (BLUE).6 When the parameter estimates obtained from an experiment are close to the true process parameters p*, the Fisher information matrix F (eq 2) evaluated in these estimates can be seen as an approximation of the matrix F evaluated in the true process parameters p*. Consequently, its inverse matrix provides an approximation of the best attainable estimation error covariances. This is the basis for the experimental design methodology. For models which are nonlinear in the parameters, the application of this methodology involves the assumption that the output can be approximated by a first-order Taylor series expansion in the vicinity of the true parameters p*. In practice, the true process parameters are unknown. Therefore, during experimental design for parameter estimation, a so-called nominal parameter set p0 is to be used instead of the unknown true process parameters p*. Commonly, these nominal parameters p0 are obtained from preliminary experiments or the literature. Scalar functions φ of the Fisher information matrix F (eq 2) evaluated in the nominal parameters p0 are used as OED criteria for increasing the practical identifiability of the model parameters from experimental data. Because the nominal parameters p0 are generally not coinciding with the true process parameters p*, an iterative design scheme must be pursued to obtain eventually a truly optimal experiment with respect to the true set p*. In each iteration cycle, the parameter values that are estimated in the preceding experiment are used as the nominal parameter set p0 during experimental design. Thus, the mathematical formulation of the OED problem can be written as follows: Find u(t) to minimize J ) φ(F) subject to

f[x3 ,x,p,u,t] ) 0, x(t0) ) x0

(3)

h[x,y,p,u,t] ) 0

(4)

g[x,y,p,u,t] e 0

(5)

xL e x e xU

(6)

uL e u e uU

(7)

where J is the performance index (cost function), x is the vector of state variables, u is the vector of control (input) variables, y is the vector of output variables (measurements), eq 3 is the system of ordinary-differential equality constraints with its initial conditions, eqs 4 and 5 are the equality and inequality algebraic constraints, and eqs 6 and 7 are the upper and lower bounds for the state and control variables. As mentioned above, during optimization, the model parameters p are fixed at a nominal parameter set p0. The performance index is a scalar function φ(F) of the Fisher information matrix with output sensitivities evaluated in the nominal parameters p0. Because no analytical solutions for the sensitivity functions are available, they result from the numerical integration of an associated set of differential equations. Examples of common functions are the so-called D criterion (maximizing the determinant of F), where J ) -|F|, which measures the global accuracy of the estimated parameters, and the modified E criterion (minimizing the condition number of F, i.e.,

the ratio between the maximum and minimum eigenvalues), where J ) Λ(F) ) λmax(F)/λmin(F), which measures the parameter decorrelation. Thus, OED of these processes implies the solution of a class of very challenging nonlinear, highly constrained dynamic optimization problems. Besides, additional difficulties may arise because of the nonsmooth nature of the abovementioned cost functions. 3. Solution Methods Dynamic optimization methods can be classified into two major groups: indirect approaches, based on the minimum principle of Pontryagin, and direct approaches, which transform the original problem into a nonlinear programming (NLP) problem via discretization of the control and/or the state variables. In the case of the problem stated above, indirect approaches are not applicable because of the nondifferentiable nature of the costs. Thus, one must resort to direct approaches, reducing the problem to the solution of an NLP, following either the control vector parametrization (CVP) strategy7,8 or the complete (controls and states via collocation) parametrization strategy.9 In any case, standard gradient-based NLP solvers may fail to converge, or converge to a local solution, because of the nonsmooth and highly nonlinear nature of many problems.10,11 In this study, we have followed the CVP approach, where the control variables are approximated by some type of discretization, usually based on low-order polynomials (e.g., piecewise constant or piecewise linear). This discretization results in an NLP main problem, where the decision variables are the coefficients used to build the approximation. This main NLP problem requires the solution of an inner initial value problem (IVP), formed by eqs 3 and 4, for each function evaluation. The most popular methods for solving the main NLP are gradient-based local optimization techniques, usually variants of the powerful sequential quadratic programming (SQP) method. However, as was already mentioned, because this NLP may be nonconvex due to the nonsmooth nature of the performance indexes, SQP solvers may converge to local solutions or even fail to converge. Global optimization techniques are needed in order to surmount these difficulties. The use of stochastic or hybrid stochastic-deterministic methods have been suggested as efficient and robust alternatives for similar difficult optimal control problems.12,13 Here, we will apply two stochastic techniques, namely, integrated controlled random search (ICRS) and differential evolution (DE), which have been found particularly suitable for these types of problems.14 ICRS is an adaptive stochastic method developed by Banga and Casares,15 while DE is a population-based stochastic method recently presented by Storn and Price.16 Both methods are very simple to implement and use, which is an important characteristic for many potential users. The computational implementation of this approach is, in fact, rather easy. The user only needs to set up the global optimization stochastic solvers (details and kernel codes for that purpose are given elsewhere12,15,16) coupled with an initial value solver (i.e., an ODE or DAE solver) and the process dynamic model, augmented with the first-order sensitivities needed for the computation of the Fisher information matrix. Adequate solvers are widely available in many languages and/or environ-

Ind. Eng. Chem. Res., Vol. 41, No. 10, 2002 2427

ments, ranging from low-level FORTRAN solvers (e.g., LSODA or DASSL) to high-level environments, which can be general purpose systems such as Matlab (The Mathworks Inc., Natick, MA) or more specific dynamic simulation packages such as gPROMS (PSE Ltd., London, U.K.), Abacuss II (MIT, Cambridge, MA), or EcosimPro (EA Int., Madrid, Spain). All of them allow adequate interfacing with the external global optimization solvers. 4. Case Studies The presented method is general; i.e., it is applicable to any nonlinear dynamic process model. In this section we illustrate its advantages considering an important problem in biotechnology: the estimation of kinetic parameters of unstructured microbial growth models. It should be noted that this case study was chosen as an illustrative example which has received wide attention in the literature. In particular, there are two main reasons to explain its importance: (i) This class of models is an excellent example of a “simple”, common model for which standard experimental approaches result in estimated parameters of very poor statistical quality. (ii) Fermentation experiments are usually rather time-consuming, so bioprocess engineering is an area where the benefits of the proposed approach (in terms of both reducing the experimental work and increasing the statistical quality) will be especially meaningful. 4.1. Mathematical Model. The dynamics of a fedbatch bioreactor where one biomass is growing on one limiting substrate are described by the following model:

Fin dCS ) - σ(CS) CX + (C - C S) dt V S,in

(8)

Fin dCX ) µ(CS) CX C dt V X

(9)

dV ) Fin dt

(10)

with CS [g/L] the (volumetric) concentration of the limiting substrate, CX [g of DW/L] the biomass concentration, V [L] the volume of the liquid phase, CS,in [g/L] the substrate concentration in the influent with flow rate Fin [L/h], σ [g/(g of DW)‚h] the specific substrate consumption rate, and µ [1/h] the specific growth rate. These rates are functions of CS:

σ(CS) )

1 µ(CS) + m YX/S

CS µ(CS) ) µm Kp + CS + CS2/Ki

(11)

(12)

with YX/S [g of DW/g] the biomass on the substrate yield coefficient and m [g/(g of DW)‚h] the (overall) specific maintenance demand. Equation 12 is known as the nonmonotonic Haldane growth kinetics, with parameter Kp [g/L] indicating how fast the optimum for the specific growth rate µ is reached and Ki [g/L] the inhibition parameter. The kinetics is depicted as a function of the substrate concentration CS in Figure 1 (see section 4.2 and Table 1 for values of the parameters).

Figure 1. Haldane kinetics: specific growth rate as a function of the substrate concentration (CS,µ is the substrate concentration where the specific growth rate µ reaches its maximum value). Table 1. Initial Conditions, Operating Conditions, and Other Model Parameters for the Bioreactor Case Study µm KP,0 KI,0

2.1 10 0.1

X(0) S(0) V*

1/h g/L g/L

Parameters m YX/S CS,in a

0.29 0.47 500 1500

Initial Conditions 10.5 free 7

g/(g of DW)‚h g of DW/g g/L g g of DW g L

In this case study, our objective is to estimate the parameters Kp and Ki based on measurements of the two model outputs which are defined as the substrate concentration CS and the biomass concentration CX. The specific growth rate parameter µm [1/h] is assumed to be known. The volumetric feed rate u(t) into the fedbatch bioreactor is considered as the control input. Hence, a single input-multiple output system with two unknown parameters is under study. The information matrix F (eq 1) can be computed from the sensitivities obtained by solving an extended IVP. 4.2. Initial Conditions, Operation Conditions, and Model Parameters. The values for initial conditions, operating conditions, and other model parameters are listed in Table 1. These values are adopted from the literature which deals with the optimization of the penicillin G fed-batch fermentation.17,18 The following relations hold between the initial bioreactor states (with X(0) the initial biomass amount):

V(0) )

V*CS,in CS,in - CS(0)

CX ) X(0)/V(0)

(13) (14)

with V* the initial volume without substrate. The maximum volume is fixed at VMAX ) 10 L. Observe that this physical constraint is equivalent to fixing the total amount of the substrate available to R ) 1500 g (see the differential equation for V in the bioreactor model):

R ) S(0) +

∫0t CS,inu(t) dt f

(15)

with S(0) [g] the initially supplied amount of substrate. The maximum feed rate is limited by the pump capac-

2428

Ind. Eng. Chem. Res., Vol. 41, No. 10, 2002

ity: UMAX ) 1 L/h. For calculation of the weighting matrix Q in the information matrix F, an additive, uncorrelated, zero-mean white Gaussian noise is assumed on the measurements of the substrate concentration CS and the biomass concentration CX. For the cases discussed in this paper, the following values are used as the (time-invariant) variances of the measurement errors in the weighting matrix Q: σCS2 ) 10-2 g2/L2 and σCX2 ) 6.25 × 10-4 g2/ L2, respectively. 4.3. Problem Statement. The OED problem is formulated as finding the optimal substrate feed rate so that maximum global accuracy and/or decorrelation of the estimates for Kp and Ki is obtained, subject to the differential equality constraints (model dynamics, eqs 8-12) plus bounds for the control and the states. The following two subproblems are considered: (A) Optimal control for parameter decorrelation: the cost function is the modified E cost (ratio of the maximum and minimum eigenvalues of matrix F). (B) Optimal control for parameter accuracy: the cost function is the D cost (determinant of matrix F). Note that the nominal values for KP,0 and KI,0 (10 and 0.1 g/L, respectively) differ by 2 orders of magnitude. An appropriate weighting of the output sensitivities with respect to the model parameters can be introduced into the matrix F by multiplication of each sensitivity function with the corresponding nominal parameter value (see, e.g., work by Munack19). In this way the experimental design problem is formulated in terms of the dimensionless parameters KP/KP,0 and KI/KI,0. Consequently, the minimization of the condition number of the matrix F in subproblem a corresponds to the equalization of the relative estimation errors on both parameters. In this manner, a funnel-shaped error functional Jwls in the dimensionless parameter space (resulting in circular contour lines of Jwls in the dimensionless parameter plane) is aimed at. As mentioned before, the D criterion in subproblem b is not affected by scaling of the parameters (see, e.g., work by Mehra20). 4.4. Results and Discussion. Preliminary runs of the numerical optimization algorithms indicated that path inequality constraints on CS must be imposed in order to guarantee model validity. More specifically, a freakish evolution of the substrate concentration is not acceptable because the basic assumption of balanced growth would not apply. This assumption permits the use of unstructured growth kinetics. Balanced growth implies that the conversion rates of all cell compounds are equal such that they can be represented by one single overall specific growth rate.21 Considering subproblem A, it should be noted that the modified E cost should be as close as possible to 1, which corresponds to a situation of complete decorrelation (or equal relative parameter estimation errors). Versyck et al.4 followed a heuristic procedure, based on theoretical analysis of the optimal process performance feed rate profile, to design optimal control inputs which provide the desired complete decorrelation. Here, we started from no a priori assumptions and solved the problem in a purely numerical way. Control parametrizations based on piecewise constant and piecewise linear polynomials were tested considering 10-20 time elements of varying length. The piecewise constant approximation (steps) was found to give the best results. It should be noted that this parametrization is generally preferable from the control point of view because it is easier to implement in practice.

Figure 2. Typical convergence curves of ICRS and DE for subproblem A.

Figure 3. Optimal feed rate by Versyck et al.4 for subproblem A.

Apart from the stochastic techniques mentioned in the previous section, we also tried to solve the problem by using a deterministic (SQP-based) solver. Many difficulties were encountered by the SQP solver, and in fact convergence failures or convergence to local solutions after excessive computation times was found, in agreement with the reasoning presented previously about the nonsmoothness of the cost function. In contrast, the stochastic solvers arrived to very good solutions (cost functions close to 1 within a tolerance of 10-5 in most cases) in very reasonable computation times (order of minutes in a PC Pentium-II). Typical convergence curves (performance index versus computation time) are presented in Figure 2 (note the log-log axes). The ICRS method was clearly faster in this case, arriving at a solution very close to the global in just 10 s. Note that the curve of DE starts later than that of ICRS because of the initial time spent in the generation of the initial population. These solutions, which in this case we know are globally optimal, agree quite well with those of Versyck et al.,4 obtained via a heuristic procedure based on a conjecture. For example, Figure 3 shows the optimal feed rate found by these authors, while those obtained with the stochastic methods are presented in Figure 4 (for the particular case of fixed tf ) 86.3 h and CS(0) ) 38.4 g/L).

Ind. Eng. Chem. Res., Vol. 41, No. 10, 2002 2429

5. Conclusions and Future Work

Figure 4. Optimal feed rates obtained with ICRS and DE for subproblem A.

In this paper, we compared several numerical strategies for solving two OED problems in the context of parameter estimation for (bio)chemical process models. The stochastic global optimization methods show the best performance in terms of the criterion values and the computational efficiency for these nonlinear and highly constrained dynamic optimization problems. One additional advantage of these stochastic global optimization methods is that parallel versions can be created very easily. In particular, parallel versions capable of running in standard clusters of commodity PCs (a` la Beowulf; see www.beowulf.org) are being developed. It is envisaged that these versions will allow the solution of much larger OED problems in very modest computation times. Further, this will be done using already available local networks of standard PCs, thus avoiding additional investment in hardware. Finally, it should be noted that the sequence of experimental design-experiment execution-parameter estimation constitutes one step in an iterative design cycle. Evidently, there is always a need for an initial experiment which can be designed using prior information. The outcome of this initial experiment is used further in the subsequent design. The general robustness of the technique is the subject of ongoing work, but preliminary results have been very promising.23 Acknowledgment

Figure 5. Optimal feed rates obtained with ICRS and DE for subproblem B.

Considering subproblem B, the determinant of F must be maximized. There is no a priori knowledge about the globally optimal value for this cost function, so we have to limit ourselves to compare the solutions obtained by the different stochastic solvers and the heuristic solution obtained by Versyck et al.4 These authors obtained a D cost of |F| ) 2.2143 × 1014, with CS(0) ) 68.5 g/L and tf ) 218 h. Here, the problem was solved for several time horizons with the different solvers. In this way, a compromise between two criterions (statistical quality and process time) can be reached by using the results as a pareto-optimal response curve. In agreement with the previous results, the SQP-based deterministic solver failed to arrive at good solutions. In contrast, the ICRS stochastic method arrived at a better value of |F| ) 1.9770 × 1016 for a shorter process time of 86.3 h and CS(0) ) 38.1 g/L. The DE method clearly outperformed the ICRS method in this subproblem, arriving at a value of |F| ) 2.7155 × 1018, with CS(0) ) 31.7 g/L and the same process time. The corresponding feed rates obtained with the stochastic methods are shown in Figure 5. When considering the optimization results obtained for subproblem B with respect to the modified E cost in A, and vice versa, it should be noted that the minimization of either cost function is at the expense of the alternative cost. For example, the best solution (obtained by DE) for subproblem B yields a modified E cost of 5.1580 × 1011. This confirms the results reported for heuristic solution procedures in Versyck and Van Impe.22 Current research focuses on the stochastic optimization of a novel cost combining both criteria in a weighted sum. As such, a tradeoff between the two parameter estimation quality measures is aimed at.

K.J.V. contributed to this paper as a research assistant with the Fund for Scientific ResearchsFlanders (FWO). This work was supported in part by Project OT/ 99/24 of the Research Council of the Katholieke Universiteit Leuven and the Belgian Program on Interuniversity Poles of Attraction, initiated by the Belgian State, Prime Minister’s Office for Science, Technology and Culture. The Fund for Scientific ResearchsFlanders (FWO) is acknowledged for the study journey grant that allowed us to perform part of this research at the Instituto de Investigaciones Marinas (CSIC) in Vigo, Spain. Scientific responsibility is assumed by the authors. Literature Cited (1) Munack, A.; Posten, C. Design of optimal dynamical experiments for parameter estimation. Proceedings of the American Control Conference, ACC89, Pittsburgh, PA, 1989; American Automatic Control Council; Northwestern University: Evanston, IL, 1989; pp 2011-2016. (2) Lohmann, Th. W.; Bock, H. G.; Schlo¨der, J. P. Numerical methods for parameter estimation and optimal experiment design in chemical reaction systems. Ind. Eng. Chem. Res. 1992, 31, 54. (3) Baltes, M.; Schneider, R.; Sturm, C.; Reuss, M. Optimal experimental design for parameter estimation in unstructured growth models. Biotechnol. Prog. 1994, 10, 480. (4) Versyck, K. J.; Claes, J.; Van Impe, J. Practical identification of unstructured growth kinetics by application of optimal experimental design. Biotechnol. Prog. 1997, 13, 524. (5) Bauer, I.; Bock, H. G.; Ko¨rkel, S.; Schlo¨der, J. P. Numerical methods for initial value problems and derivative generation for DAE models with application to optimum experimental design of chemical processes. Proceedings of the International Workshop on Scientific Computing in Chemical Engineering, TU Hamburg, Germany, May 26-28, 1999; Vol. II. (6) Walter, E.; Pronzato, L. Identification of parametric models from experimental data; Springer: London, 1997. (7) Vassiliadis, V. S.; Sargent, R. W. H.; Pantelides, C. C. Solution of a class of multistage dynamic optimization problems, Part I. Ind. Eng. Chem. Res. 1994, 33, 2111.

2430

Ind. Eng. Chem. Res., Vol. 41, No. 10, 2002

(8) Balsa-Canto, E.; Banga, J. R.; Alonso, A. A.; Vassiliadis, V. S. Efficient optimal control of bioprocesses using second-order information. Ind. Eng. Chem. Res. 2000, 39, 4287. (9) Cuthrell, J. E.; Biegler, L. T. Simultaneous optimization and solution methods for batch reactor control profiles. Comput. Chem. Eng. 1989, 13, 49. (10) Barton, P. I.; Allgor, R. J.; Feehery, W. F.; Galan, S. Dynamic optimization in a discontinuous world. Ind. Eng. Chem. Res. 1998, 37, 966. (11) Barton, P. I.; Banga, J. R.; Galan, S. Optimization of hybrid discrete/continuous dynamic systems. Comput. Chem. Eng. 2000, 24, 2171. (12) Banga, J. R.; Seider, W. D. Global optimization of chemical processes using stochastic algorithms. In State of the Art in Global Optimization; Floudas, C. A., Pardalos, P. M., Eds.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996; pp 563583. (13) Banga, J. R.; Alonso, A. A.; Singh, R. P. Stochastic dynamic optimization of batch and semicontinuos bioprocesses. Biotechnol. Prog. 1997, 13, 326. (14) Balsa-Canto, E.; Alonso, A. A.; Banga, J. R. Dynamic optimization of bioprocesses: deterministic and stochastic strategies. Proceedings of ACoFoP IV, Go¨teborg, Sweden, Sept 21-23, 1998. (15) Banga, J. R.; Casares, J. J. Integrated Controlled Random Search: application to a wastewater treatment plant model. IChemE Symp. Ser. 1987, 100, 183. (16) Storn, R.; Price, K. Differential Evolutionsa simple and efficient heuristic for global optimization over continuous spaces. J. Global Optim. 1997, 11, 341. (17) Lim, H. C.; Tayeb, Y. J.; Modak, J. M.; Bonte, P. Computational algorithms for optimal feed rates for a class of fed-batch

fermentation: numerical results for penicillin and cell mass production. Biotechnol. Bioeng. 1986, 28, 1408. (18) Van Impe, J. F.; Nicolaı¨, B.; Vanrolleghem, P.; Spriet, J.; De Moor, B.; Vandewalle, J. Optimal control of penicillin G fedbatch fermentation processes with growth/production decoupling. Control Eng. Pract. 1992, 3, 939. (19) Munack, A. Optimization of sampling. In Biotechnology (Measuring, modelling, and control); Rehm, H.-J., Reed, G., Eds.; VCH: Weinheim, Germany, 1991; Vol. 4; pp 252-264. (20) Mehra, R. Optimal input signal for parameter estimation in dynamic systemsssurvey and new results. IEEE Trans. Autom. Control 1974, 19, 753. (21) Roels, J. A. Energetics and kinetics in biotechnology; Elsevier: New York, 1983. (22) Versyck, K. J.; Van Impe, J. F. Trade-offs in design of fedbatch experiments for optimal estimation of biokinetic parameters. Proceedings of the 1998 IEEE Conference on Control Applications, Control Applications in Biological Systems, Trieste, Italy; IEEE Catalog No. 98CH36104; 1998; pp 51-55. (23) Versyck, K. J.; Van Impe, J. F. Robustness of the iterative optimal experiment design scheme for parameter estimation: identification of microbial heat resistance parameters. In Proceedings of the 8th International Conference on Computer Applications in Biotechnology; Dochain, D., Perrier, M., Eds.; Elsevier Science: New York, 2001; pp 31-36.

Received for review February 22, 2001 Revised manuscript received February 10, 2002 Accepted February 19, 2002 IE010183D