Dynamic Optimization of HIPS Open-Loop Unstable Polymerization

Ind. Eng. Chem. Res. 2005, 44, 2659-2674

2659

Dynamic Optimization of HIPS Open-Loop Unstable Polymerization Reactors Antonio Flores-Tlacuahuac* and Lorenz T. Biegler Department of Chemical Engineering, Carnegie Mellon University, 5000 Forbes Avenue, Pittsburgh, Pennsylvania 15213

Enrique Saldı´var-Guerra† Centro de Investigacio´ n y Desarrollo Tecnolo´ gico, Avenida de los Sauces 87-6, Lerma, Edo de Me´ xico 52000, Me´ xico

We consider dynamic optimization strategies for grade transitions for high-impact polystyrene reactors. Because our desired operating conditions are at unstable points, we apply a simultaneous dynamic optimization (SDO) approach, where state and control variables in the optimal control problem are discretized and a large-scale nonlinear programming solver is applied. For this purpose, we consider Radau collocation on finite elements and the IPOPT NLP solver. In addition, we describe the stability of the SDO strategy through the presentation of dichotomy properties for boundary-value problems. The resulting SDO approach is then demonstrated on a wide variety of operating scenarios for the high-impact polystyrene (HIPS) reactor, with highly reliable and efficient performance results. 1. Introduction The polymer industry represents an important segment of the chemical processing industry with around 100 million tons/year produced worldwide.1 Moreover, the growth of the European polymer market is reported in the range of 1-8% depending on the type of polymer.2 Polymerization reactors are used to manufacture diverse types of the same kind of polymer; such products are called grades. In a continuous plant, the term grade transition refers to the application of a set of heuristic rules, based mainly on plant experience, for carrying out product plant transitions. The use of heuristic rules leads to nonhomogeneous grade products. Therefore, to minimize the production of off-specification products, there is a clear need for systematic grade transition procedures. Given strong marketing pressures, there is also a need to operate processing systems closer to their design limits. Moreover, although the analysis and control of chemical reactors is not a new research topic,3 it continues to attract the attention of both practical and theoretical researchers.4 To improve the performance of polymerization reactors, it appears that operation around highly nonlinear dynamic regions would be economically desirable. Nonlinearities in polymer reactions stem from complicated interactions among kinetics, heat- and mass-transfer limitations, and model uncertainties. Such nonlinear dynamic behavior is characterized by the presence of input/output multiplicities, stability interchanges, oscillations, and chaotic behavior.1 However, operation around regions displaying nonlinear behavior seems to * To whom correspondence should be addressed. On leave from Universidad Iberoamericana. E-mail: antonio.flores@ uia.mx, [email protected]. Tel./fax: +52(55)59504074. Web: http://200.13.98.241/∼antonio. † Present address: Centro de Investigacio´n en Quı´mica Aplicada Blvd. Enrique Reyna 140, Saltillo, Coahuila CP 25100 Me´xico.

be avoided because such operation is often considered to be difficult to achieve. It seems to be a generally accepted idea that operation around such regions might lead to operability problems and, therefore, that operation should move to a new region where nonlinearities are less severe. Process design studies often avoid open-loop unstable regions5-7 because of the perception that adequate process operation would be hard to achieve. On the other hand, from a control point of view, instability might not be a major issue. In fact, one of the basic aims of control theory is to design control systems that are able to deliver satisfactory control performance for unstable processes. Even open-loop stable processes require control systems because of unavoidable disturbances and modeling errors. Also, operation at unstable points can be desirable and might even be required. Some recent evidence shows that, from an economic point of view, unstable operation might be well justified. For example, Blanco et al.7 recently presented a case study of the Tennessee Eastman process where the imposition of stability constraints leads to optimal operating costs that are over 19 times the operating costs of the openloop unstable process. Another interesting example dealing with unstable operation is provided by Varma8 for a serially connected train of isothermal continuous stirred tank reactors (CSTRs). For seven CSTRs in series, over 500 steady-state solutions are reported, but only eight of them are open-loop stable. It is therefore highly likely that the best steady-state operating point will be unstable. As a result, we claim that operation around unstable regions should not be avoided only on operability grounds; sometimes, open-loop unstable steady states might be one of the best ways to run a process. In particular, this is especially true in polymerization reactors that display the classic three-steady-state pattern (see Figure 1). Here, the intermediate steady state is always unstable, whereas the upper and lower

10.1021/ie049534p CCC: $30.25 © 2005 American Chemical Society Published on Web 12/23/2004

2660

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005

regions of challenging nonlinear behavior. Finally, section 5 summarizes these results and concludes the paper. 2. Dynamic Optimization Formulation

Figure 1. Classical three-states multiplicity pattern.

steady states might be either stable or unstable. Operation at the lower steady state is not profitable because of the low conversion yields. From an economic point of view, the upper steady state would be preferred if not for the presence of the gel effect. Hence, the intermediate steady state, which features open-loop instability, is better suited to carry out the polymerization reactions; yield is not as high, but the gel effect is avoided. Optimal grade transition policies have been addressed in a number of studies.9-12 However, in all of these works, only optimal transition trajectories between open-loop stable steady states were addressed. In this work, we address the problem of computing optimal grade transition policies for continuous systems operating at unstable nominal operating points embedded in a multiplicity region. We have previously reported an initial effort to tackle this problem.13 In this work, we apply a robust dynamic optimization formulation able to compute optimal transient trajectories with modest computational effort. One of our aims is to show that a simultaneous dynamic optimization (SDO) approach is the only optimal control methodology able to solve unstable transition problems without stabilizing the initial and final desired steady states beforehand. High-impact polystyrene (HIPS) polymerization processes have significant commercial importance. Styrenic polymers are among the five families of polymer commodities that have the largest world production, and in that family, general-purpose polystyrene (GPPS) and HIPS are by far those of largest consumption. Many of the commercial processes for the production of HIPS have a CSTR as a first reactor in the process, and very little has been published regarding the optimal operation of either a specific reactor or the train of reactors in the process. Here, we apply the SDO methodology for the calculation of optimal transient trajectories for a CSTR in a HIPS polymerization process operating around an open-loop unstable process region. In the next section, we describe the optimization problem formulation and the simultaneous dynamic optimization strategy to solve it. Here, we also include a discussion of the relation of this problem formulation to boundary-value problems (BVPs) and properties of BVPs that handle open-loop instability. Section 3 then presents a mathematical model for the HIPS polymerization reaction process and provides some detail on its dynamic features. Section 4 presents a number of grade transition case studies for the HIPS reactor, even in

