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 rate constant

the present and the future of chemical process control. reaction forward reverse
Its fundamentals reached popularity in the early 1970s. 1 TEG + tTPA ‚ bEG + bTPA + W k1 k1/K1
Currently, the availability of powerful computers has 2 TEG + tEG ‚ bEG + EG k2 k2/K2
permitted the application of such model predictive
controllers to real systems, and in the future, an ever- by the long dead time of the plug-flow reactor and by
increasing interest in this field is expected. More details the weak effects of the manipulated variables, which
on the history and trends of this technology can be found also cause a low efficiency for traditional PID control-
in Morari and Lee.1 The MPC strategy adopts a process lers.
model, which can be the same as that used for process The next section deals with the analytical reduction
design, with a high degree of detail, to predict the future of the detailed model, which was necessary to obtain a
behavior of the controlled system and to perform the practical MPC controller for the reactor. The original
calculation of the optimal set of control actions, dictated model of the reactor (see part I of this work3), consisting
by the minimization of a common objective function. of a system of 16 time-dependent partial differential
This function is structured to increase in value when equations in two spatial coordinates, was simplified by
the system moves away from its steady state, in which contracting the radial direction into a unique grid point
case the controlled variables differ significantly from and substituting the corresponding derivatives with
their set points. Bounds and constraints can also be simplified linear relationships accounting for the mass
taken into account for specific requirements on both and thermal diffusion into the solid. The model was
manipulated and controlled variables. This feature solved with the VLUGR2 library (by Blom et al.4), and
make MPC capable of successfully treating multivari- the resulting savings of CPU time was around 2 orders
able systems with critical characteristics such as strong of magnitude, dropping from 1 h needed by the complete
nonlinear behavior, inverse responses, long dead times, model to a few seconds required by the reduced one.
or marked interactions among variables. The degree of In a subsequent step, this model was incorporated
detail adopted for the translation of the real process into a model predictive control structure, and the overall
behavior into mathematical terms synthesizes the chal- control loop was implemented, through the use of the
lenge represented by MPC, which does not impose any detailed model as a “virtual plant”, to demonstrate the
limits on its model structure, thus improving on many effectiveness of both the designed control scheme and
traditional control techniques that, today, might be the model reduction technique.
rendered obsolete. A practical restriction derives only
from the ability to obtain detailed models with a Model Reduction
computing time that satisfies the limits needed for the The model simplification was realized following simi-
control to be implemented in a real-time system. The lar assumptions as highlighted by Yao et al.5 in the third
solid-state polymerization (SSP) process is exactly such part of their work on SSP for the nylon process.
a system, as already observed by Krishnan et al.,2 for Substantially, it consists of the removal of the radial
which the choice of model size becomes critical for MPC dimension by substituting it with a single node and
applications but, at the same time, conventional control translating the description of mass and heat diffusion
shows its weakness. Such poor controllability is caused within the polymer into simple linear relationships.
Moreover, the kinetic scheme was simplified by taking
* To whom correspondence should be addressed. E-mail: into account only the two main reactions (esterification
[Link]@[Link]. Fax: +39.02.7063.8173. 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) ∂C̃j 3

Terms ) -k0,j(p̃j - pj,g) +
∂t R
Generation Rate Terms Gj(t)
GEG ) R2 GtEG ) -R1 - 2R2 m̆p ∂C̃j ∂2C̃j
GW ) R 1 GtTPA ) -R1 Gj(t) - + Db 2 (1)
(1 - )acFp ∂z ∂z
Reaction Rate Terms Rj(t)
CbEG The parameter k0,j can be considered a global mass-
R1 ) k1CtEGCtTPA - 2(k1/K1)CWCbTPA transfer coefficient of component j, and it depends on
CtEG + CbEG
two primary coefficients: kg,j corresponding to the
R2 ) k2CtEGCtEG - 4(k2/K2)CEGCbEG interphase mass transfer and kp,j related to the internal
diffusion resistance. The former can easily be calculated
Table 3. Legend of Abbreviations for Molecular Species 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

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

where Hj is the Henry’s law pseudo-equilibrium con-

stant calculated as the ratio between p̃j, the equilibrium
vapor pressure obtained from the Flory-Huggins equa-
tion, and the volume fraction φj of component j in the
solid particle

