Optimal Control and Optimal Time Location Problems of Differential

The differential equations are used to describe the slower dynamic behaviors of the processes, ..... ]T, where The crossover factor CR ∈ [0,1] is se...
0 downloads 0 Views 327KB Size
Download PDF
5348

Ind. Eng. Chem. Res. 1997, 36, 5348-5357

Optimal Control and Optimal Time Location Problems of Differential-Algebraic Systems by Differential Evolution Feng-Sheng Wang* Department of Chemical Engineering, National Chung Cheng University, Chiayi 621, Taiwan

Ji-Pyng Chiou Department of Electrical Engineering, National Chung Cheng University, Chiayi 621, Taiwan

An efficient method is introduced for solving optimal control and optimal parameter selection problems of nonlinear differential-algebraic systems involving general constraints. These infinite dimensional problems are converted into the finite dimensional optimization problems by using the control parametrization technique. Moreover, the transformed time locations for the control variables are included in such problems. As a result, optimal control problems become acausal optimal parameter selection problems. A modified version of differential evolution is introduced to solve these problems. The integration of the penalty functions is used to ensure the solution onto the feasible domain of the problems. Two industrial systems, a chemical reactor and a robotics arm, are used to illustrate the robustness and effectiveness of the method. The control effects can be improved by selecting the proper time locations as observed from the simulated results. Introduction Evolutionary algorithms (EAs) have received considerable and increasing interest over the past decade. Such algorithms have been used to solve many practical problems (Androulakis and Venkatasubramanian, 1991; Seywald and Kumar, 1995; Ma and Wu, 1995; Moros et al., 1996). Evolutionary algorithms are a class of stochastic search and optimization methods based on the mechanics of natural biological evolution that includes genetic algorithms (Holland, 1975), evolutionary programming (Fogel, 1994), evolution strategies (Schwefel, 1995), and genetic programming (Koza, 1994). EAs operate on a population of potential solutions by applying the principle of survival of the fittest to produce successively better approximations to an optimal solution. These algorithms are robust and global when compared to traditional optimization methods, such as gradient-based and enumerative strategies. While many excellent EA packages are available, users generally demand that a practical optimization technique should fulfill three requirements. First, the method should find the true global minimum, regardless of the initial guess. Second, convergence should be fast. Third, the program should have a minimal number of free factors so that it is easy to use. Storn and Price (1996) have developed a differential evolution method that not only is simple but also performs well on a wide variety of test problems. Differential evolution (DE) turned out to be the best genetic type of algorithm for solving the real-valued test function suite of the first International Contest on Evolutionary Computation that was held in Nagoya, 1996. Although DE has shown promising results (Storn and Price, 1996), little is known about applying DE to solve optimal control problems in industrial systems. Many practical computing techniques for solving optimal control problems have drawn great attention over the last two decades. Most of these methods successfully solve simple constrained problems, but the presence of, for example, algebraic equality or inequality * Author to whom correspondence should be addressed. E-mail: [email protected]. S0888-5885(97)00248-0 CCC: $14.00

constraints often results in both analytical and computational difficulties. Moreover, many dynamic optimization problems in industrial processes are described by a set of coupled differential and algebraic equations called differential-algebraic equations (DAEs). The differential equations are used to describe the slower dynamic behaviors of the processes, while the algebraic equations are used to described the faster phases or the constraints of the problems. The modeling of chemical processes (Butala et al., 1992; Soroush and Kravaris, 1993) and the dynamic analyses of nonlinear circuits (Newcomb, 1981) are typical DAEs. Dynamic optimization of such systems has recently been of interest in the literature. Pontryagin minimum principle method has been used to transform the problems into nonlinear twopoint boundary value problems (TPBVP) (Bryson and Ho, 1975). Even though the order of the system is small, the analytical solution for the TPBVP is difficult to obtain. The well-documented gradient method in function space is iterative in nature and is extended for solving the problems (Jones and Finch, 1984). The method requires the integration of two sets of differential equations, the system state equations and their co-state equations, to provide the directional gradient for control vector iteration. Methods of state and control parametrization (Logsdon and Biegler, 1989) have been proposed to convert dynamic optimization problems into nonlinear programming problems so that a standard mathematical programming technique could be employed to solve these approximate problems. Control parametrization techniques (Chew and Goh, 1989) have been applied to convert optimal control problems into finite dimensional optimization problems. A unified computational technique (Teo et al., 1991; Chen and Hwang, 1990) has been used to provide the descent gradients in these finite dimensional optimization problems. The time partitions for the control variables have to be selected first in solving these finite dimensional optimal control problems. Little discussion has appeared concerning the selection of time partitions in the optimal control problems. In this study, the differential evolution algorithm is applied for determining the optimal control solution for © 1997 American Chemical Society

Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 5349

