Optimal Start-Up and Product Transition Policies of a Reactive

2092

Ind. Eng. Chem. Res. 2007, 46, 2092-2111

Optimal Start-Up and Product Transition Policies of a Reactive Distillation Column Rodrigo Lo´ pez-Negrete de la Fuente and Antonio Flores Tlacuahuac* Departamento de Ingeniera y Ciencias Quı´micas, UniVersidad Iberoamericana, Prolongacio´ n Paseo de la Reforma 880, Me´ xico DF, 01210, Me´ xico

In the processing industry, there is the need to approach the operation of industrial equipment so they increase their energy efficiency, leading to more-economical and environmentally oriented processes. A feasible way to achieve these purposes lies in the optimal dynamic operation of industrial operations. In this work, we address the dynamic optimal operation of a highly energy-intensive industrial operation (distillation) under two common operating scenarios: start-up and steady-state transition operations. The start-up operations and steady-state transitions of a reactive distillation column were addressed in this paper by approaching the problem as a dynamic optimization problem. A detailed tray-by-tray model that considers internal tray hydraulics, but ignores vapor dynamics, was derived and used for calculations. Several manipulated variables were considered, besides the reboiler heat duty. The large scale nonlinear programming (NLP) problems generated from the application of the simultaneous dynamic optimization method were successfully solved with the use of an interior point optimizer (IPOPT). It was determined that, with the use of dynamic optimization, large time reductions can be obtained, when compared to empirical ramplike changes in the manipulated variables, thus reducing the amount of waste and energy consumption. Overall, when using optimal, rather than empirical dynamic operation policies, energy and raw material savings on the order of one order of magnitude, which clearly demonstrates that significant economic and environmental advantages can be achieved by approaching the dynamic operation of industrial processes as a formal dynamic optimization problem. 1. Introduction Reactive distillation (RD) is the process in which vaporliquid separation and one or more chemical reactions occur simultaneously. In this way, only one piece of equipment (the RD column) is used, possibly reducing investment and operation costs.1,2 There are several other advantages associated with this process, such as (i) chemical equilibrium limitations can be eliminated, (ii) parallel reactions can be reduced, and/or (iii) azeotropes can be reacted away.3 Selectivity can be increased by maintaining low concentrations of one of the reactants, thus reducing the reaction rate of parallel reactions. If the reaction that is occurring is exothermic, the heat generated could be used to reduce the reboiler heat requirement. Finally, in packed columns, the amount of catalyst required may be reduced. Unfortunately, RD processes have some disadvantages that are related to the nonlinearities introduced by chemical reaction kinetics, coupled with the thermodynamic equilibrium relationships. Such nonlinearities could result in the emergence of input/ output multiplicities,1,4 and because of this, RD columns have a tendency to be difficult to control. To operate distillation columns, the design of the column will likely be based on a mathematical model of the process. Because there are strong nonlinearities (trough mass, energy and chemical kinetic couplings), the modeling and design problem could be difficult to handle, which means that simulation and design tasks will be rather complicated. Of course, difficulties associated with the numerical solution of simulation and design tasks do not necessarily represent a disadvantage to the use of reactive columns, but operability problems clearly do. Initial research in RD focused on single steady-state calculations,5 simulation,6,7 * To whom correspondence should be addressed. Tel./Fax: (+52) 55 59504074. E-mail: [email protected]. URL: http://200.13.98.241/ ∼antonio.

and dynamic modeling.3,8 However, recent studies have been performed to address the calculation of multiple steady states, process synthesis, and open- or closed-loop operation. Several papers have been written about the emergence of multiplicity of steady states in RD. In the production of ethylene glycol, Ciric and Miao4 found multiple steady states by homotopy continuation techniques. They used the volume holdup as their continuation parameter, and found the existence of output multiplicities. Using the same system but with the boil-up ratio, which is a variable likely to be used as a manipulated variable in a closed-loop control system, as the bifurcation parameter, Monroy-Loperena and Alvarez-Ramı´rez9 also observed both input and output multiplicities. Other processes, such as the production of ethyl tert-butyl ether10 and methyl tert-butyl ether,11 have been shown to present output and input multiplicities. In the production of ethyl acetate, Vora and Daoutidis1 determined that only input multiplicities existed. This type of nonlinear behavior presents challenging problems for process synthesis, operation, and control. For process synthesis purposes, some techniques use graphical methods (for example, residue curve maps) to analyze ideal and nonideal mixtures.12 Other authors proposed a method using RD line diagrams based on component concentration profiles,13 or the use of a method based on the McCabe-Thiele method used in nonreactive distillation.14 On the other hand, the use of optimization and iterative methods applied to RD synthesis has gained popularity. For instance, mixed-integer programming has been used for equilibrium-limited reactions.15 Also, for nonequilibrium RD, there has been some work using tray-by-tray, rigorous models,2 and, more recently, the use of disjunctive programming has been applied to solve the synthesis problem in a more systematic way.16 Closed-loop operation of RD columns has also been studied by different research groups, mainly focusing on controller

10.1021/ie061312z CCC: $37.00 © 2007 American Chemical Society Published on Web 03/03/2007

Ind. Eng. Chem. Res., Vol. 46, No. 7, 2007 2093