Simultaneous dynamic optimization (SDO) provides a way to compute optimal dynamic policies such as grade transitions, even in the presence of challenging nonlinear behavior. These include transitions to unstable points, optimizations with chaotic systems,14,15 and systems with limit cycles and bifurcations. In SDO, the computation of optimal transitions policies reduces to the solution of a nonlinear optimization problem.16 The solution of such problems provides values of the decision variables (i.e., the manipulated variables) that minimize transition time and off-specification product. To solve dynamic optimization problems, three main approaches have been suggested:16 (a) iterative methods based on variational conditions, (b) feasible path nonlinear programming methods, and (c) simultaneous nonlinear programming (SDO) methods. Iterative methods based on variational conditions use Pontryagin’s maximum principle. Algorithms in this family have also been termed control vector iteration algorithms; they involve the integration of a model forward in time and adjoint equations backward in time and updating of the control profile. This approach is well suited only for small-scale problems without states and final constraint problems.17 In the feasible path approach (also termed the control vector parametrization approach), the control variable is represented as a piecewise polynomial function.18 In this approach, the model is solved in an inner loop, and the parameters representing the control variable are updated outside the inner loop. A major problem in this approach refers to the computation of gradients with respect to the states; this has been addressed by using sensitivity equations for the DAE model.19 Despite the success of this technique in solving large-scale optimal control problems, it still requires the repeated numerical integration of the model, and this is a time-consuming task for large-scale problems. Moreover, it is well-known that sequential approaches have properties of single shooting methods that cannot handle open-loop instability.20,16 In the SDO approach, both the manipulated and controlled variables are discretized along the expected solution trajectory. This leads to a nonlinear programming problem whose solution provides the full manipulated and controlled time variable profiles (early SDO formulations were advanced by Tsang et al.,21 Oh and Luus, 22 and Neuman and Sen23). Because fast variations in the controlled variables might arise, the whole solution space is commonly divided into time intervals called finite elements. Inside each finite element, the differential-algebraic equations are satisfied at Radau collocation points. This approach corresponds to a fully implicit Runge-Kutta method with high-order accuracy and stability properties. Other methods based on different collocation discretizations24 and backward difference formulas25 have also been used. 2.1. Formulation. A common requirement in polymerization reactors is that the grade transition policy features a minimum transition time because this implies minimum waste material and minimum utility consumption; it also means that the reactor can be used earlier for manufacturing the desired products. The

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2661

minimum transition policy requires setting of the following optimization problem

min s.t.

(1)

ncol

dz(t) ) F(z(t),y(t),u(t),t,p) dt

(2)

0 ) G(z(t),y(t),u(t),t,p)

(3)

initial conditions z(0) ) z0 U

z e z(t) e z

yL e y(t) e yU

(5)

uL e u(t) e uU pL e p e pU

where F is the vector of right-hand sides of differential equations in the DAE model; G is the vector of algebraic equations, assumed to be index 1; z is the differential state vector; z0 represents the initial values of z; zˆ is the new desired transition state; y is the algebraic state vector, u is the control profile vector; p is a timeindependent parameter vector; and θ is the transition horizon. The DAE optimization problem is converted into a nonlinear program (NLP) by approximating state and control profiles by a family of polynomials on finite elements (t0 < t1 < . . . < tne ) θ). Here, we use a monomial basis representation for the differential profiles, as follows

∑ Ωq q)1

( )

dz

hi

dti,q

Ω′q(Fr) ) δq,r

for q ) 1, ..., ncol for q, r ) 1, ..., ncol

ncol

∑ q)1

Ωq

( ) t - ti-1

dz

hi

dti,q

hi

(8)

ui,q

(9)

for q, r ) 1, ... , ncol

ne

min Φ )

ncol

ωj||z(ti,j) - zˆ ||2 ∑ ∑ i)1 j)1 hi

(10)

Finally, substitution of eqs 6-9 into eqs 1-5 leads to the following NLP

min f(x)

(11)

x∈Rn

c(x) ) 0

(12)

xL e x e xU

(13)

where

where Fr is the location of the rth collocation point within each element. Continuity of the differential profiles is enforced by

zi ) zi-1 + hi

t - ti-1

yi,q

From eq 6, the differential variables are required to be continuous throughout the time horizon, whereas the control and algebraic variables are allowed to have discontinuities at the boundaries of the elements. It should be mentioned that, with representation 6, the bounds on the differential variables are enforced directly at element boundaries; however, they can be enforced at all collocation points by writing appropriate point constraints. In addition, the integral objective function is approximated with Radau quadrature with ne finite elements and ncol quadrature points in each element. This leads to the following objective function

(6)

where zi-1 is the value of the differential variable at the beginning of element i, hi is the length of element i, dz/ dti,q is the value of its first derivative in element i at the collocation point q, and Ωq is the polynomial of order ncol satisfying

Ωq(0) ) 0

