92

Ind. Eng. Chem. Res. 2008, 47, 92-104

Integrated Multilevel Optimization in Large-Scale Poly(Ethylene Terephthalate) Plants Flavio Manenti* and Maurizio Rovaglio CMIC Department “Giulio Natta”, Politecnico di Milano, Piazza Leonardo Da Vinci 32, I-20133 Milano, Italy

The present paper discusses the vertical integration of advanced predictive control, businesswide optimization and markets laws, based on a large-scale first-principles model. The case study is the poly(ethylene terephthalate) (PET) plant, with esterifier sections, finisher reactors, a crystallizer, and a solid-state polymerizer, and the overall mathematical model consists of more than 1500 differential and algebraic equations (DAEs). The complex integration among supply chain layers is reached using a Windows-Unix hybrid system; the methodology adopted is described, and a real-time strategy is implemented by using a secure shell connection for linking a PentiumIV, which simulates the plant, and a cluster Opteron, where the optimizer operates. 1. Introduction The life-motive to study a new optimization tool is the need of finding an efficient solution for the supply chain management problem. There is the need of developing suitable integrated tools, which allow the linking between a first process optimization layer, the model predictive control (MPC) level, and an upper level, the businesswide, real-time dynamic optimization (RTDO), so to solve the global problem in a reasonable computational time. In this perspective, the scope of the present work is to investigate consequences of the multilevel dynamic integration in the process control hierarchy. The present industrial state of art in process control is represented by the linear MPC,1 for which several practical applications are known. Instead, the plant optimization is still realized by steady models,2 even if, during this last period, some dynamic optimal application has been implemented in the industrial world,3 demonstrating the profit increase in the order of $12.5 million/year and 1 month payback. However, in order to improve process flexibility and efficiency, nonlinear MPC and dynamic optimization tools, both based on the same mathematical model, seem to be appealing for a global simultaneous solution.4 The model developed in this manuscript regards a poly(ethylene terephthalate) (PET) synthesis plant. In fact, the polymer industry is particularly concerned in changes of product amount, which are very frequent because of market dynamics and demand portfolio, and, further, the same polymer can be required with different specifications in physical properties and quality (grade variation of the end product). Moreover, the process management of such a plant is hard to perform because of the complex and nonlinear behaviors and, generally, to poor (traditional) control schemes. Numerous studies on MPC and businesswide RTDO are available in the recent polymers literature,5-9 and moreover, previous research in the PET process optimization is already developed.10,11 Here, real-time methodologies for achieving the dynamic optimization will be applied on an industrial largescale system, and its mathematical modeling and numerical integration will be described in Section 3. Section 4 will show the plantwide control system of the PET plant and the discussion about the online implementation of the integrated multilevel * Corresponding author. Phone: +39-02-2399-3287. Fax: +39-027063-8173. E-mail: [email protected].

optimization. Numerical methods and solvers adopted will be discuss in Section 5, together with results for a typical gradechange problem and the optimized case. Finally, the summary and conclusions will be presented in Section 6. Nomenclature will be explained at the end of the paper. 2. Closed-Loop Technique Particularly in the polymer industry, continuous processes assume a strongly nonlinear behavior, and they are frequently subject to large changes in their operating conditions, dictated by market terms. Sales and purchasing conditions motivate grade transitions and the need of minimizing the off-spec product during each transient and, at the same time, of satisfying other process constraints. There is an evident need of integration techniques based on nonlinear models, able to better predict, control, and optimize unsteady conditions.12 The closed-loop technique could be a solution: it allows the integration between the predictive control, the dynamic optimization, and the market laws, by evaluating the best values for manipulated variables (MVs) and, at the same time, by the economic estimation of optimal set points, with respect to actual operating conditions of the plant. To obtain it, the former problem of businesswide RTDO, with a typical time-scale in the order of hours/days, is lowered to the MPC time-scale, in the order of minutes, and both setpoints and MVs are continually evaluated and implemented in the real plant, by using a realtime strategy tool, an efficient numerical library for integrating model equations, and a simultaneous optimizer. The methodology employed follows ideas discussed in the literature about path-constrained problems.13-15 The real-time strategy consists of a continuous interaction between the plant site, or its mathematical model, and the optimizer. Generally, it can be described through three different steps: (i) the reading of the actual operating conditions of the process (snapshot of current operating conditions and storage of the process data); (ii) the resolution of the nonlinear programming (NLP) problem, by adopting an optimization tool (usually a sequential quadratic programming, SQP); and (iii) the implementation of optimal values in the real plant, until the next sampling time. The linking between the plant and the optimizer must be guaranteed by an object linking and embedding (OLE) for process control (OPC) or a secure shell for file transfer (SSH), respectively, if a plant or its mathematical model are available. By continually repeating these steps every sampling time, it is

10.1021/ie070355s CCC: $40.75 © 2008 American Chemical Society Published on Web 11/30/2007

Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 93 Table 1. Molecular Structure of Kinetic Scheme Components

possible to speak about real-time strategy, if the selected time span is in the order of minutes or, at most, a few hours. 3. Large-Scale Mathematical Model Only continuous industrial processes that do not present a strongly nonlinear behavior are usually good candidates for the implementation of linear families of algorithms. This kind of linear advanced control does not require an optimization for their solution; in fact, they can be solved through analytical procedures. An optimizer becomes necessary when the process presents a nonlinear nature, which is much more difficult to solve in an analytical way. In this instance, the main reasons for the selection of the PET are the typical nonlinear behavior of polymer plants and the need of operating under frequent grade transitions, with traditional and inadequate control systems. A lot of studies can be quoted in this field during the last few years,8,9,16 and each of them applies the dynamic optimization for overcoming grade-change transition problems in polymeric processes. 3.1. PET Process Configuration. The PET can be synthesized through two different processes: (1) Dimethylterephthalate (DMT): The reactants are the ethylene glycol (EG) and the DMT. Some problematic issues of DMT processes concern the availability of terephthalic acid (TPA) with a high degree of purity and its utilization: in fact, the DMT is solid at normal operating conditions and it could lead to higher operating temperatures. (2) Pure terephthalic acid (PTA): It is the most widespread. EG and TPA are fed to the plant. PTA processes allow better yields and faster reaction rates at the same operating conditions of DMT processes. Moreover, the absence of methanol and the 20% reduction in variable costs render the PTA technology the almost sole solution for new PET plant realization. The DMT technology can be considered overcome, especially by considering actual plant production capacity17 around 200500 t/d, with a final intrinsic viscosity (IV) in the range of 0.500.70 dL/g and an uncertainty of 0.005 dL/g. The PTA configuration will be considered in the present study. 3.2. PTA Process. Practical problems are the solid state of TPA at ambient conditions and its low solubility in EG, even at high temperatures. Moreover, initial esterification and polycondensation processes are controlled by mixing, so the