design and configuration. Working with the production of ethylene glycol, Kumar and Daoutidis17 addressed the dynamic modeling and control of such RD columns. They compared the use of a nonlinear inversion-based controller with a simple PI controller, and found that nonlinear control worked best for setpoint changes. On the other hand, Al-Arfaj and Luyben,18 who were working with the same system and using Aspen Dynamics, concluded that only simple PI controllers were needed. This could be attributed to the fact that both groups used different modeling techniques, software, and assumptions. In the production of ethyl acetate, Vora and Daoutidis,1 using nonlinear feedback controllers, effectively addressed the problem of setpoint changes and disturbance rejection. Furthermore, they also compared their results with a simple PI controller, concluding that nonlinear control worked better. There are several other papers in the open literature that have addressed the design and configuration of different types of controllers for RD that discuss set-point changes and/or disturbance rejection. (Please refer to Al-Arfaj and Luyben19 and Luyben et al.20 for more-detailed reviews.) The startup of RD columns has also been addressed, but to a lesser degree. Alejski and Duprat21 described a dynamic model for modeling kinetically controlled RD processes. In their work, they described the column start-up and disturbance rejection behavior during continuous operation. They found that the calculated dynamic temperature profiles were similar to their experimental data, but the concentration profiles differed from them. They also suggested that tray hydraulics be taken into consideration in dynamic models. Monroy-Loperena and Alvarez-Ramı´rez9 analyzed the multiplicity maps for the ethylene glycol column and suggested a possible start-up procedure, based on the shape of the bifurcation diagram. Bisowarno and Tade´22 used SpeedUp to simulate several start-up scenarios for the ethyl tert-butyl ether column, while analyzing the potential problems that were due to input multiplicities embedded in the system. Scenna and Benz,23 using HYSYS, simulated the startup of the closed-loop ethylene glycol column. Although they do not give details on the controllers used, they successfully analyzed different start-up policies and suggested the introduction of a feed below the reactive trays to reduce start-up time. Reepmeyer et al.24 described “time optimal” strategies for column startup, although they did not solve any optimization problems. In reality, using gPROMS, and a rigorous tray-by-tray model, they simulated different start-up heuristic strategies to determine which one required less time. They do not recommend the use of start-up strategies described for conventional distillation for RD, because they may lead to longer start-up times. Compared to the number of publications related to the design, simulation, and control of RD columns, the number of works devoted to compute optimal start-up and steady-state transition trajectories is rather scarce. Ruiz et al.25 were among the first to use optimization numerical techniques to address the computation of optimal operating policies. Cervantes and Biegler26 used the simultaneous dynamic optimization formulation to minimize the heat required to start a RD column. Recently, Raghunathan et al.27 formulated and solved a hybrid optimal start-up problem where modeling switches were used to take care of different types of models, because of different operating conditions emerging during dynamic transition operations. Dynamic optimization of continuous, batch, and semi-batch processes has been proven to be a useful tool for the enhancement of the dynamic performance of complex chemical processes. Time reductions in the transition between steady states obtained with this method can be very large, thus resulting in

Figure 1. Reactive distillation (RD) column configuration.

a reduction of waste production and energy consumption. Although some operating issues have been addressed, mainly during start-up operations, there are some aspects that have not received consideration. Most of the published works only address the effect of one manipulated variable on the optimal transition trajectory. In this work, we have considered the use of several manipulated variables, therefore increasing the number of available degrees of freedom, as a way to enhance the startup and steady-state transition period measured in terms of reducing the underlying transition times. Moreover, bifurcation diagrams have been used to get a clear idea about the nonlinearities embedded in the model that could lead to potential control problems. In this paper, the dynamic optimization problem is transformed to a nonlinear programming (NLP) problem, using the simultaneous dynamic optimization technique. This approach is also called direct transcription, because both state and control variables are discretized, leading to a large-scale NLP problem.28 This has the benefit that the solution of the NLP problem is coupled with the solution of the algebraic and differnetial equations (DAEs), and, therefore, the DAE is solved only once, at the optimal point. This method has been successfully applied to different continuous, batch, and semibatch processes in the works of Cervantes and Biegler,26 Cervantes et al.,29 and Lemoine-Nava et al.30 2. Mathematical Model The RD column design was approached using the mixedinteger nonlinear and disjunctive programming techniques devised by Jackson and Grossmann.16 The tower consists of 27 trays with multiple feeds (trays 12, 13, 17, and 19; see Figure 1) in which the reaction that occurs is the metathesis of 2-pentene to form 2-butene and 3-hexene:31

2C5H10 T C4H8 + C6H12

(1)

The reaction occurs at atmospheric pressure, and its thermodynamic behavior can be represented with good accuracy, using ideal vapor-liquid equilibrium.32 Table 1 presents a list of the kinetic and physical properties data of the system. In this section, the mathematical model that governs the dynamic behavior of this process is described. The model

2094

Ind. Eng. Chem. Res., Vol. 46, No. 7, 2007

Condenser liquid split fraction:

Table 1. Kinetic and Physical Properties for the Pentene-Butene-Hexene System property rate

(

(h-1)

value/expression

0.5kf xC5 2

)

xC4xC6

Keq 27633 kJ/kg-mol 213216 exp RT -∆G exp RT

(

kf (h-1)

(

Keq

)

)

consists of a tray-by-tray application of the so-called MESH equations. They consist in the mass balance (M), thermodynamic equilibrium (E), composition summations (S), and enthalpy balance (H) equations. The following modeling assumptions are made: negligible vapor-phase dynamics, nonconstant molar holdups, liquid hydraulics modeled by the Francis weir equation, thermodynamic equilibrium between phases, and a total condenser. Because, during start-up operations, the concentration of reactants is high all over the column, chemical reaction was considered in all the trays, the reboiler, and the condenser. If the trays are numbered from bottom to top, then the model is represented by eqs 2-7, where Nc, Nt, and Nr are the set of components, stages (including the reboiler and condenser), and the number of chemical reactions, respectively. Component mass balance (M):

dni,k dt

) zi,kFk +

ni,k+1 ni,k-1 Lk+1 + Ki,k-1Vk-1 Nk+1 Nk-1

NR ni,k (Lk + Ki,kVk) + Wkl νj,iξj,k Nk j)1

∑

i ∈ Nc, k ∈ Nt (2)

Thermodynamic equilibrium (E):

yi,k ) xi,kKi,k )

ni,k K Nk i,k

D L0 + D

R)

i ∈ N c, k ∈ Nt

(3)

(8)

