Using dynamic optimization technique to study the operation of batch reactors.

Sun, Daim-Young ; Lin, Pi-Min ; Lin, Shu-Ping 等

INTRODUCTION

Many chemicals, including bio-reagents, polymer additives, and pharmaceuticals, are essentially difficult to produce by continuous processes, because of varied species and small production quantities. For such multiproduct plants, batch processes may provide needed flexibility. However, with the increase in demand for the above chemicals, how to improve the operating efficiency of batch reactors becomes a very meaningful concern. Unlike the continuous reactors using the heating or cooling flow rate (the manipulated variable) to regulate the deviation from a fi xed reaction temperature (the set-point), the batch reactors with variant reaction conditions are almost impossible to operate under the framework of conventional PID loops. Consequently, the operating policies of a batch reactor are often determined off-line, based on the operating purposes and the imposed constraints. For most batch reactor operations, two distinct purposes are often considered:

1. Acquire the maximum yields of the desired product in a fi xed operating period under some constraints;

2. Reach the specified yield for the desired product(s) in the minimum operating period under some constraints.

The operating constraints in the above expressions usually arise from the physical bounds on the control variable(s), safety limitations, and environmental regulations. In optimal control studies, the first expression can be formulated as a maximum yield (or conversion) problem, while the second one as a minimum operating time problem. By comparing with the solutions from the above optimal control problems (OCPs), the operators can easily assess whether or not the real production qualities are satisfactory and further consider how to modify the practical operating policies to make the processes more better.

Basically, to solve the above OCPs is a challenging task due to the highly non-linear, non-convex and multimodal nature derived from chemical systems. In the past decades, most developed solution techniques are based on one of the following theories:

* Pontryagin's minimum principle;

* Dynamic programming (DP);

* Non-linear programming (NLP).

For the method based on the minimum principle (Pontryagin et al., 1962), an OCP is converted into a two-point-boundary-value problem with state and co-state equations. The optimal control policy is then determined by iteratively integrating these equations. However, it is well-known that such integrations are notoriously diffi cult, especially when there are path constraints in the problems. Arising from the principle of optimality, the recent DP-based methods (Luus, 1989, 1990) still execute the computations in backward manner. Because of the lack of system's historical information, such kinds of methods are thus forced to use a great deal of time grids and state grids to refine the solution. Although the exhaustive computations make DP-based methods have the best numerical quality, however, they also confront the situation of too slow convergence. In recent years, NLP-based techniques almost become the mainstream in solving optimal control problems. The methods are basically composed of two primary sections: parameterization (or discretization) and optimization. In parameterization, an OCP with infinite dimensions is converted into a NLP problem with finite dimensions by using a trial function with unspecified coefficients. The decision variables in the coming optimization are therefore in terms of these coefficients. According to the portions to be parameterized in an OCP, the NLP-based methods may be further classified into two categories: (1) complete parameterization (CP) (Biegler, 1984; Cuthrell and Biegler, 1987); and (2) control vector parameterization (CVP) (Vassiliadis et al., 1994a,b, 1999; Carrasco and Banga, 1997; Pham, 1998; Huang et al., 2002; Sun and Chen, 2004). State variables and control input(s) in CP are both parameterized. As a result, the state equations and the related constraints are approximated by a set of algebraic equations. The converted problem is easy to optimize by any gradient-based NLP technique, such as successive quadratic programming (SQP). Although CVP is a very straightforward notion, its major difficulty comes from the dissimilarity between the original system and its converted counterpart. Such difficulty becomes even clearer for highly non-linear systems. For resolving the problem, the most general way is to increase the order of trial function. Because the size of the discretized problem is expanded in this way, therefore, the solution is very easy to converge to a local region. Furthermore, since the path constraints are only satisfied at the discretized variables, the validity inside the solutions will be hard to verify. In CVP, the control input is the only portion to be parameterized. Gradient-based techniques (Vassiliadis et al., 1994 a, b, 1999; Huang et al., 2002) and direct search techniques (Carrasco and Banga, 1997; Pham, 1998; Sun and Chen, 2004) are used to optimize the unspecified parameters. By integrating the sensitivity equations to renew the decision variables, gradient-based methods are in general quickly converge the solution, especially when the initial guesses are fortunately close to the optimum. However, besides the local convergent problem, these techniques also confront the problem that some cases are indifferentiable. As for the direct search methods, the new decision variables are sequentially generated by a pseudo-random code. As well, the preservation of a decision variables will be merely dependent on whether or not the corresponding performance index has improvement. So, the use of the direct search techniques is simple and straightforward. Moreover, randomly producing vast amount of decision vectors also increases the opportunity for the solution to approach the global optimum.

In the past decade there have been many works about the optimization of batch operations primarily concerned with developing brand-new numerical methods and the related robustness and efficiency (Aziz and Mujtaba, 2002). In engineering aspects, Garcia et al. (1995) used the generalized reduced gradient (GRG) algorithm coupled with the golden search method to determine the optimal operating profile for some classical cases by applying piecewise linear input with fixed time interval in the framework of CVP. Luus and Okongwu (1999) consider how to operate the system with two inputs (the flow rates of heating and cooling fluids) with the help of iterative dynamic programming (IDP). Even if the results from IDP have very nice numerical quality, large computation burden often discourages its use. In addition, using excessive time and state grids also make the operating policy to be sophisticated and impractical. To acquire an implementable input from SQP, Aziz and Mujtaba (2002) suggested using the piecewise constant input with variable time interval in CVP. Like other gradient-based techniques, however, GRG and SQP still confront the local convergent problem. By applying the algorithm based on simulated annealing (SA), the authors have successfully solved some typical optimal control problems in the framework of CVP. Using piecewise linear input with variable time interval to rebuild the input policy, our numerical strategy (Sun and Chen, 2004; Sun and Lin, 2006) provides the solutions with very nice qualities, even if SA is just a convenient but old-fashioned technique. In this paper, we further apply the algorithm to explore the batch reactor operations.

In the rest of the article, we will first formulate the batch reactor operations as two optimal control problems based on the operating purposes. These problems are further parameterized via the concept of CVP. Subsequently, an introduction is provided to brief simulated annealing. We also discuss how to handle the path (equality and inequality) constraints that occurred in the OCPs. The model describing a non-linear continuous stirredtank reactor (CSTR) with the known global and local optima is first used to present the efficiency and robustness of the proposed algorithm. Finally, we will discuss the batch reactor operation based on two typical batch systems by applying the proposed method. The detailed computational steps of SA may be referred to in the Appendix.

PROBLEM FORMULATION

A batch reactor system is described by a set of ordinary differential equations (ODEs):

dx/dt = f(x(t), u(t), t) (1)

with the initial conditions:

x(0) = [x.sub.0] (2)

where x(t) and u(t), respectively, stand for an n x 1 state vector and m x 1 control vector over the period from the starting time, [t.sub.0], to the ending time, [t.sub.f]. Notably, Equation (1) results from the material and energy balance relations. The control inputs are also restricted by their physical lower and upper bounds

u [less than or equal to] u(t) [less than or equal to] u

In general, the reaction temperature and the coolant flow rate are the inputs that are often used to control batch reactors. Furthermore, the state variables could be restricted by L equality constraints deriving from thermodynamic and kinetic relations:

[h.sub.l](x(t), x(t), (t) = 0 l=1, ... , L

and/or K inequality constraints:

[g.sub.k](x(t), x(t),t) [less than or equal to] [g.sub.k] k=1, ... , K

where [g.sub.k] denotes the upper bound for [g.sub.k], the k-th inequality constraint. Because of being prescribed by safety considerations and/or environmental regulations, [g.sub.k] is thus never changed over sequent searches. According to the mentioned operating purposes, the above equations, from Equation (1) to Equation (5), can be further organized as maximum yield problem and minimum operating time problem as follows:

* Maximum yield problem

[min.sub.u(t) J = [-x.sub.desired](x([t.sub.f))

subject to:

dx/dt = d(x(t), u(t), t) (6b)

x(0) = [x.sub.0] (6c)

[h.sub.l](x(t), x(t), t) = 0 l=1, ... , L (6d)

[g.sub.k](x(t), x(t), t) [less than or equal to] [g.sub.k] k=1, ... , K (6e)

u [less than or equal to] u(t) [less than or equal to] u (6f)

where [x.sub.desired](x([t.sub.f])) denotes the yields of the desired product at the specified ending time [t.sub.f]. Without loss the generality, the performance index in this paper is presented as the form of Equation (6a). Other types of performance index can be easily recast as this form via augmenting the state space. Furthermore, Equation (6a) is modified as minimization by changing its sign to unify the following formulations.

* Minimum operating time problem

[min.sub,u(t)] J = [t.sup.*.sub.f] (7a)

subject to:

dx/dt = f(x(t), u(t), t) (7b)

x(0) = [x.sub.0] (7c)

[h.sub.l](x(t), x(t), t) = 0 l=1, ... , L (7d)

[x.sub.k](x(t), tx(t), t) [less than or equal to] [g.sub.k] k=1, ... , K (7e)

[x.sub.desired]([t.sup.*.sub.f]) = [x.sup.*] (7f)

u [less than or equal to] u(t) [less than or equal to] u (7g)

[t.sup.*.sub.f] [member of] [[t.sub.f, min] [t.sub.f, min]] (7h)

where x* denotes the yields of the desired product to be reached at the end of the operation.

CONTROL VECTOR PARAMETERIZATION

To optimize the above problems under the framework of control vector parameterization, the control portions, in Equations (6a) to (6f) and Equations (7a) to (7h), are needed to be parameterized. Here, to simplify the use of the notation, we only present the formulations for the single input systems, namely m = 1. The extension to the systems with multiple inputs is very straightforward. Thus, Equations (6a) to (6f) and Equations (7a) to (7h) are parameterized as follows:

* Maximum yield problem

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8a)

subject to:

dx/dt = f(x(t), [u.sub.i](t, [xi], [omega]), t) (8b)

x(0) = [x.sub.0] (8c)

[h.sub.l](x(t), x(t), t) = 0 = 1, l=1, ... , L (8d)

[g.sub.k](x.(t), x(t), t) [less than or equal to] [g.sub.k] k=1, ... , K (8e)

u [less than or equal to] [u.sub.t](t, [xi], [omega]) [less than or equal to] u (8f)

[t.sub.f] = [[omega].sub.N+1] (8g)

[for all]t [member of] [[omega].sub.t], [[omega].sub.l+1]) i=1, ... , N (8h)

* Minimum operating time problem

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9a)

subject to:

dx/dt = f(x(t), [u.sub.i](t, [xi], [omega], t) (9b)

x(0) = [x.sub.0] (9c)

[h.sub.l](x(t), x(t), t) = 0 l= 1, ... , L (9d)

[g.sub.k](x(t), x(t), t) [less than or equal to] [g.sub.k] k=1, ... , K (9e)

[x.sub.desired]([t.sup.*.sub.f] = x* (9f)

u [less than or equal to] [u.sub.t](t, [xi], [omega]) [less than or equal to] u (9g)

[for all] t [member of] [[omega].sub.t], [[omega.sub.t+1]) i=1, ... , N (9h)

[t.sup.*.sub.f] = [[omega].sub.N+1] [member of] [[t.sub.f, min] [t.sub.f, max] (9i)

where N([greater than or equal to] 1) denotes the number of control (time) stage from [t.sub.0] to [t.sub.f], and [[xi].sub.i] denotes the control input applied at the i-th time grid, [[omega].sub.i]. Therein, [[omega].sub.1] and [[omega].sub.N+1] are corresponding to [t.sub.0] and [t.sub.f], respectively. [u.sub.i](t, [[xi].sub.i], [[xi].sub.i+1], [[omega].sub.i], [[omega].sub. i+1]) represents the control input applied within the duration from [[omega].sub.i] to [[omega].sub.i+1]. For facilitating the selection of [[omega].sub.N+1] in Equations (9a) to (9i), a range [[t.sub.f,min], [t.sub.f,max]] that [t.sub.f] probably locates can be assigned in advance by referring to physical insight of the operator to the system. Moreover, Carrasco and Banga (1997) and Aziz and Mujtaba (2002) report that applying the input in variable time interval is good for the numerical qualities. The same conclusion is also reached in previous investigations (Sun and Chen, 2004; Sun and Lin, 2006). Thus, this policy is still adopted in this study. Consequently, other time grids are considered as the variables to be optimized.

So far, the above OCPs have been completely converted into two non-linear programming problems. Therein, Equations (8a) to (8h) have 2N dimensions, since the starting time and the ending time are both specified. Equations (9a) to (9i) have 2N + 1 dimensions due to the unknown ending time. Furthermore, many functions (Vassiliadis et al., 1994a: Carrasco and Banga, 1997) can be used to interpolate the value of [u.sub.i](t), t [member of] [[[omega].sub.i], [[omega].sub.i+1]) in the optimization. Therein, the piecewise constant function is widely used by many workers for the numerical simplicity. However, such function is in general hard to implement for practical uses, because of the intrinsic discontinuities occurring at the time grids. Increasing the grid numbers may ease the discontinuity, however, the dimension of the converted problems will be inevitably enlarged. To diminish the discontinuous effect, piecewise linear function seems to be more feasible. A piecewise linear input applied in the interval of [[[omega].sub.i], [[omega].sub.i+1]) is expressed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Since studies (Carrasco and Banga, 1997; Sun and Chen, 2004; Sun and Lin, 2006) have pointed that the piecewise linear inputs may give better results in most cases, therefore, Equation (10) is further adopted in this study. Thus, from practical point of view, the continuous input in Equations (6) and Equations (7) are approximated by a series of ramp inputs, [u.sub.i](t, [[xi].sub.i], [[xi].sub.i+1], [[omega].sub.i], [[omega].sub.i+1]), i = 1, 2, ... ,N, with different actuating periods. Even the resultant optimal policy has very abrupt changes such as bangbang input, the piecewise linear inputs with variable time interval are very effective in refining the solution by properly adding the grid numbers (Carrasco and Banga, 1997; Sun and Chen, 2004).

After having finished the previous formulations, we can further apply simulated annealing to select the optimum decision vector. A brief introduction about simulated annealing is provided in the next section.

OVERVIEW OF SIMULATED ANNEALING

The so-called annealing is a slow-paced cooling scheme. When a structured solid, composed of atoms, molecules or ions, is heated to some extent, the binding strength between composed particles starts to weaken. Further being heated, the solid will melt. Once the thermal energy is slowly released, these particles will gradually solidify toward a new but stable configuration. However, even at the same thermal state, the established configurations might be very diverse, which is the so-called meta states. Theoretically, the configuration with the highest stability can be eventually reached, provided the cooling scheme is slow enough.

Kirkpatrick et al. (1983) formulated the above procedure as a metaheuristics, called simulated annealing, by analogizing the melted particles and their associated energy states to the decision variables and the corresponding performance indexes. Suppose the performance index to be minimized is improved, it takes for granted that the corresponding decision variables are automatically accepted. However, even the performance index has no improvement, the corresponding decision variables are not directly abandoned but conditionally accepted as the seeds for escaping from the local regions. In SA, the acceptance of the unimproved solution is related to the following probability:

P = exp[[J.sup.i+1]([x.sub.i+1]) - [J.sup.i]([x.sub.i]/[T.sub.a]] (11)

where [T.sub.a] is the artificial annealing temperature, [J.sub.i+1] is the unimproved performance index from the current decision vector [x.sub.i+1], [J.sup.i] is that from xi, the previously accepted decision vector. P is then compared to an uniformly random number, P, ranging form [0, 1]. Once P is greater than P, [x.sub.i+1] is accepted to replace [x.sub.i] for further searches. Otherwise, [x.sub.i+1] is rejected. Such procedure is named as Metropolis criterion. Obviously, SA accepts the solutions trapped in local regions as the seeds so that the direction of convergence may go toward the global optimum. Of course, the cost is to prolong the time to converge.

Modified from the algorithm proposed by Goffe et al. (1994), the detailed computational steps of this study are in the Appendix. Therein, proportionally reducing the temperature is adopted as the annealing schedule, namely:

[T.sup.new.sub.a] = [gamma][T.sup.old.sub.a] (12)

where [gamma](< 1) is reducing factor specified by the user. The exploration about the effects caused by the parameters used in the algorithm, such as initial annealing temperature and the reducing factor and iterations before reducing temperature, is provided in the first example.

THE TREATMENT OF EQUALITY CONSTRAINTS

Because the time derivative terms of the state variables remain intact after parameterizing the control variable, the ordinary differential equations and the equality constraints in DOPs consist of a set of differential-algebraic equations (DAEs). In general, such systems are highly stiff. Whether or not a set of DAEs could be integrated by a numerical code is related to its index. The index of a set of DAEs is defined as the minimal number of times that the algebraic equations must be differentiated to obtain the standard ODEs. Thus, for an optimal control problem described by DAEs, the index is equivalent to the number of differentiation required to obtain a standard ODEs for the control variable (Dadebo and McAuley, 1995). Therefore, a system consisting of pure ODEs is obviously a special set of DAEs whose index is 0. The so-called high-index DAE means its index value is at least 2. For the high-index DAEs, inconsistency of initial conditions, and fail (or drift) of integration might happen in computing. So far, no general codes are available for the high index systems. Currently, DASSL (Petzold, 1982; Brenan et al., 1989) and DVERK (Hull et al., 1976) are two popular solvers for the systems with index 1 and 0. For applying these integrators, therefore, the index of the DAEs in this study is presumed to be less than or equal to 1. Furthermore, since the equality constraints can be simultaneously integrated, the path constraints in CVP are automatically satisfied by the state paths. As a result, the DAE solver provides an unified and convenient approach to handle the equality constraints.

THE TREATMENT OF INEQUALITY CONSTRAINTS

To ensure complete satisfaction of the inequality constraints over the entire time horizon, Mekarapiruk and Luus (1997) introduce additional state equations:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with initial condition [x.sub.n+k](0) = 0 where k = 1, ... , K, to measure the violations. As a result, the K inequality constraints are transferred into K end-point constraints. Using the similar idea, Aziz and Mujtaba (2002) develop the differentiable square-type function to accumulate the violations.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Sun and Lin (2006) replace the right-hand side of Equation (13) by the normalized function to increase the sensitivity.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For the cases with [g.SUB.k] = 0, Equation (15) is modified as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [alpha] is the parameter with a very small positive value. In this study, Equation (15) or (16) is still adopted to measure the violation of the path constraints. K end-point constraints and other end-point constraints, such as Equation (7f), are afterward augmented in the original performance index as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [[theta].sub.k] > 0 is the factor for penalizing the violation of the k-th inequality constraints. Thus, these end-point constraints may gradually shrink to a acceptable tolerance with the optimization.

CASE STUDIES

In this section, the non-linear CSTR model is first used to demonstrate the efficiency and robustness of the algorithm in acquiring the global solution for non-linear chemical systems. Then, this method is applied to explore the operations of batch reactors based on pure kinetic model and rigorous reactor model. By referring the temperature profile acquired from the optimization of the pure kinetic model, the engineers may improve the design of batch reactors. The optimization of the rigorous reactor model provides the information to evaluate the operating conditions of an already designed reactor. Notably, for a discretized OCP, the solution acquired by any NLP technique is merely suboptimal, since it is just very close to the global solution. For completing the above optimizations, all the computations are executed on a Pentium/4 personal computer with 3.0 GHz speed.

Non-Linear CSTR Model

This model is governed by the following differential equations:

d[x.sub.1]/dt = -(2 +_ u)([x.sub.1] + 0.25) + ([x.sub.2] + 0.5) exp [25[x.sub.1]/ [x.sub.1] + 2] (18a)

d[x.sub.2]/dt = 0.5 - [x.sub.2] - ([x.sub.2] + 0.5) exp [25[x.sub.1]/[x.sub.1] + 2] (18b)

d[x.sub.3]/dt = [x.sup.2.sub.1] + [x.sup.2.sub.2] + 0.1[u.sup.2] (18c)

with the initial condition

x(0) = [0.09 0.09 0]' (19)

where [x.sub.1] and [x.sub.2] represent the deviations from dimensionless steady-state temperature and concentration. The control input in this case is unbounded. Therefore, the problem is to fi nd the optimal input u to minimize the performance index [x.sub.3]([t.sub.f]) at [t.sub.f] = 0.78.

Using control vector iteration based on the minimum principle, Luus and Cormack (1972) have reported that there exist one local minimum, J = 0.24425, and one global minimum, J = 0.133094. Using IDP with 78 time stages and 19 iterations, Luus (1990) got an optimal value of 0.13321. The solution is then modified to 0.13312 by linear extrapolation. Luus (2000) further reports that the value may be further refined to 0.133096, which is within 0.0015% of the known best by IDP with 20 stages and 200 iterations. Luus (1990) also noticed that even if a large number of testing points are used in IDP, the possibility that gets a local optimum cannot be ruled out. And, to rerun the problem with different grid specifications is required for minimizing the chance of obtaining a local optimum. These opinions hint that the optimization of this case is sensitive to the initial conditions of the optimizers.

As previously indicated, simulated annealing conditionally accepts the unimproved solutions for departing from the local regions. The acceptance of an unimproved solution is dependent on whether or not the probability obtained from Equation (11) is larger than or equal to an uniform random number in the range of [0, 1]. In such a procedure, the initial value of [T.sub.a] will definitely affect the convergent speed. By setting the initial input as 2.5 and time stages as 5, the changes of the performance index for [T.sub.a] initiated at different values and cooled at the rate of 10% are shown in Figure 1. Obviously, the unimproved solutions are easy to be accepted as [T.sub.a] is in high level. The chance that escapes from the local region is thus greatly increased. However, excessive computations also make the convergence to be slow. In other respect, a small probability caused by a low [T.sub.a] will reject most unimproved solutions. Under such condition, simulated annealing and quasi Monte Carlo method will become very much alike. The changes of performance index with different annealing schedule as [T.sub.a] is started at 10 is presented in Figure 2. We can see that although all these solutions gradually get close to each other, the slow cooling rate will take much longer time in converging the solution. Figure 3 demonstrates the effect of different [N.sub.T] , the iteration numbers before reducing the temperature, on the performance index. It implies that a moderate [N.sub.T] is essentially required to quickly converge the solution to an numerical accuracy. Basically, the above shows that SA never trap these solutions to around 0.24425, but quickly converge to the vicinity of 0.1331. Because of being insensitive to the initial conditions of the optimization, therefore, simulated annealing behaves more robust than dynamic programming in this case. By taking [N.sub.T] = 30, [gamma] = 0.5 and [T.sub.a] = 0.1, the performance index acquired by different time stages are list in Table 1. All these values are very close each other. The trajectory of [x.sub.1](t) and [x.sub.2](t) and the control profile by SA with 8 stages are shown in Figures 4 and 5. Next, we further apply this method to investigate the operation of the batch reactors.

Pure Kinetic Model

The target of this model is to formulate the following consecutive-competitive reaction:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where A and B are raw materials, P and S denote the desired product and the waste or the low-priced product that should be inhibited. Obviously, the reaction temperature and the operating time are two factors that affect the conversion to P. As proposed by Burghardt and Skrzypek (1974) and studied by Garcia et al. (1995), the material balances for this system are given by:

d[x.sub.1]/dt = -([k.sub.1][x.sub.1][x.sub.2)V (20a)

d[x.sub.2]/dt = (-[k.sub1][x.sub.1][x.sub.2] - [k.sub.2][x.sub.2][x.sub.3])V (20b)

d[x.sub.3]/dt = (-[k.sub.1][x.sub.1][x.sub.2] - [k.sub.2][x.sub.2][x.sub.3])V (20c)

d[x.sub.4]/dt = ([k.sub.2][x.sub.2][x.sub.3])V (20d)

with the initial conditions

x(0) = [ 1. 1. 0. 0. ]' (20e)

[FIGURE 1 OMITTED]

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

where [x.sub.1], [x.sub.2], [x.sub.3], and [x.sub.4], respectively, denote the concentrations of A, B, P and S. The total reaction volume, V, is equal to 1l . The reaction rate constants, [k.sub.1] and [k.sub.2], in [mole.sup.-1][s.sup.-1] are expressed as:

[k.sub.1] = 1.667 x [10.sup.3][e.sup.-6.688 x 4/R(T+273)] (20f)

[k.sub.2] = 1.667 x [10.sup.3][e.sup.-8.360 x 4/R(T+273)] (20g)

where T in [degrees]C, the reaction temperature, is restricted by:

302 [les than or equal to]T(t) [less than or equal to] 352 (20h)

R denotes the gas constant. This model is further treated as two optimal control problems based on the previous formulations.

Maximum Yield Problem

The aim of maximum yield problem for this model is to optimize the temperature profile so that not only the yields of P is maximized at the end of the operating time [t.sub.f] = 6000 s, but Equations (20a) to (20h) are simultaneously satisfied.

Using GRG algorithm coupled with golden section method, Garcia et al. (1995) acquire the optimal temperature input and the performance in good agreement with the result from using the minimum principle (Burghardt and Skrzypek, 1974), where the maximum yields of P are 0.8664 for 5 stages and 0.8665 for 10 stages. By SA with [T.sub.a] = 50, [N.sub.S] = 30 and [N.sub.T] = 30, the maximum yields of P for the stages from 1 to 5 as shown in Table 2 are 0.8655, 0.8663, 0.8665, 0.8665 and 0.8666, respectively. Therein, when the stage numbers are increased to 3 or above, the yields of P exhibit almost the same values as reported in Garcia et al. (1995). The optimal temperature profiles for the above cases are plotted in Figure 6. Notably, the optimum yield of P is merely 0.13% lower than the best reported even just when one control stage is used. From the engineering point of view, such difference is almost negligible. On the other hand, when the control the stages are refined to 10, the conversion to P as shown in Table 3 has no further improvement. From the above facts, they seem to imply the possibility to approach the optimal conditions of this system by using small stages to change the reaction temperature. Thus, the operating complexity of this system can be greatly reduced. These data are plotted in Figures 7 and 8. Furthermore, Figures 6 and 7 both reveal a trend that the reaction will be initiated at a higher temperature level, when the stage numbers are increased. Then, the reaction rate has to be decelerated by steeply reducing the temperature to suppress the conversion of P. Finally, the temperature is gradually risen to another level, and kept constant to the end. This temperature variation appears when the stage numbers are increased to 3 or above. Similar observations are also reported by Burghardt and Skrzypek (1974) and Garcia et al. (1995).

[FIGURE 4 OMITTED]

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

Minimum Operating Time Problem

More P will be converted into S if the consecutive reaction is ended too late. It is reasonable to produce P with the specified yield as fast as possible in order to reduce the cost for sequently treating the low-priced by-product. Therefore, the minimum operating time problem for this model is to fi nd the optimum temperature profile so that not only Equations (20a) to (20h) are satisfied, but the specified yields of P is generated in the minimum operating time [t.sup.*.sub.f].

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

In the following discussions, the yields of P will be specified as 0.70, 0.75, 0.80, 0.85 and 0.866, respectively. With the increase of the specified yields of P, the system is needed to be given a much longer time to react. Thus, the ranges that the optimal operating time could be located are specified as [600, 1000] for 0.70, [600, 1300] for 0.75, [600, 1500] for 0.80, [2500, 3500] for 0.85 and [5000, 6100] for 0.866. As mentioned above, the operator may specify the ranges by referring to his physical insight to the systems. The related parameters used in SA are set as [T.sub.a] = 50, [N.sub.T] = 30 and [N.sub.S] = 30. Most computations in this case are completed within 5 s. The minimum operating times for the different yields of P with the different operating stages are list in Table 4. Therein, the time grid and the associated control are respectively shown in the third and the fourth column. The last grid marked with " * " in the third column denotes the time that the reaction is stopped. The conversions to S corresponding to the different conversions of P are list in the last column. The minimum operating time could be around 625 s, 862 s, and 1340 s, as the yields of P are specified as 0.70, 0.75 and 0.80. If the control stages are refined to 10, the operating time can be slightly improved to 623.16 s for 0.70, 861.8 s for 0.75, and 1337.5 s for 0.80. Notably, as shown Figure 9, the isothermal policy is favourable to the operation in the above cases. Once the yield of P is specified as 0.85, the operating time is minimized to 3190.4 s, 3189.0 s and 3186.2 s for 1, 3 and 5 stages. It can be further improved to 3185.3 s and 3182.0 s, as the control stages are refined to 8 and 10. Because the yield of P at 0.85 is favourable to the generation of S, the temperature profiles shown in Figure 10 are very different from those in Figure 9. We can fi nd that constant temperature profile and variable temperature profile both can be used to generate such concentration of P. When 8 stages are used, the shape of the temperature profile is very closed to on-off type. The temperature variation appear when 10 stages are used. If the specified yield of P is further increased to 0.866, which is a value that is very close to the solution from the previous maximum yield problem, the operating time is minimized to 6053 s for 3 stages, 6034 for 5 stages and 6037 for 10 stages. The above temperature profiles are plotted in Figure 11. The temperature variation appears as the stage numbers are increased to 5 or above.

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

Rigorous Reactor Model

Studied by Aziz and Mujtaba (2002), this model describes a first-order consecutive exothermic reaction in a jacketed batch reactor shown in Figure 12.

[FIGURE 11 OMITTED]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The material balances around the reactor are given as:

d[x.sup.1]/dt = -[k.sup.1][x.sub.1] (21a)

d[x.sup.2]/dt = [k.sup.1][x.sub.1] - [k.sub.2][x.sub.2] (21b)

d[x.sup.3]/dt = [k.sub.2][x.sub.2] (21c)

The energy balances around the reactor, the reactor wall and the jacket are described as:

d[x.sup.4]/dt = 193.4524[k.sub.1][x.sub.1] + 35.7143[k.sub.2][x.sub.2] - 8.8923([x.sub.4] - [x.sub.5] (21d)

d[x.sub.5]/dt = 33.1978([x.sub.4] - [x.sub.5] - 38.7940([x.sub.5] - [x.sub.6]) (21e)

d[x.sub.6]/dt = u/0.53(298 - [x.sub.6]) + 19.2925([x.sub.5] - [x.sub.6]) (21f)

[FIGURE 12 OMITTED]

where [x.sub.1], [x.sub.2] and [x.sub.3] denote the concentration of A, P, and S, respectively; [x.sub.3], [x.sub.4] and [x.sub.5] are the temperatures in K for the contents, the reactor wall and the jacket. The initial conditions of this system are:

x(0) = [ 0.975 0.025 0 350 373 300 ]' (21g)

The reaction rate constants in material balance relations are estimated as follows:

[k.sub.1] = 4.38 x [10.sup.4][e.sup.-3.49 x [10.sup.7]/(8314x4) (21h)

[k.sub.2] = 3.94 x [10.sup.5][e.sup.-4.65 x [10.sup.7]/(8314x4) (21i)

u(t), the coolant flow rate, in [m.sup.3]/h is bounded by:

0 [less than or equal] u(t) [less than or equal to] 9 (21j)

Four different constraints will be individually imposed on the above model.

C1: [x.sub.4]([t.sub.f]) [less than or equal to] 320 K.

C2: [x.sub.4](t) [less than or equal to] 370 K and [x.sub.4]([t.sub.f]) [less than or equal to] 320 K

C3: [x.sub.4]([t.sup.f]) [less than or equal to] 320 K and [x.sub.3]([t.sub.f]) = 0.1.

C4: [x.sub.4](t) [less than or equal to] 370 K, [x.sub.4]([t.sub.f]) [less than or equal to] 320 K and [x.sub.3]([t.sub.f]) = 0.1

Maximum Yield Problem

The aim of maximum yield problem for this model is to seek the optimum coolant flow rate which satisfies Equations (21a) to (21j) and the imposed constraints so that the yields of P is maximized at the end of the operation, [t.sub.f] = 3.5 h.

Table 5 summarizes the solutions for this model restricted by the above constraints by taking 2, 3 and 5 control stages and setting [T.sub.a] as 50. Most conversions to P are very closed to those from SQP shown in Table 6. Furthermore, the refined solutions in which the stages are increased to 10 are tabulated in Table 6 to compare with those from SQP. The profiles for the coolant flow rate and the reaction temperature for different constraints in this problem can be seen in Figures 13 and 14.

For C1 case, it is an end-point constraint that the reaction temperature should be ended at 320 K. To maximize the yields of P, the reaction temperature during the operation is thus allowable to be as high as possible. In practice, C1 could be the situation that the price of A is costly or the high temperature in reactor is considered to be endurable. In the computation, the highest temperature appears around 383.4 K, and the maximum yields of P is 0.6534. Furthermore, the yields of S in C1 is the largest among the four constraints, since the high reaction temperature simultaneously stimulate its generation.

The next consideration in C2 is to prevent the reaction temperature being raised too high so as to endanger the reactor. Consequently, the reaction temperature not only ends at 320 K, but has no permission to go up to 370 K during the operating period (a path constraint). Using square-type function to accumulate the violations, Aziz and Mujtaba (2002) partition a discretized DOP into more sub-problems, the so-called path-division technique, to ensure the satisfactions of the path constraints. In principle, the total violations of path constraints could be reduced by asking the constraints to be satisfied in each sub-problem. They report that the optimum conversion to P and the corresponding yield of S in one-path-interval method (without partition) are 0.6253 and 0.0991, while increasing to 0.6401 and 0.1161 in two-path-interval method. Because the results by 2 stages are impossible to converge to a satisfactory accuracy, only the results from 3 and 5 stages are shown in Table 5. Spending 15 s of CPU time, SA converges the yield of P to 0.6341 and 0.6350, respectively. Their differences with respect to the best value are within 1%. Even when the stages are refined to 10, the yield of P and the corresponding yield of S are improved to 0.6421 and 0.1214 within 19 s of CPU time.

[FIGURE 13 OMITTED]

[FIGURE 14 OMITTED]

[FIGURE 15 OMITTED]

In contrast to C2, the yield of S in C3 case is directly specified at a given level, namely S = 0.1. Such constraint might derive from the situation that the operator directly prescribe the amount of S, because of its costly expenditure in the sequent treatments or of the potential risk brought by S to the rector. By SA with 2, 3, 5 and 10 stages, the values, including 0.6227, 0.6276, 0.6270 and 0.6274 shown in Tables 5 and 6, are very close to 0.6249 from SQP. In Figure 14, we see the reaction temperature will go beyond the safety margin 350 K for rapidly producing P. As soon as the highest temperature, about 374.5 K, is reached, the reaction heat should be greatly removed to inhibit the conversion to S. By restricting the production of S, the extent over the safety margin for the temperature profile in C3 is smaller than that in C1.

All the above considerations are finally imposed on the model in C4 case. By using 3, 5 and 10 stages, the yield of P may be converged to 0.6270, 0.6280 and 0.6297. These results are very closed to 0.6253 from SQP. Although the results are only increased by 0.27%, 0.43% and 0.7%, using two or three stages may reduce the operating complexity.

Minimum Operating Time Problem

The aim of this problem is to fi nd the optimum coolant flow rate which satisfies Equations (21a) to (21j) and the imposed constraints so that the specified yields of P is generated in the minimum operating time [t.sup.*.sub.f]. In this problem, the yields of P at the end of the operation will be specified as 0.6. C1 and C2 will be respectively considered in the following discussions.

For C1 case, because of no path constraint imposed on the reaction temperature, the operation is expected to reach the specified yields of P as fast as possible, provided the reaction temperature is ended at 320 K. It is reasonable for an exothermic reaction to turn off the coolant flow in the initial stage. Thus, the reaction will be accelerated with the increase of the reaction temperature. Then, the coolant flow is fully opened to remove the reaction heat till the operating time is up. Such operations indeed create an on-off coolant flow. Setting [T.sub.a] as 50 and assuming [t.sup.*.sub.f] [member of] [0.5, 2.8], we acquire the minimum operating time to be 2.592 h and 2.578 h for 3 and 5 stages. Once the stages are refined to 10, the operating time can be further minimized to 2.40 h. In Figure 15, the changes of the coolant flow rate for 10 stages appears almost on-off form. The reaction temperature shown in Figure 16 initially ascends to a highest point, about 403 K, then is all the way down to 320 K. Additionally, because the reaction is allowed to be higher than 370 K, the productions of S is certainly increased. These data are shown in Table 7. Compared with the results by SQP (Aziz and Mujtaba, 2002), the optimal operating time by SA has been reduced by 2.4%, but the conversion to S is slightly increased by 1.0%.

In C2 case, we first solve the operating time for this problem to be 2.98 h by using 8 stages. The stages numbers are then increased to 10, the operating time is further minimized to 2.89 h. Because there exists the constraint on the highest reaction temperature, part of the reaction temperature profile in Figure 16 is carefully maintained at 370 K to increase the production of P. The operating time of this problem is thus longer than that in C1 due to the lower temperature section. Although the operating time compared with C1 is increased by 20%, however, the yield of S is reduced by 24.5%.

CONCLUSION

In this study, we explored the batch reactors operations by global optimization technique, simulated annealing. Based on the considerations in practical operation, the operations are first formulated two optimal control problems, including: a maximum yield problem and a minimum operating time problem. By discretizing the control inputs, the optimal control problems are converted into non-linear programming problems. The embedded parameters in discretized input can thus be optimized by the proposed algorithm. How to handle the equality and inequality constraints is also included into this article.

[FIGURE 16 OMITTED]

By using the CSTR model with the well-known local and global optimum, we first explore the effects caused by changing the parameters, including: the artificial annealing temperature [T.sub.a], the reducing factor in annealing schedule [gamma] and the iteration numbers before reducing the temperature [N.sub.T], on the convergence. The results show that SA is effective in excluding the local solutions, as well as less sensitive to the initial conditions in optimization. In conclusion, the basic strategy used by simulated annealing is to accept the unimproved solutions as the seeds for escaping the local regions. Therefore, the cost is certain to increase the computation efforts. On the other hand, once the degree of acceptance for the unimproved solutions becomes very small, SA and pseudo Monte Carlo method behave very much alike.

Based on the above formulations and discretization, then, four batch operation problems derived from the pure kinetic model and the rigorous reactor model are considered in this study. Therein, the pure kinetic model is only composed of material balance and the related kinetic relations. The solutions from this model may provide the optimum temperature trajectories that are helpful to the detailed design of the batch reactor. The rigorous reactor model consists of kinetic relations and material/energy balance. Four limitations on the operation are formulated as the path constraints, and also considered into the problems. The solutions from this model provide the optimal coolant flow rate which can be referred by the operation of an already designed reactor.

The above results verify that using simulate annealing to determine the optimal operating policy of batch reactors is simple, convenient and robust. Although SA has the critique of spending much longer time in exhaustive searches, its ability in search of global solution has been shown in this study. Therefore, the authors believe that the advancement of computers will gradually diminish such critique, and direct search techniques in the foreseeable future will be very promising in the solution of optimal control problems.

APPENDIX

Based on the Goffe's algorithm (Goffe et al., 1994), the detailed steps for determining optimal operating policy of batch reactors are stated as follows:

Step 1. Input problem's data and assign the parameters related to the computation, including:

1. Numbers of control (time) stages, N.

2. Artificial annealing temperature, [T.sub.a].

3. Numbers of the time grids, [N.sub.P] = P + 1.

4. Numbers of cycles at a given temperature, [N.sub.S].

5. Numbers of iteration before reducing temperature, [N.sub.T]. Noted that [N.sub.P] in original algorithm is set as the dimensions of the optimization problem. To increase the convergent speed, [[xi[.sub.i] and [[omega].sub.i] in our calculations are simultaneously updated. Furthermore, [N.sub.S] and [N.sub.T] are adequate for setting as 30, unless special stiff conditions are found in the problems, such as state path constraints.

Step 2. Initiate [N.sub.P] time grids and their associated control inputs by:

[[omega].sub.i] = (i - 1) * [t.sub.f]/P (A-1a)

[[xi].sub.i] = u + u/2 (A-1b)

where i = 1, ... , [N.sub.P]. The initial lower bound [[omega].sub.i] and the upper bounds [[omega].sub.i] of [[omega].sub.i] are given by:

[[omega].sub.i] = [[omega].sub.i-1] + [[omega].sub.i]/2 (A-2a)

[[omega].sub.i] = [[omega].sub.i] + [[omega].sub.i+1]/2 (A-2b)

Thus, the initial searching region for each [[omega].sub.i] is defined as:

[DELTA][[omega].sub.i] = [[omega].sub.i] - [[omega].sub.i] (A-2C)

Notably, [DELTA][[omega].sub.1] is always set as 0 since [[omega].sub.1] is correspondent to [t.sub.0]. Meanwhile, [DELTA][[omega].sub.P+1] is also fixed as 0 for the maximum yield problems whose [t.sub.f] specified. For the minimum operating time problems, the initial guess of [[omega].sub.P+1] is given as:

[t.sub.f,min] + [t.sub.f,min]/2 (A-3)

and [DELTA][[omega].sub.P+1] = [t.sub.f,max] - [t.sub.f,min.]

Step 3. Set counter l for [N.sub.T] as 1.

Step 4. Set acceptance indices [C.sub.i] as 0 where i = 1, ... , [N.sub.P].

The function of acceptance index is to record the numbers of the accepted performance indices in each grid during [N.sub.S] cycles at a given T.

Step 5. Set counter k for [N.sub.S] as 1.

Step 6. Set counter j for [N.sub.P] as 1.

Step 7. Simultaneously update the j-th time grid [[omega].sub.j] and the associated control input [[xi].sub.j] and the time grid by:

[[xi]'.sub.j] [left arrow] + [[GAMMA].sub.j] * [DELTA][u.sub.j] (A-4A)

[[omega]'.sub.j] [left arrow] [[omega].sub.j] + [[GAMMA].sub.j] * [DELTA][[omega].sub.j] (A-4b)

where [DELTA][u.sub.j](= u-u) means the physical bound for [[xi].sub.j] and [[GAMMA].sub.j] represents a pseudo-random number from Gaussian probability distribution with zero mean and one standard deviation. Notably, for the situation that [[xi].sub.j] or [[omega].sub.j] is out of its bound, [DELTA][u.sub.j] or [DELTA][w.sub.j], its new value is needed to regenerate by:

[[xi]'.sub.j] [left arrow] + [Y.sub.j] * (u - u) (A-5a)

or

[[omega]'.sub.j] [left arrow] [[omega].sub.j] + [Y.sub.j] * [DELTA][[omega].sub.j] (A-5b)

where [Y.sub.j] is a uniform distributed random number ranging from 0 to 1.

Step 8. Use the decision variables from Step 7 to calculate the performance index J by integrating the dynamic equations.

Step 9. Determine whether a new J and its corresponding decision variables would be accepted as follows:

* If J is better than J', the performance index from the accepted decision variables, J' with the corresponding decision variables is accepted as the new solution. Make [C.sub.j] [left arrow] [C.sub.j] + 1,

* Once J is also better than the current optimum [J.sup.old], replace [J.sup.new] by J, and record the associated decision variables as new optimum ones. Check the stopping criterion:

|J.sup.new] - [J.sup.old|/|[J.sup.old]| < [epsilon] (A-6)

where [epsilon] is the specified tolerance factor. If Equation (A-6) is satisfied, stop the computation. Otherwise, make [J.sup.old] [left arrow] [J.sup.new] go to Step 10. Note that J' and [J.sup.old] both have very huge initial values.

* If J is inferior to the current J', calculate the probability P as follows:

P = exp[- J - J'/[T.sub.a] (A-7)

If P is greater or equal to an uniformly random number P [member of] [0, 1], then J and the corresponding decision variables are retained for the next step. Otherwise, they are rejected. Make [C.sub.j] [left arrow] [C.sub.j] + 1. Then, go to Step 10.

Step 10. Make j [left arrow] j + 1. If j [less than or equal to] [N.sub.P] , then go back to Step 7, or go to

Step 11. Notably, for j = [N.sub.P], the current [t.sub.f] should be modified by:

[t.sub.f] [left arrow] [t.sub.f] + [GAMMA] * [DELTA][[omega].sub.P+1]

Being out of the range of [[t.sub.f,max], [t.sub.f,min]], the new [t.sub.f] may be regenerated by:

[t.sub.f] [left arrow] [t.sub.f,min] + Y * [DELTA][[omega].sub.P+1] (A-9)

Step 11. k [left arrow] k + 1. If k [less than or equal to] [N.sub.S], then go to Step 6, or go to Step 12.

Step 12. Modify the bounds for the control inputs and time grids as follows:

Step 12a. Use Equations (A-2a) to (A-2c) to rebuild the new bound for each time grid.

Step 12b. Calculate j-th frequency indices [IF.sub.j] = [C.sub.j]/[N.sub.S], and recast the j-th new searching regions as:

* IFj > 0.6

[DELTA][u.sub.j] [left arrow] [DELTA][u.sub.j] x [1. + [F.sub.j] - 0.6/0.4] (A-10a)

[DELTA][[omega].sub.j] [left arrow] [DELTA][[omega].sub.j] x [1. + [F.sub.j] - 0.6/0.4] (A-10b]

* 0.4 [less than or equal to] [IF.sub.j] [less than or equal to] 0.6

[DELTA][u.sub.j] [left arrow] u - u (A-11a)

[DELTA][[omega].sub.j] [left arrow] [omega].sub.j] - [[omega].sub.j] (A-11b)

* [IF.sub.j] < 0.4

[DELTA][u.sub.j] [left arrow] [DELTA][u.sub.j]/1. + (0.4 - [F.sub.j])/0.4 (A-12a)

[DELTA][[omega].sub.j] [left arrow] [DELTA][[omega].sub.j]/1. + (0.4 - [F.sub.j]/0.4 (A-12b)

The values of 0.6 and 0.4 in Equations (A - 10a) to (A - 10b) and (A - 12a) to (A - 12b) are made possible to maintain the average percentage of accepted changes at about one-half of the total number of the cycles (Goffe et al., 1994; Li et al., 2000).

Step 13. l left arrow] 1, + 1 if l [less than or equal to] [N.sub.T], go to Step 4. Otherwise, go to

Step 14.

Step 14. Make [T.sub.a] [left arrow] [gamma][T.sub.a], where [gamma] is the temperature reduction factor with default value 0.9. If the above scheme has iterated 30 times, print out the current optimum [J.sup.new] and its corresponding control vector, then stop. Otherwise, Make J [left arrow] [J.sup.opt], [[xi].sub.j] [left arrow] [[xi].sub.j] [left arrow] [[xi].sup.opt.j] and ([[xi].sub.j], [[omega].sub.j]) [left arrow] ([[xi].sup.opt], [[omega].sup.opt.sub.j]), j = 1, ... , [N.sub.p]. Go back to Step 3.

ACKNOWLEDGEMENT

The authors would like to thank the anonymous reviewers for their constructive comments and suggestions.

Manuscript received May 21, 2006; revised manuscript received September 6, 2006; accepted for publication February 6, 2007.

REFERENCES

Aziz, N. and I. M. Mujtaba, "Optimal Operation Policies in Batch Reactors," Chem. Eng. J. 85, 313-325 (2002).

Bellman, R. E., "Dynamic Programming," Princeton University Press, Princeton, NJ, U.S.A. (1957).

Brenan, K. E., S. L. Campbell and L. R. Petzold, "Numerical Solution of Initial Value Problems in Differential-Algebraic Equations," North-Holland, NY, U.S.A. (1989).

Biegler, L. T., "Solution of Dynamic Optimization Problems by Successive Quadratic Programming and Orthogonal Collocation," Comput. Chem. Eng. 8, 243-248 (1984).

Burghardt, A. and J. Skrzypek, "Optimal Temperature Profiles in A Tubular Reactor for a System of Consecutive-Competing Reaction," Chem. Eng. Sci. 29, 1315-1331 (1974).

Carrasco, E. F. and J. R. Banga, "Dynamic Optimization of Batch Reactors Using Adaptive Stochastic Algorithms," Ind. Eng. Chem. Res. 36, 2252-2261 (1997).

Cuthrell, J. J. and L. T. Biegler, "On the Optimization of Differential-Algebraic Process Systems," AICHE J. 33, 1257-1270 (1987).

Dadebo, S. A. and K. B. McAuley, "Dynamic Optimization of Constrained Chemical Engineering Problems Using Dynamic Programming," Comput. Chem. Eng. 19, 513-525 (1995).

Garcia, V., M. Cabassud, M. V. Le Lann, L. Pibouleau and G. Casamatta, "Constrained Optimization for Fine Chemical Productions in Batch Reactors," Chem. Eng. J. 59, 339-241 (1995).

Goffe, W. L., G. D. Ferrier and J. Rogers, "Global Optimization of Statistic Functions with Simulated Annealing," Econometrics 60, 65-99 (1994).

Huang, Y. J., G. V. Reklaitis and V. Venkatasubramanian, "Model Decomposition Based Method for Solving General Dynamic Optimization Problems," Comput. Chem. Eng. 26, 863-873 (2002).

Hull, T. E., W. D. Enright and K. R. Jackson, "User Guide to DVERK - A Subroutine for Solving Nonstiff ODE's," Dept. of Computer Science, University of Toronto, ON, Canada 1000 (1976).

Kirkpatrick, S., C. D. Gelatt, Jr. and M. P. Vecchi, "Optimization by Simulated Annealing," Science 220, 671-680 (1983).

Li, P., K. Lowe, H. Arellano-Garcia and G. Wozny, "Integration of Simulated Annealing to a Simulation Tool for Dynamic Optimization of Chemical Processes," Chem. Eng. Process. 39, 357-363 (2000).

Luus, R. and D. E. Cormack, "Multiplicity of Solutions Resulting from the use of Variational Methods in Optimal Control Problems," Can. J. Chem. Eng. 50, 309-312 (1972).

Luus, R., "Optimal Control by Dynamic Programming Using Accessible Grids Points and Region Reduction," Hung. J. Ind. Chem. 17, 523-544 (1989).

Luus, R., "Optimal Control by Dynamic Programming Using Systematic Reduction in Grid Size," Int. J. Control 51, 995-1013 (1990).

Luus, R. and O. N. Okongwu, "Towards Practical Optimal Control of Batch Reactors," Chem. Eng. J. 75, 1-9 (1999).

Luus, R., "Iterative Dynamic Programming," Chapman & Hall/ CRC (2000).

Mekarapiruk, W. and R. Luus, "Optimal Control of Inequality State Constrained Systems," Ind. Eng. Chem. Res. 36, 1686-1694 (1997).

Petzold, L. R., "Differential/Algebraic Equations Are Not ODEs," SIAM J. Sci. Stat. Comput. 3, 367-384 (1982).

Pham, Q. T., "Dynamic Optimization of Chemical Engineering Processes by an Evolutionary Method," Comput. Chem. Eng. 22, 1089-1097 (1998).

Pontryagin, L. S., V. G. Boltyanskii, R. V. Gamkrelidze and E. F. Mishchnkov, "The Mathematical Theory Of Optimal Process," Inter-Science Publisher, Inc., NY, U.S.A. (1962).

Sun, D. Y. and C. L. Chen, "An Algorithm for Optimizing the Unsteady Chemical Processes by Simulated Annealing," J. Chem. Eng. Japan 37, 711-722 (2004).

Sun, D. Y. and P. M. Lin, "The Solution of Time Optimal Control Problems by Simulated Annealing," J. Chem. Eng. Japan 39, 753-766 (2006).

Vassiliadis, V. S., R. W. H. Sargent and C. C. Pantelides, "Solution of a Class of Multistage Dynamic Optimization Problems. 1. Problems without Path Constraints," Ind. Eng. Chem. Res. 33, 2111-2122 (1994a).

Vassiliadis, V. S., R. W. H. Sargent and C. C. Pantelides, "Solution of a Class of Multistage Dynamic Optimization Problems. 2. Problems with Path Constraints," Ind. Eng. Chem. Res. 33, 2123-2133 (1994b).

Vassiliadis, V. S., E. Balsa-Canto and J. R. Banga, "SecondOrder Sensitivities of General Dynamic Systems with Application to Optimal Control Problems," Chem. Eng. Sci. 54, 3851-3860 (1999).

Daim-Yuang Sun, (*) Pi-Min Lin and Shu-Ping Lin

Department of Chemical and Materials Engineering, National Chinyi Institute of Technology, Taiping City, Taichung County, 411 Taiwan R.O.C.

* Author to whom correspondence may be addressed. E-mail address: dysun@chinyi.ncit.edu.tw

Table 1. The performance indexes of the non-linear CSTR with
different control stages

Time stages J

5 0.133381
8 0.133129
10 0.133361
20 0.133133

Table 2. The maximum yield of P for the pure kinetic model where the
best known value of P is 0.8665

Control Time grids Associated control inputs
stages

1 (0, 6000) (317.3, 352.0)
2 (0, 3887, 6000) (302.9, 351.9, 352.0)
3 (0, 277, 3723, 6000) (339.3, 302.4, 351.7, 352.0)
4 (0, 304, 3685, 3842, 6000) (337.2, 302.0, 351.3, 352.0,
 352.0)
5 (0, 154, 1133, 3481, 3674, (351.7, 304.5, 309.7, 350.7,
 6000) 352.0, 352.0)

Control P (l [mole.sup.-1][s.sup.-1]) S (l [mole.sup.-1][s.sup.-1])
stages

1 0.8655 5.571 x [10.sup.-2]
2 0.8663 5.678 x [10.sup.-2]
3 0.8665 5.675 x [10.sup.-2]
4 0.8665 5.672 x [10.sup.-2]
5 0.8666 5.658 x [10.sup.-2]

Table 3. The distributions of reaction temperature and product
concentrations for the pure kinetic model in maximum yield
problem with 10 stages

t T A B P S

 0 352.0 1.0000 1.0000 0.0000 0.0000
 131 306.9 0.7354 0.7339 0.2631 0.0014
 874 308.6 0.3932 0.3829 0.5965 0.0102
1463 315.8 0.2816 0.2644 0.7011 0.0172
1647 318.0 0.2569 0.2374 0.7237 0.0194
2391 329.4 0.1840 0.1559 0.7878 0.0281
2702 334.9 0.1623 0.1307 0.8059 0.0317
3326 350.6 0.1279 0.0886 0.8329 0.0392
3630 352.0 0.1153 0.0726 0.8420 0.0426
4824 352.0 0.0883 0.0365 0.8599 0.0517
6000 352.0 0.0768 0.0203 0.8666 0.0565

Table 4. The optimization results for the pure kinetic model to achieve
the different specified yields of P in minimum operating time

P Control Time grids Associated inputs
(mole/l) stages

0.70 1 0, 628.0 * 351.2, 352.0
 3 0, 21.7, 544.7, 624.9 * 351.5, 351.8, 351.9,
 351.7
 5 0, 145.2, 162.0, 162.8, 351.9, 351.9, 351.7,
 592.8, 625.4 * 351.8, 351.8, 351.2

0.75 1 0, 861.8 * 349.1, 350.5
 3 0, 16.4, 386.2, 861.9 * 351.7, 352.0, 352.0,
 352.0
 5 0, 4.6, 23.9, 94.1, 184.1, 346.2, 351.6, 352.0,
 862.5 * 352.0, 352.0, 352.0

0.80 1 0, 1342.1 * 351.8, 351.8
 3 0, 1194.5, 1330.4, 1340.4 * 352.0, 351.9, 351.0,
 347.7
 5 0, 8.3, 68.4, 754.4, 770.0, 349.7, 351.8, 352.0,
 1339.8 * 351.8, 351.8, 352.0

0.85 1 0, 3190.4 * 351.9, 352.0
 3 0, 905.5, 930.9, 3189.0 * 348.0, 351.8, 351.9,
 351.9
 5 0, 629.6, 638.8, 763.0, 347.7, 349.6, 351.9,
 943.3, 3186.2 * 352.0, 352.0, 352.0

0.8666 3 0, 1972.9, 3454.7, 6053.9 * 309.0, 318.8, 351.8,
 351.9
 5 0, 705.3, 2676.3, 3397.7, 314.7, 304.3, 333.4,
 3457.5, 6034.7 * 350.2, 351.8, 352.0

P Control S
(mole/l) stages (mole/l)

0.70 1 0.0219
 3 0.0219
 5 0.0219

0.75 1 0.0282
 3 0.0286
 5 0.0285

0.80 1 0.0383
 3 0.0383
 5 0.0383

0.85 1 0.0565
 3 0.0560
 5 0.0560

0.8666 3 0.0566
 5 0.0564

Table 5. The optimization results for the rigorous reactor model to
acquire the maximum yield with different constraints

Constraint Control Time grids Associated inputs
type stages

C1 2 0, 1.71, 3.50 0.369. 0.027, 5.195
 3 0, 2.02, 2.29, 3.50 0, 0.009, 0.657, 0.096,
 7.634
 5 0, 1.50, 2.11, 2.45, 0, 0.515, 0.665, 0.372,
 2.60, 3.50 1.370, 8.869

C2 3 0, 1.73, 2.59, 3.50 0.080, 0.7316, 0.6377,
 8.5410
 5 0, 0.92, 1.42, 1.49, 0, 0.467, 0.649, 0.712,
 2.62, 3.50 0.653, 8.823

C3 2 0, 1.66, 3.50 0.5989, 0.0679, 4.5707
 3 0, 0.55, 2.88, 3.50 0.2958, 0.1087, 1.5862,
 8.2114
 5 0, 1.03, 1.34, 2.38, 0.1329, 0.4866, 0.6275,
 2.74, 3.50 0.6614, 1.6833, 8.1289

C4 3 0, 1.35, 2.33, 3.50 0.0218, 0.6856, 0.6085,
 5.9695
 5 0, 0.84, 1.62, 1.90, 0.0174, 0.2290, 1.4050,
 2.75, 3.50 0.0309, 2.1591, 5.9059

Constraint Control P (mole/l) S (mole/l)
type stages

C1 2 0.6457 0.1707
 3 0.6494 0.1700
 5 0.6500 0.1664

C2 3 0.6341 0.1093
 5 0.6350 0.1102

C3 2 0.6227 0.1000
 3 0.6276 0.1000
 5 0.6270 0.1000

C4 3 0.6270 0.1000
 5 0.6280 0.1000

Table 6. The yields of P and S for the rigorous reactor model
in maximum yeild problem with different operating constraints

 P S

Constraints SA SQP SA SQP T([t.sub.f])

C1 0.6534 0.6426 0.1674 0.1690 320.0
C2 0.6421 0.6401 0.1214 0.1161 320.0
C3 0.6274 0.6249 0.1000 0.1000 320.0
C4 0.6297 0.6253 0.1000 0.1000 320.0

Table 7. The minimum operating time [t.sub.f] * and the yields of S of
the rigorous reactor model with different operating constraints

 S [t.sub.f] *

Constraints SA SQP SA SQP

C1 0.1131 0.1071 2.404 2.46
C2 0.0854 0.0845 2.888 2.91