problems described by differential-algebraic equations. A control parametrization method is used to convert such infinite dimensional systems to finite ones. Time partitions in this discretization must be selected by the user. Equal time intervals are typically used. In this study, a method for selecting the time partitions within the optimization problem is introduced. With this scheme, the user need not provide any initial guess for the control policy and time partitions. Instead, initial guesses that lie within the domain of convergence for a gradient-based method is generated by the DE algorithm. This evolved solution is then used as the initial starting point for sequential quadratic programming in order to find a more accurate solution.

problem is to find the control policy u(t) and the parameter vector θ so that the performance index J is minimized without violating constraints (4)-(7) of the problem. Approximation of the Problem Control Parametrization. In an optimal control problem, a control u(t) as a function of time is to be chosen. This is an infinite dimensional problem. To solve this problem efficiently, the control is represented by a finite set of control parameters u(j). The continuous time is first partitioned as

t0 ) 0 < t1 < t2 (zmax)k, k ) 1, ..., ni; i ) 1, ..., Np

5352 Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997

Step 6. Repeat step 3 until the maximum iteration is reached or the best function value has achieved the desired level. As discussed in step 3, many, but not all, DAEs can be solved conveniently using numerical codes for ordinary differential equations (ODEs) (Audry-Sanchez, 1988). Historically, a set of DAEs was solved as ODEs often destroying the natural structure of the system (Pantelides et al., 1988; Unger, 1995). This was done primarily to cast such problems in the form required by the particular integration method or software package being used. Today, it is becoming more common to deal with the original DAEs. One important reason for this (Brenan, 1989) is that the variables in the original DAEs typically have some physical significance, whereas those that result after manipulation into ODEs may not. Thus, by solving the original DAEs directly, the effect of parameter changes or model changes is more easily studied. In this study, a modified collocation method (Wang and Chiou, 1995; Wang and Chiou, 1997) is used to solve DAE in (28).

dx5 ) {kty + kp[β1 + β2 + cmm0(1 - x3)] + kid1x1 + dt kid2x2}λ1 - kp(1 + cm)y(1 - x3)x5 dx6 ) {kty + kp[β1 + β2 + cmm0(1 - x3)] + kid1x1 + dt kid2x2}λ2 + ktλ21 - kp(1 + cm)y(1 - x3)x6 (36) where the expressions for primary radical concentration are

β1 )

β2 )

(2kd1x1 + kid1x1y)K1 m0(1 - x3) + K1kpy (2kd2x2 + kid2x2y)K2 m0(1 - x3) + K2kpy

and the volume contraction factor is given as

 ) Fm(1/Fp - 1/Fm) Simulation Examples The computational schemes were tested on a wide variety of optimal control problems. Two examples are presented here. The original version of the DE in Matlab code was modified to implement the proposed algorithm and to test these problems. Example 1. An optimal temperature control of batch styrene polymerization initiated with mixed initiator system is considered in this example. Reducing the batch reaction time and obtaining the desired polymer properties are simultaneously required in this batch free radical polymerization process. Although one can use both high-temperature and high-initiator concentration to achieve a near-complete monomer conversion in a short reaction time, it is often difficult to achieve both high monomer conversion and desired high polymer molecular weight simultaneously by doing so. It has been well known that merely increasing the amount of initiator for rapid polymerization results in lower polymer molecular weight, gelation, and discoloration of polymer product. The polymer molecular weight decreases with an increase in reaction temperature. Therefore, it is desirable to select the proper concentration of initiator and polymerization temperature profile to achieve both high monomer conversion and desired high polymer molecular weight, simultaneously. The free radical polymerization of styrene in a batch reactor catalyzed by mixed initiators modeled by a set of implicit formulas of differential-algebraic equation is employed in this example. The modeling equations discussed in Butala et al. (1992) are as follows:

dx1 ) -(kd1 + kid1y)x1 - kp(1 + cm)x1(1 - x3)y dt dx2 ) -(kd2 + kid2y)x2 - kp(1 + cm)x2(1 - x3)y dt dx3 ) kp(1 + cm)[1 + (1 - x3)](1 - x3)y dt

{

kt dx4 ) y + kp[β1 + β2 + cmm0(1 - x3)] + kid1x1 + dt 2

}

kid2x2 y - kp(1 + cm)y(1 - x3)x4

Here x1 and x2 denote concentration of the first and second initiator, respectively; x3 is conversion of monomer; and x4, x5, and x6 are the corresponding zeroth, first, and second moment of dead polymer distributions. The physical data in the above equations discussed in Butala et al. (1992) are shown in Table 2. The modeling equations for the zeroth, first, and second moments λ0, λ1, and λ2 of live polymer radicals are based on a quasi-steady-state approximation as in the forms

λ0 ) y λ1 ) λ2 )

a(y + γ) 1-a

a(y + γ + 2λ1) 1-a

(37)

where

a ) m0(1 - x3)/[m0(1 - x3) + β1 + β2 + kid1x1/kp + kid2x2/kp + m0cm(1 - x3) + (ktc + ktd)y/kp] γ)

2kd1x1 kp[m0(1 - x3) + K1kpy]

+

2kd2x2 kp[m0(1 - x3) + K2kpy] cm )

+

2kdmm20(1 - x3)2 + cmy kp

kfm kp

Here the total concentration of live polymer y is governed by the following nonlinear algebraic equation:

2kd1x1 + 2kd2x2 + 2kdmm30(1 - x3)3 - 2β1kpy 2β2kpy - (ktc + ktd)y2 ) 0 (38) Since the parameters β1 and β2 in (38) are also in terms of y, the concentration of live polymer cannot be explicitly obtained from this constrained equation. The molecular weight properties, e.g., number average chain length and molecular weight distribution, are

Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 5353 Table 2. Definition of Physical Parameters for the Batch Polymerization of Styrene

0 mol/L e x1(0) e 0.005 mol/L

parametersa

0 mol/L e x2(0) e 0.005 mol/L

kd1 ) 1.44 × 1013 exp(-14700/T) kd2 ) 8.315 × 1013 exp(-16162/T) kdm ) ki ) 2.19 × 105 exp(-13810/T) kp ) kp0 ) 1.06 × 107 exp(-3569/T) Cm ) kfm/kp ) kfm0/kp0 + B1x3 kfm0 ) 2.31 × 106 exp(-6377/T)

B1 ) -1.013 × 10-3 log10

(473.12-T 202.5 )

kid1 ) 3.116 × 109 exp(-6359/T) kid2 ) 1.62 × 1013 exp(-9518/T) kt ) ktc ) kt0gt kt0 ) 1.25 × 106 exp(-847/T) gt ) exp[-2(A1x3 + A2x23 + A3x33)] A1 ) 2.57 - 5.05 × 10-3T A2 ) 9.56 - 1.76 × 10-2T A3 ) -3.03 + 7.85 × 10-3T

K1 )

ktp1 ) 3.44 × 10-14 exp(14022/T) ki1kp

K2 )

ktp2 ) K1 ki2kp

Fm ) 924.0 - 0.918(T - 273.1) Fp ) 1084.8 - 0.605(T - 273.1) a

Units: g, mol, L, s, cal, K.

computed using the method of moments. The number average chain length or the degree of polymerization XN is therefore calculated as

XN )

x5 + λ1 x4 + y

(39)

As mentioned earlier, we are interested in the minimization of batch time while the desired polymer properties are achieved. These properties are expressed in the forms

(

)

XNd - XN(tf) XNd

2

e 0.05

(40)

x3(tf) e 0.99

(41)

x1(tf) + x2(tf) e 10 ppm

(42)

where the desired value of number average chain length is set to be XNd ) 1000 in this example. The first constraint requires that the number average chain length at the final time XN(tf) has to be within 5% of the desired level XNd. The final conversion x3(tf) has to be greater than or equal to 99% as expressed in (41). The last constraint denotes that the initiator residue concentration at the final time has to be less than 10 ppm. The performance index is defined as