where ni,k is the individual molar holdup of component i in stage k; Nk is the total molar holdup in stage k; zi,k is the composition of component i in feed stream k; Fk is the molar flow rate of feed stream k; Lk and Vk are, respectively, the liquid and vapor molar flow rates of stage k; D is the distillate flow rate; L0 is the refluxed molar flow rate; R is the condenser liquid split fraction (hereafter, this will be called the “reflux split fraction”); τk is the time constant of the linearized Francis weir equation for stage k; Ki,k is the thermodynamic equilibrium constant for component i in stage k; Wjl is the reaction volume at stage k; Wkw is the volume capacity of stage k; νj,i is the stoichiometric coefficient of component i in reaction j; ξj,k is the reaction rate of reaction j at stage k; hi,kl and hi,kV are, respectively, the liquid and vapor enthalpies of component i at stage k; hi,kF,l is the enthalpy of component i in feed stream k; and, finally, ∆Hj,kR is the heat of reaction j at stage k. Also, βk is introduced to model the opening percentage of the valve that allows the flow of bottoms from the reboiler. This variable is useful to model batch operations by setting R ) 0 (total reflux), Fk ) 0, and β1 ) 0. The model applies to all stages with the following modifications: if 1 < k < Nt, then Qk ) 0 and βk ) 1; if k ) 1, then Qk is the reboiler heat duty (QR) and 0 e βk e 1; and, finally, if k ) Nt, then Qk is the condenser heat duty (QC), Lk ) L0 + D, Vk ) 0, and βk ) 1. Note that this model is a system of DAEs of high index (actually, with an index of 2). This can be quickly noticed from the fact that the vapor flow rate is an algebraic variable that does not appear in any of the algebraic equations. Because most solvers for dynamic simulation and the formulation for dynamic optimization require model with an index of 1 (or lower), the index reduction scheme used by Cervantes and Biegler26 was applied here. The first step in the index reduction procedure is to write the left-hand side (LHS) of the energy balance as33

Liquid and vapor-phase composition summations (S): NC

ni,k ∑ i)1

NC

ni,k ) Nk ∑ i)1

k ∈ Nt

(4)

NC

ni,kKi,k ) Nk ∑ i)1

k ∈ Nt

(5)

dhi,kl

∑ i)1

ni,k

dt

Nk+1 i)1

NC

d

i,k+1

F,l

i,k(hi,k

NC

(

xi,k

Vk

i)1

Vk-1 NC

∑n

Nk-1 i)1

i,k-1Ki,k-1(hi,k-1

V

-

l

i,kKi,k(hi,k

- hi,kV) +

k i)1 NR Wkl ξj,k(-∆Hj,kR) j)1

∑

dTk + Qk

k ∈ Nt (6)

Francis weir equation:

Lk )

dxi,k

βk (N - WkwFˆ k) τk k

k ∈ Nt

+

dt

∂Tk

)

d

(1) ) 0 (10)

dt

Solving for dTk/dt, one obtains

NC

∑n N

(9)

dt

k

∂Ki,k dTk

NC

- hi,kl) +

(hi,k+1l - hi,kl) +

hi,kl) +

ni,k ∑ ∂T i)1

( ) ∑( )

NC

i)1

∑n

∂hi,kl dTk

xi,kKi,k) ) ∑ Ki,k ∑yi,k) ) dt(∑ dt i)1 dt i)1 i)1

∑z

Lk+1 NC

)

dt

NC

) Fk

NC

The derivative of the temperature can be obtained from the summation of the vapor compositions as follows:

d

Energy balance (H): NC

dhi,kl

(7)

dt

)-

NC

dxi,k

NC

∂Ki,k

Ki,k ∑ dt i)1

(11) xi,k ∑ ∂T i)1

k

Combining the LHS of the energy balance in eq 9 with its righthand side (RHS) (eq 6) with eqs 11 and 2, one obtains a morecomplex modified algebraic energy balance (eq 12):

[

∂hi,kl

NC

ni,k ∑ ∂T i)1

k

where

-

∑ i)1

(

dni,k

NC

Ki,k

dt

NC

Nk dt

]

)

ni,k dNk

∂Ki,k

ni,k ∑ ∂T i)1

k

Table 2. Mathematical Model Scaling Factors scaling factor

) RHS

Vk-1 NC

∑ni,k-1Ki,k-1(hi,k-1V - hi,kl) +

Nk-1 i)1

Vk NC

variable

NR

ni,kKi,k(hi,kl - hi,kV) + Wkl∑ξj,k(-∆Hj,kR) + Qk ∑ N i)1 j)1 k

Finally, the model comprised by eqs 2, 3, 4, 5, 7, 8, and 12 is the one used for simulations and dynamic optimization, where the new algebraic equation (given by eq 12) is used, instead of the original energy balance, to calculate the vapor flow rates. 3. Simultaneous Dynamic Optimization Formulation The general dynamic optimization formulation for the DAE system can be stated as follows:

min

∫0θ {|| z(t) - zˆ (t) ||2 + || u(t) - uˆ (t) ||2} dt

10 kmol/h 10 kmol 340.3 K 250 kmol/h 250 kmol/h 250 kmol/h 106 kJ/h 40 kmol/h 106 kJ/h

Table 3. Manipulated Variables (QR, R, β, F12, F13, F17, F19), Temperature (T), and Compositions (xC5, xC4, xC6) at the Condenser and Reboiler Used To Track the Column Behavior for Each of the Analyzed Steady States

k+1

hi,kl) +

ni,k Nksfactor Tksfactor Lksfactor Vksfactor Dksfactor Qcksfactor Fksfactor QRsfactor

(12)

zi,k(hi,kF,l - hi,kl) + ∑ ∑ni,k+1(hi,k+1l N i)1 i)1

value

sfactor

Lk+1 NC

NC

RHS ) Fk

-


(13)

s.t. the DAE model:

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

(14)

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

(15)

z(0) ) z0

(16)

Initial conditions:

