Global optimization of bioprocesses using stochastic

0 downloads 0 Views 344KB Size Report
Keywords: bioprocess engineering, global optimization, nonlinear dynamic systems, ... solution strategies, both deterministic and stochastic, which we have used and compared ... The main reason is that, for many bioprocesses, optimal ..... Clustering methods were derived from the initial concepts of multi-start methods,.
Banga, J. R., Moles, C. G., Alonso, A. A. (2003) Global Optimization of Bioprocesses using Stochastic and Hybrid Methods. In "Frontiers in Global Optimization", C. A. Floudas and P. M. Pardalos (Eds.), Nonconvex Optimization and Its Applications, vol. 74, pp. 45-70, Kluwer Academic Publishers, ISBN 1-4020-7699-1.

Frontiers In Global Optimization C. A. Floudas and P. M. Pardalos, Editors c °2003 Kluwer Academic Publishers

Global optimization of bioprocesses using stochastic and hybrid methods Julio R. Banga Process Engineering Group IIM-CSIC C/Eduardo Cabello 6, 36208 Vigo (SPAIN) [email protected] Carmen G. Moles Process Engineering Group IIM-CSIC C/Eduardo Cabello 6, 36208 Vigo (SPAIN) [email protected] Antonio A. Alonso Process Engineering Group IIM-CSIC C/Eduardo Cabello 6, 36208 Vigo (SPAIN) [email protected]

Abstract In this contribution, we will focus on problems arising in the context of biochemical process engineering. Many of these problems can be stated as the optimization of non-linear dynamic systems. Relevant classes in this domain are (i) optimal control problems (dynamic optimization), (ii) inverse problems (parameter estimation), and (iii) simultaneous design and control optimization problems. Most of these problems are, or can be transformed to, nonlinear programming problems subject to differentialalgebraic constraints. It should be noted that their highly constrained and non-linear nature often causes non-convexity, thus global optimization methods are needed to find suitable solutions. Here, we will present our experiences regarding the use of several stochastic, deterministic and hybrid global optimization methods to solve those problems. Several parallel versions of the most promising methods, which are able to run on standard clusters of PCs, will also be presented. Results for selected challenging case studies will be given. Keywords: bioprocess engineering, global optimization, nonlinear dynamic systems, optimal control, integrated design

1

1

Introduction

Systems engineering approaches are increasingly applied in the bioprocess industries (i.e., biotechnological, food, pharmaceutical, environmental, etc.). In order to increase the productivity, profitability and/or efficiency of bioprocesses, considerable research effort has been devoted to their improvement via computer aided process engineering methods. In this way, mathematical modeling, optimization and control have become fundamental tools to optimally design and operate production facilities in these sectors [77] [2] [9] [8]. During the last decade, our group has been especially interested in robust and efficient optimization techniques which can be used, in combination with suitable models, to obtain optimal or near-optimal solutions regarding the design, operation and control of bioprocesses. Since most bioprocess are operated in batch or semi-continuous modes, they have an inherent dynamic nature. Thus, there is a need to use methods designed for the optimization of dynamic systems. In this context, there are three types of optimization problems which are especially relevant: • Optimal operation: given a process dynamic model and a set of specifications (constraints, goals), the objective is to compute the optimal operating conditions which lead to maximum performance as measured by some pre-defined criteria. These problems belong to the domain of dynamic optimization (or open loop optimal control) • Integrated process design: to find simultaneously the static design variables (e.g. sizes and number of units), the operating conditions (e.g. flows) and other design issues (e.g. the controllers) which minimize capital and operation costs while optimizing certain characteristics of the dynamics of the process (e.g. maximizing controllability) • Model calibration: this is the well known problem of parameter estimation (inverse problem), i.e., to find the parameters of a nonlinear dynamic model which give the best fit to a set of experimental data. This is one of the mandatory steps in dynamic model development, but unfortunately many modelers are not aware of the dangers of applying standard optimization methods for its solution These problems are, or can be transformed to, nonlinear programming problems subject to dynamic (usually, differential-algebraic) constraints. Their highly constrained, non-linear and sometimes non-smooth nature often causes non-convexity, thus global optimization methods are needed to find suitable solutions. This paper is structured as follows. In the next sections, the above three classes of problems will be examined, briefly reviewing the state of the art regarding solution methods, and highlighting the need of using global optimization. Then, we will discuss several GO solution strategies, both deterministic and stochastic, which we have used and compared based on their results for a number of problems. Several hybrid (stochastic-deterministic) approaches will also be presented and evaluated, since our experience indicates that they can bring significant advantages in terms of robustness and computational effort. Results for selected challenging case studies will be outlined, finalizing with a set of conclusions. 2

2

Dynamic optimization

The dynamic optimization (open loop optimal control) of bioprocesses has been, and continues to be, a hot research topic. The main reason is that, for many bioprocesses, optimal operating policies can improve productivity and profitability in a very significant way. A good example is fermentation, which can be carried out in continuous, batch and fed-batch modes. In fed-batch fermentation, cells or micro-organisms are grown in a bioreactor where nutrient(s) are provided along the process using a controlled (time-varying) feed. Fed-batch bioreactors have a number of well known advantages over batch or continuous fermentors. For example, fed-batch can be the best (or even the only) alternative due to its effectiveness in overcoming undesired effects like substrate inhibition and catabolite repression. Besides, it provides better control of deviations in the organism’s growth pattern, and production of high cell densities can be possible due to extension of process time, which is useful for the production of substances associated with growth. Typical product revenues for many pharmaceutical products obtained through fed-batch fermentation are in the order of billions of USD [19]. Not surprisingly, the dynamic optimization of fed-batch bioreactors is a class of problems that has received major attention during the last decades [9]. Basically, dynamic optimization allows the computation of the optimal operating policies for these units, i.e. the best time-varying feed rate(s) which ensure the maximization of a pre-defined performance index (usually, a productivity, or an economical index derived from the operation profile and the final concentrations). Dynamic optimization has also been successfully used for computing optimal operating policies in food process engineering. In particular, many studies have illustrated the benefits of using optimal controls in several important operations like thermal processing (sterilization, pasteurization), drying, contact cooking or microwave heating [8].

2.1

Dynamic optimization: problem statement

The general dynamic optimization (optimal control) problem of a bioprocess, considering a free terminal time, can be stated as finding the control vector(s) u{t} and the final time tf to minimize (or maximize) a performance index J [x, u]: J [x, u] = Θ [x {tf }] +

Z tf t0

Φ [x {t} , u {t} , t] dt

(1)

subject to the dynamics of the system, usually expressed as a set of ordinary differential equality constraints, Eqns. (2): dx = Ψ [x {t} , u {t} , t] dt

(2)

where x is the vector of state variables, with initial conditions x{t0 } = x0 , and also subject to sets of algebraic equality and inequality constraints, Eqns. (3) and (4): h [x {t} , u {t}] = 0

(3)

g [x {t} , u {t}] ≤ 0

(4)

3

An additional set of algebraic inequality constraints are the upper and lower bounds on the state and control variables, Eqns. (5) and (6): xL ≤ x {t} ≤ xU u

L

(5)

U

≤ u {t} ≤ u

(6)

The above formulation assumes that the process is modelled as a lumped system (i.e., described by ordinary differential equations, so no spatial distributions are taken into account). If the process is modelled as a distributed system (e.g., most processes in the food industry, or bioreactors where state variables are function of both time and spatial position), the corresponding governing partial differential equations (PDEs) are introduced as an additional set of equality constraints: p (x, xκ , xκκ , . . . , κ) = 0 a (x, xκ , . . . , κ) = 0