Figure 1. PE section consists of the initial mixer, the primary esterifier, and a distillation column for recycling EG. Table 2. Reaction Scheme and Rate Constants19 rate constants no.

reactions

forward

reverse

1 2 3 4 5 6 7 8 9 10 11

EG + TPA T tEG ≈ tTPA + W EG + tTPA T tEG ≈ bTPA∼ + W ∼tEG + TPA T ∼bEG ≈ tTPA + W ∼tEG + ∼tTPA T ∼bEG ≈ bTPA∼ + W ∼tEG + ∼tEG T ∼bEG∼ + EG ∼bTPA ≈ tEG + ∼tEG f ∼tDEG + ∼tTPA ∼bTPA ≈ tEG f ∼tTPA + AA ∼bTPA ≈ tEG + EG f ∼tTPA + DEG ∼tEG + ∼tDEG T ∼bEG∼ + DEG ∼bTPA ≈ bEG∼ f ∼tVIN + ∼tTPA ∼tVIN + ∼tEG f ∼bEG∼ + AA

k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11

k1/K1 k2/K2 k3/K3 k4/K4 k5/K5

k9/K9

thermodynamics and the kinetics of chemical reactions will be limited, even if clear data is not available in literature. A PET production train usually consists of six sections:10 (1) primary esterification ) PE; (2) secondary esterification ) SE; (3) low polymerization ) LP; (4) intermediate polymerization ) IP; (5) high polymerization ) HP; and (6) solid-state polymerization ) SSP. 3.2.1. Primary Esterifier. The first reactor is the primary esterifier, which is directly fed by raw materials, EG and TPA, previously mixed in metered amounts in a dedicated unit; a possible scheme is reported in Figure 1. The pressure is usually 2-8 atm, and the temperature is 530-555 K. The relatively high temperature allows a higher solubility of TPA inside the oligomer; nevertheless, high pressures are needed in order to limit the EG evaporation. The model of the first continuous stirred tank reactor (CSTR) has been constructed by following a consolidated kinetic scheme.18,19 With a two-phase reactor, flash equilibrium relations are needed because the vapor flow rate has to be determined for rendering solvable the differential and algebraic equations (DAEs) system. From Tables 1-4, the segment approach useful for representing the polymer/solvent properties and polymerization kinetic scheme is shown. In these tables, Fmix is useful for evaluating the concentrations of chemical compounds involved in recycles. Fmix is useless in

94 Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 Table 3. Terms of Reaction and Generation Rates reaction rates

generation rates

R1 ) 2k1C1C2 - (k1/K1)C4C7 R2 ) 2k2C1C4 - (k2/K2)C6C7 R3 ) 2k3C2C3 - (k3/K3)C4C7

G1 ) -R1 - R2 + R5 - R8 G2 ) -R1 - R3 G3 ) R1 + R2 - R3 - R4 - 2R5 2R6 - R7 - R8 - R9 - R11 G 4 ) R1 - R2 + R 3 - R4 + R7 + R8 +R10 G5 ) R3 + R4 + R5 + R9 - R10 + R11 G6 ) R2 + R4 - R10 G 7 ) R1 + R2 + R 3 + R4 G 8 ) R8 + R9 G9 ) R7 + R11 G10 ) R6 + R8 - R9 G11 ) R10 - R11

R4 ) 2k4C3C4 - (k4/K4)C6C7 R5 ) k5C3C3 - 4(k5/K5)C1C5 R6 ) k6C3C3 R7 ) k7C3 R8 ) k8C1C3 R9 ) k9C3C11 - 4(k9/K9)C5C9 R10 ) k10C5C6/(C3 + C5) R11 ) k11C3C12

balance in the reactor. The relative algebraic or differential variable is reported in the last column of the table, for each equation. The vapor-liquid equilibrium is calculated with the Flory-Huggins equations:

PEG ) γEGPpEGxEG PW ) γWPpWxW PDEG ) γDEGPpDEGxDEG PAA ) γAAppAAxAA ln

Table 4. Kinetic ki ) Ai exp(-Ei/(RgasT)), by Supposing a Catalyst Concentration (Sb2O3) in the Order of 0.05 wt % Data;19

reaction no.

Ai

Ei (kcal/kmol)

Ki

1 2 3 4 5 6 7 8 9 10 11

2.08 × 106 (L mol-1 min-1) 2.08 × 106 (L mol-1 min-1) 2.08 × 106 (L mol-1 min-1) 2.08 × 106 (L mol-1 min-1) 1.36 × 106 (L mol-1 min-1) 8.32 × 107 (L mol-1 min-1) 8.32 × 107 (min-1) 8.32 × 107 (L mol-1 min-1) 1.36 × 106 (L mol-1 min-1) 7.2 × 109 (min-1) 1.36 × 106 (L mol-1 min-1)

17.6 17.6 17.6 17.6 18.5 29.8 29.8 29.8 18.5 37.8 18.5

2.5 2.5 1.25 1.25 0.5

Pj Ppj

(

) ln Φpj + 1 -

Cj ) Cout j

Fout Fout out (j ) (1,3)/11); C ) C 2 2 Fout - Fsol Fout - Fsol

all the other cases. It is possible to evaluate the main properties of the polymer, such as the IV and the average molecular weight MWn as follows:

IV ) 2.1 × 10-4(192.17 × MWn)0.82

(2)

Finally, the model equations are summarized in Table 5, where 12 ordinary differential equations (ODEs), related to mass and energy balances, are listed, together with 3 algebraic equations, which define the pressure equilibrium, the overall EG mass balance in the distillation column, and the overall mass

)

1 Φp + χpj (Φpj )2 X hN j

The existence of undissolved TPA is taken into account by calculating the solubility into EG/oligomer solution.20 The TPA concentration in the liquid solution (Cliq,out ) RFout) and the liq,out) are easily retrievsolid TPA residue (Csol,out ) Cout 2 - C able. Concentration values used in the kinetic scheme are obviously referred to the liquid volume only, and hence, the following corrections are needed:

1.0

MWn ) CtEG + CtTPA + CbEG + CbTPA + CbDEG + CtVIN + CtDEG CtEG + CtTPA + CtVIN + CtDEG (1)

(3)

Fsol ) Fout

Csol,outMW2 FTPA

(4)

3.2.2. Secondary Esterifier. The oligomer produced in the PE section is fed to the secondary esterifier that operates close to atmospheric conditions (pressure ∼1-2 atm and temperatures in the range of 535-575 K). To enhance effectiveness, this unit is often divided into several equal-sized chambers, as reported in Figure 2, and each one is modeled as a CSTR. In order to model this unit, the same equations of the PE are adopted under the assumption that each chamber is independent from the others and no back-mixing is considered. The backmixing has the capability to lower the conversion of functional end groups and, hence, to change the polymer molecular weight and properties. Nevertheless, this assumption leads to negligible errors thanks to small gas holdups and operating regimes far away from the mass transfer limitation.10 Subsequently, the SE

