1948
Ind. Eng. Chem. Res. 1999, 38, 1948-1955
Optimization of Fed-Batch Reactors by the Luus-Jaakola Optimization Procedure Rein Luus* and Denis Hennessy Department of Chemical Engineering, University of Toronto, Toronto, Ontario M5S 3E5, Canada
The presence of numerous local optima makes optimization of fed-batch reactors a challenging problem. In addition, the low sensitivity of the performance index on the control policy makes it very difficult to determine the optimal control policy accurately. As an alternative to iterative dynamic programming, we consider the application of the direct search optimization procedure of Luus and Jaakola (AIChE J. 1973, 19, 760-766) with the recent refinements involving the use of a quadratic penalty function with a shifting term to handle equality constraints. In using the optimization procedure in multipass fashion, we use the extent of variation of the variables in a pass to provide the region size for the beginning of the subsequent pass. In establishing the optimal control policies for three fed-batch reactor models, we show that this approach provides a good way of solving such optimization problems. 1. Introduction For design and operation of chemical engineering systems, it is important to determine the maximum possible economic benefit that can be realized under physical constraints. In fed-batch reactors, we are concerned with determining the feed rate to the reactor that will give the maximum amount of the desired product without violating any constraints. Although having a single control variable in the form of the feed rate may appear to present a simple optimal control problem, considerable difficulties have been reported in the determination of the optimal feed-rate policy for fedbatch reactors.2-6 Advantages of using iterative dynamic programming7 for establishing the optimal control of fed-batch reactors were pointed out by Hartig et al.8 and by Bojkov and Luus.9 Because there is always a possibility of obtaining a local optimum, it is useful to crosscheck results with a different optimization procedure. One such means is to use direct search optimization, where the optimal control policy is obtained by optimizing the policy over all of the time stages simultaneously, rather than using stage to stage optimization as in iterative dynamic programming (IDP). The availability of very fast computers at low cost makes such a procedure feasible for solving very difficult optimal control problems such as the cancer chemotherapy scheduling problem.10 By using the Luus-Jaakola (LJ) direct search optimization procedure1 in a multipass fashion with the search region determined according to the extent of the change in the variables during a pass, a nonseparable optimal control problem involving three state variables and three control variables over 100 stages was successfully solved by Luus11 as a 300-variable nonlinear optimization problem. Using the same approach, very accurate parameter estimates for the Lotka-Volterra problem were obtained.12 Because the LJ optimization procedure exhibits a high reliability of obtaining the global optimum13 and it is easy to program, the goal of * To whom correspondence should be addressed. E-mail:
[email protected]. Fax: 416-978-8605. Phone: 416-9785200.
this paper is to examine the viability of this approach for establishing the optimal control of fed-batch reactors. 2. Models of Typical Fed-Batch Reactors In fed-batch reactors, the reactors are operated in fedbatch mode, where the feed rate into the reactor is used for control. Because there is no outflow, the feed rate must be chosen so that the batch volume does not exceed the physical volume of the reactor. There is also an upper limit for the feed rate. For the present work we consider models of three fed-batch reactors that have been used recently for optimal control studies. Model 1. For the first model we choose the fermentation process involving ethanol fermentation as formulated and used for optimal control by Hong,14 and used for optimal control studies by Chen and Hwang,15 Luus,5 Hartig et al.,8 and Bojkov and Luus.9 The differential equations describing the ethanol fermentation are
x1 dx1 ) Ax1 - u dt x4
(1)
(
)
150 - x2 dx2 ) -10Ax1 + u dt x4
(2)
x3 dx3 ) Bx1 - u dt x4
(3)
dx4 )u dt
(4)
with
A) B)
(
(
)( )(
)
x2 0.408 1 + x3/16 0.22 + x2
(5)
)
x2 1 1 + x3/71.5 0.44 + x2
(6)
where x1 represents the cell mass concentration, x2 is the substrate concentration, x3 is the desired product
10.1021/ie980731w CCC: $18.00 © 1999 American Chemical Society Published on Web 04/10/1999
Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999 1949
concentration, and x4 is the volume of the liquid inside the reactor. The initial state is specified as
xT(0) ) [1 150 0 10]
give a value reasonably close (within 0.4%) to the optimal performance index.23 The bioreactor is described by five differential equations
(7)
with the feed rate to the reactor u constrained by
0 e u e 12
(8)
The liquid volume is limited by the 200 L vessel size, so that at the final time tf
h ) x4(tf) - 200 e 0
(9)
The performance index to be maximized by choosing u is the total mass of the desired product at the final time tf:
I ) x3(tf) x4(tf)
(13)
dx2 u ) g2x3 - x2 dt x5
(14)
dx3 u ) g3x3 - x3 dt x5
(15)
dx4 u ) -7.3g3x3 + (20 - x4) dt x5
(16)
dx5 )u dt
(17)
(10)
In the original formulation of the problem, the final time tf is not specified. Luus5 showed that with 20 time stages of equal length the maximum yield can be obtained with tf ) 62.9 h. By using flexible stage lengths, Bojkov and Luus9 showed that tf should be in the range 61.1-61.3 h for maximum yield. Hartig et al.8 by using Hong’s equations14 concluded that the optimum batch time is 62 h but carried out a comparison of optimization procedures with tf ) 63.0 h. Therefore, we consider here the fed-batch optimization problem with tf ) 63.0 h to provide direct comparison to the results reported in the literature. For the maximum amount of desired product, it is expected that the final volume must be at its upper limit. Therefore, the inequality expressed by eq 10 becomes equality, and we formulate the augmented performance index to be maximized as 2
J ) x3(tf) x4(tf) - θ(h - s)
where
g3 )
21.87x4 (x4 + 0.4)(x4 + 62.5)
(18)
x4e-5x4 g2 ) 0.1 + x4
(19)
4.75g3 0.12 + g3
(20)
g1 )
with the initial condition
x(0) ) [0 0 1 5 1]T
(21)
(11)
where a penalty term is added to the performance index in eq 10. In the penalty function we have included the shifting term s. The advantages of such a penalty function including a shifting parameter to handle equality constraints have been recently demonstrated.7,16-18 The shifting parameter s is initially set to zero and then is updated after every pass of the optimization procedure according to
sq+1 ) sq - (x4(tf) - 200)
dx1 u ) g1(x2 - x1) - x1 dt x5
(12)
where q denotes the pass number in the multipass optimization procedure. At the optimum, the factor -2θs gives the Lagrange multiplier that yields the sensitivity information, showing the effect of changing the volume constraint x4(tf) ) 200. Model 2. A very challenging optimal control problem is the optimization of the fed-batch reactor model as presented by Park and Ramirez.3 This model was used by Luus19,20 to show early applications of iterative dynamic programming to singular control problems and by Yang and Wu21 to test their extension of simulated annealing. The difficulties in establishing the optimal control policy are partly due to the low sensitivity of the performance index on the control policy.22 This low sensitivity can be used to one’s advantage, because a crude approximation to the optimal control policy can
Here x1 denotes the concentration of secreted SUC2-s2 in culture and x2 denotes the concentration of total SUC2-s2 in culture, both expressed in arbitrary units; x3 is the culture cell density (g/L), x4 is the culture glucose level (g/L), and x5 denotes the culture volume (L). There is no upper limit imposed on the culture volume, but the feed rate to the reactor, which is the control u, is bounded by
0eue2
(22)
The performance index to be maximized is the product of the concentration of the protein x1 and the culture volume x5 at the final time tf
I ) x1(tf) x5(tf)
(23)
where the final time tf is specified as 15 h. Model 3. A very challenging optimal control problem is determining the optimal feed rate to a fed-batch fermentor in which the maximum production of penicillin is required, as studied by Lim et al.2, Cuthrell and Biegler,4 Luus,6 Dadebo and McAuley,24 Gupta,25 and Mekarapiruk and Luus.26 The system is described by four differential equations
1950 Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999
( ) ( )
dx1 x1 u ) h1x1 dt 500x4
(24)
x2 dx2 u ) h2x1 - 0.01x2 dt 500x4
(25)
(
)
dx3 h1x1 h2x1 0.029x3x1 x3 u )+ 1dt 0.47 1.2 0.0001 + x3 500 x4 (26) dx4 u ) dt 500
(27)
where
h1 )
0.11x3 0.006x1 + x3
(28)
{ {
x1 - 40 dx5 ) 0 dt 1000
if x1 > 40 if 0 e x1 e 40 if x1 < 0
(37)
x3 - 25 dx6 ) 0 dt 1000
if x1 > 25 if 0 e x3 e 25 if x3 < 0
(38)
and then form the augmented performance index to be maximized
J ) I - θ1(x4(tf) - 10 - s)2 - θ2x5(tf) - θ3x6(tf) (39) For flexibility we have chosen here different penalty function factors for each constraint, but, as is shown in the numerical results, nothing is lost if the factors are taken to be equal. As in model 1, the shifting parameter s is initially put to zero and is then updated according to
and
h2 )
0.0055x3 0.0001 + x3(1 + 10x3)
sq+1 ) sq - (x4(tf) - 10) (29)
The first three state variables are the concentration of biomass, product, and substrate, respectively, and x4 denotes the volume. The initial state is given as
x(0) ) [1.5 0 0 7]T
(30)
The constraints on the feed rate are
0 e u e 50
(31)
The constraints on the state variables are
0 e x1 e 40
(32)
0 e x3 e 25
(33)
0 e x4 e 10
(34)
and
The performance index to be maximized is the total amount of penicillin produced in the batch time tf ) 132 h, given by
I ) x2(tf) x4(tf)
(35)
The state constraints imposed by eqs 32-34 make the optimization problem very difficult. As suggested by Mekarapiruk and Luus,26 because to get the maximum amount of penicillin, the liquid volume must be at its maximum value, we replace the inequality in eq 34 by the final state equality constraint
x4(tf) - 10 ) 0
(36)
To handle inequality constraints in eqs 32 and 33, we use a penalty function approach as suggested by Luus19 and Vassiliadis et al.27 and used by Mekarapiruk and Luus26 for optimal control of this fed-batch reactor. We first introduce the auxiliary variables through the differential equations
(40)
where q denotes the pass number in the optimization procedure. At the optimum the product -2θ1s gives the sensitivity of the performance index to the volume constraint and tells us by how much the performance index can be increased if the volume constraint is relaxed. 3. LJ Optimization Procedure Let us outline the LJ direct search optimization procedure1 for this set of problems where stages of equal length are used and where the final time is specified. The direct search optimization is straightforward: 1. Divide the time interval into P stages, each of equal length L ) tf/P. 2. Take an initial estimate for the control u at each stage, accumulating the P values in the P × 1 vector u, and select the initial region size for the control used at the stages denoted by rin. Choose the number of iterations N to be used in each pass, the number of passes, and the region contraction factor γ. 3. Set the pass index q to 1. 4. Set the iteration index j to 1, and put rj ) rinq. 5. Select R sets of values for the variables
uj ) u*j + Drj
(41)
where D is a diagonal matrix with random numbers between -1 and 1 as diagonal elements. Each set consists of different random numbers. 6. For each of the R sets, evaluate the appropriate performance index. 7. Compare the values of the performance index, and save the set of variables u that yields the maximum value for the performance index; denote this by u*j. 8. Reduce the size of the region by the contraction factor γ
rj+1 ) γrj
(42)
9. Increment the iteration index j by 1, go to step 5, and continue for the specified number of iterations. This concludes one pass of the optimization procedure. 10. Choose a new initial region size for the next pass. For the first few passes, choose a new region size to be a fraction η of the initial region size used in the present
Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999 1951
pass. After a small number of passes (such as 10), choose the initial region size to be equal to the change in the variables during the current pass that yielded the best performance index, i.e. q+1 rin,i ) |uj*(N) - ui*(1)|, i ) 1, 2, ..., P
Table 1. Number of Passes, Each Consisting of 20 Iterations, Required To Reach within 0.001% of the Optimum (First Number) and the Optimum I ) 20 746.69 (Second Number) with P ) 10 Stages as a Function of the Number of Random Points R and the Penalty Function Factor θ for Model 1
(43)
*(N)
denotes the best value for ui at the end of where ui the current pass consisting of N iterations and ui*(1) denotes the best value for ui at the beginning of the current pass. 11. Update the shifting term s according to eq 12 for model 1 and according to eq 40 for model 3. 12. Increment the pass number index q by 1, go to step 4, and continue until the total number of passes has been reached or some convergence criterion is satisfied. 13. Interpret the results. To improve the probability of getting to the vicinity of the global optimum, for the first 10 passes, instead of using eq 43, the region size was restored to η times the region size used at the beginning of the previous pass. Here η was chosen to be 0.70. To avoid region collapse, if a particular region becomes zero, because there was no change in the parameter during the pass, then arbitrarily a value of 10-5 is assigned to that region size at the beginning of the subsequent pass.
θ ) 0.5
θ ) 1.0
50 100 200 300 400 500
a 17, 37 24, 33 13, 17 15, 24 15, 27
17, 40 a 13, 20 16, 16 14, 29 15, 19
s
-102.989
-51.495
R
a
All calculations were done in double precision on Pentium/120 and Pentium/166 personal computers. Throughout we used γ ) 0.95 as the region contraction factor for reducing the region size after every iteration. For each pass we chose N ) 20 iterations. To obtain an accurate integration, the subroutine DVERK28 was used with a local error tolerance of 10-6 for models 1 and 3 and of 10-7 for model 2. Model 1 with tf ) 63.0 h. The initial value for control at each time stage was chosen to be 3.0. The initial region size in the first pass was chosen to be equal to 12 and for the second pass the initial region size was chosen as 8.4 (i.e., η × 12), continuing until the 10th pass For the 11th pass and all subsequent passes, the initial region size was determined according to the extent of variation of the variables, as suggested by Luus11 and given in the above algorithm in eq 43. Because the difficulty of obtaining the global optimum increases as the number of stages P is increased, we consider three cases. P ) 10 Stages. By using 400 different starting points, Hartig et al.8 obtained convergence to the global optimum with successive quadratic programming (SQP) 14 times, showing a success rate of 3.5%. With the present algorithm, the success rate was 75% for the cases tried. As is shown in Table 1, the penalty function factor θ is not very important if at least 100 random points are chosen per iteration. In about half of the cases, fewer than 16 passes were required to get to within 0.001% of the optimum (I ) 20 746.48), and fewer than 25 passes were required to reach the optimum (I ) 20 746.69) to within seven figures. The last line in Table 1 gives the shifting parameter s at the optimum. It is seen that the product θs is constant (-51.5) over the entire range of θ that was used. This means that if the volume constraint were relaxed by 0.01, then the performance index would be increased by 1.03, so that a volume constraint of 200.01
θ ) 100
a 16, 22 a 14, 20 15, 22 15, 22
a 16, 36 14, 20 16, 25 a 15, 24
-5.1495
-0.514 94
Denotes the local optimum.
Table 2. Number of Passes, Each Consisting of 20 Iterations, Required To Reach within 0.001% of the Optimum (First Number) and the Optimum I ) 20 832.1 (Second Number) with P ) 15 Stages as a Function of the Number of Random Points R and the Penalty Function Factor θ for Model 1 R 50 100 200 300 400 500 s
4. Numerical Results
θ ) 10
a
θ ) 0.5
θ ) 1.0
θ)2
θ)5
θ ) 10
θ ) 100
19, 19 15, 19 28, 29 22, 24 16, 16 18, 18
23, 65 19, 21 21, 24 16, 29 16, 16 17, 19
28, 30 20, 22 20, 21 15, 16 a 16, 20
60, 85 a 17, 18 18, 20 17, 17 17, 19
a 37, 55 a a a 38, 39
a 33, 37 a a 28, 36 20, 21
-106.185 -53.086 -26.521 -10.614 -5.3083 -0.53039
Denotes the local optimum.
will give a performance index of I ) 20 746.69 2(-51.5)(0.01) ) 20 747.72. To illustrate this aspect, a computer run was performed, where the final volume constraint in eq 9 was changed to 200.01. The new optimum of I ) 20 747.717 was obtained. P ) 15 Stages. The global optimum was again easily obtained. However, as is shown in Table 2, when the penalty function factor θ is increased beyond 5, the global optimum is more difficult to reach. With θ ) 10 it takes 38 passes to reach within 0.001% of the optimum and 39 passes to reach the optimum (I ) 20 832.1), whereas with θ ) 0.5, even with R ) 50, the global optimum is reached within 19 passes. With SQP Hartig et al.8 obtained the global optimum only three times from 400 starting points. The computation time for 50 passes, each consisting of 20 iterations, with R ) 100 was 815 s, which is about double the time required for iterative dynamic programming by Hartig et al.8 P ) 20 Stages. With 20 stages, Hartig et al.8 with SQP found 34 local optima, the best being I ) 20 805, but could not get the global optimum of I ) 20 841.1. To obtain the global optimum, they used iterative dynamic programming with 23 grid points and between 21 and 41 allowable values for control. They were able to get the global optimum in 18 out of 28 runs. The computation time (converted to Pentium/120 and WATCOM compiler) for a typical run was about 20 min. Here, as is shown in Table 3, the overall success rate is only 25%, but for θ ) 0.5 and 1, the success rate is 40%. The computation time for R ) 200 for 100 passes, each consisting of 20 iterations, was 3453 s on Pentium/120, and that for R ) 400 for 50 passes was 3710 s. Although iterative dynamic programming is computationally more efficient for P ) 20, the LJ optimization procedure still provides a very good way of crosschecking the results.
1952 Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999 Table 3. Number of Passes, Each Consisting of 20 Iterations, Required To Reach within 0.001% of the Optimum (First Number) and the Optimum I ) 20 841.1 (Second Number) with P ) 20 Stages as a Function of the Number of Random Points R and the Penalty Function Factor θ for Model 1 θ ) 0.5
θ ) 1.0
θ)2
θ)5
100 200 300 400 500
a 79, 94 a 34, 47 a
a a a 19, 33 32, 39
a a a 19, 23 a
a a a 19, 25 84, 89
s
-107.972
-54.018
-26.993
-10.797
R
a
Denotes the local optimum.
Table 4. Optimal Control Policies for tf ) 63 h for Model 1 Obtained by the LJ Optimization Procedure stage no.
P ) 10, I ) 20 746.69
P ) 15, I ) 20 832.1
P ) 20, I ) 20 841.1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
0 0 1.149 21 1.501 63 2.186 39 3.081 36 4.471 87 5.878 88 11.889 39 0
0 0 0 0.880 53 1.231 93 1.523 21 1.935 69 2.443 28 3.089 15 3.888 86 4.918 24 6.253 53 7.504 24 11.569 45 0
0 0 0 0 0.767 00 1.119 70 1.278 57 1.541 24 1.835 41 2.180 98 2.610 97 3.102 43 3.694 81 4.399 76 5.237 58 6.236 95 7.359 36 8.934 72 10.017 99 0
The optimal control policies with tf ) 63 h obtained by LJ optimization procedure as outlined here for these three cases are given in Table 4. Model 2 with tf ) 15.0 h. To obtain accurate integration of the differential equations, we used the subroutine DVERK28 with a local error tolerance of 10-7. The initial value for the control was taken as 1.0, and the initial region size for the first pass was taken as 2.0. The region contraction factor γ ) 0.95 was used after every iteration, and 20 iterations were used in every pass. For the first 10 passes, the region was restored to 0.70 of the value used at the beginning of the previous pass (η ) 0.70), and starting with the 11th pass, the algorithm given in this paper was used. If the initial region size at the beginning of any pass was calculated to be less than 10-5, then the region size was arbitrarily put equal to 10-5. This prevented the region from collapsing for any particular variable. P ) 31 Stages. When the time interval was divided into 31 stages of equal length, the optimization problem involves a 31-dimensional search. With the given algorithm no difficulties were experienced in getting the maximum value of I ) 32.613 392 for the performance index, giving the control policy in Figure 1. As can be seen in Table 5, the use of 100 random points per iteration is adequate, and the optimum is obtained within 43 min of computation time on the Pentium/120 digital computer. It was also interesting to note that with R ) 400 a local optimum with I ) 32.418 11 was obtained, so the global optimum is not always obtained if a very large number of points are used.
Figure 1. Optimal control policy for model 2 with P ) 31, yielding I ) 32.613 392. Table 5. Computational Effort Required by the LJ Optimization Procedure for Model 2 with P ) 31 To Reach I ) 32.613 392 no. of random points per iteration, R
no. of passes required to reach I ) 32.613 392
computation time, s on a Pentium/120
100 200 300 400 500
143 127 75 a 72
2570 4785 4571 a 4762
a
Denotes the local optimum.
Table 6. Computational Effort Required by Iterative Dynamic Programming for Model 2 with P ) 31 To Reach I ) 32.613 392 no. of random points per iteration, R
no. of passes required to reach I ) 32.613 392
computation time, s on a Pentium/120
3 5 7 9 11 13
a 36 44 39 38 35
1527 1187 2916 3330 3952 4279
a
After 50 passes, I ) 32.613 380.
The results obtained with the use of iterative dynamic programming (using a single grid point with γ ) 0.95, η ) 0.85, with 40 iterations per pass, an initial control of 0.1, and an initial region size of 1.0) are given in Table 6. It is seen that with R ) 5 the computation time to reach I ) 32.613 392 is 20 min and with R ) 7 it is 49 min. Therefore, the computation times for the two methods are quite comparable. P ) 45 Stages. When P ) 45 was used, a more accurate switching was obtained in the control policy, which resulted in a better performance index of I ) 32.686 867. It is interesting to note that this value of the performance index is marginally better than can be obtained with P ) 80 stages (I ) 32.686 016), showing that the performance index is very sensitive to the time of switching. Although 80 stages give a more refined and smoother control policy, this refinement does not make up for the inability of switching at the “proper” moments. The optimal control policy is given in Figure
Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999 1953
Figure 2. Optimal control policy for model 2 with P ) 45, yielding I ) 32.686 867.
Figure 4. Convergence profile for IDP showing the deviation from the optimum performance index (32.686 867 - I) as a function of the pass number: (b) R ) 7, (2) R ) 9, (9) R ) 13. Table 7. Performance Index Obtained after 100 Passes with P ) 15 for Model 3 with tf ) 132 h, Showing the Effects of the Number of Random Points R and the Penalty Function Factor θ θ)5
50 100 150 200 250 300 400
87.668 87.957 87.920 87.970 87.957 87.920 87.970
87.956 87.957 87.970 87.782 87.970 87.970 87.957
87.331 87.970 87.956 87.472 87.957 87.913 87.970
87.920 87.934 87.924 87.472 87.970 87.956 87.920
-2.8846
-1.4422
-0.7211
-0.2885
s
Figure 3. Convergence profile for the LJ optimization procedure showing the deviation from the optimum performance index (32.686 867 - I) as a function of the pass number: (b) R ) 100, (2) R ) 200, (9) R ) 300.
2. To get this optimal control policy by the LJ optimization procedure required a large number of passes, as is shown in Figure 3. With the use of R ) 100, after 500 passes the performance index I ) 32.686 866 was reached, whereas I ) 32.686 867 was reached in 350 passes with R ) 300 and in 400 passes with R ) 200. A large number of figures in the performance index is required to yield a relatively smooth control policy, because of the very low sensitivity of the control policy on the performance index. This very low sensitivity was pointed out by Luus,22 who showed that with P ) 80 stages the control policy with I ) 32.686 02 was relatively smooth, whereas the control policy giving I ) 32.686 01 exhibited oscillatory sections. The convergence profile for IDP is shown in Figure 4. It is seen that the number of random points per iteration is not that important. To obtain an accurate optimal control policy, about 50 passes, each consisting of 40 iterations, are required. With R ) 7, the computation time for 50 passes was 3584 s, and with R ) 13, the computation time for 50 passes was 7391 s on a Pentium/120. The computation time with R ) 300 with the LJ optimization procedure for 500 passes, each consisting of 20 iterations, was 36 968 s, so IDP becomes
penalty function factor θ ) 10 θ ) 20
no. of random points
θ ) 50
noticeably faster relative to the LJ optimizaton procedure as the number of stages is increased. This effect is expected, because in IDP a sequential 1-dimensional optimization is carried out over the stages, but in the LJ optimization procedure, a simultaneous optimization involving a P-dimensional search is required. Model 3 with tf ) 132.0 h. To obtain sufficient accuracy in the integration of the six differential equations, we used the subroutine DVERK28 with a local error tolerance of 10-6. The problem was solved with the use of P ) 10, 15, and 20 stages. The initial control policy u(k) ) 11.9, k ) 0, 1, ..., P - 1, was used in each case, and the initial region size for the first pass of 10.0 was used. We took γ ) 0.95 and η ) 0.70 with 20 iterations per pass. After the 10th pass, the region sizes were obtained from the extent of the variables during the previous pass as outlined in the algorithm in eq 43. If the region size calculated was less than 10-5, then the value 10-5 was used, as in the previous examples. The same penalty function factor θ was used for all of the constraints; i.e., θ ) θ1 ) θ2 ) θ3. When P ) 10 was used, the maximum value of the performance index I ) 87.9164 was readily obtained when θ was chosen in the range 5-50 and R ) 100 random points per iteration were used. With R ) 37 the local optimum having the performance index I ) 87.8165 was obtained. P ) 15 Stages. With P ) 15 stages, numerous local optima were found, as is shown in Table 7. The global optimum I ) 87.970 was obtained 8 times out of 28 runs, and the second best local optimum I ) 87.957, which is 0.015% less than the global optimum, was obtained 7 times. Other local optima were obtained less frequently.
1954 Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999 Table 8. Control Policies with tf ) 132 h for Model 3 Obtained by the LJ Optimization Procedure, Yielding the Global Optimum and the Second Best Local Optimum P ) 10
P ) 15
P ) 20
stage no.
I ) 87.9164
I ) 87.8165
I ) 87.9699
I ) 87.9572
I ) 87.9964
I ) 87.9959
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
4.9327 36.5413 8.2841 8.6942 8.9132 9.0514 9.1611 9.2597 9.3537 9.4451
3.3850 7.3479 36.6555 9.0184 9.2464 9.3885 9.5014 9.6025 9.6987 9.7920
3.7421 6.9730 13.1833 44.6235 8.6796 8.9169 9.0584 9.1608 9.2439 9.3170 9.3848 9.4494 9.5123 9.5744 9.6352
4.3298 9.1617 48.6433 8.2643 8.5211 8.7570 8.8737 8.9950 9.0641 9.1528 9.2192 9.2794 9.3388 9.4004 9.4542
3.5099 5.7735 9.4129 50.0000 21.2620 8.4138 8.6887 8.8430 8.9509 9.0339 9.1035 9.1639 9.2188 9.2700 9.3204 9.3679 9.4153 9.4625 9.5084 9.5538
3.4617 5.6719 9.1663 21.2892 50.0000 8.4611 8.7170 8.8693 8.9757 9.0577 9.1241 9.1862 9.2410 9.2925 9.3418 9.3900 9.4372 9.4837 9.5295 9.5748
Table 9. Performance Index Obtained after 100 Passes with P ) 20 for Model 3 with tf ) 132 h, Showing the Effects of the Number of Random Points R and the Penalty Function Factor θ penalty function factor θ ) 10 θ ) 15
no. of random points
θ)5
50 100 150 200 250 300 400
87.8801 87.5321 87.6255 87.9218 87.9344 87.9594 87.9964
87.9964 87.6177 87.9964 87.9959 87.8737 87.9959 87.4467
87.7516 87.9336 87.9220 87.5326 87.7280 87.9964 87.7514
s
-2.8906
-1.4454
-0.9635
There appears to be no clear pattern with respect to the number of random points used per iteration and the probability of getting the global optimum. There is a wide range for the penalty function factor θ for which the global optimum is obtained, but it should not be chosen larger than 50. The range 5 e θ e 20 appears to be good for this problem. The control policies corresponding to I ) 87.970 and to I ) 87.957, as given in Table 8, show that the control policies are quite different yet give a difference of only 0.015% in the performance index. P ) 20 Stages. Of special interest is the optimal control policy with P ) 20, giving I ) 87.9964, because it is quite different from the optimal control policies reported in the literature.4,6,25,26 The second best local optimum with I ) 87.9959 is within 0.001% of the global optimum but is quite different where the maximum feed rate occurs at stage 5 rather than at stage 4, as is shown in Table 8. The third best local optimum with I ) 87.9594 corresponds to the literature value that had been thought to be the global optimum. Although the improvement in the performance index over the best reported value is only 0.04%, it is clear that the global optimal control policy had been missed in the past. This points out the need to have a different optimal control method to cross-check the results. In highly nonlinear and constrained systems where there are numerous local optima, it is easy to miss the global optimum if local optima give performance indices that are close to the global optimum.
Figure 5. Trajectories of the three constrained state variables.
As is shown in Table 9, the global optimum is not easily obtained with P ) 20. Of the 21 runs, the global optimum was obtained 4 times. With R ) 50 and θ ) 10, the global optimum was obtained in 81 passes requiring 17 726 s of computation time on a Pentium/ 166 digital computer. At the optimum, the Lagrange multiplier is given by the product 2θs ) -28.9. This means that since the final volume is 10.000 00, the maximum error related to reporting the maximum performance index is (5 × 10-6)28.9 ) 1.5 × 10-4. In Figure 5, obtained by integrating the system equations with 200 time steps, it is seen that all of the constrained state variables are inside the constraint boundaries. In fact, the maximum value of x3 is 22.81, occurring at t ) 26.6 h, and is well below the upper constraint 25.0. 5. Conclusions The proposed method of determining the region sizes at the start of iterations in every pass enables the LJ optimization procedure to establish the optimal control policy very accurately. The method therefore extends the usefulness of the LJ optimization procedure for determining the optimal control of complex systems. The computation time required for the optimal control of
Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999 1955
typical fed-batch reactors is quite reasonable, so the approach outlined here provides an attractive means of solving such optimal control problems if the number of stages is not too large. If the number of stages is greater than 15, then IDP is noticeably faster, but the present method still provides a very good means of crosschecking the results obtained by other methods, as was illustrated in the last example. Acknowledgment Financial support from the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged. Nomenclature D ) diagonal matrix with random numbers in the range -1 to 1 along the diagonal gi ) nonlinear functions defined in eqs 18-20 hi ) nonlinear function defined in eqs 28 and 29 h ) inequality written as an equality I ) performance index to be maximized J ) augmented performance index L ) length of a stage P ) number of stages q ) pass number r ) region size vector R ) number of randomly chosen points used per iteration s ) shifting term used in the quadratic penalty function t ) time, h tf ) batch time u ) control, feed rate into the reactor u ) vector consisting of control values for each stage Greek Letters γ ) region contraction factor by which the region size is reduced after every iteration η ) region restoration factor, used in the first 10 passes to determine the size of the initial region size in the LJ optimization procedure and for all passes with IDP θ ) penalty function factor Abbreviations A ) expression defined by eq 5 B ) expression defined by eq 6 LJ ) optimization method of Luus and Jaakola1 IDP ) iterative dynamic programming
Literature Cited (1) Luus, R.; Jaakola, T. H. I. Optimization by Direct Search and Systematic Reduction of the Size of Search Region. AIChE J. 1973, 19, 760-766. (2) Lim, H. C.; Tayeb, Y. J.; Modak, J. M.; Bonte, P. Computational Algorithms for Optimal Feed Rates for a Class of FedBatch Fermentation: Numerical Results for Penecillin and Cell Mass Production. Biotechnol. Bioeng. 1986, 28, 1408-1420. (3) Park, S.; Ramirez, W. F. Optimal Production of Secreted Protein in Fed-Batch Reactors. AIChE J. 1988, 34, 1550-1558. (4) Cuthrell, J. E.; Biegler, L. T. Simultaneous Optimization and Solution Methods for Batch Reactor Control Profiles. Comput. Chem. Eng. 1989, 13, 49-62. (5) Luus, R. Application of Dynamic Programming to Differential-Algebraic Process Systems. Comput. Chem. Eng. 1993, 17, 373-377. (6) Luus, R. Optimization of Fed-Batch Fermentors by Iterative Dynamic Programming. Biotechnol. Bioeng. 1993, 41, 599-602. (7) Luus, R. Iterative Dynamic Programming: From Curiosity to a Practical Optimization Procedure. Control Intell. Syst. 1998, 26, 1-8.
(8) Hartig, F.; Keil, F. J.; Luus, R. Comparison of Optimization Methods for a Fed-Batch Reactor. Hung. J. Ind. Chem. 1995, 23, 141-148. (9) Bojkov, B.; Luus, R. Optimal Control of Nonlinear Systems with Unspecified Final Times. Chem. Eng. Sci. 1996, 51, 905919. (10) Luus, R.; Hartig, F.; Keil, F. J. Optimal Drug Scheduling of Cancer Chemotherapy by Direct Search Optimization. Hung. J. Ind. Chem. 1995, 23, 55-58. (11) Luus, R. Determination of the Region Sizes for LJ Optimization Procedure. Hung. J. Ind. Chem. 1998, 26, 281-286. (12) Luus, R. Parameter Estimation of Lotka-Volterra Problem by Direct Search Optimization. Hung. J. Ind. Chem. 1998, 26, 287-292. (13) Wang, B. C.; Luus, R. Reliability of Optimization Procedures for Obtaining Global Optimum. AIChE J. 1978, 24, 619626. (14) Hong, J. Optimal Substrate Feeding Policy for Fed Batch Fermentation with Substrate and Product Inhibition Kinetics. Biotechnol. Bioeng. 1986, 27, 1421-1431. (15) Chen, C. T.; Hwang, C. Optimal Control Computation for Differential-Algebraic Process Systems with General Constraints. Chem. Eng. Commun. 1990, 97, 9-26. (16) Luus, R. Handling Difficult Equality Constraints in Direct Search Optimization. Hung. J. Ind. Chem. 1996, 24, 285-290. (17) Luus, R.; Storey, C. Optimal Control of Final State Constrained Systems. Proceedings of IASTED International Conference on Modelling, Simulation and Control, Singapore, Aug 1113, 1997; IASTED/Acta Press: Anaheim, CA, 1997; pp 245-249. (18) Luus, R.; Mekarapiruk, W.; Storey, C. Evaluation of Penalty Functions for Optimal Control. Proceedings of the International Conference on Optimization Techniques and Applications ICOTA’98, Perth, Australia, July 1-4, 1998; Curtin Printing Services: Curtin, Australia, 1998; pp 724-731. (19) Luus, R. Application of Iterative Dynamic Programming to State Constrained Optimal Control Problems. Hung. J. Ind. Chem. 1991, 19, 245-254. (20) Luus, R. On the Application of Iterative Dynamic Programming to Singular Optimal Control Problems. IEEE Trans. Autom. Control 1992, 37, 1802-1806. (21) Yang, R. L.; Wu, C. P. Global Optimal Control by Accelerated Simulated Annealing. Paper presented at the First Asian Control Conference, Tokyo, Japan, 1994. (22) Luus, R. Sensitivity of Control Policy on Yield of a FedBatch Reactor.Proc. IASTED Int. Conf. on Modelling and Simulation, Pittsburgh, PA, April 27-29, 1995, 224-226. (23) Banga, J. R.; Irizarry-Rivera, R.; Seider, W. D. Stochastic Optimization for Optimal and Model-Predictive Control. Comput. Chem. Eng. 1998, 22, 603-612. (24) Dadebo, S. A.; McAuley, K. B. Dynamic Optimization of Constrained Chemical Engineering Problems Using Dynamic Programming. Comput. Chem. Eng. 1995, 19, 513-525. (25) Gupta, Y. P. Semiexhaustive Search for Solving Nonlinear Optimal Control Problems. Ind. Eng. Chem. Res. 1995, 34, 38783884. (26) Mekarapiruk, W.; Luus, R. Optimal Control of Inequality State Constrained Systems. Ind. Eng. Chem. Res. 1997, 36, 16861694. (27) Vassiliadis, V. S.; Sargent, R. W. H.; Pantelides, C. C. Solution of a Class of Multistage Dynamic Optimization Problems. 2. Problems with Path Constraints. Ind. Eng. Chem. Res. 1994, 33, 2123-2133. (28) Hull, T. E.; Enright, W. D.; Jackson, K. R. User Guide to DVERKsa Subroutine for Solving Nonstiff ODE’s; Report 100; Department of Computer Science, University of Toronto: Toronto, Canada, 1976.
Received for review November 19, 1998 Revised manuscript received March 2, 1999 Accepted March 3, 1999 IE980731W
