Optimal Non-Linear Design and Control of Bioprocesses via Linear Programming

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

ELSEVIER

IFAC PUBLICATIONS www.elsevier.comlloc:atelifac

OPTIMAL NON-LINEAR DESIGN AND CONTROL OF BIOPROCESSES VIA LINEAR PROGRAMMING Julio Vera\ Carmen G. MoIes 2, Julio Banga2 and Nestor V. Torres 1,*

Grupo de Tecnologia Bioquimica, Departamento de Bioquimica y Biologia Molecular. Facultad de Biologia .2Instituto Universitario de Bioorgimica Antonio Gonzalez. Avda. Astfeo. Feo. Sanchez. sin. 38206 La Laguna. Tenerife. Islas Canarias. Espaiia [email protected] J Process Engineering Group. Instituto de Investigaciones Marinas (CSIC). Eduardo Cabe//o. 6. 36308- Vigo. Espaiia. J

Abstract. An optimization approach for bioprocesses based on linear programming is developed and applied to a wastewater treatment plant. After defining an objective function reflecting both investment costs and "paracosts" (stability, flexibility and controllability), we define a set of constraints determined by the system components and technical and economical factors . Results obtained reveals that significant improvements in controllability and cost reduction are achieved. We conclude that the incorporation of this method to a dynamic process simulator would lead to a fast and interactive tool for obtaining near-optimal integrated designs for bioprocess plants. Copyright © IFA C 2004. Key words: Integrated design, bioprocesses, control and optimization, linear programming, power law formalism, waste-water treatment.