p̃j
Hj ) (3)
φj

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 m̆p ∂T̃s h0 3

)- + (T - T̃s) +
∂t Fp(1 - )ac ∂z FpCp g R
corresponding reaction rates reported in Table 2. There- ∂2T̃s Gc(t)(-∆HCR) ∆Hev,W 3
fore only nine of the original 16 state variables involved Db + - k0,W (p̃ - pW,g) -
∂z 2 Cp,p FpCp W R
are still needed, namely, the ethylene glycol (EG), water
(W), EG end group (tEG), and terephthalic acid end ∆Hev,EG 3
k0,EG (p̃EG - pEG,g) (4)
group (tTPA) concentrations within the solid particle; FpCp R
the EG and W concentrations in the gas phase; the solid
and gas temperatures, Ts and Tg, respectively; and the The thermal exchange between solid and gas can be
crystallinity χc. The concentrations of the TPA and EG evaluated as a linear function of the driving force and,
repeat units (bTPA and bEG, respectively) can be using the lumped coefficient h0, obtained in a complete
obtained algebraically, and they do not lead to further analogy with coefficient k0,j for the mass balance.
differential equations. (See Table 3 for the abbreviations Similarly to eq 2, the interphase heat-transfer coefficient
for the molecular species considered in this work.) hg can be coupled with the internal heat conduction
Moreover, for the evaluation of the gas-phase proper- coefficient hp, which primarily depends on the thermal
ties, the hypothesis of a pure nitrogen gas flow was conductivity κp, giving the global thermal coefficient as
employed because the concentrations of byproducts are a sum of the two resistances
negligible, being always below 0.1%.
The model reduction method adopted here, suggested 1 1 1
) + (5)
by Yao et al.,5 substitutes Fick’s law with a simple linear ho hp hg
correlation that defines the interphase and internal
transfers by coupling their effects in a unique constant Obviously, the mass and thermal balances for the gas
coefficient. A reduction in the radial coordinate implies phase are not affected by the model reduction because
that the computation of the concentration and temper- their structure does not depend on the particle radial
ature gradients within the solid particle is no longer coordinate. The only observable changes are those
needed. Obviously, all state variables y now have the related to the new transfer coefficients adopted, k0,j and
property of being averaged, ỹ, along the particle radius, h0, which must also appear in gas-phase equations.
corresponding to a general representation given as ỹ ) Consequently, the terms reported in the driving force
ỹ(z,t) rather than y ) y(r,z,t). Therefore, the material for the mass and heat flux must be p̃j and T̃s, implicitly
balances of the diffusing compounds j (i.e., EG and W) referring to the average values of the flux into the
in a plug-flow reactor scheme assume the following particle, and not to pj,R and Ts,R which are related to
structure the surface values.
 Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 4269

Finally, the mass and thermal balances for the gas Tin - TR π2
phase assume the following form ) (12)
T̃ - TR ∞
1
∂Cj,g m̆g ∂Cj,g ∂2Cj,g Dg,j ∂ac ∂Cj,g
6 ∑
n)1n
2
e -Rn2π2t/R2

) + Dg,j + +
∂t Fgac ∂z ∂z2 ac ∂z ∂z
3(1 - ) which, once substituted into eq 11, gives rise to the
k0,j(p̃j - pj,g) (6) following new relationship
R
∞
∂2Tg ∑ e-Rn π t/R
2 2 2
∂Tg m̆g ∂Tg κg κg dac ∂Tg κpπ2
) + + - n)1
∂t Fgac ∂z FgCp,g ∂z2 FgCp,gac dz ∂z hp ) (13)
∞
1
2vxacπ 3R ∑ e -Rn2π2t/R2

hw(Tg - TW) +
FgCp,gac j)EG,W
∑
k0,j(p̃j - pj,g) n)1n
2