κ∈Ω

(7)

κ ∈ δΩ

(8)

where κ are the independent variables (time and spatial position), Eqns. (7) are the system of governing PDEs within the domain Ω, where xκ = ∂x/∂κ, and Eqns. (8) are the auxiliary conditions of the PDEs (boundary and initial conditions) on the boundary δΩ of the domain Ω. In many cases, especially in food process engineering, the dynamics can be much more complex, like e.g. coupled heat and mass transfer involving dimensional changes (moving fronts) and phase change, multiphase flows, etc.

2.2

Dynamic optimization: local solution methods

The dynamic optimization of bioprocesses can be a very challenging problem due to several reasons. Frequently, the control variable (e.g. feed rate to a bioreactor) appears linearly in the system differential equations, so the problem is singular, creating additional difficulties for its solution (especially using indirect methods). For this type of problems, the optimal operating policy will be either bang-bang, or singular, or a combination of both. Further, most bioprocesses have highly nonlinear dynamics, and constraints are also frequently present on both the state and the control variables. These characteristics introduce new challenges to the existing solution techniques. Therefore, efficient and robust methods are needed in order to obtain the optimal operating policies. Numerical methods for the solution of optimal control problems are usually classified under two categories: indirect and direct approaches. Indirect (classical) approaches are based on the transformation of the original optimal control problem into a two point boundary value problem using Pontryagin’s necessary conditions [24]. Although many researches have followed this approach for the optimization of fed-batch reactors (e.g.[72]), the resulting boundary value problems (BVPs) can be very difficult to solve, especially when state constraints are present. Alternatively, direct approaches transform the original optimal control problem (which is infinite dimensional) into a non-linear programming (NLP) problem, either using control vector parameterization (CVP) ([83]) or complete (control and state) parameterization [30]. The CVP approach transforms the original dynamic optimization problem into a nonlinear programming (NLP) problem. Gradient-based local methods are the best option 4

for solving NLPs, provided that these problems are unimodal (i.e. single optimum) and smooth. Sequential Quadratic Programming (SQP) methods are usually recognized as the state of the art in this domain. In this scheme, gradients are usually estimated using either finite differences, adjoints or first order sensitivities, the latter being the preferred approach [83]. The simultaneous integration of the system’s dynamics with the first order sensitivities provides both the states and the needed gradients accurately and with a low computational cost. Recently, this approach has been extended using second order sensitivities [84], so the estimation of the Hessian can be done efficiently. Using second order information of high quality, SQP methods can solve the master NLP more efficiently. Following these ideas, Balsa-Canto et al [4] have presented a CVP method which makes use of restricted second order information and a mesh refinement procedure in order to solve these problems in a very efficient way, even for high levels of control discretization. These authors have demonstrated how this new method allows a rapid and robust solution of several challenging bioreactor optimization problems. However, solving the NLPs arising from direct approaches like CVP is not always trivial. On the one hand, these NLPs are frequently multimodal (nonconvex, i.e. presenting multiple local optima), due to the highly nonlinear and constrained nature of the dynamics [12], or to the presence of discontinuities [16]. On the other hand, an inner initial value problem must be solved iteratively within the master NLP. If the integration tolerances are not tight enough (or are similar to the optimization ones), then the resulting ”numerical noise” will make the NLP non-smooth, leading the NLP solver to bad convergence and/or early stops. Therefore, deterministic (gradient based) local optimization techniques may converge to local optima, especially if they are started far away from the global solution. In order to surmount these difficulties, global optimization methods must be used.

3

Integrated process design

During the last decade, the importance of a simultaneous (integrated) design approach, considering operability together with the economic issues, has been recognized (e.g. see [60] [75] [63] [14] [35] and the references cited therein). It should be noted that the optimization problems arising from these formulations are very challenging. The multimodal (non-convex) nature of these problems has been highlighted by e.g. Schweiger and Floudas [75] and Bansal et al. [14], among others.

3.1

Integrated process design: problem statement

A general statement can be considered taking into account process and control superstructures which indicate the different design alternatives (e.g. [75] [14]). These statements result in mixed integer optimal control problems (MIOCPs). A simpler case is often considered, where it is assumed that the process flowsheet is given, as well as the process specifications. Although this problem statement is obviously simpler than the above mentioned, it has been shown to be challenging enough for many optimization methods. Besides, it is a case often encountered in the real world, where many 5

bioprocesses have well established process flowsheets, so the process and control superstructures are not an issue for the integrated design problem. For this latter case, the objective is to simultaneously find the static variables of the process design, the operating conditions and the parameters of the controllers which optimize a combined measure of the plant economics and its controllability, subject to a set of constraints which ensure appropriate dynamic behavior and process specifications. Mathematically, the statement is: Find v to minimize C=

P

αi φi

(9)

f (x, x, p, v) = 0

(10)

x(t0 ) = x0

(11)

h(x, p, v) = 0

(12)

g(x, p, v) = 0

(13)

subject to ·

L

v ≤v≤v

U

(14)

where v is the vector of decision variables (equipment static design variables, operating conditions and tuning parameters of controllers), C is the cost (objective function) to minimize (typically, a weighted a combination of capital, operation and controllability costs, φ1 , φ2 , φ3 ), f is the set of differential and algebraic equality constraints describing the system dynamics (i.e. the nonlinear process model), and h and g are possible equality and inequality path and point constraints which express additional requirements for the process performance.

3.2

Integrated process design: local solution methods

The formulation above is that of a non-linear programming problem (NLP) with differentialalgebraic (DAEs) constraints, i.e. similar to the one resulting from applying the CVP approach to dynamic optimization problems. Again, due to the nonlinear and constrained nature of the system dynamics, these problems are very often multimodal (non-convex). Further, it is known that using standard controllability measures, such as the Integral Square Error (ISE), in the objective function often causes non-convexity [75]. Therefore, if this NLP-DAEs is solved via standard NLP methods, such as Sequential Quadratic Programming (SQP), it is very likely that the solution found will be of local nature. It is often argued that the non-convexity of most process optimization problems can be surmounted by using a local method repeatedly, starting from a number of different initial decision vectors (this is the so-called multi-start strategy). However, as we will show in the following, even the use of a large number of trials within a multi-start SQP can miss the global solution of integrated process design problems. Therefore, global optimization (GO) methods must be used. 6

4

Parameter estimation in dynamic bioprocess models

In parameter estimation (or inverse problems) the objective is to find the parameters of a nonlinear dynamic model which give the best fit to a set of experimental data. An excellent review and practical introduction, which includes an impressive and well documented list of applications (including bio related ones) has been recently presented by Schittkowski [73]. This is a topic of major importance, and it has received great attention in the particular domain of bioprocesses and biosystems engineering, since the calibration of nonlinear dynamic models can be rather challenging. For example, Mendes and Kell [53] provide several examples of the challenging nature of parameter estimation in biochemical pathways. Moreover, modeling large biological systems from functional genomic data is an emerging area which requires proper parameter estimation methods [52]. A related class of problems is the optimal experimental design of dynamic experiments, where the objective 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. A good example in bioprocess engineering is the practical identification of growth kinetics [85]. The problem can be posed as a dynamic optimization, and certain global optimization methods have already demonstrated their usefulness to ensure proper solutions [13]. However, for the sake of briefness, in this contribution we will only consider the parameter estimation problem.

