Nonlinear model predictive control of a packed distillation column

control was used to control the distillation column at an operating condition where the process gain ... linearize the model equations about the nomin...
0 downloads 0 Views 1MB Size
Ind. Eng. Chem. Res. 1993,32, 2345-2356

2345

Nonlinear Model Predictive Control of a Packed Distillation Column Ashutosh A. Patwardhant and Thomas F. Edgar' Department of Chemical Engineering, The University of Teras at Austin, Austin, Texas 78712

A rigorous dynamic model based on fundamental chemical engineering principles was formulated for a packed distillation column separating a mixture of cyclohexane and n-heptane. This model was simplified to a form suitable for use in on-line model predictive control calculations. A packed distillation column was operated at several operating conditions to estimate two unknown model parameters in the rigorous and simplified models. The actual column response to step changes in the feed rate, distillate rate, and reboiler duty agreed well with with dynamic model predictions. One unusual characteristic observed was that the packed column exhibited gain-sign changes, which are very difficult to treat using conventional linear feedback control. Nonlinear model predictive control was used to control the distillation column a t an operating condition where the process gain changed sign. An on-line, nonlinear model-based scheme was used to estimate unknownhimevarying model parameters. 1. Introduction

Distillation is one of the most frequently used unit operations in the chemical industry. Distillation columns are of two types-tray columns, in which vapor and liquid are contacted on discrete stages, and packed columns, in which contact occurs continuously through the column. There are several instances when the .use of packed columns is preferred to that of plate columns (Maddox, 1973). A large amount of work dealing with the modeling, simulation, and control of tray columns has been published. In contrast, relatively little control work with packed columns has been reported. To design packed columns with a good degree of confidence, to predict their behavior under various operating conditions accurately, and to design a control system which keeps the products at their set points in an efficient manner, both steady-state and dynamic models of good accuracy are necessary. Packed columns contain either random or structured packings. The past few years have seen an increse in the use of high-efficiency structured packings, especially in the retrofitting of existing plants. Structured packings offer increased interfacial area, while reducing form drag. Structured packings are available in a variety of surface textures-sheet metal with smooth, fluted, lanced, or perforated surfaces. Packings may also be made from gauze material. Knowledge of the response of packings with different surface treatment to different operating conditions is very important. A column containing efficient packing may exhibit behavior that requires special consideration when designing a control system. Chemical processes in general, and distillation columns in particular, are inherently nonlinear in nature. Several approaches have been used to design controllers for such nonlinear processes. The most common approach is to linearizethe model equations about the nominal operating point, and then apply one of the many methods available for designing linear control systems. This approach has proved t o be satisfactory in a large number of process control applications. There are, however, some critically important control cases where process operation is not restricted to mildly nonlinear regions. In such situations, linear control proves ineffective. Nonlinear control methodologies based on differential geometric concepts have received considerable attention from chemical engineers

* To whom correspondence should be addressed.

+ Present address: Shell DevelopmentCompany,Houston, TX.

0888-588519312632-2345$Q4.QQlQ

in recent years. These strategies involve the transformation of the nonlinear process model into an equivalent linear process by means of nonlinear variable transformations. The resulting transformed process is then controlled using a linear controller, and the manipulated variables for the actual process are calculated using the inverse transformations. However, it may not always be possible to use the variable transformation method, because suitable transformations may not always exist. Nonlinear transformations are tedious to develop and may be too sensitive to modeling error (Alsop and Edgar, 1990; Patwardhan et al., 1990). Constraints on the states and manipulated variables are also difficult to implement using this approach. An attractive alternative is model predictive control using nonlinear programming techniques. Model predictive control algorithms can be applied to processes described by linear as well as nonlinear model equations. These algorithms can deal with process constraints, such as bounds on state and manipulated variables,in an explicit manner. An approach that uses orthogonal collocation to convert dynamic model equations into algebraic constraints is called nonlinear model predictive control, or NMPC (Patwardhan et al., 1990). With the incorporation of nonlinear process models, which are accurate over a broad range of operating conditions, NMPC permits process operation over a wide region without requiring controller retuning, and also allows operation in a region that is economically attractive, such as a steady-state optimum where a standard linear controller will not be able to controlthe process satisfactorily. For these reasons, nonlinear model predictive control has become a popular research topic in chemical process control. This paper reports on the experimental verification of a packed distillation column model, and the application of the model with NMPC to control the column in an operating region where the process gain changes sign.

2. Literature Survey 2.1. Dynamic Modeling of Packed Columns. Hitch et al. (1986) and Krishnamurthy and Taylor (1985) have presented comprehensive reviews of packed column models. This section outlines the unsteady-state packed column models proposed in the literature. Gray and Prados (1963) and Sakata and Pradoa (1972) investigated the dynamic response of a packed absorber using the frequency response technique and pulse testing, 0 1993 American Chemical Society

2346 Ind. Eng. Chem. Res., Vol. 32, No. 10, 1993

respectively. The model predictions were compared to experimental measurements for the ail-water-carbon dioxide system in a 0.1524-m4.d. column packedto a depth of 1.56m with0.016-m ceramicRaschig rings. The authors concluded that the steady-state gains were not well described by the model, and that the phase shift was dependent mainly upon the residence time in the column; axial diffusion did not play a major role. Tommasi and Rice (1970) modeled the dynamic response of the liquid composition profile in the enriching section of a packed distillation column to step changes in the reflux ratio. The model parameters were not related to any physically meaningful constants, but were determined by fitting the solution to experimentally determined dynamic responses. The use of field tests in modeling the unsteady-state behavior of a large-scale absorber was demonstrated by McDaniel and Holland (1970). As in their steady-state models, the transfer section approach was used to describe heat and mass transfer. Vaporization efficiencies were obtained from operating data. Bradley and Andre (1972) studied the dynamic response of a packed tower (0.076-m i.d., 2.82 m long) in which carbon dioxide was absorbed from an air-carbon dioxide mixture into aqueous monoethanolamine. The differential equations were simplified by making suitable assumptions and solved analytically usingthe method of characteristics. The model predictions were correct to within the repeatability of the experiments. Some model parameters were determined from steadystate experiments. von Stockar and Wilke (1977) emphasized the importance of including heat effects. These investigators used the relaxation method in order to obtain the steady-state composition and temperature profiles within the column, thereby avoiding problems encountered with the shooting method. Both the material and energy holdups in the gas phase were neglected. A mixing-cell type of model was employed,and mass transfer was described using the filmfactor concept. Hitch et al. (1987) extended their earlier work (Hitch et al., 1986) to include the dynamic behavior of packed absorbers by introducing accumulation terms into the material and energy balances. The finite-difference method was used to discretize the partial differential equations in the time dimension. It was found that neglecting the gas-phase holdup can introduce significant error into the predicted dynamic response, especially at high pressures, where the density of the vapor phase becomes comparable to that of the liquid phase. Srivastava and Joseph (1984) used orthogonal collocation to discretize model partial differential equations in the spatial direction, and a semiimplicit Runge-Kutta method to integrate the resulting set of ordinary differential equations in time. Using orthogonalcollocationnot only reduced the dimensionality of the problem (compared to using finite differences),which resulted in a reduction in the CPU time for steady-state simulation,but also made the incorporation of boundary conditions easier. Their simulation results agreed well with those of previous researchers (Tommasi and Rice, 1970;von Rosenberg and hadi, 1980; Treybal, 1969). The dynamics of isothermal absorption accompanied by chemical reaction (Benfieldhot-carbonate process) was studied by Suenson et al. (1985). A dispersion-typemodel was formulated to account for axial mixing. The model partial differentialequationswere discretized in the spatial coordinateusing orthogonal collocation, and the resulting ordinary differential equations were integrated in time using a fourth-order Runge-Kutta method. The authors