ψq(Fr) ) δq,r

s.t.

t - ti-1

hi

Here, yi,q and ui,q represent the values of the algebraic and control variables, respectively, in element i at collocation point q. ψq is the Lagrange polynomial of order ncol satisfying

(4)

bounds

z(t) ) zi-1 + hi

∑ ψq q)1

u(t) )

semi-explicit DAE model

ncol

ψq

q)1

∫0θ||z(t) - zˆ (t)||2 dt

L

∑

( ) ( ) t - ti-1

ncol

y(t) )

(7)

Here, Radau collocation points are used because they allow constraints to be set easily at the end of each element and to stabilize the system more efficiently if high-index DAEs are present. In addition, the control and algebraic profiles are approximated using a similar monomial basis representation that takes the form

x)

(

dz ,z ,y ,u ,t,p dti,q i i,q i,q

f: Rn f R,

and

)

T

c: Rn f Rm

2.2. Dichotomy. When using Pontryagin’s maximum principle, optimal control problems are naturally cast in terms of two-point boundary-value problems (BVPs).17 Hence, the stability of optimal control problems can be appreciated by examining the stability of BVP solvers. It is important to note that this approach applies to the problem formulation itself and does not depend on the computational algorithm for these optimal control problems. There is a fundamental reason to explain why the SDO approach is able to successfully solve unstable optimal control transition problems, and this is based on the dichotomy property20 usually invoked to guarantee stability properties for boundary-value problems. Dichotomy is the BVP equivalent of the asymptotic stability concept normally used in linear and timeinvariant initial-value problems (IVPs). It requires that all eigenvalues of the linearized Jacobian be located strictly inside the stability region, where, for a continuous-time system, the stability region is defined as the

2662


left half of the complex plane (not including the imaginary axis).26 However, in BVP problems, the Jacobian of a well-conditioned problem can have both negative and positive eigenvalues.27 Consider the following two-point BVP

z˘ ) f(z), 0 < t < 1

(14)

B0z(0) + B1z(1) ) b

(15)

subject to

The BVP has dichotomy if there exists a constant κ of moderate size such that the following inequalities hold

||G0(t,s)|| e κ, s e t

(16)

||G1(t,s)|| e κ, s > t

Here, G0(t,s) and G1(t,s) are Green’s functions given by

G0(t,s) ) Φ(t)B0Φ(0)Φ-1(s)

(17)

G1(t,s) ) -Φ(t)B1Φ(1)Φ-1(s)

(18)

and X(t) and Φ(t) are the fundamental solutions of the IVP and BVP, respectively, given by

Φ(t) ) X(t)Q-1

(19)

Q ) B0 + B1X(1)

(20)

In the following example, we show a simple case of a BVP whose stability properties change depending the way the boundary values are specified. To illustrate, we consider the following linear timeinvariant, homogeneous BVP with λ > 020

[ ] [ ][ ] z˘ 1 z1 λ 0 z˘ 2 ) 0 -λ z2

[ ] [] z1(1) 1 ) z2(0) 1

Here, X(t), B0, and B1 of the above problem are given by

[

eλt 0 X(t) ) 0 e-λt

]

[ ]

[ ]

0 1 B0 ) 0 0

0 0 B1 ) 1 0

and the Green’s functions are given by

G0(t,s) )

[

0 0 0 eλ(s-t)

[

]

λ(t-s) 0 G1(t,s) ) - e 0 0

∀t g s

]

∀t < s

In this case, the BVP is stable, as the exponential terms always decay because κ g e-λ is easily satisfied. On the other hand, if the problem is analyzed as an initial value problem, for the same X(t) with B0 ) I and B1 ) 0, then this problem is totally unstable. The above example clearly shows that asymptotic instability of the IVP does not necessarily imply instability of the BVP, because stability properties of the BVP also depend on the boundary conditions. To enforce stability, one needs to pin down unstable modes with

the appropriate boundary conditions. Hence, the solution of a linear BVP can be split into two spaces: a subspace with increasing modes and a subspace with nonincreasing modes. A well-conditioned BVP (i.e., moderate κ values) and stability follow if and only if dichotomy holds, i.e., increasing modes are bounded by final conditions and nonincreasing modes have corresponding initial conditions. This property has been analyzed in ref 28. Proving whether a problem has or does not have dichotomy properties is much harder for a larger nonlinear system, and in general, dichotomy is shown instead for linearized BVPs about a fixed trajectory. Nevertheless, solution of BVPs with finite difference or collocation methods leads to discretized linearized systems that have conditioning constants closely related to κ in inequalities 16. Therefore, dichotomy guarantees a well-conditioned discretized system as well. To extend this concept to optimization, we now consider the solution of eqs 11-13. Cuthrell and Biegler29 showed that the Karush-Kuhn-Tucker (KKT) conditions of this NLP are equivalent to collocation applied to the boundary-value problem generated by the optimality conditions of eqs 1-5. Moreover, Biegler et al.16 described an efficient matrix decomposition strategy that preserves well-conditioning based on dichotomy. As a result, the dichotomy and conditioning properties of this BVP system are preserved by the SDO method, regardless of the stability characteristics of the DAE model itself. On the other hand, competing feasible path methods rely heavily on IVP solvers that cannot handle increasing modes. This essential difference allows SDO to deal with the challenging nonlinear characteristics discussed in the next sections. 2.3. Formulation and Solution of Nonlinear Program. The dynamic optimization formulation given in eqs 11-13 was implemented using the AMPL mathematical programming language30 and solved using the IPOPT algorithm31 for large-scale nonlinear programming. This algorithm follows a barrier approach, where the bound constraints in inequality 13 are replaced by logarithmic barrier terms that are added to the objective function to give n

n