4.1

Parameter estimation: problem statement

This type of problem, also known as data fitting in dynamic systems, or dynamic model calibration, is usually stated as minimizing a cost function that measures the goodness of the fit of the model with respect to a given experimental data set, subject to the dynamics of the system (acting as a set of differential equality constraints) plus possibly other algebraic constraints. Mathematically, the formulation is that of a non-linear programming problem (NLP) with differential-algebraic constraints: Find p to minimize J=

R tf 0

(ymsd (t) − y(p, t))T W (t)(ymsd (t) − y(p, t))dt

(15)

subject to ·

f (x, x, y, p, v, t) = 0

(16)

x(t0 ) = x0

(17)

h(x, y, p, v) = 0

(18)

g(x, y, p, v) = 0

(19)

L

U

p ≤p≤p

(20)

where J is the cost function to be minimized, p is the vector of decision variables of the optimization problem. The set of parameters to be estimated, ymsd are the experimental measures of a subset of the (so-called) output state variables, y(p, t) are the model predictions for those outputs, W (t) is a weighting (or scaling) matrix, x are the differential state 7

variables, and v is a vector of other (usually time-invariant) parameters which are not estimated. Besides, f is the set of differential and algebraic equality constraints describing the system dynamics (i.e. the nonlinear process model), and h and g are the possible equality and inequality path and point constraints which express additional requirements for the system performance. Finally, p is subject to upper and lower bounds acting as inequality constraints. The formulation above is that of a non-linear programming problem (NLP) with differential and algebraic (DAEs) constraints, i.e. similar to the integrated design problem. And, as already mentioned in the previous sections, due to the nonlinear and constrained nature of the system dynamics, these problems are very often multimodal (non-convex). Therefore, if this NLP-DAEs is solved via standard local methods for parameter estimation, such as the standard Levenberg-Marquardt method, it is very likely that the solution found will be of local nature, as discussed by e.g. Mendes and Kell [53] .

5

Global optimization of dynamic systems

As already mentioned, the application of direct methods (i.e. control vector parameterization or complete parameterization) to optimal control problems frequently leads to nonconvex NLPs subject to nonlinear differential-algebraic constraints. Similarly, the latter also often arise in the framework of integrated design or parameter estimation problems. The more naive approach to surmount nonconvexity, i.e. multi-start local methods, fails for any mildly realistic problem. Thus, there is a clear need of robust and efficient global optimization problems in order to ensure proper solutions. The global optimization (GO) of nonlinear dynamic systems is receiving increased attention from engineers, mathematicians and computer scientists. In the domain of deterministic GO methods, Esposito and Floudas [32] [33] have recently presented approaches to solve nonlinear optimal control (dynamic optimization) and parameter estimation problems. This is indeed a very promising and powerful approach, but the objective function and the dynamics of the system must be twice continuously differentiable, and restrictions may also apply for the type of path constraints which can be handled. Other groups [79] [62] are also making good progress in deterministic global optimization of dynamic systems, yet several issues regarding requirements and computational performance are still present. In any case, research along these lines continues and it might result in breakthrough results in the short term, Regarding stochastic GO methods, several researches have shown that they can locate the vicinity of global solutions for nonlinear dynamic problems with relative efficiency [11] [47] [6] [12] [7] [1] [86] , but the cost to pay is that global optimality can not be guaranteed. However, in many practical situations these methods can be satisfactory if they provide us with a ”good enough” (often, the best available) solution in modest computation times. Furthermore, stochastic methods are usually quite simple to implement and use, and they do not require transformation of the original problem, which can be treated as a black box. Thus, they can handle problems with complicated dynamics (e.g. discontinuities [15], non-smoothness, etc.). 8

5.1

Stochastic methods

Since the number and types of stochastic methods being used for global optimization is increasing rapidly, we provide here a rough classification, followed by a brief review of their applications in our domain of interest, the optimization of bioprocesses. A crude taxonomy is as follows: • Random search and adaptive stochastic methods: these methods have their origins in research performed during the ’50s and ’60s [23] [51] [64], with a number of more refined and efficient methods being developed during the last decade (e.g. [89] [1] [12] [81]). • Clustering methods were derived from the initial concepts of multi-start methods, i.e. local methods started from different initial points. Clustering methods are more efficient and robust than multi-start because they try to identify the vicinity of local optima, thus increasing efficiency by avoiding the repeated determination of the same local solutions [82] [66]. • Evolutionary computation: most of these algorithms were created following the ideas of biological evolution, but in fact, they can be regarded as population-based adaptive stochastic methods. At least three different types were developed independently in the late ’60s and early ’70s: Genetic Algorithms (GAs) [39] [36], Evolutionary Programming (EP) [34] and Evolution Strategies (ES) [38] [74] [17] [65] [18]. • Simulated Annealing (SA), and variants: developed (originally, for combinatorial optimization problems) by simulating certain natural phenomena taking place at the atomic level regarding the cooling of metals [43] [44]. • Other bio-inspired methods and metaheuristics: in recent years, a number of (so called) meta-heuristics have been presented, mostly based on other biological or physical phenomena, and with combinatorial optimization as their original domain of application. Examples of these more recent methods are Ant Colony Optimization (ACO) [22] [31] and particle swarm methods [21]. A review of these and other recent techniques can be found elsewhere [28] [55] Currently, Genetic algorithms (GAs) and simulated annealing (SA) are the most popular types of methods. However, as many authors have reported during recent years, and as our own experience has confirmed, GAs and SA are usually not the most efficient and robust algorithms for global optimization. In fact, for global optimization in real variables, many other simple techniques outperform both types of methods, which were originally developed with combinatorial problems (integer variables) in mind. In any case, the literature is very fragmented, and there is a lack of sound comparative studies. This complicates the selection of methods for a given type of GO problem. Furthermore, although the topic is still the subject of great debate, it should be noted that for the general GO problem with a priori unknown structure, there is no best method (the so-called ”no free lunch”, or NFL, theorem [87]). 9

However, for the case of global optimization of nonlinear dynamic processes, different recent works [3] [58] [59] indicate that certain simple stochastic methods, namely the Differential Evolution method [80] and certain Evolution Strategies (ES) [71] consistently present the best performance in terms of efficiency and robustness. Additionally, our own recent experience indicates that ES methods exhibit particularly nice properties regarding their use for real-world GO problems: • good efficiency and robustness • good scaling properties (almost linear with problem size in some cases) • inherent parallel nature, i.e. these methods lend themselves to parallelization very easily, which means that medium-to-large scale problems can be handled in reasonable wallclock time. Thus, our recent and current research on GO methods for NLP-DAEs is focused on extending the ES paradigm with two main issues in mind: • Devising better procedures for the handling of constraints (i.e. with a good compromise between efficiency and robustness) • Designing hybrid stochastic-deterministic methods, combining the robustness and global convergence properties of ES methods with the efficiency of gradient-based local search methods when they are started close to the global solution

5.2

Handling of constraints in stochastic methods

