Ssp

Ind. Eng. Chem. Res. 2004, 43, 4267-4277

4267

Dynamic Modeling of a Poly(ethylene terephthalate) Solid-State Polymerization Reactor II: Model Predictive Control Maurizio Rovaglio,* Carlo Algeri, and Davide Manca Dipartimento di Chimica, Materiali e Ingegneria Chimica “Giulio Natta”, Politecnico di Milano, Piazza Leonardo Da Vinci 32, 20133 Milano, Italy

A detailed dynamic model for the poly(ethylene terephthalate) (PET) solid-state polymerization (SSP) process is employed for the design of an advanced control scheme based on model predictive control (MPC). The original model of the reactor, consisting of a system of 16 time-dependent partial differential equations in two spatial coordinates, is simplified by contracting the radial direction. Such a reduced model is then incorporated into an MPC methodology. The overall control scheme is simulated to verify its reliability by maintaining the reactor temperature and the final intrinsic viscosity (IV) of the outlet PET resin constant. The validation of the resulting control loop composed of the MPC controller and the SSP unit is performed through the use of a detailed model that substitutes for the real plant. Control parameters such as the prediction and control horizons and the sampling time are tested and optimized. Control performances in response to unmeasured disturbance and set-point variations are illustrated and analyzed. Introduction

Table 1. Kinetic Scheme of the Reduced Model

Model predictive control (MPC) technology represents the present and the future of chemical process control. Its fundamentals reached popularity in the early 1970s. Currently, the availability of powerful computers has permitted the application of such model predictive controllers to real systems, and in the future, an everincreasing interest in this field is expected. More details on the history and trends of this technology can be found in Morari and Lee.1 The MPC strategy adopts a process model, which can be the same as that used for process design, with a high degree of detail, to predict the future behavior of the controlled system and to perform the calculation of the optimal set of control actions, dictated by the minimization of a common objective function. This function is structured to increase in value when the system moves away from its steady state, in which case the controlled variables differ significantly from their set points. Bounds and constraints can also be taken into account for specific requirements on both manipulated and controlled variables. This feature make MPC capable of successfully treating multivariable systems with critical characteristics such as strong nonlinear behavior, inverse responses, long dead times, or marked interactions among variables. The degree of detail adopted for the translation of the real process behavior into mathematical terms synthesizes the challenge represented by MPC, which does not impose any limits on its model structure, thus improving on many traditional control techniques that, today, might be rendered obsolete. A practical restriction derives only from the ability to obtain detailed models with a computing time that satisfies the limits needed for the control to be implemented in a real-time system. The solid-state polymerization (SSP) process is exactly such a system, as already observed by Krishnan et al.,2 for which the choice of model size becomes critical for MPC applications but, at the same time, conventional control shows its weakness. Such poor controllability is caused * To whom correspondence should be addressed. E-mail: [email protected]. Fax: +39.02.7063.8173.

rate constant reaction 1 2

TEG + tTPA ‚ bEG + bTPA + W TEG + tEG ‚ bEG + EG

forward

reverse

k1 k2

k1/K1 k2/K2

by the long dead time of the plug-flow reactor and by the weak effects of the manipulated variables, which also cause a low efficiency for traditional PID controllers. The next section deals with the analytical reduction of the detailed model, which was necessary to obtain a practical MPC controller for the reactor. The original model of the reactor (see part I of this work3), consisting of a system of 16 time-dependent partial differential equations in two spatial coordinates, was simplified by contracting the radial direction into a unique grid point and substituting the corresponding derivatives with simplified linear relationships accounting for the mass and thermal diffusion into the solid. The model was solved with the VLUGR2 library (by Blom et al.4), and the resulting savings of CPU time was around 2 orders of magnitude, dropping from 1 h needed by the complete model to a few seconds required by the reduced one. In a subsequent step, this model was incorporated into a model predictive control structure, and the overall control loop was implemented, through the use of the detailed model as a “virtual plant”, to demonstrate the effectiveness of both the designed control scheme and the model reduction technique. Model Reduction The model simplification was realized following similar assumptions as highlighted by Yao et al.5 in the third part of their work on SSP for the nylon process. Substantially, it consists of the removal of the radial dimension by substituting it with a single node and translating the description of mass and heat diffusion within the polymer into simple linear relationships. Moreover, the kinetic scheme was simplified by taking into account only the two main reactions (esterification and transesterification), reported in Table 1, with the

10.1021/ie034324i CCC: $27.50 © 2004 American Chemical Society Published on Web 06/22/2004

4268 Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 Table 2. Generation Rate Gj(t) and Reaction Rate Rj(t) Terms Generation Rate Terms Gj(t) GEG ) R2 GtEG ) -R1 - 2R2 GW ) R 1 GtTPA ) -R1

∂C ˜j 3 ) -k0,j(p˜ j - pj,g) + ∂t R Gj(t) -

∂C ˜j ∂2C ˜j + Db 2 (1) (1 - )acFp ∂z ∂z m ˘p