(i) (i) ln(x(i) - x(i) ∑ L ) + ∑ln(xU - x )] i)1 i)1

min f(x) - µ[

s.t.

c(x) ) 0

(21) (22)

with a barrier parameter µ > 0. Here, x(i) denotes the ith component of the vector x. Because the objective function of this barrier problem becomes arbitrarily large as x(i) approaches either of its bounds, a local solution x/(µ) of this problem lies in the interior of this set, i.e., xU > x/(µ) > xL. The degree of influence of the barrier is determined by the size of µ, and under certain conditions, x/(µ) converges to a local solution x/ of the original problem 11-13 as µ f 0. Consequently, a strategy for solving the original NLP is to solve a sequence of barrier problems 21 and 22 for decreasing barrier parameters µl, where l is the counter for the sequence of subproblems. IPOPT follows a primal-dual approach and applies a Newton method to the resulting KKT conditions. Exact first and second derivatives for this method are provided automatically through the AMPL interface. More information on IPOPT can be found in Waechter and Biegler.32

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2663 Table 1. Analyzed Casesa

3. Mathematical Model In this section, we describe the differential algebraic equation (DAE) model of the free-radical bulk polymerization of the system styrene-polybutadiene, using a monofunctional initiator. The polymerization reactions are carried out in a nonisothermal CSTR, and the model assumes perfect mixing, constant physical properties, quasi-steady state, and the long-chain hypothesis. Constant volume in the reactor has been also assumed. Changes in density of the monomer-polymer mixture have been neglected. The kinetic mechanism involves initiation, propagation, transfer, and termination reactions. Polybutadiene is added to guarantee desired mechanical properties by promoting grafting reactions. Taking into account all of these modeling assumptions, the dynamic mathematical model of the nonisothermal CSTR is given as follows

dCi FiCoi - FCi ) - KdCi dt V

(23)

dCm F(Com - Cm) ) - KpCm(µ0r + µ0b) dt V

(24)

dCb F(Cob - Cb) ) - Cb(Ki2Cr + Kfsµ0r + Kfbµ0b) (25) dt V dCr ) 2f*KdCi - Cr(Ki1Cm + Ki2Cb) dt

(26)

dCbr ) Cb[Ki2Cr + Kfb(µ0r + µ0b)] - Cbr[Ki3Cm + dt Kt(µ0r + µ0b + Cbr)] (27) dµ0r ) 2Ki0Cm3 + Ki1CrCm + CmKfs(µ0r + µ0b) dt [KpCm + Kt(µ0r + µ0b + Cbr) + KfsCm + KfbCb]µ0r + KpCmµ0r (28) dµ0b ) Ki3CbrCm - [KpCm + Kt(µ0r + µ0b + Cbr) + dt KfsCm + KfbCb]µ0b + KpCmµ0b (29) 0 0 o dT F(T - T) ∆HKpCm(µr + µb) UA(T - Tj) ) dt V FsCps FsCpsV (30)

dTj Fj(Toj - Tj) UA(T - Tj) + ) dt Vj FjCpjVj

(31)

The definitions of the model symbols are provided in the Nomenclature section. The details of the kinetic mechanism, rate constants, and design parameters are contained in ref 33; complete nonlinear bifurcation analysis of the aforementioned polymerization model has been presented there as well. 4. Computation of Unstable Transition Trajectories Figure 2 displays the nonlinear bifurcation diagram of the HIPS CSTR using as continuation parameters the

initial point case

Fcw

A1 A2 A3 A4 A5 B1 B2 B3 C1 C2 D1 D2 D3 D4

1 1 1 1 1

1 1 0 0 0 0

Fi

0.0015 0.0015 0.0015 0.0015 0.0015 0 0 0

final point

F

0 0 0

Fcw U U U U U U U U U U SU SU SU SU

Fi

F

0.5 0.1 0.1 0.5 0.5

10 10 1 1 1 1

0.02 0.04 0.05 0.05 0.11 0.0015 0.0015 0.0015

1.14 1.14 1.14

U U S S S U U U, TP U U, TP U U S S

a

U denotes an open-loop unstable steady state, S denotes an open-loop stable steady state, TP means a turning point, and SU stands for start-up conditions. Blanks mean that the value of the given variable was kept at its nominal value. All flow rates are in L/s. The second column shows the initial values of the manipulated variables, and the third column contains the end values, when the transition is completed. Table 2. Statistics of the Grade Transition Problemsa case

n

m

NFE

iters

CPU time

A1 A2 A3 A4 A5 B1 B2 B3 C1 C2 D1 D2 D3 D4

496 496 996 746 1496 496 746 496 556 556 836 926 1116 1856

316 316 636 476 956 316 476 316 316 316 476 476 636 956

20 20 40 30 40 20 30 20 20 20 30 30 40 60

35 223 67 73 68 330 450 289 869 292 119 36 952 893

1.69 5.01 6.06 4.2 3.24 6.3 13.57 6.91 18.66 9.07 7.47 0.95 54.82 72.88

a n is the number of variables, m is the number of constraints, NFE is the number of finite of elements (two internal collocation points were always used), iters is the number of iterations. CPU time is in seconds using a Pentium 4 PC running at 1.6 GHz.

