Global Optimization for Integrated Design and Control of

Moreover, these results show that surrogate-based methods may reduce ... past decade from many researches in computer-aided process engineering.2,3 In...
4 downloads 0 Views 152KB Size
Download PDF
9148

Ind. Eng. Chem. Res. 2007, 46, 9148-9157

Global Optimization for Integrated Design and Control of Computationally Expensive Process Models Jose A. Egea Process Engineering Group, IIM-CSIC, Eduardo Cabello 6, 36208 Vigo, Spain

Dirk Vries Systems and Control Group, Wageningen UniVersity, P. O. Box 17, 6700 AA Wageningen, The Netherlands

Antonio A. Alonso Process Engineering Group, IIM-CSIC, Eduardo Cabello 6, 36208 Vigo, Spain

Julio R. Banga* Process Engineering Group, IIM-CSIC, Eduardo Cabello 6, 36208 Vigo, Spain

The problem of integrated design and control optimization of process plants is discussed in this paper. We consider it as a nonlinear programming problem subject to differential-algebraic constraints. This class of problems is frequently multimodal and “costly” (i.e., computationally expensive to evaluate). Thus, on the one hand, local optimization techniques usually fail to locate the global solution, and, on the second hand, most global optimization methods require many simulations of the model, resulting in unaffordable computation times. As an alternative, one may consider global optimization methods which employ surrogate-based approaches to reduce computation times and which require no knowledge of the underlying problem structure. A challenging wastewater treatment plant (WWTP) benchmark model1 is used here to evaluate the performance of these techniques. Numerical experiments with different optimization solvers indicate that the proposed benchmark optimization problem is indeed multimodal, and that via global optimization we can achieve an improvement of the controllers’ performance compared to the best tuned controllers’ settings available in the literature. Moreover, these results show that surrogate-based methods may reduce computation times while ensuring convergence to the best known solutions. Introduction Many optimization problems which arise during the design and/or operation of chemical and biochemical processes are frequently multimodal. Thus, it is not surprising that global optimization (GO) of such processes have received much attention during the past decade from many researches in computer-aided process engineering.2,3 In fact, many advances have been made regarding both deterministic and stochastic methods during recent years.4 However, the current state of the art is far from satisfactory, especially when we consider the global optimization of complex process models. These models are typically complex due to their nonlinear dynamic behavior and large number of states. Therefore they are typically expensive to evaluate (i.e., each simulation can take minutes or even hours in CPU computation time on an ordinary PC), and for a number of reasons they can only be treated as blackbox models. Thus, several authors5-8 have recently proposed the use of so-called surrogate models which are (computationally) cheaper to evaluate. In the area of process systems engineering, Moles et al.9 showed how stochastic global optimization methods can be applied for simultaneous design and control of process plants of medium complexity. In this contribution, our aim is to apply a similar approach to more complex (and costly) models. To keep the computational burden acceptable, the capabilities of * To whom correspondence should be addressed. Tel.: +34 986 214 473. Fax: +34 986 292 762. E-mail: [email protected].

recent global surrogate model based optimization methods are evaluated here. Such evaluation is performed considering a challenging benchmark case study: the integrated design and control of a wastewater treatment plant (WWTP) for nitrogen removal, such as the one developed by the European Cooperation in Science and Technology (COST) 624 work group.10 This paper is structured as follows: in the next section, the general statement of the integrated design problem is presented. Next, we briefly review the state of the art regarding the global optimization of such problems, and we present the methods selected for our research. In the following section, the WWTP case is outlined. Finally, sections with results, discussion, and conclusions are provided. Integrated Design: Problem Statement The general statement of the simultaneous (integrated) design and control problem takes into account the process and control superstructures indicating the different design alternatives.11,12 This general approach results in mixed integer optimal control problems. In this work, we consider a simpler, yet nontrivial, subproblem, where it is assumed that the flowsheet (i.e., plant configuration) is given. It should be noted that this subproblem is challenging enough to serve as a case study for the comparison of different (including surrogate-based) global optimization methods, which is the main objective of this work.

10.1021/ie0705094 CCC: $37.00 © 2007 American Chemical Society Published on Web 11/28/2007

Ind. Eng. Chem. Res., Vol. 46, No. 26, 2007 9149

The aim is to simultaneously find the static variables of the process design, the operating conditions, and the controllers’ parameters 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. We state our problem as follows: Find V to minimize

∑wiφi

(1)

f(x˘ ,x,V) ) 0

(2)

x(t0) ) x0

(3)

h(x,V) ) 0

(4)

g(x,V) e 0

(5)

V L e V e VU

(6)

C) subject to