As discussed in previous sections, our problems of interest are highly constrained. Unfortunately, most stochastic methods, including Evolutionary Computation techniques, were originally devised for unconstrained optimization. In most cases, ad hoc procedures for constraint handling were later added, usually by means of either direct rejection of unfeasible points, or by means of penalty functions. It is not our intention to review here the many approaches suggested in the literature. Instead, the interested reader will find this information in [54], [57] and [27]. In his thorough survey of the state of the art on constraint handling techniques, Coello [27] concludes that, if little is known about the problem, a simple static or dynamic penaltybased approach is advisable, since it is very easy to implement and is also rather efficient. This author also states that more comparisons of different approaches are needed, taking special care regarding the metrics used in these exercises. Our own experience in this topic during the last years is in agreement with Coello’s recommendations, so in most of our applications we have used penalty functions to handle the constraints. However, the recent SRES (Stochastic Ranking Evolution Strategy) method developed by Runarsson and Yao [71] has presented very good efficiency and robustness even for large nonlinear problems. As we will describe in the next section regarding case studies, this method has consistently outperformed all the other methods tested for most problems from a large suite of constrained examples. All this empirical evidence lead us to recommend SRES as possibly the best stochastic method currently available for the optimization of constrained dynamic systems. 10

In any case, as noted by Coello [27], if the NFL theorems [87] hold, then one should not expect a best general constraint handling technique. The best techniques for each type of problems will have to exploit specific knowledge about that type.

5.3

Hybrid methods

Using the above mentioned stochastic methods, refined solutions are usually obtained at a very large computational cost, especially if the dynamic systems are large (i.e. arising from a finite element discretization). Although there is always a trade-off between convergence speed and robustness in both stochastic and deterministic (local) methods, the latter usually have the opposite behavior, i.e. they converge very fast if they are started close to the global solution. Clearly, a convenient approach would be to combine both methodologies in order to compensate for their weaknesses while enhancing their strengths. Hybrid approaches have been common for the solution of many tough numerical problems, and global optimization is no exception. In the case of dynamic optimization, a simple two-phase hybrid (stochastic-deterministic) approach was suggested by Banga and Seider [12], and later considered by Balsa-Canto et al. [3] and Carrasco and Banga [25] with very good results. This approach was developed by adequately combining the key elements of a stochastic and a deterministic method (usually, an SQP method), taking advantage of their complementary features. Other authors have also used different hybrid approaches for dynamic optimization, confirming their usefulness [26]. However, in certain cases the global optimum might be situated in a region of the search space with discontinuities and/or non-smoothness, so the use of a local deterministic method for the refinement phase of the hybrid might result in lack of convergence or premature stop. In order to surmount these difficulties, we have been exploring hybrid approaches where both phases are performed by stochastic methods, which are more robust. More information about this approach is given in the next sections.

6

Case Studies

Here we present results for selected challenging case studies regarding the dynamic optimization of bioprocesses, the simultaneous design and control of a wastewater treatment plant, and the robust estimation of kinetic parameters in complex reaction systems.

6.1

Optimal control of bioreactors

Considering one of the main topics of this contribution, the dynamic optimization of bioprocesses, adaptive stochastic methods have been proposed as robust alternatives to local gradient-based methods [6] [7] [3] [5] . Other researches have used alternative stochastic algorithms, including different random search algorithms, arriving to similar conclusions [78] [48] [67]. Genetic algorithms and related evolutionary algorithms, which have been used for the solution of general optimal control problems by several authors [56] [76] [88], have also been extensively used for the optimization of bioprocesses, like fed-batch fermentation, during the last decade (e.g. [3] [50] [70] [86] [46] [61] [68]). Other metaheuristics, like Ant Colony Optimization, have also been employed for the same purpose [41]. 11

Although partial comparisons of different types of methods have been presented [69] [3] [26] [45], more work is needed in order to arrive to meaningful conclusions. In particular, a uniform set of benchmark problems should be used (one has been suggested in [9]). Also, many studies present simple comparisons of computation times and final performance index values, but this can lead to wrong conclusions. Instead, more sound procedures, like the use of convergence curves (i.e. performance index versus computation time) should be used, together with a systematic approach for the estimation of the computational effort (e.g. the use of function evaluations for comparing computational effort should be avoided, since different methods may have very different overheads). Recently, Banga et al. [5] considered the general problem of dynamic optimization of bioprocesses with unspecified final time. Several solution strategies, both deterministic and stochastic, were compared based on their results for a set of case studies. These strategies included two types of gradient-based (local) control vector parameterization (CVP) methods, CVP-fd (gradient via finite differences) and CVP-sg (gradient via sensitivities), one complete parameterization (CP) local method, three different stochastic methods, ICRS/DS [7], GA [56] and DE [80], all of them using the CVP framework, and a two-phase hybrid method, TPH. The comparative evaluation of their efficiency and robustness indicated the superiority of the hybrid approach, which presented the best compromise between robustness and efficiency, confirming the preliminary results of Carrasco and Banga [25]. In order to illustrate these results, let us consider the case study of a fed-batch reactor for the production of ethanol, as studied by a number of authors [7]. The (free terminal time) optimal control problem is to maximize the yield of ethanol using the feed rate as the control variable. This problem has been solved using CVP and gradient based methods, but convergence problems have been frequently reported, something which has been confirmed by our own experience. Table 1: Ethanol bioreactor: summary of optimization results (relative CPU times and final relative error w.r.t. best known value) for several methods

CVP-fd CVP-sg CP-col ICRS/DS GA DE TPH

Rel. CPU time 1.13 0.69 2.77 1.27 51.08 6.49 1.00

% error 7.17 0.41 0.30 1.07 4.01 0.06 0.01

In our runs, the best result was obtained by the hybrid TPH method, which arrived to J = 20839 and tf = 61.17 h) in only 35 s of computation (PC PII), which is orders of magnitude faster than results obtained by other methods like IDP [47], and very close to the best reported solution, J = 20842 . We also found that no other method was faster in arriving at within 0.5% of this J value, confirming the good efficiency of TPH. After careful initialization, CP-col arrived close to that performance index but with a CPU time three times larger. Only CVP-sg and DE were able to reach within 0.5%, the latter with a 12

rather large computation effort. Also, the CVP-sg method required a careful initialization, which is typical of local gradient methods. The CVP-fd and GA methods failed quite dramatically, the latter at a huge computational cost (50 times that of TPH), reinforcing our belief that genetic algorithms, despite their popularity, are not the most competitive stochastic method. The detailed results (relative CPU time and % of error with respect to the best known solution, J = 20842) for all methods are shown in Table 1.

6.2

Integrated design of a WWT plant

This case study was based on a real wastewater treatment (WWT) plant placed in Manresa (Spain), as described by Guti´errez and Vega [37]. The plant is formed by two aeration tanks, acting as bioreactors, and two settlers, as depicted in Figure 1. A flocculating microbial population (biomass) is kept inside each bioreactor, transforming the biodegradable pollutants (substrate), with the aeration turbines providing the necessary level of disolved oxigen. The effluents from the aeration tanks are separated in their asociated settlers into a clean water stream and an activated sludge, which is recycled to the corresponding aeration tank. Since the activated sludge is constantly growing, more is produced that can be recycled to the tanks, so the excess is eliminated via a purge stream (qp ). The objective of the control system is to keep the substrate concentration at the output (s2 ) under a given admisible value. The main disturbances come from large variations in both the flowrate and substrate concentration (qi and si ) of the input stream. Although there are several possibilities for the manipulated variable, here we have considered the flowrate of the sludge recycle to the first aeration tank, as originally considered by Guti´errez and Vega [37]. qi si xi q12 s1

sir1

x1 c1

xir1

xd1 xb1 xr1

q1 s1

x2 c2 fk2

fk1

qr1

q22 s2

qr2 sir2 xir2

q2

qsal xd2 xb2 xr2

