A Case Study in Simultaneous Design and Control Using Rigorous

760

Ind. Eng. Chem. Res. 2002, 41, 760-778

A Case Study in Simultaneous Design and Control Using Rigorous, Mixed-Integer Dynamic Optimization Models Vikrant Bansal,† John D. Perkins, and Efstratios N. Pistikopoulos* Centre for Process Systems Engineering, Department of Chemical Engineering, Imperial College, London SW7 2BY, U.K.

This article presents an in-depth case study that demonstrates (i) how the process and control design, involving both discrete and continuous decisions, can be simultaneously optimized for systems described by realistic dynamic models and that (ii) mixed-integer dynamic optimization problems, involving thousands of differential-algebraic equations, can now be solved using stateof-the-art algorithms and technology. The study in question involves a distillation column that is subject to time-varying disturbances and time-invariant uncertainty. A new, multicomponent, mixed-integer dynamic distillation model is developed and used to simultaneously determine the number of trays, feed tray location, column diameter, reboiler and condenser surface areas, 5 × 5 control structure, PI-gains, reset times, and set-points, to obtain an economically optimal system that is also guaranteed to give feasible dynamic operation over all uncertainty scenarios. 1. Introduction Almost 60 years ago, Ziegler and Nichols1 observed that the ability to control a system in the face of uncertainties and disturbances is dependent on its process design, and as such, systematic methods are needed that account for operability issues at the earliest stages of the process design procedure. Over the past 20 years or so, this realization has been reawakened in both academia and industry, leading to a glut of publications in the field. Much of the integrated design and control literature (for a description of major works, see Pistikopoulos and van Schijndel2 and Bansal3) has concentrated on the development of analysis tools that provide some measure of a system’s controllability and that allow design alternatives to be compared on a common basis. The principal, common drawback of these approaches is that it is usually unclear what implications the value of the particular controllability metric has for changes to the process design; this leads to design changes being “postulated on physical grounds in an ad hoc manner”.4 Furthermore, although these tools have been developed to give an indication of the closed-loop behavior of a system, how they actually relate is again uncertain. This is especially true, since most of the tools rely on the use of steady-state or linear dynamic models. While the use of such models may be adequate in some cases, in general it is quite unpredictable whether the conclusions drawn are correct, particularly in the face of process nonlinearities, disturbances, uncertain parameters, and constraints. That this is the case is evidenced by the large number of papers in which the authors “verify” their findings through some form of closed-loop, nonlinear dynamic simulation. A number of more systematic procedures have been developed in recent years for designing systems that are * To whom all correspondence should be addressed. Telephone: +44 (0)20 7594 6620. Fax: +44 (0)20 7594 6606. E-mail: [email protected]. † Current address: Process Systems Enterprise Ltd., Bridge Studios, 107a Hammersmith Bridge Road, London W6 9DA, U.K.

Figure 1. Binary distillation system for the design and control study.

economically optimal and operable. However, despite these advances in integrated design methodologies, limitations still exist. To illustrate this point, consider the binary distillation column shown in Figure 1, which is adapted from one presented by Viswanathan and Grossmann.5 Here, a mixture of benzene and toluene is to be separated into a top product with at least 98 mol % benzene and a bottom product with no more than 5 mol % toluene. The system is subject to uncertainty in the feed flow rate and the cooling water inlet temperature (where the latter can be described dynamically by a slow sinusoid representing diurnal, ambient variations), as well as a high-frequency sinusoidal disturbance in the feed composition. The objective is to design the distillation column and its required control scheme at minimum total annualized cost (comprising capital costs of the column, internals, and exchangers, and operating costs of the steam and cooling water utilities), capable of feasible operation over the whole of a given time horizon, where feasibility is defined through the satisfaction of constraints such as (i)

10.1021/ie010156n CCC: $22.00 © 2002 American Chemical Society Published on Web 01/22/2002

Ind. Eng. Chem. Res., Vol. 41, No. 4, 2002 761 Table 1. Comparison of Integrated Design Methods versus the Desired Features discrete decisions authors

model

typea

Lennhoff and Morari (1982)41 Palazoglu and Arkun (1986, 1987)42,43 Abbas et al. (1993)44 Luyben and Floudas (1994a,b)45,46 Brengel and Seider (1992)47 Nishida and Ichikawa (1975)48 Nishida et al. (1976)49 Walsh and Perkins (1994, 1996)6,7 Rowe et al. (1997)50

NLD LD SS SS NLD NLD NLD NLD NLD

Mohideen et al. (1996, 1997a)8,34 Schweiger and Floudas (1997)30 Bahri et al. (1997)51

NLD NLD NLD

a

uncertain parameters

disturbances

process

control

x

x x x x

x x x x x x

0-1

0-1

relaxed 0-1 relaxed 0-1 heuristics heuristics

heuristics heuristics

0-1 0-1 0-1

0-1 0-1

SS ) steady-state. LD ) linear dynamics. NLD ) nonlinear dynamics.

product quality specifications; (ii) minimum column diameter to avoid flooding; (iii) fractional entrainment limit; (iv) temperature driving forces in the reboiler and condenser; (v) limit on the heat flux in the reboiler; (vi) limit on the cooling water outlet temperature; (vii) above atmospheric pressure operation for the column; (viii) limits on the liquid levels in the reflux drum and reboiler; and (ix) limits on the flow rates of the steam and cooling water. Solution of the problem thus requires the determination of (i) the optimal process design, in terms of the number of trays and feed location (discrete decisions), and the column diameter and condenser and reboiler surface areas (continuous decisions); and (ii) the optimal control design, in terms of the pairings of manipulated and controlled variables (discrete decisions), and the tuning parameters for the given control structure (continuous decisions). While, at first glance, the task of designing this column and its associated control system may seem simple, it is apparent from the problem statement above that, in fact, any methodology that is used to solve the problem must (i) be applicable to nonlinear, dynamic systems; (ii) adequately treat and ensure feasibility in the presence of both time-varying disturbances and time-invariant (or relatively slowly varying) uncertainties; (iii) involve selection of the best process design, taking into account both discrete and continuous decisions; and (iv) involve selection of the best control scheme, again incorporating structural and continuous decisions. Table 1 cross-checks the principal integrated design methods reported in the literature with the criteria above. It can be seen that there are, in fact, only two approaches that have all the required featuressthose of Walsh and Perkins6,7 and Mohideen et al.8 These two approaches are similar in spirit; however, the approach of Walsh and Perkins uses heuristics for the discrete optimization decisions that are more tailored toward the wastewater neutralization processes that they were considering. Thus, the framework of Mohideen et al. is the only integrated design approach that has all the conceptual features required for selecting the best process design and control scheme for general, nonlinear dynamic processes where feasibility must be guaranteed in the presence of both disturbances and uncertainties. Nevertheless, there are a number of issues that can be highlighted related to the application of this framework. First, it requires the solution of mixed-integer dynamic optimization (MIDO) problems, which is, in general, a very difficult task. Second, because of the lack of

efficient methods for accomplishing this, the examples to date have used process and control model simplifications,8 or rigorous models with fixed discrete decisions and simplifications in the treatment of uncertainty (the double-effect distillation system of Bansal et al.;9 the industrial two-column system of Ross et al.10). Recently, much progress has been made in the development of efficient MIDO algorithms.3,11 The objective of this article is to revisit the distillation problem shown in Figure 1 and present a detailed case study that demonstrates how such algorithms allow the framework of Mohideen et al. to be used to solve problems where realistic models are employed. This case study also illustrates how MIDO problems can now be readily tackled that involve thousands of differential-algebraic equations. The remainder of this article is organized as follows. In the next section, a rigorous, mixed-integer dynamic model is developed that can be used for optimizing the process and control structure, design, and dynamic operation of a general, multicomponent distillation system. New features are highlighted with particular reference to the incorporation of binary variables related to the number of trays and the treatment of dynamic path constraints. The framework of Mohideen et al. is then applied in a transparent, step-by-step manner for the case study in question. Issues are described that relate to the solution of the MIDO problem, and the results of the study are discussed in detail. 2. Modeling Aspects Due to the complexity and highly constrained nature of the problem stated in the previous section, it is likely that a simplified dynamic model using “traditional” assumptions (such as constant molal overflow, constant relative volatility, constant liquid, and negligible vapor hold-ups) will be inadequate for realistically portraying the operability characteristics of the distillation system over time. We thus choose to develop a rigorous, mixedinteger dynamic model that relaxes many of the assumptions used in distillation studies that have appeared in the literature. Some of the key features that set this model apart from many previous works are (i) dynamic material and energy balances for each tray, reboiler, condenser, and reflux drum; (ii) the consideration of both liquid and vapor, material and energy holdups; (iii) the consideration of tray nonideality through use of Murphree tray efficiencies; (iv) the consideration of liquid hydraulics and level on each tray using

762

Ind. Eng. Chem. Res., Vol. 41, No. 4, 2002

with each other, and that the dynamics of the downcomers can be neglected. The material and energy balances presented below do not consider entrainment, weeping, draw-offs, or external heat inputs for the trays, though these can be easily incorporated if required. 2.2.1. Fundamental Equations

Feed conditions: NC

zi,f ) 1 ∑ i)1

(4)

hf ) hf(Tf,zf)

(5)

Feed flow rate to each tray: Fk ) F‚yfk, k ) 1, ..., N

(6)

Component molar balances: N

∑

dMi,k

) Lk+1xi,k+1 + Vk-1yi,k-1 + Fkzi,f + dt Rkxi,d - Lkxi,k - Vkyi,k, i ) 1, ..., NC, k ) 1, ..., N (7) Energy balances:

(

yrk′)

k′)k

Figure 2. Distillation column superstructure.

modified Francis weir formulas; (v) equations for the pressure drop from tray to tray by considering the pressure drop incurred by the vapor in flowing through the openings at the bottom of each tray and the hydrostatic pressure on each tray; and (vi) detailed flooding and entrainment calculations and subsequent evaluation of bottlenecks in the column and the minimum allowable column diameter. The specific details of the model are discussed in the following sections. All the relevant units and case-studyspecific parameter values can be found in Appendix A. 2.1. Process Superstructure. A similar superstructure to that presented by Viswanathan and Grossmann12 is used in order to represent a column where the number of trays and the feed tray location are to be optimized (Figure 2). It is assumed that a reasonable estimate is available of the upper bound, N, on the number of trays, while binary variables, yfk and yrk, k ) 1, ..., N, are introduced into the model to account for the locations of the feed and reflux, respectively. Thus, yfk ) 1 if all the feed enters tray k; similarly, yrk ) 1 if all the reflux enters tray k. The following constraints are used in order to ensure that only one tray receives the feed and one tray receives the reflux:

N

(

∑

k′)k

yrk′)

dUk dt

l v + Vk-1hk-1 + Fkhf + Rkhld ) Lk+1hk+1

Lkhlk - Vkhvk, k ) 1, ..., N (8) Component molar hold-ups: Mi,k ) Mlkxi,k + Mvkyi,k, i ) 1, ..., NC, k ) 1, ..., N (9) Energy hold-ups: Uk ) Mlkhlk + Mvkhvk - 0.1Pk‚voltray, k ) 1, ..., N (10) Volume constraints: Mlk Flk