Reaction Rate Terms Rj(t)

CbEG R1 ) k1CtEGCtTPA - 2(k1/K1)CWCbTPA CtEG + CbEG R2 ) k2CtEGCtEG - 4(k2/K2)CEGCbEG Table 3. Legend of Abbreviations for Molecular Species

The parameter k0,j can be considered a global masstransfer coefficient of component j, and it depends on two primary coefficients: kg,j corresponding to the interphase mass transfer and kp,j related to the internal diffusion resistance. The former can easily be calculated from the Colburn analogy, whereas the latter has to be correlated to the diffusion, Dp,j, which necessitates some further insights, as discussed later. Once the single coefficients are calculated, k0,j can be simply obtained by expressing the sum of two resistances

HjMj 1 1 ) + k0,j kp,j kg,j

(2)

where Hj is the Henry’s law pseudo-equilibrium constant calculated as the ratio between p˜ j, the equilibrium vapor pressure obtained from the Flory-Huggins equation, and the volume fraction φj of component j in the solid particle

Hj )

p˜ j φj

(3)

With similar assumptions, the following modified expression for the thermal balance of the solid particles moving into the SSP reactor can be easily obtained

∂T ˜s ∂T ˜s m ˘p h0 3 (T - T ˜ s) + )+ ∂t FpCp g R Fp(1 - )ac ∂z corresponding reaction rates reported in Table 2. Therefore only nine of the original 16 state variables involved are still needed, namely, the ethylene glycol (EG), water (W), EG end group (tEG), and terephthalic acid end group (tTPA) concentrations within the solid particle; the EG and W concentrations in the gas phase; the solid and gas temperatures, Ts and Tg, respectively; and the crystallinity χc. The concentrations of the TPA and EG repeat units (bTPA and bEG, respectively) can be obtained algebraically, and they do not lead to further differential equations. (See Table 3 for the abbreviations for the molecular species considered in this work.) Moreover, for the evaluation of the gas-phase properties, the hypothesis of a pure nitrogen gas flow was employed because the concentrations of byproducts are negligible, being always below 0.1%. The model reduction method adopted here, suggested by Yao et al.,5 substitutes Fick’s law with a simple linear correlation that defines the interphase and internal transfers by coupling their effects in a unique constant coefficient. A reduction in the radial coordinate implies that the computation of the concentration and temperature gradients within the solid particle is no longer needed. Obviously, all state variables y now have the property of being averaged, y˜ , along the particle radius, corresponding to a general representation given as y˜ ) y˜ (z,t) rather than y ) y(r,z,t). Therefore, the material balances of the diffusing compounds j (i.e., EG and W) in a plug-flow reactor scheme assume the following structure

Db

∂2T ˜s ∂z

2

+

Gc(t)(-∆HCR) ∆Hev,W 3 - k0,W (p˜ - pW,g) Cp,p FpCp W R ∆Hev,EG 3 (p˜ EG - pEG,g) (4) k0,EG FpCp R

The thermal exchange between solid and gas can be evaluated as a linear function of the driving force and, using the lumped coefficient h0, obtained in a complete analogy with coefficient k0,j for the mass balance. Similarly to eq 2, the interphase heat-transfer coefficient hg can be coupled with the internal heat conduction coefficient hp, which primarily depends on the thermal conductivity κp, giving the global thermal coefficient as a sum of the two resistances

1 1 1 ) + ho hp hg

(5)

Obviously, the mass and thermal balances for the gas phase are not affected by the model reduction because their structure does not depend on the particle radial coordinate. The only observable changes are those related to the new transfer coefficients adopted, k0,j and h0, which must also appear in gas-phase equations. Consequently, the terms reported in the driving force ˜ s, implicitly for the mass and heat flux must be p˜ j and T referring to the average values of the flux into the particle, and not to pj,R and Ts,R which are related to the surface values.

Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 4269

Finally, the mass and thermal balances for the gas phase assume the following form

∂Cj,g ∂2Cj,g Dg,j ∂ac ∂Cj,g m ˘ g ∂Cj,g + ) + Dg,j + ∂t Fgac ∂z ac ∂z ∂z ∂z2 3(1 - ) (6) k0,j(p˜ j - pj,g) R ∂Tg ∂t

)

m ˘ g ∂Tg Fgac ∂z

+

∂2Tg

κg

FgCp,g ∂z2

+

FgCp,gac dz ∂z

6(TR - Tin)

∞

∑ n)1

2

1

n

e-Rn π t/R 2 2

2

2

dT ˜ 4 ) 4πR2hp(T ˜ - TR) - πR3FpCp 3 dt

(8)

(9)

where hp is the laminar transfer coefficient representing the polymer-side resistance. Differentiating the expression of T ˜ in eq 8 gives

dt

6(TR - Tin)R R

2

∑

∞

e-Rn π t/R ∑ n)1 2 2

2

(10)

