Nonlinear Multi-Criteria Optimization for the Integrated Design and Control of Bioprocesses

Copyright © IFAC Computer Applications in Biotechnology, Nancy, France, 2004

ELSEVIER

IFAC PUBLICATIONS www.elsevier.comllocatelifac

NONLINEAR MULTI-CRITERIA OPTIMIZATION FOR THE INTEGRATED DESIGN AND CONTROL OF BIOPROCESSES

Oscar H. Send in, Carmen G. Moles, Antonio A. A1onso and Julio R. Banga·

Process Engineering Group, Instituto Investigaciones Marinas (CSIC) Eduardo Cabello, 6, 36208 Vigo, SPAIN - e-mail: [email protected]

Abstract: The simultaneous design and control of continuous bioprocesses is considered here as a multi-objective optimization problem subject to non-linear differential-algebraic constraints. This type of problem is very challenging to solve due to its non-convexity (multi-modality). We present two novel solution approaches based on extensions of an stochastic global optimization method. The robustness and efficiency of these approaches are illustrated with a wastewater treatment plant case study. Copyright © 2004 IFA C. Keywords: Multiple-criterion optimisation, global optimisation, integrated plant control, biotechnology, waste treatment.

1. INTRODUCTION

Schweiger and Floudas (1997). Further, as also noted by these authors, in these problems one must find the best alternative for several (often conflicting) objectives, so a multi-objective optimization framework should be considered. Recently, Moles, et aI., (2003) considered the global optimization of a composite function including weighted economic and controllability objectives. In fact, this is a common approach for solving multi-objective optimization problems since it provides a compromise solution between all the objectives. However, choosing suitable values for the weights can be difficult.

Most bioprocesses present highly nonlinear dynamics which can lead to challenging controllability problems. The optimal design of bioprocesses operated in continuous mode (around a certain steady state) is usually carried out in twosteps. First, the optimal steady state with respect to some economic measure is found using an ad hoc procedure or standard mathematical optimization techniques. Second, controllers are placed and tuned in order to obtain a satisfactory response. However, this common practice of considering process control issues as a sequential second step after process design ignores the interaction of design and control. In other words, it does not take into account the economic consequences of the process dynamics arising from such a design. In the area of bioprocessing, wastewater treatment plants are good examples of how the traditional two step method can lead to designs with very poor controllability.

In this paper we consider the real multi-objective formulation for the optimization problems arising from integrated design and control, and we present two alternative methods for their solution. These methods are ultimately based on extensions of a recent stochastic global optimisation strategy for constrained non-linear programming problems. The usefulness and efficiency of these novel approaches are illustrated considering the integrated design of a wastewater treatment plant model.

In

order to surmount these difficulties, a simultaneous approach, considering operability together with the economic issues, has been suggested by a number of researches in process systems engineering (e.g. Morari and Perkins, 1994; Schweiger and Floudas, 1997; Pistikopoulos and van Schijndel, 1999). However, the solution of the associated optimization problems is not trivial. To begin with, these problems are frequently multimodal (non-convex), as described by e.g.

2. PROBLEM STATEMENT For most real applications, the optimisation problem should be a multi-objective one, i.e., it should consider the simultaneous optimisation of two or more objectives, which are usually conflicting. In general, there is no feasible solution which is

571

simultaneously optimal for all the objectives, but instead a set of solutions that constitute the relatively best alternatives can be found. This is known as the Pareto-optirnal (or non-dominated) solutions set. All the points in this set are optimal in the sense that it is not possible to improve one objective without degrading one or more of the others.