+

Mvk Fvk

) voltray, k ) 1, ..., N

Definition of Murphree tray efficiencies: / yi,k ) yi,k-1 + effi,k(yi,k - yi,k-1), i ) 1, ..., NC, k ) 1, ..., N (12)

N

∑ yfk ) 1 k)1

(1)

N

∑ yrk ) 1

(2)

k)1

k

effi,k ) ai,k + (1 - ai,k)

∑ yrk′, k′)1

i ) 1, ..., NC,

k ) 1, ..., N (13) Equilibrium vapor phase composition: v / l Φi,k yi,k ) Φi,k xi,k, i ) 1, ..., NC, k ) 1, ..., N

The reflux is also constrained so that it can only be placed at or above the feed tray. Hence

(11)

(14)

Mole fractions normalization: NC

xi,k ) 1, ∑ i)1

N

yfk -

∑ yrk′ e 0, k′)k

k ) 1, ..., N

(3)

2.2. Trays. The trays are modeled as sieve trays in a manner similar to that presented by Pantelides13 and Bansal et al.,9 but with the introduction of binary variables in keeping with the superstructure in Figure 2. It is assumed that the liquid and vapor phases are well-mixed and in thermal and mechanical equilibrium

k ) 1, ..., N

(15)

k ) 1, ..., N

(16)

, k ) 1, ..., N

(17)

NC

yi,k ) 1, ∑ i)1 Liquid levels: levelk )

Mlk FlkAtray

Ind. Eng. Chem. Res., Vol. 41, No. 4, 2002 763

Liquid outlet flow rates (modified Francis formula for liquid flow over a rectangular notch weir: Coulson et al.14):

{

0, if levelk e heightweir Lk ) 1.84Flk‚lengthweir(levelk - heightweir)1.5 × 60, k ) 1, ..., N, otherwise (18) Pressure driving force for vapor inlet15: Pk-1 - Pk ) 10-5(

2 v yrk′)(R‚velk-1 F˜ k-1 + βF˜ lkg‚levelk), ∑ k′)k

(

)

Vk-1 1 , k ) 1, ..., N v 60 F A k-1 holes

(20)

2.2.2. Geometry

Free volume between trays: voltray ) space‚Atray

(21)

Active area of a tray:15 Atray ) 0.8Acol

(22)

Cross-sectional area of the column: π Acol ) D2col 4

(23)

Total area of all the active holes:16 Aholes ) 0.12Atray

(24)

Weir length (chord of a circle): lengthweir ) 0.7266Dcol

(25)

2.2.3. Flooding and Entrainment Correlations

Fractional entrainment for 80% flooding factor:17 entk ) 0.224 × 10-2 + 2.377 exp(-9.394‚FLV0.314 ), k k ) 1, ..., N (26) Sherwood flow parameter: FLVk )

()

L ˜ k F˜ vk 0.5 , V ˜ k F˜ lk

k ) 1, ..., N

(27)

Mass flow rates: NC

xi,k‚MWi, ∑ i)1

L ˜ k ) Lk

k ) 1, ..., N

(28)

k ) 1, ..., N

(29)

NC

yi,k‚MWi, ∑ i)1

V ˜ k ) Vk

Flooding velocity:15 ) 60

() ( ) σlk 0.2

20

), k ) 1, ..., N (31) exp(-1.463‚FLV0.842 k Minimum allowable column diameter and area: min Dcol,k )

( ) min 4Acol,k π

0.5

, k ) 1, ..., N

(32) (33)

N

Vapor velocities:

velflood k

K1k ) 0.0105 + 0.1496‚space0.755 ×

min min ) 0.9Acol,k , k ) 1, ..., N Anet,k

k ) 1, ..., N (19)

velk-1 )

Empirical coefficient:17

‚K1k‚

F˜ lk

-

F˜ vk

F˜ vk 0.5

Minimum net area for vapor-liquid disengagement: min ) Anet,k

Vk , v Fk‚floodfac‚velflood k

k ) 1, ..., N

(34)

2.2.4. Remarks. At this stage it is worth remarking on the presence of the binary variables, yrk, in several of the tray equations above. To ensure that the properties of the vapor outlet from the top tray N are not different from those of the vapor outlet from the tray receiving the reflux (i) binary terms are included in the material and energy balances (eqs 7 and 8), so that these equations reduce to their steady-state equivalents for the “shadow” trays above reflux; (ii) the Murphree efficiencies of the reflux tray and above are set to 1 via the binary term in eq 13; (iii) the binary term in eq 19 ensures that there is no pressure drop across the “shadow” trays; and (iv) the flow rate of liquid, LN+1, into the top tray is assigned to zero (see Section 2.7). Hence, the dynamic model allows no physical mechanism for any separation to be effected in the vapor stream flowing up through the “shadow” trays, thus making it suitable for determining the optimal number of trays. 2.3. Reboiler. A kettle-type reboiler is considered. The reboiler has a liquid inlet from tray 1 of the column which it partially vaporizes using an external heat input. Vapor is returned to tray 1 while the remaining liquid forms the bottom product from the column. The reboiler effectively acts as another tray, and so its modeling equations are very similar to those presented in Section 2.2. The principal difference is that additional equations are required to model the heating medium.

Component molar balances: dMi,0 ) L1x1 + V0yi,0 - Bxi,0, i ) 1, ..., NC dt Energy balance:

(35)

dU0 ) L1hl1 + V0hv0 - Bhl0 + 60Qreb dt

(36)

Component molar hold-ups: Mi,0 ) Ml0xi,0 + Mv0yi,0, i ) 1, ..., NC

(37)

Energy hold-up: U0 ) Ml0hl0 + Mv0hv0 - 0.1P0‚volreb

(38)

Volume constraint: , k ) 1, ..., N (30)

Ml0 Fl0

+

Mv0 Fv0

) volreb

(39)

764


Volume constraint:

Phase equilibrium: v l Φi,0 yi,0 ) Φi,0 xi,0, i ) 1, ..., NC

Mlc

(40)

Flc


+

Mvc Fvc

) volcond

(56)

Phase equilibrium:

xi,0 ) 1 ∑ i)1

(41)

NC

yi,0 ) 1 ∑ i)1

(42)

v l Φi,c yi,c ) Φi,c xi,c, i ) 1, ..., NC

(57)


xi,c ) 1 ∑ i)1

Cylindrical geometry:

(58)

NC

π volreb ) D2reb‚lengthreb 4

(43)

yi,c ) 1 ∑ i)1

(44)

Lc ) 0, if levelc e heightcondweir 1.84Flc‚lengthcondweir(levelc - heightcondweir)1.5 × 60, otherwise

Liquid level:

(59)

Liquid outlet flow rate: Ml0

)

π Fl0 level20‚lengthreb 4

Heat transfer equations for saturated steam condensing and then subcooling, with negligible hold-up in the coil: Qreb ) Qcondsn reb

)

Qcondsn reb

+

Qsub reb

Ucondsn Scondsn (Ts,in reb reb

Qcondsn reb

(45) - T0)

1 ) Fs∆hvapo (Ts,in) s 60

Ts,in - Ts,out sub sub Qsub reb ) Ureb Sreb Ts,in - T0 ln Ts,out - T0

(

Qsub reb )

(46) (47)

)

(48)

1 F [h (T ) - hw(Ts,out)] 60 s w s,in

(49)

+ Ssub Sreb ) Scondsn reb reb

(50)

dMs ) Fs dt

(51)

2.4. Condenser. The condenser, like a tray, is treated as a vapor-liquid equilibrium system with a weir. It receives a vapor inlet from the top tray, N, which it totally condenses by means of a coolant such as cooling water flowing through a coil. The resulting (saturated) liquid is then fed to the reflux drum.

dMi,c ) VNyi,N - Lcxi,c, i ) 1, ..., NC dt

(52)

Energy balance: dUc ) VNhvN - Lchlc - 60Qcond dt

(53)

Component molar hold-ups: Mi,c )

+

Calculation of vapor inlet flow using equations related to minor head loss due to friction:18

( )( )[

VN ) 60

2 Kloss

0.5

]

1e5(PN - Pc) π 2 Dpipe FvN 4 F˜ vN

0.5

(61)

Cylindrical geometry: π volcond ) D2cond‚lengthcond 4

(62)

π Mlc ) Flc level2c ‚lengthcondweir 4

(63)

Liquid level:

i ) 1, ..., NC

Qcond ) UcondScond

(TN - Tw,out) - (Tc - Tw,in) TN - Tw,out ln Tc - Tw,in

Qcond )

(

)

1 - Tw,in) F cj (T 60 w pw w,out

(54)

Energy hold-up: Uc ) Mlchlc + Mvc hvc - 0.1Pc‚volcond, k ) 1, ..., N (55)

(64)

(65)

Total cooling water consumption: dMw ) Fw dt

Component molar balances:

Mvc yi,c,

(60)

Heat-transfer equations, with negligible hold-up of cooling water:

Total steam consumption:

Mlcxi,c

{

(66)

2.5. Reflux Drum. Like the sieve trays and other units, the reflux drum is a fixed volume vessel with vapor-liquid equilibrium.

Component molar balances: dMi,d ) Lcxi,c - (R + D)xi,d, i ) 1, ..., NC (67) dt Energy balance: dUd l ) Lchc - (R + D)hld dt

(68)


Component molar hold-ups: Mi,d ) Mldxi,d + Mvdyi,d, i ) 1, ..., NC

(69)

Energy hold-up: Ud ) Mldhld + Mvdhvd - 0.1Pd‚voldrum,

k ) 1, ..., N

(70)

Volume constraint: Mld Fld

+

Mvd Fvd

) voldrum

(71)

Phase equilibrium:

Figure 3. Sinusoidal variations for the distillation study.

v l yi,d ) Φi,d xi,d, i ) 1, ..., NC Φi,d

(72)


xi,d ) 1 ∑ i)1

(73)

NC

yi,d ) 1 ∑ i)1

(74)

Reflux flow rate to each tray: Rk ) R‚yrk, k ) 1, ..., N

(75)

Vessel volume (cylinder with hemispherical ends): π π voldrum ) D2drum‚lengthdrum + D3drum 4 6

(76)

Liquid level: Mld ) Fld‚volld volld )

(77)

19 π level2d‚lengthdrum 15 4

hl∪ ) hl(P∪,T∪,x∪), ∪ ) {f, k ) 0, ..., N, c, d}

(78)

) h (P∪,T∪,y∪), ∪ ) {k ) 0, ..., N, c, d}

(79)

v

(2πt 60 )

z1,f ) 0.7 + 0.07 sin

2.6. Physical Properties. The process model eqs 4-78 are independent of the types of models used to describe the physical properties. Hence, the required physical properties can be generically represented by the following equations:

hv∪