Then, replacing dT ˜ /dt in eq 9 and solving for hp, one obtains

(Tin - TR) ∞ 2 2 2 e-Rn π t/R hp ) 2κp R(T ˜ - TR)n)1

∑

(11)

A further simplification can be realized rearranging eq 8 as

e

(12)

-Rn2π2t/R2

2

which, once substituted into eq 11, gives rise to the following new relationship

hp )

where R ) κp/FpCp is the thermal diffusivity inside the sphere. Therefore, defining the rate of heat transfer from the particle body to its outer surface through a linear correlation, the particle heat balance can be formulated as follows

)

6

1

n)1n

-

The coefficients hp and kp,j previously defined in eqs 2 and 5 and representing the resistances to mass and heat diffusion into the solid particle, need some further clarification about their evaluation. The method for estimating hp, reported below, was proposed by Yao et al.,5 and for sake of brevity, the description of the corresponding procedure for the mass-transfer coefficients is omitted because it is completely analogous to that for hp. The main idea is to adapt the analytical solution of the heat conduction problem of a sphere with uniform initial temperature Tin and with its outer surface held at a constant temperature TR. Carslaw and Jaeger6 provided the following solution for such a problem

dT ˜

∞

κpπ2

∑

π

T ˜ - TR

π2

)

∞

dac ∂Tg

κg

2vxacπ + k0,j(p˜ j - pj,g) hw(Tg - TW) FgCp,gac j)EG,W (T ˜ s - Tg)3(1 - ) 3(1 - ) Cpv,j ˜ s - Tg) + h0(T (7) FgCp,gR FgCp,gR

T ˜ ) TR -

Tin - TR

e-Rn π t/R ∑ n)1

∞

3R

∑

2 2

1

n)1n

e

2

(13)

-Rn2π2t/R2

2

This last equation allows for the evaluation of the polymer-side heat-transfer coefficient that has the intrinsic characteristic of changing with time: as t f 0, hp f ∞, and as t f ∞,hp f κpπ2/3R. Moreover, hp also depends on the pellet dimension and its physical thermal properties. In their work, Yao et al.5 also suggested the convenient introduction of a time-invariant value of this lumped parameter equal to the infinite-time heattransfer coefficient. Assuming such a hypothesis for hp, eq 9 can be solved analytically to obtain

T ˜ ) TR + (Tin - TR)e-3hpt/FpCpR

(14)

Results given by the analytical solution and the simplified one can be compared and an adjusting factor ch can be introduced as hp ) chh∞p . In this study, this correction factor has been used instead as a fitting parameter between models, as suggested by Yao et al.5 The lumped parameter kp,j for mass transfer can be determined following exactly the same method. Obviously, the value of kp,j depends on the diffusivity of component j and on the pellet dimensions. Therefore, in complete analogy with hp, when t f 0, kp,j f ∞, whereas when t f ∞, kp,j f Dw,jFpπ2/3R. Neglecting the time dependency and choosing the constant infinite-time value of kp,j ) Dw,jFpπ2/3R, the dynamic material balance of component j can be solved analytically to give

C ˜ j ) Cj,R + (Cj,in - Cj,R)e-3kp,jt/RFp

(15)

Even in this case, an adjustment factor of ck,j can been adopted as a correction between the simplified and analytical solutions or extended as a fitting parameter between the detailed and simplified models (kp,j ) ck,j ∞ kp,j ). The system equations and initial/boundary conditions of the resulting reduced model are summarized in Tables 4 and 5, respectively. The system of timedependent one-dimensional partial differential equations was numerically solved through the same library, VLUGR2,4 as adopted for the detailed model. Although this numerical solver is designed for time-dependent two-dimensional partial differential equations, it was possible to solve the one-dimensional system by imposing a fictitious coordinate (replacing the radial one) constituted by a unique grid node. The CPU time saved is remarkable: the reduced model requires a computing time 2 orders of magnitude smaller than the detailed

4270 Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 Table 4. Reduced Model Equations

∂C ˜j ∂C ˜j ∂2C ˜j m ˘p 3 ) -k0,j(p˜ j - pj,g) + Gj(t) + Db 2 ∂t R (1 - )acFp ∂z ∂z

(19)

m ˘p ∂C ˜j ∂C ˜j ∂2C ˜j ) Gj(t) + Db 2 ∂t (1 - )acFp ∂z ∂z

(20)

∂χc ∂2χc m ˘p ∂χc )+ kc(χmax - χc) + Db 2 ∂t Fpac(1 - ) ∂z ∂z

(21)

∂Cj,g ∂2Cj,g Dg,j ∂ac ∂Cj,g m ˘ g ∂Cj,g 3(1 - ) ) + Dg,j + k0,j(p˜ j - pj,g) + ∂t Fgac ∂z ac ∂z ∂z R ∂z2

(22)

∂T ˜s ∂T ˜s ∂2T ˜ s Gc(t)(-∆HCR) m ˘p h0 3 (Tg - T ˜ s) + Db 2 + )+ ∂t FpCp R Cp,p Fp(1 - )ac ∂z ∂z

