Simultaneous optimization of dynamic bioprocesses

SIMULTANEOUS OPTIMIZATION OF DYNAMIC BIOPROCESSES

C.A.M.Riascos and J.M.Pinto^*

Department of Chemical Engineering, University of São Paulo

Av. Prof. Luciano Gualberto, t.3, 380 - São Paulo SP 05508-900 - Brazil.

E-mail: [e-mail]

(Received: March 5, 2002 ; Accepted: May 27, 2002)

Abstract - It is known that the optimal control may introduce significant economical benefits into production processes, thus being an important and challenging research area with practical relevance. The modeling and optimization of biotechnological processes has been object of research and their related results have generated improvements in operating conditions and strategies, however, the inherent features of dynamical bioprocesses prevent the application of conventional optimization algorithms, hence making necessary the development of tailored methods and strategies. The objective of this work is to develop mathematical programming strategies for simultaneous optimization of dynamic systems and evaluate their computational performance. Simultaneous optimization with orthogonal collocation is applied to a simplified model for biosynthesis of penicillin from glucose, which was studied by Cuthrell and Biegler (1989). The results show that discretization of differential equations systems (DAE) by orthogonal collocation in finite elements efficiently transforms dynamic optimization problems into nonlinear programming (NLP) problems, enabling to solve complex problems with several control variables and minimizing the approximation error.

Keywords: optimal control, orthogonal collocation, nonlinear programming.

INTRODUCTION

Dynamic processes represent a great variety of operations, in particular all processes in batch or fed-batch mode. Their optimal control requires special effort, due to the intrinsic dynamics and mainly by the presence of path constraints in the state and control variables. Moreover, adding the possible existence of discontinuities in the variable profiles and in the differential equations, a complex dynamical optimization problem is generated whose resolution strategies must be investigated. Cuthrell and Biegler (1987) discuss that the differential algebraic optimization problems (DAOP) cannot be solved by straightforward nonlinear programming (NLP) techniques or by optimal control methods.

Two popular approaches to solve DAOP are based on the discretization of variables. The first discretizes only the control variables, u(t), (sequential method or control vector parameterization) and the DAE systems is integrated using standard integration algorithms; the optimization is carried out in the space of the decision variables. The second approach relies on the discretization of all variables that converts the dynamic optimization problem into an NLP (simultaneous method); in this approach the optimization is carried out in the full space of the discretized variables, thus enabling the solution of problems with constraints on state and control variables.

The present work aims to develop and to evaluate mathematical programming techniques based on the orthogonal collocation method for the simultaneous optimization of dynamical processes, taking the fed-batch biochemical reactor as an example.

METHODOLOGY

Discretization and Approximation by Orthogonal Collocation on Finite Elements

Consider the optimal control problem (OCP) stated as follows

where F = objective function;

z^j(t) = state variable j profile;

u(t) = control profile;

h and g = equality and inequality constraints;

= initial value for state variable j;

NV = number of state variables;

U^L, U^U, Z^jL and Z^jU = control and state bounds.

To convert the OCP into an NLP, the differential algebraic equations (DAE) are discretized by orthogonal collocation on finite elements. Defining the normalized time q e [0, 1], yields:

where a_n = element limits and q_i = collocation point i. Note that a more compact notation is used in place of that in Cuthrell and Biegler (1987, 1989). The corresponding Lagrange polynomials are:

Thus, the state profile j is approached by a (K+1)^th degree polynomial, and the control by a degree K polynomial ( and u_n,i = polynomial coefficients). Functions j_i and y_i (polynomial building blocks) depend on the location of the collocation points that are roots of the Legendre polynomial (Cuthrell and Biegler, 1989; Rice and Do, 1995).

The derivatives of the Lagrange polynomials at the collocation points are given as follows

In (4), is a square matrix of dimension K+1 that is calculated from q_i (Rice and Do, 1995), and is the vector of values for variable j in the n^th-element. Then, the DAE can be written in the form of a residual for every point i and element n as follows

Constraints must be added to guarantee the continuity of the state profiles at the element limits. Therefore, the polynomials are extrapolated to generate the initial point of the next element. The control variable values on the element limits are obtained from the extrapolation of each approximation polynomial. With the orthogonal collocation method, the OCP turns into an algebraic nonlinear programming problem:

Control of the Approximation Error by Direct Enforcement

The algebraic approximation of the state profiles brings errors to the DAE resolution. Vasantharajan and Biegler (1990) present two strategies for error control: equidistribution and direct enforcement. The former reformulates the problem to find a solution with alternation in error sign, whereas in the latter, constraints are added to reduce the error in noncollocation points. The second approach is more straightforward and allows a direct control of the approximation error, but one must choose a sufficient number of elements to satisfy the error tolerance. Nevertheless, both formulations provide criteria for partitioning the process time on elements.

Denoting the noncollocation points as q_nc, the added constraints are as follows:

where x is the permissible error, R^j(t_n,nc) is the residue that corresponds to variable j at time t_n,nc, which represents the noncollocation point nc in the element n. The constant C (which penalizes the error by its proximity to the collocation points) is estimated as follows:

Differentiating the interpolation polynomials at the noncollocation points, and calculating the analytical derivatives with the DAE, the approximation residue is obtained according to:

THE PENICILLIN BIOSYNTHESIS PROBLEM