Table 5. Balance Equations of the CSTR in the PE Section; The DAE System no. 1 2 3 4 5 6 7 8 9 10 11 12

13 14 15 16

equation out out out dt ) 1 / FmixV (FinFinCin 1 - F F C1 in in out out / dt ) 1 V (F C2 - F C2 ) + G2(t) out out dt ) 1 / V (FinCin 3 - F C3 ) + G3(t) in in dt ) 1 / V (F C4 - FoutCout 4 ) + G4(t) out out dt ) 1 / V (FinCin 5 - F C5 ) + G5(t) out out dt ) 1 / V (FinCin 6 - F C6 ) + G6(t) in in in dt ) 1 / FmixV (F F C7 - FoutFoutCout 7 out out out dt ) 1 / FmixV (FinFinCin 8 - F F C8 out out out dt ) 1 / FmixV (FinFinCin 9 - F F C9 out out dt ) 1 / V (FinCin 10 - F C10 ) + G10(t) out out dt ) 1 / V (FinCin 11 - F C11 ) + G11(t) in in / dT dt ) 1 / FmixCpV [FmixCpF (T - Tref) NCvap Cvap(∆Hev + FmixCpFout(T - Tref) - Fvap ∑j)1 j j rec FEGCpF (T - Tref) - UArea(T - Th)] p p p p Ptot ) P + P + P +P

dCout 1 dCout 2 dCout 3 dCout 4 dCout 5 dCout 6 dCout 7 dCout 8 dCout 9 dCout 10 dCout 11

/ / / / / / / / / / /

EG

F

vap

W

DEG

AA

) CHV‚Y‚P xP - P / TPtot; tot

tot

d

rec Frec ) Fvap(Cvap 1 )/C1 Fout ) (FinFin + FEGFrec - FvapFvap)/Fout

variable

FvapFvapCvap 1

+

FrecFrecCrec 1 )

+ G1(t)

FvapFvapCvap + FrecFrecCrec 7 7 ) + G7(t) FvapFvapCvap + FrecFrecCrec 8 8 ) + G8(t) FvapFvapCvap + FrecFrecCrec 9 9 ) + G9(t)

Cout 1 Cout 2 Cout 3 Cout 4 Cout 5 Cout 6 Cout 7 Cout 8 Cout 8 Cout 10 Cout 11 T

Cpgas j (T - Tref)) +

Y ) 1 - 2 3 xP - P / TPtot /

tot

d

Ptot Fvap Frec Fout

Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 95

Figure 2. Multichambers SE section and the associated distillation column for the EG recycle.

those of previous stages, adding a single constraint: no recycle stream exists in this unit. 3.2.4. Intermediate and High Polymerization. The polycondensation stage of the reactor train consists of IP and HP, the so-called finishers, where the melt-phase polymerization reaches its maximum feasible degree. In fact, these units are characterized by a marked influence of the mass transfer, and they significantly differ from previous stages. The most common configuration adopts disks-ring reactors (Figure 4) in which a number of disks, attached to a rotating shaft, push toward the melt polymer by enhancing the mass transfer through the generation of a thin surface film. The evaporation of condensed byproducts, EG and W, is also promoted by lower vacuum pressures, 1.5-4 mmHg in the IP section and 0.5-2 mmHg in the HP section, in order to favor the lengthening of polymer chains. To overcome the high polymer viscosity, a high power to the agitator and a high temperature, ∼555-575 K, are required, but this also leads to an undesired degradation due to side reactions, especially through the formation of AA and DEG byproducts. Therefore, the final IV obtainable with these units is limited to 0.5-0.7 dL/g, and higher qualities are denied here, so to avoid higher residual contents in secondary specs (DEG and AA). However, higher IV can be obtained by continuing the polymerization in the SSP section that, directly operating on the solid polymer at lower temperatures, can overcome constraints about byproducts residuals, and the bottle grade can be easily reached. The typical textile fiber grade is ∼0.55-0.65 dL/g, whereas the bottle grade requires an IV ≈ 0.80 dL/g. With the SSP section, this may be easily obtained, and the IV can reach >1.0 dL/g, by satisfying also the tire-cord resins market (very high IV). In any case, for modeling a disks-ring reactor, the proposed path21,22 is to model the finisher as a two-phase plug-flow reactor. The mass transfer is kept into account by using a coefficient kja that includes the real coefficient of mass transfer kj and the specific interphase area a. This unique parameter allows one to easily fit the model results to specific equipment behavior, and it must be emphasized that a marked sensitivity of reactor performances to this constant is attended. Therefore, a generic balance equation will have the following form,

∂y F ∂y )+ f(t) ∂t Area ∂z

(5)

where Area is the reactor section and f(t) is a well-defined function. It is a partial differential equation (PDE) system, and in order to solve this model coupled with previous parts of the overall DAE system, it is convenient to rewrite it in finite differences terms,

Figure 3. LP section. The EG recovered by the spray condenser is recycled to the beginning of the process.

∂y yn-1 - yn ≈ ∂z ∆z

(6)

that once substituted gives section will be modeled by adopting the same DAE system listed in Table 5, which is repeated three times, one for every chamber. 3.2.3. Low Polymerizer. A first polymerizer is placed in the LP section, just after the esterification train, as reported in Figure 3. In this kind of reactor, which is typically a CSTR operating at 555-580 K and at a medium vacuum pressure (50-500 mmHg), the exceeding EG and W produced are stripped off and the influence of the mass transfer in the reaction equilibrium gradually grows. The model of this unit is still analogous to

∂y F yn-1 - yn )+ f(t) ∂t Area ∆z

(7)

and finally, substituting a fictitious volume of the n finite elements V ) Area‚∆z, the following form is obtained:

∂y F ) - (yn-1 - yn) + f(t) ∂t V

(8)

96 Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008

Figure 4. Disks-ring reactors for IP and HP sections. The recovered EG is recycled to the LP section. Table 6. Balance Equations for Finisher Sections; DAE System no.