where V is the vector of decision variables (e.g., design variables, operating conditions, parameters of controllers, setpoints, etc.); C is the cost (objective function) to minimize (normally a weighted combination of capital, operation, and controllability costs, φi); f is a functional for the system dynamics (i.e., the nonlinear process model); x is the vector of the states; t0 is the initial time for the integration of the ODE’s (and, consequently, x0 is the vector of the states at that initial time); h and g are possible equality and inequality path and/or point constraint functions which express additional requirements for the process performance; and, finally, VL and VU are the lower and upper bounds for the decision variables. The dependence of φi upon the decision variables, V, is defined by the problem formulation. In some cases this dependence is simple and straightforward (i.e., when minimizing the cost of a chemical process one of the φi can be equal to a reactor size, which may also be a decision variable), whereas in others there might not be an explicit expression for this dependence (i.e., φi can be the integral square error (ISE) of a controller which, in general, does not depend explicitly on the decision variables). The formulation above is that of a nonlinear programming problem (NLP) with differential-algebraic (DAEs) constraints. Due to the nonlinear and constrained nature of the system dynamics, these problems are very often multimodal. Further, it is known that using standard controllability measures such as the ISE in the objective function often causes nonconvexity.11 Therefore, if this NLP-DAEs is solved via standard local NLP methods, such as sequential quadratic programming (SQP), it is very likely that the solution found will be nonglobal. To avoid this, GO should be used. Optimization Methods Need of Global Optimization. In principle, model-based optimization can be successfully used to improve the design and/or operation of single units or full process plants. Typically, most of the problems in process engineering applications are highly constrained and exhibit nonlinear dynamics. These properties often result in nonconvexity and multimodality. Furthermore, in many complex process models some kind of noise and/or discontinuities (due either to numerical methods or to intrinsic properties of the model) is present. Therefore, there is a great need of robust global optimization solvers which

can (i) locate the vicinity of the global solution in a reasonable number of iterations, (ii) handle noise and/or discontinuities, and (iii) use some sort of approximation or reduction of the original process model, in order to keep the computational effort acceptable. In general, (iterative) gradient-based local methods for constrained NLP problems are very efficient, but they can only handle differentiable objective and constraint functions based on a set of continuous variables.13 Additionally, only convergence to nonglobal solutions is guaranteed. Therefore, one must use the so-called GO methods9,14,15 in order to provide an approximation of the global optimum. GO methods can be roughly classified as being deterministic16,17 or stochastic.18 Deterministic GO methods ensure convergence to the global optimum for certain problems, although no algorithm can solve general GO problems with certainty in a finite time.18 For these methods, the computational effort usually increases exponentially with the problem size. Further, they have requirements (e.g., smoothness, differentiability) which are rarely met in realistic applications. Stochastic GO methods can usually find solutions in the vicinity of the global solution in relatively short computation times compared to deterministic GO methods. Note that, with these stochastic methods, the convergence to global optimality (in finite time) is not guaranteed either, but many empirical studies have shown that they are often the best methods for many classes of problems. Other advantages of these methods are that they are easy to implement and they can treat the objective function as a black box (i.e., a simple connection between inputs and outputs, with no derivative information needed). This feature is especially appealing in the case of complex dynamic systems where the objective function is the result of, e.g., a simulation provided by a third party software with restricted access for the user (a common case in real industrial applications). Surrogate-Based Global Optimization. In general, all these GO approaches require a significant number of evaluations of the objective and constraint functions. In the case of realistic problems, these models are costly to evaluate, posing a major challenge to the application of GO methods. In recent years, a number of approaches have been proposed to obtain surrogate models which are cheaper to evaluate than the original ones and which imitate the original model based on a reduced number of sampled points (simulations). Hence, these so-called surrogate model based solvers try to approximate the original model over a region by a model that is cheaper to evaluate. Provided the surrogate model is accurate enough, the computation times can be dramatically reduced. In addition, surrogate models go beyond simply reducing computation time (e.g., kriging-based methods even provide statistical information about decision variables). Surrogate model based methods can be classified into two groups:19 non-interpolating (e.g., quadratic polynomials and other regression models) and interpolating methods (e.g., basis functions and kriging). At the same time, both methods can be one-stage or two-stage methods. Two-stage methods fit first a response surface using sample points from the real model and then optimize an auxiliary function based on the fitted surface. A potential disadvantage of these methods is that the initial surface may not accurately fit the real model, which can provoke the optimization to stop prematurely or converge to nonglobal solutions. On the other hand, one-stage methods evaluate hypotheses about the location of the optimum. This is done by examining the best-fitting response surface passing through the

9150

Ind. Eng. Chem. Res., Vol. 46, No. 26, 2007

Figure 1. Benchmark layout.

