Semiexhaustive Search for Solving Nonlinear Optimal Control Problems

As reported in the literature, iterative dynamic programming (IDP) found better optimums for ... profiles using nonlinear programming. Luus (1989)has ...
0 downloads 0 Views 722KB Size
I n d . Eng. C h e m . Res. 1995,34,3878-3884

3878

Semiexhaustive Search for Solving Nonlinear Optimal Control Problems Yash P. Gupta* Department of Chemical Engineering, Technical University of Nova Scotia, Halifax, Nova Scotia B3J 2x4, Canada

It is difficult to obtain optimums, especially global optimums, for optimal control problems that are highly nonlinear and involve state constraints. The recommended practice in the literature

is to use different methods and then to implement the best solution found. To aid this process, a method that uses a semiexhaustive search to solve such problems is proposed in this paper. Its performance is presented on six example problems from the chemical engineering literature. As reported in the literature, iterative dynamic programming (IDP) found better optimums for several nonlinear problems. The proposed method not only found the IDPs optimums, but also, in some cases, it found even better optimums. 1. Introduction

+(XI 5

Nonlinear optimal control problems are frequently encountered in chemical engineering design and operation. The determination of an optimal control trajectory is important to maximize the yield of high valued final product. Although a number of methods are available for the solution of optimal control problems, the choice is limited when the equations are highly nonlinear. The problem becomes more difficult if state constraints are present andlor the optimal control profiles have singular arc portions. Recently, Cuthrell and Biegler (1989) have presented a method to convert these problems into a form which can be solved by using nonlinear programming. Chen and Hwang (1990) have proposed another method for finding piecewise constant optimal control profiles using nonlinear programming. Luus (1989) has proposed iterative dynamic programming (IDP) for solving such problems. Using IDP, better optimal solutions were found for several problems than were obtained using other methods. In some cases, IDP found a better optimum than was found earlier using IDP due to the difference in the way it was used. Since one cannot be sure if a global optimum has been reached for a problem of this type, one needs to cross-check the results using different optimization methods as recommended by Bojkov and Luus (1992). Although a global optimum still cannot be guaranteed, one can implement the best solution found by using the different methods. To support the above approach, a method that uses semiexhaustive search for solving nonlinear optimal control problems is proposed in this paper. 2. Optimal Control Problem

The following optimal control problem is considered in this study. The system is described by a set of nonlinear differential equations dxldt = flx,u,t)

(1)

where x is a state vector and u is a control vector bounded by asus/l

(2)

The initial state x(0) is given. There may be inequality constraints on the states of the type

* E-mail: [email protected]

0

(3)

The final time performance index is given by I[X(O),tfI= @(X(tf))

(4)

where the final time (tf)is given. The time interval (0, tf) can be divided into P time stages of equal length L where L = tflP. The optimal control problem is t o find a piecewise constant control trajectory over P time stages such that the performance index given in eq 4 is minimized. It is assumed that qNx(t))changes monotonically with time, which is a common occurrence. 3. Methodology

The proposed method is an iterative procedure that is started from an assumed control trajectory. A coarse grid consisting of piecewise constant paths over P stages is chosen around the assumed control trajectory, and a semiexhaustive search is conducted to seek the optimal trajectory in the grid. This completes one iteration. At each subsequent iteration, a narrower grid is formed around the optimal control trajectory found at the previous iteration by reducing the span (range) of the grid. The semiexhaustive search is then conducted in the narrower grid. The procedure is carried for a specified number of iterations, and this completes one pass. A number of passes may be made by restarting the search on a coarse grid around the optimal control trajectory found in the previous pass to see if the solution improves. The concept of finding an adequate solution by narrowing the grid around the optimal control trajectory found at the previous iteration is the same as used in the IDP method. However, the search method used to find an optimal trajectory in a given control grid is entirely different, and the grids on state variables used in IDP are not required. To explain the semiexhaustive search, let us consider a simplified case where we have only one control variable and four state variables and the optimal control trajectory is to be found over 10 time stages. Consider a coarse grid shown in Figure 1 where there are only two values of u on either side of the optimal trajectory found at the previous iteration. In other words, there are five values of u t o be considered over any stage. Let tf = 10 and the performance index t o be minimized be x&), that is, @(x(tf)) = x&) = Xd10).