andlor discontinuities. A large population size is usually required, which is translated into a large Pareto-optimal set in such a way that it is difficult to distinguish between two adjacent solutions. Additionally, they are very expensive in terms of computational effort. In order to surmount the limitations of both approaches, and after an initial screening of a number of global optimisation approaches (Moles et aI, 2003), here we propose two extensions of the SRES (Stochastic Ranking Evolutionary Strategy) algorithm of Runarsson and Yao (2000). These extensions consider the solution of the set of NLPs generated from two different frameworks: the traditional E-constraint technique, and the recent Normal Boundary Intersection (NBD method of Das and Dennis (1998). It should be noted that the SRES algorithm was developed as a novel global technique for mono-objective optimisation problems, so the methods we are presenting here can be seen as extensions of SRES for multi-objective problems. In the following sections they are described briefly.

When the integration of process design and control is considered, the problem is translated to finding the equipment static design variables, the operating conditions and the controllers' parameters which optimise the plant economics and, simultaneously, a measure of the plant controllability, subject to a set of constraints which ensure appropriate dynamic behaviour and process specifications (here it is assumed that the process flowsheet is given). The general mathematical formulation of the problem is:

min x

F( t,z,p,x ) = [Fl (t, Z, p, X)] F2 (t,z,p,x)

(1)

subject to:

f(t,z,p,x) = 0

(2)

z(to) = Zo

(3)

h(z,p,x) =0 g(z,p,x) $ 0 xL $x $x u

(4)

3.1 E-constraint SRES method Our strategy is based on the combination of the well known E-constraint method and the recent SRES method for the solution of the resulting NLPs. In the E-constraint technique, the multi-objective problem is reduced to a single NLP by minimizing one of the objectives and converting the rest to inequality constraints by means of a characteristic parameter, E. Mathematically, the problem statement is as follows:

(5) (6)

where x is the vector of decision variables, z is the vector of dynamic state variables, F is the vector of objective functions (FJ is a combination of capital and operation costs, and F2 is the controllability measure), f is the set of differential and algebraic equality (DAEs) constraints describing the system dynamics, and h and g are possible equality and inequality path andlor point constraints.

mm x

F2 (t,z,p,x)

$

E

(8)

and the same set of constraints given by Eq. (2) to (6). The solutions to the E-constraint problem is expected to be a Pareto-optirnal solution, so a vector of different E values produces the set of Paretooptimal solutions. It should be noted that prior information is required in order to choose appropriate values of E for generating a good approximation of the entire Pareto set.

During the last decades, many algorithms have been proposed for finding the set of Pareto-optimal solutions. A review can be found in Deb (2001) and Ehrgott and Gandibleux (2002), and the references cited therein. Basically, there exist two main groups:

•

(7)

subject to:

3. MULTIOBJECTIVE OPTIMIZATION METHODS

•

Fl (t,z,p,x)

'Classical' approaches transform the original non-linear multi-objective optimisation problem into a single-objective non-linear programming (NLP) problem, whose solution via a local method is expected to be Pareto-optirnal. Multiple solutions can be obtained by changing the characteristic parameter of the method and solving the resulting NLP. Although easy to implement, their main drawback is that these methods usually fail in finding Pareto points in non-convex problems. Genetic and evolutionary algorithms can find multiple Pareto-optimal solutions in one single optimization run, instead of solving a set of NLPs. They can easily deal with non-convexities

3.2 Stochastic NBI (NBI-SRES) The recent NBI method developed by Das and Dennis (1998) produces an even spread of points on the Pareto front by transforming the original nonlinear multi-objective optimization problem into a set of NLPs which are solved sequentially. The characteristic parameter of the method is a weighting vector ~, which defines a point in the so-called Convex Hull of Individual Minima (CHIM). Starting from that point, the normal to the CHIM will intersect the boundary of the objective space in a

572

point that will be a Pareto-optimal solution if the Pareto front is convex and the individual minima of the objectives are the global ones. Thus, the NLP subproblem is defined as maximizing the distance along the normal direction, subject to a set of additional equality constraints which ensure that solution is actually on its corresponding normal to the CHIM (called here NBI-constraints).

The aim of the control system is to keep the substrate concentration at the output (S1) under a certain admissible value. The main disturbances come from large variations in both the flowrate and substrate concentration (qj and SI) of the input stream. The flowrate of the sludge recycle to the first aeration tank is taken here as the manipulated variable (Gutierrez and Vega, 2000);

This method has the advantage that an equally distributed set of fJ produces an equally distributed set of points on the Pareto front, which is an useful feature for helping the decision-making process. But, as a drawback. it can only handle unimodal problems, since the original MATLAB implementation works by solving the set of NLPs by means of SQP (Sequential Quadratic Programming), which is a local gradient-based method. Thus, it usually fails with non-convex problems.

The dynamic model consists ofa set of33 DAEs (14 of them are ODEs) and 44 variables. The value of three flow rates (qr1, qrJ and qp) are fixed at their steady-state values corresponding to a given nominal operational conditions. Therefore, this leaves 8 design variables for the integrated design problem, namely the volume of the aeration tanks (vJ and V1), the areas of the settlers (adJ and ad1), the aeration factors (fTeJ and fk1), and the gain and the integral time of a PI controller.

In order to improve the robustness of the technique, we have replaced the SQP solver by SRES. The additional NBI-constraints are added to the objective function as a quadratic penalty function. This technique works well for our purposes, since an approximate solution close to the normal suffices.

The integrated design problem is formulated to minimize a weighted sum of economic terms (FJ) and simultaneously a measure of control performance taken here as the Integral Square Error (lSE): FI

4. CASE STUDY: A WASTEWATER TREATMENT PLANT

(9) (10)

This case study, described by Gutierrez and Vega (2000), represents an alternative configuration of a real wastewater treatment plant located in Manresa (Spain). Fig. I shows the process flowsheet, which consists of two aeration tanks, acting as bioreactors, and two settlers. The biodegradable pollutants (substrate) are transformed by a flocculating microbial population (biomass) inside each bioreactor, in which the necessary level of dissolved oxygen is provided by the aeration turbines. The effluents from the aeration tanks go to their associated settlers where the activated sludge is separated from the clean water and recycled to the corresponding reactor. The activated sludge in excess is eliminated via a purge stream (qp).

The weighting vector Wj in the economic objective is chosen in such a way that it implies a similar contribution for each term in the objective function (Moles, et al., 2003, Gutierrez and Vega, 2000). The ISE is evaluated considering a step disturbance to the input substrate concentration, Sj, whose behaviour is taken from the real plant. The objectives above are subject to several sets of equality and inequality constraints. Besides the 33 model DAEs (system dynamics), which acts as equality constraints, there exist a set of 32 inequality constraints which impose limits on residence times and biomass loads in the bioreactors, the hydraulic

ql SI Xi

T

ST

qp

=(wI .v[ )+(W2 'vi )+(w3·adl2)+(w4.adi)+ (ws .fkl2 )+ (W6 .fki )

XT

Fig. I. Process flowsheet of the wastewater treatment plant.

573

capacity in the settlers, the sludge ages in the decanters, and the recycle and purge flow rates. Also, an additional set of 120 double inequality constraints which impose limits on the state variables.

~~ .....:... , .. ..'

· ISO .:.. .. -1 25 ., .....

~ 5. RESULTS AND DISCUSSION

~_ -I

The Pareto-optimal solutions set was generated considering the scenario described in Moles et al. (2003). The corresponding settings for the methods proposed here were as follows: . •

•

.'

I) Values of the ISE between 0.28 and 1.0: It can be observed that the cost can be greatly reduced while maintaining a good controllability. Since cell recycle permits that the bioreactor can be operated at dilution rates higher than the specific growth rate, it is obvious that a reduction in the plant economics is achieved by an increase in the flow rate of the recycling sludge (qr,), which is mainly translated to a reduction in the volume of the first tank. 1

0

le

14

12

I:

0

0.8

• Il) Values of the ISE between 1.0 and 9.5: In this region, there is a point from which no further improvement can be achieved in the cost function. Only the E-constraint SRES technique was able to find feasible solutions in this region. It should be noted that this strategy minimises the economic function while imposing the constraint to the ISE measure, i.e., we are forcing the controller to have a low lSE, but at the expense of increasing the oscillatory behaviour, and even the instability of the system. It is important to note that, for the majority of designs in this region, their corresponding E-constraint (Eq. 8) is not active, which only happens when concave parts and/or discontinuities are present in the Pareto front.

,.-con. traint SRL<;

NBI·SRES

OX

IC

0

0.6

le le

I

0

oX

..

ox

OA 0 .2 woo

1100

x 1200

X

1300

0

o

.

CD

~ /

•

ocxxx!XXIeXX

'~~~~~~~~~~--~--~------~ ~ I~ IM I~ I@ I ~~~~ ~

Fe...

Fig. 2. Pareto-optimal solutions wastewater treatment plant.

set

for

~_l

Consequently, the biomass concentration rise up to a level in which some constraints regarding biomass concentration (mainly that of the top layer of first decanter, xd,) become almost active (see Fig. 3). This behaviour is more obvious in the solutions obtained with the E-constraint SRES technique. Additionally, certain constraints are always active (or very close to being active). They are referred to the lower limits in the sludge age in the second decanter at steady-state, and in both settlers after disturbance. An important difference between both Pareto-optimal sets is found in the dissolved oxygen concentrations inside the bioreactors. They are always in their lower bounds in the solutions obtained with the E-constraint SRES technique, while those of NBI-SRES present higher levels (consequently, the aeration factors are different too).

the

)(

2000 . 1500 1750 - . Fc.. 1000 I<'SO .

Fig. 3. Variation, along the Pareto set, of biomass concentration in the top layer of first settler after disturbance.

NBI-SRES technique, the shadow minimum (or utopia point, i.e., the vector containing the individual global minima of the objectives) was found by minimizing each objective function with SRES, and refining the solutions with SQP. Then, 24 NLPs were solved with the SRES algorithm for a population size of 200 which was evolved for 200 generations.

16

::::j::::::~:. ...:. ..

\SE

The Pareto-optimal solutions found with these strategies are shown in Fig. 2. Regarding ISE values, we can distinguish several regions in the Pareto front.

•

....

t· ..·

= ii

In the E-constraint SRES approach, the population size and the number of generations were fixed at 200 and 300, respectively. The values of E were chosen taking into account several preliminary runs in order to find a good distribution of points. The number of NLPs solved was 26. In

·r:l. · ·. ·.·. . ~ "."'1'"' +

the

574

IIl) Values of the ISE greater than 9.5: This region is a virtually vertical line where all designs have the same cost. It is clear that it is possible to obtain different dynamic responses (and thus different ISE values) for the same design depending on the parameters of the controller, which have no effect on the economic objective.