QR (kJ/h) R β F12 (kg-mol/h) F13 (kg-mol/h) F17 (kg-mol/h) F19 (kg-mol/h) reboiler temperature (K) xC5 at reboiler xC4 at reboiler xC6 at reboiler condenser temperature (K) xC5 at condenser xC4 at condenser xC6 at condenser

steady state 1 3.433 × 106 0.2 1 36.6 34.3 25 24

steady state 2 4.400 × 106 0.2 1 36.6 34.3 25 24

steady state 3 2.694 × 106 0.227 1 36.6 34.3 25 24

338.72

339.58

339.17

0.0187 4.4 × 10-5 0.9813

0.0003 1.5 × 10-8 0.9997

0.0090 1.1 × 10-5 0.9910

276.87

280.02

277.06

10-5 0.9999 1011

0.1491 0.8466 0.0043

0.0099 0.9901 1.8 × 10-5

the value of the derivative in element i at the collocation point q, and Ωq is the polynomial of order Ncol, satisfying

Ωq(0) ) 0

(for q ) 1, ..., Ncol)

Ω′q(Fq,r) ) 0

(for q, r ) 1, ..., Ncol)

where Fr is the location of the collocation point r within each element. Continuous profiles are enforced by

Bounds:

(

Ncol

L

z e z(t) e z

zi ) zi-1 + hi

U

∑

Ωq

q)1

xL e x(t) e xU uL e u(t) e uU pL e p(t) e pU

(17)

where F is the vector of the RHS of the differential equations, G is the vector of algebraic equations (of index 1), z is the differential state vector, z0 are the initial values of z, x is the algebraic state vector, u is the control profile vector, and p is a time-independent vector of parameters. The DAE optimization problem is converted to an NLP by approximating the state and control variable profiles by a family of polynomials on finite elements (t0 < t1 < ...< tne ) θ). Here, a monomial basis representation for the differential profile is Ncol

z(t) ) zi-1 + hi

(

∑ Ωq

q)1

t - ti-1 dz hi

dt

i,q

)

(18)

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

t - ti-1 dz hi

dt

)

i,q

(19)

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. In addition, the control and algebraic profiles are approximated using a similar monomial basis representation, which takes the form Ncol

x(t) )

∑Ψq q)1 Ncol

u(t) )

∑

q)1

( ) ( )

Ψq

t - ti-1 hi

xi,q

(20)

ui,q

(21)

t - ti-1 hi

where xi,q and ui,q represent the values of the algebraic equations and control variables, respectively, in element i at collocation point q. Also, Ψq is the Lagrange polynomial of order Ncol that satisfies

Ψq(Fr) ) δq,r

(for q, r ) 1, ..., Ncol)

2096


Figure 2. Continuation diagrams using the reboiler heat duty as the main bifurcation parameter: (a and c) condenser conditions and (c and d) reboiler conditions. The abbreviations “ss1”, “ss2”, and “ss3” represent the first, second, and third steady states, as defined in Table 3.

From eqs 18 and 19, it can be seen that the differential state 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. In addition, the integral objective function is approximated with a Radau quadrature with Ne finite elements and Ncol quadrature points in each element. This leads to the following objective function: Ne

min Φ )

Ncol

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

NC Nt

min

∫0θ{∑ ∑ Ri,kn

||(ni,k)(ni,ksfactor) - ni,kdes||2 +

i,k

i)1 k)1

Nt

∑ k)1

Nt

RkTk||(Tk)(Tksfactor) - Tsdes||2 +

X)

Xunscaled Xsfactor

k

Nt

Nkdes||2 +

RkF ||(Fk)(Fksfactor) - Fkdes||2 + ∑ k)1 k

RQR||(QR)(QRsfactor) - QRdes||2 + RR ||R - Rdes||2 + Rβ ||β - βdes||2} (24)

(22)

In this work, the objective function used is shown below (eq 24). Note that all variables in the model were substituted by scaled versions of them (e.g., eq 23), so that each term in the objective function has the same order of magnitude. The scaling factors used are shown in Table 2.

RkN ||(Nk)(Nssfactor) ∑ k)1

Here, θ is the transition horizon, all variables with the superscript “des” refer to the unscaled desired values (i.e., the new set point), and the weight factors for each deviation are denoted by R. The values of the weight factors were selected to give importance to certain terms of the objective function during the solution of the optimization problem, and also, to obtain smooth solution profiles of the manipulated and state variables.

(23) 4. Results and Discussion Xunscaled

where X represents any variable to be scaled and is the unscaled variable, whereas Xsfactor is the scaling factor, as given in Table 2. Therefore, in terms of scaled variables, the objective function reads as follows:

In this section, continuation diagrams, the transition between steady states, and start-up operations are shown and discussed. In all cases, the feed streams considered are pure pentene at its boiling point with the following values: F12 ) 36.6 kg-mol/h,


Figure 3. Continuation diagrams using the reflux split fraction as the mainbifurcation parameter: (a and c) condenser conditions and (c and d) reboiler conditions. Table 4. Objective Function Weights Used for Start-Up Calculations weight

i

Rk

Tk

RkNk Rkni,k RkFk RQR RR Rβ

C4 C6

k

value

26 4 26 4 26 4 12, 13, 17, 19

105 105 105 105 105 105 103 103 103 103

weight

i

Rk

Tk


C4 C6

F13 ) 34.3 kg-mol/h, F17 ) 25 kg-mol/h, and F19 ) 24 kgmol/h. The steady states considered for the transitions and startups are characterized by the manipulated variable and output tracking variables values that are shown in Table 3. Although most of the system states were part of the objective function, the dynamic tracking results are only shown in terms of the two main product streams (i.e., mole fractions and temperature at both the condenser and the reboiler). In the first subsection, the continuation diagrams are analyzed; in the second subsection, the start-up operations are examined; and in the third subsection, the transitions between steady states are addressed. IPOPT34 was used to obtain the solutions of the large-scale NLPs that resulted from applying the simultaneous dynamic optimization (SDO) formulation for each of the startups and transitions.