min

T(t),x1(0),x2(0)

tf

(43)

the goal of this objective is to find the optimal reaction temperature T(t) and the initial concentration of the mixed initiators, x1(0) and x2(0), to achieve the desired polymer properties (40)-(43). These decision variables are determined within the bounded region as

343 K e T(t) e 413 K

By choosing 10 time partitions for the temperature profile, 23 parameters (i.e., 10 constant temperature actions, 10 time interval parameters, two initial concentrations of initiators, and final time) have to be determined in this parameter selection problem. In this work, the scaling factor F ) 0.8, crossover ratio CR ) 0.8, and the second strategy of mutation operator were used in the DE algorithm. A modified collocation method (Wang and Chiou, 1995; Wang and Chiou, 1997) was used to solve the differential-algebraic equations (36)-(38). The dynamic-state variables and steadystate variables were then provided for evaluating the fitness function, which was the performance index (43) of the problem including the penalty functions. By choosing population size Np ) 20, we were able to obtain convergence to 233.24 after 100 iterations. Histories of the best and worst fitness values are shown in Figure 1 where it can be seen that the best fitness values converge very quickly. To obtain more accurate results, sequential quadratic programming was then employed to solve this optimal control problem again. The best results were used as the initial starting point for a modified SQP solver which was an extended version of the modified reduced gradient method for solving optimal control problems (Wang and Chiou, 1995; Wang and Chiou, 1997). The optimal results are represented in Figures 2-5. From these results, three inequality constraints, (40), (41), and (42), were satisfied. The initiator residue at final time tf ) 164.438 min was less than 10 ppm, which achieved the desired level. The first initiator was almost exhausted at t > 85 min, as observed in Figure 2. The monomer conversion at final time achieved 99%, as observed in Figure 3. The number average chain length achieved 776.39, as shown in Figure 4. The optimal stepwise temperature profile was continuously increasing from 389.1 to 413.0 K, as observed in Figure 5. Although use of the modified SQP yielded rapid convergence to the refined optimal solution, it required that very good initial starting points be used in order to ensure that the global optimum was obtained. The modified SQP alone was used to solve this example. However, the optimal solution was not easy to obtain if the time interval parameters νi were also considered as the decision parameters in this optimization problem. To test the robustness of DE, three different scaling factors, F ) 0.2, 0.4, and 0.6, respectively, were used to run this problem again. The best fitness value related to each scaling factor was 209.96, 190.97, and 237.84, respectively. Each best solution was then used as the initial starting point for the modified reduced gradient method to obtain a more accurate solution. The optimal solutions by using these four initial guesses gave almost the same results. In this example, five time partitions were used to discretize the control variable of temperature for the same problem. The proposed DE algorithm was employed to find a converged solution. This solution was then used as the initial starting points for the modified SQP solver to find the more accurate solution. The optimal results are shown in Figures 2-5 as dashed line for comparison. Although the control actions were different from these discussed in the first result, as shown in Figure 5, nearly identical control effects, as observed in Figures 2-4, could be obtained using this

5354 Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997

Figure 1. History of the best and worst fitness values (b, best fitness value; O, worst fitness value).

Figure 2. Optimal concentration of the first and second initiators (solid lines, 10 time partitions; dashed lines, 5 time partitions).

smaller number of time partitions. Less control actions for the reaction temperature were more convenient to implement this optimal temperature control of the batch polymerization. Example 2. Here, we consider a robotic arm that consists of two identical links that is assumed to be rigid and to move in the horizontal plane (Weinrb and Bryson, 1985; Lee, 1992). The joints are assumed to be frictionless and to have unlimited angular freedom. The equations of motion given by Lee (1992) are expressed as follows:

dx1 4 9 ) sin(x3) cos(x3)x21 + 2x22 + (u1 - u2) dt 4 3 31 9 3 cos(x3)u2 + sin2(x3) 2 36 4

[

(

)

]/[

dx2 9 7 ) sin(x3) x21 + cos(x3)x22 + dt 2 4 7 3 cos(x3)(u1 - u2) - u2 2 3

[

(

]

)