(T̃s - Tg)3(1 - ) 3(1 - ) This last equation allows for the evaluation of the
Cpv,j + h0(T̃s - Tg) (7) polymer-side heat-transfer coefficient that has the
FgCp,gR FgCp,gR 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
The coefficients hp and kp,j previously defined in eqs thermal properties.
2 and 5 and representing the resistances to mass and In their work, Yao et al.5 also suggested the conve-
heat diffusion into the solid particle, need some further nient introduction of a time-invariant value of this
clarification about their evaluation. The method for lumped parameter equal to the infinite-time heat-
estimating hp, reported below, was proposed by Yao et transfer coefficient. Assuming such a hypothesis for hp,
al.,5 and for sake of brevity, the description of the eq 9 can be solved analytically to obtain
corresponding procedure for the mass-transfer coef-
ficients is omitted because it is completely analogous
to that for hp.
T̃ ) TR + (Tin - TR)e-3hpt/FpCpR (14)
The main idea is to adapt the analytical solution of
Results given by the analytical solution and the simpli-
the heat conduction problem of a sphere with uniform
fied one can be compared and an adjusting factor ch can
initial temperature Tin and with its outer surface held
at a constant temperature TR. Carslaw and Jaeger6 be introduced as hp ) chh∞p .
provided the following solution for such a problem In this study, this correction factor has been used
instead as a fitting parameter between models, as
suggested by Yao et al.5
6(TR - Tin) ∞
1 The lumped parameter kp,j for mass transfer can be
∑ e-Rn π t/R
2 2 2
T̃ ) TR - (8) determined following exactly the same method. Obvi-
2 2
π n)1 n ously, the value of kp,j depends on the diffusivity of
component j and on the pellet dimensions. Therefore,
where R ) κp/FpCp is the thermal diffusivity inside the in complete analogy with hp, when t f 0, kp,j f ∞,
sphere. Therefore, defining the rate of heat transfer whereas when t f ∞, kp,j f Dw,jFpπ2/3R.
from the particle body to its outer surface through a Neglecting the time dependency and choosing the
linear correlation, the particle heat balance can be constant infinite-time value of kp,j ) Dw,jFpπ2/3R, the
formulated as follows dynamic material balance of component j can be solved
analytically to give
4 dT̃
- πR3FpCp ) 4πR2hp(T̃ - TR) (9)
3 dt C̃j ) Cj,R + (Cj,in - Cj,R)e-3kp,jt/RFp (15)

where hp is the laminar transfer coefficient representing Even in this case, an adjustment factor of ck,j can been
the polymer-side resistance. Differentiating the expres- adopted as a correction between the simplified and
sion of T̃ in eq 8 gives analytical solutions or extended as a fitting parameter
between the detailed and simplified models (kp,j ) ck,j
∞
dT̃ 6(TR - Tin)R ∞ kp,j ).
∑ e-Rn π t/R
2 2 2
) (10) The system equations and initial/boundary conditions
dt R 2
n)1 of the resulting reduced model are summarized in
Tables 4 and 5, respectively. The system of time-
Then, replacing dT̃/dt in eq 9 and solving for hp, one dependent one-dimensional partial differential equa-
obtains tions was numerically solved through the same library,
VLUGR2,4 as adopted for the detailed model. Although
this numerical solver is designed for time-dependent
(Tin - TR) ∞ two-dimensional partial differential equations, it was
∑
e-Rn π t/R
2 2 2
hp ) 2κp (11) possible to solve the one-dimensional system by impos-
R(T̃ - TR)n)1 ing a fictitious coordinate (replacing the radial one)
constituted by a unique grid node. The CPU time saved
A further simplification can be realized rearranging eq is remarkable: the reduced model requires a computing
8 as 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 3 m̆p ∂C̃j ∂2C̃j (19)

) -k0,j(p̃j - pj,g) + Gj(t) - + Db 2
∂t R (1 - )acFp ∂z ∂z

∂C̃j m̆p ∂C̃j ∂2C̃j (20)

) Gj(t) - + Db 2
∂t (1 - )acFp ∂z ∂z

∂χc m̆p ∂χc ∂2χc (21)

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

∂Cj,g m̆g ∂Cj,g ∂2Cj,g Dg,j ∂ac ∂Cj,g 3(1 - ) (22)

) + Dg,j + + k0,j(p̃j - pj,g)
∂t Fgac ∂z ∂z2 ac ∂z ∂z R

∂T̃s m̆p ∂T̃s h0 3 ∂2T̃s Gc(t)(-∆HCR)

)- + (Tg - T̃s) + Db 2 + -
∂t Fp(1 - )ac ∂z FpCp R ∂z Cp,p
(23)
∆Hev,W 3 ∆Hev,EG 3
k0,W (p̃ - pW,g) - k0,EG (p̃EG - pEG,g)
FpCp W R FpCp R

∂Tg m̆g ∂Tg κg ∂2Tg κg dac ∂Tg 2vxacπ

) + + - hw(Tg - TW) +
∂t Fgac ∂z FgCp,g ∂z2 FgCp,gac dz ∂z FgCp,gac
(24)
j)EG,W (T̃s - Tg)3(1 - ) 3(1 - )
∑ j
k0,j(p̃j - pj,g)Cpv,j
FgCp,gR
+ h0(T̃s - Tg)
FgCp,gR

Table 5. Initial and Boundary Conditions of the