k

value

29 1 29 1 29 1 12, 13, 17, 19

105 105 105 105 105 105 103 103 103 103

weight

i

Rk

Tk


C4 C6

k

value

22 8 22 8 22 8 12, 13, 17, 19

105 105 105 105 105 105 103 103 103 103

Because the solution of dynamic optimization problems by simultaneous techniques have a tendency to be sensitive to the initialization strategy, it deserves to be discussed. Each one of the dynamic optimization problems solved were initialized in the same manner. The first optimization problem solved for each case was initialized using the dynamic ramp simulation results, taking the value of each of the states at the times corresponding to each of the collocation points used, and a linear interpolation of the manipulated variables. Whenever any of the calculations done required adjustments (i.e., if the transition time needed to be reduced to soften the curves), then the previous optimal solution was used to initialize the new problem. 4.1. Nonlinear Steady-State Behavior. Figure 2 shows the nonlinear steady-state diagrams of the light component (butene)

2098


Table 5. Problem Sizes, Number of Iterations, and CPU Time Required to Attain an Optimal Solution steady state transition

Ne/Ncol

no. of variables

no. of restrictions

no. of iterations

CPU time (h)

start-up to 1 start-up to 2 start-up to 3 1 to 2 2 to 1 2 to 3 3 to 2 1 to 3 3 to 1

20/3 20/3 20/3 20/3 20/3 20/3 20/3 20/3 20/3

19533 19533 19533 19173 19173 19233 19233 19233 19233

19253 19253 19253 19113 19113 19113 19113 19113 19113

153 594 189 85 159 35 52 52 250

9.99a 90.01b 13.75a 7.82a 18.05a 2.48a 3.36a 3.13a 14.80a

a Problems are initialized from a previous optimal solution. b Problems are initialized from a dynamic simulation.

mole fraction and the temperature at the condenser (Figures 2a and 2c), and the heavy component (hexane) mole fraction and the temperature of the reboiler (Figures 2b and 2d), using the reboiler heat duty as the main continuation parameter. The nominal operating point is denoted by the symbol “O”, and the labels “ss1”, “ss2”, and “ss3” denote the different steady states (steady states 1, 2, and 3) used in the optimal transitions and startups that have been calculated. The first noticeable aspect of the diagrams is that, near the examined operating region, there are no output multiplicities present and no unstable steady states. An important fact is that, at the condenser, all the nominal

operating points are very close to high-sensitivity regions. This could bring operating problems if the controller is not properly tuned; for example, if the main continuation parameter is increased to 4 × 106 kJ/h, the composition of the light component in the condenser will be reduced to ∼0.9, which is 10% below the nominal value. Also, if the heat duty is reduced to 2.5 × 106 kJ/h, the composition of the heavy component in the bottoms flow rate will be reduced to ∼0.8, almost 20% below the nominal setpoint. On the other hand, input multiplicities do exist. In the case of the condenser temperature, it can be seen that, with a high reflux split fraction (curve 3), input multiplicities exist at low values of the main continuation parameter, and the same is true for the other two values of the reflux split fraction for values of QR below the nominal point. This behavior could lead to operating problems if this state is used to track the progress of this process. At the reboiler, the temperature behaves as expected: it increases monotonically with the heat duty until it reaches the boiling temperature of the hexene. Also, in the internal trays, input multiplicities exist for the compositions of the different components as well as for the temperatures; these are not shown, because of space limitations. In Figure 3, the multiplicity maps of the butene mole fraction (Figure 3a) and the temperature (Figure 3c) at the condenser, as well as the hexene mole fraction (Figure 3b) and the temperature (Figure 3d) at the reboiler, are shown. Note that no input/

Figure 4. Optimal start-up transition to steady state 1. Reboiler data: (- - -) ramp dynamic simulation and (s + s) Optimal transition trajectory. Condenser data: (‚ ‚ ‚) ramp dynamic simulation and (s * s) optimal transition trajectory.


Figure 5. Manipulated variables for start-up to steady state 1.

output multiplicities were obtained at the condenser and reboiler. However, the three addressed nominal operating points are located around highly nonlinear operating regions, which could create potential control problems for a control system that has not been tuned properly. For example, if the reflux split fraction is increased only 10%, the composition of the butene in the condenser is decreased from 1 to 0.55, which is far from the desired purity. Similar behavior is also observed at the reboiler. 4.2. Optimal Start-Up Operations. In the industry, two different column start-up methods are used, and are differentiated on the basis of whether the column is filled or empty at the start-up procedure. Ganguly and Saraf35 suggested that, by starting the column filled with the feed mixture, then the total reflux constraint is reduced, as is the start-up time. Ruiz et al.25 described three steps of a column start-up procedure. The first step is called “discontinuous”, meaning that the column is dry and being filled up. The second step begins exactly at the end of the first one, and is called the “semi-continuous” step. During this period, the column is filled with some mixture and is equivalent to the behavior of a column subjected to large disturbances. Finally, the third step is the “continuous” phase and is similar to maintaining the operation of the column around a desired setpoint. From this observation, it could be said that, by starting the process with a partially filled column, the discontinuous step is skipped, and only the semi-continuous and continuous steps occur. In this work, the discontinuous step is not considered, and the model used requires that the column be filled with a mixture