(23)

∆Hev,W ∆Hev,EG 3 3 k0,W (p˜ - pW,g) - k0,EG (p˜ EG - pEG,g) FpCp W R FpCp R 2vxacπ ∂Tg dac ∂Tg m ˘ g ∂Tg κg ∂2Tg κg ) + - hw(Tg - TW) + + ∂t Fgac ∂z FgCp,g ∂z2 FgCp,gac dz ∂z FgCp,gac j)EG,W

∑

k0,j(p˜ j - pj,g)Cpv,j

j

(T ˜ s - Tg)3(1 - ) FgCp,gR

(24)

3(1 - )

˜ s - Tg) + h0(T FgCp,gR

Table 5. Initial and Boundary Conditions of the Reduced Model Cj|t)0 ) C0j Cj|z)H ) C0j

|

∂Cj ∂z

z)0

∀ z ∈ {0 e z e H}

)0

j ) EG, W, tTPA,tEG 0 Cj,g|t)0 ) Cj,g ∀ z ∈ {0 e z e H} j ) EG, W, AA 0 Cj,g|z)H- ) Cj,g |z)H -

|

∂Cj,g ∂z

z)0

|

Dg,jac ∂Cj,g m ˘g ∂z

z)H-

)0

j ) EG, W, AA χc|t)0 ) χ0c ∀ z ∈ {0 e z e H} χc|z)0 ) χ0c

∂χc ∂z

|

z)H

)0

Ts|t)0 ) T0s Ts|z)0 ) T0s

|

∂Ts ∂z

z)H

)0

Tg|t)0 ) T0g

∂Tg ∂z

|

z)0

∀ z ∈ {0 e z e H}

∀ z ∈ {0 e z e H}

)0

Tg|z)H- ) T0g|z)H -

Figure 1. IV profiles along the reactor height calculated with the detailed (solid lines) and reduced (dashed lines) models at different reaction times.

information returned by the reduced model. Therefore, some checks are needed to verify the feasibility of the derived tool. To carry out a comparison between model performances, it becomes necessary to introduce the molecular average number and the related intrinsic viscosity, IV

Xn )

|

κj,gac ∂Tg m ˘ gCp,g ∂z

z)H-

model, mainly because of the simplifications that lighten the numerical iterations. Conversely, the reduction in terms of the number of state variables and the other assumption adopted give rise to minor amounts of

CtEG + CtTPA + CbEG + CbTPA CtEG + CtTPA

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

(16) (17)

In Figure 1 is shown a comparison between the dynamic evolutions of the IV along the reactor height when calculated with the detailed model (solid lines) and with the reduced model (dashed lines). A good agreement between the two models was obtained by setting

Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 4271

Figure 2. IV profiles at the reactor outlet predicted by reduced and detailed models.

Figure 4. EG concentration profiles at the reactor outlet predicted by the reduced and detailed models.

Figure 3. Particle temperature profiles along the reactor height predicted by the reduced and detailed models at different reaction times.

Figure 5. W concentration profiles at the reactor outlet predicted by the reduced and detailed models.

the adjustment factors, ck,j, to 3.5 for ethylene glycol and to 1.5 for water. These relatively high values are mainly due to the reduction of the kinetic scheme. Moreover, it is important to verify the analogy of the dynamic behaviors for the two models and, hence, to determine whether the variable profiles at the reactor outlet as a function of time can provide reasonably the same information. Figure 2 shows a comparison of the outlet IV predicted by the models, and the presented curves confirm a very good agreement. In addition, Figure 3 compares the particle temperature profiles obtained from the two models. The corresponding curves essentially overlap, even if the corresponding lumped parameter hp is not adjusted by any fitting parameter ch. The last comparison, illustrated in Figures 4 and 5, is of the two diffusing components EG and W. Their profiles at the reactor outlet, as determined by the detailed and reduced models, show a good agreement for the predicted dynamics even though the water concentration has a significant error at the steady state because of the simplification of the kinetic scheme. However, with the aim of developing an MPC controller, the presence of a tolerable deviation between the predictions of the simplified model and the detailed model (representing the real plant) can be accepted,

given that a feedback error corrector term is always incorporated in the controller structure (for details, see, for example, Morari and Lee1). Model Predictive Control Configuration From its original theory, the MPC algorithm calculates, during each iterative cycle repeated after a time step equal to the sampling time ts, the optimal profiles of the manipulated variables along a predefined temporal horizon, the so-called control horizon. The set of best possible actions is evaluated by minimizing an objective function, which is determined mainly by the difference between the output of the model and the corresponding imposed set points, predicted along a prediction horizon. Usually, these two temporal horizons are represented as multiples of the sampling time, and their lengths can be indicated with nHC and nHP, respectively. A control horizon of nHC means that the objective function will be minimized by a control action at each time step between 1 and nHC with a total number of degrees of freedom equal to nHC times the number of controlled variables. The prediction horizon is longer than, or at least equal to, the control horizon, and therefore, in any control actions between nHC and nHP, the manipulated variables will maintain the last value achieved in nHC. At the end of this optimizing