The specific physical property models used for the benzene-toluene case study are described in Appendix B. 2.7. Specifications. Equations 4-86 constitute a system of [N(7NC + 27) + 15NC + 58] differentialalgebraic equations (DAEs) in [N(7NC + 27) + 17NC + 71] variables. For simulation purposes, this means that (2NC + 13) variables need to be specified in order for there to be a square system. These correspond to (i) the feed conditions F, Tf, and zi,f for i ) 1, ..., NC - 1; (ii) the inlet temperatures of the utilities Tw,in and Ts,in; (iii) the design variables Dcol, Sreb, and Scond; (iv) the manipulated variables R, D, Fw, Fs, and B; and (v) the spurious superstructure stream LN+1 ) 0 and arbitrary l and xi,N+1 for i ) 1, ..., NC. values for hN+1 In the context of the benzene-toluene case study: (i) The feed flow rate F is an uncertain parameter with nominal value 6 and range 6-6.6 kmol min-1. (ii) The feed temperature Tf ) 363.29 K. (iii) The feed composition of benzene is a highfrequency sinusoidal disturbance with a 60 min time period

(87)

(iv) The cooling water inlet temperature is subject to diurnal ambient variations. At steady-state a designer might consider it as an uncertain parameter with nominal value 298.15 and range 293.15-303.15 K. Dynamically, it can be represented by a slow-moving sinusoid with a 1440 min (24 h) time period

2πt (1440 )

Tw,in ) 298.15 + 5 sin

(88)

l

Fl∪ ) F (P∪,T∪,x∪), ∪ ) {k ) 0, ..., N, c, d} (80) Fv∪

) F (P∪,T∪,y∪), ∪ ) {k ) 0, ..., N, c, d} (81) v

NC

F˜ lk ) Flk

xi,k‚MWi, ∑ i)1

k ) 0, ..., N

(82)

k ) 0, ..., N

(83)

NC

F˜ vk ) Fvk

yi,k‚MWi, ∑ i)1

l Φi,∪ ) Φli(P∪,T∪,x∪), i ) 1, ..., NC, ∪ ) {k ) 0, ..., N, c, d} (84) v Φi,∪ ) Φvi (P∪,T∪,y∪), i ) 1, ..., NC, ∪ ) {k ) 0, ..., N, c, d} (85)

σlk ) σl(Pk,Tk,xk), k ) 0, ..., N

(86)

The feed composition (eq 87) and the cooling water inlet temperature (eq 88) are shown graphically in Figure 3. (v) Ts,in ) 420 K. (vi) The diameter of the column, Dcol, and the areas of the exchangers, Sreb and Scond, are optimization variables. (vii) The manipulated variables, R, D, Fw, Fs, and B, are determined by the control scheme and its tunings. Note that these manipulated variables correspond to the “standard” five for distillation control.19 Also, for fixed feed flow rates, at steady-state it is not possible to independently specify the reflux, distillate, and bottoms flow rates, R, D, and B. In this case, one would specify the values of the liquid levels in the reboiler and the reflux drum, level0 and leveld, respectively, since these correspond to pure integrators. This leaves three manipulated variables at steady-state: R or D, Fw, and Fs.

766


2.8. Initial Conditions. There are [N(NC + 1) + 3NC + 5] differential states in eqs 4-86. If the process is initially at steady-state, then the corresponding [N(NC + 1) + 3NC + 5] initial conditions are

|

dMi,∪ dt

t)0

) 0, i ) 1, ..., NC, ∪ ) {k ) 0, ..., N, c, d} (89)

|

dU∪ dt

t)0

) 0, ∪ ) {k ) 0, ..., N, c, d}

(90)

Ms(t)0) ) 0

(91)

Mw(t)0) ) 0

(92)

2.9. Control Superstructure. In this study, the selection of the optimal multiloop proportional-integral (PI) control scheme is considered using the superstructure in Figure 4 proposed by Narraway.20 A set of binary variables, ycm,n, m ) 1, ..., Nin, n ) 1, ..., Nmeas, is introduced so that ycm,n takes a value of 1 if the mth manipulated variable, um, is paired with the nth measured variable, measn, and is 0 otherwise. The following constraints apply: Nmeas

um ) bm,1 +

∑ Km,n(em,n + Im,n),

m ) 1, ..., Nin (93)

n)1

em,n ) setm,n - measn, m ) 1, ..., Nin, n ) 1, ..., Nmeas (94) dIm,n em,n , m ) 1, ..., Nin, n ) 1, ..., Nmeas ) dt τm,n dbm,n ) 0, m ) 1, ..., Nin, n ) 1, ..., Nmeas dt

(95)

(96)

setm,n ) setm-1,n, m ) 2, ..., Nin, n ) 1, ..., Nmeas (97) L U ‚ycm,n e Km,n e Km,n ‚ycm,n, m ) 1, ..., Nin, Km,n n ) 1, ..., Nmeas (98) Nmeas

∑ ycm,n ) 1,

m ) 1, ..., Nin

(99)

n ) 1, ..., Nmeas

(100)

n)1 Nin

∑ ycm,n ) 1, m)1 with the initial conditions

Im,n(t)0) ) 0, m ) 1, ..., Nin, n ) 1, ..., Nmeas (101) em,n(t)0) ) 0, m ) 1, ..., Nin, n ) 1, ..., Nmeas (102) It should be noted that the modeling equations of Narraway20 have been modified here by the addition of eq 96. This is required in order to enable the distillation

Figure 4. Multiloop PI control superstructure.

system to initialize at steady-state, since the steadystate and dynamic degrees-of-freedom are different, as discussed in Section 2.7. In the context of the distillation model described by eqs 4-86, the number of manipulated variables (inputs) Nin ) 5. From Section 2.7, these correspond to

u1 ) R

(103)

u2 ) D

(104)

u3 ) Fw

(105)

u4 ) Fs

(106)

u5 ) B

(107)

In principle, any number of measured variables, nmeas, could be considered by including, for example, tray temperatures and so on. However, in this study the five common control objectives for a binary system19 are considered. Hence

meas1 ) x1,d

(108)

meas2 ) leveld

(109)

meas3 ) Pc

(110)

meas4 ) x1,b

(111)

meas5 ) level0

(112)

Thus, eqs 93-97 and 103-112 constitute an additional 110 equations to the process in eqs 4-86, with 160 new unknown variables (of which 50 are differential states). The required specifications (that replace the manipulated variables in Section 2.7) are the gains, reset times, and set-points of the controllers, Kmn, τmn, and set1n for m ) 1, ..., 5 and n ) 1,..., 5.


2.10. Feasibility Path Constraints. 2.10.1. Description

Distillate purity: x1,d g 0.98

(113)

x1,b e 0.05

(114)

Bottoms purity:

Prevent subatmospheric operation: Pc g 1.013 25

Figure 5. Tracking the maximum value of a variable over time.

(115)

Liquid level in the reflux drum:19 0.1Ddrum e leveld e 0.9Ddrum

(116)

convert it into a number of point constraints. However, this does not guarantee that the constraint will not be violated between the time points used. The constraint can be more stringently enforced by converting it into an end-point constraint:23

Liquid level in the reboiler (keep the tubes covered but avoid flooding): 0.6Dreb e level0 e 0.9Dreb

dZ ) [max(0,X (t)-X max)]2 dt

(126)

Z(t0) ) 0

(127)

(118)

Z(tf) ) 0

(128)

(119)

As pointed out by Walsh,24 the end-point constraint (eq 128) often leads to slow convergence if implemented within a dynamic optimizer, and so a more practical alternative is to relax the constraint (eq 128) to

(117)

Minimum column diameter due to flooding: min - Dcol e 0, k ) 1, ..., N Dcol,k

Fractional entrainment limit:17 N

(

∑ yrk′)‚entk e 0.1, k′)k

k ) 1, ..., N

Z(tf) e

Temperature driving force for the reboiler: Ts,out - T0 g 0

(120)

Limit on the heat flux to prevent film boiling in the reboiler (The critical heat flux is predicted from the Zuber equation for a single tube (ref 21, p 592), with a scaling factor of 0.5 to account for tube bundles.):

Qreb i)1 - 6.55 × 10-8 x ∆hv (T ) × Sreb NC i,0 i,0 0

[

]

[σl0g(F˜ l0 - F˜ v0)(F˜ v0)2]0.25 e 0 (121)

where is a user-specified tolerance. The potential difficulty in applying eq 129 to a problem with several path constraints is that of how to choose a suitable value of for each constraint, particularly since this will be dependent on the scaling of the variables involved and the length of the time horizon. In this study we employ a different approach, developed by Swartz25 and stemming from some initial work by Cobden.26 The basic idea is to define a variable Xh that tracks the maximum value of X over time, and to then enforce an end-point constraint on this new variable. Hence

Limit on the outlet cooling water temperature to prevent corrosion (ref 22, p 437): Tw,out e 323.15

(129)

Xh (t) ) max X (th)

(130)

Xh (tf) e X max

(131)

ht ∈[t0,t]

(122)

Limits on the flow rates of the utilities: Fs e 10

(123)

Fw e 650

(124)

2.10.2. Mathematical Treatment. Consider a path constraint on a variable X of the form

X (t) e X

max

, ∀t ∈ [t0, tf]

(125)

A common way of dealing with such a constraint is to apply it only at various discrete time instants, that is

Figure 5 shows a graphical interpretation of Xh (t). At any time t > 0, Xh (t) can be mathematically expressed by the following conditions:

IF [(Xh < X ) AND (dX /dt > 0)] THEN dXh dX ) dt dt

(132)

dXh )0 dt

(133)

ELSE

768


Table 2. Summary of Distillation Model Constraints constraints

a

description

hd ) 0 ha ) 0

DAE system

h0 ) 0

initial conditions

ge e 0 gq e 0 gy e 0

end-point inequalities time-invariant inequalities pure integer constraints

expressions

number

4-88 93-97, 103-112 134 or 139, 135, 137 for each inequality 113-124 141-149 89-92 101-102 136, 138 for each inequality 113-124 131 or 140 for each inequality 113-124 98 1-3 99, 100

1655a

295

72 50 42

94 time-invariant.

The discontinuous eqs 132 and 133 are smoothly approximated by

dXh 1 dX ) [1 + tanh(106‚float)] 1 + tanh 106 dt 4 dt

[

(

dX dt (134)

)]

Annualized installed cost of a column shell (ref 22, p 574) assuming a 3-year payback period: Cshell ) 1.4 ×

( )(

)

3.18 M&S 100Dcol 1.066 × 101.9 3 280 12 × 2.54 (heightstack + 15)0.802 (142) N

float ) X - Xh

(135)

heightstack ) space‚[(

∑ k‚yrk) - 1]

(143)

k)1

float(t0) ) 0

(136)

If X corresponds to a differential state, then its derivative is explicitly available for use in eq 134. Otherwise, dX/dt can be very accurately approximated by dXˆ /dt, where

|

t)t0

(137)

)0

( ) (

1 M&S × 4.7 3 280 100Dcol 12 × 2.54

)

1.55

‚heightstack (144)

(138)

The advantages of the new path constraint formulation defined by eqs 131 and 134-138 are that the path constraints are rigorously enforced; experience shows that it is generally scaling-independent (an exception may be for very badly scaled systems, in which case the constants of 106 in eq 134 and 10-8 in eq 137 need to be adjusted); and it is applicable to any length of time horizon. Compared to the approach described by eqs 126-128, this comes at the expense of defining one extra state variable, Xˆ , and this is only the case if the path constraint is on a nondifferential state. Note that, for a path constraint of the form X g X min, the minimum value of X is tracked via