for mode ling. analysing and optimizing biochemical systems, the so called Indirect Optimization Method, IOM (Torres et al. 1997; Alvarez-Vasquez et a1. 2(00). The IOM approach is based on the fact that, although S-system models are noniinear, the steadystate equations are linear when the variables are expressed in logarithmic coordinates. The key argument now is that the logarithmic transformation does not change the locations of maxima and minima. In other words, the two functionsj(x) and In j(x), for positive x-values, have their maxima and minima for the same values of x. This paper is focused on demonstrating how IOM can be applied to integrated design problems ofbioprocess plants. In this context the selection of a wastewater treatment plant as a case study is not arbitrary: these plants do present many controllability problems in real practice due to bad designs. Thus, the optimization approach must consider the optimal steady state with respect to

I.INTRODUCTION The proper design and operation of any microbialbased biotechnological process requires the quantitative description of the variables relevant for the system"s kinetics. However, the nonlinear nature of these processes has impaired such optimizations (Vagners 1983). One approach that fully acknowledges the nonlinearities, yet ultimately leads to linear optimization tasks is based on a mode ling framework called Biochemical Systems Theory (BST; Voit 2(00). The hallmark of this theory is the approximation of rate laws and other processes with products of power-law functions. The BST framework offers some choices for the formulation of biochemical systems, among which the most relevant for our present purposes is the S-system (Voit 2000). Within this formalism, it has been presented and applied to different systems an integrated approach

469

some economic measure. However, also important from the economic point of view is the process dynamics arising from the given design. In this vein, a suitable objective function should reflect some monetary costs or investment costs and also some "paracosts" such as stability, flexibility or controllability (G6rak et al. 1987). The constraints are dictated by the system components and by technical and economical factors. Moreover, we will compare the optimum profiles and the optimized responses with the ones obtained for the same problem by using selected global optimization methods, like the GLOBAL method (Csendes, 1988; Moles et aI2001).

qr,

2. CASE STUDY: WASTE-WATER TREATMENT PLANT MODEL

q.

i

x, .• •.c,

.J0..' .• r-~ ~ ._ r ~ ' tank i ~.-. " -" "",

In order to test the applicability and reliability of the proposed optimization approach to process design and control we have used a waste-water treatment plant model previously presented by Gutierrez and Vega (2000). In Figure I the plant scheme is presented. The plant is composed by two aeration tanks working in series acting as bioreactors and two settlers. The settler model is inspired on that by Takacs et al (1991). A flocculating microbial population (biomass) that transforms the biodegradable pollutants (substrate) is kept inside each bioreactor. Oxygen is provided by aeration turbines. Effluents coming from the aeration tanks are separated in their associated settlers into a clean water stream and an activated sludge, which is recycled to the corresponding aeration tank. Due to the excess of activated sludge production, the amount that can not be transferred to the recycle tanks is eliminated via a purge stream. The objective of the control system is to keep the substrate concentration at the output (qsal) under a given admissible value.

ar.

x.

q.

q.~;T.:....:-~=···~.·:.· . __.;:.q,_. _._-:q,:.-__ tN PUT

Figure 1. Plant scheme of the activated sludge process. Table 1. Model's S-System equations dX 1 =0.213 XIO.OSSX30.049XllO.908XI20.014X130.989 dt

X I4-O·944 X22.o·IS8 X24-O.06S _ 0.213 XIl.Ol6 X13 0.293 X3 -0.016 XI4·0.983

dX 2 = 0.203 X20.033 ~0.031 X70.171 X120.794X220.1 71 dt

X23 0.794 XIS·O.966 _ 0.203 X/ 003 X220·8S7X230.138 -0.003 XISO.996 dX, = 1.572 X IO.06 X30.006~0.001 X130.039 X14.o·969 ~

dt

xt

90S X13 0.039 X22 -0.003 X24-0.001 _ 1.572 X IO.866 X I40. 133 dX. = 0.345 X 0.OO9 X3 0.912 ~0.077 X22 0.912 X23 0.081 2

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.

dt

XIS -0.994 _ 0.345 X20.4123 ~0.973 X22O.SOS X23 0.081 XIS -O.S87 dX s =0.757 X 18 - 0.757 X IO.OOI X30.OO1 X sO.998 X130.082 dt

The mathematical model consists of 14 ordinary differential equations and 19 algebraic equations, with 14 dependent and II independent variables. In Table I the classification, definitions, basal values and equivalencies among the original and S-system nomenclature are shown.

XI80.722 X I4-0.2764

dX6 =0.737XsO.I~190.809X220. 19

X IS.o·19

dt

_ 0.737 X 20.0000S ~o.oooos ~ XI90.724 X22o.236 X23 0.038 XIS -O.27S dX, = 1.648 X9 X22 X 16·1 - 1.648 X70. 9644 XI6·0.097 dt

X 220. =1.756XIOX221.04XI7·IX24.o·04 097

The S-System model version was derived through aggregation of the input (output) fluxes in a.r (T) flux function and then expanding this in Power-Law formalism for each equation (for details see Voit,

dX. dt

2000).

dX9

dt

470

-1.756Xs0·~220. 108 X17-O·I04X24-0.004 = 0.735XI o.884xl ll OX 130.264 X I6-0.884

- 0.735 X 9o.669 X I3 0.033 XI6:O·18 dX 'O = 1.004X20.906 X 80.09 X 22 0.18 X230.126 X I1.o·906 dl

_ 1.004 XIO0.916 X 0.108 X23 0.011 X .o·12 22 I1 dX 027X 0.621 X O.06IX .o·06X .0.146 _ " = ' 9 \3 16 22

XI

Biomass conc. bioreactor I

X2

Biomas conc. bioreactor 2

SI

Organic matter conc. bioreactor 1

XI X2 X3

S2

Organic matter conc. bioreactor 2

X.

dJ

Cl

Oxygen conc. bioreactor 1

Xs

dl

Oxygen conc. bioreactor 2 C2 xd l Suspended org. matter conc. settler 1

X1

xdz Suspended org. matter conc. settler 2

X8

xb l Settling organic matter conc. settler 1

X9

- 0.27 X l1 X\31.02 X 16-1X 22 -2.43 dX 0 153X 100.91X 23 o·02X 24O·OO4X 11.0.024 _'2 = ' -0. 153X I2X 23 0.8~240. 19 X I1-1 0.992 X210.118 X220.142X230.02 dXll 0.045 X 20. 1l5

-=

dl

x.

xb2 Settling organic matter conc. settler 2

XIS .o.16S X .o·118 _ 0.045 X 0.OO2 X30.2S6 X. 0.021 20 2 X 21 0.118 X 220.2S6 X 23 0.023 X IS .o.219 X 20.0.118

xrl

2.1. Simultaneous Optimization ofthe design and control of the waste-water treatment plant model

XIO Sedimen.organic matter conc. settler 1 X l1

xr2

Sedimen.organic matter conc. settler 2

X 12

qrl

Activ. sludge net flux to bioreactor 1

X\3 X I4

VI Bioreactor 1 volume V2 Bioreactor 2 volume

XIS

ad 1 Settler 1 surface

One significant advantage of representing a complex system with an S-system rather than some other nonlinear model is the fact that its steady state can be expressed explicitly in the form of a system of linear equations. One important consequence is that the steady state itself, as well as derived features like sensitivities, are easily computed. The numerical values for the system variables at steady state coincides with those previously defined (Gutierrez and Vega 2000). The local stability is confirmed by exammmg the eigenvalues for the local representation of the system showing that the real parts of all the eigenvalues are negative (results not shown). Moreover, there are imaginary parts, indicating that dumping or oscillatory behavior is to be expected, which coincides with the observed experimental behavior. The quality of the representation was subsequently analyzed in order to check whether the model represented a realistic picture of the system and was thus worth optimizing (results not shown). Once we are provided with an Ssystem model version, and the steady state analysis and quality assessment have been carried out, we are in conditions to perform the logarithmic transformation and subsequently, to define the linear programming problem. The following features guarantied comparison of the IOM results with those obtained by Moles et al. X 22 , X23 and X 24 were kept at the same value (all equal I) of the original model as well as the control configuration of the plant (fixed qr1 value). In the objective function the only paracost that we have considered was controllability (ISE). The optimization program must define a feasible space:: where large disturbances (e.g. the washout of organism) are avoided.

X I6

adz Settler 2 surface

X I1

fkl

Bioreactor 1 aeration factor

X I8

fk2 ,i

Bioreactor 2 aeration factor

X I9

Integral time of the controller

X 20

Is,

Gain of the controller

ql

Flux settler 1 ~ bioreactor 2

X 21 X 22

qr2

Flux settler 2 ~ bioreactor 2

X 23

qr3 Flux settler 2 ~ bioreactor 1

X 24

The implementation of this step is as follows : Objective function. The optimization problem consists of evaluating the physical and control parameters of the plant simultaneously, while the investment and the operation costs are minimized. In addition, the dynamics of the system should be within desired ranges. These desired dynamic characteristics of the system are good for rejection of the disturbances at the output variable S2 and to ensure the stability of the resulting system. Accordingly an objective function, C, was originally defined as a weighted sum of controllability (ISE) and economic (41) cost terms as: Minimize C = 0)1 . ISE + 41 The economic term, 41 has the following expression: ; = (m, '\1~) + (m. '\1;) + (m.ad,' ) + (m, 'ad:) + (m. ft,') + (m, ·ft;)

where VI, V2, ad l and ad2 are related to the bioreactor volumes and the surface of the settler tanks respectively and fkl and fk2 to the aeration devices. In all cases 0) represents the vector of weighting factors, 0) = [103,2'IO-s, 2'IO- s , lO- s , lO- s , 12, 12]. If we translate the objective function to Power-Law formalism, the function takes the form:

Table 2.Variables. and S-system equivalence for the model optimization of the waste-water treatment plant.

",s-s

'f'econ

where : Definition

~

S-S

471

=r~a

econ

.II~9 x~econ'i] 1=14

1

organic matter present in the fluxes coming from the mixing flux nodes.

.i =14.15 ... .19

Table 4. Mixing coefficient constrains

The objective of the control system is to keep the substrate concentration of the output (S2) close to a predefined setpoint value. ISE is the integral square error of the output (S2) evaluated considering a step disturbance to the input substrate concentrations, Si j • The integral ISE function can be substituted by an equivalent one as the following:

Variable Linear constraint Xirl -0.781<0.961Yll+O.015YI2+O.748YI3 - 0.167Yn -0.069y24< 0.014 Sirl -0.728<0.037y3+O.OO2Y4- 0.257 Yl3 - 0.003 Yn - 0.001 Y24 < 0.271 Xir2 -0.083 <0.177 Y7 + 0.8222 Y12 - 0.6827 Yn + 0.6827 Y23 < 0.9161 Sir2 -0.061<0.9176 Y3 + 0.0824 Y4 + 0.0571 Y22 - 0.0571 Y23 < 1.1607 Xr -0.597<0.984Yll+O.015yI2+O.072 Y13 - 0.171 Y22 - 0.07 Y24 < 0.044