s2

q3

qr sr xr qp

Figure 1: Scheme of the wastewater treatment plant A detailed description of the process dynamic model, derived from a first principles approach, is given in [37]. The overall model consists of 33 DAEs (14 of them are ODEs) and 44 variables. The value of three flowrates (qr2 , qr3 and qp) are fixed at their steady-state values corresponding to a certain nominal operational conditions. Therefore, this leaves 8 design variables for the integrated design problem, namely the volume of the aeration tanks (v1 and v2 ), the areas of the settlers (ad1 and ad2 ), the aeration factors (f k1 and f k2 ), the gain and the integral time of a PI controller. More complex formulations are possible, but 13

our objective is to illustrate how this problem of medium size is already very challenging for many GO methods. The integrated design problem is formulated as an NLP-DAEs, where the objective function to be minimized is a weighted sum of economic (φecon ) and controllability cost terms (measured here as the ISE): C = w1 · ISE + φecon = w1 · ISE + (w2 · v12 ) + (w3 · v22 ) + (w4 · ad21 ) + (w5 · ad22 ) + (w6 ·

f k12 )

+ (w7 · R

(21)

f k22 )

where the ISE is the integral square error, ISE = 0∞ e2 (t) · dt. The ISE is evaluated considering a step disturbance to the input substrate concentration, si , whose behaviour is taken from the real plant (similarly to [75]). The minimization is subject to several sets of constraints. On the one side, there are 33 model DAEs (system dynamics), acting as differential-algebraic equality constraints. On the other, there are 9 inequality constraints (see [37]) which impose limits on the residence times and biomass loads in the aeration tanks, the hydraulic capacity in the settlers, the sludge ages in the decanters, and the recycles and purge flow rates respectively. Finally, upper and lower bounds were considered for the state variables. A weighting vector wi = [103 , 2 · 10−5 , 2 · 10−5 , 1 · 10−5 , 1 · 10−5 , 12, 12] was considered for the optimization runs, which implies a similar contribution of each term in the objective function. Alternatively, the ²-constraint method could be used to obtain a pareto-optimal solution set, as in e.g. [75]. The starting point for the optimization runs had an objective function of C0 = 1.374 × 104 . Some of the bounds from the original model formulation [37] were slightly relaxed since initial runs showed that it was not possible to obtain feasible solutions for the above mentioned controller location. This problem was solved with several selected stochastic and deterministic methods which can handle black-box models. The selection has been made based on their published performance and on our own experiences considering their results for a set of GO benchmark problems. This selection included two adaptive stochastic methods, ICRS [10] and LJ [49], the GLOBAL clustering method [29] [20] and the GCLSOLVE deterministic method [40], which is a version of the DIRECT algorithm [42]. Two evolutionary computation methods were also considered, the Differential Evolution (DE) method [80], and SRES [71]. For the sake of fair comparison, we have considered Matlab implementations of all these methods, except for the case of GLOBAL, where only a Fortran implementation, difficult to translate to Matlab, was available. In order to speedup the solution of the DAEs required by the optimizers, the system’s dynamic simulation plus the objective function and constraints evaluations were implemented in compiled Fortran modules, which are callable from the solvers via simple gateways. Since most stochastic methods use 90% (or more) of the computation time in system simulations this procedure ensures good efficiency while retaining the main advantages of the Matlab environment.

14

Table 2: Summary of optimization results (weighted objective function and its components, plus relative CPU times) for best methods

C∗ φecon ISE Rel.CPU t

ICRS 1551.02 1165.14 0.3858 2.38

DE 1538.41 1130.02 0.408 4.98

GLOBAL 1544.54 1145.87 0.3986 1.00

ICRS DE LJ GLCSOLVE GLOBAL SRES

4

Objective function

SRES 1538.01 1133.20 0.405 4.74

10

3

10 0 10 4

2

1

10

2

10

3

10

CPU time,s

x 10

1.8

Objective function

1.6

Figure 2: Convergence curves

1.4

The best result was obtained with the SRES method, which arrived at C ∗ = 1538.01 after 1800 s of computation (PC Pentium III/450 MHz), followed by DE, which arrived at C1∗ = 1538.41 after 1900 s of computation. GLOBAL arrived at C ∗ = 1544.54 after 385 s of 0.8 computation (a computation time similar to that of SRES if a Matlab implementation were used). The stochastic ICRS algorithm converged to C ∗ = 1551.02 in about 900 s of 0.6 computation. A summary of these results, showing the components of the objective function 0.4 relative CPU times, are presented in Table 2. and ∗ 0.2 The other methods converged to bad objective function values, above C = 2500. The problem was also solved using a multi-start SQP approach, considering several hundred 0 0 1 2 3 4 random initial10vectors, which satisfying 10 10 were generated 10 10 the decision variables bounds. Only CPU time,s 20% of the SQP runs converged, and the corresponding results correspond to a large number of local solutions. It is very significant that despite the huge computational effort associated with all these runs, the best value found (C ∗ = 1644.65) was still far from the solutions of 1.2

15

several of the GO methods, obtained with much smaller computation times. These results illustrates well the inability of the multi-start approach to handle this type of problems.

4

10

Objective function

HyBRID

SRES

3

10 −2 10

−1

10

0

10

1

10

2

10

3

10

4

10

CPU time,s

Figure 3: Convergence curves of SRES and ICRS-SRES

The comparison of single figures of final objective function values and the associated computation times can be misleading. In order to provide with a more fair comparison of the different methods, a plot of the convergence curves (objective function values versus computation times) is presented in Figure 2 (note the log-log scales). Leaving GLOBAL apart, it can be seen that the ICRS method presented the most rapid convergence initially, but was ultimately surpassed by DE and SRES. It can also be seen that for a computation time of 200 seconds, several methods had arrived to objective function values reasonably close to the best solution. This is a common feature of stochastic GO methods, and in fact slightly relaxing the tolerance for the desired solution is exponentially efficient. Although the results obtained with GLOBAL for this problem seem to indicate that the clustering approach can be a competitive GO strategy, in fact we have found that for larger problems (i.e. more than 20 decision variables) the method becomes inefficient and usually fails. In other words, it does not scale well with problem size. In contrast, our experience indicates that SRES and DE scale up rather well. We also tested several hybrid approaches. Regarding stochastic-deterministic hybrids, like DE-SQP or SRES-SQP, they arrived to good results, but they were not able to improve the optimum reached by SRES alone. Based on the fact that the ICRS algorithm usually exhibits an initial period of fast convergence, we also experimented with ICRS-SRES hybrids, since SRES was the most efficient method in the vicinity of the best solution. Typical convergence curves for SRES and the ICRS-SRES method are shown in Figure 3, where 16

it can be clearly seen that the use of the hybrid allows faster convergence at earlier computation times. Furthermore, this method finally converged to the best known objective function for this problem (C ∗ = 1537.8). In fact, for CPU times as low as 10 seconds, rather good objective function values were already reached. Although these early solutions are of medium quality, they might be interesting for applications which do not demand very stringent convergence tolerances (for example, an interactive dynamic simulation software package with capabilities for fast near-optimal integrated design). In any case, these results are promising and we continue research in this direction.

7

6

Speedup

5

4

3

2

1 1

2

3

4

5

6

7

Number of processors

Figure 4: Speedup as a function of the number of processors for the parallel version of ICRS