31 9 + sin (x )] ] / [36 4 2

3

Figure 3. Optimal trajectory of monomer conversion (solid lines, 10 time partitions; dashed lines, 5 time partitions).

Figure 4. Optimal trajectory of number average chain length (solid lines, 10 time partitions; dashed lines, 5 time partitions).

dx3 ) x2 - x1 dt dx4 ) x1 dt

(44)

where x1 and x2 are the inertial angular velocities of the shoulder and elbow links, respectively, and x3 and x4 represent angles. The control variables u1 and u2, representing the torque, are bounded by

|uj(t)| e 1, j ) 1, 2 The aim of this optimization problem is to find the control policy of the robotic arm such that the arm starts from the initial state x(0) ) [0.0 0.0 0.5 0.522]T and reaches the final state x(tf) ) [0.0 0.0 0.5 0.0]T in minimum time tf. The performance index, including the penalty terms at the final time, is therefore expressed as

Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 5355

Figure 5. Optimal temperature control policy (solid lines, 10 time partitions; dashed lines, 5 time partitions).

min J ) tf + R1x1(tf)2 + R2x2(tf)2 + R3(x3(tf) - 0.5)2 +

u1,u2,tf

Figure 6. Optimal control action of u1 for 13 time partitions (solid line for the optimal time location; dashed line for locating times at roots of Chebyshev polynomial; dotted line for equal time interval).

R4(x4(tf) - 0.522)2 This optimization problem is a singular/bang-bang control problem. The control policy should relate to switching times. Accordingly, we focus our attention on strategies of time partitions in the control parametrization approach. Three strategies, the optimal time location, equal time interval, and locating times at roots of Chebyshev polynomial (Boor, 1978; Wang and Chiou, 1997), respectively, are applied to discretize the control variables. In this example the modified SQP alone was applied to solve these optimal parameter selection problems by using many different sets of randomly distributed starting points. However, a converged solution could not be obtained by the modified SQP alone. This failure could be overcome by the proposed approach. DE was first used to solve these parameter selection problems to find a converged suboptimal solution. The modified SQP used this suboptimal solution as the starting points to solve the parameter selection problems again to find a more accurate solution. The second strategy of mutation operator of DE with population size Np ) 20, scaling factor F ) 0.8, and crossover factor CR ) 0.8 were used in this work. By choosing 13 time partitions for the control variables, a total of 40 parameters, i.e., 26 constant control actions for u1 and u2, 13 time interval parameters, and final time tf, had to be determined in this problem. However, when equal time interval or locating times at roots of Chebyshev polynomial was used to discretize these control variables, a total of only 27 parameters had to be determined. The computational results of the control policy for these three different time partitioned strategies are represented in Figures 6 and 7 for comparison. The optimal final time was 2.9839 min for the strategy of optimal time locations. However, the optimal final time of tf ) 3.2305 or 3.0092 min was respectively obtained by the strategies of equal time interval or locating times at roots of Chebyshev polynomial. The normalized time partitions are expressed in Table 3. In this work, 20 time partitions were also used to discretize the control variables in this example. The

Figure 7. Optimal control action of u2 for 13 time partitions (solid line for the optimal time location; dashed line for locating times at roots of Chebyshev polynomial; dotted line for equal time interval).

computational results of the control policy for these three different time partitioned strategies are represented in Figures 8 and 9 for comparison. The optimal final time was 2.9832 min for the strategy of optimal time location. However, the optimal final time of tf ) 2.9870 or 2.9893 min was respectively obtained by the strategies of equal time interval or locating times at roots of Chebyshev polynomial. The normalized time partitions are expressed in Table 3. From Figures 6-9, we observed that the switching times in this singular optimal control problem could be determined by the proposed strategy so that the control actions for 13 and 20 time partitions were nearly identical ones at such times. On the other hand, time partitions were fixed by the strategies of equal time interval or locating times at roots of Chebyshev polynomial so that the control

5356 Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 Table 3. Normalized Time Partitions in Example 2a 13 time partitions

20 time partitions

location

OT

ET

CT

OT

ET

CT

t0 t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 t11 t12 t13 t14 t15 t16 t17 t18 t19 t20