dX dXh 1 ) [1 - tanh(106‚float)] 1 - tanh 106 dt 4 dt

[

Annualized installed cost of a column’s trays and internals (ref 22, p 575): Ctrays ) 1.4 ×

dXˆ X ) Xˆ + 10-8 dt dXˆ dt

where a scale-up factor of 1.4 has been used in eq 142 after comparison with some industrial cost data.

(

)] dX dt (139)

(140)

2.11. Cost Functions

Total annualized capital costs: Ccap ) Cshell + Ctrays + Creb + Ccond

( ) ( )

(

)

104‚Sreb 1 M&S × 101.3 Creb ) 0.6 3 280 144 × 2.542

(

(141)

4

0.65

× 3.29 ×

)

1.35 (145)

0.65

10 ‚Scond 1 M&S × 101.3 × Ccond ) 0.6 3 280 144 × 2.542 3.29 × 1.0 (146) Annual operating costs: Cop ) Csteam + Cwater

(147)

Average annual cost of steam and cooling water for 8150 h of operation per year: Csteam ) MWw‚Mscs

instead of eq 134, with the end-point inequality

Xh (tf) g X min

Annualized installed costs of the reboiler and condenser using a scaling factor of 0.6 (ref 22, p 572):

60 × 8150 tf

(148)

60 × 8150 tf

(149)

Cwater ) MWw‚Mwcw

2.12. Summary of Model Statistics. The constraints and variables for the complete process and control model are classified in Tables 2 and 3, respectively. The tables also show the number of each constraint and variable type for the case study, where the

Ind. Eng. Chem. Res., Vol. 41, No. 4, 2002 769 Table 3. Summary of Distillation Model Variables symbol

description

xd(t)

differential states

u(t) z ν(t) θ d yp yc xa(t)

manipulated variables operating variables disturbances uncertain parameters design variables process binaries control binaries algebraic variables

variables

number

Mi,∪, i ) 1, ..., NC, ∪ ) ∀k, c, d U∪, ∪ ) ∀k, c, d Ms, Mw, Im,n, bm,n, ∀m,∀n X h and X ˆ for ge e 0 R, D, Fw, Fs, B none zi,f, i ) 1, ..., NC - 1, Tw,in F Dcol, Sreb, Scond, Km,n, τm,n, set1,n, ∀m, ∀n yfk, yrk, k ) 1, ..., N ycm,n, ∀m, ∀n all others

295

5 0a 2 1 58 60 25 1353b (134 time-invariant)

a In this study, the set-points are treated as design rather than operating variables. b Does not include the six specified variables T ; f l ; xi,N+1; i ) 1, ..., NC; and Ts,in. LN+1; hN+1

maximum number of trays N ) 30 and the number of components NC ) 2. 3. Simultaneous Design and Control In this section, the application of the simultaneous design and control framework of Mohideen et al.8 (see Appendix C) is described in detail for the distillation case study. 3.1. Step 1: Initialization. An initial set of two scenarios, θ1 ) 6 and θ2 ) 6.6, is chosen with weights wT ) [0.75, 0.25]. These correspond to the nominal and upper values, respectively, of the feed flow rate, F. 3.2. Step 2: Multiperiod Design and Control. For this case study, using the nomenclature in Tables 2 and 3, the multiperiod design and control problem shown in Step 2 of Appendix C has the form ns

(expected) cost ) min {Ccap + d,yp,yc

s.t.

wjCop,j} ∑ j)1

do for this case study), the master problem does not require any explicit dual information with respect to the DAE system and so no intermediate adjoint problem is required for its construction. The master problem also has a very simple form compared to those of algorithms where adjoint variables are required, as will be illustrated later for this case study. Furthermore, the MIDO approach is independent of the type of method used for solving the dynamic optimization primal problems. It should be noted, however, that, like other MIDO algorithms that have been proposed in the literature, this approach guarantees local but not global optimality. The algorithm was implemented with gPROMS/gOPT v1.727 for solving the dynamic optimization primal problems and GAMS v2.50/CPLEX28 for the MILP master problems. Issues relating to the formulations of the primal and master problems are discussed below. 3.2.1. Primal Problems. In applying Step 2 of the MIDO algorithm described in Appendix D, the pth primal problem corresponds to ns

hd(x3 jd(t),xjd(t),xja(t),uj(t),ν(t),θj,d,yp, t) ) 0

cost ) min{Ccap + d

(150) ha(xjd(t),xja(t),uj(t),ν(t),θj,d,yp,

t) ) 0

h0(x3 jd(0),xjd(0),xja(0),uj(0),ν(0),θj,d,yp,t)0) ) 0 ge(x3 jd(tf),xjd(tf),xja(tf),uj(tf),ν(tf),θj,d,yp,tf) e 0, j ) 1, ..., ns gq(d,yc) e 0 gy(yp,yc) e 0 Note that, in this study, a time horizon tf ) 1440 min (1 day) was chosen because it is long enough to allow the system to reach a “cyclic steady-state” and so be representative of a whole year’s dynamic operation, while not leading to unacceptably long computation times. Equation 150 corresponds to a MIDO problem, which, with ns ) 2, consists of 3212 DAEs, 590 differential states, 144 end-point inequalities, 50 time-invariant inequalities, 42 pure integer constraints, 58 timeinvariant search variables, and 85 binary variables. This formidably sized problem was solved using the MIDO algorithm of Bansal,3 the steps of which are outlined in Appendix D. The main advantage of this algorithm over other MIDO algorithms is that, even when the binary variables participate within the DAE system (as they

s.t.

wjCop,j} ∑ j)1

hd(x3 jd(t),xjd(t),xja(t),uj(t),ν(t),θj,d,ypp,t) ) 0 (151) ha(xjd(t),xja(t),uj(t),ν(t),θj,d,ypp,t)

)0

h0(x3 jd(0),xjd(0),xja(0),uj(0),ν(0),θj,d,ypp,t)0) ) 0 ge(x3 jd(tf),xjd(tf),xja(tf),uj(tf),ν(tf),θj,d,ypp,tf) e 0, j ) 1, ..., ns Note that, with the control structure fixed, the timeinvariant inequalities, gq e 0, corresponding to eq 98, become redundant, while the controllers’ gains and reset times corresponding to loops not in the structure under consideration are all set to zero. This reduces the number of time-invariant search variables from 58 to 18. The size of the DAE system in eq 151 was also reduced by removing superfluous equations corresponding to the “shadow” trays above the reflux tray of the particular primal problem. In addition, the number of end-point inequalities was reduced from 144 to 28 by observing that the bottlenecks for the minimum column diameter due to flooding and the fractional entrainment limit are at the bottom and the top of the column, respectively, for this particular system. Clearly, the pure integer constraints, gy e 0, are not required for the primal problems, either.

770


For Step 3 of the MIDO algorithm, the problem was re-solved at the optimal solution from Step 2, with all of the constraints and variables included. Since the control structure binaries, yc, do not appear in the DAE system, the following formulation was used: ns

(expected) cost ) min{Ccap + d,yp

wjCop,j} ∑ j)1

hd(x3 jd(t),xjd(t),xja(t),uj(t),ν(t),θj,d,yp,t) ) 0

s.t.

(152) ha(xjd(t),xja(t),uj(t),ν(t),θj,d,yp,t)

)0

h0(x3 jd(0),xjd(0),xja(0),uj(0),ν(0),θj,d,yp,t)0)

where wfk are the Lagrange multipliers for the constraints yfk - yfpk> ) 0 and k ) 1, ..., 30; wrk are the Lagrange multipliers for the constraints yrk - yrpk > ) U,p L,p and λm,n are the Lagrange 0 and k ) 1, ..., 30; and λm,n p multipliers for the control structure inequalities, Km,n U L p - Km,n‚ycm,n e 0 and Km,n‚ycm,n - Km,n e 0, respectively. The Lagrange function when not using this MIDO algorithm is made up of terms relating to all the DAEs and inequalities that contain binary variables, namely eqs 6-8, 13, 19, 75, 98, 119, and 143. Hence 30

)0

ge(x3 jd(tf),xja(tf),uj(tf),ν(tf),θj,d,yp,tf) e 0, j ) 1, ..., ns

min

yfk(1,...,30),yrk(1,...,30),ycm,n(m)n-1,...,5),η

p l,p hk+1 etc.) (Lk+1

(153)

η

∑ yfk ) 1

dt +

30

∫0 ∑ k)1 1440

-

dt

{ [∑ 30

p ν2,k (

yrk′)

∑

-

dt

k′)k

30

dt +

dUpk

k

p µ2,k (effpk

- 0.8 - 0.2

k)1

∑ yrk′) +

k′)1

30

30

k)1

k′)k

30

30

∑ k)1

N

yrk ) 1

k)1

p (Rpk - Rp‚yrk)] dt + ∫01440 ∑ [ν4,k

30

p λ1,k [(

∑ k′)k

30

yrk′)‚entpk - 0.1] +

30

N

∑

yrk′)

p dMi,k

k)1

k)1

yfk -

(

k′)k

p 2 (R‚velk-1 etc.)]} dt +

N

∑

{ [∑ 30

p ν1,i,k

p p [Pk-1 - Ppk - 10-5( ∑ yrk′) × ∫01440 ∑ {ν3,k

p p Lp(dp,yjfk(‚) ,yrk(‚) ,ycm,n(‚)) e η, p ∈ Pfeas

s.t.

]} ]}

p p xi,k+1 etc.) (Lk+1

and the Lagrange multipliers obtained for hy ) 0 and gq e 0. 3.2.2. Master Problem. The explicit form of the master problem (for feasible primal problems) is

30 2

∫0 ∑ ∑ k)1i)1 1440

hy ) yp - ypp ) 0 gq(d,ycp) e 0

p µ1,k (Fpk - Fp‚yfk) + ∑ k)1

Lp ) costp +

yrk′ e 0, k ) 1, ..., 30

2[(

∑ k)1

5

k‚yrk) - 1]} +

p p µ3,k {heightstack ∑ k)1

5

U,p p U [λm,n (Km,n - Km,n ‚ycm,n) + ∑ ∑ n)1m)1

k′)k L,p L p (Km,n ‚ycm,n - Km,n )] (155) λm,n

5

∑ ycm,n ) 1,

m ) 1, ..., 5

n)1 5

∑ ycm,n ) 1,

n ) 1, ..., 5

m)1

Earlier, we remarked that the master problem when using the MIDO algorithm of Bansal3 is of a much simpler form compared to that when using other MIDO algorithms such as those of Mohideen et al.29 and Schweiger and Floudas.30 This can be clearly seen for this example by comparing the Lagrange function, Lp, that results in both cases. For the MIDO algorithm used in this article, with the primal problem formulation eq 152, the Lagrange function in eq 153 is given by 30

Lp ) costp + 5

∑ k)1

30

wf pk(yf pk - yfk) +

wf pk(yrpk - yrk) + ∑ k)1

5