cooling water and initiator flow rates; the only type of bifurcations shown are of static nature (saddle points). As usual, open-loop stable steady states are represented by the continuous line, and unstable steady states are represented by the dashed line. Under nominal operating conditions, the reactor exhibits three steady-state values (denoted by N1, N2, and N3). The lower stable steady-state branch represents extinction conditions because temperature is so low there that polymerization reactions cannot proceed to any meaningful rate. Along the intermediate unstable steady-state branch, monomer conversion rises to meaningful values (around 30%). This is the zone in which our HIPS CSTR is operated; the selected nominal point is denoted by N2 and corresponds to an open-loop unstable operating point. The upper stable steady-state branch denotes a highmonomer-conversion operation region. To compute the optimal transition unstable trajectories, several cases were formulated (all of the analyzed cases are summarized in Table 1. The grade transition problems were solved on an Intel Pentium 4, 1.6-GHz PC running the Linux operating system. Problem sizes and performance of the optimizer are shown in Table 2). Except in four cases, all of the transition trajectories always started from the same steady-state conditions

2664


Figure 2. HIPS reactor multiple-steady-state diagrams using the (a) cooling water and (b) initiator flow rates as the main bifurcation parameters.

(the N2 unstable steady state in Figure 2). The A, B, and C transitions always end up around a region of lower monomer conversion where the polymerization reactions tend to extinguish, whereas D transitions were formulated for seeking optimal transitions for start-up conditions. D transitions end up around either the intermediate or high-monomer-conversion region but use different sets of manipulated variables. A1 Transition. For this transition, the cooling water flow rate was taken as the manipulated variable (see Figure 2a). From Figure 3, we notice that the dynamical optimal transition trajectory implies complete closing of the cooling water flow rate for most of the transition time; the control valve remains saturated, hitting its lower bound except at the end of the transition when it goes to its final desired value. The transition time is short (around 20 min) compared to the reactor time constant. The behaviors of the rest of the state trajectories shown in Figure 3 are what we normally expect to see in a polymerization reaction system. Because of a decrease in reactor temperature, the propagation rate step also reduces, leading to an increase in the monomer

concentration. Together with this falling in the polymerization rate, less initiator is used, and the polymer concentration starts decreasing. However, toward the end of the transition, conversion increases, ending up at a value higher than the initial one. A2 Transition. For this transition, the cooling water flow rate was also taken as the manipulated variable; see Figure 2a. As can be seen in this figure, the A2 case might mean a difficult-to-achieve transition because it is closer to the turning point than the A1 transition was. Case A2 was set up for comparison purposes against the A1 case; we wanted to have a clear idea about potential computational difficulties when operating near turning points. As expected, the A2 transition trajectory did not represent a problem from a convergence point of view. In Figure 4, we can see that all of the states, except the cooling jacket temperature, have reached the desired end values in about 1 h. However, because the jacket temperature still has a small offset, the cooling water flow rate recovers from saturation and finally goes to its final value; the complete optimal transition trajectory takes around 1.3 h. The dynamic behaviors


Figure 3. Case A1.

Figure 4. Case A2.

of the rest of the states follow a similar pattern previously discussed under the A1 trajectory. A3 Transition. The aim of performing this transition is two-fold. On one side, we wanted to verify that our dynamic optimization formulation is also able to track trajectories ending in an open-loop stable steady state. On the other side, we wanted to carry out a comparison of the shapes of the profiles of the manipulated and controlled time-variant variables when a transition is sought to two different steady-state solutions. The A3 and A2 steady-state operating points (see Figure 2a) are two of the three steady-state solutions present at a 0.1

L/s cooling water flow rate; A3 represents an open-loop stable steady state. Figure 5 displays the shapes of the transition profiles. The shapes of the profiles of the A3 and A2 transitions are qualitatively similar, except that, now, the cooling water flow rate takes initially a large value. This is easy to understand because the A3 steady state operates at a lower reactor temperature, so the cooling system has a greater demand. This is also the reason the transition trajectory time increases to 3 h. A4 Transition. Similarly, the A4 transition was setup as a way to compare it against the A1 transition (see Figure 6). A1 and A4 represent two of the three

2666


Figure 5. Case A3.

Figure 6. Case A4.

steady-state solutions that the reactor shows at a 0.5 L/s cooling water flow rate. The qualitative transition profiles in these two transitions are similar except for the jacket reactor temperature and cooling water flow rate profiles. Once again, we found that a larger initial variation in cooling water flow rate was needed because of the large reactor temperature decrease imposed on the A4 transition; this also explains the increase in transition time. A5 Transition. Until now, all of the addressed transitions have featured low monomer conversions. We followed this approach because, in HIPS manufacturing,

a train of series-connected reactors is normally used to achieve higher conversions, so that only modest monomer conversions are obtained in the first reactor. Despite this observation, we carried out a transition to a high-conversion steady state. In fact, as displayed in Figure 2a, the final desired steady state is the highermonomer-conversion steady state when Fw ) 0.5 L/s. As can be noted from Figure 7, because of the monomer conversion’s high sensitivity to increasing temperature, the monomer concentration drops quickly. On the other hand, the reactor temperature features a slow dynamic response. The amount of heat released by the polym-


Figure 7. Case A5.

Figure 8. Case B1.

erization reaction system is enough to reach the final desired reactor temperature; excess heating is withdrawn by the cooling system. Notice that, because the cooling water temperature must lie below the 100 °C physical upper bound, the dynamic optimization formulation finds that, to meet this constraint, the cooling water flow rate should be increased to 0.866 L/s. The results show that, as expected, our dynamic optimization formulation is able to track transitions ending around high-monomer-conversion regions. B1 Transition. For this transition the initiator feed flow rate was taken as the manipulated variable (see Figure 2b). The idea behind using different manipulated variables was to confirm whether potential difficulties

associated with the dynamic optimization transition computations between open-loop unstable steady states depended on the manipulated input. We actually do not claim that the initiator flow rate should be considered the best manipulated variable for carrying out the transition. In Figure 8, we see that the best input policy is just a simple step change from the initial to the final value. The total global transition time is around 2 h. Because the transition implies a decrease in reactor temperature, the polymerization reaction propagation step stops almost completely until both the reactor and jacket temperatures reach their final values. When this happens, around 1.5 h, the propagation step restarts; this is what allows both the monomer and initiator