As it was mentioned previously, the computation times of several of these methods can also be reduced by exploiting parallelization. Parallel implementations of the DE and ICRS methods have been tested with this problem using a local network of several low-cost platforms (Pentium III). An almost linear speed up in the number of processors was achieved (see Figure 4 for the case of ICRS), arriving to equivalent solutions, thus confirming its great potential for the solution of larger problems. The corresponding parallel implementation of the hybrid approach is being carried out.

6.3

Parameter estimation of bioprocesses

Parameter estimation problems, or inverse problems, are probably the most common type of optimization problems in engineering and science, since they are a mandatory step in any model development process. In the case of dynamic bioprocesses, the main focus of this contribution, this is also true. In this section we will mention a particular type of inverse 17

problems, originated in the area of systems biology and bioinformatics, which is especially challeging. Mendes and Kell [53] considered the parameter estimation of several rate constants of the mechanism of irreversible inhibition of HIV proteinase. This problem has a total of 20 parameters to estimate, and these authors obtained the best fit using the simulated annealing (SA) method. However, they highlighted the huge computational effort associated with this method, noting that the topic deserves more detailed study. Mendes [52] considered a larger inverse problem regarding a 3-step pathway, finding that gradient methods could not converge to the solution from any arbitrary starting vector. After comparing a number of stochastic methods, evolutionary programming (EP) was found to be the best performing algorithm, with a final refined solution which replicated the true dynamics rather well, although not perfectly. It should be noted that pseudoexperimental data was generated by simulation from a known set of parameters, so the solution of the inverse problem should lead to this original parameter set. As a drawback, the needed computation time was again excessive. More recently, Moles et al. [59] studied the same problem using several state of the art (deterministic and stochastic) methods, concluding that Evolution Strategies (ES) were the only methods to successfully solve the problem, achieving an almost perfect fit with a reduced computation time (although several hours of a Pentium III are still needed). The results of the best three methods, SRES [71] (the only one to provide a satisfactory solution), DE [80] and ICRS [10], are presented in Table 3. It can be seen that both DE and ICRS converged to local solutions which are not very far from the global one (in contrast, many other GO methods tested converged to much worse local solutions). Note that the small computation time of ICRS is due to a premature stop due to stagnation. Several on-going tests carried out by our group indicate that the use of a hybrid approach allows the reduction of the above computation time by one order of magnitude or more, while converging to essentially the same parameters as the original set. Since modeling large biological systems from functional genomic data is already an emerging area which requires proper parameter estimation methods, we envision hybrid approaches, implemented in parallel hardware, as one of the most promising solution strategies.