Table 1. Representative designs from the Pareto-optimal set.

VI V2

adI ad2 flcI flc2

Kp T; F cort ISE Method

Design A 4669.7182 9998.7969 2853.5886 2211.8980 0.7510 0.8219 -100.0000 1.0832 2580.8750 0.2822 NBI-SRES

Design B 6428.7663 4280.1057 2154.4923 3894.7295 0.0560 0.0123 -99.9999 1.1517 1391.1138 0.3202 E-constraint

Design C 5415.4870 4212.0830 2341.6368 3994.3647 0.3909 0.0665 -99.9986 0.7453 1157.6515 0.4001 NBI-SRES

In this region, the constraints related to the dissolved oxygen concentrations are active for both sets of Pareto-optimal solutions. In addition to these constraints and those mentioned previously, the biomass concentration in the inlet stream to the first bioreactor (xirJ) is active too. Thus, the system is so constrained that no further reduction in the cost is possible. In fact, in region II it is possible to find designs as economic as these ones, but the constraint to the ISE and the restrictions of the system lead to unstable dynamics.

DesignD 5512.8773 3986.5565 2291.7307 3998.6114 0.06742 0.0112 -99.9999 0.7400 1138.1542 0.4000 E-constraint

Design E 4813.1861 4007.8199 2471.5581 3990.9076 0.0794 O.oI13 -99.9957 0.6411 1005.0243 1.4999 E-constraint