0888-5885/95/2634-3878$09.OQIQ 0 1995 American Chemical Society

Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995 3879 ................. ..... ........ . . .. .. . . . ... . .. ....................

.

1

2

3

Stage Number I

1

2

3

4

5 6 7 Stage Number

8

9

10

Figure 1. Control grid over 10 time stages where five values of control are to be evaluated at each stage. The dark arrows show the piecewise constant optimal trajectory found at the previous iteration.

In an exhaustive search, we would start from do)and integrate eq 1 five times over stage 1 using the five values of u. Then, starting from each of the five d l ) vectors obtained over stage 1, we would integrate eq 1 five times over stage 2 using the five values of u. In other words, we would integrate eq 1 25 times over stage 2. It can be seen that the total number of such integrations for the 10 stages is given by Ci1_015i= 12 207 030. We see that even for a modest grid configuration, the exhaustive search requires too many integrations t o be carried out in a reasonable time. This is because the number of x vectors from which the integrations need to be started over stage i is given by 5i.1,and this number increases sharply as i increases. Since n4(t) increases monotonically, the optimal trajectory in a given control grid may be found by carrying out integrations over stage i starting from a relatively small number of x(i - 1) vectors that contain one of the lowest values of 24(i - 1). To explain this, let us assume that we want to carry integrations over any stage from a maximum of 20 x vectors. As in the exhaustive search, start from do)and integrate eq 1 five times over stage 1 using the five values of u. Store the resulting five x(1) vectors and the corresponding u value. Then starting from each of the five dl)vectors, integrate eq 1 five times over stage 2 using the five values of u. At any time, store a maximum of 20 x(2) vectors that contain the 20 lowest values of xq(2). For each of the 20 vectors, store the correspondingu value and the startingx(1) vector. Note that out of the 25 x(2) vectors generated over stage 2, five vectors that contain the five highest values of xd2) are discarded. Now starting from each of the 20 stored vectors, integrate eq 1 five times over stage 3 using the five values of u. Again, store a maximum of 20 x(3) vectors that contain the 20 lowest values of ~ ( 3 ) For . each of these 20 vectors, store the corresponding u value and the starting x(2) vector. Note that out of the 100 x(3) vectors generated over stage 3, 80 vectors that contain the 80 highest values of xq(3) are discarded. Continue this procedure over stages 4 through 10. Then pick up the control trajectory that results in the 0 )moving backward through minimum value of ~ ~ ( 1by the stored results. In the above example, five integrations would be carried over the first stage, 25 integra-

Figure 2. Trajectories of @(dt)) over a three-stage system where two values of control are considered at each stage. The dark line shows the optimal trajectory.

tions over the second stage, and 100 (20 x 5 = 100) integrations over each of the remaining stages. Therefore, the total number of integrations would be 830 which is a greatly reduced number from 12 207 030. For any stage i, the decision about storing or discarding a currently calculated vector can be facilitated by keeping the stored vectors in an ascending order of x4(i)’s. We find that the above semiexhaustive search is enough in many cases t o find the optimal trajectory. Since 4(x(t)) increases monotonically, the minimum value of 6(dtf))for a control grid generally results from a group of x vectors that contain one of the lowest values of $(dt)). For an illustration of this point, see Figure 2 where the plots of 4(x(t)) for a three-stage system with only two values of u are shown. Here it is assumed that t f = 3. In this case the minimum value of $(x(tf))can be found by carrying out integrations over stage 3 from only two of the four x vectors that contain the two lowest values of 4(x(2)). The integration paths represented by dotted lines from the other two x vectors over stage 3 need not be carried. Moreover, it is not necessary t o determine the optimal trajectory for each and every grid configuration. Even if we find a suboptimal trajectory for a particular grid, it is possible t o recapture the optimal trajectory through the search at subsequent iterations and passes. Other features of the method, and the strategies used t o reduce the computational effort, are given below. (a) The starting control trajectory is taken at the midpoint of urnaxand umin unless a better guess is available. The starting span is usually (urnax - umin)/2 unless the search is started from a suboptimal control trajectory where a reduced span may be used. The span at subsequent iterations is reduced by multiplying the previous span by a parameter y . Initially, an approximate optimal control trajectory is found by dividing the time interval (0, tf) into a few time stages (3-5) and carrying out a couple of iterations. This approximate control trajectory is then used as the starting trajectory in the subsequent iterations where the required number of time stages is considered. (b) During the first iteration, the optimal trajectory in the starting grid is found by using the semiexhaustive search procedure described above. The search procedure for second and subsequent iterations is slightly different to make use of the optimal trajectory found. After an integration is carried over a stage using one of the five values of u, the integration is continued over the remaining stages by using the optimal control trajectory obtained at the previous iteration. The value