2668


Figure 9. Case B2.

Figure 10. Case B3.

concentration to reach their final values. However, because the value of the lowest reactor temperature remains around 70 °C, a more likely reason for the stopping of the polymerization reaction has to do with the small and sluggish increase in the initiator concentration. B2 Transition. In this transition, the initiator feed flow rate was taken again as the manipulated variable (see Figure 2b). Case B2 represents a slight variation of the case B1 conditions. This time, the total transition time amounts to 2.7 h (see Figure 9), but the same qualitative dynamic behavior patterns remain. B3 Transition. Again, we kept the initiator flow rate as the input manipulated variable (see Figure 2b). The

main objective of this case was again to verify whether our dynamic optimization approach can handle end transitions located at turning points. Once again, using the simultaneous dynamic optimization approach, no problems were observed when the computations were carried out. This time, the total transition time is around 3 h (see Figure 10). The qualitative dynamic response of case B3 is essentially the same as that observed for the B1 transition. C1 Transition. In this case, the manipulated variables are the cooling water and initiator flow rates. As depicted in Figure 2b, this transition implies large variations in both the cooling water and initiator flow rates; the total transition time is about 1.7 h, as


Figure 11. Case C1.

Figure 12. Case C2.

displayed in Figure 11. The best transition policy found for both input variables is a step change from their initial values to their final values. An important difference with respect to the rest of the analyzed cases is the fast jacket temperature dynamics, a predictable result because of the large increase in cooling water flow rate. However, the temperature reactor dynamics is much slower, as it depends on the polymerization reaction as well as heat-transfer characteristics. The slow observed reactor temperature response is caused by a low polymerization rate that is just large enough to keep the reactor temperature moving slowly. Once again, the best policy implies the extinction of the polymerization reaction during most of the transient

time. Only near the end does the reaction ignite, allowing the monomer concentration to reach its final transition value. C2 Transition. In this case, the manipulated variables are also the cooling water and initiator flow rates. According to Figure 2b, this transition implies changing the initiator flow rate 77 times from its nominal value, a large increase in terms of its relative value. To complicate matters further, the C2 point turns out to be a turning point; however, no computational problems were observed for the computation of the dynamic optimal transition trajectory using the simultaneous dynamic optimization approach. As noted from Figure 12, the total transition time is about 2.5 h. The same

2670


Figure 13. Case D1.

Figure 14. Case D2.

qualitative dynamic response pattern as discussed before for the C1 transition is observed.

D1 Transition. To test our dynamic optimization formulation, we decided to compute optimal transition


Figure 15. Case D3.

policies for the start-up of the reactor so it finally reaches the N2 open-loop unstable nominal operating point. To do so, the initiator and cooling water flow rates were chosen as the manipulated variables. At the start, the reactor was assumed to be operated with a continuous monomer feed stream whose initial concentration and temperature were taken as those of the feed stream. The initial reactor jacket temperature was taken as the cooling water feed stream temperature. Starting from these values, our dynamic optimization formulation was able to find an optimal control policy featuring a minimum start-up time; the time-variant profiles of both states and manipulated variables are shown in Figure 13. The dynamic start-up behavior seems interesting. The optimal policy dictates a step change in the initiator flow rate. However, the cooling water flow rate is not used as a manipulated variable until essentially both the reactor and jacket temperatures have reached their final design values, and then the cooling water flow rate has a ramp-like variation. At the beginning, heat for the observed temperature rise is provided by the highly exothermic polymerization reaction. The transition time is short, around 1 h, taking into account the reactor residence time. An interesting point to note is that, during the 0.1-0.3-h time interval, the monomer concentration vanishes completely. D2 Transition. In this case, we again decided to compute an optimal starting policy but using, in addition to the initiator and cooling water flow rates, the monomer flow rate to check whether a better control policy might be obtained. Initially, the reactor was assumed to contain a monomer mass holdup equivalent to one reactor volume. Therefore, at the start, the monomer concentration and reactor temperature were

assumed to be those of the feed stream; the initial jacket temperature was set to the cooling water feed stream temperature. In Figure 14, the time-variant profiles of the optimal states and manipulated variables are displayed. Compared to the D1 case, no significant improvement was observed in terms of either the transition time or the reactant consumption; in fact, the D2 transition policy seems not to be as good as that found for the D1 case, because more monomer seems to be required to carry out the same transition and a longer transition time, around 2 h, is required. However, it is interesting to remark that, this time, the monomer concentration does not vanish because of the high initial monomer injection; in fact, the dynamic response of the monomer is the fastest one. Also, we note that computation of the optimal transition policy for this case demanded the use of an upper bound on the monomer flow rate (2 L/s); otherwise, an optimal solution of the dynamic optimization problem could not be found. D3 Transition. In all of the transition cases analyzed to this point (except for the A5 transition), we have always imposed changes aiming to reduce the reactor temperature, thus leading to a reduction in monomer conversion and an increase in the molecular weight. To keep a realistic optimization scenario, we did not try to raise monomer conversion beyond its nominal design value since, from an industrial point of view, the first reactor in a HIPS manufacturing plant reaches only about 30% monomer conversion. However, optimal start-up policies were sought assuming that the desired end point was given by the upper stable steady state (point N3 in Figure 2). The aim of running this case was to verify the robustness of our dynamic optimization approach for grade transitions between arbitrary oper-

2672


Figure 16. Case D4.