found that in order to use orthogonal collocation, it was necessary to include some dispersion in the model; otherwise the discontinuities in the liquid and gas flow rate profiles caused by step changes in the inputs could not be adequately described. 2.2. Nonlinear Model Predictive Control. Several model predictive control schemes have been reported in the literature. Those employing linear models include dynamicmatrix control (Cutler and M e r , 1980),model algorithmic control (Richalet et al., 1978),and simplified model predictive control (Arulalanand Deshpande, 1987). Nonlinear model-based control algorithms can be applied to processes described by a wide variety of model equations, such as nonlinear ordinary differential equations, algebraic equations, partial differential equations, integro-differentialequations, and delay-differentialequations. The models may be based on fundamental material and energy balances, rate expressions, thermodynamic principles, and empirically-derived relationships. Such models are accurate over a broad range of operating conditions. Model predictive control strategies based on these nonlinear models allow tight control of the process, improved constraint handling, and greater robustness in terms of handling unusual dynamics and time delays. In model predictive control, the controlproblem is posed as a nonlinear programming problem (NLP): optimize some objective function of the inputs and outputs such that (1)the model equations are satisfied and (2) other constraints (if any) on the states and manipulated variables are met. The approaches to solving this NLP fall into two broad categories. The sequential solution and optimization strategy was originated over 20 years ago and employs separate algorithmsto solve the differential equations and to carry out the optimization. First, the manipulated variable profile is guessed, and the differential equations are solved numerically to obtain the controlled variable profile. The objective function is then calculated. The gradient of the objective function with respect to the manipulated variable can be found either by numerical perturbation or by solving sensitivity equations. The control profile is then updated using some optimization algorithm, and the process repeated until the optimal profiles are obtained. Various versions of this strategy have been reported by Hicks and Ray (19711, Hasdorff (19761, Sargent and Sullivan (19781, Asselmeyer (1985), Morshedi (1986),Economouetal.(19861,Jangetal. (1987), Kiparissides and Georgiou (19871, and Peterson et al. (1989). The availability of accurate and efficient integration and optimization packages permits implementation of this method with little programming effort. However, it is difficult to incorporate constraints on state variables into this procedure. Also, the algorithm requires the solution of differential equations at each iteration. Jones and Finch (1984) found that such methods spend about 86% of the time integrating the model equations in order to obtain gradient information. This can make them prohibitive in terms of computation time and thus unattractive for use in real-time applications. The simultaneous solution and optimization strategy is a more recent and attractive alternative. This method was used for feedback control of nonlinear processes by Patwardhan et al. (1988,1989,1990),Eaton et al. (1990), and Eaton and Rawlings (1990). Examples of model predictive control based on nonlinear process models may also be found in Brengel and Seider (19891, Li and Biegler (1988, 1989, 19901, and Li et al. (Li and Biegler, 1990).

Ind. Eng. Chem. Res., Vol. 32, No. 10, 1993 2347 CwheadCmdenser

The reboiler is of the single shell, two tube pass BKU type. The tubes have an exterior surface area of 18.58m2, and the design duty is 6.95 X los W. Steam at up to 8.26 X lo5Pa (absolute) can be supplied to the tubes. A single shell pass, two tube pass BEM type heat exchanger is used as the overhead condenser. The area available for heat transfer is 14.31 m2,and the design duty is 6.95 X 106 W. Cooling water temperatures as low as 280 K can be achieved. Liquid stream samples can be withdrawn from the column a t several locations, and these are analyzed using a Gow-Mac gas chromatograph and Shimadzu integrator. The distillation test system is monitored and controlled by a Fischer-Porter DCI-5000 distributed control system. Input signals include signals from pressure transmitters, level transmitters, flow transmitters, and thermocouples. Output signals are mainly controlvalve signals,particularly for column pressure control, cooling water flow control, process stream flow control, and steam flow control. The 16-bit word length of a Micro-PDP-11/73 computer purchased in 1986imposes a 32K-word limit on the virtual address range of a task. Larger tasks must be built using overlays, so that, at any given time, the memory requirement does not exceed 32K. As a result of this limitation, and in order to use existing nonlinear programming software without extensive modification, a steady-state process model was used for state/parameter estimation and process control. It should be mentioned that many of these modeling limitations would not exist for currently available distributed control systems.

6@--i Reflux Drum

Reflux Heater

1

&

-w

Reflux Pump

F e d Pump

-------------

I

1

Bottom Product Bottoms Pump

Figure 1. Flowsheet of the distillation test system.

Biegler and Rawlings (1991) recently reviewed optimization-based approaches to nonlinear model predictive control. The most serious obstacle to using NMPC is that nonlinear dynamic process models are usually not available in the process industries, since their development and maintenance is both difficult and expensive. As the availability and fidelity of dynamic modeling platforms increase, it will become easier to obtain dynamic process models for use in design, control, and optimization. However, steady-state process models are easier to develop and are frequently available from the design stage, and it is sometimes possible to utilize these for process control. NMPC using a steady-state process model is described in this paper. 3. Distillation Test System

The distillation test system is located at the Separations Research Program's facility a t The University of Texas at Austin and is mainlyused to test and evaluate mass transfer and hydraulic characteristics of various column internals. The system can be operated at total reflux, as well as at finite reflux with feed and product streams. The flowsheet of the distillation test system in the material-balance control configuration is shown in Figure 1. The heart of the system consists of a distillation column made from an 0.457-m nominal diameter, schedule 40 carbon steel pipe. The inside diameter is 0.429 m. In the present study, the column was packed with Sulzer BX (Sulzer Brothers, Ltd.), a structured packing made of stainless steel gauze. The rectifying and stripping sections each contained 2.74 m of packing. In addition to nozzles for connections to the condenser, reboiler, and feed and reflux preheaters, there are ports at various locations for temperature measurement and liquid sample withdrawal. The operating pressure can be varied from full vacuum to 4.46 X lo5Pa (absolute) and the temperature from ambient up to 803 K.

4. Dynamic Modeling of the Distillation Test System 4.1. Rigorous Model. A dynamic model was formulated for the distillation test system, based on the following assumptions: 1. The chemical system of interest, cyclohexane + n-heptane, is an ideal system at low and moderate pressures, and heat-of-mixing effects are negligible. The largest parametric uncertainty in this model was introduced via the liquid and vapor film heat- and mass-transfer coefficients (Patwardhan, 1991). The liquid and vapor properties (density, heat capacity, and latent heat of vaporization) were assumed to be independent of composition and temperature, which did not introduce large additional error. The liquid density, heat of vaporization, and liquid and vapor heat capacities were obtained from Daubert and Danner (19851, and the Antoine coefficients were from Gmehling et al. (1980). 2. Equimolar counterdiffusion was assumed between the phases. This implied that the liquid and vapor flow rates were constant within each packed section, which decoupled the material and energy balances in the packed sections. The energy balance equations within the packing were thus omitted from the model. 3. The liquid and vapor film mass-transfer coefficients were combined into an overall vapor-phase mass-transfer coefficient. It is also possible to use the overall liquidphase mass-transfer coefficient instead. However, in distillation applications, the controlling resistance to mass transfer is usually in the vapor film, and the use of the overall vapor-phase mass-transfer coefficient is more common. The partial differential-algebraic equations formulated using the overall vapor-phase mass-transfer coefficient are listed below. This rigorous dynamic model will be denoted by RM. The term "rigorous" is used relative to the simplified model (SM), developed in section 4.2.

2348 Ind. Eng. Chem. Res., Vol. 32, No. 10, 1993

1.liquid-phase material balance for the light component:

2. vapor-phase material balance for the light component:

Apv4v

=

I62- Aqa,(y

- Y*)

(2)

3. vapor-liquid equilibrium at the interface:

y*U- x ) (3) x ( l - Y*) The unsteady-state equations for the auxiliary equipment and the boundary conditions for the packing partial differential equations are listed below. Perfect level control was assumed in the overhead accumulator and sump. 1. material balance at the feed location: a=

FL+ LBR= LTS Fv + VTS = VBR F r f L + LBRXBR LTSXTS

(4) (6)

~XTR dt = V T R ~ T-RXD) VTR = D + L T R

3. distillate rate specification:

D = distillate rate 4. material balance around the sump: dxB

Hsurnpp~ dt = LBS(XBS - XB) 5. energy balance around the sump: HsumppLcPL

+

(5)

Fvfv + VTSYTS= VBRYBR (7) 2. material balance around the overhead accumulator: %-PL

number may be used if necessary. The collocation points were chosen to be the zeros of Jacobi polynomials on the interval [0,11 in addition to the two end points. WithNc collocation points on each element (including the two end points), the dynamic model contained 6Nc 8 ordinary differential and algebraic equations in 6Nc + 8 dependent variables. These were integrated in time using DASSL (Petzold, 1982), a differential-algebraic system solver. 4.2. Simplified Model. The rigorous model formulated in section 4.1 was not suitable for use in on-line model predictive control calculations due to the memory limitation imposed by the Micro-PDP-11. This section describes the simplification of RM to a form suitable for use in on-line calculations. 1. RM used the reboiler heat duty as a manipulated variable. Specifying the vapor boilup rate instead of the reboiler heat duty decoupled the material and energy balances for the sump and reboiler. The energy balances around the sump and the reboiler were then omitted from the set of model equations, and

dt= LBSCPL(TL, - TB) dTB

was used as the vapor-liquid equilibrium relationship for the reboiler. 2. In order to conform to the memory limitation, it was necessary to use a steady-state process model in on-line computations. The time derivatives in all the RM equations were set to zero. This reduced the partial differential equations to ordinary differential equations, for which an analytical solution is available (Patwardhan, 1991). 3. When the material balance control configuration is employed, the product compositions are not sensitive to the thermal condition of the feed (Bequette, 1986). Thus, regardless of its actual condition, the feed was assumed to be completely vaporized, which implied xm = XBR. This simplified model (SM) was found to be suitable for use in on-line controller calculations. The simplified model equations are listed below: 1. mass transfer in the rectifying section:

6. material balance around the total reboiler: dXreb

HregL

dt= VB&B

- YBS)

7. energy balance around the total reboiler:

(18) 2. mass transfer in the stripping section:

dTVm

HreeLCP,

dt= VBSCP~(TB- T r e J +

8. vapor-liquid equilibrium in the reboiler:

The model equations for the auxiliary equipment serve as boundary conditions for the packing equations. For the sake of numerical implementation, the distillation column was divided into two finite elements, one finite element corresponding to the rectifying section and the other corresponding to the stripping section. The model partial differential equations were discretized in space on each of the finite elements using orthogonal collocation (Villadsen and Stewart, 1967; Holland and Liapis, 1983). In this study, an equal number of collocation points was used on each finite element, but an unequal

3. material balances:

Here,

Ind. Eng. Chem. Res., Vol. 32, No. 10, 1993 2349 0.994

of allowable values. This is especially true when the controlled and manipulated variables have "leveled out." In practical situations,therefore, one must estimate either p or n on-line, while the other is fixed at some nominal value. It should be mentioned that extrema of the type observed in Figure 2 are not generally observed in columns containing trays or random packing, but seem to be restricted to columns containing certain types of structured packing (Bravo et al., 1992).

0.990

0.986

5. On-Line State and Parameter Estimation

0.982

0.978 0

4 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0 . 0 Reboiler duty (million W)

0.018 1

I

For effective control, it is necessary to know the current values of the controlled variables and to have sufficiently accurate estimates of the model parameters. In distillation, it is difficult to obtain quick measurements of the product compositions. This makes it necessary to estimate product compositions and model parameters from the available process measurements-flow rates, temperatures, heat duty, and pressure. The use of nonlinear programming techniques for state and parameter estimation, or data reconciliation, has been previously investigated (Liebman, 1991; Wright and Edgar, 1991). State and parameter estimation was performed by solving the following nonlinear programming problem (NLP) at each sampling instant: NS

min

W,(S,

- in)'

n=1

0.04 0.06 0.08 0.10 0 . 1 2 0.14 0.16 0.18 0.20 Reboiler d u t y (million W)

such that eqs 18 and 19 are satisfied, and the estimated values of the distillate and boilup rates are accurate:

b=d

Figure 2. Predicted %D and %B versus heat duty for F = 0.044 kg/s. The constants in eq 24 were n = 0.01 and p = 50.

+

+

and X2R are the roots of x2 ~ R X ER = 0. CIS, (22.9, xls, and xsare defiied in a similar fashion. The volumetric X1R

mass-transfer coefficient was expressed in terms of the vapor flow rate: q 2 m = p ( 10"VJA)

(24)

K$zp was calculated in the rectifying and stripping sections by using the appropriate value of the vapor flow rate in eq 24. For a given set of inputs to the distillation column, the location of the extremum in a steady-state plot of bottoms composition versus the reboiler duty (Figure 2) was governed predominantly by p and n. These parameters were estimated from steady-state operating data. In order to estimate both p and n with reasonable confidence, the vapor rate should be varied over the widest possible range. The upper limit of this range is dictated by the flooding point, while heat losses to the surroundings determine the lower limit. The reboiler duty should be high enough that vapor reaches the top of the column despite condensation caused by heat loss. During normal column operation, the reboiler duty does not usually cover the entire range

(26)

vs = vs The NLP was solved using GRG2, a software package which employs the generalized reduced gradient method (Lasdon and Warren, 1986). During the solution of this NLP,D, Vs,andFwere held fixed at their measured values. On the basis of sensitivity analysis, n was fixed at a nominal value, while p was allowed to vary. In the objective function, i n is the measured composition and Sn is the estimated composition at the location denoted by the subscript n. 2n was assumed to be the bubble point composition corresponding to the temperature at the location denoted by the subscript n. In binary distillation, the bubble-point composition can be readily calculated if the temperature and pressure are known. In the present study, temperatures at four locations within the column were used to estimate the product compositions, i.e., Ns = 4. The objective function weights, Wn, were adjusted periodically (when composition measurements were available) to obtain good agreement between the estimated and measured compositions. The solution to the NLP provided estimates of the current distillate and bottoms compositions and two model parameters-feed compositions cf, and one of the parameters @) in eq 24. These estimates were filtered using a first-order lag before being used by the controller. The estimation program was executed once every 60 s, and required between 15 and 30 s for execution. Composition analyzer readings were available every 60 min. 6. NMPC

This section describes the use of NMPC to control the packed distillation column by inverting the steady-state process model described in section 4.2. The control

2350 Ind. Eng. Chem. Res., Vol. 32, No. 10, 1993 x.CLD

D ,Q R

, *"et B

* PROCESS

COSTROLLER

-

zD,XW

c

Figure 3. NMPC block diagram.

objective was to maintain the distillate and bottoms compositions at their set pointa by manipulating the distillate rate and the reboiler duty (material-balance scheme). NMPC may be represented using the IMC structure (Garcfa and Morari, 1982) as shown in Figure 3. At every sampling instant, the additive disturbance (dD, dg) is computed as the difference between the measured (XD, X B ) and predicted output ( 3 ~ZB) , values. As in IMC, a control move is calculated which takes the output to the set point biased by the additive disturbance. The control move is used to predict the output values at the next sampling instant. The control move and output prediction is part of a NLP formulated as shown below, and solved at every samplinginstant. For convenience, the subscript k denoting the current sampling instant has been omitted. It is understood that, unless explicitly denoted, all variable values are for the current sampling instant. Values from the previous sampling instant are clearly identified with the subscript k - 1,

such that eqs 18 and 19 are satisfied, and

D,, I D I D , ,

v,,

I

(29)

vsIv,

m u

(30)

The additive disturbances, dD and dg, were estimated as in DMC (Cutler and Ramaker, 1980):

dD = ZD - XD, dg = Zg - Xg,, When NMPC is implemented using a steady-stateprocess model, the additive disturbances, dD and dg, assume large values following step changes in the manipulated variables, and then decay to smaller values representative of steadystate model error. This is true even when the steady-state model is "perfect." Using dD and d g as feedback to the controller can therefore lead to poor controller performance. PerforFance may be improved by using filtered values (a, and dg) as feedback to the controller. A firstorder filter was used in the present investigation. Implementing manipulated variable constraints via explicit bounds in the NLP is not a good strategy when a feasible-path algorithm is used to perform the nonlinear optimization. Explicit bounds on the change in the manipulated variables can prevent the algorithm from achieving feasibility, thus causing it to terminate at an infeasible point. This is an advantage of using an infeasible-path algorithm for solving the NLP, because the constraints have to be satisfied only at the final iteration. Since GRG2 (a feasible-path algorithm) was used to solve the NLP, bounds on AD and AVS were

implemented by penalizing these changes in the objective function. F,f , p, and n were held constant while solving the control NLP. The solution provided new values for the manipulated variables, which became the set points for lower-level controllers; Le., D became the set point for the distillate flow controller and Qmb became the set point for the reboiler duty controller. The reboiler duty set point was obtained by multiplying the vapor boilup set point by the latent heat of vaporization. The solution also provided predictions of the product compositions for use in estimating the additive disturbances, dD and dg, at the next sampling instant. The NMPC algorithm was executed once every 3 min. 7. Open-Loop Tests

The dynamic response of the distillation column to changes in the feed rate, distillate rate, and reboiler duty was determined experimentally by conducting step testa. These step testa served to validate the dynamic process model, and totest the state/parameter estimator described in section 5. During a test, a step change was made in one of the inputs, while the other two inputs were held constant. Step increases and step decreases were made in both the feed rate and the distillate rate, while the heat duty was changed to cover the entire range of operating values. The upper bound on the heat duty was governed by the flooding point, while the lower bound depended upon heat losses to the surroundings. The reboiler duty was always kept large enough that vapor reached the top of the column despite condensation caused by heat loss. While the effect of feed composition changes was not studied explicitly, f changed during the course of each step test. Since the distillate and bottom products were returned to the feed tank forrecirculation, the totalamount of cyclohexane and n-heptane and hence the average composition in the system remained constant. Because the sump liquid holdup was comparable to that in the feed tank, a change in the bottoms composition caused the feed composition to change in the opposite direction. The influence of the distillate composition on the fed composition, while similar, was not large in magnitude, because the holdup in the reflux drum was much smaller than the sump holdup. Therefore, any change in operating conditions caused the feed composition to change in a manner that counteracted the change in bottoms composition. For example, an increase in feed rate increased both X D and XB. Since the total sump holdup was constant, an increase in X B led to an increase in the amount of cyclohexane contained in the sump. Since the total amount of cyclohexane in the system was fixed, there was less cyclohexane and more n-heptane in the feed tank; i.e., f decreased, which counteracted the increase in X D and X B . Therefore, with the present configuration (distillate and bottoms products returned to the feed tank, comparable sump and feed tank holdups, and small reflux drum holdup), changes in the bottoms compositionswere smaller than if the feed composition were kept constant. 7.1. Step Changes in the Reboiler Duty. Steadystate data were used to estimate two model constanta-p and n. The distillation column was operated at different heat duties for fixed feed and distillate flow rates. Steadystate compositions were measured at several locations in the column for each value of the reboiler duty. This was repeated at other values of feed and distillate flow rates, and the two constants p and n in eq 24 were determined for each (F, D ) pair by minimizing the least-squares objective function

Ind. Eng. Chem. Res., Vol. 32,No. 10,1993 2351 0.16F=0.108 kg/s, D=0.063 kg/s F=0.055 kg/s, D=0.032 kg/s

----

0.12-

m

0.08

-

0.04

-

n x

+

0.007

0

Figure 4. Bottoms composition versus reboiler heat duty for f = 0.4565. The constants in eq 24 were p = 7.43 and n = 0.0265. 0.020

F=0.108 kg/s, Dz0.063 kg/s F=0.055 kg/s, D=0.032 kg/s 0.015

$

-

5000

10000 15000 Time ( s )

20000

25000

Figure 6. Distillate composition response when the feed rate was increased from 0.056 to 0.066 kg/s (t = 4079 8 ) and then decreased to 0.044 kg/s (t = 14 700 8 ) . The distillate rate was constant at 0.032 kg/s and the reboiler duty at 1.47 X 10'5 W. Estimator Model RM Measured

-

-

0.010-

0.005-

-.""" . 0.04

0.08 0.12 0.16 Reboiler duty (million W )

0.00

Figure 5. Bottoms composition versus reboiler heat duty for f = 0.4565. The constants in eq 24 were p = 34.9 and n = 0.0117. Table I. Experimental Conrtants in Eq 24 for Various F and D Values

0.108 0.055 0.028

0.063 0.032 0.016

7.43 34.9 27.3

'.-----_______

0.20

0.0265 0.0117 0.0178

Note that the error between the measured and the predicted cyclohexane compositions was minimized, and not the relative error. This was because n-heptane was considered to be the bottom product, and the mole fraction of n-heptane in the bottoms is comparable to the mole fraction of cyclohexane in the distillate. As a result, the fitting error in XD is comparable to that for XB. The values of p and n at different operating conditions are shown in Table I. These constants vary with feed and distillate rates, because eq 24 is a simplification; in reality, %urn depends not only upon the vapor loading, but also upon the liquid loading and the vapor and liquid distribution. A change in loading alters the liquid distribution in the packing. After each reboiler duty step test, measured feed and product compositions revealed variations in the feed rate. This variation made it difficult to compare dynamic model predictions with the estimated response, because the change in XD and XB caused by variations in F was larger than the change caused by step changes in Qreb. Figures 4 and 5 show the predicted bottom product

0

5000

10000 15000 Time 1s)

20000

25000

Figure 7. Bottoms composition response when the feed rate was increased from 0.055 to 0.066 kg/s (t = 4079 8 ) and then decreased to 0.044 kg/s (t = 14 700e). The distillate rate was constant at 0.032 kg/s and the reboiler duty at 1.47 X 10'5 W.

composition as a function of the reboiler duty predicted by the model for two @, n) pairs (see eq 24) at different feed flow rates. For the lowp-high n combination (Figure 4), the bottoms composition decreased as the heat duty increased for both values of the feed rate; i.e., the process gain was always negative. In Figure 5 (high p, low n),on the other hand, the sign of the process gain depends upon the feed rate and the heat duty. This shows that the process gain changes sign for large p-small n combinations, which emphasizes the need for accurate knowledge of p and n. Also,since these parameters changed with varying feed rate, it became necessary to estimate them on-line in order to have values that closely reflected conditions existing at any given time. 7.2. Step Changes in the Feed Rate. The response to step changes in the feed rate was determined by making both positive and negative changes in F. Figures 6 and 7 show the distillate and bottoms composition response respectively, when the feed rate was increased from 0.055 to 0.066 kg/s and then reduced to 0.044 kg/s. In a separate run, the feed rate was f i t decreased and then increased. Responses for this run are shown in Figures 8 and 9.The compositionmeasurementsdenotsd by the diamonds were available at infrequent intervals via gas chromatographic analysis. In all cases, RM predicted the composition to within 0.02 mole fraction of the measured value, and the steady-state composition estimates also agreed well with the measurements. The process gain predicted by the

2362 Ind. Eng. Chem. Res., Vol. 32, No. 10, 1993 1.00

Zstimator Model RM Measured

Estimator Model RM Measured

8-

- - -- - - -

-’

0.90

a 0.80

0.70

0.70 5000

10000 15000 20000 25000 3C Time I s )

Figure 8. Distillate composition response when the feed rate was decreased from 0.055 to 0.044 kg/s (t = 5220 8) and then increased to 0.066 kg/s (t = 16 200 8). The distillate rate was constant at 0.032 kg/s and the reboiler duty at 1.47 X lob W. 0.04

Zstimator Model RM Measured

0

IO

10000 15000 Time (si

20000

25000

Figure 10. Distillate composition response when the feed rate was decreased from 0.032 to 0.025 kg/s (t = 4979 s) and then increased to 0.038 kg/s (t = 15 120 8). The feed rate was constant at 0.055 kg/s and the reboiler duty at 1.47 X l@W.

-

Estimator Model RM

0.03

-

----

0.03

*

2

5000

*

““?;,-.;I

c.02

II I

0.01 I

0.00 0

5000

10000 15000 20000 25000 Time (si

30000

Figure 9. Distillate composition response when the feed rate was decreased from 0.055 to 0.044 kg/s (t = 5220 s) and then increased to 0.066 kg/s (t = 16 200 8). The distillate rate was constant at 0.032 kg/s and the reboiler duty at 1.47 X l@W.

model also agreed with the experimental value. The dynamic response prediction was in good agreement with the estimated composition. The model appears to predict the distillate composition more accurately than the bottoms composition. This is because the composition range in Figures 6 and 8 is much larger than the composition range in Figures 7 and 9. The prediction error for XD is comparable to that for XB due to the reason given in section 7.1. 7.3. Step Changes in the Distillate Rate. Figures 10 and 11 show the distillate and bottoms composition response respectively, when the distillate rate was decreased from 0.032 to 0.025 kgls and then increased to 0.038 kg/s. In a separate run, the distillate rate was first increased and then reduced. Responses for this run are shown in Figures 12 and 13. As with step changes in F, RM predicted the composition to within 0.02 mole fraction of the experimental value. The predicted and observed process gain and time constant were also in good agreement. As before, the difference in the scale makes distillate composition predictions appear more accurate than the bottoms composition predictions. Values for p and n were estimated from data obtained at several values Of &be Thus, the p and n values in Table 1 represent “average” values over the range of &b. The open-loop responses in Figures 6-13 are all at a constant +b. Consequently, predicted values of XD and XB are biased with respect to the measured values. As a result,

0

5000

15000 10000 Time ( s )

20000

25 00

Figure 11. Bottoms composition response when the distillate rate was decreased from 0.032to0.025 kg/s (t = 4979 s) and then increased to 0.038 kg/s (t = 15 120 8). The feed rate was constant at 0.065 kg/s and the reboiler duty at 1.47 X locl W.

0

X

Estimator Model RM Measured 0

5000

----

10000 15000 20000 25000 30000 Time (si

Figure 12. Distillate composition response when the distillate rate was increased from 0.032to 0.038 kg/s (t = 5760 8) and then decreased to 0.025 kg/s (t = 16 020 8). The feed rate was constant at 0.055 kg/s and the reboiler duty at 1.47 X 1V W.

the measured bottoms composition is always higher than the composition predicted by RM, while the predicted distillate composition is always higher than the measured value. For some Qreb values different from the value in Figures 6-13, the predicted XB values will be higher than measured XB values, and predicted XD’S will be lower than the measured XD’S.

Ind. Eng. Chem. Res., Vol. 32, No. 10,1993 2353 0.03

0.9

0.02 0.8

X I

0.01

I

II

0.7

II

Measured

6

0.6

0.00 5000

10000

15000

20000

25000

31

1500

0

Time

Time i s )

3000

4500

3000

4500

3000

4500

(6)

Figure 13. Bottoms composition response when the distillate rate was increased from 0.032 to 0.038kg/s (t = 5760 8 ) and then decreased to 0.025 kg/s (t = 16 020 e). The feed rate was constant at 0.055 kg/s

and the reboiler duty at 1.47 X I@ W.

8. Closed-Loop Tests

Since the process gain changes sign at low feed flow rates, a linear controller (e.g., PID, DMC) cannot be used to control the column. Patwardhan et al. (1992)compared, via simulation, the performance of NMPC (using a dynamic nonlinear model) and DMC with gain scheduling to control this process. The sign of the controller gain was adjusted to coincide with the sign of the process gain, based on the region of operation. They found that the closedloop reponse exhibited sustained oscillations, and the controller performance was not satisfactory. NMPC, on the other hand, was found to provide satisfactory closedloop control in simulation despite changes in the sign of the process gain. To test the performance of the NMPC scheme using the steady-state model in the presence of gain-signchanges, closed-loop tests were performed with different values for the feed rate. Figure 14shows the top and bottom product composition and reboiler duty responses for a feed rate of 0.133 kg/s when the column was taken to the set-point (xFt = 0.9410,x r t = 0.0338) from an arbitrary initial condition. The parameter n was held fixed at a nominal value of 0.01,while p was estimated on-line. The reboiler duty increased until it reached its upper bound (slightly below the flooding point). As the controlled variables approached their set points, the heat duty dropped slowly to a steady-state value. In this example, the process gain did not change sign. For a lower feed rate (0.044kg/s), the situation was entirely different. Figure 15 shows the input and output trajectories from an arbitrary initial condition for xFt = 0.8832 and = 0.0338. The reboiler duty reached its lower bound in order to drive the outputs to the set point at a fast rate. As soon as the outputs approached their set points, the reboiler duty increased to a steady-state value between the upper and lower bounds. On the basis of operating experience, the column was deemed to have reached steady state at the end of the trace in Figure 15, and the system was shut down. At this steady state, neither XD nor XB was at its set point, because at steady state, p = 50. Figure 2 shows predicted steady-state plots of X D and XB versus the reboiler duty for p = 60 and n = 0.01. Each of these plots shows an extremum at a reboiler duty of roughly 1.6 X lo5 W. They are show that the set point was infeasible according to the process model. A linear controller would have been unstable due to the incorrect sign of its gain, and would have saturated the reboiler

xrt

0.004 0

1500

Time

( 6 )

0.20

0.16

0.12

0.08 0

1500

Time ( s )

Figure 14. Closed-loop ZD, XB, and Qmb trajectories using NMPC in response to a set-point change (F = 0.133 kg/s).

duty either at its lower bound or its upper bound, therefore driving the outputs even farther from their set points. Whe NMPC was used with an additive disturbance estimator, and the controller model exhibited an extremum, an infeasible set point caused the inputs to achieve steady-state values corresponding to the manipulated variable value at the extremum. The steady-state reboiler duty was very close to the value at which the extremum in Figure 2 occurs. The present example demonstrates one of the major advantages of NMPC over linear control, viz., NMPC's ability to satisfy infeasible set points in a least-squares sense without necessarily saturating the inputs. 9. Conclusions

The distillation test system was operated at several values of the feed rate, distillate rate, and reboiler duty.

Ind. Eng. Chem. Res., Vol. 32, No. 10, 1993 iL . n"

I

,

Time ( s i

1

0.127

0

1

n n n Jr V . " "

2500

5000

7500

10000

rigorous model via reasonable assumptions, and its use was found to satisfy the memory requirement. The simplified model captured the gain-sign change exhibited by the process, which was the nonlinearity of greatest concern in this work. The steady-state model was found to provide good control of the process. Thus, in situations where gain nonlinearities lead to poor performance by linear controllers, steady-state nonlinear models may provide adequate control within the NMPC framrework. This is an important result, since the development of detailed nonlinear dynamic models may be unnecessary. It will be worthwhile to investigate the use of steady-state models for control of other industrial process units, such as tray towers, reactors, etc. When the process gain changes sign, the steady-state characteristic exhibits an extremum. This means that certain set points are infeasible. If an infeasible set point was specified as the control objective, a linear controller makes the closed loop unstable and saturates the manipulated variables. NMPC, on the other hand, recognizes the infeasible nature of the set point, and takes the process as close to the set point, in a least-squares sense, as possible. This capability of NMPC was demonstrated experimentally (see also Patwardhan and Edgar (1991)) on a packed distillation column.

Acknowledgment The authors acknowledge support from NSF Grant CBTE-8420001, the Eastman Kodak Company, and the Separations Research Program a t The University of Texas at Austin.

Time ( s i 0.2c

0.15

0.10

0.05 2500

5000

7500

10000

Time ( s )

Figure 16. Closed-loop XD, XB, and Qreb trajectories using NMPC in response to a set-point change (F = 0.044 kg/s).

The parameters p and n in the mass-transfer-coefficient expression were found to be functions of the feed and distillate flow rates. A t high F and D, the bottoms composition decreased with increasing reboiler duty. At lower F andD values, the A.%B/AQ,b and A.%D/A&b process gains changed sign as the reboiler duty was varied over the range of permissible values. Thus, the column was found to exhibit input multiplicity. The response to step changes in the feed and distillate rates predicted by the rigorous model was compared to the experimentally observed response. The rigorous model predicted XD and XB to within 0.02 mole fraction. The predicted process gains and time constants agreed well with the experimental values. The rigorous model was not suitable for on-line use, due to memory limitations on the process control computer. A simplified steady-state model was obtained from the

Nomenclature A = cross-sectional area of the column (m2) a = interfacial area per unit volume of packing (m2/m3) B = bottoms product flow rate (g-mol/@ C p = molar heat capacity (J/(g-moEK)) D = distillate rate (g-mol/s) d = additive disturbance estimate e = error between the set point and the controlled variable F = feed flow rate (g-mol/@ f = light component mole fraction in the feed H = liquid holdup in auxiliary equipment (m3) h = height of packing (m) A",, = latent heat of vaporization (J/g-mol) K = vapor-liquid equilibrium constant Ky = overall vapor-phase mass-transfer coefficient (g-mol/ (m2.s-mole fraction)) L = liquid flow rate (g-mol/s) M = number of collocation points on each finite element in time n = constant in eq 24 NC= number of collocation points in space including the two end points on each finite element N s = number of composition measurement locations p = constant in eq 24 Q = heat duty (W) 'T = sampling time T = temperature (K) t = time ( 8 ) V = vapor flow rate (g-mol/s) W = objective function weight for the estimation NLP x = light component liquid mole fraction y = light component vapor mole fraction z = distance (m) Greek Letters a= p

relative volatility of the light component

= molar density (g-moVm3)

Ind. Eng. Chem. Res., Vol. 32, No. 10, 1993 2366 $t

= holdup in packing (mVm3)

Q = objective function weight for the control NLP Subscripts A = overhead accumulator a = location in the packing B = bottom product b = location in the packing BR = bottom of the rectifying section BS = bottom of the stripping section c = location in the packing D = distillate drum = reflux drum F = feed I = between packed beds i = finite element j = collocation point k = sampling instant L = liquid 1 = lower bound m = mass transfer max = maximum min = minimum n = composition measurement location 0 = output R = rectifying section reb = reboiler ref = reference S = stripping section sump = sump TR = top of the rectifying section TS = top of the stripping section u = upper bound V = vapor 1 = component 1 2 = component 2

Superscripts

* = vapor-liquid

interface = measured value = estimated value a = location a b = location b set = set point

A

Literature Cited Alsop, A. W.; Edgar, T. F. Nonlinear Control of a High Purity Distillation Column by the use of Partially Linearized Control Variables. Comput. Chem. Eng. 1990,14,665. Arulalan, G. R.; Deshpande, P. B. Simplified Model Predictive Control. Znd. Eng. Chem. Res. 1987,26, 347. Asselmeyer, B. Optimal Control for Nonlinear Systems Calculated with Small Computers. J. Optim. Theory Appl. 1985,45,533. Bequette, B. W. Measurement Selection and Control System Design for Multivariable Interacting Processes. Ph.D. Dissertation, The University of Texas at Austin, 1986. Biegler,L. T.; Rawlings,J. B. Optimization Approaches to Nonlinear Model Predictive Control. Presented at the Fourth International Conference on Chemical Process Control, South Padre Island, TX, 1991. Bradley, K. J.; Andre, H. A Dynamic Analysis of a Packed Gas Absorber. Can. J. Chem. Eng. 1972,50,528. Bravo, J. L.; Patwardhan, A. A,; Edgar, T. F. Influence of Effective Interfacial Areas in the Operation and Control of Packed Distillation Columna. Znd. Eng. Chem. Res. 1992, 31, 604. Brengel, D. D.; Seider, W.D. Multistep Nonlinear Predictive Control. Znd. Eng. Chem. Res. 1989,28, 1812-1822. Culter, C. R.; b a k e r , B. L. Dynamic Matrix Control-A Computer Control Algorithm. Joint Automatic Control Conference Proceedings: 1980, paper WPdB. Daubert, T. E.; Danner, R. P. Data Compilation of Tables of Properties of Pure Components; AIChE: New York, 1985.

Eaton, J. W.; Rawlings,J. B. Feedback Controlof ChemicalProceeses Using On-line Optimization Techniques. Comput. Chem. Eng. 1990,87,469-479. Eaton, J. W.; Rawlings, J. B.; Edgar, T. F. Model-Predictive Control and Sensitivity Analysis for Constrained Nonlinear Processes. In ZFAC Workshop on Model-based Control; Pergamon Press: 1990; p 129. Economou, C. G.; Morari, M.; Palsson, B. 0. Internal Model Control. 5. Extension to Nonlinear Systems. Znd. Eng. Chem. Process Des. Dev. 1986,25,403. Garcia, C. E.; Morari, M. Internal Model Control. 1. A Unifying Review and Some New Results. Znd. Eng. Chem. Process Des. Dev. 1982,21, 308. Gmehling, J.; Onken, U.; Ark, W. Vapor Liquid Equilibrium Data Collection. Aliphatic Hydrocarbons, C&. In Chemistry Data Series; Behrene, D., Eckermann, R., Eds.; DECHEMA: Frankfurt, Federal Republic of Germany, 1980; Vol. I, Part 6a, pp 300-306. Gray, R. I.; Prados, J. W. The Dynamics of a Packed Gas Absorber by Frequency-Response Analysis. AIChE J. 1963,9,211. Hasdorff, L. Gradient Optimization and Nonlinear Control; Wiley-Interscience: New York, NY, 1976. Hicks, G. A.; Ray, W. H.Approximation Methodsfor Optimal Control Synthesis. Can. J. Chem. Eng. 1971, 49, 522. Hitch, D. M.;Rouseeau, R. W.; Ferrell, J. K. Simulation of ContinuousContact Separation Processes: Multicomponent Adiabatic Absorption. Znd. Eng. Chem. Process Des. Dev. 1986,25, 699. Hitch, D. M.;Rousseau, R. W.; Ferrell, J. K. Simulation of ContinuousContact Separation Processes: Unsteady-State, Multicomponent, Adiabatic Absorption. Znd. Eng. Chem. Res. 1987,26, 1092. Holland, C. D.; Liapis, A. I. Computer Methodsfor Solving Dynamic Separation Problems; McGraw-Hill: New York, NY, 1983. Jang, S.; Joseph, B.; Mukai, H. On-line Optimization of Constrained Multivariable Processes. AZChE J. 1987,33 (l),26. Jones, D. I.; Finch, J. W. Comparison of Optimization Algorithms. Znt. J. Control 1984,40, 747. Kiparissides, C.; Georgiou, A. Finite-Element Solution of Nonlinear Optimal Control Problems with a Quadratic Performance Index. Comput. Chem. Eng. 1987,11,77. Krishnamurthy, R.; Taylor, R. Simulation of Packed Distillation and Absorption Columns. Znd. Eng. Chem. Process Des. Dev. 1985,24, 513. Lasdon, L. S.; Warren, A. D. ‘GRG2 User’s Guide”; Department of General Business, School of Business Administration, The University of Texas at Austin, Austin, TX, 1986. Li, W. C.; Biegler, L. T. Process Control Strategies for Constrained Nonlinear Systems. Znd. Eng. Chem. Res. 1988,27,1421-1433. Li, W. C.; Biegler, L. T. Multistep, Newton-Type Control Strategies for Constrained, Nonlinear Processes. Chem.Eng. Res. Dev. 1989, 67,562-577. Li, W. C.; Biegler, L. T. Newton-Type Controllers for Constrained Nonlinear Processes with Uncertainty. Zng.Eng. Chem. Res. 1990, 29,1647-1657. Li, W. C.; Biegler, L. T.; Economou, C. G.; Morari, M. A Constrained Pseudo-NewtonControl Strategy for Nonlinear Systems. Comput. Chem. Eng. 1990,14 (4/5), 451-468. Liebman, M. J. Reconciliation of Process Measurements using Statistical and Nonlinear Programming Techniques. Ph.D. Dissertation. The University of Texas at Austin, Austin, TX, 1991. Maddox, R. N. Gas Absorption. In Chemical Engineers’ Handbook; Perry, R. H., Chilton, C. H., Eds., McGraw-Hik New York, NY, 1973. McDaniel, R.; Holland, C. D. Modeling of Packed Absorbers at Unsteady State Operation-IV. Chem. Eng. Sci. 1970,25,12831296. Morshedi, A. M. Universal Dynamic Matrix Control. Chemical Process Control Conference ZZZ, Aailomar, CA; CACHE Corp.: Austin, TX, 1986. Patwardhan, A. A. Modeling and Control of a Packed Distillation Column. Ph.D. Dissertation, The University of Texas at Austin, Austin, TX, 1991. Patwardhan, A. A.; Edgar, T. F. Dynamics of Packed Columns with Low Holdups. Presented at the AIChE Spring Meeting, Houston, TX, 1989; paper 35a. Patwardhan, A. A.; Edgar, T. F. Nonlinear Model-Predictive Control of a Packed Distillation Column. Presented at the American Control Conference, Boston, MA, 1991. Patwardhan, A. A.; Rawlings, J. B.; Edgar, T. F. Model Predictive Control of Nonlinear Processes in the Presence of Constraints. Presented at the AIChE Annual Meeting, Washington, DC, 1988.

2356 Ind. Eng. Chem. Res., Vol. 32,No. 10, 1993 Patwardhan, A. A.; Rawlings, J. B.; Edgar, T. F. Nonlinear Model Predictive Control. Chem. Eng. Commun. 1990,87, 123-141. Patwardhan, A. A.; Wright, G. T.; Edgar, T. F. Nonlinear Model Predictive Control of Distributed-Parameter Systems. Chem.Eng. S C ~1992,47 . (4), 721-735. Peterson, T. E.; Arkun, Y .;Schork, F. J. Nonlinear Predictive Control of a Semi Batch Polymerization Reactor by an Extended DMC. Presented at the American Control Conference, Pittsburgh, PA, 1989; paper TP4. Petzold, L. F. ‘A Description of DASSL A DifferentWAlgebraic System Solver”;Technical Report SAND82-8637;Sandia National Laboratories, Albuquerque, NM, and Livermore, CA, 1982. Richalet, J.; Testud, J. L.; Papon, J. Model Predictive Heuristic Control: Applications to Industrial Processes. Autornatica 1978, 14,413. Sakata, N.; Prados, J. W. The Dynamics of a Packed Gas Absorber by the Pulse Response Technique. AZChE J. 1972,18, 572. Sargent, R. W. H.; Sullivan, G. R. The Development of an Efficient Optimal Control Package. In Optimization Techniques, Proceedings of the Eighth ZFIP Conference on Optimization Techniques, Wurzburg, 1977; J., Stoer, Ed.; Springer-Verlag, Berlin, 1978. Srivastava, R. K.; Joseph, B. Simulation of Packed Bed Separation Processes using Orthogonal Collocation. Comput. Chem. Eng. 1984, 8, 43. Suenson,M. M.; Georgakii,C.; Evans, L. B. Steady-state and Dynamic Modeling of a Gas Absorber-stripper System. Znd. Eng. Chem. Fundam. 1986,24, 288.

Tommasi, G.; Rice, P. Dynamics of Packed Tower Dietillation. Znd. Eng. Chem. Process Des. Dev. 1970,9,9. Treybal, R. E. Adiabatic Gas Absorption and Stripping in Packed Towers. Znd. Eng. Chem. 1969.61, 36. Villadaen,J. V.; Stewart, W. E. Solution of Boundary-value Problems by Orthogonal Collocation. Chem. Eng. Sci. 1967,22,1483. von Rosenberg, D. U.; Hadi, M. S. Numerical Solution of Multicomponent, Packed Tower Distillation Problems. Chem. Eng. Commun. 1980,4, 313. von Stockar, U.; Wilke, C. R. Rigorous and Short-cut Design Calculations for Gas Absorption Involving Large Heat Effects. 1. A New Computational Method for Packed Absorbers. Znd. Eng. Chem. Fundam. 1977,16,88. Wright, G. T.; Edgar, T. F. Adaptive Nonlinear Control of a Laboratory Water-Gas Shift Reactor. Preeentedat ITAC91, IFAC International Symposium on Intelligent Tuning and Adaptive Control, Singapore, 1991. Received for review August 10, 1992 Accepted February 9, 1993.

* Abstract published in Advance ACS Abstracts, August 15, 1993.