equation

) 1 V ) 1/V ) 1/V ) 1/V ) 1/V ) 1/V ) 1/V ) 1/V ) 1/V ) 1/V ) 1/V dT / dt ) 1 / FmixCpV FmixCpF (T dCout 1 dCout 2 dCout 3 dCout 4 dCout 5 dCout 6 dCout 7 dCout 8 dCout 9 dCout 10 dCout 11

1 2 3 4 5 6 7 8 9 10 11 12

/ dt / dt / dt / dt / dt / dt / dt / dt / dt / dt / dt

/

(FinCin 1 (FinCin 2 (FinCin 3 (FinCin 4 (FinCin 5 (FinCin 6 in in (F C7 (FinCin 8 (FinCin 9 (FinCin 10 (FinCin 11

variable

+ G1(t) + G2(t) + G3(t) + G4(t) + G5(t) + G6(t) + G7(t) - kja(Cout 7 + G8(t) - kja(Cout 8 + G7(t) - kja(Cout 9 + G10(t) + G11(t) - Tref) - FmixCpFout(T

FoutCout 1 ) FoutCout 2 ) FoutCout 3 ) FoutCout 4 ) FoutCout 5 ) FoutCout 6 ) FoutCout 7 ) FoutCout 8 ) FoutCout 9 ) FoutCout 10 ) FoutCout 11 ) in in

kja(Cout 1

-

Cout 1 Cout 2 Cout 3 Cout 4 Cout 5 Cout 6 Cout 7 Cout 8 Cout 9 Cout 10 Cout 11 T

- C/7) - C/8) - C/9) - Tref) -

V kja(cout - c/j )(∆HEV + Cpgas j j j (T - Tref)) - UArea(T out / NC k a(Cout - C/) yEG ) kja(CEG - CEG) ∑j)1 j j j / NC k a(Cout - C/) yW ) kja(Cout C ) ∑ j)1 W W j j j / NC k a(Cout - C/) yDEG ) kja(Cout ∑j)1 DEG - CDEG) j j j NC k a(Cout - C/)MW ) Fout ) 1 / Fout (FinFin - V ∑j)1 j j j j NC ∑j)1

C/1)

- Th)

/

13

yEG

/

14

yW

/

15 16

Practically, the plug-flow model has been translated into a series of CSTRs, and the equation can be included in the overall DAE system. In this instance, a series constituted by 10 CSTRs has been employed for both finishers. Discretized mass and energy balance equations are reported in Table 6. The first three algebraic equations are simplified estimations of gas-phase molar fraction that are necessary to evaluate the interphase mass transfer. Interfacial equilibrium concentrations C/j of the volatile compound j are calculated with the Flory-Huggins equations through the following,

C/j )

( ) CPET NC

1-

x/j ∑ j)1

x/j ; x/j )

Pyj Ppj γj

(9)

yDEG Fout

and the value of molar fraction related to the overall existing mass transfer with the following:

yj )

/ kja(Cout i - Ci )

(10)

NC

∑ j)1

kja(Cout j

-

C/j )

This equation hides the hypothesis that no back-mixing between blocks exist. This leads to negligible errors with gas holdups being small and operating regimes being far away from the mass transfer limitation, as already discussed in the literature.10 Through the back-mixing exclusion, a notable advantage in the tailored structure of the problem is obtained. The molar fraction of the component AA is retrievable, trivially imposing the normalization (yAA ) 1 - yEG - yW - yDEG). The last algebraic equation needed is the overall mass balance:

Fout )

1 F

out

(

)

NC

FinFin - V

/ kja(Cout ∑ j - Cj )MWj j)1

Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 97 Table 9. Reaction Rate Rj(t) Terms Referred to Kang’s Kinetic Scheme20

(11)

3.2.5. Solid-State Polymerizer. The last section of the reactor train is the SSP (see Figure 5). In this unit, the solid polymer is extruded in pellets and flows into a moving-bed reactor that is a cylindrical unit terminating with a hopper discharge with a countercurrent flow of inert gas (usually 0.5-4 (kg of N2)/(kg of PET)). The operating temperature is in the range of 475505 K, far away from the melting temperature (530-536 K) and also from the sticking temperature (518 K). The proposed model for the SSP is substantially analogous to the Nylon-6,6 SSP reactor model,20,23 in order to describe the degradation and polycondensation reactions; our model adopts a reduced kinetic model with 10 species and 9 reactions,11,24 where the proposed system of PDE is able to account for both the diffusion inside the polymer particle, considered as an additional mass transfer resistance, and the evolution along the axial coordinate z of the reactor, in the gas-phase and the solid-phase concentrations. Tables 7-10 detail Kang’s kinetic scheme for the SSP. In the SSP unit, another important piece of equipment is the preheater and crystallizer (PHCR). In industrial practice, they

Figure 5. Pelletizer in the SSP section. Table 7. Kang’s Kinetic Scheme20 for the SSP no.

reaction

forward rate constant

reverse rate constant

1 2 3 4 5 6 7 8 9

EG + TPA T tEG + tTPA + W EG + tTPA T tEG + bTPA + W tEG + TPA T bEG + tTPA + W tEG + tTPA T bEG + bTPA + W tEG + tEG T bEG + EG tEG + tEG f bDEG + W bEG + bTPA f tVIN + tTPA tEG + bTPA f AA + tTPA tEG + tVIN f bEG + AA

k1 k2 k3 k4 k5 k6 k7 k8 k9

k1/K1 k2/K2 k3/K3 k4/K4 k5/K5

Table 8. Kinetic Parameters of Kang’s Scheme20 Ai A1 ) A2 ) 2A3 A3 ) A4 ) 6.8 × 1012 (L/mol)/min A5 ) 5.4 × 1012 (L/mol)/min A6 ) 1.8 × 1015 (L/mol)/min A7 ) 3.6 × 109 (1/min) A8 ) 2.3 × 109 (1/min) A9 ) A 5

Ei (kcal/kmol)

Ki

E1 ) E 2 ) E 3 K1 ) K2 ) 2.5 E3 ) E4 ) 17.6 (kcal/mol) K3 ) K4 ) 1.25 E5 ) 18.5 (kcal/mol) E6 ) 29.8 (kcal/mol) E7 ) 37.8 (kcal/mol) E 8 ) E6 E9 ) E 5

K5 ) 0.50

R1 ) 4k1CEGCTPA - (k1/K1)CWCtTPA(CtEG)/CtEG + CbEG R2 ) 2k2CEGCtTPA - 2(k2/K2)CWCbTPA(CtEG)/CtEG + CbEG R3 ) 2k3CtEGCTPA - (k3/K3)CWCtTPA(CbEG)/CtEG + CbEG R4 ) k4CtEGCtTPA - 2(k4/K4)CWCbTPA(CbEG)/CtEG + CbEG R5 ) k5CtEGCtEG - 4(k5/K5)CEGCbEG R6 ) k6CtEGCtEG R7 ) k7CbEG(CbTPA)/CtEG + CbEG R8 ) k8CtEG(CbTPA)/CtEG + CbEG R9 ) k9CtEGCtVIN Table 10. Generation Rate Gj(t) Terms Referred to Kang’s Kinetic Scheme20 GEG(t) ) -R1 - R2 + R5 GTPA(t) ) -R1 - R3 GtEG(t) ) R1 + R2 - R3 - R4 - 2R5 - 2R6 - R8 - R9 GtTPA(t) ) R1 - R2 + R3 - R4 + R7 + R8 GbEG(t) ) R3 + R4 + R5 - R7 + R9 GbTPA(t) ) R2 + R4 - R7 - R9 GtEG(t) ) R1 + R2 + R3 + R4 + R6 GbDEG(t) ) R6 GtVIN(t) ) R7 - R9 GAA(t) ) R8 + R9 Table 11. CVs and MVs in the PE Section CVs (or bounded variables)