Reduced Model
Cj|t)0 ) C0j ∀ z ∈ {0 e z e H}
Cj|z)H ) C0j

∂Cj
∂z| z)0
)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 -
Dg,jac ∂Cj,g
m̆g ∂z | z)H-

∂Cj,g
∂z |
z)0
j ) EG, W, AA
)0

χc|t)0 ) χ0c ∀ z ∈ {0 e z e H}
χc|z)0 ) χ0c Figure 1. IV profiles along the reactor height calculated with
the detailed (solid lines) and reduced (dashed lines) models at
∂χc
∂z | z)H
)0
different reaction times.

information returned by the reduced model. Therefore,

Ts|t)0 ) T0s ∀ z ∈ {0 e z e H}
Ts|z)0 ) T0s some checks are needed to verify the feasibility of the
derived tool.
∂Ts
∂z | z)H
)0
To carry out a comparison between model perfor-
mances, it becomes necessary to introduce the molecular
average number and the related intrinsic viscosity, IV
Tg|t)0 ) T0g ∀ z ∈ {0 e z e H}

∂Tg
∂z | z)0
)0 Xn )
CtEG + CtTPA + CbEG + CbTPA
CtEG + CtTPA
(16)

Tg|z)H- ) T0g|z)H -
κj,gac ∂Tg
m̆gCp,g ∂z | z)H-
IV ) 2.1 × 10-4(192.17 × Xn)0.82 (17)

In Figure 1 is shown a comparison between the

model, mainly because of the simplifications that lighten dynamic evolutions of the IV along the reactor height
the numerical iterations. Conversely, the reduction in when calculated with the detailed model (solid lines) and
terms of the number of state variables and the other with the reduced model (dashed lines). A good agree-
assumption adopted give rise to minor amounts of ment 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 Figure 4. EG concentration profiles at the reactor outlet predicted
and detailed models. by the reduced and detailed models.

Figure 3. Particle temperature profiles along the reactor height Figure 5. W concentration profiles at the reactor outlet predicted
predicted by the reduced and detailed models at different reaction by the reduced and detailed models.
times.
given that a feedback error corrector term is always
the adjustment factors, ck,j, to 3.5 for ethylene glycol and incorporated in the controller structure (for details, see,
to 1.5 for water. These relatively high values are mainly for example, Morari and Lee1).
due to the reduction of the kinetic scheme.
Moreover, it is important to verify the analogy of the
Model Predictive Control Configuration
dynamic behaviors for the two models and, hence, to
determine whether the variable profiles at the reactor From its original theory, the MPC algorithm calcu-
outlet as a function of time can provide reasonably the lates, during each iterative cycle repeated after a time
same information. Figure 2 shows a comparison of the step equal to the sampling time ts, the optimal profiles
outlet IV predicted by the models, and the presented of the manipulated variables along a predefined tem-
curves confirm a very good agreement. In addition, poral horizon, the so-called control horizon. The set of
Figure 3 compares the particle temperature profiles best possible actions is evaluated by minimizing an
obtained from the two models. The corresponding curves objective function, which is determined mainly by the
essentially overlap, even if the corresponding lumped difference between the output of the model and the
parameter hp is not adjusted by any fitting parameter corresponding imposed set points, predicted along a
ch. prediction horizon. Usually, these two temporal horizons
The last comparison, illustrated in Figures 4 and 5, are represented as multiples of the sampling time, and
is of the two diffusing components EG and W. Their their lengths can be indicated with nHC and nHP,
profiles at the reactor outlet, as determined by the respectively. A control horizon of nHC means that the
detailed and reduced models, show a good agreement objective function will be minimized by a control action
for the predicted dynamics even though the water at each time step between 1 and nHC with a total
concentration has a significant error at the steady state number of degrees of freedom equal to nHC times the
because of the simplification of the kinetic scheme. number of controlled variables. The prediction horizon
However, with the aim of developing an MPC controller, is longer than, or at least equal to, the control horizon,
the presence of a tolerable deviation between the and therefore, in any control actions between nHC and
predictions of the simplified model and the detailed nHP, the manipulated variables will maintain the last
model (representing the real plant) can be accepted, 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

{ [ ]}
yc(i) - yc,set(i)
k+nHP nc 2

F) ∑ c)1
∑ ωc yc,set(i)
+