observed data and another point in which the optimum is presumed to be located. Each hypothesis is evaluated and the surface is constructed by evaluating the new points where this credibility is maximum. Jones19 presented a taxonomy of these methods with a complete overview of the different approaches in which it is concluded that interpolating methods are more suitable than noninterpolating methods to locate the global optimum, specifically mentioning kriging and basis functions. However, noninterpolating quadratic fitting methods have been widely used in response surfaces based engineering design, as surveyed in Simpson et al.20 Thus, to compare the most used techniques (i.e., quadratic methods) and the most recent and promising ones (i.e., radial basis functions and/or kriging) the different types of response surface methodologies have been tested and evaluated, considering a realistic and challenging process engineering problem. Selected Optimization Methods. Regarding surrogate-based GO, we have considered three recent solvers which are Matlab implementations of the strategies mentioned above: (a) rbfSolve. Included in the Tomlab toolbox,21 it solves costly global optimization problems using a radial basis functions (RBF) interpolation algorithm.5,6 It fits a response surface (based on splines) to data collected by evaluating the objective function at some points and then applies a global optimization algorithm (a DIRECT22 implementation called glcFast from the same toolbox) and a local algorithm (npsol23) using different initial points over that surrogate model. The first set of points to create the response surface may be given by the user or selected by the algorithm based on different strategies. (b) ego. Also included in the Tomlab toolbox,21 it solves costly global optimization problems using the Efficient Global Optimization (EGO) algorithm,7 based on kriging interpolation. The idea of the EGO algorithm is to first fit a response surface to data collected by evaluating the objective function at a few points. Then, EGO balances between finding the minimum of the surface and improving the approximation by sampling where the prediction error may be large. (c) snobfit. (Stable noisy optimization by branch and fit) builds local quadratic models from selected sampled points and combines local and global search over these models. It is specially designed for costly black-box functions where some noise can be present.8 To critically evaluate the performance of these three surrogate-based strategies, we have also considered selected local and global solvers which only rely on evaluations of the original (costly) model:

(1) fminsearch. This is a direct search local method included in the Matlab Optimization Toolbox24 based on the simplex method by Nelder and Mead. Although generally less efficient than gradient-based methods, the simplex method may be more robust if the problem is highly discontinuous or contains noise. (2) fmincon. This local solver, which is also part of the Matlab Optimization Toolbox,24 finds a local minimum of a constrained multivariable function by means of a SQP algorithm. The method uses numerical or, if available, analytical gradients. (3) NOMADm. (Nonlinear optimization for mixed variables and derivatives-Matlab) is a Matlab code that runs various generalized pattern search (GPS) algorithms to solve nonlinear and mixed variable optimization problems.25 (4) Differential Evolution (DE). This method26 is a heuristic population-based stochastic approach to global optimization for minimizing possibly nonlinear and nondifferentiable continuous space functions. (5) Stochastic Ranking with Evolution Strategy (SRES). This method uses a (µ, λ) evolution strategy combined with an approach to balance objective and penalty functions stochastically. This feature makes it specially appealing for the case of constrained problems.27 (6) GLOBAL. This is a hybrid GO method28 based on the algorithm by Boender et al.29 It uses iteratively a random search followed by a local search routine. It carries out a clustering phase regularly. Next, two different local search procedures can be selected for a second step. The first one is an algorithm of quasi-Newton type. The second, more appropriate for problems with discontinuous objective functions or derivatives, is a robust random search method (UNIRANDI) by Ja¨rvi.30 Case Study: Optimization of Controllers Settings in a WWT Plant A number of control strategies have been proposed to meet the strict standards that WWTPs must comply with, while also trying to reduce costs.31-34 Relevant examples from the recent literature of attempts to optimize the controllers of these plants are as follows: (1) ad hoc extensive simulation studies35,36 (strictly speaking these may not be called optimizations, because there is no evidence that a locally or globally best solution is found); (2) dynamic optimization design strategies using local gradient-based37 or evolutionary38 optimization methods, often based on simplified models and without use of any control strategy; (3) global optimization methods for simultaneously optimizing operation and design;9,39 (4) an

Ind. Eng. Chem. Res., Vol. 46, No. 26, 2007 9151