min

max

reactor temperature (K) reactor pressure (atm)

515 1.0

565 8.0

MVs

min

max

EG feed flow (L/min) EG-TPA flow ratio wall temperature (K) discharge pressure (atm)

2.0 1.0 500 3.0

12.0 1.2 580 8.0

exist both as separated units or as a unique apparatus. Characteristics of such a unit are not so clear: bubble beds or agitated vessels are the most common configurations, but further details are hidden by patents of plant vendors. In this study, the SSP unit is modeled as a plug-flow reactor where the flow of solid particles raises the temperature in the range of 475 K. During the heating, another phenomenon occurs: the crystallization of the polymer that is initially amorphous and during the process reaches its natural semicrystalline state. Equations of the detailed mathematical model have been already reported in a previous work. 24 At the moment, no probing confront is possible because of the lack of experimental data.19 Only some literature confront10 and the general know-how about the process normal conditions allow us to sustain a good agreement of the overall model with the real behavior. 4. Formulation Of The Multilevel Dynamic Optimization 4.1. Plantwide Control in PET Plants. Besides 1500+ state variables, the control problem consists of 15 MVs, 2 controlled variables (CVs), and 16 upper and lower constrained variables. Loops used in the IV control are shown in Figure 6: an inferential, based on the dissipation energy of the electric motor that moves the disks of the reactor, is adopted. In the short range of interest of the present paper, a piecewise linear correspondence between the IV of the polymer and the electric consumption can be considered. Figure 7 reports controls implemented in the SSP section: a flow-rate controller manages the feed of warm nitrogen in the bottom of the pelletizer, whereas a temperature controller manages the utility, positioned on the nitrogen recycle stream. In the other sections, the reactor temperature is regulated by wall temperatures and pressures by

98 Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 Table 13. CVs and MVs in the LP Section CVs (or bounded variables)

min

max

reactor temperature (K) reactor pressure (atm)

525 0.05

585 1.5

MVs

min

max

wall temperature (K) discharge pressure (atm)

520 0.1

600 1.0

Table 14. CVs and MVs in the IP Section CVs (or bounded variables)

min

max

reactor temperature (K) reactor pressure (Torr)

525 1

585 4

MVs

min

max

wall temperature (K) discharge pressure (atm)

520 2.0

600 6.0

Table 15. CVs and MVs in the HP Section

Figure 6. IP control scheme: an inferential (energy consumption of the electric motor) is adopted for estimating the polymer IV.

CVs (or bounded variables)

min

max

reactor temperature (K) reactor pressure (Torr) AA concentration (mol/L) DEG concentration (mol/L) intrinsic viscosity (dL/g)

525 0.5 0.0 0.0 unconstrained

585 2.0 1 × 10-4 0.1

MVs

min

max

wall temperature (K) discharge pressure (Torr)

520 0.5

600 2.5

Table 16. CVs and MVs in the SSP Section CVs (or bounded variables)

min

max

reactor temperature (K) reactor pressure (Torr) AA concentration (mol/L) DEG concentration (mol/L) intrinsic viscosity (dL/g)

453 0.5 0.0 0.0 unconstrained

514 2.0 1 × 10-4 0.1

MVs

min

max

inlet gas temperature (K) gas-solid flow ratio polymer mass flow (kg/min)

453 1.0 0.5

514 4.5 2.0

4.2. Integrated Objective Function. The solution of a multilevel optimization problem requires a single (integrated) objective function: it has to include MPC tasks and businesswide RTDO objectives. The MPC has the aim to control the plant, by keeping CVs as close as possible to their setpoints. So, the MPC will evaluate the best set of MVs that drives the plant along optimal trajectories. According to Tables 11-16, the MPC objective function could assume the following form, Figure 7. SSP control scheme: a flow controller manages the inlet nitrogen and a temperature controller regulates the utility open command of the nitrogen recycle valve.

ΦMPC ) min { u∈U

∫0t [ω1(IVHP(t) - IVHPset)2 + ω1 f

MV

(IVSSP(t) - IVSSPset)2 + Table 12. CVs and MVs in the SE Section

2 ||MVn,(t) - MVtar ∑ n ||ω ] dτ} n)1 2

(12)

CVs (or bounded variables)

min

max

reactor temperature (K) reactor pressure (atm)

525 1.0

585 3.0

MVs

min

max

wall temperature (K) discharge pressure (atm)

520 0.5

600 2.0

the opening of discharge valves. In Tables 11-16, CVs and MVs for each session are reported. PE and SE distillation columns and LP, IP, and HP spray condensers are not taken into consideration in the plantwide control problem.

where ω1 ) 100 is the weight of squared terms related to CVs, whereas ω2 ) I is the weighted matrix for MVs terms. The MPC tracking problem has the aim to keep CVs as close as possible to their set points, with starting points, respectively, equal to

IVHPset (0) ) 0.575 dL/g IVSSPset (0) ) 0.782 dL/g

(13)

set set and IVSSP(t) values along the producBy moving IVHP(t) tion time, a grade-change problem is set. Nevertheless, the MPC

Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 99

Figure 8. Part of the Jacobian’s Boolean (in black, the no null terms). The first block represents the PE section, blocks 2-4 are referred to chambers of the SE section, the block 5 is referred to the LP unit, and blocks 6-15 concern the IP section.

optimizes MVs around nonoptimal set points, and this is pushing toward the integrated multilevel optimization, where some set points can become degrees of freedom of the businesswide problem and the optimizer should estimate their best value, by taking into consideration the market laws. In this way, the RTDO problem is lowered to the MPC time-scale and simultaneously solved. The optimal set points will be employed in the MPC solution, for calculating real optimal MVs and further improving the plant performance. The main goal of this paper is to optimize the production through the maximization of the end product flow, in the observance of grade specifications. The integrated objective function could be written as follows,

Φ ) ΦMPC + ΦRTDO Φ)

min u∈U,MVtar 15 ∈T

{

(14)

∫0t [ω1(IVHP(t) - IVHPset)2 + ω1 f

MV

(IVSSP(t) - IVSSPset)2 +

2 ||MV(t),n - MVtar ∑ n ||ω n)1 2

tar ] dt} (15) µ15‚MV15,(t,iter)