Table 3: Results of the GO methods (performance index J, number of evaluations and relative CPU time

J Neval CPU time, h

SRES 0.0013 28E5 1.00

18

DE 151.779 22.5E5 1.17

ICRS 183.579 16515 0.01

7

Conclusions

In this contribution, we have shown how several problems related with the optimization of dynamic bioprocesses can be efficiently solved by using stochastic and hybrid global optimization (GO) methods. Relevant examples of these problems are nonlinear dynamic optimization (open loop optimal control), simultaneous design and control optimization, and parameter estimation in nonlinear dynamic models. Empirical evidence seems to indicate that certain Evolution Strategies, like the SRES method [71], are the most promising stochastic solvers for these types of problems. Despite their main drawback (i.e., inability to guarantee global optimality) these methods have been shown to be capable of reaching very good solutions in moderate computation times, while adequately handling constraints. Parallel prototypes capable of running on existing networks of low-cost personal computers have been tested, presenting a good speedup. Several hybrid strategies have also been tested, with some of them arriving at better solutions than any other method, and with the minimum (or close to the minimum) computational effort. These results confirm the usefulness of stochastic methods in order to ensure robustness, and the superiority of the hybrid approach to obtain the best trade-off between such robustness and efficiency. Using these hybrid solvers and parallel implementations, we expect to be able to solve more challenging and larger problems with reasonable computation efforts.

Key for abbreviations ACO: ant colony optimization BVP: boundary value problem CP: complete (state and control) parameterization CVP: control vector parameterization DAEs: differential algebraic equations DE: differential evolution EC: evolutionary computation EP: evolutionary programming ES: evolution strategy GA: genetic algorithm ICRS: Integrated Controlled Random Search ISE: integral square error LJ: Luus-Jaakola MIOCP: mixed integer optimal control problem NFL: no free lunch (theorem) NLP: nonlinear programming ODEs: ordinary differential equations PDE: partial differential equation PI: proportional integral (controller) SA: simulated annealing SQP: sequential quadratic programming SRES: Stochastic Ranking Evolution Strategy 19

Acknowledgements This work was financially supported in part by the Spanish Ministry of Science and Technology (MCyT project AGL2001-2610-C02-02).

References [1] Ali M., Storey C., and T¨orn A. (1997), “Application of stochastic global optimization algorithms to practical problems,”J. Optim. Theory Appl. Vol. 95, 545–563. [2] Bailey J. E. (1998), “Mathematical modeling and analysis in biochemical engineering: Past accomplishments and future opportunities,” Biotechnol. Prog. Vol. 14, 8–20. [3] Balsa-Canto E., Alonso A. A. and Banga J. R. (1998), “Dynamic optimization of bioprocesses: deterministic and stochastic strategies,” In ACoFop IV, Automatic Control of Food and Biological Processes, G¨oteborg, Sweden, 21-23 September. [4] Balsa-Canto E., Banga J. R., Alonso A. A., and Vassiliadis V. S. (2000),“Efficient optimal control of bioprocesses using second-order information,” Ind. Eng. Chem. Res. Vol. 39, 4287–4295. [5] Banga J. R., Alonso A. A., Moles C. G. and Balsa-Canto E. (2002), “Efficient and robust numerical strategies for the optimal control of non-linear bio-processes,” In Proceedings of the Mediterranean Conference on Control and Automation (MED2002), Lisbon, Portugal. [6] Banga J. R., Alonso A. A. and Singh R. P. (1994), “Stochastic Optimal Control of Fed-Batch Bioreactors,”Presented at the AIChE Annual Meeting, San Francisco. [7] Banga J. R. , Alonso A. A., and Singh R. P. (1997),“Stochastic dynamic optimization of batch and semicontinuous bioprocesses,” Biotechnol. Prog. Vol. 13, 326–335. [8] Banga J. R. , Balsa-Canto E., Moles C. G. and Alonso A. A. (2003), “Improving food processing using modern optimization methods,” Trends Food Sci. Technol., Vol. 14 , 131–144. [9] Banga J. R. , Balsa-Canto E., Moles C. G. and Alonso A. A. (2003), “Optimization of bioreactors and bioprocesses: A review,” Proc. Indian Nat. Sci. Acad. Part A (Physical Sci.), in press. [10] Banga J. R. and Casares J. (1987),“ICRS: Application to a wastewater treatment plant model,” In IChemE Symp. Ser. No. 100, 183–192, Pergamon Press, Oxford, UK. [11] Banga J. R. , Mart´ın R. P., Gallardo J. M. and Casares J. J. (1991),“Optimization of thermal processing of conduction-heated canned foods: study of several objective functions,” J. Food Eng. Vol. 14, 25–51. 20

[12] Banga J. R. and Seider W. D. (1996),“Global optimization of chemical processes using stochastic algorithms,” In “State of the Art in Global Optimization”, C. A. Floudas and P. M. Pardalos (Eds.), Kluwer Academic Pub., Dordrecht, pages 563–583. [13] Banga J. R. , Versyck K. J. and Van Impe J. F. (2002),“Computation of optimal identification experiments for nonlinear dynamic process models: an stochastic global optimization approach,” Ind. Eng. Chem. Res. Vol. 41, 2425–2430. [14] Bansal V., Perkins J., Pistikopoulos E., Ross R. and Schijndel J. V. (2000),“Simultaneous design and control optimization under uncertainty,” Comput. Chem. Eng. Vol. 24, 261–281. [15] Barton P. I., Allgor R., Feehery W. and Gal´an S. (1998),“Dynamic optimization in a discontinuous world,” Ind. Eng. Chem. Res. Vol. 37, 966–981. [16] Barton P. I., Banga J. R. and Galan S. (2000),“Optimization of hybrid discretecontinuous dynamic systems,” Comput. Chem. Eng. Vol. 24, 2171–2191. [17] Beyer H. G. (1996),“Toward a theory of evolution strategies: Self-adaptation,” Evolutionary Computation Vol. 3, 311–347. [18] Beyer H. G. and Schwefel H. P. (2002),“Evolution strategies - a comprehensive introduction,” Natural Computing Vol. 1, 3–52. [19] Bibila T. A. and Robinson D. K. (1995),“In pursuit of the optimal fed-batch process for monoclonal antibody production,” Biotech. Progr. Vol. 11, 1–15. [20] Boender C., Kan A. R., Timmer G. and Stougie L. (1982),“A stochastic method for global optimization,” Math. Programming Vol. 22, 125–145. [21] Bonabeau E., Dorigo M. and Theraulaz G. (1999),“Swarm intelligence: from natural to artificial isystems,” Oxford University Press, New York. [22] Bonabeau E., Dorigo M. and Theraulaz G. (2000),“Inspiration for optimization from social insect behaviour,” Nature Vol. 406, 39–42. [23] Brooks S. (1958),“A discussion of random methods for seeking maxima,” Operations Res. Vol. 6, 244–251. [24] Bryson A. and Ho Y.-C. (1975),“Applied Optimal Control,” Hemisphere Publishing Corporation, New York. [25] Carrasco E. F. and Banga J. R. (1998),“A hybrid method for the optimal control of chemical processes,” IEE Conference Publication Vol. 2, 925–930 [26] Chiou J. P. and Wang F. S. (1999),“Hybrid method of evolutionary algorithms for static and dynamic optimization problems with application to a fed-batch fermentation process,” Comput. Chem. Eng. Vol. 23, 1277–1291. 21

[27] Coello C.A. (2002), “Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: a survey of the state of the art,” Comput. Methods Appl. Mech. Eng. Vol. 191, 1245–1287. [28] Corne D., Dorigo M.and Glover F. (1999),“ New Ideas in Optimization,” McGraw-Hill, New York. [29] Csendes T. (1988), “Nonlinear parameter estimation by global optimization - Efficiency and reliability,” Acta Cybernetica Vol. 8, 361–370. [30] Cuthrell J.E. and Biegler L.T. (1989),“ Simultaneous Optimization and Solution Methods for Batch Reactor Control Profiles.,” Comput. Chem. Eng. Vol. 13, 49–60. [31] Dorigo M., Maniezzo V.and Colorni A. (1996),“ The Ant System: Optimization by a Colony Cooperaing Agents,” IEEE Transac. Systems Man Cyber. B Vol. 26, 29–41. [32] Esposito W. R. and Floudas C. A. (2000), “Deterministic Global Optimization in Nonlinear Optimal Control Problems,” J. Global Optim. Vol. 17, 97–126. [33] Esposito W. R. and Floudas C. A. (2000), Global Optimization for the Parameter Estimation of Differential-Algebraic Systems,” Ind. Eng. Chem. Res. Vol. 39, 1291– 1310. [34] Fogel L.J., Owens A.J. and Walsh M.J. (1966),“ Artificial intelligence through simulated evolution,” Wiley, New York. [35] Fraga E.S., Hagemann J., Estrada-Villagrana A.and Bogle I.D.L. (2000),“Incorporation of Dynamic Behaviour in a Automated Process Synthesis System,” Comput. Chem. Eng. Vol. 24, 189–205. [36] Goldberg D.E. (1989),“Genetic Algorithms in Search, Optimization and Machine Learning,” Addison Wesley Longman, London. [37] Gutirrez G. and Vega P. (2000),“Integrated Design of Activated Sludge Process Taking in to Account the Closed Loop Controllability,” Proc. ESCAPE-10, 63–69. [38] Heckler R. and Schwefel H.-P. (1978),“Superimposing direct search methods for parameter optimization onto dynamic simulation models. In Proc. Winter Simulation Conf. Vol. 1, 173–181, Miami Beach FL. [39] Holland J.H. (1992),“Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence,” MIT Press, Cambridge, MA. [40] Holmstr¨om K. (1999),“The TOMLAB Optimization Environment in Matlab,” Adv. Model. Optim. Vol. 1, 47. [41] Jayaraman V.K., Kulkarni B.D., Gupta K., Rajesh J. and Kusumaker H.S. (2001),“Dynamic Optimization of Fed-Batch Bioreactors Using the Ant Algorithm,” Biotech. Progr. Vol. 17, 81–88. 22

[42] Jones D.R. (2001),“DIRECT algorithm,” In Encyclopedia of Optimization, C.A. Floudas and P.M. Pardalos eds., Kluwer Academic Publishers, Dordrecht. [43] Kirkpatrick S., Gellatt C.D.and Vecchi M.P. (1983),“Optimization by Simulated Annealing,” Science Vol. 220, 671–680. [44] Laarhoven P.J.M. and Aarts E.H.L. (1987),“Simulated annealing: theory and applications,” Reidel Publishing Company, Dordrecht. [45] Lee J. H. (1999),“Comparison of various optimization approaches for fed-batch ethanol production,” Appl. Biochem. Biotech. Vol. 81, 91–105. [46] Lopez Cruz I.L., van Willigenburg L.G. and van Straten G. (2000),“Evolutionary algorithms for optimal control of chemical processes,” In Proceedings of the IASTED International Conference on Control Applications, Cancun (Mexico). [47] Luus R. (1993),“ Optimization of Fed-Batch Fermentors by Iterative Dynamic Programming,” Biotechnol. Bioeng. Vol. 41, 599–602. [48] Luus R. and Hennessy D. (1999),“Optimization of Fed-Batch Reactors by the LuusJakola Optimization Procedure,” Ind. Eng. Chem. Res. Vol. 38, 1948–1955. [49] Luus R. and Jaakola T.H.I. (1973),“Optimisation of Non-Linear Functions Subject to Equality Constraints,” IEC Proc. Des. Dev. Vol. 12, 380–383. [50] Matsuura K., Shiba H., Nunokawa Y.and Shimizu H. (1993),“Calculation of optimal strategies for fermentation processes by genetic algorithm,” J. Soc. Ferment. Bioeng. Vol. 71, 171–178. [51] Matyas J. (1965),“Random Optimization,” Automation Remote Control Vol. 26, 246– 253. [52] Mendes P. (2001),“Modeling Large Biological Systems from Functional Genomic Data: Parameter Estimation,” In Foundations of Systems Biology, H. Kitano (ed.), 163–186, MIT Press, Cambridge, MA. [53] Mendes P. and Kell D.B. (1998), “Non-Linear Optimization of Biochemical Pathways: Applications to Metabolic Engineering and Parameter Estimation,” Bioinformatics Vol. 14, 869–883. [54] Michalewicz Z., Dasgupta D., Le Riche R.G. and Schoenauer M. (1996),“Evolutionary algorithms for constrained engineering problems,” Comput. Ind. Eng. Vol. 30, 851–870. [55] Michalewicz Z. and Fogel D.B. (2000),“How to solve it: modern heuristics,” Springer, Berlin. [56] Michalewicz Z., Janikow C.Z. and Krawczyk J.B. (1992),“A modified genetic algorithm for optimal control problems,” Comput. Math. with Appl. Vol. 23, 83–94. [57] Michalewicz Z. and Schoenauer M. (1996),“Evolutionary Algorithms for Constrained Parameter Optimization Problems,” Evolutionary Comput. Vol. 4, 1–32. 23

[58] Moles C. G., Gutierrez G., Alonso A.A. and Banga J. R. (2001),“ Integrated Process Design and Control via Global Optimization: a Wastewater Treatment Plant Case Study,” In European Control Conference (ECC) 2001, Porto (Portugal). [59] Moles C.G., Mendes P. and Banga J. R. (2001),“Global optimization of biochemical pathways: the parameter estimation problem,” In Proceedings of the Third International Conference in Systems Biology ICSB 2002, Stockholm. [60] Morari M. and Perkins J.D. (1994),“Design for Operations,” In Foundations of Computer-Aided Process Design (FOCAPD). [61] Na J. G., Chang Y. K., Chung B. H. and Lim H. C. (2002),“Adaptive optimization of fed-batch culture of yeast by using genetic algorithms,” Bioproc. Biosyst. Eng., Vol. 24, 299–308. [62] Papamichail I. and Adjiman C. S. (2002),“A Rigorous Global Optimization Algorithm for Problems with Ordinary Differential Equations,” J. Global Optimization, Vol. 24, 1–33. [63] Pistikopoulos E.N., Ross R. and Schijndel J.M.G. (1999),“ Towards the Integration of Process Design, Process Control and Process Operability. Current Status and Future Trends,” In Foundations of Computer-Aided Process Design (FOCAPD), Snowmass, CO. [64] Rastrigin L.A. and Rubinstein Y. (1969),“The Comparison of Random Search and Stochastic Approximation while Solving the Problem of Optimization,” Automatic Control, Vol. 2, 23–29. [65] Rechenberg I. (2000),“Case studies in evolutionary experimentation and computation,” Comput. Methods Appl. Mech. Eng., Vol. 186, 125. [66] Rinnooy-Kan A. H. G. and Timmer G. T. (1987),“Stochastic Global Optimization Methods. Part I: Clustering Methods,” Math. Progr., Vol. 39, 27. [67] Rodriguez-Acosta F., Regalado C. M. and Torres N. V. (1999), “Non-linear optimization of biotechnological processes by stochastic algorithms: Application to the maximization of the production rate of ethanol, glycerol and carbohydrates by Saccharomyces cerevisiae,” J. Biotechnol., Vol. 68, 15–28. [68] Ronen M., Shabtai Y.and Guterman H. (2002),“Optimization of feeding profile for a fed-batch bioreactor by an evolutionary algorithm,” J. Biotechnol., Vol. 97, 253–263. [69] Roubos J.A., Gooijer C.D. de, Van Straten G. and Van Boxtel A.J.B. (1997), “Comparison of Optimization Methods for Fed-Batch Cultures of Hybridoma Cells,” Bioproc. Eng., Vol. 17, 99–102. [70] Roubos J.A., Van Straten G. and Van Boxtel A.J.B. (1999), “An Evolutionary Strategy for Fed-Batch Bioreactor Optimization; Concepts and Performance,” J. Biotechnol., Vol. 67, 173–187. 24

[71] Runarsson T.P. and Yao X. (2000), “Stochastic Ranking for Constrained Evolutionary Optimization,” IEEE Transac. Evol. Comp., Vol. 4, 284–294. [72] San K.Y. and Stephanopoulos G. (1984), “A note on the optimality criteria for maximum biomass production in a fed-batch fermentor.,” Biotechnol. Bioeng., Vol. 26, 1261–1264. [73] Schittkowski K. (2003),“Numerical data fitting in dynamical systems: a practical introduction with applications and software,” Kluwer Academic Publishers, Dordrecht. [74] Schwefel H.-P. (1995),“Evolution and optimum seeking,” Wiley, New York. [75] Schweiger C.A. and Floudas A. (1997), “Interaction of Design and Control: Optimization with Dyamic Models,” In Optimal Control Theory, Algorithms and Applications, (Hager, W.W. and Pardalos, P.M., eds.), Kluwer Academic Publishers, Dordrecht. [76] Seywald H., Kumar R.R. and Deshpande S.M (1995), “Genetic Algorithm Approach for Optimal Control Problems with Linearly Appearing Controls,” J. Guid. Contr. Dyn., Vol. 18, 177–182. [77] Shimizu K. (1996), “A Tutorial Review on Bioprocess Systems Engineering,” Comput. Chem. Eng., Vol. 20, 915–941. [78] Simutis R. and Lubbert A. (1997), “A comparative study on random search algorithms for biotechnical process optimization,” J. Biotechnol., Vol. 52, 245–256. [79] Singer A. B., Bok J. K.and Barton P. I. (2001), “Convex Underestimators for Variational and Optimal Control Problems,” Computer Aided Chem. Eng., Vol. 9, 767–772. [80] Storn R. and Price K. (1997), “Differential Evolution - a Simple and Efficient Heuristic for Global Optimization over Continuous Spaces,” J. Global. Optim., Vol. 11, 341–359. [81] T¨orn A., Ali M. and Viitanen S. (1999), “Stochastic Global Optimization: Problem Classes and Solution Techniques,” J. Global Optim., Vol. 14, 437. [82] T¨orn A. A. (1973),“Global optimization as a combination of global and local search,” In Proceedings of Computer Simulation Versus Analytical Solutions for Business and Economic Models, Gothenburg. [83] Vassiliadis V. S. (1993), “ Computational Solution of Dynamic Optimization Problems with General Differential-Algebraic Constraints,” PhD Thesis, Imperial College, University of London. [84] Vassiliadis V. S., Balsa-Canto E.and Banga J. R. (1999), “Second Order Sensitivities of General Dynamic Systems with Application to Optimal Control Problems.,” Chem. Eng. Sci., Vol. 54, 3851–3860. [85] Versyck K., Claes J. and Van Impe J. (1997), “Practical identification of unstructured growth kinetics by application of optimal experimental design.,” Biotechnol. Prog., Vol. 13, 524–531. 25

[86] Wang F.S. and Chiou J.P. (1997), “Optimal Control and Optimal Time Location Problems of Differential-Algebraic Systems by Differential Evolution,” Ind. Eng. Chem. Res., Vol. 36, 5348–5357. [87] Wolpert D.H. and Macready W.G. (1997), “No Free Lunch Theorems for Optimization,” IEEE Transac. Evol. Comp., Vol. 1, 67–82. [88] Yamashita Y. and Shima M. (1997), “Numerical Computational Method Using Genetic Algorithm for the Optimal Control Problem with Terminal Constraints and Free Parameters,” Nonlinear Anal. Theory Methods Appl. Vol. 30, 2285–2290. [89] Zabinsky Z.B. and Smith R.L. (1992), “Pure Adaptive Search in Global Optimization,” Math. Progr., Vol. 53, 323–338.

26