integrated approach for the optimization of control strategies, where a small selection of global and local optimization methods was used.40 Evaluation of these and similar strategies, either in practice or by simulation, is a real problem due to the lack of a standard with respect to evaluation criteria, process complexity, and large variations in plant configuration. To enhance the development and acceptance of new control strategies, the International Water Association (IWA) Task Group on Respirometry, together with the European COST work group, proposed a standard simulation benchmarking methodology for evaluating the performance of activated sludge plants. The COST 624 work group defines the benchmark as “a protocol to obtain a measure of performance of control strategies for actiVated sludge plants based on numerical, realistic simulations of the controlled plant”. According to this definition, the benchmark consists of a description of the plant layout, a simulation model, and definitions of (controller) performance criteria. The layout of this benchmark plant combines nitrification with pre-denitrification by a five-compartment reactor with an anoxic zone (see Figure 1). A secondary settler separates the microbial culture from the liquid being treated. A basic control strategy consisting of two PI controllers is proposed to test the benchmark. Its aim is to control the dissolved oxygen level in the final compartment of the reactor (AS unit 5) by manipulation of the oxygen transfer and to control the nitrate level in the last anoxic compartment (AS unit 2) by manipulating the internal recycle flow rate.41 In this work, a Simulink implementation of the benchmark model by Jeppsson was used for the simulations.1,10,41 Each function evaluation consists of an initialization period of 100 days to achieve steady state, followed by a period of 14 days of dry weather and a third period of 14 days of rainy weather. Calculations of the controller performance criterion are based on data from the last 7 rain days. As far as we know, global optimization methods have not been used for the integrated design and control of WWTPs models of this complexity. Further, since each simulation of this benchmark model takes a significant time on a standard PC (about 90 s in a PC-PIV 2,4 GHz), it is an illustrating example to evaluate the surrogate-based strategies mentioned previously. The control performance is tested by using the ISE as a controller performance criterion. Both the nitrate level and oxygen level controllers (further referred as N- and O-controller, respectively) are optimized with respect to their controller parameters, that is, the gain K, integral time constant τi, and antiwindup time constant τt. The problem is formulated as follows:

min J(V,tf) ) WTISE

(7)

controllers; f describes the system dynamics, and h and g are possible equality and inequality trajectory and/or endpoint constraint functions which express additional requirements for the process performance. The weighting vector WT, the ISE, and the decision parameter vector are as follows:

WT ) [w1 w2] )

1000 1 [1001 1001]

[ ]

ISE ISE ) ISEO N ISE(.) )

[]

(13) (14)

∫t t (τ)2(.) dτ

(15)

K(O) τi(O) τt(O) ) K (N) τi(N) τt(N)

(16)

[ ]

V V ) VO N

f

0

The weighting vector is chosen such that the ISE(O) equals the ISE(N) part when using the benchmark default settings provided by the COST project.10 The system dynamics is described by algebraic mass balance equations, ordinary differential equations for the biological processes in the bioreactors as defined by the ASM1-model,42 and the double-exponential settling velocity function43 as a fair presentation of the settling process, with x ∈ ∇13×1 the state vector, p the system parameters, and d the influent disturbance. A detailed model of the benchmark process is provided by the IWA Task Group on Benchmarking of Control Strategies for WWTPs.44 Due to the complexity of the system dynamics, the problem cannot be solved analytically, and the optimal values for the decision variables V ∈ ∇6×1 (i.e., the PI controller parameters) have to be retrieved by optimization techniques or open loop controller tuning. However, the problem remains that most tuning techniques are intended for linear systems so that model approximation by linearization around different operating points would be inevitable. In the particular case considered here, no further constraints apart from the systems dynamics are laid on the problem. Boundaries on the decision variables (VL and VU) are chosen such that the process dynamics would not show (exceptional) unstable behavior: VL ) [100.00 7.00 × 10-4 1.00 × 10-4 100 1.00 × 10-2 1.00 × 10-4]T

(17) VU ) [1000.00 7.00 × 10-1 7.00 × 10-1 50 000 1.00 × 100 7.00 × 10-2]T

s.t.

(18) The objective function values are normalized with respect to the performance obtained with the tuned controller settings provided by the COST project, which are

x˘ ) f(x,V,p,d)

(8)

x(t0) ) x0

(9)

h(x,V,p) ) 0

(10)

VCOST )

g(x,V,p) e 0

(11)

[500.00 1.00 × 10-2 2.00 × 10-4 15 000 5.00 × 10-2 3.00 × 10-2]T

V L e V e VU

(12)

Results and Discussion

(19)

where WT ∈ ∇1×2 contains the weighting coefficients, ISE ∈ ∇2×1 contains the integral squared errors of the two PI

All runs were carried out on a Pentium-IV 2.4 GHz computer. To illustrate the nonconvexity of the problem, we investigated

9152

Ind. Eng. Chem. Res., Vol. 46, No. 26, 2007

Figure 2. Histogram of solutions obtained by a local method (fmincon) using 15 different initial decision vectors.

presented in Figure 2. It is interesting that, despite the large computational effort (about 8 h of CPU time per run), only two of the runs were able to improve the controllers nominal settings (which we have normalized to J ) 1.0), and the best of them, with an objective function value of J ) 0.64, is worse than the best solutions provided by the global solvers based on surrogate models, presented below. As expected, the global solvers performed much better due to their capabilities to escape from local solutions. Typical convergence curves (objective function value versus CPU time) from the same initial point, randomly chosen within the bounds and denoted by x0 with a corresponding objective function value of J ) 36.14, are presented in Figure 3 (note the log scale for J). x0 ) [750.62 5.07 × 10-1 4.62 × 10-1 27 831 9.32 × 10-2 2.47 × 10-2] Figure 3. Convergence curves for the different optimization algorithms. Table 1. Solutions Provided by Stochastic Solvers Using the Same Initial Point in 10 Runs solver