where the optimizer will simultaneously calculate optimal MVs, degrees of freedom of the MPC problem in every sampling time, and the maximum reachable for the MVtar 15 , which represents the SSP production target, economic degree of freedom of businesswide RTDO problem for each SQP iteration. Weights employed in the numerical solution are ω1 ) 100, ω2 ) 1, and µ15 ) 5. 5. Results 5.1. Numerical Integration of the DAE System. The complete PET process model consists of 1500+ differential and algebraic equations; nevertheless, assumptions made in the modeling of PET stages lead to a particular and favorable structure of the resulting DAE system. In the model, each representative equation block is composed of 16 differential or algebraic equations except some stage that have only 15 equations: by adding a fictitious algebraic equation to these blocks, the system assumes a tailored equal-sized block structure, and efficient DAE solvers for structured-sparsity problems can be used. Moreover, each block in the melt-polymerization section is sequentially solvable: once inputs are known, the

100

Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008

Figure 10. Grade-change transitions in the IV of the HP unit.

Figure 9. Flow chart of the real-time optimization strategy using a hybrid (Windows/Unix) environment

remaining unknown variables of the block are retrievable from the equation subsystem itself, when the stationary solution is searched; in this particular case, all the equations are treated as algebraic. This solution strategy is necessary because the assumption of the phase equilibrium in the first stages made the system index two and the steady solution is achieved with difficulties, if a good guess is not supplied. This property is no longer true for the SSP equations, with these being interconnected by the finite difference approximation of second-order derivatives. Therefore, the overall steady system is practically obtained by independently solving each equation block, by starting from the feed inlet and ending necessarily at the last block of the HP section. The solution of this last block gives initial conditions of the SSP unit that, being represented by a purely ODE system, can be solved by simulating a complete startup. The overall system can now be initialized with this steady solution in order to simulate the dynamic effect of disturbances acting on the process. The tailored structure has clear advantages in the numerical integration of the DAE system, too. In Figure 8, a small part of the Jacobian’s Boolean is shown, in order to demonstrate that the system has a banded and equal-size block diagonal structure. This property allows integrating the system by using the tailored solver of the BzzMath5.0.25 This strategy employs a minimal CPU elapsed time (no more than 3 s on a PentiumIV 2.6 GHz), in order to simulate a transient after a disturbance in the input conditions. For merging the BzzMath5.0 in C++ and the PET model, written in FORTRAN 90, a mixed language is adopted. The resolution of the integrated optimization is performed by using the MUSCOD-II (multiple shooting code for optimization), which is a simultaneous tool based on a direct multiple-shooting method,26,27 which can be extended to the real-time dynamic optimization problem, through the ReTiST28 (real-time strategy tool). Since BzzMath5.0 library works in MS Windows operating system only, whereas MUSCOD-II runs in Unix systems, the real-time application requires a secure shell for file transfer (SSH) technology, and the simplified algorithm of ReTiST is reported in Figure 9. The mathematical model is initialized, fully discretized, and integrated by a numerical solver. Then, the model integration is frozen at time τ, and all the stored data is sent to the UNIX environment, where the optimizer operates. The MUSCOD-II evaluates the best MVs, if we are solving a MPC; otherwise (in a two-level integration problem), optimal set points are

Figure 11. Grade-change transitions in the IV of the SSP unit.

Figure 12. Selected manipulated variables during the simulation.

calculated, too. Result files are sent back to the MS Windows environment, and new optimized values are implemented in the model: the integration of the DAE model restarts for another time interval. The procedure is reiterated till the simulation ending. The integrated multilevel optimization problem can be solved every 2 h, even if, in the presence of stochastic disturbances only or near to stationary conditions, the solution could be obtained in <30 min. In any case, the adopted sampling time allows for speaking about real-time strategy. At every sampling time, the structured NLP problem consists of 9 discretization points, 13 995 variables, 12 480 equality constraints, 28 188 inequality constraints, and selected parameters for the SQP method are as follows: 10-5, the final accuracy of the NLP solution; 10-7, the initial integration tolerance; 10-7, the final integration tolerance; a single-gradient evaluation for each SQP iteration; 4 maximum SQP iteration numbers. In the step immediately after the grade change, one of the hardest time-consuming optimizations, CPU time statistics may be resumed in the following points: • CPU time for sensitivity generation: 2730.27 s (87.00%) • CPU time for state integration: 192.82 s (6.14%) • CPU time for online graphics: 213.39 s (6.80%) • CPU time for constraint reduction: 1.26 s (0.04%) • CPU time for condensing: 0.04 s (0.00%) • CPU time for solution of condensed QPs: 0.18 s (0.01%) • CPU time for remaining calculations: 0.23 s (0.01%)

Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 101

Figure 13. Inlet nitrogen temperature at the bottom of the SSP pelletizer.

Figure 14. Effective production, 0.5 MMt/y of pellets, at the end of the PET process (SSP unit).

The effective CPU without online graphics goes down to 50 min, allowing one to solve the integrated optimal problem in real time. 5.2. Polymer-Grade Change and Integrated Solution. Finally, the above-mentioned methods, techniques, and models are adopted in the polymer-grade-transition case. Let us adopt a market that roughly imposes a quality variation every 3 days and where the grade change is required in both HP and SSP end products, as reported in Table 17. The MPC is able to drive the PET process during the whole simulation time, as illustrated in the following trends. In the HP unit, the process is able to rapidly answer to grade transitions, as reported in Figure 10, with a transient period of 15 h in the first step and shorter transients, <2 h, in the second and third changes, considering an uncertainty in the end product equal to 0.0015 dL/g. The IVHP is around the set point (and not perfectly overlapped) because penalties forced among MVs and their targets avoid the complete achievement of control objectives. Instead, transients in SSP, reported in Figure 11, appear longer above all for the characteristic time of the solid-state unit, which can require a residence time of 20 h, if the desired IVSSP is very high. Therefore, it induces long periods for reaching onspecs, which sometimes overcome the day. In the present study, the selected uncertainty is fixed at 0.005 dL/g, but customers could also require a smaller region of confidence, and off-spec times become inevitably larger.