U,p p U L,p L (Km,n - Km,n ‚ycm,n) + λm,n (Km,n ‚ycm,n ∑ ∑ [λm,n n)1m)1 p )] (154) Km,n

where ν refers to time-dependent dual (adjoint) multipliers for the DAE system, µ refers to Lagrange multipliers for time-invariant equalities, and λ refers to Lagrange multipliers for the end-point and time-invariant inequalities. The Lagrange function (eq 154) consists of the primal objective value plus 85 simple linear terms with multipliers whose values are all obtained directly from the primal solution of eq 152. In sharp contrast, the Lagrange function (eq 155) consists of the primal objective plus 295 terms. One hundred and fifty of these terms occur within an integral, a large number of primal solution values need to be stored (e.g. the values of dMi,k/ dt, dUk/dt, etc.), and the values of the dual multipliers ν are difficult to obtain, since a complex, intermediate adjoint problem must be solved. Note that, when solving the master problem (eq 153), convergence of the overall problem was significantly aided by placing some restrictions on the values of the binary variables, based on engineering knowledge of the process. For example, the composition of the feed means that the feed tray cannot be located very near the top

Ind. Eng. Chem. Res., Vol. 41, No. 4, 2002 771 Table 4. Progress of Iterations for the Multiperiod MIDO Design and Control Problem iteration number primal solutions: discrete decisions no. of trays feed tray control schemea process design Dcol (m) Sreb (m2) Scond (m2) controllers’ gains K∪,1 (x1,d) K∪,2 (leveld) K3,3 (pc)b K∪,4 (x1,b) K∪,5 (level0) reset times τ∪,1 τ∪,2 τ∪,3 τ∪,4 τ∪,5 set-points set1,1 set1,2 set1,3 set1,4 set1,5 costs ($100k yr-1) capital operating(1)c operating(2)c expected UB master solutions: no. of trays feed tray control schemea LB UB - LB e 10-4?

1

2

3

4

Table 5. Dynamic Feasibility Test Results for the Optimal System Found in Table 4

5

25 15 2

24 14 1

23 13 1

22 13 1

24 13 1

2.03 127.6 91.45

1.99 134.2 85.03

1.99 140.0 84.13

2.00 138.9 84.02

2.00 138.0 85.78

6.70 -105.0 -28.00 9.71 -1042

33.74 -41.29 -31.44 -2.22 -600.0

48.85 -18.64 -29.24 -2.38 -560.1

70.00 -24.55 -26.57 -3.37 -580.5

32.10 -25.39 -36.16 -0.93 -550.0

160.0 530.0 9935 2325 663.6

87.3 568.2 3483 59.8 693.7

100.0 568.9 3615 66.3 662.1

143.2 568.9 5032 61.7 664.2

77.2 684.5 2809 150.6 695.2

0.9883 0.5368 1.1944 0.0182 0.6002

0.9849 0.0668 1.2800 0.0223 0.8995

0.9843 0.0773 1.3022 0.0250 0.8994

0.9835 0.0746 1.3164 0.0293 0.8980

0.9853 0.0703 1.2694 0.0179 0.8995

1.941 6.367 7.220 8.521 8.521

1.883 6.268 7.122 8.364 8.364

1.858 6.287 7.136 8.357 8.357

1.823 6.334 7.194 8.372 8.357

1.894 6.269 7.097 8.370 8.357

24 14 1 8.242 No

23 13 1 8.282 No

22 13 1 8.341 No

24 13 1 8.355 No

22 11 1 8.357 Yes: STOP

Control scheme 1: R-x1,d, D-leveld, Fw-Pc, Fs-x1,b, B-level0. Control Scheme 2: R-x1,d, D-leveld, Fw-Pc, B-x1,b, Fs-level0.b For K3,3, the cooling water flow rate is scaled (0.01Fw). c Period 1: nominal feed flow rate, F ) 6 kmol min-1. Period 2: feed flow rate at upper bound, F ) 6.6 kmol min-1. a

constraint

Xl

critical feed flow rate

eq 113: 0.98 - x1,d e 0 eq 114: x1,b - 0.05 e 0 eq 115: 1.013 25 - Pc e 0 eq 116a: 0.06 - leveld e 0 eq 116b: leveld - 0.54 e 0 eq 117a: 0.6 - level0 e 0 eq 117b: level0 - 0.9 e 0 eq 118: Dmin col,1 - Dcol e 0 eq 119: ent23 - 0.1 e 0 eq 120: T0 - Ts,out e 0 eq 121: Qreb/Sreb e 0 eq 122: Tw,out/1e4 - 0.032 315 e 0 eq 123: Fs - 10 e 0 eq 124: Fw/6.5e4 - 0.01 e 0

6 × 10-6 5 × 10-6 -0.24 2 × 10-5 -0.45 -0.30 1 × 10-6 4 × 10-5 4 × 10-5 -0.15 -17.0 -1 × 10-6 -2.2 -2 × 10-6

6.6 6.6 6.6 6.6 6.6 6.6 6.6 6.6 6.6 6 6.6 6 6.6 6.6

the optimal solution is actually found on the third iteration. 3.3. Step 3: Dynamic Feasibility Test. In this step, the optimal process and control design from Step 2 is tested for feasibility over the whole range of feed flow rates, 6-6.6 kmol min-1. Because there are no operating variables for this problem and all the path inequality constraints have been converted into end-point inequalities, the dynamic feasibility test (eq 173) becomes

χ(d,yp,yc) )

max

6eθe6.6,l∈L

ge,l(x3 d(tf),xd(tf),xa(tf),u(tf), ν(tf),θ,d,yp,tf)

hd(x3 d(t),xd(t),xa(t),u(t),ν(t),θ,d,yp,t) ) 0

s.t.

(156) ha(xd(t),xa(t),u(t),ν(t),θ,d,yp,t) ) 0 h0(x3 d(0),xd(0),xa(0),u(0),ν(0),θ,d,yp,t)0) ) 0 This can be solved in the following manner:32

Step 1 or bottom of the column if the products’ quality specifications are to be met. This was modeled by setting yfk ) 0 for k ) 1, ..., 6 and k ) 25, ..., 30. In terms of the control structure, it is well-established in industrial distillation control practice that in order to provide effective pressure control, for a liquid distillate, the pressure is almost always paired with the condensation rate, which is determined by the cooling water flow rate (ref 19, p 490). Hence, yc3,3 ) 1. In addition, on the basis of the fact that a very long time response is undesirable between the action of a manipulated variable and its effect on a controlled variable (ref 31, p 525), pairings such as R-x1,b, D-x1,b, Fs-x1,d, and R-level0 can be eliminated; that is, yc1,4 ) yc2,4 ) yc4,1 ) yc1,5 ) 0. 3.2.3. Results. With the restrictions described above, the problem is still of a highly combinatorial nature, with 1953 discrete alternatives. [There are 7 possible control structures; there are 18 possible feed tray locations from trays 7-24 inclusive; and for the feed located on tray k, there are (31 - k) alternatives for the reflux tray location. Hence, the total number of discrete 24 (31 - k) ) 1953.] Nevertheless, alternatives is 7∑k)7 with a termination tolerance ) 10-4 the MIDO algorithm converged in just five iterations. The primal and master design and control solutions at each iteration are shown in detail in Table 4. It can be seen that

For each constraint, l ∈ L, solve the dynamic optimization problem χ(d,yp,yc) ) max ge,l(‚) 6eθe6.6

(157)

hd(‚) ) 0 ha(‚) ) 0 h0(‚) ) 0 Step 2 Set χ(d,yp,yc) ) max {χl} l∈L

The results of applying the above strategy for the dynamic feasibility test are summarized in Table 5. It can be seen that there are no additional critical scenarios to the two feed flow rates already considered and that, within the accuracy used for the dynamic optimization (3 × 10-4), χ = 0. Hence, the simultaneous design and control algorithm terminates. 3.4. Remarks. 1. The third primal solution in Table 4 is an economically optimal process and control structure and design that is guaranteed to give feasible dynamic operation for the disturbances considered, over all feed flow rates in the range 6-6.6 kmol min-1.

772


Figure 6. Simultaneous design and control result: controlled distillate composition at F ) 6 kmol min-1.

2. The optimal control structure consists of the following pairings: (i) reflux flow rate, R, controlling the distillate composition, x1,d; (ii) distillate flow rate, D, controlling the liquid level in the reflux drum, leveld; (iii) cooling water flow rate, Fw, controlling the pressure of the condenser, Pc, and hence the pressure in the column; (iv) steam flow rate, Fs, controlling the bottoms composition, x1,b; and (v) the bottoms flow rate, B, controlling the liquid level in the reboiler, level0. This corresponds to one of the three or four principal control structures recommended in the practical distillation control literature (e.g. see ref 19). 3. Refer to Table 4. Although there were seven different control structures under consideration in the study, the MIDO algorithm immediately found the optimal structure with the first master problem. From iteration 2 onward, the algorithm kept this control structure but modified the number of trays and the feed location of the column. Given the flatness of the expected total costs of primal solutions 2-5 compared to those of primal solution 1, it is apparent that, for this example, there are a number of alternative distillation designs that give similarly good dynamic performance provided that the “correct” control structure is used. 4. Unlike the many integrated design approaches that rely on the use of steady-state or approximate, linear dynamic models, the simultaneous design and control framework used in this study does not require an a posteriori nonlinear dynamic simulation to be performed to verify the resultsthis is given as part of the solution of the mixed-integer dynamic optimization problem. Figures 6 and 7 show the controlled distillate and bottom product compositions, respectively, for the optimal system at the nominal feed flow rate of 6 kmol min-1, while those for the second scenario of 6.6 kmol min-1 are shown in Figures 8 and 9. Notice how one of the compositions, in this case the distillate, is tightly controlled relative to the other. This effect of controlling both compositions with one tight loop and one loose loop is due to the negative interaction of the two loops and is again a common feature of distillation control.33 It is also interesting to note that Figures 6-9 all show the dynamic operation over 1 week, which is seven times longer than the time horizon that was used for solving the primal dynamic optimization problems. Although no rigorous stability test has been included in this study,34 the cyclical behavior and the fact that none of the

Figure 7. Simultaneous design and control result: controlled bottoms composition at F ) 6 kmol min-1.

Figure 8. Simultaneous design and control result: controlled distillate composition at F ) 6.6 kmol min-1.

Figure 9. Simultaneous design and control result: controlled bottoms composition at F ) 6.6 kmol min-1.