best solution

worst solution

mean

SRES DE snobfit rbfSolve GLOBAL

0.9026 0.6434 0.5496 0.5344 0.8913

2.2191 1.1670 1.8886 0.8958 3.5262

1.6514 0.8407 0.8593 0.6037 1.7004

the performance and robustness of the local solvers fminsearch and fmincon by applying a multistart procedure. Hence, initial values for the decision variables were randomly selected within their valid range, and optimization runs are started from there. Both solvers converged to nonglobal solutions in all the runs, confirming the need of using global optimization methods. The histogram of the solutions obtained by fmincon in 15 runs is

(20) The best final solutions were consistently found by rbfSolve, which also showed a good convergence rate. In fact, this method provided a good solution in less than 1 h of computation time. snobfit, another method based on response surfaces, achieved also an excellent final result as well as a very good convergence rate. Other solvers such as DE, NOMADm, SRES, and GLOBAL finally reached good values of the performance index (better than the one obtained in the literature), but not as good as those provided by rbfSolve and snobfit, although their convergence rates (specially in the case of NOMADm) were quite good. The other surrogate-based method (ego) showed a rather good convergence rate but stopped prematurely in almost all runs. Thus, it seems that surrogate models based on RBFs perform better than kriging for this problem. The convergence curve for a local solver such as fminsearch presented in Figure 3 shows how local strategies are not suitable for GO problems. In this case fminsearch obtained a result far from the best known values for this problem.

Ind. Eng. Chem. Res., Vol. 46, No. 26, 2007 9153

Figure 4. Convergence curves for (a) SRES, (b) DE, (c) snobfit, (d) rbfSolve, and (e) GLOBAL using x0 as starting point (except for GLOBAL that does not use initial given points).

For the case of stochastic solvers (e.g., SRES, DE, and GLOBAL) and surrogate-based algorithms (e.g., rbfSolve or snobfit) with stochastic initial sampling set, 10 runs were performed using x0 as starting point (except for GLOBAL, which does not use any initial point as input) in order to check the variability of each algorithm. Table 1 shows the best and the worst solution obtained by each solver, as

well as the mean value of all runs. Parts a-e of Figure 4 show the convergence curves for each solver. fminsearch, fmincon, NOMADm and snobfit were applied using their default parameters. For the rest of the solvers some tuning was performed as follows (we maintain the same notation as in their original codes):

9154

Ind. Eng. Chem. Res., Vol. 46, No. 26, 2007

Figure 5. Evolution of ISE(O) for nominal and optimized (by rbfSolve) controller parameters. Table 2. Best Optimization Results Obtained by the Different Optimization Algorithms (Starting from the Same Initial Point) PI controller parameter solver NOMADm fmincon fminsearch DE GLOBAL SRES rbfSolve snobfit ego a

K(O) 750.62 661.39 819.08 357.51 345.05 606.97 375.31 548.01 444.59

τi(O) 10-3

6.91 × 3.98 × 10-2 4.87 × 10-1 1.73 × 10-3 8.88 × 10-4 3.74 × 10-3 7.00 × 10-4 a 1.40 × 10-3 7.00 × 10-4 a

τt(O)

K(N)

10-1

6.30 × 6.59 × 10-1 4.51 × 10-1 1.00 × 10-4 2.67 × 10-1 3.21 × 10-1 1.00 × 10-4 4.55 × 10-1 1.00 × 10-4

a

a a

21 175 28 478 26 508 17 876 25 181 26 426 21 018 18 278 6 129

τi(N)

τt(N)

10-2

3.07 × 3.74 × 10-1 9.47 × 10-2 2.00 × 10-2 8.82 × 10-2 5.94 × 10-2 2.71 × 10-2 2.38 × 10-2 1.00 × 10-2

J

10-3

a

2.97 × 5.85 × 10-2 2.55 × 10-2 1.00 × 10-4 1.85 × 10-2 3.90 × 10-2 1.46 × 10-2 4.92 × 10-2 7.00 × 10-2

a

a

0.8170 8.6229 30.1660 0.6434 0.8913 0.9026 0.5333 0.5496 1.4722

Active bound constraint.

• DE: Crossover probability (CR) ) 0.9, Binomial crossover: Crossover is applied to every decision variable with probability defined by CR (strategy ) 2). • rbfSolVe and ego: Initial sampling points to fit the surface chosen by the user. In every case a latin hypercube sampling plus x0 was chosen (GO.Percent ) 1000), Change the seed for generating random numbers in every optimization. This allows one to generate different initial surfaces when applying the latin hypercube sampling (RandState ) -10). • SRES Population size (λ) ) 150, Number of parents (µ) ) 15. • GLOBAL Numbers of points sampled in every iteration (nsampl) ) 60, Number of points added for the clustering phase in every iteration (nsel) ) 20.