at its boiling point. Ruiz et al.25 were among the first to formulate and solve start-up problems in RD columns using numerical optimization techniques. The authors used a sequential optimal control procedure. This would require that the DAE system be solved at each iteration of the optimization problem, thus consuming much CPU time. Another problem with this type of formulation could arise if startup or transitions are performed at unstable steady states, because the sequential dynamic optimization method cannot handle them. On the other hand, it has been shown36 that the SDO strategy is well-suited to handle dynamic transitions from and/or to open-loop unstable steady states. The initial conditions used during start-up calculations consider that the column is filled with only pentene at its boiling point. Therefore, the individual molar holdups are all zero, except for the pentene, which is calculated from the tray volume (see the work of Jackson and Grossmann16) and the molar density of the pure component. The total molar holdups are equal to the individual holdup of the pentene, and the temperature of each tray is the bubble temperature of the pure component (310.09 K). The reboiler and condenser holdups were 257.9 and 711.23 kg-mol, respectively. Finally, all the internal flow rates (liquid and vapor) are set to zero. The manipulated variables that were used as decision variables to take care of the start-up procedure were as follows: the reflux split fraction (R), reboiler thermal duty (QR), opening percentage of the valve that allows the flow of bottoms from the reboiler (β), and the flow rate of each feed stream (Fi). As shown in eq 24, all the system states were made part of the objective function and,

2100


Figure 6. Optimal start-up transition to steady state 2. Reboiler data: (- - -) ramp dynamic simulation and (s + s) Optimal transition trajectory. Condenser data: (‚ ‚ ‚) ramp dynamic simulation and (s * s) optimal transition trajectory. Table 6. Comparison of Control Energy Using Both Ramp and Optimal Operating Policies during Column Startup

R [h] ramp optimal β [h] ramp optimal QR [kJ] ramp optimal F12 [kg-mol] ramp optimal F13 [kg-mol] ramp optimal F17 [kg-mol] ramp optimal F19 [kg-mol] ramp optimal ΣF [kg-mol] ramp optimal

to steady state 1

to steady state 2

to steady state 3

5.3130 5.9004

5.8162 5.6075

1.1101 × 101 1.1152 × 101

2.7834 × 101 2.9467 × 101

2.6865 × 101 2.9467 × 101

4.8866 × 101 4.9112 × 101

8.4622 × 108 9.3994 × 107

1.0846 × 109 1.2013 × 108

6.6407 × 108 1.3002 × 108

9.0951 × 103 1.0996 × 103

9.0951 × 103 1.0757 × 103

9.0951 × 103 1.7716 × 103

8.5235 × 103 1.0269 × 103

8.5235 × 103 1.0149 × 103

8.5235 × 103 1.6826 × 103

6.2125 × 103 7.6485 × 102

6.2125 × 103 7.5076 × 102

6.2125 × 103 1.3068 × 103

5.9640 × 103 7.4248 × 102

5.9640 × 103 7.9689 × 102

5.9640 × 103 1.2873 × 103

2.77951 × 104 3.6338 × 103

2.97951 × 104 3.6383 × 103

2.97951 × 104 6.0483 × 103

therefore, were the desired values (after startup is over) of the manipulated variables.

To compare the performance of optimal start-up strategies, the column was also started-up using a heuristic approach. Both



the reboiler thermal duty (QR) and the flow rate of the feed streams (F12, F13, F17, F18) were changed in a ramp-like form. Hence, it was assumed that QR changed from its initial start-up value up to its final steady-state value within 7 h; similarly, all the feed stream flow rates reached their final steady-state values within 2 h. On the other hand, the reflux split fraction (R) and the opening percentage of the bottoms flow rate valve (β) were changed in step-like form. To get smooth and fast start-up profiles, a set of weighting values were used for each one of the cases to be discussed afterward; these values are shown in Table 4. All weights are set to zero, except for those shown. Similarly, Table 5 shows the systems statistics for the addressed start-up and steady-state transitions. 4.2.1. Steady State 1. The optimum start-up profiles to steady-state 1 are shown in Figures 4 and 5. In these figures, the evaluated optimal transition is compared with a ramp startup simulation in the manipulated variables, as it has been previously described. The arrowheads on each curve show which y-axis to use to read them, and, as usual, the lower x-axis is used with the left y-axis; therefore, the upper x-axis is to be used with the right y-axis. As one may notice, the compositions in the condenser depict faster dynamics, compared to the composition dynamics in the reboiler. This behavior can be partially explained recalling that the holdup in the reboiler is greater than that in the condenser. The manipulated variables profiles (Figure 5) show that the column should be operated in a semi-batch mode for ∼2 h. This is observed with feeds at trays 12, 13, and 19, with total reflux, and with no bottoms flow rate. Interestingly, the feed trays that are closer to the

reboiler (12 and 13) are the ones that get the flow rate of the reactant first. Table 6 shows a comparison of the control energy used for starting-up the addressed RD column, using both ramplike and optimal operating policies. As it can be noticed, optimal operating policies lead to an order-of-magnitude energy savings when compared against heuristic control policies. The savings in the amount of raw material are also approximately an order of magnitude. These savings should lead to a more-profitable process and less energy demand (therefore reducing the amount of green house gases (GHGs)), and they clearly show the advantages of using optimal control policies rather than empirical procedures. Overall, the optimal start-up policies seem to require shorter transition times, compared to the heuristic startup policies. 4.2.2. Steady State 2. The optimal start-up profiles for the case in which the set point is steady state 2 are displayed in Figures 6 and 7. In this case, a large reduction in transition time is shown. For example, the composition of the light component in the condenser reaches steady state after ∼5 h of operation, whereas the ramp calculation requires much longer than 25 h. The same can be noticed for the other compositions and the temperature at the condenser. On the other hand, the compositions and temperature at the reboiler require about the same amount of time as the ramp calculation (∼5 h). This could be explained by the fact that the internal liquid hydraulics are modeled, and, therefore, a time decoupling effect between the condenser and the reboiler exists. A fast change in the reboiler heat duty (Figure 7c) causes a quick response in the reboiler. This requires a longer time to be displayed in the condenser

2102


Figure 8. Optimal start-up transition to steady state 3. Reboiler data: (- - -) ramp dynamic simulation and (s + s) Optimal transition trajectory. Condenser data: (‚ ‚ ‚) ramp dynamic simulation and (s * s) optimal transition trajectory.