4272 Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004

A general-purpose objective function, also adopted in this study, can be given as follows

{ [ ]} ∑ (∑{ [ [

k+nHP

F)

nc

∑ c)1 ∑ ωc

yc(i) - yc,set(i)

i)k+1

k+nHC-1 nm l)k

Figure 6. Schematic illustration of an MPC computing cycle.

procedure, the first control action, among the nHC mentioned above, is implemented and maintained the same throughout the sampling time, after which a new measurement becomes available, and the calculation cycle can restart. During the optimization computation, a feedback correction for the model prediction must also be taken into account. In fact, the error between the controlled variables computed by the model and the real measurements is estimated at the beginning of the evaluation cycle, and this error is employed to correct all future model predictions along the new optimization step. Such a strategy avoids a model prediction of the process behavior that is not in agreement with the real state of the system. Therefore, this correction method allows one to reject and control unmeasured disturbances. To better illustrate the MPC computation procedure mentioned above, Figure 6 shows schematically the main steps, which can be summarized as follows: (1) At time t0 (iHC ) 1), the measurement y0 of the current process status becomes available. (2) The optimization, constrained or not, is executed in real time, minimizing the objective function F ) F(y,u) along the prediction horizon nHP and supplying the optimal profile u j (t) of manipulated variables along the control horizon nHC. (3) The first control action, u j 0, of the optimal profile, determined at step 2, is implemented on the process, and the cycle restarts at t0 ) t0 + ts returning to step 1. It is clear that, to have a feasible control solution, the optimization problem must be solvable in a computing time less than or equal to the sampling time ts. Consequently, this requirement can also be applied to the model that has to be fast enough to allow the optimization procedure to be performed within the feasible time limit. The core of the MPC control philosophy is the objective function, with regard to which some interesting features can be highlighted. In fact, the MPC control actions are obtained not by evaluating only the absolute error (y - yset), as in conventional control, but rather by adopting a more detailed minimization objective function, which can be structured by taking into account several constraints and different operating policies. In a general formulation, not only the absolute error on controlled variables can be penalizing but also the excessive deviation of the corresponding manipulated variables from their steady-state values. Soft and hard operating constraints can also be imposed, allowing the real systems to be fully managed while the predictive model takes into account the coupling between variables and actions, thereby making easier an effective multivariable control structure.

m)1

2

+ yc,set(i) um(l) - um(l - 1) 2 ωm,M + um(l - 1) um(l) - um,SST(l) 2 ωm,T um,SST(l)

]

] })

(18)

The structure is a quadratic function in terms of error for each controlled variable yc with respect to the corresponding set points yc,set along the prediction horizon hp and also in terms of the deviation for each of manipulated variables um along the control horizon hc. These last terms can be determined in two different forms. The first is the difference between two subsequent values of the manipulated variable um(l) and um(l - 1), which is include to avoid control actions that are too strong. The second is the difference between the actual value of the manipulated variable um and the corresponding value of its desired steady-state condition um,SST (steady-state target). Generally, the importance of each single contribution of the objective function is managed through the weights ωc, ωm,M, and ωm,T. SSP Control Scheme Defining which variables have to be controlled and which have to be manipulated is clearly the first step of any control system setup. In the previous paper (see part I of this work3), the analysis pointed out that only the purge gas temperature and byproduct concentrations can be used efficiently to control the process outputs. The more obvious controlled variable is the intrinsic viscosity IV, which is the primary qualitycontrol parameter related to the average polymer chain length, but also a temperature control is needed because this parameter is limited by a maximum value dictated by the sticking polymer temperature. In the literature, similar conclusions can easily be found. In fact, in their simulations, Mallon and Ray7 showed that an increased concentration of condensates in the gas phase can influence significantly the product quality, whereas Leffew et al.8 emphasized the role of the purge gas flow and its temperature to better control the solid-bed temperature profile and consequently the polymer intrinsic viscosity. Direct control of the IV at reactor outlet requires a continuous measurement of this variable. For the similar nylon SSP process, Krishnan et al.2 adopted an industrial scheme in which the polymer molecular weight, and consequently the IV, is inferred through the pressure difference between the discharging point of the extruder (located at reactor outlet) and the corresponding spinning manifold. Therefore, because the two processes are very similar, this type of measurement is also assumed to be available in this work. Moreover, Krishnan et al.2 adopted an additional outlet water concentration control. In summary, all of the possible controlled variables can be listed as follows: (1) the intrinsic viscosity, IV, of PET at reactor outlet (IVout); (2) the temperature of the outlet inert gas (Tg,out); and (3) the water concentration in the outlet inert gas (CW,out). The corresponding

Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 4273 Table 7. Tuned Parameters of 3 × 2 Controller parameter IVout set point Tg,out set point controlled variable weights ωc controlled variable soft-constraint weights ωi,soft manipulated variable weights ωm,M manipulated variable steady-state target weights ωm,T um,SST - CW,in um,SST - Tg,in um,SST - m ˘g minimum soft constraints controlled variable 1 controlled variable 2 maximum soft constraint controlled variable 1 controlled variable 2 minimum hard constraint controlled variable 1 controlled variable 2 controlled variable 3 maximum hard constraint controlled variable 1 controlled variable 2 controlled variable 3