The best solutions found by each solver are presented in Table 2, where it can be seen that several solvers were able to improve the controllers nominal settings. Simulations using the best solution obtained by rbfSolve and snobfit, denoted by Vrbf and Vsnobf, respectively, were performed to compare other criteria such as effluent quality and pollutant violations with those using the nominal parameter values, denoted by VCOST, provided by the COST project. Results are presented in Table 3, indicating that Vrbf and Vsnobf provide better results than VCOST for those criteria (except in the case of pumping energy). Hence, although we only minimized an ISE-based objective function in this work, the optimized plant also performed better in almost all the other relevant indexes. In the future we will consider multiobjective global optimization to check whether those criteria can be further improved. The ISE evolution for the two controllers is shown in Figures 5 and 6 (for dissolved oxygen and nitrates levels, respectively). It is clear how the optimized (by rbfSolve) controller’s parameters with respect to the default parameters provided by the COST project mainly affect the ISE of the N-

Ind. Eng. Chem. Res., Vol. 46, No. 26, 2007 9155

Figure 6. Evolution of ISE(N) for nominal and optimized (by rbfSolve) controller parameters.

Figure 7. Response for oxygen concentration in tank 5 (solution provided by rbfSolve). Table 3. Values of WWTP Performance Criteria for the Reference and the Optimized Decision Vectors criterion effluent quality (EQ) index violations max. total N-level violations max. total ammonia N-level aeration energy av pumping energy ISE(O) ISE(N)

results with VCOST

results with Vsnobf

results with Vrbf

9032 11.16 25.89 7173 1919 1.223 × 10-5 0.8335

8984 9.97 25.15 7165 1984 2.404 × 10-5 0.4426

8980 9.97 25.15 7164 1985 1.285 × 10-5 0.4378

controller. Indeed, ISE(O) remains almost the same using both nominal parameters and optimized parameters. The major influence in the normalized ISE is due to the improvement in ISE(N).

units (kg poll‚units)/day % of operating time % of operating time kWh/day kWh/day (mg COD/L)2‚day (mg N/L)2‚day

In Figures 7 and 8 the dynamic behavior of the controlled variables (i.e., total N concentration in tank 2 and dissolved oxygen in tank 5) using the optimized parameters provided by rbfSolve are presented. As shown, the control strategy is

9156

Ind. Eng. Chem. Res., Vol. 46, No. 26, 2007

Figure 8. Response for N concentration in tank 2 (solution provided by rbfSolve).

effective to maintain the concentration of oxygen in tank 5, whereas the concentration of N in tank 2 presents rather large fluctuations around the set point. This can be explained because the oxygen controller is assumed to be ideal with no delay or noise, while the N controller has a time delay of 10 min with white, normally distributed zero-mean noise (with standard deviation of 0.1 g m-3). Thus controlling the N concentration in tank 2 is a more difficult task than controlling the oxygen concentration in tank 5. The implementation of the two PI controllers is a simple control strategy, and the benchmark is built for different designers to test their own control strategies.45-48 Our aim was not testing different control strategies for the WWT COST Benchmark but showing that global optimization techniques can efficiently be applied to the case of integrated design of a plant. Conclusions Optimization problems which arise during an integrated design and control approach are frequently multimodal. Thus, it is not surprising that global optimization (GO) has received much attention during the past decade from many researchers in computer-aided process engineering. However, the current state of the art is far from satisfactory, especially when we consider the global optimization of complex process models which are typically expensive to evaluate. In this contribution, we have considered a number of recent global optimization methods, including several surrogate-based approaches, and we have evaluated their performance based on the results for a benchmark problem in the management of wastewater treatment systems. Numerical experiments with the different GO solvers indicate that the proposed optimization problem is indeed multimodal and that, as expected, standard (local) solvers converge to local optima. In contrast, the results obtained show that surrogate-based methods, in particular RBF’s and quadratic interpolation, can indeed reduce computation times significantly while ensuring convergence to the best known solutions.