than in the ramp simulations. In the calculated optimum operating conditions, this is countered by quickly increasing the reflux split fraction above the nominal point, thus increasing the amount of mass flowing through in the upper trays and incrementally increasing the heat transfer and temperature of these stages. This is aided by not allowing the flow of bottoms by fixing β ) 0 at the condenser, and by setting high feed flow rates during the first hours of operation. It should be stressed that, similar to the past case, as shown in Table 6, energy and raw material savings of ∼1 order of magnitude are achieved using optimal start-up control policies. 4.2.3. Steady State 3. Finally, the last case considered is the startup to steady state 3, which is shown in Figures 8 and 9. In this case, the condenser responds quickly enough in the ramp calculation, such that the optimum solution closely resembles it. The reboiler requires a longer time to reach the steady state in the ramp simulation, because of the tray liquid hydraulics. The optimum transition for the reboiler manages to reach the setpoint in a time period almost 10 h less than the ramp evaluation. This is accomplished by completely opening the feed valves allowing high feed rates and increasing the reboiler heat duty. The other manipulated variablessthe opening percentage of the bottoms flow rate and the reflux split fractionsare changed in an almost steplike fashion. As shown in Table 6, energy savings obtained using optimal control profiles are not

as large as they were in the past two cases. However, the energy savings are still acceptable. Moreover, the raw material savings are again ∼1 order of magnitude. From Table 3, we see that start-up times turn out to be similar, except for starting the column to the second steady state. Startup times to the first and third steady states are similar, because the operating conditions for these steady states are not very different. The reason the start-up time to the third steady state is larger than that for the first steady state is because the hexane composition at the reboiler is slightly greater for the third steady state, when compared to the same composition at the first steady state. However, the situation for the second steady state appears different: here, the hexane composition at the reboiler is similar, compared to the same composition at the third steady state. The startup to the second steady state is complicated because the C4 condenser purity is rather low. This means that not enough hexane is formed to meet reboiler specifications and, therefore, a long CPU time is needed to determine an optimal start-up policy able to meet such purity reboiler specifications. Moreover, as the results discussed in this part show, the optimal start-up time to the third steady state is larger, when compared to similar start-up times to the first and second steady states. Overall, the results of this part clearly demonstrate the economic advantages of using optimal control policies when addressing the startup of the examined RD column. The desired



steady states were chosen to meet product purity. Around the steady states, different patterns of dynamic behavior during the start-up period were observed. For instance, the dynamic behavior up to steady state 1 shows that the difference between optimal control and empirical operation policies is not large, although the energy and raw material savings shown in Table 6 clearly demonstrate that optimal control policies lead to important savings. The optimal dynamic start-up policies up to steady state 2 show a sluggish response in the condenser dynamic response and a relative faster response in the reboiler. The opposite behavior is observed when analyzing the column start-up behavior up to steady state 3. The partial conclusion of this subsection is that optimal control policies lead to important energy and raw material savings, when compared to industrially used empirical start-up procedures. 4.3. Transitions between Steady States. Optimal transition policies between the desired steady-state operating points were sought. As observed previously, to get a clear idea about the potential advantages of using optimal transition trajectories, these results were compared with those obtained from using ramp and steplike control policies. The manipulated variables are the thermal duty in the reboiler and the reflux split fraction. The objective function weights used for the below discussed transitions between the steady states are presented in Table 7. As mentioned previously, the reboiler thermal duty was implemented using a 3.5-h ramp. Moreover, the changes in the reflux split fraction were implemented in a steplike manner. 4.3.1. Steady State 1 S Steady State 2. In Figure 10, the transition between steady state 1 to steady state 2 is shown. In

Table 7. Objective Function Weights Used for Transitions between Steady States Transition 1 S 2 weight

i

RkTk RkNk Rkni,k RQR RR

C4 C6

Transition 2 S 3

k value weight 22 7 22 7 22 7

103 103 103 103 103 103 1

i

k

value

weight

5 × 107 5 × 107 5 × 107 5 × 107 5 × 107 5 × 107 5 × 104 5 × 104

RkTk

C4 C6

22 7 22 7 22 7

RkTk RkNk Rkni,k RQR RR

Transition 3 S 1 i

k

value

C4 C6

22 7 22 7 22 7

5 × 107 5 × 107 5 × 107 5 × 107 5 × 107 5 × 107 5 × 104 5 × 104

RkNk Rkni,k RQR RR

this case, the only manipulated variable considered is the reboiler heat duty, and the optimum transition is compared with dynamic simulation, considering a ramp change in the manipulated variable. Figure 10 shows that a great improvement in transition time is obtained for the variables in the condenser, whereas the ramp change requires ∼45 h to achieve the same transition. In contrast, the reboiler (Figure 10) manages the optimal transition within ∼5 h, whereas the step change is completed within ∼15 h. The time-response difference between the condenser and the reboiler is partially explained by the fact that tray hydraulics is modeled. The optimal heat duty policy shown in Figure 10e consists of increasing the heat duty above the nominal point for ∼9 h and then reducing it to the set-point value, thus increasing the temperatures and vapor flow rates in the upper trays quickly, which reduces the residence time in each tray and causes the lower purity of light component in the condenser.

2104


Figure 10. Compositions and temperature dynamic response during transition from steady state 1 to steady state 2. Reboiler data: (- - -) ramp dynamic simulation and (s + s) Optimal transition trajectory. Condenser data: (‚ ‚ ‚) ramp dynamic simulation and (s * s) optimal transition trajectory.

The return path from steady state 2 to steady state 1 is shown in Figure 11. First, note that, for nonlinear systems, the transition from one steady state to another does not necessarily require the same time and trajectory as those of the return transition. The condenser dynamic response for the first transition (from steady state 1 to steady state 2) required ∼45 h, whereas the reversed transition required ∼30 h. From Figure 11, it can be

noticed that the condenser requires less time than the reboiler to complete the change. This could be explained by looking at the reboiler heat duty (Figure 11e), which is decreased initially, causing a lower flow rate of vapor to the top of the column. Because the vapor dynamics is not modeled, this change passes through the remaining trays very fast, quickly reducing the condenser temperature, thus reaching the setpoint. On the other