constraints are violated over the 7 day period suggest that the optimal solution found is indeed stable. 5. Table 6 compares the process design resulting from the application of the simultaneous design and control framework (labeled “III”) with the process designs that would have resulted if design and control had been considered sequentially (labeled “I” and “II”; both were obtained through solution of mixed-integer, nonlinear programs (MINLPs) using the steady-state analogue of the mixed-integer dynamic distillation model presented

Ind. Eng. Chem. Res., Vol. 41, No. 4, 2002 773 Table 6. Comparison of Optimal Steady-State and Dynamic Process Designs design variable no. of trays feed location Dcol (m) Sreb (m2) Scond (m2) capital cost ($100k yr-1) expected cost

I: nominala II: operable atb III: operablec steady-state steady-state dynamically 23 12 1.82 112.9 83.47 1.687 7.596

23 12 1.91 131.6 83.33 1.785 7.908

23 13 1.99 140.0 84.13 1.858 8.357

a

Inoperable at steady-state. b Inoperable dynamically; that is, no PI-control scheme can be found for which all the constraints can be satisfied. c Dynamically operable design resulting from application of simultaneous design/control framework; see Table 4.

Table 7. Solution of Dynamic Infeasibility Minimization Problem (Eq 175) for the Flexible Steady-State Design II in Table 6 Process Design (Fixed) no. of trays feed tray Dcol (m) Sreb (m2) Scond (m2)

23 12 1.91 131.6 83.33 Control Designa

K1,1 K2,2 K3,3 K4,4 K5,5 τ1,1 τ2,2 τ3,3 τ4,4 τ6,5 set1,1 set2,2 set3,3 set4,4 set5,5

250 -81.72 -37.50 -63.25 -351.4 271.3 549.1 429.6 91.5 599.8 0.9810 0.0660 1.2500 0.0469 0.6013

Inequality Violations (κl) at F ) 6 kmol min-1 eq 119:b 0.1 ent23 e 0.01 9.78 × 10-4 eq 122:c 10-4Tw,out e 3.2315 × 10-2 3.10 × 10-4 Inequality Violations (κl) at F ) 6.6 kmol min-1 1.50 × 10-2 eq 118:d Dmin col,1 - Dcol e 0 eq 119:e 0.1 ent23 e 0.01 1.02 × 10-3 eq 124: Fw/(6.5 × 104) e 10-2 6.68 × 10-5 sum of violations (obj) 1.74 × 10-2

Figure 10. Solution of dynamic infeasibility minimization problem (175) for the flexible steady-state design: violation of the minimum column diameter constraint at F ) 6.6 kmol min-1.

in Section 2 of this article). Design I corresponds to the optimal design at nominal conditions. This system is inoperable at steady-state over the ranges of feed flow rate, cooling water temperature, and feed composition. Design II corresponds to the optimal “flexible” design. This was obtained by applying the design algorithm of Halemane and Grossmann35 where the expected cost was minimized over nine different scenarios (nominal point with a weight of 0.5 and eight vertex points, each with a weight of 0.125) with the reflux, cooling water, and steam flow rates used as operating variables. Design II is operable at steady-state but inoperable dynamically when disturbances are considered, as shown in Appendix E. By comparing design II with design III in Table 6, it can be seen that to accommodate feed flow rates, cooling water temperatures, and feed compositions that are away from nominal conditions requires far more overdesign compared to the nominal design I when the dynamic behavior of the system is accounted for than when only steady-state effects are considered. Design II has a smaller column diameter than design III, leading to flooding problems (Figure 10 in Appendix E) and excessive fractional entrainment (Table 7 in Appendix E), and a smaller condenser area, leading to overheating of the cooling water (Figure 11 in Appendix E). This illustrates a weakness of considering design and control in a sequential manner. 4. Conclusions In this article, a new, rigorous, mixed-integer dynamic model for a multicomponent distillation system has been

a Control structure fixed at that of Scheme 1: R-x , D-level , 1,d d Fw-Pc, Fs-x1,b, B-level0. b Corresponds to the fractional entrainc ment reaching 0.11. Corresponds to the outlet cooling water reaching 326.25 K. d Corresponds to the column diameter being 1.50 cm above the minimum allowable. e Corresponds to the fractional entrainment reaching 0.11.

Figure 11. Solution of dynamic infeasibility minimization problem (175) for the flexible steady-state design: violation of the cooling water outlet temperature constraint at F ) 6 kmol min-1.

developed that relaxes many of the assumptions used in distillation studies that have appeared in the literature. This model has been utilized in a study where the number of trays, feed location, column diameter, reboiler and condenser surface areas, 5 × 5 control structure, PI-gains, reset times, and set-points of a benzenetoluene system subject to disturbances and uncertainty were simultaneously determined to give an economically optimal system that is guaranteed to give feasible operation. This required the solution of a mixed-integer dynamic optimization problem involving approximately

774


3000 DAEs, 600 differential states, 150 end-point inequalities, and 85 binary variables. As such, we have demonstrated how the framework of Mohideen et al.8 can be used for the simultaneous design and control of processes described by realistic models, and illustrated that mixed-integer dynamic optimization problems involving thousands of differential-algebraic equations can now be solved using the latest algorithms. Acknowledgment The authors gratefully acknowledge financial support from the EPSRC, the Centre for Process Systems Engineering Industrial Consortium, and DETR/ETSU. V.B. would also like to thank ICI (Ph.D. Scholarship) and the Society of the Chemical Industry (SCI Messel Scholarship). An earlier version of the study presented in this article was presented at the PSE 2000 Conference, Keystone, CO, as briefly described in Bansal et al.36 Appendix A: Modeling Parameters and Variables A.1. Process Parameters. See Table 8. Table 8. Process Parameters symbol ai cjpw cs cw Dcond Ddrum Dpipe Dreb Floodfac g heightcondweir heightweir Kloss lengthcond lengthcondweir lengthdrum lengthreb M&S MWw N NC Nin Nmeas space Ucond Ucondsn reb Usub reb R β ∆hvapo (Ts,in) s

definition Murphree efficiency parameter for component i on tray k average molar heat capacity of water unit cost of steam unit cost of cooling water diameter of the condenser diameter of the reflux drum diameter of the condenser inlet pipe diameter of the reboiler flooding factor acceleration due to gravity condenser weir height tray weir height minor loss coefficient length of the condenser length of the condenser weir length of the reflux drum length of the reboiler Marshall & Swift index molar mass of water upper bound on number of trays no. of components no. of manipulated variables no. of measured variables tray spacing heat-transfer coefficient for the condenser heat-transfer coefficient for the condensing section of the reboiler heat-transfer coefficient for the subcooling section of the reboiler dry pressure drop coefficient aeration factor molar enthalpy of vaporization of steam at its inlet temperature

value

units

0.8a 0.0753

MJ kmol-1 K-1

0.01 10-5 0.5 0.6 0.1731

$ kg-1 $ kg-1 m m m

1.0 0.8 9.81 0.35 0.05b 0.04c 3 2 1.5 4 1100d 18.015 30

m m s-2 m m

pairing

index, m, n

lower bound on the controller L gain, Km,n

R-x1,d R-leveld R-Pc R-x1,b R-level0 D-x1,d D-leveld D-Pc D-x1,b D-level0 Fw-x1,d Fw-leveld Fw-Pc Fw-x1,b Fw-level0 Fs-x1,d Fs-leveld Fs-Pc Fs-x1,b Fs-level0 B-x1,d B-leveld B-Pc B-x1,b B-level0

1, 1 1, 2 1, 3 1, 4 1, 5 2, 1 2, 2 2, 3 2, 4 2, 5 3, 1 3, 2 3, 3 3, 4 3, 5 4, 1 4, 2 4, 3 4, 4 4, 5 5, 1 5, 2 5, 3 5, 4 5, 5

5 -100 1 0.4 5 -50 -150 -104 20 200 -1000 10 -100 -500 -1000 -100 5 1 -10 -1500 100 -104 -104 1 -650

kmol-1

Since benzene and toluene are very similar compounds and the distillation system operates at close to atmospheric pressure, their liquid and vapor thermophysical properties can be accurately described by ideal models. B.1. Process Parameters. See Table 10. B.2. Liquid NC

hl∪ )

l xi,∪hi,∪ ∑ i)1

0.002

MW m-2 K-1

0.001

MW m-2 K-1

1.0b 0.6b 38.2531f MJ kmol-1

a

For binary systems, ai,k is the same for both components. Theoretical and empirical methods are available for its calculation (see ref 21 for a review). Here an 80% efficiency is assumed for all the trays. b Source: ref 15. c Corresponds to rounded pipe entrances where curvature/pipe diameter > 0.15.18 d Source: Chemical Engineering. e Source: ref 22. f Corresponds to Ts,in ) 420 K. Source: ref 37.

A.2. Control Structure Parameters. See Table 9.

(158)

l is the molar enthalpy of pure component i and is hi,∪ given by

vapo , ∫TT cp,i dT) - ∆hi,∪

ft MW m-2 K-1

100 -1 100 40 500 50 -1 -10 104 105 -10 1000 -10 -10 -50 -5 100 100 0 -10 104 -100 -100 1000 -5

Appendix B: Physical Properties

form )( hi,∪ - href,i

2 5 5 2e 0.001

upper bound on the controller U gain, Km,n

U a The values of KL m,n and Km,n were found by spanning the inverse of the steady-state process gain between the mth manipulated variable and the nth measured variable.20 For K3,n, n ) 1, ..., 5, the cooling water flow rate is scaled (0.01Fw).

Molar enthalpy:

m m m m kg

Table 9. Control Structure Parametersa

∪

ref

i ) 1, ..., NC (159)

1 ) 0.001 cpA,i(T∪ - Tref) + cpB,i(T2∪ 2 1 T2ref) + cpC,i(T3∪ - T3ref) + 3

[

1 vapo c (T4 - T4ref) - ∆hi,∪ , i ) 1, ..., NC (160) 4 pD,i ∪

]

The enthalpy of vaporization of component i is determined from38 vapo ) ∆hi,∪

10

-6

{ [ hA,i

T∪ 1Tcrit,i

]

(hB,i+hC,i(T∪/Tcrit,i)+hD,i(T∪2/Tcrit,i2))

}

,

i ) 1, ..., NC (161)

Note that eq 49 requires the calculation of the molar enthalpy of pure water, hw. Equation 160 can be used


obtained from39

Table 10. Physical Properties value symbol

benzene (i ) 1)

toluene (i ) 2)

Constants for Calculating the Molar Heat Capacity of the Ideal Gas of Component i cpA,i -33.92 -24.35 cpB,i 0.4739 0.5125 -4 cpC,i -3.017 × 10 -2.765 × 10-4 cpD,i 7.130 × 10-8 4.911 × 10-8

form href,i

Standard Enthalpy of Formation of Component i at Tref ) 298.2 K 82.98 50.03

hA,i hB,i hC,i hD,i

Constants for Calculating the Molar Heat of Vaporization of Component i 0.475 × 108 0.50144 × 108 0.452 38 0.385 90 0.0534 0 -0.1181 0

MWi

Molar Mass of Component i 78.114 92.141

Pcrit,i

Parachor of Component i 205.1 245.1

R

Universal Gas Constant 8.314

Tcrit,i

Critical Temperature of Component i 562.16 591.80

vpA,i vpB,i vpC,i vpD,i

Constants for Calculating the Saturated Vapor Pressure of Component i -6.982 73 -7.286 07 1.332 13 1.380 91 -2.628 63 -2.834 33 -3.333 99 -2.791 68

FA,i FB,i FC,i FD,i

Constants for Calculating the Liquid Molar Density of Component i 1.0162 0.8488 0.265 50 0.265 55 562.16 591.80 0.282 12 0.287 80

1

0 Fi,∪ )

0 Fi,∪

(165)

FA,i

, i ) 1, ..., NC

(166)

NC

pari(Fl∪xi,∪ - Fv∪yi,∪)]4 ∑ i)1