In this equation x.,.,.,t represents the control variable and Ppert the perturbation parameter, that represents the step disturbance. So expressing the maximization of dXa.n/dPpert is equivalent to maximizing the sensitivity of the flux through the control variable with regard to the perturbation parameter. And due to the fact that the flux changes are controlled by changes in the perturbation parameter, maximizing the system's controllability is just maximizing the sensitivity of the flux with respect to the perturbation parameter. Using the power-law fonnalism it is possible to transfonn the above expression into the following, which will be used in our optimization set up:

Dependent variable constraints. The dependent variables (Table 1) were allowed to vary between different limits ranging from 0.75 to 10 times the basal values. As an exception, an additional constraint was imposed in order to assure that at any optimized solution found the value of J4 was the predefmed "set point" (S2)0. Table 5. Dependent variables constrains -1.3862 -0.2876 -0.2876 -1.3862 -1.3862 -0.6931 -0.2876 -1.3862 -0.2876 -0.6931 -0.2876 -1.3862

Max [0.116 X 2 + 0.993 J4 - X I4 - 0.965 XIS - 0.1 X I6 - 0.08 X I7 - 0.005 XIS - 0.005 X 19 - 0.719 X20 + 0.719 X 21 + 0.142 X 22 + 0.023 X 23 ]

Constraints. In addition to the steady state constraints there are four sets of constraints to be imposed in the fonnulation of the optimization problem. Flux constraints. These constraints were already imposed in the original presentation of the wastewater treatment system. They impose upper and lower limits to some dependent fluxes, which are related to the tank's dimensions. The constraints can be straightforwardly translated into power law expression.

::;; ::;; ::;; ::;; ::;; ::;; ::;; ::;; ::;; ::;; ::;; ::;;

YI Y2 Y3 Ys Y6 Y7 Ys Y9 YIO Yll YI2 Y13

::;; ::;; ::;; ::;;

::;; ::;; ::;; ::;; ::;; ::;; ::;; ::;;

0.2231 2.3025 1.6094 0.2231 0.2231 1.6094 2.3025 0.9162 2.3025 0.2231 1.6094 0.9162

Independent variables. We imposed limits to the independent variables ranging from 0.25 and 10 times the base-line steady-state values except for X22, X23 and X 24 that were kept at the base-line steadystate values.

Table 3. Flux constrains. yi stands In(Xi) Table 6. Indcmendent variables constrains Fluxes Q2 Q3 Q12 Q22 Qsal Qr

Linear constraints -0.4312 < 1.0250 Y13 - 2.4335 Yn < 0.7449 -0.1202 < 0.8061 Y23 + 0.1931 Y24 < 1.0559 -1.5690 < 0.2985 Y13 < 0.2761 -1.4847 < 0.8605 Yn + 0.1394 Y23 < 0.3604 -1.1013 <1.0403 Yn - 0.0403 Y24 < 0.3759 -1.072 < 0.936 YI3 - 2.223 Yn+ 0.086 Y24 < 0.529

Mixing coefficient constraints. As in the previous case these constraints were defined in the original model. They serve to limit the biomass and other

472

-1.3862 -1.3862 -1.3862 -1.3862

::;; ::;; ::;; ::;;

Yl4 YIs YI6 Y17

-2.3025 -2.3025 -2.3025 -2.3025

::;; ::;; ::;; ::;;

YI8 YI9 Y20 Y21

::;;

::;; ::;; ::;; ::;;

::;; ::;; ::;;

0 0 0 0 0 0 2.3025 2.3025

Steady state constraints. This set of constraints, derived from the S-system model formulation, assures that the optimized system is in a steady state:

3. RESULTS AND DISCUSSION As the . result of the optimization procedure above presented we obtained the optimum solution vector profile shown in Table 8. The optimized steady state is stable and does not violate any of the imposed constraints. This solution has a robust behaviour (results not shown) and better perturbation dynamics than the original set up. It can be shown that the time course of the substrate concentration X. (S2) after a disturbance in Si (approx. +3%) shows less and smaller oscillations than the basal steady state.