hand, the reboiler, which is full of liquid, requires a longer period to reach its setpoint, because the increase in the manipulated variable is completed slowly. Almost 40 h are required for it to attain its nominal value. The ramp change simulation is barely visible, because the time period required to attain the setpoint is almost 10 times longer is that for the optimal trajectory. For example, in Figure 11d, a dashed, almost-horizontal, line can

be observed at a temperature of 339.5 K, which corresponds to the ramp dynamic simulation. It appears to be straight, because it has a very small negative slope, and it will eventually reach the desired setpoint. From Table 8, we notice that significant energy savings can be achieved using optimal control policies rather than empirical procedures. Actually, the energy savings are greater for the steady state 2 f steady state 1 transition

2106


Figure 12. Compositions and temperature dynamic response during transition from steady state 1 to steady state 3. Reboiler data: (- - -) ramp dynamic simulation and (s + s) Optimal transition trajectory. Condenser data: (‚ ‚ ‚) ramp dynamic simulation and (s * s) optimal transition trajectory. Table 8. Comparison of Control Energy Using Both Ramp and Optimal Operating Policies during Steady-State Transitions Transition 1 w 2 R [h] ramp optimal QR [kJ] ramp optimal

4.3831 × 108 2.0342 × 108

Transition 2 w 1

1.7182 × 109 1.8392 × 108

Transition 1 w 3

Transition 3 w 1

Transition 3 w 2

Transition 2 w 3

4.9172 4.5433

3.5806 4.0009

2.9925 2.9517

6.5240 6.6935

1.0789 × 109 5.3692 × 107

1.7152 × 109 6.6691 × 107

2.6101 × 108 6.5208 × 107

6.7649 × 108 7.8999 × 107



than for the steady state 1 f steady state 2 transition, stressing the highly nonlinear behavior of the RD column. 4.3.2. Steady State 1 S Steady State 3. Figure 12 depicts the case when the transition is done from steady state 1 to steady state 3. In this case, two manipulated variables are used: the reflux split fraction and the reboiler heat duty. Again, here, a large reduction in transition time can be achieved, in comparison

to the ramp dynamic simulation results. In the condenser (Figure 12), the optimal transition is completed within ∼10 h, whereas the ramp dynamic simulation results display a sluggish response. Also, the temperature change in this section of the column is almost negligible, because of the column design. The trays just below the condenser have compositions similar to that of the distillate and, therefore, similar operating temperatures. Simi-

2108



larly, in the reboiler (Figure 12), the optimal transition is achieved well before the ramp simulation; however, in this case, the ramp calculation can barely be observed, because it shows an inverse response (i.e., the temperature is first reduced before it reaches the nominal value above the starting point). Although the increase in temperature in the reboiler with a lower heat duty seems counterintuitive, it can be explained, because of the

increase in the reflux split fraction. This can be observed when moving from curve 1 to curve 3 in the continuation diagram (see Figure 2d). The optimal profiles for these variables first increase the heat duty and the reflux split fraction, and then they reduce them to their new operating values. This behavior is not obvious from the bifurcation diagrams, where the new steady state has lower values for both manipulated variables



than those for the starting point, and this can be expected of optimal trajectories. The optimal operating trajectory between steady state 3 and steady state 1 is shown in Figure 13. In this figure, the behavior of the output variables at the condenser can be observed, and a great reduction in transition time is displayed. The ramp dynamic simulation requires >20 h, whereas the optimal

transition only requires ∼2 h. The reboiler dynamic response displayed in Figure 13c shows that the optimal transition time again is much smaller than the time required for the ramp simulation. In this section of the column, the setpoint is reached after 8 h by the optimal trajectory, whereas the ramp change requires a time period much longer than 20 h. In fact, in Figure 13c, the horizontal dashed line at the top represents the

2110


ramp dynamic simulation response, which has a very small negative slope that will reach the setpoint after a very long time. In Figures 13a-c, the ramp change cannot be observed, because the outputs present an inverse response (i.e., in Figure 13c, the composition of hexene first increases to 1 before it decreases to the setpoint of ∼0.98). The manipulated variables shown in Figure 13 indicate that (i) the column should be operated at total reflux (Figure 13f) for ∼2 h, and (ii) the reboiler heat duty (Figure 13e) should be first reduced, allowing the temperatures to decrease quickly along all the stages, therefore allowing the condenser to reach its setpoint fast, and, finally, it is increased until the reboiler reaches its new operating temperature. Table 8 shows that the 1 S 3 steady-state transitions feature the largest energy savings. As a matter of fact, energy savings of ∼2 orders of magnitude can be achieved using optimal control trajectories, instead of ramplike empirical control policies. 4.3.3. Steady State 2 S Steady State 3. Figure 14 depicts the transition between steady state 2 to steady state 3. Here, the reflux split fraction and the reboiler heat duty were both used as manipulated variables. The optimal transition at the condenser requires ∼5 h, whereas the ramp simulation response requires >30 h. At the reboiler, similar dynamic responses are achieved using both optimal and empirical control policies. In this case, the reboiler (Figure 14) responds quickly, because the change in its setpoint is very small, and the decrease in temperature is achieved by a fast reduction in the reboiler heat duty. In contrast, the changes in the compositions and temperature of the condenser (Figure 14) are larger and require a longer time period. This could be attributed to the fact that the variation in the reboiler heat duty is greater than that of the reflux split fraction, and, because the liquid tray hydraulics are modeled, the changes in the first parameter require a longer time period to reach the condenser. The transition from steady state 3 to steady state 2 is shown in Figure 15. Again, in this case, the transition is achieved faster at the reboiler, because the change in temperatures and compositions in this stage is small, compared to the changes in the condenser. In fact, Figure 15 shows that the changes in the reboiler temperature are

Optimal Start-Up and Product Transition Policies of a Reactive

Recommend Documents