(167)

B.3. Vapor

Molar enthalpy: NC

hv∪ )

v yi,∪hi,∪ ∑ i)1

(162)

v l vapo hi,∪ ) hi,∪ + ∆hi,∪ , i ) 1, ..., NC

(163)

where the saturated vapor pressure of component i is

(169)

Fugacity coefficient of component i: (170)

Molar density (perfect gas equation): Fv∪ )

The fugacity coefficient of component i

(168)

Molar enthalpy of pure component i:

v Φi,∪ ) 1, i ) 1, ..., NC

0.38

Psat,i , i ) 1, ..., NC P∪

i,∪

F FB,i[1+(1-(T∪/FC,i)) D,i]

σl∪ ) [0.001

where Tcrit,w ) 647.3 K and, for example, Tref,w ) 373.15 K with the corresponding ∆href,w ) 40.6581 MJ kmol-1.

l ) Φi,∪

NC x

∑ i)1

The surface tension using the Macleod-Sugden correlation39

for this purpose with constants cpA,w ) 32.24, cpB,w ) 1.924 × 10-3, cpC,w ) 1.055 × 10-5, and cpD,w ) -3.596 × 10-9. The molar enthalpy of vaporization of water at Ts,in or Ts,out can be estimated using the Watson correlation:39

( )

)

where the molar density of pure component i is calculated from38

Source for h∪,i and F∪,i: ref 38. All other data were obtained from ref 39.

vapo ∆hw,∪ ) ∆hvapo ref,w

)

Fl∪

a

T∪ 1Tcrit,w Tref,w 1Tcrit,w

[ ( ) ( ) ( ( )]

Tcrit,i T∪ Psat,i ) vpA,i 1 + Pcrit,i T∪ Tcrit,i T∪ 1.5 T∪ 3 + vpC,i 1 + vpB,i 1 Tcrit,i Tcrit,i T∪ 6 vpD,i 1 , i ) 1, ..., NC (164) Tcrit,i

The molar density

Critical Pressure of Component i 48.98 41.06

pari

( )

ln

100P∪ RT∪

(171)

Appendix C: The Framework of Mohideen et al. (1996)8 The steps of the iterative, decomposition algorithm proposed by Mohideen et al.8 can be summarized as follows (see Tables 2 and 3 for nomenclature): Step 1. Choose an initial set of points θi, i ) 1, ..., ns, representing alternative scenarios for the uncertainties. Step 2. For the current set of scenarios, determine the optimal set of process structure and design variables, control structure, and tuning parameters. This

776


is achieved by solving the multiperiod, mixed-integer dynamic optimization (MIDO) problem

ha(xid(t),xia(t),ui(t),zi,νi(t),θi,d,y,t) ) 0

Step 3. Re-solve the primal problem at the optimal solution found in Step 2 with the full set of search variables and constraints included, and additional constraints of the form y - y j ) 0, where y is now a set of continuous search variables and y j is the set of (complicating) binary variables, fixed at y j ) yp. Convergence will be achieved in one iteration. Obtain the multipliers ωp corresponding to these new constraints. Step 4. Construct the pth relaxed master problem from the p primal solution, Jp, and the Lagrange multipliers, ωp. This corresponds to the mixed-integer linear program (MILP)

h0(x3 id(t0),xid(t0),xia(t0),ui(t0),zi,νi(t0),θi,d,y,t0) ) 0

minη

ns

min

∑wiC(x3 id(tf),xid(tf),xia(t),ui(tf),zi,νi

{

d,y,z1,z2,...,zns i)1

(tf),θi,d,y,tf)} hd(x3 id(t),xid(t),xia(t),ui(t),zi,νi(t),θi,d,y,t) ) 0 (172)

s.t.

g(x3 id(t),xid(t),xia(t),ui(t),zi,νi(t),θi,d,y,t) e 0, i ) 1, ..., ns where wi are discrete probabilities for the selected ns wi ) 1. parameter points such that ∑i)1 Step 3. Test the design and control system from Step 2 for feasibility (and stability, if required: see Mohideen et al.34) over the whole range of the uncertain parameters, over the entire time horizon of interest. This corresponds to solving the dynamic feasibility test problem

χ(d,y) ) max min max gl(x3 d(t),xd(t),xa(t), θ

z

l∈L,t∈[0,tf]

u(t),z,ν(t),θ,d,y,t) s.t.

hd(x3 d(t),xd(t),xa(t),u(t),z,ν(t),θ,d,y,t) ) 0 (173) ha(xd(t),xa(t),u(t),z,ν(t),θ,d,y,t) ) 0

h0(x3 d(t0),xd(t0),xa(t0),u(t0),z,ν(t0),θ,d,y,t0) ) 0 If χ(d,y) e 0, then the system is feasible and the algorithm terminates; otherwise, the solution of eq 173 defines a critical uncertainty realization, θc, that is added to the current set of scenarios, before returning to Step 2. Note that if the active set formulation of Grossmann and Floudas32 is used, as proposed by Dimitriadis and Pistikopoulos40 and Mohideen et al.,8 then eq 173, like eq 172, corresponds to a MIDO problem. Appendix D: The MIDO Algorithm of Bansal (2000)3 The steps of the algorithm are briefly summarized below: Step 1. Set the termination tolerance, . Initialize the iteration counter, p ) 1; lower bound, LB ) -∞; and upper bound, UB ) ∞. Step 2. Fix the values of the binary variables, y ) yp, and solve a standard dynamic optimization problem (pth primal problem) to obtain a solution Jp. (For computational speed, it may be useful to omit the search variables and constraints that are superfluous due to the current choice of values for the binary variables.) Set UB ) min(UB,Jp) and store the values of the continuous and binary variables corresponding to the best primal solution so far. Note that if the primal problem is infeasible for the particular choice of binary variables, then an infeasibility minimization problem must be solved.

(174)

y j ,η

s.t.

Jp + (ωp)T(yp - y j ) e η, p ∈ Pfeas (ωp)T(yp - y j ) e 0, p ∈ Pinfeas

where Pfeas denotes the set of iteration numbers for which feasible primal problems are found and Pinfeas denotes the set for which infeasible primal problems are found; the pure binary constraints of the original system are added, as are integer cuts in order to exclude previous primal integer solutions. Solve eq 174 to obtain a solution ηp. Update the lower bound, LB ) ηp. Step 5. If UB - LB is less than or equal to a specified tolerance , or the master problem is infeasible, then stop. The optimal solution is given by UB and the values stored in Step 2. Otherwise, set p ) p + 1 with yp+1 equal to the integer solution of the pth master problem, and return to Step 2. Appendix E: Testing the Dynamic Operability of Design II The dynamic operability of the “flexible” process design II shown in Table 6 was tested considering the two critical scenarios and optimal PI-control structure found during application of the simultaneous design and control framework in Section 3. An infeasibility minimization problem was solved of the form (cf. eq 151 for nomenclature) 28

obj ) min dc

s.t.

κl ∑ l)1

hd(x3 jd(t),xjd(t),xja(t),uj(t),ν(t), θj,dc,dp,yp,t) ) 0 (175) ha(xjd(t),xja(t),uj(t),ν(t),θj,dc,dp,yp, t) ) 0

h0(x3 jd(0),xjd(0),xja(0),uj(0),ν(0),θj,dc,dp,yp,t)0) ) 0 ge,l(x3 jd(tf),xjd(tf),xja(tf),uj(tf),ν(tf),θj,dc,dp,yp,tf) e κl, l ) 1, ..., 28 j ) 1, ..., 2 κl g 0, l ) 1, ..., 28 For a fixed process design (dp,yp) and control structure, eq 175 corresponds to a dynamic optimization problem where the search variables are the gains, reset times, and set-points of the controllers (denoted by the vector dc). The solution of eq 175 is shown in Table 7. Since the objective function is greater than zero, this indicates


that no set of values for the tuning parameters can be found for which all the inequality constraints are satisfied. The principal bottlenecks for feasibility are the minimum column diameter requirement due to flooding (eq 118) at F ) 6.6 kmol min-1 (see Figure 10); the fractional entrainment limit (eq 119) at both F ) 6 and F ) 6.6 kmol min-1; and the limit on the outlet cooling water temperature (eq 122) at F ) 6 kmol min-1 (see Figure 11). Nomenclature Acol ) cross-sectional area of the column, m2 min Acol,k ) minimum allowable column cross-sectional area based on flooding calculations for tray k, m2 Aholes ) total area of all the active holes, m2 min Anet,k ) minimum net area for vapor-liquid disengagement based on flooding calculations for tray k, m2 Atray ) active area of the tray, m2 bmn ) bias of the mth manipulated variable B ) bottoms flow rate, kmol min-1 C∪ ) annual costs (total, capital, operating, column shell trays, reboiler, condenser, steam, water), $ yr-1 D ) distillate flow rate, kmol min-1 D∪ ) diameter (column, reflux drum), m min Dcol,k ) minimum allowable column diameter based on flooding calculations for tray k, m em,n ) error for the control loop between the mth manipulated variable and the nth measured variable effi,k ) Murphree efficiency for component i on tray k entk ) fractional entrainment for tray k F ) total feed flow rate, kmol min-1 Fk ) feed flow to tray k, kmol min-1 FLVk ) Sherwood flow parameter for tray k Fs,w ) flow rate of steam or cooling water, kmol min-1 hl,v ∪ ) molar liquid or vapor enthalpy (feed and all units), MJ kmol-1 hw ) molar enthalpy of water, MJ kmol-1 heightstack ) total height of the trays, ft i ) index set for the components Im,n ) integral error for the control loop between the mth manipulated variable and the nth measured variable k ) index set for the reboiler (k ) 0) and trays (k ) 1, ..., N) Km,n ) controller gain for the control loop between the mth manipulated variable and the nth measured variable K1k ) empirical flooding velocity coefficient for tray k L∪ ) liquid molar flow rate (trays and condenser), kmol min-1 L ˜ k ) liquid mass flow rate leaving tray k, kg min-1 lengthweir ) tray weir length, m level∪ ) liquid level (all units), m m ) index set for the manipulated variables Mi,∪ ) molar hold-up of component i (all units), kmol Ml,v ∪ ) molar liquid or vapor hold-up (all units), kmol Ms,w ) steam or water cumulative molar consumption, kmol measn ) nth measured variable n ) index set for the measured variables P∪ ) pressure (feed and all units), bar Qcond ) external cooling rate for the condenser, MW Qreb ) external heating rate for the reboiler, MW Qcondsn ) external heating rate for the condensing section reb of the reboiler, MW Qsub reb ) external heating rate for the subcooling section of the reboiler, MW R ) total reflux flow rate, kmol min-1