The penicillin biosynthesis from glucose, shown in Cuthrell and Biegler (1989) and previously studied by Lim et al. (1986), is developed in a fed-batch reactor with volume (V), biomass (X), product (P) and substrate (S). The control variable is the feed rate of substrate (U), and the concentration of the feed is constant (S_F). The problem is defined by:

where m(X,S) = biomass growth rate (h^-1) and r(S) = penicillin production rate (g P/g X h). The parameter values are shows in Table 1.

COMPUTATIONAL RESULTS

The methodology was applied to several cases, which are shown in Table 1. Table 2 shows the corresponding results and statistics that were obtained with GAMS/CONOPT (Brooke et al., 2000).

In cases CG1 and CG2 the end time is constant, in order to evaluate the capacity of global collocation for reproducing discontinuities in the control profile. Note that despite the objective value, F, being next to the analytical solution, the state and control profiles are far from the analytical ones (Figure 1). In case EF1, the final time and the element limits were fixed; the results and the profiles are very close to the analytical ones. These parameters were free in case EF2; in the optimization it occurred the collapse of two elements, which reduced the number of constraints, thus damaging the approach quality.

Observing the statistics (Table 2) it can be seen that the number of iterations and the computational effort increase with the number of collocation points and elements. In Figure 1 it is observed that the profiles are closer to the analytical solution when these parameters are increased.

Since there was a large increase on the processing time (t_f) with small effect in the objective function, productivity was maximized in cases EF3, CE1 to CE4. Note in Figures 1b and 1c that in case EF3 biomass continued to grow exponentially after 28h (when the growth stage finishes in the analytical solution) and the product concentration profile surpasses the analytical one.

To verify the need to control the approximation errors, these were estimated for case EF3. In Table 2 it is observed that the error values are significant. Another evidence is that the biomass profile surpasses the analytical solution (in the interval [0-28h]), which is unrealistic since the analytical solution attains maximum growth. Case CE1 presents error control with 0.1 absolute tolerance for all variables, while in case CE2 the tolerances were reduced to 1%. In CE2, despite the collapse of one element, the resulting control profile has similar features of those from the analytical solution (high feed rate, followed by reduction to zero and finally almost constant at an intermediate value).

Figure 2b shows that error control prevents the biomass concentration from surpassing the analytical profile (in the interval [0-28h]). Finally, in cases CE3 and CE4 the conditions of case CE2 were repeated, for different initial estimates, to verify the robustness of the approach. In Figure 2 and Table 2 it is observed that in case CE4 the profiles and the final results are identical to the ones of CE2. Nevertheless, case CE3 presents a different solution, thus implying that convergence to local solutions is an issue in such nonlinear and nonconvex optimal control problems.

Table 3 presents the results and statistics for the sequential approach and Figure 3 presents the optimal control and state profiles, which were obtained with MATLAB (Grace, 1995). Figure 3b-d show that the state profiles are very similar to the analytical ones, whereas the control profiles (Figure 3a) are not; furthermore, the objective function values for these cases (Table 3) are very near to the analytical one. From the statistics of Table 3 it can be seen that the number of iterations and the computational effort increase enormously with the number of intervals; in this Table the CPU times were normalized because the absolute values depend on the software employed.

CONCLUSION

The simultaneous approach for solving optimal control problems was investigated in this paper for a simplified model of a biorreator. The impossibility to estimate discontinuities in the control profile with global collocation motivates the incorporation of finite elements. Moreover, sequential strategies are unsuitable to solve complex problems due to the large increment in computational effort.

The approximation errors are significantly reduced when error enforcement is used and the increment of computational effort is small. The distribution of noncollocation points must cover the complete batch time to guarantee that the new constraints improve the complete approximation profile.

ACKNOWLEDGEMENTS

The authors would like to thank CAPES and PADCT/CNPq (grant 62.0239/97) for their financial support.

• Brooke, A., Kendrick, D., Meeraus, A. A. and Raman, R. (2000). GAMS - A User's Guide (Release 2.50). The Scientific Press. Redwood City, USA.
• Cuthrell, J.E. and Biegler, L.T. (1987). On the Optimization of Differential-Algebraic Process Systems. AIChE J., 33, (8), 1257-1270.
• Cuthrell, J.E. and Biegler, L.T. (1989). Simultaneous Optimization and Solution Methods for Batch Reactor Control Profiles. Comput. Chem. Engng., 13, (1), 49-62.
• Grace, A. (1995) Optimization Toolbox: for Ùse with Matlab, User's Guide. Math Works, Natick (MA).
• Lim, H.C., Tayeb, Y.J., Modak, J.M. and Bonte, P. (1986). Computational Algorithms for Optimal Feed Rates for a Class of Fed-Batch Fermentation. Biotech. Bioeng., 28, (9), 1408-1420.
• Rice, R.G. and Do, D.D. (1995). Applied Mathematics and Modeling for Chemical Engineers. Wiley, New York (NY).
• Vasantharajan, S. and Biegler, L.T. (1990). Simultaneous Strategies for Optimization of Differential-Algebraic Systems with Enforcement of Error Criteria. Comput. Chem. Engng., 14, (10), 1083-1100.

*

To whom correspondence should be addressed

Publication Dates

• Publication in this collection
  20 Jan 2003
• Date of issue
  Dec 2002

History

• Accepted
  27 May 2002
• Received
  05 Mar 2002