Design F 4737.1086 4067.2123 2493 .9900 3969.9916 0.0809 0.0117 -99.9982 13.3338 999.5367 9.4172 Both

10.8 20.7 20.6 10.5

~

. 21!.4

!. ~

203 ~O. ~

20.1 20.8r----..----.------.'""i===::===;--, - Design A I 1- DcsignB 20.7

2°0

600 ri_Ibr)

8(K)

HliM)

We have analysed the dynamic behaviour of several representative designs of these regions. The values of the decision variables and the objective functions are given in Table 1. Designs A and B correspond to the designs with minimum ISE found with both methods, whose dynamic responses for the outlet substrate concentration are plotted in Fig. 4. Designs C and D are compromise solutions between the two objectives. Note that both designs have similar values for the objectives, but the decision variables are rather different, specially regarding the aeration factors. Their corresponding dynamics are very similar (not shown). Design E belongs to region 11, in which the concentration of biomass in the top layer of the first decanter is at its upper limit after the disturbance. Although the ISE is relatively low, the system presents a highly oscillatory behaviour, as it can be seen in Fig. 5. Finally, design F corresponds to the most economic design with an ISE below 10, i.e. a probably acceptable dynamic response (Fig. 6).

20.5

203

o

4(~)

Fig. 6. Dynamic response of the outlet substrate concentration for design F.

20.6

'it=