Figures 12 and 13 report temperatures of reactors in the melt section and the inlet nitrogen temperature at the bottom of the pelletizer, respectively. The optimizer imposes the largest moves in sampling times immediately successive to the grade change and then slowly regulates optimal values for MVs in the following intervals. By increasing the number of SQP iterations, the optimization may result faster and the off-spec period can be reduced; nevertheless, the real-time solution is no more guaranteed. Even if starting and arrival IV values are the same, initial and ending reactor temperatures are different: the PET process is free to choose a new steady state, dictated by the history of the last grade change and in accordance with boundaries and path constraints. Figure 14 illustrates the solid PET production status. The target of 0.5 MMt/y is usually preserved, with the exception of transient times, where the production has a decrease until the MPC is near the optimum. The trend is strictly connected to the inlet gas temperature, which regulates the correct grade change. A first comparison between the classical grade change and the integrated multilevel optimization is reported in Figures 15 and 16. In the melt section, improvements in the final IVHP, dictated by the integrated optimization, are apparently small (Figure 15), even if every variation imposed inside these 5 sections becomes fundamental to obtain the best initial conditions for the SSP unit, where the exiting polymer is characterized

Table 17. Sequence of Three Grade-Change Transitions in Both HP and SSP Sections (9 Day Period) and the Required Uncertainty in the Quality Specification (d Indicates the day and h Indicates the Hour of Each Starting-Point Transition) set point

init. condition (d ) 1; h ) 0)

grade change (d ) 1; h ) 2)

grade change (d ) 4; h ) 0)

grade change (d ) 7; h ) 0)

uncertainty (dL/g)

IVHPset (dL/g) IVSSPset (dL/g)

0.575 0.782

0.598 0.813

0.587 0.798

0.575 0.782

(0.0015 (0.005

102

Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008

Figure 15. Grade-change transitions in the IV of the HP unit: common grade change (light solid line), integrated multilevel optimization (dark solid line), and desired product quality (dotted line).

Figure 16. Grade-change transitions in the IV of the SSP unit: common grade change (light line), integrated multilevel optimization (dark thick line), and desired product quality (dark thin line). Banded areas represent quality regions.

Figure 17. Reactor temperature comparisons: common grade change (light solid line) and integrated multilevel optimization (dark solid line).

Figure 18. Inlet gas temperature at the bottom of the SSP pelletizer: common grade change (light line) and integrated multilevel optimization (dark line) comparison.

by the highest added value. In this instance, apart from small differences in IVHP trends, benefits are clearly shown in the solid-state section. Since large improvements in the overall revenues of the production site can be reached by increasing

the SSP production and/or minimizing off-spec times in the solid-state unit, the dynamic optimization has a reduced impact on the melt section, where added values for low-IV polymers are significantly smaller than those for the high-IV polymer

Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008 103

Figure 19. Increase of the overall production: during the steady state (left side) by optimizing setpoint trajectories and during the transient (right side) by reducing off-spec times.

exiting by the SSP. Note that, in both cases, penalties on MVs avoid the complete achievement of the desired set point. Nevertheless, the final grade obtained is composed in the quality region. In this context, IVSSP trends, illustrated in Figure 16, clearly show the gap between the MPC and the integrated optimization: grade-change transitions appear strongly faster, guaranteeing the achievement of quality regions some hours beforehand, by decreasing off-specs. This gap is enhanced when the quality region is reduced. This peculiarity of the integrated optimization is also confirmed by other simulations, with different disturbances and different entities. Moreover, in this particular case, on-spec times are improved by one-half day. Key comparisons of MVs are shown in Figures 17 and 18. Generally speaking, trends present similar behaviors, especially in the SE, IP, and HP units. Instead, in the low polymerizer and the solid-state pelletizer, the gap can overcome 2 K. All these smaller than significant variations allow the increase in the PET production, as qualitatively reported in Figure 19. With the selected series of product grades, the integrated multilevel optimization allows higher production without any variation in the inlet feed flows but simply managing units in an optimal way. On the left side of Figure 19, the gap between the integrated multilevel optimization and the MPC about the final production of the PET plant is reported. By considering the same time of on-spec production during steady-state conditions, the net profit margin is inevitably increased by using the dynamic optimization. Moreover, by considering that faster transients are obtained by using the businesswide RTDO, the final production results furthermore higher. A lot of series of grade changes have been considered: in any case, a significant increase in the onspec zone, and then in the overall polymer production, has been detected. 6. Conclusions The integration of two (or more) optimal control levels via closed-loop technique allows industrial plants to operate inside their stability region and, at the same time, to increase the production (rising efficiencies of unit operations) through the definition of optimal economic trajectories. Moreover, the integrated optimization is able to expand the on-spec time since transients are usually faster than a typical grade change managed by the MPC. Future challenges will concern the introduction of a detailed market demand and the development of MILP/ MINLP problems for solving the production scheduling together with the integrated operational supply chain. Acknowledgment The authors thank the Interdisciplinary Centre for Scientific Computing, IWR department “Ruprecht-Karls”, University of Heidelberg, for the MUSCOD-II package and Prof. M. Diehl

from ESAT Department, Catholic University of Leuven, for his insightful advice. The authors also wish to thank the technicians of “Sistemi e Tecnologie Infromatiche”, CMIC department, Politecnico di Milano, for their support in the cluster room setting up and ReTiST implementation. Nomenclature Acronyms CHV ) characteristic of valve CPU ) central processing unit CSTR ) continuous stirred tank reactor CV(s) ) controlled variables DAE ) differential algebraic equation(s) HP ) high polymerizer IP ) intermediate polymerizer IV ) intrinsic viscosity IVHP ) intrinsic viscosity HP IVSSP ) intrinsic viscosity SSP LP ) low polymerizer MPC ) linear/nonlinear model predictive control MV(s) ) manipulated variable(s) MW ) molecular weight NC ) number of components NLP ) nonlinear programming ODE(s) ) ordinary differential equation(s) OLE ) object linking and embedding OPC ) OLE for process control PDE(s) ) partial differential equation(s) PE ) primary esterifier PHCR ) preheater crystallizer RTDO ) real-time dynamic optimization SE ) secondary esterifier SQP ) sequential quadratic programming SSH ) secure shell file transfer SSP ) solid-state polymerizer Chemical Species AA ) acetaldehyde DEG ) diethylene glycol DMT ) dimethylterephthalate EG ) ethylene glycol PET ) poly(ethylene terephthalate) PTA ) pure terephthalic acid TPA ) terephthalic acid VIN ) vinyl W ) water Prefixes b ) repeat group t ) end group

104

Ind. Eng. Chem. Res., Vol. 47, No. 1, 2008