ating points. In the present case, only the monomer feed stream and cooling water flow rates were considered as manipulated variables. Initially, a mass reactor holdup consisting of both initiator and monomer concentrations was assumed. The initial reactor and cooling water temperatures were taken as the monomer and water feed stream temperatures, respectively. As seen in Figure 15, the computation of the optimal transition policies did not present any additional complexity, although the monomer flow rate remained saturated during a significant duration of the transition period. D4 Transition. To investigate whether a better optimal transition policy might be found, we decided to free the initiator feed stream flow rate and to consider it as an additional manipulated variable. As seen in Figure 16, no significant improvement over case D3 was observed. The reason using the initiator flow rate as an extra manipulated variable did not improve the transition policy becomes clear because, right from the beginning, the optimizer decided to keep the initiator flow rate constant at 1.5 × 10-3 L/s, rendering cases D3 and D4 almost equivalent. 5. Conclusions Frequently, when chemical processes are forced to operate close to optimal steady-state conditions, the desired operating point can end up being open-loop unstable. Also, dynamic optimal set-point changes are sometimes required between unstable operating regions;

this situation emerges particularly during grade transition in polymerization reactors. However, it is clear that the same situation can occur in many other science and engineering fields. Until now, the dynamic optimization of systems exhibiting instabilities has not drawn enough attention in the process systems engineering (PSE) research community. This is partially due to the fact that the use of dynamic optimization approaches based on the open-loop integration of model equations gives rise to unbounded states leading to model numerical integration difficulties. Although this problem can partially be fixed by using a stabilizing controller, the dynamic optimal solution trajectory would depend on the control law structure and its tuning parameters. Moreover, incorporating a controller can induce undesired closed-loop behavior. On the other hand, the SDO approach efficiently handles unstable transition problems because of the fact that the dynamic optimization problem is cast as an NLP without stability worries. Given that the main aim of this work was to propose a novel dynamic optimization approach featuring openloop unstable systems, we did not intend to compute the best transition trajectory for a given polymer grade. In fact, because of the presence of nonconvexities in the formulation of the dynamic optimization problem, it is likely that better transition trajectories can be found. On the other hand, it is interesting to note that the computational demand tends to increase as the optimizer is forced to track trajectories crossing through the


complete nonlinear high-sensitivity regions (see Figure 2a,b and Table 2). This is specially clear for both the D3 and D4 transitions, as they demanded an increase of approximately 2 orders of magnitude in CPU time compared to the rest of transitions. The greater computational load is due to the nonlinearity embedded in the model. Because of the number of states of the HIPS model and the numbers of both internal collocation points and finite elements, the size of the underlying dynamic optimization formulation did not represent any problem for finding the optimal solutions. In this work, we have shown the theoretical foundations that allow the SDO approach to solve unstable transition problems. The solution of the resulting NLP is addressed using state-of-the-art numerical optimization algorithms. Such algorithms easily allow an extension to larger-scale dynamic systems. Finally, our unstable transition dynamic optimization strategy was applied to a realistic industrial HIPS model that has been previously validated against plant data. As future work, we will apply the same SDO approach to the complete HIPS plant composed of seven series-connected CSTRs operating around highly nonlinear and open-loop unstable regions. Acknowledgment Antonio Flores-Tlacuahuac acknowledges the financial support from the Fullbright Commission, Universidad Iberoamericana-Santa Fe, and the Department of Chemical Engineering at Carnegie-Mellon University. We also thank three anonymous reviewers for helping us to improve our presentation. Nomenclature A ) heat-transfer area (m2) Cb ) butadiene concentration (mol/L) Cob ) feed stream butadiene concentration (mol/L) Cbr ) branched radical concentration (mol/L) Ci ) initiator concentration (mol/L) Coi ) feed stream initiator concentration (mol/L) Cm ) monomer concentration (mol/L) Com ) feed stream monomer concentration (mol/L) Cps ) monomer heat capacity [J/(kg K)] Cpj ) cooling water heat capacity [J/(kg K)] Cr ) radical concentration (mol/L) F ) feed stream volumetric flow rate (L/s) Fj ) cooling water volumetric flow rate (L/s) f* ) initiator efficiency, dimensionless kd ) chemical initiation kinetic constant (L/s) kfb ) polybutadiene transfer kinetic constant [L/(mol s)] kfs ) styrene monomer transfer kinetic constant [L/(mol s)] ki0 ) thermal initiation kinetic constant [L2/(mol2-s)] ki1 ) chemical initiation kinetic constant (size 1) [L/(mol s)] ki2 ) polybutadiene activation chemical initiation kinetic constant [L/(mol s)] ki3 ) grafted chemical initiation kinetic constant (size 1) [L/(mol s)] kp ) propagation kinetic constant [L/(mol s)] kt ) coupling termination kinetic constant [L/(mol s)] T ) reactor temperature (K) To ) feed stream reactor temperature (K) Tj ) jacket temperature (K) Toj ) feed stream cooling water (K) U ) heat-transfer coefficient [J/(s K m2)]

V ) reactor volume (L) Vj ) jacket volume (L) Greek Symbols ∆H ) heat of reaction (J/mol) µ0b ) zeroth moment of dead butadiene µ0r ) zeroth moment of dead radicals Fj ) water density (kg/L) Fs ) monomer density (kg/L)