200

50

150

100

200

250

Time (hr)

Fig. 4. Dynamic responses of the outlet substrate concentration for designs A and B.

1(l( ~ I

1000

3000

-lrol

It should be noted that both solution approaches found difficulties when computing solutions in region 11. In Fig. 7, we show the set of solutions found with NBI-SRES and the lines that joins them with their corresponding starting points at the ClllM. Only a few set of solutions are feasible with respect to the NBI-constraints. The solutions of the NLPs whose normal should intersect the Pareto front in region 11 are moved up to the region rn, where they have approximately the same value of the cost

5000

Tim
Fig. 5. Dynamic response of the outlet substrate concentration for design E.

575

function but a more stable dynamic (and higher values of the ISE). Regarding computational effort, the majority of NLPs were solved in similar CPU times, which varied from 1 to lA hours (pentium IV1.7 GHz).

Their efficiency and robustness was illustrated by solving a multi-criteria optimisation problem related with a wastewater treatment plant model. It should be noted that this problem could not be successfully solved with standard gradient-methods. In contrast, the techniques presented in this paper allowed the computation of the Pareto-optimal solution set with reasonable effort.

On the other hand, the E-constraint SRES technique was able to find solutions in region n, but the cost to pay was a significant increase in the computational effort (Fig. 8). Furthermore, for these NLPs, the majority of individuals in the final population were not feasible. This is a consequence of the highly nonconvex and constrained nature of the search space in this region.

Although both methods found several difficulties when dealing with designs which were very close to dynamic instability, they exhibited complimentary features in that region, thus suggesting certain research directions in order to further increase their robustness.

6. CONCLUSIONS AND FUTURE WORK 7. ACKNOWLEDGEMENTS Here we have presented two novel strategies, based on multi-objective optimisation, for the integrated design and control of bioprocesses, modelled as nonlinear dynamic systems. These techniques are ultimately based on an stochastic global optimisation method, so they can deal adequately with the typical non-convexities of these problems.

Author Oscar H. Sendin acknowledges a pre-
16 14 I~

criteria optimization: State of the Art Annotated Bibliographic Surveys. Kluwer Academic

6

1200

1400

1600

1800

2000

2200

~400

Publishers, The Netherlands. Gutierrez, G. and P. Vega (2000). Integrated design of activated sludge process taking into account the closed loop controllability. Proceedings of ESCAPE-lO, pp. 63-69. Firenze. Moles, C.G., G. Gutierrez, A.A. Alonso and 1.R Banga (2003). Integrated design and control via global optimization: A wastewater treatment plant case study. Chemical Engineering Research & Design, 81, pp. 507-517. Morari M., and lD. Perkins (1994). Design for operations. In: Foundations of Computer-Aided

2600

Fe.... Fig. 7. Performance of the NBI-SRES approach in terms of feasibility with respect to the normal directions to the CHIM.

.........

·······r ;..... ,

x 10'

....:;....••.. .

:::".:r: :: .:·.

Process Design (FOCAPD). Pistikopoulos, E.N., R Ross and 1.M.G. Schijndel (1999). Towards the integration of process design, process control and process operability. Current status and future trends. In: Foundations

:.:.:f... :'.'..... .'.[..::::....':. .......

~ 2.5

....... "',

=

t

, . i"·

2 ·"

~L5

: ••••. . . ••

. ..... .....

.... ........ ]

!

•. •...

·····t'·· .. ... ~...........

..

".

~

; ......

. -' ·i : .•... ..... i .......;, . . ..r:::.· .· ··:.~........ . ... . ,J

>. .

of Computer-Aided Process Design (FOCAPD) .

5,; .- :::· r;~: :i"t:;'J .: to

....

::. _;: •. "

".:.:

Snowmass, Colorado. Runarsson, T.P. and X. Yao (2000). Stochastic ranking for constrained optimisation. IEEE

Transactions on Evolutionary Computation, 4:3,

.... .::.:.. ......

pp. 274-283. Schweiger, C.A. CA Floudas (1997). Interaction of design and control optimization with dynamic models. In: Optimal control: Theory. algorithms and applications. (W.W. Hager and P.M. Pardalos (eds.)). Kluwer Academic Publishers, The Netherlands.

1000

Fig. 8. Performance, along the Pareto set, of the Econstraint SRES approach in terms of computation time.

576