7392
Ind. Eng. Chem. Res. 2004, 43, 7392-7400
Optimal Temperature-Time Programming in a Batch Copolymerization Reactor D. Salhi, M. Daroux, C. Gentric, J. P. Corriou, F. Pla, and M. A. Latifi* Laboratoire des Sciences du Ge´ nie Chimique, CNRS-ENSIC, B.P. 451, 1 rue Grandville, 54001 Nancy Cedex, France
This paper addresses the optimization of a batch copolymerization reactor of styrene and R-methylstyrene. The reactor is equipped with a jacket for heat exchange and, thus, for temperature control. The performance index is the batch period, and the decision variable is the jacket inlet temperature. The constraints involved include the process model, the polymer quality, the final monomer conversion rate, the heating/cooling rate of the reactor and upper and lower bounds of the variables. The optimization method is a direct method with a sequential strategy known as control vector parametrization method. The optimal profiles are determined with great care in terms of treatment and fulfilment of constraints. The results show a substantial improvement with respect to the complete discretization method previously used. 1. Introduction The optimization problem considered here concerns a batch emulsion copolymerization reactor of styrene and R-methylstyrene. The reactor is equipped with a jacket for heat exchange and, thus, for temperature control. In the previous experimental work of Gentric (1997),1 a three-way valve allowed for manipulation of the temperature at the jacket inlet. The typical load of the reactor consists of water (dispersion medium), monomers that are slightly soluble in water, a watersoluble initiator, and an emulsifier. The contents of both the reactor and jacket are assumed to be perfectly mixed, and the control or decision variable is the jacket inlet temperature (Figure 1). The optimization problem is represented in Figure 2, where the objective is to maximize the productivity, which is expressed in terms of the determination of jacket inlet temperature profiles that minimize the batch period, tf, under specified constraints. The latter are of two types, namely, process model and operating constraints, and will be detailed later, as will be the optimizer developed here. 1.1. Previous Works. In our previous works,1,2 for the optimization stage, the control variable was the reactor temperature, and therefore, the heat balance equations in the reactor and in the jacket were not taken into account in the process model; rather, only the mass balance equations were considered. On the other hand, the model was described by a system of ordinary differential equations (ODEs) that was correct, but we will see a more appropriate representation in this work. The optimization method (open-loop control) was a direct method with a simultaneous strategy known as the complete discretization method, i.e., it is based on the parametrization of the state and the control variables by means of Lagrange polynomials. In addition, the method of orthogonal collocation on finite elements was mainly used. For control purposes, the heat balance equations were incorporated into the model. A control * To whom correspondence should be addressed. Tel: +33 (0)3 83 17 52 34. Fax:+33 (0)3 83 17 53 26. E-mail: latifi@ ensic.inpl-nancy.fr.
Figure 1. Schematic representation of the batch polymerization reactor.
Figure 2. Diagram of the optimization problem.
law was then designed in closed-loop control using nonlinear geometric control and a state observer based on an extended Kalman filter (EKF). The experimental implementation involved the tracking of the optimal reactor temperature profiles determined off-line in openloop control, i.e., set-point trajectories, and the control variable was the jacket inlet temperature. 1.2. Present Work. In the present work, the optimization objective is the same as before,2 but a different approach is employed. It is based on the use of the jacket inlet temperature, Tjin, as the control variable instead of the reactor temperature T. In this case, it is necessary to take the heat balance equations for the reactor and the jacket into account in the model. On the other hand, the total set of model equations is described by a hybrid system of ODEs that is more appropriate in this case. The optimization method used is a direct method with a sequential strategy known as control vector parametrization (CVP). A comparison is made between the optimization results obtained previously by the method of complete
10.1021/ie0498549 CCC: $27.50 © 2004 American Chemical Society Published on Web 09/25/2004
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7393
discretization2 and those obtained in the present work by the CVP method.
Radical entry into active particles
2. Process Model
Propagation
The first type of constraints involves those required by the process model. The overall process is described by a kinetic model, a molecular weight distribution model, and heat balance equations in the reactor and in the jacket. 2.1. Kinetic Model. The copolymerization of styrene and R-methylstyrene has been previously studied both in the bulk3 and in solution.4-9 Similarities with the homopolymerization behavior of styrene or R-methylstyrene are often emphasized. The kinetic model of an emulsion polymerization usually divided into three stages10 is detailed in ref 1, where different assumptions and subsequent models are discussed. The experimental validation was carried out by Castellanos11 and Gentric.1 The three stages of the model are as follows: (i) In the first stage, free radicals are produced in the aqueous phase by initiator decomposition. They are captured by the micelles swollen with monomer. The polymerization begins in these micelles. This stage, corresponding to particle heterogeneous nucleation, stops when all of the micelles have disappeared. (ii) Particle growth occurs during the second stage. Monomer diffuses rapidly from monomer droplets toward the particles, which are saturated with monomer as long as monomer droplets exist. This stage ceases when all of the monomer droplets have disappeared. (iii) The third and final stage is characterized by the decrease of the monomer concentration in particles. The following assumptions are used in the modeling of the kinetic process: (A1) Styrene and R-methylstyrene are both hydrophobic monomers; thus, only micellar nucleation is taken into account. (A2) Because of this hydrophobic character, the reactions of propagation, transfer to monomer, termination in the aqueous phase, and radical desorption are neglected. (A3) Termination in particles is considered to be very rapid compared to radical entry into particles; thus, it can be assumed that no more than one radical is present in each polymer particle (zero-one system).12 This allows us to write P•j ) N•. (A4) The maximum conversion rate considered is generally around 65%, so that gelation has not yet occurred.13 Thus, the gel effect is not included in the model to avoid unnecessary complexity.14 (A5) Coagulation is assumed to be negligible. The kinetic mechanism is then written as follows:
Ra ) 2fkdA
(1)
Particle formation R • + m f N•
Rn ) kcmmR •
(2)
Ri ) kcpNR •
Rt ) kcpN •R •
Rp ) kpMpP •j ) kpMpN •
(3)
(4)
(5)
Transfer to monomer P•j + M f M• + Pj
RtrM ) ktrMMpP •j ) ktrMMpN • (6)
The rate of particle formation can be derived from this mechanism. The quasi-steady-state approximation allows for the determination of the initiator radical concentration in the water phase as
R• )
Ra kcmm + kcpNp
(7)
The total number of particles is deduced from eq 7 and the rate of radical entry into the micelles (eq 2)
dNp RaNa ) kcmm dt kcmm + kcpNp
(8)
In the same way as Harada et al.,15 who studied the emulsion homopolymerization of styrene, we introduce a capturing efficiency of the particles with respect to the micelles
)
kcpns kcm
(9)
where ns is the aggregation number of a micelle, defined as
ns )
SNa m
(10)
Thus, during step 1
dNp ) dt
RaNa Np 1+ SNa
(11)
where, with the classical assumption that emulsifier molecules are adsorbed in monomolecular layers at the polymer particles surface, we can write
S ) S0 - kv(XM0)2/3Np1/3 with
[
36πMM2
ωP2Fp2(asNa)3
]
(12)
1/3
(13)
In the above expressions, X is the conversion rate, and M0 is the overall initial monomer concentration. The rate of monomer consumption is classically expressed as
Np dM ) -Rp ) -kpMp n j dt Na
Radical entry into inactive particles N + R • f N•
• P•j + M f Pj+1
kv )
Initiator decomposition A f 2R•
N + R• f N•
(14)
Here, because the objective is temperature control, only
7394
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004
the overall monomer concentration is described. Castellanos11 studied the copolymerization of styrene and R-methylstyrene at three temperatures (50, 65, and 85 °C) and four different initial monomer compositions (10, 25, 35, and 50% by mass of R-methylstyrene). He found experimentally that no composition change occurred during the polymerization under these conditions and showed that the global propagation constant could be written as
kp ) kpsty exp(afMS)
(15)
where kpsty is the styrene homopolymerization propagation rate constant, a is a constant, and fMS is the initial R-methylstyrene mole fraction. The monomer concentration in the particles, Mp, corresponds to the particle saturation during steps 1 and 2 and then decreases with increasing conversion rate during step 3
Mp ) Mpc )
(1 - Xc)FM [(1 - Xc) + XcFM/FP]MM X e Xc, steps 1 and 2 (16)
Mp )
(1 - X)FM
X > Xc, step 3
[(1 - X) + XFM/FP]MM
(17)
Rudin and showed experimentally that n j is equal to 0.5 for styrene/R-methylstyrene copolymerization. Castellanos11 confirmed this value, even at high conversion rates. 2.2. Molecular Weight Distribution Model. Here, the molecular weight distribution model is based on the tendency model developed by Villermaux and Blavier.16 It consists of three ODEs where the unknowns are the different moments of the molecular weight distribution. The integration of these equations leads to the numberand weight-average molecular weights and the polydispersity index given by Samanta7
dP dQ0 ) ) Rt + RtrM dt dt dQ1 ) L(Rt + RtrM) ) Rp dt
(18)
dQ2 ) 2L2(Rt + RtrM) dt
Ran j Np S Np +
RtrM ) ktrMMp
Np n j Na
(19)
Q2 M h w ) MM Q1 Ip )
(22)
M hw M hn
Remark. The rate constant for transfer to monomer is computed similarly to the propagation rate constant. Because no composition change occurs during this copolymerization under the conditions considered, a global constant is thus assumed. 2.3. Heat Balance Equations. Finally, the heat balance equations in the reactor and in the jacket, the contents of which are assumed to be perfectly mixed, are given, respectively, as
UA dT V∆H R + (T - T) )dt mrCp p mrCp j
(23)
dTj Fj UA ) (Tjin - Tj) (T - T) dt Vj FjVjCj j
(24)
The initial temperatures can be considered as known and equal to room temperature or as unknown parameters to be optimized. Here, they are included in the optimization process as unknown parameters. 2.4. Hybrid Representation of the Process Model. The resulting model is then a hybrid system of ODEs with two state modes denoted as 1 and 217
mode 1 x3 (1) ) f(1)(x(1),u(1),p) mode 2 x3 (2) ) f(2)(x(2),u(2),p)
(1) t ∈ [t(1) 0 , tf ]
(25)
(2) t ∈ [t(2) 0 , tf ]
(i) Each mode is characterized by state and control variables, parameters, a set of ODEs, and a time interval in which these equations are defined. Note that the upper bound of the time interval of mode 2 is the batch period tf ) t(2) f . Note also that a possible discontinuity in control is taken into account in these equations. (ii) The switching time from mode 1 to mode 2 is determined by a transition condition
(26)
which, in our case, is the time when the monomer conversion rate equals its critical value (see eqs 16 and 17). (iii) The variables of the two modes at the transition time are linked through the following transition function (1) (1) T(1)(x(1),u(1),x(2),u(2),p) ) 0 w x(2)(t(2) 0 ) ) x (tf ) (27)
(20)
and
Rp L) Rt + RtrM
Q1 M h n ) MM Q0
(1) L(1)(x(1),u(1),p) ) 0 w t(2) 0 ) tf
where
Rt )
(ODEs) provides the number- and weight-average molecular weights and the polydispersity index
(iv) A special case of this transition condition is the set of initial conditions of the initial mode, which is mode 1 (1) (1) (1) T0(1)(x(1), x(1) 0 ) ) 0 w x (t0 ) ) x0
(21)
The integration of these ordinary differential equations
(28)
3. Operating Constraints The second type of constraints is the operating constraints. They consist of (i) the lower and upper
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7395
bounds of the seven state variables, the single control variable, and the two parameters
form similar to the cost form (eq 34)
Ji ) Gi[x(tf),p] +
xLi e xi(t) e xU i , i ) 1, ..., 7 xT ) [M, Np, Q0, Q1, Q2, T, Tj] uL e u(t) e uU, u ) Tjin pLi e pi(t) e pU i , i ) 1, 2 pT ) [T0, Tj0]
(29)
(ii) the final monomer conversion rate
X(t(2) f ) ) 1 -
(2) x(2) 1 (tf )
x(1) 10
) Xf
(30)
(iii) the polymer quality expressed in terms of the number- or weight-average molecular weight
MM
(2) x(2) 4 (tf ) (2) x(2) 3 (tf )
MM
(2) x(2) 4 (tf )
5.1. Treatment of Constraints. To transform the constrained optimization problem (eqs 34-36) into a problem constrained only by the process model eqs 35, the following variables are defined i)Nc
(37)
i)Nc
νi G i ∑ i)1
G ) G0 +
)M h wf
(38)
i)Nc
| | | |
dx6 dT ) e hcr dt dt
(32)
Note that the final conversion rate and polymer quality constraints are terminal-point equality constraints, whereas the heating/cooling rate is an infinite-dimensional inequality constraint. The resulting optimization problem is then given by
min {J ) tf} u(t),p
(33)
subject to constraints 25-32. 4. Optimization Method The above dynamic optimization problem involves the performance index, which has to be optimized; the process model equations combined with initial state conditions; and different types of constraints. Without loss of generality, a typical formulation of a dynamic optimization problem can be stated as
u(t),p
νiJi ∑ i)1
J ) J0 +
∫tt F0(x,u,p,t) dt f
0
(34)
The other types of constraints to be taken into account typically include infinite-dimensional, interior-point and terminal-point equality and inequality constraints.18,19 Moreover, they can be written in the following canonical
(39)
where νi represents the multipliers that are determined at the optimum. The resulting global performance index is given by
min J ) G[x(tf),p] + u(t),p
∫tt F(x,u,p,t) dt f
0
(40)
5.2. Optimality Conditions. For the resulting dynamic optimization problem, eqs 40 and 35, the firstorder necessary conditions for an optimum can be derived using the calculus of variations.20,21 It consists of forming an augmented performance index J h by appending the ODE constraints (eqs 35) to the cost (eq 40) using the adjoint variables λ(t), defined by i)Nc
λ(t) ) λ0(t) +
νiλi(t) ∑ i)1
(41)
The development of variations of the augmented performance index is straightforward. By setting δJ h equal to zero, the first-order necessary conditions can be derived as
∂H ∂λ
x3 (t) )
subject to
x3 (t) ) f(x(t),u(t),p,t) with x(0) ) x0(p) (35)
νiFi ∑ i)1
F ) F0 +
and (iv) finally the heating/cooling rate of the reactor
min J0 ) G0[x(tf),p] +
(36)
5. Optimality Conditions
(31) (2) x(2) 5 (tf )
i
0
where ti e tf, i ) 1, 2, ..., Nc, and Nc is the number of constraints. For example, if an infinite-dimensional inequality constraint φ(x,u,p,t) e 0 is involved, then Gi ) 0 and Fi ) ω(max{φ(x,u,p,t), 0]}2, where ω is an adjustable weight. The constraint is then expressed as Ji e 0. In the case of a terminal-point equality constraint φ(x,u,p,tf) ) 0, Gi ) ω[φ(x,u,p,tf)]2 and Fi ) 0. The constraint is then expressed as Ji ) 0. The derivation of Gi and Fi for the other types of constraints is straightforward.18
)M h nf
or
∫tt Fi(x,u,p,t) dt
(42)
with
x(0) ) x0(p)
(43)
∂H ∂x
(44)
λ4 (t) ) -
7396
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004
with
|
(45)
∂H )0 ∂u
(46)
H(x,λ,u,p,t) ) F(x,u,p,t) + λTf(x,u,p,t)
(47)
λ(tf) )
∂G ∂x t)tf
where
Conditions 44, referred to as adjoint or costate equations, allow for the determination of the profiles of adjoint or costate variables λ(t). The corresponding boundary conditions are provided by conditions 45. It is important to notice that the set of conditions (eqs 4246) constitutes a two-point boundary value problem (TPBVP) because initial conditions 43 are specified for the state equations (eq 42), whereas for costate eq 44, only terminal conditions 45 are available. Because TPBVPs are known to be very difficult to solve numerically, appropriate computational methods are therefore needed.
Figure 3. Control variable approximation by piecewise-constant functions.
6.3. Computational Method Used. In this paper, the CVP method is used, and the gradients are evaluated using the adjoint equations. More specifically, the control vector u(t) is approximated by N piecewiseconstant variables over N intervals Ij of the time horizon [t0, tf] as (see Figure 3) N
u(t) )
6. Computational Methods The numerical methods used for the solution of dynamic optimization problems can be grouped into two categories: indirect and direct methods. 6.1. Indirect Methods. The objective of indirect methods is to solve the TPBVP (eqs 42-46), thus indirectly solving the dynamic optimization problem defined by eqs 40 and 35. The most well-known methods in this category are boundary condition iteration (BCI) and control vector iteration (CVI).22 However, these methods can exhibit numerical instabilities or slow convergence rates for many problems. 6.2. Direct Methods. In this category, there are two strategies: sequential and simultaneous. The sequential strategy, often called control vector parametrization (CVP), consists of approximating the control trajectory by a function of only a few parameters and leaving the state equations in the form of the original ODE system.18 Mostly, piecewise-constant control profiles are used. Hence, an infinite-dimensional optimization problem in continuous control variables is transformed into a finite-dimensional nonlinear programming (NLP) problem that can be solved by any gradient-based method [e.g., a sequential quadratic programming (SQP) method23]. The performance index evaluation is carried out by solving an initial value problem (IVP) of the original ODE system, and gradients of the performance index as well as the constraints with respect to the parameters (p and approximations of u) can be evaluated by solving either the adjoint equations (eqs 44 and 45)18 or the sensitivity equations.24 Moreover, this method is a feasible-type method, i.e., the solution is improved at each iteration. In the simultaneous strategy, both the control and state variables are discretized using polynomials (e.g., Lagrange polynomials) whose coefficients become the decision variables in a much larger NLP problem.25,26 Unlike the CVP method, the simultaneous method does not require the solution of IVPs at every iteration of the NLP. The method is, however, of infeasible type, i.e., the solution is available only once the iterative process has converged.
ujχj(t) ∑ j)1
(48)
with
χj(t) )
{
1 if t ∈ Ij 0 otherwise
The intervals Ij, which are not necessarily of the same length, are defined as Ij ) [tj-1, tj]; j ) 1, 2, ..., N. Their lengths are therefore additional parameters to be optimized. Finally, the optimization problem is set as a finitedimensional nonlinear programming problem N
min {J ) tf ) θ
∆ti} ∑ i)1
(49)
subject to the hybrid process model (eq 25) and constraints 26-32. The resulting vector of parameters θT ) [u1T, ..., uNT, ∆t1, ..., ∆tN, pT] is then of dimension (1 + Nu) × N + Np. The gradients of both the performance index and the constraints with respect to θij are derived from the variations of the augmented performance index27 as
∂Gi ∂Ji ) Hi(tN) + ∂tN ∂tN ∂Gi ∂Ji + ) Hi(t, j ) 1, ..., N - 1 j ) - Hi(tj ) + ∂tj ∂tj
(50) (51)
∂x0 ∂Ji ∂Gi ) - Jp,i(0) + λiT(t+ 0) ∂p ∂p ∂p
(52)
∂Ji ) Ju,i(tj-1) - Ju,i(tj), i ) 1, ..., Nc ∂uj
(53)
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7397
Figure 4. Algorithm of the CVP method. Table 1. Influence of the Number of Control Segments on the Performance Indexa
a
number of control segments
tf (s)
1 2 3 4 5 6 7 8 9 10
4605 4361 4027 4007 3966 3873 3848 3848 3848 3848
χf ) 0.7, M h nf ) 2 × 106 g‚mol-1, hcr ) -3 °C‚min-1.
The computational method developed here focuses on the accurate computation of gradients. It is noteworthy that two methods are currently used in the process optimization community: (i) the sensitivity method, which is implemented for example in gPROMS28 or in ABACUSS II,29 and (ii) the adjoint system method, which is implemented for example in DYNO.27 In addition to the adjoint system method developed within our optimization package DYNO, the derivatives of the process model equations and constraints with respect to the state vector/parameters are performed using available automatic differentiation packages (e.g., Adifor, Tapenade, ...). The resulting gradients are very accurate and greatly improve the convergence (if there is any) and its rate to the optimum.
where 7. Results and Discussion
∂Hi λ4 i ) ∂x
(54)
∂Hi ∂u
(55)
∂Hi i ) 1, ..., Nc ∂p
(56)
J4 u,i ) J4 p,i )
The CVP algorithm consists of the following steps: (i) estimation of the vector of parameters; (ii) forward integration of the state equations (eqs 42); (iii) backward integration of the adjoint equations (eqs 44); (iv) computation of gradients by means of eqs 50-53; (v) estimation of a new vector of parameters θ by means of an NLP solver; go to step ii if not converged, go to step vi otherwise; (vi) end. A schematic representation of this algorithm is given in Figure 4.
7.1. Influence of the Number of Control Segments. The first result concerns the influence of the number of control segments. Table 1 reports the variations of the resulting performance index for different values of the number of control segments and for fixed values of the final monomer conversion rate, the final number-average molecular weight, and the heating/ cooling rate. It can be seen that tf decreases with increasing number N of control segments. For values of N greater than 6, tf no longer changes. However, we found that the number of four control segments provides a good compromise between the accuracy of the computed optimal profiles and the computational speed in the implementation of on-line optimization. Four control segments are therefore used in the remainder of this work. 7.2. Influence of the Polymer Quality. Two variables can be used to characterize the polymer quality: the final number-average molecular weight, M h nf, or the final weight-average molecular weight, M h wf. Here, the
7398
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004
Figure 5. Control or decision variable (χf ) 0.7, hcr ) -3 °C‚min-1).
Figure 6. Reactor temperature obtained with time-varying jacket inlet temperature (χf ) 0.7, hcr ) -3 °C‚min-1).
Table 2. Influence of the Final Number-Average Molecular Weight on the Duration of the Control Segmentsa 10-6 M h nf (g‚mol-1)
∆t1 (s)
∆t2 (s)
∆t3 (s)
∆t4 (s)
1.0 2.0 3.0
653 675 621
615 1054 1796
603 1119 1871
613 1127 1338
a
χf ) 0.7, hcr ) -3 °C‚min-1.
influence of M h nf is analyzed. The influence of M h wf does not pose any specific problem and is not investigated here only for brevity. Figure 5 presents the variations of the control or decision variable, i.e., the jacket inlet temperature, versus the reduced time for fixed values of the final monomer conversion rate and the heating/cooling rate and for different values of M h nf. The corresponding batch periods are also reported in the figure. Table 2 lists the control segment durations in each case. In all cases, the optimal profile roughly consists of two temperature levels that depend on the value of M h nf. Moreover, it consists of operating at the high value for about 10 min in all three cases and at the low value for the rest of the batch period. This confirms the results obtained previously.1,2 On the other hand, the decrease in polymer quality requirement results in the decrease of temperature levels and, consequently, in longer batch period values. Finally, it is interesting to notice that Franc¸ ois et al.30 also obtained four arcs of temperature as the optimal solution for the same problem using the run-to-run optimization approach. However, in their physical interpretation of the solution, the optimal profiles are also reduced to two levels of temperature as in our case. Figure 6 presents the optimal profiles of the reactor temperature. It can be seen that the reactor temperature remains in the prescribed range and that no thermal runaway occurs. 7.3. Comparison of the Strategies. The operating constraints highly influence the value of the performance index and the shape of the optimal profiles. Thus, it is of primary importance to check their fulfilment. To this end, Figure 7 presents the constraint on the heating/cooling rate versus reduced time. This constraint, which is set to -3 °C‚min-1, is satisfied at all times and for all cases. For the rest of constraints, Table 3 reports that the fixed values of the final monomer conversion rate and the final number-average molecular weight are exactly the same as those obtained from the
Figure 7. Reactor heating/cooling rate (χf ) 0.7, hcr ) -3 °C‚min-1). Table 3. Influence of the Final Number-Average Molecular Weight on the Final Time with Computed Final Conversion Rate, Number- and Weight-Average Molecular Weights, and Polydispersity Indexa 10-6 M h nf (g‚mol-1)
tf (s)
Xf,c
10-6 M h nf,c (g‚mol-1)
h wf 10-6 M (g‚mol-1)
IPf
1.0 2.0 3.0
2483 3975 5627
0.70 0.70 0.70
1.0 2.0 3.0
2.11 5.63 8.82
2.11 2.82 2.94
a
χf ) 0.7, hcr ) -3 °C‚min-1.
calculation. The predicted values of the weight-average molecular weight and polydispersity index are also presented. In the complete discretization method used previously,2 the constraints were not perfectly fulfilled, hence leading to optimal profiles that exhibited the right overall variations but did not reach their final shape. This is mainly due to the fact that the optimization method used was not appropriate given that the computed optimal profiles exhibit a sort of stiff behavior. The consequence is that the optimal performance indices are larger than those obtained in the present work under the same operating conditions. This confirms the efficiency of the CVP method with respect to the complete discretization method for the optimization problem considered here. 7.4. Benefits of Time-Varying Optimal Control. To quantify the benefits realized upon the use of timevarying optimal control, a first comparison is performed between time-varying optimal control and the best uniform jacket inlet temperature (Table 4), and the second comparison is realized between time-varying optimal control and the best uniform reactor tempera-
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004 7399 Table 4. Influence of the Final Number-Average Molecular Weight on the Final Time, the Jacket Inlet Temperature, and the Polydispersity Index in the Case of the Best Uniform Controla 10-6 M h nf (g‚mol-1)
tf (s)
Tjin (K)
IPf
tf/tuniform f
1.0 2.0 3.0
2701 4605 6450
336 326 321
2.03 2.13 2.21
0.92 0.86 0.87
a
χf ) 0.7, hcr ) -3 °C‚min-1.
Table 5. Influence of the Final Number-Average Molecular Weight on the Final Time, the Inlet Jacket Temperature, and the Polydispersity Index in the Case of the Isothermal Reactora 10-6 M h nf (g‚mol-1)
tf (s)
Tjin (K)
IPf
tf/tisotherm f
1.0 2.0 3.0
2827 5588 8396
339 329 323
2.02 2.02 2.01
0.88 0.71 0.67
a
χf ) 0.7, hcr ) -3 °C‚min-1.
Figure 8. Reactor temperature obtained with the best uniform jacket inlet temperature (χf ) 0.7, hcr ) -3 °C‚min-1).
ture (Table 5). It is important to notice that both best uniform temperatures (jacket inlet and reactor) are deduced from the same optimization process where only a single control segment is used. It can be seen that up to 15% reduction in the batch period can be obtained when operating at the best uniform control. The corresponding reactor temperature profiles are presented in Figure 8 and show here also that the temperature lies in the range between the upper and lower limits with no thermal runaway. The reductions in batch period are even higher, i.e., up to 30%, when compared to the best uniform reactor temperature (Table 5). 8. Concluding Remarks The optimal profiles in open-loop control have been determined with great care in terms of the treatment and fulfilment of operating constraints. The results obtained show that significant benefits can be realized when operating at time-varying optimal profiles with respect to the best uniform control or the best reactor temperature. However, in the iterative optimization process described in the computational algorithm, we determined that, when we started with different initial estimates of the vector of parameters (a multistart method), in some cases, different values of the performance index were obtained and, consequently, different optimal control profiles. This means that the optimization problem exhibits multiple optima. On the other
hand, only first-order necessary conditions of optimality are used in this work, and the global optimality of the computed profiles is not guaranteed. The dynamic optimization problem should then be solved to global optimality prior to any experimental implementation. Therefore, the resulting profiles will provide an improvement of not only the performance index but also the optimal value within the specified constraints. In addition, global optimization methods have recently been developed that can be adapted to our problem.31 Note Added after ASAP Posting. This article was released ASAP on September 25, 2004, without all of the requested galley changes. The correct version was posted on October 13, 2004. Nomenclature A ) initiator concentration (mol‚L-1), cooling surface (m2) as ) surface area occupied by an emulsifier molecule (dm2) Cj ) cooling fluid heat capacity (J‚kg-1‚K-1) f ) initiator efficiency Fj ) cooling fluid flow rate (L‚s-1) fMS ) R-methylstyrene molar fraction in the initial load H ) Hamiltonian function hcr ) maximum heating/cooling rate of the reactor (K‚s-1) Ip ) polydispersity index J ) objective function kcm ) rate constant for initiator radical entry into micelles (L‚micelle-1‚s-1) kcp ) rate constant for initiator radical entry into particles (L‚part-1‚s-1) kd ) rate constant for initiator decomposition (s-1) kp ) rate constant for propagation (L‚mol-1‚s-1) ktrM ) rate constant for transfer to monomer (L‚mol-1‚s-1) L ) kinetic chain length (g‚mol-1) m ) number of micelles per unit volume (micelle‚L-1) M ) global monomer concentration (mol‚L-1) mrCp ) reactor total heat capacity (J‚K-1) MM ) monomer molecular weight (g‚mol-1) M h n ) number-average molecular weight (g‚mol-1) Mp ) monomer concentration in the particles (mol‚L-1) M h w ) weight-average molecular weight (g‚mol-1) n j ) average number of radicals per particle N ) number of inactive particles per unit volume (particle‚L-1) N • ) number of active particles per unit volume (particle‚L-1) Np ) total number of particles per unit volume (particle‚L-1) Na ) Avogadro’s number ns ) aggregation number of micelles P ) dead polymer concentration (mol‚L-1) p ) parameter Qi ) ith moment of the molecular weight distribution R • ) initiator radical concentration (mol‚L-1) Ra ) initiator decomposition rate (mol‚L-1‚s-1) Ri ) initiation rate (mol‚L-1‚s-1) Rn ) particle formation rate (mol‚L-1‚s-1) Rp ) polymerization rate (mol‚L-1‚s-1) Rt ) termination rate (mol‚L-1‚s-1) RtrM ) transfer to monomer rate (mol‚L-1‚s-1) S ) emulsifier concentration (mol‚L-1) t ) time (s) T ) reactor temperature (K) Tj ) jacket temperature (K) Tjin ) cooling fluid inlet temperature (K) u ) control variable U ) heat-transfer coefficient (J‚K-1‚s-1‚m-2) V ) reactor contents volume (L) Vj ) reactor jacket volume (L) x ) state vector
7400
Ind. Eng. Chem. Res., Vol. 43, No. 23, 2004
X ) conversion rate Xc ) critical conversion rate Greek Symbols ∆H ) polymerization reaction enthalpy (J‚mol-1) ) constant describing the efficiency of the particles relative to the micelles in collecting an initiator radical λ ) adjoint variable Fp ) polymer particle density (g‚L-1) Fj ) cooling fluid density (g‚L-1) FM ) monomer density (g‚L-1) FP ) polymer density (g‚L-1) θ ) parameter vector ωP ) polymer weight fraction in the particles Subscripts f ) final 0 ) initial
Literature Cited (1) Gentric, C. Optimisation Dynamique et Commande non Line´aire d’un Re´acteur de Polyme´risation en Emulsion. Ph.D. Thesis, INPL, Nancy, France, 1997. (2) Gentric, C.; Pla, F.; Latifi, M. A.; Corriou, J. P. Optimization and nonlinear control of a batch emulsion polymerization reactor. Chem. Eng. J. 1999, 75, 31-46. (3) O’Driscoll, K. F.; Dickson, J. R. Copolymerization with depropagation. II. Rate of copolymerization of styrene-R-methylstyrene. J. Macromol. Sci. A: Chem. 1968, 2, 449-457. (4) Branston, R. E.; Plaumann, H. P.; Sendorek, J. J. Emulsion copolymers of R-methylstyrene and styrene. J. Appl. Polym. Sci. 1990, 40, 1149-1162. (5) Johnston, H. K.; Rudin, A. Free radical copolymerization of R-methylstyrene and styrene. J. Paint. Technol. 1970, 43, 435448. (6) Rudin, A.; Chiang, S. S. Kinetics of free-radical copolymerization of R-methylstyrene and styrene. J. Polym. Sci. A: Polym. Chem. 1974, 12, 2235-2254. (7) Rudin, A.; Samanta, M. C. Emulsion copolymerization of R-methylstyrene and styrene. J. Appl. Polym. Sci. 1979, 24, 16651689. (8) Rudin, A.; Samanta, M. C. Emulsion copolymers of R-methylstyrene and styrene. J. Appl. Polym. Sci. 1979, 24, 1899-1916. (9) Rudin, A.; Samanta, M. C.; Van der Hoff, B. M. E. Monomer chain-transfer constants from emulsion copolymerization data: Styrene and R-methylstyrene. J. Polym. Sci. A: Polym. Chem. 1979, 17, 493-502. (10) Smith, W. V.; Ewart, R. H. Kinetics of emulsion polymerization. J. Chem. Phys. 1948, 16, 592-599. (11) Castellanos, J. R. Contribution a` la mode´lisation du proce´de´ de copolyme´risation en e´mulsion de l′R-methylstyre`ne et du styre`ne. Ph.D. Thesis, INPL, Nancy, France, 1996. (12) Asua, J. M.; Adams, M. E.; Sudol, E. D. A new approach for the estimation of kinetic parameters in emulsion polymerization systems. J. Appl. Polym. Sci. 1990, 39, 118-1213. (13) Gerrens, H. Uber den Geleffekt bei der Emulsions Polymerisation von Styrol. Z. Elektrochem. 1956, 60, 400-404. (14) Qin, J.; Li, H.; Zhang, Z. Modeling of high-conversion binary copolymerization. Polymer 2003, 44, 2599-2604.
(15) Harada, M.; Nomura, M.; Kojima, H.; Eguchi, W.; Nagata, S. Rate of emulsion polymerization of styrene. J. Appl. Polym. Sci. 1972, 16, 811-833. (16) Villermaux, J.; Blavier, L. Free radical polymerization engineeringsI A new method for modeling free radical homogeneous polymerization reactions. Chem. Eng. Sci. 1984, 39, 8799. (17) Galan, S.; Feehry, W. F.; Barton, P. I. Parametric sensitivity function for hybrid discrete/continuous systems. Appl. Numer. Math. 1999, 31, 17-47. (18) Goh, C. J.; Teo, K. L. Control parametrization: A unified approach to optimal control problems with general constraints. Automatica 1988, 24, 3-18. (19) Teo, K. L.; Goh, C. J.; Wong, K. H. A Unified Computational Approach to Optimal Control Problems; John Wiley and Sons: New York, 1991. (20) Pontryagin, L. S.; Boltyanskii, V. G.; Gamkrelidze, R. V.; Mishchenko, E. F. The Mathematical Theory of Optimal Processes; Pergamon Press: New York, 1964. (21) Bryson, A. E.; Ho, Y. C. Applied Optimal Control; Hemisphere Publishing Corporation: New York, 1975. (22) Ray, W. H.; Szekely, J. Process Optimization with Applications in Metallurgy and Chemical Engineering; John Wiley and Sons: New York, 1973. (23) Schittkowski, K. NLPQL: A FORTRAN subroutine solving constrained nonlinear programming problems. Ann. Oper. Res. 1985, 5, 485-500. (24) Feehery, W. F. Dynamic optimization with path constraints. Ph.D. Thesis, MIT, Cambridge, MA, 1998. (25) Biegler, L. T. Solution of dynamic optimization problems by successive quadratic problems and orthogonal collocation. Comput. Chem. Eng. 1984, 8, 243-248. (26) Cuthrell, J. E.; Biegler, L. T. On the optimization of differential-algebraic process systems. AIChE J. 1987, 33, 12571270. (27) Fikar, M.; Latifi, M. A. User’s Guide for FORTRAN Dynamic Optimisation Code DYNO; Technical Report: Laboratoire des Sciences du Ge´nie Chimique, CNRS: Nancy, France, 2002. (28) Barton, P. I.; Pantelides, C. C. gPROMS: A combined discrete/continuous environment for chemical processing systems. Simul. Ser. 1993, 25, 25-34. (29) Tolsma, J. E.; Clabaugh, J.; Barton, P. I. ABACUSS II: Advanced Modeling Environment and Embedded Simulator; MIT: Cambridge, MA, 1999 (see http://yoric.mit.edu/abacuss2/ abacuss2.html). (30) Franc¸ ois, G.; Srinivasan, B.; Bonvin, D. Run-to-run optimization of batch emulsion polymerization. In Proceedings of the IFAC World Congress; Elsevier Science Ltd.: Amsterdam, 2002; pp 1258-1263. (31) Chachuat, B.; Latifi, M. A. A new approach in deterministic global optimisation of problems with ordinary differential equations. In Nonconvex Optimization and Its Applications; Floudas, C. A., Pardalos, P. M., Eds.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 2004; pp 83-108.
Received for review February 20, 2004 Revised manuscript received September 8, 2004 Accepted September 13, 2004 IE0498549