Superscripts * ) interface d ) discharge in ) entering liq ) liquid phase out ) exiting p ) partial rec ) recycle set ) set point sol ) solid phase tar ) target tot ) total vap ) vapor Subscripts iter ) number of iteration j ) j-esime component mix ) mixture conditions ref ) reference conditions Greek Letters R ) solid TPA solubility in oligomer solution Φ ) general objective function F ) process flow density Other Ai ) pre-exponential factor Area ) area/surface Ci ) concentration of species Cp ) specific heat Ei ) activation energy Gi ) generation rate ki ) forward rate constant Ki ) equilibrium constant P ) pressure Ri ) reaction rate T ) temperature U ) overall heat exchange coefficient V ) volume ω ) generic weight Literature Cited (1) Morari, M.; Lee, J. H. Model predictive control: Past, present and future. Comput. Chem. Eng. 1999, 23, 667-682. (2) Biegler, L. T.; Grossmann, I. E. Retrospective on Optimization. Comput. Chem. Eng. 2004, 28, 1169-1192. (3) Vettenranta, J.; Smeds, S.; Yli-Opas, K.; Sourander, M.; Vanhamaki, V.; Aaljoki, K.; Bergman, S.; Ojala, M. Dynamic Real Time Optimization Increases Ethylene Plant Profits. Hydrocarbon Process. 2006, 10, 59-66. (4) Bassett, M. H.; Dave, P.; Doyle, F. J.; Kudva, G. K.; Pekny, J. F.; Reklaitis, G. V.; Subrahmanyam, S.; Miller, D. L.; Zentner, M. G. Perspectives on Model Based Integration of Process Operations. Comput. Chem. Eng. 1996, 20, 821-844. (5) Wang, Y.; Seki, H.; Ohyama, S.; Akamatsu, K.; Ogawa, M.; Ohshima, M. Optimal grade transition control for polymerization reactors. Comput. Chem. Eng. 2000, 24, 1555-1561. (6) Kadam, J. V.; Schlegel, M.; Marquardt, W.; Tousain, R. L.; van Hessem, D. H.; van der Berg, J.; Bosgra, O. H. A Two-level Strategy of Integrated Dynamic Optimization and Control of Industrial ProcessessA Case Study. Presented at ESCAPE-12, The Hague, The Netherlands, 2002; pp 511-516.

(7) Chatzidoukas, C.; Perkins, J. D.; Pistikopoulos, E. N.; Kiparissides, C. Optimal grade transition and selection of closed-loop controllers in a gas-phase olefin polymerization fluidized bed reactor. Chem. Eng. Sci. 2003, 58, 3643-3658. (8) Cervantes, A. M.; Tonelli, S.; Brandolin, A.; Bandoni, J. A.; Biegler, L. T. Large-scale dynamic optimization for grade transitions in a low density polyethylene plant. Comput. Chem. Eng. 2002, 26, 227-237. (9) Benamor, S.; Doyle, F. J.; McFarlane, R. Polymer grade transition control using advanced real-time optimization software. J. Process Control 2004, 14, 349-364. (10) Tremblay, D. A. Using Simulation Technology to Improve Profitability in the Polymer Industry. Presented at AIChE Spring 1999 Meeting, Houston, TX, 1999. (11) Rovaglio, M.; Algeri, C.; Manca, D. Dynamic modeling of a poly(ethylene terephthalate) solid-state polymerization reactor. II. Model predictive control. Ind. Eng. Chem. Res. 2004, 43, 4267-4277. (12) Lu, J. Z. Challenging Control Problems and Emerging Technologies in Enterprise Optimization. Presented at DYCOPS-6, Chejudo Island, 2001. (13) Kameswaran, S.; Biegler, L. T. Simultaneous dynamic optimization strategies: Recent advances and challenges. Comput. Chem. Eng. 2006, 30, 1560-1575. (14) Flores-Tlacuahuac, A.; Biegler, L. T.; Saldivar-Guerra, E. Optimal grade transitions in the high-impact polystyrene polymerization process. Ind. Eng. Chem. Res. 2006, 45, 6175-6189. (15) Lang, Y. D.; Biegler, L. T. A Software Environment for Simultaneous Dynamic Optimization. Comput. Chem. Eng. 2007, 31, 931-942. (16) Joly, M.; Pinto, J. M. Optimal Control of Product Quality for Batch Nylon-6,6 Autoclaves. Chem. Eng. J. 2004, 97, 87-101. (17) HPI Construction Boxscore. Supplement to Hydrocarbon Processing; Kane, L. A., Romanow, S., Eds.; Hydrocarbon Processing: Houston, TX, 2005. (18) Kang, C. K.; Lee, B. C.; Ihm, D. W.; Tremblay, D. A. A simulation study on continuous direct esterification process for poly(ethylene terephthalate) synthesis. J. Appl. Polym. Sci. 1997, 63, 163-174. (19) Kim, I. S.; Woo, B. G.; Choi, K. Y.; Kiang, C. Two-phase model for continuous final-stage melt polycondensation of poly(ethylene terephthalate). III. Modeling of multiple reactors with multiple reaction zones. J. Appl. Polym. Sci. 2003, 90, 1088-1095. (20) Kang, C. K. Modeling of solid-state polymerization of poly(ethylene terephthalate). J. Appl. Polym. Sci. 1998, 68, 837-846. (21) Choi, K. Y.; Laubriet, C.; Lecorre, B. Two-phases model for continuous final-stage melt polycondensation of poly(ethylene terephthalate). I. Steady-state analysis. Ind. Eng. Chem. Res. 1991, 30, 2. (22) Choi, K. Y.; Saint Martin, H. C. Two-phases model for continuous final-stage melt polycondensation of poly(ethylene terephthalate). II. Analysis of dynamic behavior. Ind. Eng. Chem. Res. 1991, 30, 1712. (23) Yao, K. Z.; McAuley, K. B.; Berg, D.; Marchildon, E. K. A dynamic mathematical model for continuous solid-phase polymerization of nylon 6,6. Chem. Eng. Sci. 2001, 56, 4801-4814. (24) Algeri, C.; Rovaglio, M. Dynamic modeling of a poly(ethylene terephthalate) solid-state polymerization reactor. I. Detailed model development. Ind. Eng. Chem. Res. 2004, 43, 4253-4266. (25) Buzzi Ferraris, G. BzzMath: Numerical libraries in C++; Politecnico di Milano: Milan, Italy, 2003. (26) Bock, H. G.; Plitt, K. J. A Multiple Shooting Algorithm for Direct Solution of Optimal Control Problem; In Proceedings 9th IFAC World Congress; Pergamon Press: Elmsford, NY, 1984. (27) Leineweber, D. B. Efficient Reduced SQP Methods for the Optimization of Chemical Process Described by Large Sparse DAE Models; Fortschr. Ber. VDI Reihe 3, Verfahrenstechnik; VDI Verlag: Dusseldorf, Germany, 1999. (28) Manenti, F.; Rovaglio, M.; Ferreau, H. J.; Diehl, M. Dynamic Optimization in Polymer Plants. Presented at MCCEE, Venezia, Italy, 2006; pp 252-256.

ReceiVed for reView March 8, 2007 ReVised manuscript receiVed September 13, 2007 Accepted September 21, 2007 IE070355S