3880 Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995

of x4(10) reached is recorded. As before, a predetermined number of x ( i ) vectors for stage i are stored. However, these vectors now correspond to the lowest values of X&O) and not x&). (c) The number of x vectors from which integrations are carried is increased somewhat as the stage number increases. This is to partially compensate for the exponential increase in number of x vectors with stage number in the exhaustive search. In the above explanation, this number was assumed to be 20 for every stage to keep the explanation simpler. (d)The number of grid points (values of u) considered on either side of the optimal trajectory is decreased as iterations proceed. In the first pass, a larger integration step size (dtinitial)is first used depending upon what the system would allow without becoming unstable. This is then decreased linearly with iteration number to a final value (dtfind)that would provide sufficient accuracy at the last iteration. In subsequent passes, the integration step size is kept a t its low value dtfinal. (e) The constraints on control variables are satisfied by keeping the control grids within the allowable bounds. In this method, the state constraints can be handled simply by checking the values of the state variables at each integration step and terminating those integration paths that violate the state constraints. (0 When there is more than one control variable, integrations are done such that only one of the control variables is varied at a time. For example, if there are two control variables and we consider five values for each of the variables, then integrations would be done for 25 sets (Ei2 = 25) of control variables. The maximization problems are handled by minimizing the negative of the performance index. 4. Testing

A computer program was developed t o test the performance of the above method. Since IDP found better optimums for several problems, the performance of the method is presented on example problems for which the IDP results are available. The standard fourth-order Runge-Kutta method was used for integrations. The Runge-Kutta method was also used in IDP for most of these example problems. The performance indices were calculated using the same integration step sizes that were used in IDP. For all of the example problems, the performance indices were checked by running the integrations a t half of the final step size (dtfinal)and were found to be the same. The integration step sizes and the values of other parameters used in the program are given in Appendix A. The computations were done on an IBM-type 486DX 33 MHz personal computer, and the computation times for the search are reported. The computation times for the IDP method are also reported wherever available. Since the IDP computations were generally done on a super computer (Cray X-MP/24),approximate equivalent times on a 486/33 MHz personal computer are also indicated by using a conversion factor of 7.6. This conversion factor was calculated from the computation times given by Bojkov and Luus (1992)for the Cray X-MP/24 and a 486/ 33 MHz personal computer. To conserve space, the optimal control trajectories are presented only for the cases where these were different from those obtained using IDP and resulted in better values of performance index. Example Problem 1. Let us consider the problem studied by Ye0 (1980) and Luus (1990)where the system is described by

dx,ldt = x 2 dxdat = -x3u

+ i6t - 8

dx,ldt = u

dx41dt =

xs + x i

(5)

(6)

(7)

+ 0.0005(x2 + 16t - 8 - O . ~ X , U ~(8)) ~

x(0) = [O -1

-&

0IT and t f = 1

The control is bounded by -4 Iu I10 We wish to find a piecewise constant control trajectory over 10 time stages that minimizes x&). For this problem, Ye0 (1980) reported a minimum value of 0.140 52 using quasilinearization. Using IDP, minimum values of 0.120 11(Luus, 1989)and 0.120 12 (Luus, 1990)were found. The semiexhaustive search was started from the middle of the control range, that is, u(i)= 3 where i = 1,2, ..., 10. This initial trajectory was also used by Luus (1990). The convergence of the search to the minimum value of 0.120 11is shown in Figure 3. The minimum was reached in 24 iterations, and the computation time for 25 iterations on the 486/33 PC was 76 s. The lowest time reported by Luus (1989) for 20 iterations of IDP was 11 s on the Cray computer (equivalent to approximately 84 s on a 486/33 PC) where the minimum was reached in 16 iterations. Example Problem 2. This example considers a reaction system consisting of a series of CSTRs studied by Bojkov and Luus (1992). The following difference equations represent the concentrations of the three key components at stage 12. (9)

x,(k) = x,(k

-

+

1) 0.01u,(k)3c2(K)

(11)

where

d o ) = [l 0 O F The control ul(k)is related to the temperature T a t stage

K through the relation u,(b) = 104e -3000/Rk)