Acknowledgment We acknowledge the support of the EU Marie Curie programme (PRISM Project, No. MRTN-CT-2004-512233). The authors thank Dr. Ulf Jeppsson (IEA, Lund University of Technology, Sweden) for providing his MATLAB/ SIMULINK implementation of the COST Benchmark. J.A.E. gratefully acknowledges financial support (FPU fellowship) from the Spanish Ministry of Education and Science. Literature Cited (1) Alex, J.; Beteau, J. F.; Copp, J.; Hellinga, C.; Jeppsson, U.; MarsiliLibelli, S.; Pons, M. N.; Spanjers, H.; Vanhooren, H. Benchmark for Evaluating Control Strategies in Wastewater Treatment Plants, ECC’99, Karlsruhe, Germany, Aug. 31-Sep. 3, 1999. (2) Floudas, C. A.; Akrotirianakis, I. G.; Caratzoulas, S.; Meyer, C. A.; Kallrath, J. Global Optimization in the 21st Century: Advances and Challenges. Comput. Chem. Eng. 2005, 29 (6), 1185-1202. (3) Biegler, L. T.; Grossmann, I. E. Retrospective on Optimization. Comput. Chem. Eng. 2004, 28 (8), 1169-1192. (4) Floudas, C. A.; Pardalos, P. M., Frontiers in Global Optimization; Kluwer Academic: Boston, MA, 2004. (5) Bjo¨rkman, M.; Holmstro¨m, K. Global Optimization of Costly Nonconvex Functions Using Radial Basis Functions. Optim. Eng. 2000, 1 (4), 373-397. (6) Gutmann, H. M. A Radial Basis Function Method for Global Optimization. J. Global Optim. 2001, 19 (3), 201-227. (7) Jones, D. R.; Schonlau, M.; Welch, W. J. Efficient Global Optimization of Expensive Black-Box Functions. J. Global Optim. 1998, 13 (4), 455-492. (8) Huyer, W.; Neumaier, A. SNOBFIT-Stable Noisy Optimization by Branch and Fit, Universita¨t Wien. (9) Moles, C. G.; Gutierrez, G.; Alonso, A. A.; Banga, J. R. Integrated Process Design and Control via Global Optimization-A Wastewater Treatment Plant Case Study. Chem. Eng. Res. Des. 2003, 81 (A5), 507517. (10) Copp, J., The COST Simulation BenchmarksDescription and Simulator Manual; Office for Official Publications of the European Communities: Luxembourg, 2002. (11) Schweiger, C. A.; Floudas, C. A., Interaction of Design and Control: Optimization with Dynamic Models. In Optimal Control: Theory, Algorithms, and Applications; Hager, W. W., Pardalos, P. M., Eds.; Kluwer Academic: Dordrecht, The Netherlands, and Boston, MA, 1997.

Ind. Eng. Chem. Res., Vol. 46, No. 26, 2007 9157 (12) Bansal, V.; Perkins, J. D.; Pistikopoulos, E. N.; Ross, R.; van Schijndel, J. M. G. Simultaneous Design and Control Optimisation under Uncertainty. Comput. Chem. Eng. 2000, 24 (2-7), 261-266. (13) Gill, P. E.; Murray, W.; Saunders, M. A.; Wright, M. H. Constrained Nonlinear Programming. Optimization; Elsevier North-Holland: New York, 1989. (14) Banga, J. R.; Balsa-Canto, E.; Moles, C. G.; Alonso, A. A. Improving Food Processing Using Modern Optimization Methods. Trends Food Sci. Technol. 2003, 14 (4), 131-144. (15) Banga, J. R.; Balsa-Canto, E.; Moles, C. G.; Alonso, A. A. Dynamic Optimization of Bioreactors-A Review. Proc.-Indian Acad. Sci. 2003, 69A (3-4), 257-265. (16) Floudas, C. A.; Pardalos, P. M. Recent Developments in Deterministic Global Optimization and Their Relevance to Process Design. Fifth International Conference on Foundations of Computer-Aided Process Design, 2000. (17) Grossmann, I. E. Global Optimization in Engineering Design; Kluwer Academic: Dordrecht, The Netherlands, and Boston, MA, 1996. (18) Guus, C.; Boender, E.; Romeijn, H. E. Stochastic Methods. In Handbook of Global Optimization; Horst, R., Pardalos, P. M., Eds.; Kluwer Academic: Dordrecht, The Netherlands, 1995. (19) Jones, D. R. A Taxonomy of Global Optimization Methods Based on Response Surfaces. J. Global Optim. 2001, 21 (4), 345-383. (20) Simpson, T. W.; Poplinski, J. D.; Koch, P. N.; Allen, J. K. Metamodels for Computer-Based Engineering Design: Survey and Recommendations. Eng. Comput. 2001, 17 (2), 129-150. (21) Holmstro¨m, K.; Edvall, M. M. The Tomlab Optimization Environment. In Modeling Languages in Mathematical Optimization; Kallrath, J., Ed.; Kluwer Academic: Boston, MA, 2004. (22) Jones, D. R. DIRECT Global Optimization Algorithm. In Encyclopedia of Optimization; Floudas, C. A., Pardalos, P. M., Eds.; Kluwer Academic: Dordrecht, The Netherlands, 2001. (23) Gill, P. E.; Murray, W.; Saunders, M. A.; Wright, M. H. User’s Guide for npsol 5.0: A FORTRAN Package for Nonlinear Programming; SOL 86-1; Systems Optimization Laboratory, Stanford University: Stanford, CA, 1998. (24) Optimization Toolbox for Use with Matlab, User’s guide, Version 2; The MathWorks Inc. (25) Abramson, M. A. Pattern Search Algorithms for Mixed Variable General Constrained Optimization Problems; Rice University: Houston, TX, 2002. (26) Storn, R.; Price, K. Differential Evolution-A Simple and Efficient Heuristic for Global Optimization over Continuous Spaces. J. Global Optim. 1997, 11 (4), 341-359. (27) Runarsson, T. P.; Yao, X. Stochastic Ranking for Constrained Evolutionary Optimization. IEEE T. EVolut. Comput. 2000, 4 (3), 284294. (28) Csendes, T. Nonlinear Parameter Estimation by Global Optimization-Efficiency and Reliability. Acta Cybernet. 1988, 8, 361-370. (29) Boender, C. G. E.; Rinnooy Kan, A. H. G.; Timmer, G. T.; Stougie, L. A Stochastic Method for Global Optimization. Math. Program. 1982, 22 (2), 125-140. (30) Ja¨rvi, T. A Random Search Optimizer with an Application to a Maxmin Problem; Institute of Applied Mathematics, University of Turku: Turku, Finland, 1973. (31) Jeppsson, U.; Pons, M.-N. The COST Benchmark Simulation Model-Current State and Future Perspective. Control Eng. Pract. 2004, 12 (3), 299-304.