Literature Cited (1) Ray, W. H.; Villa, C. Nonlinear dynamics found in polymerization processessA review. Chem. Eng. Sci. 2000, 55, 275. (2) Moen, O. Machine means for design, control and economic optimization of polyolefine production processes and their products. In Escape-14; Barbosa, A., Matos, H., Eds.; Elsevier: New York, 2004; pp 117-122. (3) Aris, R. Chemical reactors and some bifurcarion phenomena. Ann. N. Y. Acad. Sci. 1979, 7, 314. (4) Antonelli, R.; Astolfi, A. Continuous Stirred Tank Reactors: Easy to Stabilise? Automatica 2003, 39, 1817. (5) Kokosis, A.; Floudas, C. A. Stability in Optimal Design: Synthesis of Complex Reactor Networks. AIChE J. 1994, 40, 849861. (6) Seider, W.; Brengel, D.; Provost, A.; Widagdo. S. Nonlinear Analysis in Process Design. Why Overdesign To Avoid Complex Nonlinearities? Ind. Eng. Chem. Res. 1990, 29, 805. (7) Blanco, A. M.; Bandoni, J. A.; Biegler, L. Re-design of The Tennessee Eastman Challenge Process: An Eigenvalue Optimization Approach. In Proceedings of Foundations of Computer Aided Process Design, manuscript accepted. (8) Varma, A. On the Number and Stability of Steady States of a Sequence of Continuous Stirred Tank Reactors. Ind. Eng. Fundam. 1980, 19, 316. (9) Cervantes, A., Tonelli, S.; Brandolin, A.; Bandoni, A.; Biegler, L. T. Large-Scale Dynamic Optimization of a Low-Density Polyethyle Plant. Comput. Chem. Eng. 2000, 24, 983-989. (10) McAuley, K. B.; McGregor, J. F. Optimal Grade Transitions in a Gas-Phase Polyethylene Reactor. AIChE J. 1992, 38, 10, 1564-1576. (11) Takeda, M.; Ray, W. H. Optimal-Grade Transition Strategies for Multistage Polyolefin Reactors. AIChE J. 1999, 45, 8, 1776-1792. (12) BenAmor, S.; Doyle, F. J.; McFarlane, R. Polymer Grade Transition Control Using Adavanced Real-Time Optimization Software. J. Process Control 2004, 14, 349-364. (13) Silva, A.; Flores-Tlacuahuac, A.; Arrieta J. J. Dynamic Trajectory Optimization Between Unstable Steady States of a Class of CSTRs, In Escape-12; Elsevier: Amsterdam, 2002; pp 547-552. (14) Schloeder-H, J. P.; Bock Kallratch, J. G. Least squares parameter estimation in chaotic differential equations. Celestial Mech. Dynam. Astron. 1993, 56, 353-371. (15) Schloeder V. J. P.; Schultz Bock, H. G. Numerik grosser differentiell-algebraischer gleichungen. simulierung und optimierung. In Prozess Simulation; Schuler, H., Ed.; VCH Verlag: Berlin, 1994. (16) Biegler, L. T. Optimization strategies for complex process models. Adv. Chem. Eng. 1992, 18, 197-256. (17) Bryson, A. E.; Ho, Y. C. Applied Optimal Control; Hemisphere: New York, 1975. (18) Ray, H. W. Advanced Process Control; McGraw-Hill: New York, 1981. (19) Vassiliadis, V. Computational Solution of Dynamic Optimization Problems with General Differential Algebraic Constraints. Ph.D. Thesis, University of London, London, U.K., 1993. (20) Ascher, U. M.; Petzold, L. R. Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations; SIAM: Philadelphia, PA, 1998. (21) Tsang, T. H.; Himmelblau, D. M.; Edgar, T. F. Optimal control via collocation and nonlinear programming. Int. J. Control 1975, 24, 763-768. (22) Oh, S. H.; Luus, R. Use of Orthogonal Collocation Method in Optimal Control Problems. Int. J. Control 1977, 26, 657-673. (23) Neuman, A.; Sen, C. P. A suboptimal control algorithm for constrained problems using cubic splines. Automatica 1973, 9, 601-613.

2674


(24) Betts, J. Practical Methods for Optimal Control Using Nonlinear Programming; SIAM Series on Advances in Design and Control. SIAM: Philadelphia, PA, 2001. (25) Jockenho¨vel, T.; Biegler, L. T.; Wa¨chter, A. Dynamic optimization of the tennessee eastman process using the optcontrolcentre. Comput. Chem. Eng. 2003, 27, 11, 1513-1531. (26) Glad, T.; Ljung. L. Control Theory: Multivariable and Nonlinear Methods; Taylor & Francis: London, 2000. (27) Ascher, U.; Bader, G. Stability of Collocation at Gaussian Points. SIAM J. Numer. Anal. 1986, 23, 412-422. (28) De Hoog, F. R.; Mattheij, R. M. M. On Dichotomy and Well Conditioning in BVP. SIAM J. Numer. Anal. 1987, 24 (1), 89105. (29) Cuthrell, J. E.; Biegler, L. T. Simultaneous optimization and solution methods for batch reactor control profiles. Comput. Chem. Eng. 1989, 13 (1-2), 49. (30) Fourer, R.; Gay, D. M.; Kernighan, B. W. AMPL: A Modeling Language for Mathematical Programming; Thomson Publishing Company: Danvers, MA, 1993.

(31) Wa¨chter, A.; Biegler, L. T. Global and Local Convergence of a Reduced Space Quasi-Newton Barrier Algorithm for LargeScale Nonlinear Programming; CAPD Technical Report B-00-06; Carnegie Mellon University: Pittsburgh, PA, 2000. (32) Waechter, A.; Biegler, L. T. On the implementation of an interior point filter line search algorithm for large-scale nonlinear programming, manuscript accepted. (33) Varazaluce-Garcı´a, J. C.; Flores-Tlacuahuac, A.; Saldı´varGuerra, E. Steady-State Nonlinear Bifurcation Analysis of a HighImpact Polystyrene Continuous Stirred Tank Reactor. Ind. Eng. Chem. Res. 2000, 39, 1972.

Received for review May 28, 2004 Revised manuscript received September 10, 2004 Accepted September 30, 2004 IE049534P

Dynamic Optimization of HIPS Open-Loop Unstable Polymerization

Recommend Documents