value(s) 0.716 dL/g 485.75 K 2.0, 0.1 0, 1, 0 0.1, 0.5, 0.4 0, 0, 0 1.0 × 10-6 mol/L 493 K 1.94 kg/min 0.70 dL/g 473 K 0.73 dL/g 498 K 1.0 × 10-8 mol/L 485 K 0.50 kg/min 1.0 × 10-4 mol/L 503 K 3.50 kg/min

Table 8. SSP Steady-State Parameters

Figure 7. Sketch of the SSP control scheme.

parameter

Table 6. Controller Parameters parameter

value(s)

IVout set point Tg,out set point CW,out set point controlled variable weights ωc controlled variable soft-constraint weights ωi,soft manipulated variable weights ωm,M manipulated variable steady-state target weights ωm,T um,SST - CW,in um,SST - Tg,in um,SST - m ˘g minimum soft constraints controlled variable 1 controlled variable 2 controlled variable 3 maximum soft constraint controlled variable 1 controlled variable 2 controlled variable 3 minimum hard constraint controlled variable 1 controlled variable 2 controlled variable 3 maximum hard constraint controlled variable 1 controlled variable 2 controlled variable 3

0.716 dL/g 485.75 K 4.27 × 10-6 mol/L 1.5, 0.1, 0.1 0, 1, 0 0.3, 0.4, 0.2 0, 0, 0 1.0 × 10-6 mol/L 493 K 1.94 kg/min 0.70 dL/g 473 K 1.0 × 10-7 mol/L 0.73 dL/g 498 K 1.0 × 10-4 mol/L 1.0 × 10-8 mol/L 485 K 0.50 kg/min 4.0 × 10-6 mol/L 503 K 3.50 kg/min

manipulated variables are (1) the inert gas mass flow (m ˘ g), (2) the temperature of the gas flow entering the reactor (Tg,in), and (3) the water concentration in the gas flow entering the reactor (CW,in). The related MPC control scheme, with all of the controlled and manipulated variables listed, is sketched in Figure 7. One of the main MPC controller features is the capability of dealing with MIMO (multiple inputmultiple output) systems both squared problems, as in the aforementioned 3 × 3 case, and nonsquared problems, as in the 3 × 2 configurations also tested in the section. For the sake of simplicity, sensors and valve

residence time polymer feed (m ˘ p) purge N2 rate (m ˘ g) inlet polymer temperature (T0s ) inlet gas temperature (T0g) inlet polymer IV pellet average radius reactor height reactor inner diameter at the entrance hopper discharge height reactor inner diameter at the outlet inlet polymer crystallinity (χ0c ) voidage () inlet concentrations 0 tEG (CtEG ) 0 tTPA (CtTPA ) 0 tVIN (CtVIN ) bEG (C0bEG) bTPA (C0bTPA) bDEG (C0bDEG) TPA (C0TPA) water (C0W) ethylene glycol (C0EG) acetaldehyde (C0AA) gas-phase water (C0W,g) gas-phase ethylene glycol (C0EG,g) 0 gas-phase acetaldehyde (CAA,g ) pressure (P)

value 6h 0.97 kg/min 1.94 kg/min 210 °C 220 °C 0.55 dL/g 0.133 cm 400 cm 40 cm 90 cm 7 cm 0.30 0.40 0.115 mol/L 0.054 mol/L 0 7.087 mol/L 7.147 mol/L 0 0 2.2 × 10-4 mol/L 5.5 × 10-4 mol/L 1.0 × 10-6 mol/L 1.0 × 10-6 mol/L 3.5 × 10-7 mol/L 1.0 × 10-8 mol/L 1 atm

dynamics are disregarded in the analyzed examples, whose control parameters are listed in Tables 6 and 7. The simulated controllers were tested with disturbances in the initial IV of the polymer mixture entering the reactor, the steady-state conditions for which appear in Table 8. Other kinds of disturbances in the reactor inputs have similar effects on the output, and therefore, in the following discussion, tests will performed only by IV input variations. The IV input value disturbance was simulated as an unmeasured disturbance, thus becoming evident to the controller only when the polymer exits the SSP reactor. Therefore, the MPC controller must demonstrate its effectiveness in feedback correction. First, several different configurations were tested to

4274 Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004

Figure 8. Comparison between different control configurations.

Figure 9. Comparison of the IV output for different sampling times of the tuned 3 × 2 control system.