(32) Lindberg, C. F. Control and Estimation Strategies Applied to the Activated Sludge Process; Uppsala University: Uppsala, Sweden, 1997. (33) Lukasse, L. J. S. Control and Identification in Activated Sludge Processes; Wageningen University: Wageningen, The Netherlands, 1999. (34) Weijers, S. Modeling, Identification and Control of Activated Sludge Plants for Nitrogen Removal; Eindhoven University of Technology: Eindhoven, The Netherlands, 2000. (35) Carucci, A.; Rolle, E.; Smurra, P. Management Optimisation of a Large Wastewater Treatment Plant. Water Sci. Technol. 1999, 39 (4), 129136. (36) Ladiges, G.; Gunner, C.; Otterpohl, R. Optimisation of the Hamburg Wastewater Treatment Plants by Dynamic Simulation. Water Sci. Technol. 1999, 39 (4), 37-44. (37) Chachuat, B.; Roche, N.; Latifi, M. A. Dynamic Optimisation of Small Size Wastewater Treatment Plants Including Nitrification and Denitrification Processes. Comput. Chem. Eng. 2001, 25 (4-6), 585-593. (38) Balku, S.; Berber, R. Dynamics of an Activated Sludge Process with Nitrification and Denitrification: Start-up Simulation and Optimization Using Evolutionary Algorithm. Comput. Chem. Eng. 2006, 30 (3), 490499. (39) Sendı´n, O. H.; Moles, C. G.; Alonso, A. A.; Banga, J. R. MultiObjective Integrated Design and Control using Stochastic Global Optimization Methods. In The Integration of Process Design and Control; Seferlis, P., Georgiadis, M. C., Eds.; Elsevier: Amsterdam, The Netherlands, and Boston, MA, 2004. (40) Schu¨tze, M.; Butler, D.; Beck, M. B. Optimisation of Control Strategies for the Urban Wastewater System-An Integrated Approach. Water Sci. Technol. 1999, 39 (9), 209-216. (41) Jeppsson, U. Modelling Aspects of Wastewater Treatment Processes, Lund, Sweden, 1996. (42) Henze, M.; Grady, C. P. L. J.; Gujer, W.; Marais, G. V. R.; Matsuo, T. ActiVated Sludge Model No. 1; IAWQ: London, 1987. (43) Taka´cs, I.; Patry, G. G.; Nolasco, D. A Dynamic-Model of the Clarification Thickening Process. Water Res. 1991, 25 (10), 1263-1271. (44) http://www.ensic.inpl-nancy.fr/benchmarkWWTP/Bsm1/ Benchmark1.htm. (45) Cadet, C.; Beteau, J. F.; Carlos Hernandez, S. Multicriteria Control Strategy for Cost/Quality Compromise in Wastewater Treatment Plants. Control Eng. Pract. 2004, 12 (3), 335-347. (46) Gernaey, K. V.; Jorgensen, S. B. Benchmarking Combined Biological Phosphorus and Nitrogen Removal Wastewater Treatment Processes. Control Eng. Pract. 2004, 12 (3), 357-373. (47) Pons, M.-N.; Potier, O. Control Strategy Robustness with Respect to Hydraulics. Control Eng. Pract. 2004, 12 (3), 349-356. (48) Zarrad, W.; Harmand, J.; Devisscher, M.; Steyer, J. P. Comparison of Advanced Control Strategies for Improving the Monitoring of Activated Sludge Processes. Control Eng. Pract. 2004, 12 (3), 323-333.

ReceiVed for reView April 10, 2007 ReVised manuscript receiVed September 18, 2007 Accepted September 26, 2007 IE0705094