0.0 0.10359 0.15954 0.18864 0.28774 0.31592 0.38589 0.41221 0.49985 0.61251 0.72517 0.82837 0.88733 1.0

0.0 0.07692 0.15384 0.23077 0.30769 0.38462 0.46154 0.53846 0.61538 0.69231 0.76231 0.84615 0.92308 1.0

0.0 0.02507 0.07396 0.14441 0.23230 0.33382 0.44366 0.55634 0.66618 0.76770 0.85579 0.92604 0.97493 1.0

0.0 0.07473 0.13175 0.14931 0.16659 0.22255 0.24259 0.31364 0.32481 0.35252 0.42520 0.49993 0.57466 0.64939 0.72412 0.79885 0.81611 0.87870 0.91441 0.92558 1.0

0.0 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.0

0.0 0.01117 0.03326 0.06577 0.10799 0.15896 0.21755 0.28245 0.35221 0.42527 0.50 0.57473 0.64779 0.71755 0.78245 0.84704 0.89201 0.93423 0.96674 0.98883 1.0

a

OT for optimal time location, ET for equal time interval, and CT for locating times at roots of Chebyshev polynomial.

Figure 9. Optimal control action of u2 for 20 time partitions (solid line for the optimal time location; dashed line for locating times at roots of Chebyshev polynomial; dotted line for equal time interval).

used to ensure that solutions fell within the feasible domain of the problems. A modified collocation method was used to solve differential-algebraic equations. The dynamic-state variables and steady-state variables were then provided for evaluating the performance index including penalty functions in DE. Control parameters as well as time interval parameters could be efficiently determined by this modified differential evolution as observed from the simulated examples. Such decision parameters were then provided as a good initial guess to the standard gradient-based algorithm to obtain a more accurate solution. Acknowledgment Financial support from the National Science Council of the ROC under grant NSC85-2214-E194-004 is gratefully acknowledged. Nomenclature Figure 8. Optimal control action of u1 for 20 time partitions (solid line for the optimal time location; dashed line for locating times at roots of Chebyshev polynomial; dotted line for equal time interval).

actions were different when the locating times were not near the switching times. This fact could be observed from Figures 6 and 8. The first control action for these three time partitioned strategies was u1(t) ) 1 during the normalized time t/tf ) 0.5. After the switching time, this control action became -1, as observed in Figure 8. However, for 13 time partitions, this switching time could not be collocated by the strategy of equal time interval or locating time at roots of Chebyshev polynomial so that the control actions were different near this switching time, as shown in Figure 6.

Djk, Djklm ) difference vector f ) a vector function of nonlinear differential equations F ) scaling factor g ) inequality constraints h ) equality constraints J ) objective function Nu ) time stages nq ) number of system parameters nt ) total number of control parameters and system parameters u ) control variables u(j) ) control action used at the stage j X ) a vector of state variables x ) dynamic-state variables y ) steady-state variables z ) a vector of control parameters and system parameters ZG ) a population at the Gth generation

Conclusions

Greek Letters

In this paper, the differential evolution method has been used to solve optimal control problems described by differential-algebraic systems with nonlinear constraints. The integration of the penalty functions was

R1,R2,R3 ) penalty constants ηi ) variable length for the ith time interval of normalized time domain νi ) variable length for the ith time interval

Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 5357 θ ) parameter variables τ ) normalized time variable