∑ (∑{ [ ]
i)k+1
k+nHC-1 nm um(l) - um(l - 1) 2
ωm,M +
um(l - 1)

[ ] })
l)k m)1
um(l) - um,SST(l) 2
ωm,T (18)
um,SST(l)
Figure 6. Schematic illustration of an MPC computing cycle.
The structure is a quadratic function in terms of error
procedure, the first control action, among the nHC for each controlled variable yc with respect to the
mentioned above, is implemented and maintained the corresponding set points yc,set along the prediction
same throughout the sampling time, after which a new horizon hp and also in terms of the deviation for each of
measurement becomes available, and the calculation manipulated variables um along the control horizon hc.
cycle can restart. During the optimization computation, These last terms can be determined in two different
a feedback correction for the model prediction must also forms. The first is the difference between two subse-
be taken into account. In fact, the error between the quent values of the manipulated variable um(l) and um-
controlled variables computed by the model and the real (l - 1), which is include to avoid control actions that
measurements is estimated at the beginning of the are too strong. The second is the difference between the
evaluation cycle, and this error is employed to correct actual value of the manipulated variable um and the
all future model predictions along the new optimization corresponding value of its desired steady-state condition
step. Such a strategy avoids a model prediction of the um,SST (steady-state target). Generally, the importance
process behavior that is not in agreement with the real of each single contribution of the objective function is
state of the system. Therefore, this correction method managed through the weights ωc, ωm,M, and ωm,T.
allows one to reject and control unmeasured distur-
SSP Control Scheme
bances.
To better illustrate the MPC computation procedure Defining which variables have to be controlled and
mentioned above, Figure 6 shows schematically the which have to be manipulated is clearly the first step
main steps, which can be summarized as follows: of any control system setup. In the previous paper (see
part I of this work3), the analysis pointed out that only
(1) At time t0 (iHC ) 1), the measurement y0 of the
the purge gas temperature and byproduct concentra-
current process status becomes available.
tions can be used efficiently to control the process
(2) The optimization, constrained or not, is executed outputs. The more obvious controlled variable is the
in real time, minimizing the objective function F ) intrinsic viscosity IV, which is the primary quality-
F(y,u) along the prediction horizon nHP and supplying control parameter related to the average polymer chain
the optimal profile u j (t) of manipulated variables along length, but also a temperature control is needed because
the control horizon nHC. this parameter is limited by a maximum value dictated
(3) The first control action, u j 0, of the optimal profile, by the sticking polymer temperature. In the literature,
determined at step 2, is implemented on the process, similar conclusions can easily be found. In fact, in their
and the cycle restarts at t0 ) t0 + ts returning to step 1. simulations, Mallon and Ray7 showed that an increased
It is clear that, to have a feasible control solution, the concentration of condensates in the gas phase can
optimization problem must be solvable in a computing influence significantly the product quality, whereas
time less than or equal to the sampling time ts. Leffew et al.8 emphasized the role of the purge gas flow
Consequently, this requirement can also be applied to and its temperature to better control the solid-bed
the model that has to be fast enough to allow the temperature profile and consequently the polymer
optimization procedure to be performed within the intrinsic viscosity.
feasible time limit. The core of the MPC control phi- Direct control of the IV at reactor outlet requires a
losophy is the objective function, with regard to which continuous measurement of this variable. For the
some interesting features can be highlighted. In fact, similar nylon SSP process, Krishnan et al.2 adopted an
the MPC control actions are obtained not by evaluating industrial scheme in which the polymer molecular
only the absolute error (y - yset), as in conventional weight, and consequently the IV, is inferred through the
control, but rather by adopting a more detailed mini- pressure difference between the discharging point of the
mization objective function, which can be structured by extruder (located at reactor outlet) and the correspond-
taking into account several constraints and different ing spinning manifold. Therefore, because the two
operating policies. In a general formulation, not only processes are very similar, this type of measurement is
the absolute error on controlled variables can be penal- also assumed to be available in this work. Moreover,
izing but also the excessive deviation of the correspond- Krishnan et al.2 adopted an additional outlet water
ing manipulated variables from their steady-state val- concentration control.
ues. Soft and hard operating constraints can also be In summary, all of the possible controlled variables
imposed, allowing the real systems to be fully managed can be listed as follows: (1) the intrinsic viscosity, IV,
while the predictive model takes into account the of PET at reactor outlet (IVout); (2) the temperature of
coupling between variables and actions, thereby making the outlet inert gas (Tg,out); and (3) the water concentra-
easier an effective multivariable control structure. tion 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 value(s)
IVout set point 0.716 dL/g
Tg,out set point 485.75 K
controlled variable weights ωc 2.0, 0.1
controlled variable soft-constraint 0, 1, 0
weights ωi,soft
manipulated variable weights ωm,M 0.1, 0.5, 0.4
manipulated variable steady-state 0, 0, 0
target weights ωm,T
um,SST - CW,in 1.0 × 10-6 mol/L
um,SST - Tg,in 493 K
um,SST - m̆g 1.94 kg/min
minimum soft constraints
controlled variable 1 0.70 dL/g
controlled variable 2 473 K
maximum soft constraint
controlled variable 1 0.73 dL/g
controlled variable 2 498 K
minimum hard constraint
controlled variable 1 1.0 × 10-8 mol/L
controlled variable 2 485 K
controlled variable 3 0.50 kg/min
maximum hard constraint
controlled variable 1 1.0 × 10-4 mol/L
controlled variable 2 503 K
controlled variable 3 3.50 kg/min