determine the best. In Figure 8 are reported the simulated responses to a disturbance in the entering polymer IV equal to +0.02 dL/g. The squared 3 × 3 scheme (CW,in, Tg,in, and m ˘ g are manipulated to control

IVout, Tg,out, and CW,out) is compared with two nonsquared 3 × 2 control schemes, the first without CW,out control and the second without IVout control. The results of the simulations (all performed with nHC ) 2, nHP ) 5, and a sampling time of 15 min) show that the IVout control is needed; otherwise, the effect of the disturbance can be only partially readsorbed. Moreover, the 3 × 3 scheme performances are slightly worse than those of the 3 × 2 scheme without CW,out control (the third controlled variable, CW,out, seems to not be useful in assisting IV control), and therefore, only this configuration is considered and further examined in the following. Once the control scheme was selected, it had to be tuned to obtain better control performances. The main parameters are nHC, nHP, and ts, and their optimal values can be selected through a sensitivity analysis that, for the sake of brevity, is partially omitted here. The chosen values are nHP ) 3, nHC ) 2, and ts ) 5 or 10 min. For a disturbance on the entering polymer IV equal to +0.02 dL/g, Figure 9 compares the tuned controller performances. Both sampling times guarantee an adequate response to the disturbance, reaching a stable state after a time equal to one nominal reactor residence time, even if the longer sampling time (ts ) 10 min) has a slightly more regular action. This is also highlighted in Figure 10, where the corresponding manipulated variable profiles for the two cases are also compared. From a practical point of view, the application of a shorter sampling time is always desired because it guarantees a more frequent realignment of the control model to the real plant data, being capable of a more efficient response to any sort of unmeasured disturbance. Figure 11 shows the effect of a variation of the inlet polymer IV by -0.04 dL/g, which is in the opposite direction and of greater amplitude than the previously examined disturbance, thus testing the reliability of the implemented control system. It can be noted that the controller acts by increasing the gas temperature to reject the disturbance, thereby enhancing the product

Figure 10. Comparison between profiles of manipulated variables for two different sampling times.

Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 4275

Figure 11. Profiles of controlled and manipulated variables in the test with a disturbance in IVin of -0.04 dL/g.

quality. However, at the new steady state, the outlet gas temperature is subject to a very slight change because of the decreased gas flow. This confirms the capability of a multivariable control scheme in managing each single manipulated variable to optimize a common objective function subject to all of the process constraints. Finally, Figure 12 shows a simulated servo problem in which the IV set point is increased from 0.716 to 0.735 dL/g. Also in this case, the controller action is effective, as the new steady state is achieved in less than onehalf of the nominal residence time, and it is compatible with the temperature control, which is obtained by a decrease of the inlet gas flow, a minimum increase of the inlet gas temperature, and a significant reduction of the water concentration in the inert gas flow. Let us conclude with some comments about the computing time needed by the MPC algorithm to define the optimal control action. As already mentioned, the simulated control action can be implemented on a real-time basis only if the computations are completed in a time less than the corresponding sampling time. An MPC controller similar to the one simulated here is unusual: the model is a complex nonlinear system of partial differential equations, and it is obviously time-expensive. Nevertheless, for a process such as SSP, with a high residence time compatible with a control action of 5 or

10 min, even a model with a high CPU time can be adopted. In Figure 13, the CPU times of a general transient are shown, as measured on an Intel Pentium 4 2.4GHz processor, together with the number of calls of the model for each time step. On the graph are also indicated the three critical thresholds of 5, 10, and 15 min. The lower threshold is generally respected, with a limited number of violations occurring. This lack of robustness of the controller, due to the high number of model calls needed by the optimizer, can be solved with further work on the optimizer settings (a robust simplex method was utilized for the simulations). In fact, by adding a limit to the number of model calls, the achievement of a good control action can be guaranteed even if this leads to a suboptimal solution. More refined and modern dynamic optimization techniques, such as the multiple shooting and simultaneous approach as reviewed by Biegler et al.,9 can strongly reduce the CPU time needed, solving the model only once at optimality. However, this is beyond scope of this work, and it is matter of current research. Conclusions This study has succeeded by including a firstprinciples large-scale model into an effective MPC

4276 Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004

Figure 12. Profiles of controlled and manipulated variables for a variation of the IVout set point.

Figure 13. Typical CPU times and numbers of calls to the model for all control actions in a standard test.

algorithm by appropriate model reduction. The proposal of a multivariable control scheme for the SSP process is another original part of this study. The choice of gas temperature and gas water concentration as manipulated variables has revealed interesting aspects, highlighting the great effect of the former on process performance and the possible importance of the latter

in polymer quality control. The servo and regulatory performance of the controller was tested on a detailed SSP reactor simulator, demonstrating that it is effective and of practical value. The CPU time analysis suggests that a real-time application is feasible, although some optimization algorithm constraint are needed to improve the controller robustness.

Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 4277