Literature Cited Androulakis, I. P.; Venkatasubramanian, V. A genetic algorithmic framework for process design and optimization. Comput. Chem. Eng. 1991, 15 (4), 217-228. Audry-Sanchez, J. On the numerical solution of differential algebraic equations. Can. J. Chem. Eng. 1988, 66, 1031-1035. Boor, C. D. A Practical Guide to Splines; Springer-Verlag: Berlin, 1978. Brenan, K. E.; Campbell, S. L.; Petzold, L. R. Numerical Solution of Initial-Value Problems in Differential-Algebraic Equation; North-Holland: Amsterdam, 1989. Bryson, A. E.; Ho, Y. C. Applied Optimal Control; Hemisphere: New York, 1975. Butala, D. N.; Liang, W. R.; Choi, K. Y. Multiobjective dynamic optimization of batch free radical polymerization process catalyzed by mixed initiator systems. J. Appl. Polym. Sci. 1992, 44, 1759-1778. Chen, C. T.; Hwang, C. Optimal control computation for differential-algebraic process systems with general constraints. Chem. Eng. Commun. 1990, 97, 9-26. Chew, E. P.; Goh, C. J. On minimum time optimal control of batch crystallization of sugar. Chem. Eng. Commun. 1989, 80, 225231. Fogel, D. B. Evolutionary programming: an introduction and some current directions. Stat. Comput. 1994, 4, 113-129. Holland, J. H. Adaptation in Natural and Artificial Systems; University of Michigan Press: Ann Arbor, 1975. Jones, D. I.; Finch, J. W. Comparison of optimization algorithms. Int. J. Control 1984, 40, 747-761. Koza, J. R. Genetic Programming II: Automatic Discovery of Reusable Programs; MIT Press: Cambridge, MA, 1994. Lee, A. Y. Solving Constrained Minimum-Time Robot Problems Using Sequential Gradient Restoration Algorithm. Optim. Control Appl. Methods 1992, 13, 145-154. Logsdon, J. S.; Biegler, L. T. Accurate solution of differentialalgebraic optimization problems. Ind. Eng. Chem. Res. 1989, 28, 1628-1639. Ma, J. T.; Wu, Q. H. Generator parameter identification using evolutionary programming. Electr. Power Energy Syst. 1995, 17 (6), 417-423. Michalewicz, Z. Genetic Algorithms + Data Structures ) Evolution Programs; Springer-Verlag: New York, 1996. Michalewicz, Z.; Schoenauer, M. Evolutionary algorithms for constrained parameter optimization problems. Evol. Comput. J. 1996, 4 (1), 1-32. Moros, R.; Kalies, H.; Rex, H. G.; Schaffarczyk, St. A genetic algorithm for generating initial parameter estimations for

kinetic models of catalytic processes. Comput. Chem. Eng. 1996, 20, 1257-1270. Newcomb, R. W. The semistate description of nonlinear timevariable circuits. IEEE Trans. Circuits Syst. 1981, CAS-28, 6271. Pantelides, C. C.; Gritsis, D.; Morison, K. R.; Sargent, R. W. H. The mathematical modeling of transient systems using differential-algebraic equations. Comput. Chem. Eng. 1988, 12, 449-454. Price, K. V. Differential evolution: a fast and simple numerical optimizer. Proc. of 1996 North American Fuzzy Information Processing Society, Berkeley, CA, 1996, ISBN:0-7803-3225-3. Schwefel, H. P. Evolution and Optimum Seeking; John Wiley: New York, 1995. Seywald, H.; Kumar, R. R. Genetic algorithm approach for optimal control problems with linearly appearing controls. J. Guid., Control Dyn. 1995, 18 (1), 177-182. Soroush, M.; Kravaris, C. Optimal design and operation of batch reactors. 1. Theoretical framework. Ind. Eng. Chem. Res. 1993, 32, 866-881. Storn, R. On the usage of differential evolution for function optimization. NAFIPS, Berkeley, 1996, 519-523. Storn, R.; Price, K. Minimizing the real functions of the ICEC’96 contest by differential evolution. IEEE Conf. on Evolutionary Computation, 1996, Nagoya, 842-844. Teo, K. L.; Goh, C. J.; Wong, K. H. A Unified Computational Approach to Optimal Control Problems; Bath Press: New York, 1991. Unger, J.; Krıˆner A.; Marquardt, W. Structural analysis of differential-algebraic equation systems-theory and applications. Comput. Chem. Eng. 1995, 19, 867-882. Wang, F. S.; Chiou, J. P. Computation of optimal control for integral and differential-algebraic systems. Proceedings of the American Control Conference, June 1995, 3933-3934. Wang, F. S.; Chiou, J. P. Nonlinear optimal control and optimal parameter selection by a modified reduced gradient method. Engineering Optimization, 1997, in press. Weinrb, A.; Bryson, A. E. Optimal control of systems with hard control bounds. IEEE Trans. Autom. Control 1985, AC-30, 1135-1138.

Received for review March 26, 1997 Revised manuscript received August 15, 1997 Accepted August 20, 1997X IE9702486

X Abstract published in Advance ACS Abstracts, November 1, 1997.