dX .. Ni --;t = 0 => a,n jXr - Pin JX) = 0

Reconfiguring equations and taking logarithms these constraints take the following form: 0= -0.96Yr+Q.065Y3+O.908Yll+O.014Yl2+O.695YI3 +O.039y14- 0.158yu -0.0652Y24 0=-0.9691 Y2+ 0.0348 Y4 +0.1718Y7 +0.7943YI2 + 0.0303 YIS - 0.6856 Yu + 0.6554 Y23 0=-0.8062 YI - 0.8986 Y3 +0.0019 Y4 -O.OOOIY130.8359 YI4 - 0.0034 Yu - 0.0014 Y24 O=-O.402Y2 +O.912Y3 - 0.896 Y4 - 0.4071 YIS +0.4071 Yu 0=-0.00 12YI-0.OO 1Y3-0.998Ys-0.082Y13+O.276y 14 +O.27YIs 0=-5·1 0·sYr5·10·sY4+O.191Ys-Y6+O.084YIs+O.084YI9 -0.045 Yu-0.038 Y23 0=-0.9645 Y7 + Y9 - 0.9022 YI6 + 0.9022 Yu 0=-0.9975 Ys + YIO - 0.8958 YI7 + 0.9319 Yu -0.0361 Y24 0=0.8847 YI + 0.1107Y7-0.6698Y9+O.2103 Y13 -0.7045 YI6 0=0.906Y2+O.0931 Ys-0.976yw-0.78y 17+O.67yu +O.108Y23 0=0.6214 Y9 - Yll - 0.9633 YI3 + 0.9398 YI6 + 2.2869 Yu 0=0.9733 YlO - YI2 + 0.9757 YI7 - 0.7873 Y23 - 0.1884 Y24 0=0.1 132Y2 -0.2566 Y3 +0.9708 y4+O.1144YIs -0.1144Y22

A second observation refers to the level of improvement both in cost and controllability obtained by using this method. In fact, examination of Table 8 shows that a significant improvement in both controllability and cost can be attained through the use of the IOM approach here applied. The ISE value changes downn to 5.34% of its initial value. Moreover, the cost function, changes down to 57.9% of the original value. Finally, the C function is 83.22% better than the original value. As can be seen for the three responses considered, lSE, and C, the solution obtained with the GLOBAL method is better than the IOM, with improvements of about 30% in three cases.

+,

+

Table 8. Optimum parameter profile obtained through the application of the IOM and the GLOBAL clustering method.

Other constraints. This is a set of constraints of different nature that forces the plant to operate within the suitable operative range. These constrains limit some flux ratios (i and ii); the residence times (iii), the settler rates (iv and v) and the ratio of organic matter input flux over the biomass (vi, vii, viii and ix).

Optimum solutions (X.)opJ(X.)bUol

Variables X 14 (VI)

IOM solution 1.0 0.282 X15 (VI) 0.495 Xl' (ad l ) X 17 (adz) 1.0 0.1 XII (fkl ) 0.1 Xl' (fk2) 0.1 Xl. (d) 10.52 Xli (k,) 0.43 ISE 1812.5 2242.5 C Improvement 83.7

Table 7. Other constrains -0.2808 <-0.2985 Y13 + YI4 < 0.2244 -1.336<-XI+O.038Y3+O.002Y4+O.042YwYI4 -0.003yu- 0.0015Y24 < 0.964 iii -1.44<-Y2+O.91Y3+O.082Y4-YIs+O.91yu +O.082Y23 < 0.33 iv 0.298 YI3 - YI6 < 0.509 v -X 17 + 0.8606 Yu+ 0.1395 Y23 < 0.529 vi -0.233 < 0.566 YI - 0.566 Yll +0.566 YI4 + 0.433 YI6 + 35.088 Yu - 1.360 Y24 < 0.289 vii -0.052<0.432Y2-O.432YI2+O.432YIS+O.56YI6 +0.35yu L36Y24 <0.47 viii -O.092<0.688yw 1.634Yu+O.264Y23+O.063Y24 < 0.163 ix -0.191 <-0.689yw33.45Yu-0.264Y23 + 1.2975Y24 < 0.176 11

+

GLOBAL 0.621 0.53 0.575 1.16 0.132 0.081 0.072 10.52 0.3986 1133.2 1538.205 88.8