The aim of this work has been to demonstrate that a first-principles model, even if appropriately sized, can be of real value for control purposes. In other words, a reduced nonlinear model, integrated into the control algorithm, is the topic subject of our work, and hopefully, in the future, it will be compared with a more traditional approach (step response) based on industrial field data. Acknowledgment The authors are grateful to Professor Kevin Yao for his courtesy.

r ) radial distance from the center of a spherical polymer particle, cm R ) ideal gas constant, L‚atm/(mol‚K) Re ) Reynolds number Rj ) rate of reaction j, mol/(L‚min) t ) Time, min Ts ) temperature of the polymer, K Tg ) temperature of the gas phase, K TR ) superficial temperature of the polymer particle, K TW ) temperature of the reactor wall, K vg, vp ) gas and polymer flow rates, respectively, cm/min yj ) mole fraction of component j in the gas phase z ) axial distance, measured from the top of the reactor, cm

Nomenclature

Greek Letters

ac ) cross-sectional area of reactor, cm2 C ˜ j ) average concentration of species within the solid particles, mol/L Cj,g ) concentration of volatile species in the gas phase, mol/L κp, Cp,g ) heat capacities of the polymer and purge gas, respectively, kcal/(kg‚K) Cpv,j ) heat capacity of volatile vapors, kcal/(kg‚K) Db ) dispersion coefficient, cm2/min dp ) particle diameter, cm Dg,j ) diffusivity of component j in the gas phase, cm2/min Dp,j ) diffusivity of component j in the solid phase, cm2/ min H ) Henry’s law coefficient, atm H ) reactor height, cm hg ) gas-side heat-transfer coefficient, kcal/(cm2‚K‚min) hp ) polymer-side heat-transfer coefficient, kcal/(cm2‚K‚ min) ho ) overall heat-transfer coefficient, kcal/(cm2‚K‚min) hw ) gas-reactor wall heat-transfer coefficient, kcal/(cm2‚ K‚min) ∆HCR ) heat of crystallization, kcal/kg ∆Hev,j ) heat of vaporization of component j, kcal/mol jD ) Colburn factor for mass transfer jH ) Colburn factor for heat transfer kj ) kinetic rate constant for reaction j, (L/mol)/min Kj ) equilibrium constant for reaction j kc ) kinetic constant of crystallization, 1/min kg,j ) mass-transfer coefficient of component j, mol‚cm/(L‚ atm‚min) ko,j ) overall mass-transfer coefficient of component j, mol‚ cm/(L‚atm‚min) kp,j ) Polymer-side mass-transfer coefficients, mol‚cm/(L‚ atm‚min) m ˘ g, m ˘ p ) mass flow rates of gas and polymer, respectively, kg/s Mj ) molecular weight of component j, kg/mol P ) reactor pressure, atm pj,g ) partial pressure of component j in the gas phase, atm p˜ j ) vapor pressure of component j in equilibrium with the polymer average concentration, atm

χc ) cystalline fraction χmax ) maximum crystalline fraction  ) voidage Φj ) volumetric fraction of component j in the polymer κg, κp ) thermal conductivities of the gas phase and polymer, respectively, kcal/(min‚cm‚K) ηg ) viscosity of the gas phase, kg/(min‚cm) Fg, Fp ) densities of the gas phase and polymer, respectively, kg/cm3

Literature Cited (1) Morari, M.; Lee, J. H.; Model Predictive Control: Past, Present and Future. Commun. Chem. Eng. 1999, 23, 667. (2) Krishnan, A.; Kosanovich, K. A.; DeWitt, M. R.; Creech, M. B. Proceedings of the American Control Conference; IEEE: Philadelphia, PA, June 24-26, 1998; pp 3386-3390. (3) 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. (4) Blom, J. G.; Trompert, R. A.; Verwer, J. G. Algorithm 758: VLUGR2: A Vectorizable Adaptive-Grid Solver for PDEs in 2D. ACM Trans. Math. Software 1996, 22, 302. (5) Yao, K. Z.; McAuley, K. B.; Marchildon, E. K. Simulation of continuous solid-phase polymerization of nylon 6,6 (III). J. Appl. Polym. Sci. 2003, 89, 3701. (6) Carslaw, H. S.; Jaeger, J. C. Conduction of Heat in Solid, 2nd ed.; Oxford University Press: Oxford, U.K., 1959. (7) Mallon, F. K.; Ray, W. H. Modeling of solid-state polycondensation (II), Reactor design issues. J. Appl. Polym. Sci. 1998, 69, 1775. (8) Leffew, K. W.; Yerrapragada, S. S.; Deshpande, P. B. Six Sigma and Solid-State Polymerization. Chem. Eng. Commun. 2001, 188, 109. (9) Biegler, L. T.; Cervantes, A. M.; Wachter, A. Advances in Simultaneous Strategies for Dynamic Process Optimization. Chem. Eng. Sci., 2002, 57, 575.

Received for review December 18, 2003 Revised manuscript received May 5, 2004 Accepted May 14, 2004 IE034324I