(12)

and ug(k) is the residence time in stage K . The constraints are as follows

T(k)I 3 9 4

0

u,(k) I 2 1 0 0

The problem is to find ul(k)and u,&) for K = 1,2, ..., 10 such that ~ ~ ( 1is0maximized. ) The semiexhaustive search was started from the initial control trajectory used by Bojkov and Luus (1992)which is as follows: u1(12) = 0.3 fork = 1,2, ...,8; ul(9) = ul(10) = 4.0; u&) = 0.45 for K = 1, 2, ..., 8; u2(9) = ~ ~ ( 1= 0 )1900. The convergence of the search to the maximum is shown in Figure 4. The maximum value of 0.6314 was reached

Ind. Eng. Chem. Res., Vol. 34, No. 11,1995 3881 0.130

0.128

U C

,

0.438

9 1

a.437

1r

- 0.436 -

U

0.126

C

0

2

.to 0.435 --

a"

h

0.434 -b

0.4331

0.120 10 15 Iteration Number

5

0

20

25

Figure 3. Convergence to the minimum for example problem 1.

0

" ' I

"

5

'

"

:

'

'

'

.

"

10 15 20 Iteration Number

~

' "

25

30

Figure 5. Convergence to the maximum for example problem 3.

We wish to find a piecewise constant temperature

O"

profile over 15 time stages such that x&) is maximized. i For this problem, Luus (1989) reported a maximum

1

0.3

" " "

0

'

' "

10

' " ' I '

" '

' " : '

20

30

" " "

1 40

Iteration Number

Figure 4. Convergence to the maximum for example problem 2.

in 28 iterations, and the computation time for 40 iterations on the 486/33 PC was 41 s. Bojkov and Luus (1992) found the same maximum value and reported a time of 118 s for 30 iterations on a 486/33 PC where the maximum was reached in 16 iterations. Example Problem 3. This example considers the tubular reactor problem that has been studied by Rosenbrock and Storey (1966) and Luus (1989). The mass balance is described by the following differential equations.

value of 0.437 65 for a piecewise constant control trajectory over 15 time stages. Starting from the initial control trajectory used by Luus which was 723.16 throughout the time interval, the semiexhaustive search found a maximum value of 0.437655 which is essentially the same as that obtained by Luus. The convergence of the semiexhaustive search to the maximum is shown in Figure 5. The value of 0.437 65 was reached in five iterations, and the computation time for 30 iterations on the 486133 PC was 77 s. The computation time for 20 iterations of IDP on the Cray computer was 48 s (equivalent to approximately 365 s on a 4861 33 PC) where the maximum was reached in only four iterations. Example Problem 4. This example considers a fedbatch reactor that has been studied by Park and Ramirez (1988) and Luus (1991a). The system is described by the following differential equations.

dxildt = gi(x, - xi) - ( ~ / ~ 5 ) 3 ~ 1 dxddt = g g 3 - (u/x5)xZ dx3/dt = gg3- (u/x5)x3 dx4/dt = -7.3gg3

+ (u/x5)(20- x4)

dx,/dt = u where

+

+

g, = 21.87x4/((x4 o.4)(x4 62.5))

g , = ~,e-&~/(O.l -I-x4)

+

where