It should be stated that while IOM optimizes over the

S-System model, the global methods optimize over the original kinetic modeL The S-System model version, derived from another kinetic model, implies a certain approximation (i.e. flux aggregation). The point is that after translating the S-System solution to the original kinetic model, we obtain a good solution quite similar to that provided by the global methods. However, with the GLOBAL method the

473

compu~tion time required was about 250 seconds (Fortran implementation, PC Pentium Ill). On the other hand the IOM approach can solve this type of problems in almost real time, since solving a LP can be done very fast and the symbolic manipulation needed by the IOM approach can be fully automated by means of suitable software (e.g. Maple). This means that we envision that the IOM procedure can be incorporated into a modern process simulator (e.g. gPROMS, Abacuss). By doing this, users will be able to define a bioprocess model and a tentative control structure and, almost immediately, obtain via the IOM the optimal integrated design ensuring low capital and operating costs plus maximum controlability. Further, the IOM scales up very well, since huge LPs can be solved very rapidly. In contrast, methods like GLOBAL do not scale well with problem size, with computational cost often increasing exponentially. Another advantage of the IOM refers to the optimization issue itself. Linear programming allows us to tackle multiobjective optimization problems that are not readily accessible by using direct, global optimization methods. In the present case we have approached the multiobjective optimization by transforming a multi-objective optimization problem into a mono-objective problem through the definition of a combined objective function that measures the plant's economics and its controllability. This approach has however some drawbacks since it constitutes a simplification that may mask important aspects of the problem. In this regard, the Multiobjective Indirect Optimization Method (Vera et al. 2003) is an alternative approach to this problem using S-system formalism. In this algorithm, we profit the properties of S-system models to solve a multi objective non-linear problem using multiobjective linear programming with low computational requirements.

G6rak, A., Kraslawski, A., and Vogelpohl, A. Simulation und Optimierung der MehrstoffRektifikation. Chem. Ing. Tech. 59: 95-106 (1987). Gutierrez, G., and Vega, P. Integrated design of activated sludge process taking into account the closed loop controllability. Proceedings of ESCAPE-IO. 63-69 (2000). Moles, e.G., Gutierrez, G., Alonso, A.A., and Banga, J.R. Integrated process design and control via global optimization: a wastewater treatment plant case study. Proceedings of European Control Conference (ECC). 4-7 September, 2001 , Porto (portugal) (2001) Takacs I., Patry G.G. and Nolasco D. A dynamic model of the clarification thickening process. Wat. Res 25(10):1263-1271. Torres, N.V., Voit, E.O., Glez-Alc6n, c., and Rodriguez, F. An indirect optimization method for biochemical systems: Description of method and application to the maximization of the rate of ethanol, glycerol, and carbohydrate production in Saccharomyces cerevisiae. Biotech. Bioeng. 55(5), 758-772 (1997). Vagners, J. . Optimization techniques. In: C.E. Pearson (Ed.): Handbook of Applied Mathematics. Selected Results and Methods, 2nd Ed., pp. 1140-1216 (1983). Vera, J., De Atauri, P., Cascante, M. and Torres, N.V. Multicriteria Optimization Of Biochemical Systems By Linear Programming. Application To The Ethanol Production By Saccharomyces Cerevisiae. Biotechnololy and Bioengineering. In Press. (2003) Voit, E.O. Computational Analysis of Biochemical Systems. A Practical Guide for Biochemists and Molecular Biologist., Cambridge University Press, Cambridge, U.K (2000).

4. ACKNOWLEDGEMENTS This work was supported by a research grant from the MCyT, n° BI02002-041 57-C02-02. JulioVera was the recipient of a research grant from the MCyT, ref. PN99-4320298.

5. REFERENCES Alvarez-Vasquez, F., Gonzatez-Alc6n, C.M. and Torres, N.V. Metabolism of citric acid production by Aspergillus niger. Model definition, steady state analysis, dynamic behavior and constrained optimization of citric acid production rate. Biotechnology Bioengineering. 70(1): 82-108 (2000). Csendes, T. Nonlinear parameter estimation by global optimization. Efficiency and reliability. Acta cibemetica 8, 361-370 (1988).

474