Rk ) reflux flow to tray k, kmol min-1 Scond ) surface area of heat exchange coil in condenser, m2 Sreb ) surface area of heat exchange coil in reboiler, m2 Scondsn ) surface area of condensing section in reboiler, m2 reb sub Sreb ) surface area of subcooling section in reboiler, m2 setm,n ) set-point for the control loop between the mth manipulated variable and the nth measured variable T∪ ) temperature (feed and all units), K T∪,in ) inlet temperature of steam or cooling water, K T∪,out ) outlet temperature of steam or cooling water, K um ) mth manipulated variable U∪ ) molar internal energy hold-up (all units), MJ kmol-1 Vk ) vapor molar flow rate from reboiler and trays, kmol min-1 V ˜ k ) vapor mass flow rate leaving tray k, kg min-1 velk ) velocity of vapor leaving tray k, m min-1 velflood ) flooding velocity of vapor leaving tray k, m min-1 k vol∪ ) volume (all units), m3 volld ) volume of liquid in the reflux drum, m3 xi,∪ ) liquid mole fraction of component i (all units) yi,∪ ) vapor mole fraction of component i (all units) / yi,k ) equilibrium vapor mole fraction of component i ycm,n ) binary variable for existence of pairing between mth manipulated variable and nth measured variable yfk ) binary variable for tray k receiving feed yrk ) binary variable for tray k receiving reflux zi,f ) mole fraction of component i in the feed -3 Fl,v ∪ ) liquid or vapor molar density (all units), kmol m l,v F˜ k ) liquid or vapor mass density (trays and reboiler), kg m-3 l σk ) liquid surface tension (trays and reboiler), dyn cm-1 l,v Φi,∪ ) liquid or vapor fugacity coefficent of component i (all units) τm,n ) reset time for the control loop between the mth manipulated variable and the nth measured variable, min

Literature Cited (1) Ziegler, J. G.; Nichols, N. B. Process Lags in Automatic Control Circuits. Trans. ASME 1943, 65, 433-444. (2) Pistikopoulos, E. N.; van Schijndel, J. M. G. Towards the integration of process design, process control & process operability - current status & future trends. In Foundations of ComputerAided Process Design, Snowmass, Colorado, July 18-23, 1999. (3) Bansal, V. Analysis, Design and Control Optimization of Process Systems Under Uncertainty. Ph.D. Thesis, Imperial College, University of London, 2000. (4) Morari, M. Effect of design on the controllability of chemical plants. In Interactions Between Process Design and Process Control; Perkins, J. D., Ed.; Pergamon Press: London, September 6-8, 1992; IFAC, pp 3-16. (5) Viswanathan, J.; Grossmann, I. E. A Combined Penalty Function and Outer-Approximation Method for MINLP Optimization. Comput. Chem. Eng. 1990, 14 (7), 769-780. (6) Walsh, S. P. K.; Perkins, J. D. Application of Integrated Process and Control-System Design to Wastewater Neutralization. Comput. Chem. Eng. 1994, 18 (SS), S183-S187. (7) Walsh, S. P. K.; Perkins, J. D. Operability and control in process synthesis and design. In Advances in Chemical Engineering; Anderson, J. L., Ed.; Academic Press: 1996; Volume 23, pp 301-402. (8) Mohideen, M. J.; Perkins, J. D.; Pistikopoulos, E. N. Optimal Design of Dynamic Systems Under Uncertainty. AIChE J. 1996, 42 (8), 2251-2272. (9) Bansal, V.; Ross, R.; Perkins, J. D.; Pistikopoulos, E. N. The Interactions of Design and Control: Double-Effect Distillation. J. Proc. Control. 2000, 10 (2-3), 219-227. (10) Ross, R.; Bansal, V.; Perkins, J. D.; Pistikopoulos, E. N.; Koot, G. L. M.; van Schijndel, J. M. G. Optimal Design and Control of an Industrial Distillation System. Comput. Chem. Eng. 1999, 23 (SS), S875-S878.

778


(11) Bansal, V.; Sakizlis, V.; Ross, R.; Perkins, J. D.; Pistikopoulos, E. N. New Algorithms for Mixed-Integer Dynamic Optimization. Submitted for publication to Comput. Chem. Eng. (12) Viswanathan, J.; Grossmann, I. E. Optimal Feed Locations and Number of Trays for Distillation Columns with Multiple Feeds. Ind. Eng. Chem. Res. 1993, 32 (11), 2942-2949. (13) Pantelides, C. C. Dynamic behaviour of process systems. MSc Lecture Notes, Imperial College, London, 1996. (14) Coulson, J. M.; Richardson, J. F. with Backhurst, J. R.; Harker, J. H. Chemical Engineering Vol. 1: Fluid Flow, Heat Transfer and Mass Transfer, 4th ed.; Pergamon: Oxford, 1990; pp 215-216. (15) McCabe, W. L.; Smith, J. C.; Harriott, P. Unit Operations of Chemical Engineering, 5th ed.; McGraw-Hill: New York, 1993; pp 327 and 561-565. (16) AspenTech. Aspen Plus Release 8 User’s Guide; 1988. (17) Lygeros, A. I.; Magoulas, K. G. Column Flooding and Entrainment. Hydroc. Proc. 1986, 65 (12), 43-44. (18) Fox, R. W.; Macdonald, A. T. Introduction to Fluid Mechanics, 4th ed.; Wiley: New York, 1992; p 365. (19) Kister, H. Z. Distillation Operation; McGraw-Hill: New York, 1990. (20) Narraway, L. T. Selection of Process Control Structure Based on Economics. Ph.D. Thesis, Imperial College, University of London, 1992. (21) Sinnott, R. K. Coulson and Richardson’s Chemical Engineering Vol. 6: Chemical Engineering Design, 2nd ed.; Pergamon: Oxford, 1993. (22) Douglas, J. M. Conceptual Design of Chemical Processes; McGraw-Hill: Singapore, 1988. (23) Sullivan, G. R.; Sargent, R. W. H. Development of Feed Changeover Policies for Refinery Distillation Units. Ind. Eng. Chem. Process Des. Dev. 1979, 18, 113-124. (24) Walsh, S. P. K. Integrated Design of Chemical Wastewater Treatment Systems. Ph.D. Thesis, Imperial College, University of London, 1993. (25) Swartz, C. L. E. Maximum of trajectory via hyperbolic tan. Personal communication, 1999. (26) Cobden, A. L. Integrated process and control system design for reactor systems. MPhil to Ph.D. Transfer Report, Imperial College, University of London, 1998. (27) Process Systems Enterprise Ltd. gPROMS Advanced User’s Guide. http://www.psenterprise.com, 1999. (28) Brooke, A.; Kendrick, D.; Meeraus, A. GAMS Release 2.25: A User’s Guide; The Scientific Press: San Francisco, 1992. (29) Mohideen, M. J.; Perkins, J. D.; Pistikopoulos, E. N. Towards an Efficient Numerical Procedure for Mixed Integer Optimal Control. Comput. Chem. Eng. 1997, 21 (SS), S457-S462. (30) Schweiger, C. A.; Floudas, C. A. Interaction of design and control: Optimization with dynamic models. In Optimal Control: Theory, Algorithms and Applications; Hager, W. W., Pardalos, P. M., Eds.; Kluwer Academic Publishers B.V.: Dordrecht, 1997; pp 388-435. (31) Stephanopoulos, G. Chemical Process Control; Prentice Hall International: Englewood Cliffs, NJ, 1984. (32) Grossmann, I. E.; Floudas, C. A. Active Constraint Strategy for Flexibility Analysis in Chemical Processes. Comput. Chem. Eng. 1987, 11 (6), 675-693. (33) Rademaker, O.; Rijnsdorp, J. E.; Maarleveld, A. Dynamics and Control of Continuous Distillation Units; Elsevier: Amsterdam, 1975.

(34) Mohideen, M. J.; Perkins, J. D.; Pistikopoulos, E. N. Robust Stability Considerations in Optimal Design of Dynamic Systems Under Uncertainty. J. Proc. Control. 1997, 7 (5), 371-385. (35) Halemane, K. P.; Grossmann, I. E. Optimal Process Design under Uncertainty. AIChE J. 1983, 29 (3), 425-433. (36) Bansal, V.; Perkins, J. D.; Pistikopoulos, E. N.; Ross, R.; van Schijndel, J. M. G. Simultaneous Design and Control Optimisation under Uncertainty. Comput. Chem. Eng. 2000, 24 (27), 261-266. (37) Perry, R. H.; Green, D. W.; Maloney, J. O. Perry’s Chemical Engineers’ Handbook, 6th ed.; McGraw-Hill: New York, 1984. (38) Daubert, T. E.; Danner, R. P. Data Compilation Tables of Properties of Pure Compounds; DIPPR, AIChE: New York, 1985. (39) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The Properties of Gases & Liquids, 4th ed.; McGraw-Hill: New York, 1987. (40) Dimitriadis, V. D.; Pistikopoulos, E. N. Flexibility Analysis of Dynamic Systems. Ind. Eng. Chem. Res. 1995, 34 (12), 44514462. (41) Lenhoff, A. M.; Morari, M. Design of Resilient Processing Plants-I. Process Design under Consideration of Dynamic Aspects. Chem. Eng. Sci. 1982, 37 (2), 245-258. (42) Palazoglu, A.; Arkun, Y. A Multiobjective Approach to Design Chemical Plants with Robust Dynamic Operability Characteristics. Comput. Chem. Eng. 1986, 10 (6), 567-575. (43) Palazoglu, A.; Arkun, Y. Design of Chemical Plants with Multiregime Capabilities and Robust Dynamic Operability Characteristics. Comput. Chem. Eng. 1987, 11 (3), 205-216. (44) Abbas, A.; Sawyer, P. E.; Yue, P. L. Integrated Design and Control of a Continuous Stirred Tank Reactor. Trans. IChemE. 1993, 71 (A), 453-456. (45) Luyben, M. L.; Floudas, C. A. Analyzing the Interaction of Design and Control-1. A Multiobjective Framework and Application to Binary Distillation Synthesis. Comput. Chem. Eng. 1994, 18 (10), 933-969. (46) Luyben, M. L.; Floudas, C. A. Analyzing the Interaction of Design and Control-2. Reactor-Separator-Recycle System. Comput. Chem. Eng. 1994, 18 (10), 971-993. (47) Brengel, D. D.; Seider, W. D. Coordinated Design and Control Optimization of Nonlinear Processes. Comput. Chem. Eng. 1992, 16 (9), 861-886. (48) Nishida, N.; Ichikawa, A. Synthesis of Optimal Dynamic Process Systems by a Gradient Method. Ind. Eng. Chem., Process Des. Dev. 1975, 14 (3), 236-242. (49) Nishida, N.; Liu, Y. A.; Ichikawa, A. Studies in Chemical Process Design and Synthesis. II. Optimal Synthesis of Dynamic Process Systems with Uncertainty. AIChE. J. 1976, 22 (3), 539549. (50) Rowe, D. A.; Perkins, J. D.; Walsh, S. P. K. Integrated Design of Responsive Chemical Manufacturing Facilities. Comput. Chem. Eng. 1997, 21 (SS), S101-S106. (51) Bahri, P. A.; Bandoni, J. A.; Romagnoli, J. A. Integrated Flexibility and Controllability Analysis in Design of Chemical Processes. AIChE. J. 1997, 43 (4), 997-1015.

Received for review February 16, 2001 Revised manuscript received August 4, 2001 Accepted October 30, 2001 IE010156N

A Case Study in Simultaneous Design and Control Using Rigorous

Recommend Documents