g , = 4.75g3/(0.12 g 3 ) x(0) = [O 0 1 5

{ ;[:

k,=k,,exp - - - - -6;8]}

i = 1 , 2 , ..., 5

kio and Ei are given by Luus (1989). The final time (tf) = 1 s, and the constraint on the temperature is 623.16 IT 5 823.16.

1IT

and the constraint is Orur10 We wish to find a piecewise constant control trajectory over 20 time stages such that the total amount of protein

3882 Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995 20835 7

I!

20815

I

354 0

'

'

,

20810

4

5

10

15

20

Iteration Number Figure 6. Convergence to the maximum for example problem 4.

(xl(tf)Xg(tf))produced at tf = 15.5 h is maximized. For this problem, Luus (1991a) reported a maximum value of 37.49 for a piecewise constant control trajectory over 20 time stages. Starting from the initial control trajectory used by Luus which was 0.5 throughout the time interval, the semiexhaustive search found a maximum value of 37.62. The convergence of the search to this maximum is shown in Figure 6. Because of the small integration step size needed, the computation time for 20 iterations on the 486/33 PC was 20 min. This problem was also solved by dividing the time interval into 80 stages. For 80 stages, the IDP resulted in a maximum value of 37.88 (Luus, 1991a). Starting from the control trajectory found over 20 stages and by making several passes, the semiexhaustive search found a different optimal solution with a maximum value of 38.08. The optimal control trajectories are given in Appendix B (Tables 2 and 3). Example Problem 5. This example considers a fedbatch fermenter with a state constraint. This problem has been studied by Chen and Hwang (1990) and Luus (1991b). The system is described by the following differential equations.

&,ldt =glxl - UX~/X,

1 0

5

10

15

20

Iteration Number

Figure 7. Convergence to the maximum for example problem 5.

value of 20 073. For a piecewise constant control trajectory over 20 time stages, Luus (1991b) found a maximum value of 20 814.8. Starting from the optimal trajectory given by Luus (1991b1, the semiexhaustive search found a different optimal solution with a maximum value of 20 830. The convergence to this maximum is shown in Figure 7, and the optimal control trajectory is given in Appendix B (Table 2). The computation time for 20 iterations on the 486/33 PC was 1 h. Example Problem 6. This example considers another fed-batch fermenter for biosynthesis of penicillin where there are constraints on the state variables. This problem has been studied by Lim et al. (19861, Cuthrell and Biegler (1989), and Luus (1993). The system is described by the following differential equations. &,ldt = hlxl - x,u/(~OOX,) h2/dt = h g l - 0 . 0 1 ~ x~, u / ( ~ O O X , )

(25) (26)

dx31dt = - h,~,/0.47- hg1/1.2 - 0.029 x,x,/ (0.0001 x,) (1 - x , / ~ O O ) U / X ~(27)

+ +

(21)

&,/dt = uI500

(28)

where

+

h, = 0.1l~,/(0.006~, x,) h, = 0.0055x3/(0.0001 where

+ x,(l + lox,))

x(0) = [1.5 0 0 7IT

and the constraints are

g, = x&l

+ X3/71.5)(O.44+ x,))

x(0) = [ l 150 0 10IT

and the constraints are

0 5 u 5 12

x,(t) 5 200

We wish to find a piecewise constant feed rate (u)over 20 time stages such that the total amount of ethanol (~3(tf)X&))produced at tf = 62 h is maximized. For this problem, Chen and Hwang (1990) reported a maximum

Oru550

0 5 ~ ~ 5 4 0

0 5 ~ ~ 5 2 05 5 ~ ~ 1 1 0 We wish to find a piecewise constant feed rate (u)over 20 time stages such that the total amount of penicillin (xz(tf)x&)) produced at t f = 132 h is maximized. For this problem, Cuthrell and Biegler (1989) reported a maximum value of 87.69. For a piecewise constant trajectory over 20 time stages, Luus (1993) found a maximum value of 87.948. Starting from a control trajectory given by Luus (19941, the semiexhaustive

Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995 3883 88.02

88.00

87.92

Table 2. Optimal Control Trajectories over 20 Stages for Example Problems 4-6

1 I

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

I

0

3

6 9 Iteration Number

12

15

Figure 8. Convergence to the maximum for example problem 6. Table 1. Values of Parameters example problem no. 1

0.1 0.01 9 5 3 3 3 3 1 1 1 1 3 0.75

2

9 7 5 5 5 3 65 15 30 30 0 0.5

3

4

5

6

0.02 0.01 7 7 5 5 3 3 10 10 10 10 0 0.75

0.141 0.097 9 7 7 7 5 3 105 50 30 15 0 0.5

0.05 0.05 3 3 3 3 3 3 30 30 30 30 2 0.5

0.025 0.025 3 3 3 3 3 3 2 3 13 1 0 0.5

search found a different optimal solution with a maximum value of 88.00. The convergence to this maximum is shown in Figure 8, and the optimal control trajectory is given in Appendix B (Table 2). In this case, the main increase in the performance index did not occur until the grid was quite narrow around the starting optimal trajectory. The computation time for 15 iterations on the 486133 PC was 47 min. 5. Comments As reported in the literature, IDP found better optimums for several problems than were found using other methods. In some cases, IDP found a better optimum than was found earlier using IDP due to the difference in the way it was used. We found that the IDP results could sometimes be improved. Although the improvements may not be large, the different results obtained using different methods show that the current methods do not always find global optimums. Although the semiexhaustive search does not guarantee the determination of the global optimum either, the different approach it employs is useful for cross-checkingthe results and for seeing if the results obtained using other methods can be further improved. The chances of obtaining a global optimum can be improved by carrying out integrations from a larger number of x vectors. However, the computational effort would also increase. To handle the state constraints, some methods only ensure that there are no violations at the end of the time stages. In our tests, we found that this approach can result in control trajectories that violate state con-

example problem no.

0.000 000 0.000 000 0.000 000 0.000 000 0.716 862 1.105 983 1.296 998 1.486 797 1.837 077 2.158 562 2.538 920 3.112 755 3.453 638 4.622 585 4.124 413 7.044 401 6.801 393 8.997 809 11.992 029 0.000 000

0.161 560 0.206 116 0.260 539 0.328 152 0.411 936 0.508 358 0.682 597 0.829 926 1.052 216 1.310 852 1.770 142 2.454 956 1.140 971 0.000 000 0.862 345 0.875 773 0.880 656 0.896 525 0.941 013 1.388 045

3.647 681 6.117 642 26.526 321 46.687 023 8.420 924 8.478 114 8.706 964 8.802 229 8.897 632 8.916 842 8.897 593 9.012 148 9.012 195 9.031 152 9.183 883 9.241 056 9.412 714 9.450 746 9.489 017 9.508 100

Table 3. Optimal Control Trajectory over 80 Stages for Example Problem 4

i 1

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

u(20

u(i) 0.156 067 0.148 743 0.181 702 0.165 527 0.200 419 0.201 335 0.206 217 0.223 307 0.253 220 0.228 806 0.288 620 0.264 893 0.335 726 0.295 443 0.357 699 0.336 947 0.405 858 0.374 120 0.442 479 0.463 231

+i)

0.449 677 0.553 436 0.494 843 0.607 147 0.596 759 0.657 794 0.648 028 0.762 774 0.765 368 0.792 224 0.852 038 0.918 261 0.977 664 0.971 865 1.145 205 1.124 453 1.235 664 1.272 896 1.414 497 1.356 209

u(40

+ i)

1.497 925 1.661 804 1.636 932 1.630 524 2.178 302 1.895 709 2.683 062 3.127 093 1.569 387 2.735 463 3.056 508 0.453 796 0.000 000 0.000 000 0.000 000 0.000 000 0.000 000 0.842 941 0.923 203 0.908 249

u(60

+ i)

0.911 911 0.910 690 0.913 437 0.913 742 0.914 047 0.916 489 0.918 320 0.921 982 0.924 423 0.929 306 0.936 020 0.944 565 0.956 899 0.972 768 0.995 046 1.029 836 1.091 176 1.193 410 1.457 082 10.000 000

straints within a time stage. In the semiexhaustive search, since the violations were checked after each integration step, the optimal trajectories satisfied state constraints within and at the end of time stages. The time requirements of the IDP and the semiexhaustive search in locating the optimums were more or less the same.

6. Conclusions The semiexhaustive search provides a viable method that can be used by itself or with other methods for solving nonlinear optimal control problems.

Acknowledgment The financial support provided by the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged.

Appendix A The number of grid points (values of u ) considered at any stage were decreased with iterations. The value of

3884 Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995

the span reduction factor ( y ) was usually 0.5. For some problems, it was initially kept a t 0.75 and then decreased to 0.5. The changes in these variables with iteration number are shown in the following code: IF Iteration = 1 THEN No-Grid-Points = nG1: y = Yinitial ELSEIF Iteration = 2 THEN No-Grid-Points = nGz: y = Yinitial ELSEIF Iteration = 3 THEN No-Grid-Points = nG3: 7 = yinitial ELSEIF Iteration < (1/3 * MaximumIterations) THEN No-GrihPoints = nG4: y = 0.5 ELSEIF Iteration < (2/3 * MaximumIterations) THEN No-Grid-Points = nGs: y = 0.5 ELSE No-GrihPoints = nGs: y = 0.5 END IF The maximum number of state vectors (n-Max) that were stored at any stage for carrying out integrations over the next stage was varied with iteration number and stage number as follows: IF Iteration = 1 THEN n =m l ELSEIF Iteration < (1/3 * MaximumIterations) THEN n=mz ELSEIF Iteration < (2/3 * MaximumIterations) THEN n=m3

ELSE n=m4 END IF n X a x (at stage i) = n (i - 1)dn The initial integration step size (dtinitial), the final integration step size (dtfinal), and the values of other parameters, nGi (i = 1, 2, ..., 61, mi (i = 1-41, dn, and yinitial are given in Table 1.

+

Appendix B

Literature Cited Bojkov, B.; Luus, R. Use of Random Admissible Values for Control in Iterative Dynamic Programming. Ind. Eng. Chem. Res. 1992, 31, 1308-1314. Chen, C. T.; Hwang, C. Optimal Control Computation for Differential-Algebraic Process Systems with General Constraints. Chem. Eng. Commun. 1990,97,9-26. Cuthrell, J. E.; Biegler, L. T. Simultaneous Optimization and Solution Methods for Batch Reactor Control Profiles. Comput. Chem. Eng. 1989,13,49-62. Lim, H. C.; Tayeb, Y. J.; Modak, J. M.; Bonte, P. Computational Algorithms for Optimal Feed Rates for a Class of Fed-Batch Fermentation: Numerical Results for Penicillin and Cell Mass Production. Biotechnol. Bioeng. 1986,28, 1408-1420. Luus, R. Optimal Control by Dynamic Programming using Accessible Grid Points and Region Contraction. Hung. J . Ind. Chem. 1989,17, 523-543. Luus, R. Optimal Control by Dynamic Programming Using Systematic Reduction in Grid Size. Int. J . Control 1990, 51, 995-1013. Luus, R. Effect of the Choice of Final Time in Optimal Control of Nonlinear Systems. Can. J . Chem. Eng. 1991a, 69, 144-151. Luus, R. Application of Iterative Dynamic Programming to State Constrained Optimal Control Problems. Hung. J . Ind. Chem. 1991b, 19, 245-254. Luus, R. Piecewise Linear Continuous Optimal Control by Iterative Dynamic Programming. Ind. Eng. Chem. Res. 1993, 32, 859-865. Luus, R. University of Toronto, personal communication, 1994. Park, S.; Ramirez, W. F., Optimal Production of Secreted Protein in Fed-Batch Reactors. M C h E J . 1988, 34, 1550-1558. Rosenbrock, H. H.; Storey, C. Computational Techniques for Chemical Engineers; Pergamon Press: New York, 1966; pp 241-277. Yeo, B. P. A Modified Quasilinearization Algorithm for the Computation of Optimal Singular Control. Int. J. Control 1980, 32. 723-730.

Received for review February 27, 1995 Accepted July 11, 1995@ IE9501438

Abstract published in Advance A C S Abstracts, September 15, 1995. @

Tables 2 and 3 comprise Appendix B.