Figure 7. Sketch of the SSP control scheme. Table 8. SSP Steady-State Parameters
parameter value
Table 6. Controller Parameters
residence time 6h
parameter value(s) polymer feed (m̆p) 0.97 kg/min
IVout set point 0.716 dL/g purge N2 rate (m̆g) 1.94 kg/min
Tg,out set point 485.75 K inlet polymer temperature (T0s ) 210 °C
CW,out set point 4.27 × 10-6 mol/L inlet gas temperature (T0g) 220 °C
controlled variable weights ωc 1.5, 0.1, 0.1 inlet polymer IV 0.55 dL/g
controlled variable soft-constraint 0, 1, 0 pellet average radius 0.133 cm
weights ωi,soft reactor height 400 cm
manipulated variable weights ωm,M 0.3, 0.4, 0.2 reactor inner diameter at the entrance 40 cm
manipulated variable steady-state 0, 0, 0 hopper discharge height 90 cm
target weights ωm,T reactor inner diameter at the outlet 7 cm
um,SST - CW,in 1.0 × 10-6 mol/L inlet polymer crystallinity (χ0c ) 0.30
um,SST - Tg,in 493 K voidage () 0.40
um,SST - m̆g 1.94 kg/min inlet concentrations
minimum soft constraints 0
tEG (CtEG ) 0.115 mol/L
controlled variable 1 0.70 dL/g 0
tTPA (CtTPA ) 0.054 mol/L
controlled variable 2 473 K 0
0
tVIN (CtVIN )
controlled variable 3 1.0 × 10-7 mol/L
maximum soft constraint bEG (C0bEG) 7.087 mol/L
controlled variable 1 0.73 dL/g bTPA (C0bTPA) 7.147 mol/L
controlled variable 2 498 K bDEG (C0bDEG) 0
controlled variable 3 1.0 × 10-4 mol/L TPA (C0TPA) 0
minimum hard constraint water (C0W) 2.2 × 10-4 mol/L
controlled variable 1 1.0 × 10-8 mol/L ethylene glycol (C0EG) 5.5 × 10-4 mol/L
controlled variable 2 485 K acetaldehyde (C0AA) 1.0 × 10-6 mol/L
controlled variable 3 0.50 kg/min gas-phase water (C0W,g) 1.0 × 10-6 mol/L
maximum hard constraint 3.5 × 10-7 mol/L
gas-phase ethylene glycol (C0EG,g)
controlled variable 1 4.0 × 10-6 mol/L
gas-phase acetaldehyde (CAA,g0
) 1.0 × 10-8 mol/L
controlled variable 2 503 K
controlled variable 3 3.50 kg/min pressure (P) 1 atm

manipulated variables are (1) the inert gas mass flow dynamics are disregarded in the analyzed examples,
(m̆g), (2) the temperature of the gas flow entering the whose control parameters are listed in Tables 6 and 7.
The simulated controllers were tested with distur-
reactor (Tg,in), and (3) the water concentration in the
bances in the initial IV of the polymer mixture entering
gas flow entering the reactor (CW,in).
the reactor, the steady-state conditions for which appear
The related MPC control scheme, with all of the in Table 8. Other kinds of disturbances in the reactor
controlled and manipulated variables listed, is sketched inputs have similar effects on the output, and therefore,
in Figure 7. in the following discussion, tests will performed only by
One of the main MPC controller features is the IV input variations. The IV input value disturbance was
capability of dealing with MIMO (multiple input- simulated as an unmeasured disturbance, thus becom-
multiple output) systems both squared problems, as in ing evident to the controller only when the polymer exits
the aforementioned 3 × 3 case, and nonsquared prob- the SSP reactor. Therefore, the MPC controller must
lems, as in the 3 × 2 configurations also tested in the demonstrate its effectiveness in feedback correction.
section. For the sake of simplicity, sensors and valve First, several different configurations were tested to
4274 Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004

IVout, Tg,out, and CW,out) is compared with two non-

squared 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. More-
over, 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
Figure 8. Comparison between different control configurations. 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 applica-
tion 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 9. Comparison of the IV output for different sampling Figure 11 shows the effect of a variation of the inlet
times of the tuned 3 × 2 control system. polymer IV by -0.04 dL/g, which is in the opposite
direction and of greater amplitude than the previously
determine the best. In Figure 8 are reported the examined disturbance, thus testing the reliability of the
simulated responses to a disturbance in the entering implemented control system. It can be noted that the
polymer IV equal to +0.02 dL/g. The squared 3 × 3 controller acts by increasing the gas temperature to
scheme (CW,in, Tg,in, and m̆g are manipulated to control 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 10 min, even a model with a high CPU time can be
gas temperature is subject to a very slight change be- adopted. In Figure 13, the CPU times of a general trans-
cause of the decreased gas flow. This confirms the cap- ient are shown, as measured on an Intel Pentium 4 2.4-
ability of a multivariable control scheme in managing GHz processor, together with the number of calls of the
each single manipulated variable to optimize a common model for each time step. On the graph are also indica-
objective function subject to all of the process constraints. ted the three critical thresholds of 5, 10, and 15 min.
Finally, Figure 12 shows a simulated servo problem in The lower threshold is generally respected, with a limit-
which the IV set point is increased from 0.716 to 0.735 ed number of violations occurring. This lack of robust-
dL/g. Also in this case, the controller action is effective, ness of the controller, due to the high number of model
as the new steady state is achieved in less than one- calls needed by the optimizer, can be solved with further
half of the nominal residence time, and it is compatible work on the optimizer settings (a robust simplex method
with the temperature control, which is obtained by a was utilized for the simulations). In fact, by adding a
decrease of the inlet gas flow, a minimum increase of limit to the number of model calls, the achievement of
the inlet gas temperature, and a significant reduction a good control action can be guaranteed even if this leads
of the water concentration in the inert gas flow. to a suboptimal solution. More refined and modern dy-
Let us conclude with some comments about the com- namic optimization techniques, such as the multiple
puting time needed by the MPC algorithm to define the shooting and simultaneous approach as reviewed by
optimal control action. As already mentioned, the simu- Biegler et al.,9 can strongly reduce the CPU time
lated control action can be implemented on a real-time needed, solving the model only once at optimality.
basis only if the computations are completed in a time However, this is beyond scope of this work, and it is
less than the corresponding sampling time. An MPC matter of current research.
controller similar to the one simulated here is unusual:
the model is a complex nonlinear system of partial dif- Conclusions
ferential equations, and it is obviously time-expensive. This study has succeeded by including a first-
Nevertheless, for a process such as SSP, with a high principles large-scale model into an effective MPC
residence time compatible with a control action of 5 or
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 in polymer quality control. The servo and regulatory
of a multivariable control scheme for the SSP process performance of the controller was tested on a detailed
is another original part of this study. The choice of gas SSP reactor simulator, demonstrating that it is effective
temperature and gas water concentration as manipu- and of practical value. The CPU time analysis suggests
lated variables has revealed interesting aspects, high- that a real-time application is feasible, although some
lighting the great effect of the former on process optimization algorithm constraint are needed to improve
performance and the possible importance of the latter the controller robustness.
 Ind. Eng. Chem. Res., Vol. 43, No. 15, 2004 4277

The aim of this work has been to demonstrate that a r ) radial distance from the center of a spherical polymer
first-principles model, even if appropriately sized, can particle, cm
be of real value for control purposes. In other words, a R ) ideal gas constant, L‚atm/(mol‚K)
reduced nonlinear model, integrated into the control Re ) Reynolds number
algorithm, is the topic subject of our work, and hope- Rj ) rate of reaction j, mol/(L‚min)
fully, in the future, it will be compared with a more t ) Time, min
traditional approach (step response) based on industrial Ts ) temperature of the polymer, K
field data. Tg ) temperature of the gas phase, K
TR ) superficial temperature of the polymer particle, K
TW ) temperature of the reactor wall, K
Acknowledgment vg, vp ) gas and polymer flow rates, respectively, cm/min
The authors are grateful to Professor Kevin Yao for yj ) mole fraction of component j in the gas phase
his courtesy. z ) axial distance, measured from the top of the reactor,
cm
Nomenclature Greek Letters
ac ) cross-sectional area of reactor, cm2 χc ) cystalline fraction
C̃j ) average concentration of species within the solid χmax ) maximum crystalline fraction
particles, mol/L  ) voidage
Cj,g ) concentration of volatile species in the gas phase, Φj ) volumetric fraction of component j in the polymer
mol/L κg, κp ) thermal conductivities of the gas phase and
κp, Cp,g ) heat capacities of the polymer and purge gas, polymer, respectively, kcal/(min‚cm‚K)
respectively, kcal/(kg‚K) ηg ) viscosity of the gas phase, kg/(min‚cm)
Cpv,j ) heat capacity of volatile vapors, kcal/(kg‚K) Fg, Fp ) densities of the gas phase and polymer, respec-
Db ) dispersion coefficient, cm2/min tively, kg/cm3
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/ Literature Cited
min (1) Morari, M.; Lee, J. H.; Model Predictive Control: Past,
H ) Henry’s law coefficient, atm Present and Future. Commun. Chem. Eng. 1999, 23, 667.
H ) reactor height, cm (2) Krishnan, A.; Kosanovich, K. A.; DeWitt, M. R.; Creech, M.
hg ) gas-side heat-transfer coefficient, kcal/(cm2‚K‚min) B. Proceedings of the American Control Conference; IEEE: Phila-
hp ) polymer-side heat-transfer coefficient, kcal/(cm2‚K‚ delphia, PA, June 24-26, 1998; pp 3386-3390.
min) (3) Algeri, C.; Rovaglio, M. Dynamic Modeling of a Poly-
ho ) overall heat-transfer coefficient, kcal/(cm2‚K‚min) (ethylene terephthalate) Solid-State Polymerization Reactor I:
hw ) gas-reactor wall heat-transfer coefficient, kcal/(cm2‚ Detailed Model Development. Ind. Eng. Chem. Res. 2004, 43,
K‚min) 4253.
(4) Blom, J. G.; Trompert, R. A.; Verwer, J. G. Algorithm 758:
∆HCR ) heat of crystallization, kcal/kg VLUGR2: A Vectorizable Adaptive-Grid Solver for PDEs in 2D.
∆Hev,j ) heat of vaporization of component j, kcal/mol ACM Trans. Math. Software 1996, 22, 302.
jD ) Colburn factor for mass transfer (5) Yao, K. Z.; McAuley, K. B.; Marchildon, E. K. Simulation
jH ) Colburn factor for heat transfer of continuous solid-phase polymerization of nylon 6,6 (III). J. Appl.
kj ) kinetic rate constant for reaction j, (L/mol)/min Polym. Sci. 2003, 89, 3701.
Kj ) equilibrium constant for reaction j (6) Carslaw, H. S.; Jaeger, J. C. Conduction of Heat in Solid,
kc ) kinetic constant of crystallization, 1/min 2nd ed.; Oxford University Press: Oxford, U.K., 1959.
kg,j ) mass-transfer coefficient of component j, mol‚cm/(L‚ (7) Mallon, F. K.; Ray, W. H. Modeling of solid-state polycon-
atm‚min) densation (II), Reactor design issues. J. Appl. Polym. Sci. 1998,
69, 1775.
ko,j ) overall mass-transfer coefficient of component j, mol‚ (8) Leffew, K. W.; Yerrapragada, S. S.; Deshpande, P. B. Six
cm/(L‚atm‚min) Sigma and Solid-State Polymerization. Chem. Eng. Commun.
kp,j ) Polymer-side mass-transfer coefficients, mol‚cm/(L‚ 2001, 188, 109.
atm‚min) (9) Biegler, L. T.; Cervantes, A. M.; Wachter, A. Advances in
m̆g, m̆p ) mass flow rates of gas and polymer, respectively, Simultaneous Strategies for Dynamic Process Optimization. Chem.
kg/s Eng. Sci., 2002, 57, 575.
Mj ) molecular weight of component j, kg/mol
Received for review December 18, 2003
P ) reactor pressure, atm
Revised manuscript received May 5, 2004
pj,g ) partial pressure of component j in the gas phase, atm Accepted May 14, 2004
p̃j ) vapor pressure of component j in equilibrium with the
polymer average concentration, atm IE034324I