How to Use Simplified Dynamics in Model Predictive Control of

In this paper, issues associated with the application of model predictive control algorithms to the product quality control of superfractionator colum...
0 downloads 0 Views 375KB Size
Download PDF
Ind. Eng. Chem. Res. 2005, 44, 2687-2696

2687

How to Use Simplified Dynamics in Model Predictive Control of Superfractionators Gabriele Pannocchia* and Alessandro Brambilla Department of Chemical Engineering, University of Pisa, Via Diotisalvi, 2 - 56126 Pisa, Italy

In this paper, issues associated with the application of model predictive control algorithms to the product quality control of superfractionator columns are addressed. Full-order and reducedorder (pseudo-ramp) models are compared in different predictive control algorithms, and the effects of different disturbance models on the closed-loop performance are emphasized. The results show that the reduced-order model, which is easier to identify, can give excellent results when it is used along with an effective disturbance model and can be implemented with a relatively small sampling time. The results also show that the output disturbance model typically used by industrial algorithms, is the source of poor closed-loop performance, and although the input disturbance model represents a better choice, the rotation factor disturbance model used with the reduced model can represent a simple and relatively effective practical alternative. 1. Introduction Distillation is one of the most important separation methods used in the process industries, and it often represents a significant challenge from a control point of view, as confirmed by the substantial industrial and academic research interest. The main driving force in distillation is the boiling point difference between the species to be separated, better defined in terms of relative volatility. In some important industrial cases, e.g., the separation of propane/propylene, ethane/ ethylene, and n-butane/i-butane, this driving force can be fairly small, and in order to achieve high-purity products, columns with a large number of trays (hundreds) and high internal flows (reflux ratios greater than 10) are used. Such distillation columns are often referred to as “superfractionators”. Effective control of superfractionators can be difficult for a number of reasons. First, superfractionators are characterized by very slow dynamics regarding the product compositions. Then, the large internal flows make bottom or overhead drum level control critical when the products are their manipulated variables. Moreover, depending on which quality control structure is chosen, interaction and illconditioning problems can be very large,1,2 and a number of different control structures (including the “unusual” DB one) have been considered in the literature.2-4 Model predictive control (MPC) currently represents the most widely used advanced multivariable control technique in the process industries and, in particular, in refinery plants,5 where the number of (controlled and manipulated) variables can be large (and often not equal); therefore, traditional multiloop single-input/ single-output (SISO) controllers do not allow optimal operation of the process, especially in the presence of constraints. An advanced multivariable control system is generally organized in two hierarchical levels: (1) At the lower level, a number of SISO proportionalintegral (PI) or proportional-integral-derivative (PID) * To whom correspondence should be addressed. E-mail: [email protected]. Fax: +39 050 511266. Tel.: +39 050 511238.

controllers guarantee automatic operation of the process. These loops usually include flow rate, level, pressure, and temperature controllers. (2) At the higher level, a multivariable predictive controller manipulates (every minute or so) the set points of the lower-level controllers to achieve quality and optimization goals respecting unit constraints. Because of the presence of the low-level PI (PID) controllers, most plant dynamics are stable, and this explains (at least in part) why convolution models are still so popular in industrial MPC algorithms, even if they cannot describe open-loop unstable or integrating dynamics or a large number of coefficients is required for slow processes. Sometimes, it can be convenient to “open” some level controllers (especially when the corresponding inventories are large or when more than one controller is preferable for a plant section) and let the predictive controller directly control those levels.6 Specific ad hoc features have been developed by practitioners to handle integrating systems within the traditional convolution model framework. MPC algorithms have been developed and applied in many areas, including distillation control (see refs 7-13). The major problems associated with the control of superfractionators with MPC algorithms are due to the (very) large time constants that arise in the key controlled variables (10-20 h). First, the identification tests required to develop a process model are very long, and often, the identified model can be unreliable because of the presence of disturbances. Then, the choice of the controller’s main parameters (such as the sampling time and the prediction horizon) can be difficult because other variables with faster dynamics (such as temperatures, levels, valve positions, etc.) must also be controlled. Finally, the performance obtained with common industrial algorithms, such as DMC, can be poor in the presence of slow disturbances, as reported and criticized in the literature.14,15 The objective of this paper, therefore, is essentially two-fold: (1) to address these control issues and propose simple and effective (i.e. applicable) solutions within the framework of currently available MPC technologies and (2) to show how alternative MPC algorithms can further improve the control performance.

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

2688

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

the control inputs are bounded

umin e u e umax

To guarantee offset-free control of the outputs in the presence of unmeasured disturbances and/or plantmodel mismatch, the process model (1) is augmented with fictitious integrating disturbances that are estimated, at each sampling time, from the difference between the actual measured outputs and those predicted by the augmented model. The number of integrating disturbances must equal the number of measured outputs,21,22 independently of the number of variables in which offset-free control is to be obtained.21 The augmented model can be described by the following general form21

Figure 1. MPC scheme.

2. Model Predictive Control MPC algorithms use a (linear) process model to predict the future plant behavior and compute a sequence of future control inputs that minimize a (quadratic) performance objective function while respecting the process constraints.5,16-19 Industrial MPC algorithms such as DMC use multivariable finite-step response (FSR) models, whereas others such as RMPCT or PFC use multivariable autoregressive with external input (ARX) or transfer function (TF) models. Most industrial algorithms are organized in three basic modules, as depicted in Figure 1: (1) a disturbance estimator, which computes a correction term that is added to the model prediction; (2) a steady-state optimization module, which computes the targets for all controlled and manipulated variables; and (3) a dynamic optimization module, which computes a sequence of future control inputs in order to minimize an objective function over some future prediction horizon. In this paper, the same three modules are considered in a fairly general setting so that different algorithms can be compared. In particular, the DMC (currently DMC-plus) controller is used as a reference because it is one of the most widely used predictive controllers in refinery plants, although alternative algorithms are considered and possible improvements in performance are emphasized. 2.1. Internal Model and Disturbance Estimation. It has been predicted5 that next-generation MPC algorithms will use state-space models, which are also the most popular kind of system models considered in academic studies. Because any particular model (FSR, FIR, ARX, TF, etc.) can be converted into an equivalent (or approximated) state-space model (see ref 19), in this paper, an “internal” state-space discrete time-invariant model in the following form is considered

xk+1 ) Axk + Buk yk ) Cxk

(1)

in which x ∈ Rn is the state vector, u ∈ Rm is the manipulated-variable vector, and y ∈ Rp is the controlledvariable vector. It is assumed that (A, B) is stabilizable and (C, A) is detectable and the following condition is satisfied

rank

[

]

I-A B )n+p C 0

(3)

(2)

which ensures that offset-free control of the outputs y is possible (see refs 20 and 21). It is also assumed that

[ ] [ ][ ] [ ] [ ]

xk+1 B A Bd xk dk+1 ) 0 I dk + 0 uk x yk ) [C Cd ] dk k

(4)

in which d ∈ Rp is the additional disturbance. Such an augmented system is detectable if and only if the following rank condition holds21,22

rank

[

]

I - A -Bd Cd ) n + p C

(5)

Provided that the augmented system is detectable, the state and the unmeasured disturbance can be estimated from the plant measurement by means of a (steadystate) Kalman filter

xˆ k|k ) xˆ k|k-1 + Lxek dˆ k|k ) dˆ k|k-1 + Ldek

(6)

in which Lx ∈ Rn×p and Ld ∈ Rp×p are constant matrices that can be computed from a discrete algebraic Riccati equation and ek is the prediction error

ek ) yk - (Cxˆ k|k-1 + Cddˆ k|k-1)

(7)

By combining the augmented model (eqs 4) with the estimator (eqs 6), the so-called “innovation form” is obtained

[ ] [ ][ ] [ ] [

]

xˆ k+1|k ALx + BdLd B A Bd xˆ k|k-1 ek dˆ k+1|k ) 0 I dˆ k|k-1 + 0 uk + Ld (8)

in which it clearly appears that the disturbance dˆ integrates the prediction error. The most common industrial choice of augmented model is the so-called “output disturbance model”, obtained from the general model in eqs 4 by setting

Bd ) 0, Cd ) I

(9)

which corresponds to modeling the disturbance as a constant step acting on the output. Most industrial algorithms do not compute the filter gain matrices from a Riccati equation but rather simply use the following values

Lx ) 0, Ld ) I

(10)

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

variable systems with more complex dynamics than a pure integrating system, the rotation factor disturbance model can also be fitted in the general model in eqs 4. The specific choices for the disturbance model matrices (Bd, Cd) are described in the case study, and as shown, the disturbance model affects significantly the closed-loop performance. Indeed, a number of researchers have criticized the output disturbance model used by DMC as a source of the sluggish rejection of slow disturbances14,15,22 and poor robustness to model errors for multivariable ill-conditioned processes.15,23 Several studies22,23 have pointed out that the input disturbance model, which, for square systems, can be obtained by setting

Bd ) B, Cd ) 0 Figure 2. Interpretation of the rotation factor disturbance model.

which lead to the following innovation form

xˆ k+1|k ) Axˆ k|k-1 + Buk dˆ k+1|k ) dˆ k|k ) yk - Cxˆ k|k-1

(11)

It should be noticed that industrial controllers typically allow one to update the disturbance with a fraction of the prediction error, to reduce the effects of output noise. The choice of the output disturbance model and related estimator is acceptable only for asymptotically stable system models. In particular, for integrating systems, the output disturbance model in eqs 9 is not detectable because eq 5 does not hold and, hence, it cannot be used to remove offset for integrating systems. This limitation has led practitioners to develop a number of ad hoc fixes. One such fix is called the “rotation factor” disturbance model, which can be described for a pure-integrating SISO system y(s) ) (k/s)u(s) by the following augmented state-space system (see ref 22)

[ ] [ ][ ] [ ] [ ]

xk+1 1 R xk b dk+1 + 0 1 dk + 0 uk x yk ) [1 1 ] dk k

dˆ k+1|k ) dˆ k|k ) yk - xˆ k|k-1

improves the feedback properties of MPC regulators. It is also important to remark that, for a given disturbance model, the estimator tuning, i.e., the choice of the filter gain matrices Lx and Ld, is significant as well. In fact, because two disturbance models (with the same number of integrators) are two different realizations of the same input/output process, it is possible to show24 that there exist filter gain matrices for each disturbance model that achieve the same input/output behavior. 2.2. Steady-State Optimization. Using the estimate of the current disturbance and assuming, for simplicity of presentation, that set points are defined for all controlled variables (even though this is usually not the case for large multivariable plants), the steady-state optimization module computes the steady-state targets that possibly drive the controlled variables to their set points while respecting the constraints. The most common objective functions used in the steady-state optimization modules are linear or quadratic. In this work, given the current set point for the controlled variables yj and the current disturbance estimate dˆ k|k at each sampling time, the following optimization problem is solved25

j - us)TRs(u j - us) (xs(k), us(k)) ) argmin(u xs,us

(12)

(15a)

subject to

in which b ) k/Ts (Ts is the sampling time). In this disturbance model, the error between the measured output and the model prediction is assumed to be caused partially by a constant output disturbance and partially by an input disturbance. Similarly to the output disturbance model, the common industrial choice of estimator is still given by eqs 11. It is clear that, with the rotation factor disturbance model, the output prediction is rotated by a factor proportional to R and not simply translated by a constant bias as in the output disturbance model. In fact, the corresponding innovation form can be written as

xˆ k+1|k ) xˆ k|k-1 + buk + R(yk - xˆ k|k-1)

(14)

(13)

A graphical interpretation of the rotation factor disturbance model is shown in Figure 2. For multi-

xs ) Axs + Bus + Bddˆ k|k

(15b)

yj ) Cxs + Cddˆ k|k

(15c)

umin e us e umax

(15d)

in which u j represents the input desired set point (often not specified and frequently assumed to be u j ) 0) and Rs is a positive definite matrix. If this optimization problem is infeasible, it means that, given the current disturbance estimate and the current set point, it is not possible to obtain offset-free control in all controlled variables. In such cases, one can solve another optimization problem to find the steady-state targets that minimize the offset25

(xs(k), us(k)) ) argmin(yj - ys)TQs(yj - ys) xs,us

(16a)

2690

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

subject to

xs ) Axs + Bus + Bddˆ k|k

(16b)

ys ) Cxs + Cddˆ k|k

(16c)

umin e us e umax

(16d)

in which Qs is a positive definite matrix. The optimal vectors (xs(k), us(k)) obtained from problem 15 (or from problem 16 when problem 15 is infeasible) represent the current state and input targets, and the corresponding output target is given by

ys(k) ) Cxs(k) + Cddˆ k|k

(17)

2.3. Dynamic Optimization. Given the current targets computed by the steady-state optimization module, the goal of the dynamic optimization module is to compute a sequence of future control inputs that are feasible with respect to the constraints and minimize an objective function that includes both the transient error of the controlled variables from their target and the control effort. In this work, the following optimization problem is considered L-1

v*(k) ) argmin v(k)

∑ j)0

N-1

(vj - vj-1)TS(vj - vj-1) +

∑ j)1

(ys(k) - yˆ j)TQ(ys(k) - yˆ j) + (xs(k) - wN)TP(xs(k) - wN)

Figure 3. Column layout.

(18a)

Given the solution vector v*(k), all predictive controllers use a receding-horizon implementation, that is, the current input injected into the plant is

w0 ) xˆ k|k, v-1 ) uk-1

(18b)

uk ) v/0(k)

wj+1 ) Awj + Bvj + Bddˆ k|k

(18c)

yˆ j ) Cwj + Cddˆ k|k

(18d)

and all modules (estimator, steady-state optimization, and dynamic optimization) are re-executed at the next sampling time.

umin e vj e umax

(18e)

subject to

(19)

3. Case Study

/ ) is the optimal input in which v*(k) ) (v/0 v/1 ... vL-1 sequence; L e N represents positive integers, usually referred to as control (L) and prediction (N) horizons; and S, Q, and P are symmetric positive definite matrices. A number of different algorithms can be fitted in this general optimization problem. (1) DMC uses P ) CTQC, and hence, the last term in eq 18a is simply equal to (ys(k) - yˆ N)TQ(ys(k) - yˆ N). Typical settings are L , N. (2) Rawlings and Muske25,26 choose L ) N and P as the solution to the Lyapunov equation P ) CTQC + ATPA. This choice is known to guarantee constrained nominal stability independently of N, Q, and S for openloop stable plants. (3) Chmielewski and Manousiouthakis27 and Scokaert and Rawlings28 also choose L ) N, but chose P as the solution to the Riccati equation associated with problem 18 (see ref 29). This choice guarantees constrained nominal stability and optimality for any plant (stable or unstable) when N is sufficiently large.

3.1. Column Characteristics and Nonlinear Model. A distillation column for the separation of propylene and propane is chosen as an example, because this is the most common superfractionator present in crude refineries or petrochemical plants. The column layout is depicted in Figure 3, in which only the basic (inventory) control loops are shown, and the main column characteristics are reported in Table 1. Even though several studies (see, e.g., ref 2) suggest the (L, B) configuration as the most appropriate for superfractionators for which the top product is more important, as in the case of propylene/propane splitter, in the present work, the (D, V) configuration is chosen because it is probably the most widely used control structure in industrial applications of this type. Plant data are surrogated by means of the results obtained from a detailed nonlinear model of the column, which is described below. Because the product composition dynamics present very large time constants, to simplify the model, it is reasonable to assume that level control of the overhead drum and column bottom is perfect. A detailed nonlinear

Ind. Eng. Chem. Res., Vol. 44, No. 8, 2005 2691 Table 1. Column Characteristics parameter

symbol

value (nominal)

no. of theoretical stages feed theoretical stage feed rate feed propylene mole fraction feed liquid fraction distillate rate bottom rate reflux rate boil-up rate relative volatility top product propylene mole fraction bottom product propylene mole fraction reflux drum holdup column trays' holdup bottom sink holdup

N Nf F zF qF D B L V R xD

170 97 5.17 0.719 0 3.70 1.47 58.57 57.09 1.11 0.995

xB

0.025

MT Mi MB

253 3.64 78.2

units

kmol/min kmol/min kmol/min kmol/min kmol/min

kmol kmol kmol

model of the column has been developed under the following standing assumptions: (1) constant vapor and liquid mole rates in each column section and (2) constant relative volatility, throughout the column. These assumptions are justified by the fact that propane and propylene have similar thermodynamic properties. Thus, the nonlinear model consists of N ) 170 ordinary differential equations (ODEs), each describing the material balance of propylene in each theoretical stage, and this nonlinear model is used throughout this work only to describe the true process. All simulations were performed using Octave (freely available at http:// www.octave.org), and integration of the ODE system was accomplished with the embedded Fortran subroutine LSODE. Because the inventory control loops are assumed to be perfect, from a control point of view, the process has two manipulated inputs (distillate rate D and boil-up rate V) and two controlled outputs (top product mole fraction of propylene xD and bottom product mole fraction of propylene xB). The following input constraints are considered

[

] [] [ ]

2.90 D 4.5 e e . 42.18 V 72

(20)

Direct measurements of product compositions are required, as the temperature changes are so small (because of the close boiling points and high purity) that a property estimator would not be reliable. It is assumed that a composition measurement is available every 20 min and the composition measurements are delayed by θ ) 20 min. 3.2. Identification of Linear Models. To derive a model that represents the response of the controlled variables (CVs) with respect to variations of the manipulated variables (MVs), a number of identification tests, typically “step tests”, are carried out. Usually, each MV is stepped individually (i.e., one at a time), and a single-input/multi-output (SIMO) identification procedure is carried out using least-squares or other techniques (e.g., subspace identification).30 To capture the steady-state characteristics of the process, typically, one has to wait until the slowest CV has approximately reached a steady value. By using the nonlinear model as the plant, under the assumption that no disturbances affect the system during the identification tests, after averaging the process response between positive and negative steps in the MVs, this identification approach

Figure 4. Feed rate and composition disturbances used in simulations.

has led to the following linear model, for simplicity shown in the Laplace domain

M1:

[

[ ]

xD(s) ) xB(s)

-0.101(9.62s + 1) 5.82 × 10-4(9.09s + 1) (1551s + 1)(92.1s + 1) (104.8s + 1)(51.4s + 1) -0.418(9.05s + 1) -1.50 × 10-3(8.29s + 1) (645s + 1)(36.5s + 1) (88.0s + 1)(22.4s + 1) e-20s

[ ] D(s) V(s)

]

(21)

in which the time constants are in minutes. From this model, two well-known characteristics appear: (1) The responses, especially for the top product composition (the most pure), have very large time constants. (2) The responses to changes in the distillate flow rate are slower but with higher gain than the responses to changes in the boil-up rate. It is important to point out that such a model would be very difficult to obtain from true plant tests for the following reasons: (1) The duration of a single step is very long (at least about 2 days in the distillate and 1 day in the vapor boil-up). (2) During the long transient period necessary to reach approximately the new steady state, it is likely that disturbances will enter the plant, thus making the identified model quite different from the true plant behavior. All of these problems are emphasized by the multivariable nature of the problem, that is, more than one variable needs to be stepped, and the overall identification time can be very large. Moreover, in the presence of relevant disturbances, it is possible to obtain inconsistent models, that is, models with a wrong “multivariable gain” sign. It is interesting to note that a recent industrial MPC implementation without plant tests (with models obtained by simulation) has been reported for two superfractionators.13 For these reasons, as well as other control issues discussed in the next paragraphs, an alternative, easier-

2692

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

Figure 5. Closed-loop simulation results (outputs) obtained with sampling time Ts ) 20 min for four predictive controllers: s, DMC 1; - -, DMC 2; - ‚ -, MPC 1; ‚‚‚, MPC 2.

to-identify model is also considered. In this model, the plant response for all controlled variables is approximated as a “pseudo-ramp”, as an integrating system with delay

M2:

[ ] [

xD(s) ) xB(s)

] [ ]

-4.68 × 10-5 2.49 × 10-6 e-20s D(s) -4.59 × 10-4 -9.15 × 10-6 s V(s)

(22)

This model can be easily identified with shorter step tests than the full model, because it is not required to wait until the controlled variables have reached a steady state. In fact, the duration of the steps used for identifying this model was 3 h for the distillate rate and 2 h for the vapor boil-up rate. 3.3. Case 1: Comparison of Controllers with Large Sampling Times. Given the slow response of the process and the fact that the product composition measurements are available only every 20 min, it would be reasonable to use a relatively large sampling time, such as 20 min. In this way, commonly used prediction and control horizons (typical industrial values of the horizons are L ) 8-15 and N ) 60-120) are significant with respect to the slowest dynamics. However, it must be kept in mind that, in practical situations, the predictive controller would be used to control other variables with a much faster dynamics as well (such as flow rates, levels, valve positions, etc.). To keep these fast variables under effective control, one must use a

Figure 6. Closed-loop simulation results (inputs) obtained with sampling time Ts ) 20 min for four predictive controllers: s, DMC 1; - -, DMC 2; - ‚ -, MPC 1; ‚‚‚, MPC 2.

relatively short sampling time, such as 1 or 2 min. Nonetheless, as a benchmark, in this subsection, a sampling time of Ts ) 20 min is chosen, whereas in the next subsection, results will be compared with those obtained for a sampling time of Ts ) 2 min. Four predictive controllers are considered: (1) DMC 1: controller based on the (state-space) full model M1 with output disturbance model eq 11 and horizons of L ) 10 and N ) 90. (2) DMC 2: controller based on the (state-space) pseudo-ramp model M2 with rotation factor disturbance model 13 (R ) 0.2 for both integrating outputs) and horizons of L ) 10 and N ) 90. (3) MPC 1: controller based on the (state-space) full model M1 with input disturbance model 14 and horizons of L ) N ) 10. (4) MPC 2: controller based on the (state-space) pseudo-ramp model M2 with input disturbance model 14, deadbeat tuning, and horizons of L ) N ) 10. All controllers use the same Q and S weight matrices in the optimization problem (eq 18)

Q)

[ ]

[

16 0 312.7 0 × 104, S ) 0 4 0 0.888

]

The terminal penalty is P ) CTQC for the DMC 1 and DMC 2, whereas it is chosen as the solution to the corresponding Riccati equation for MPC 1 and MPC 2. The estimator gain matrices for MPC 1 and MPC 2 are computed as the steady-state Kalman filter for the corresponding augmented system in eqs 4 assuming a state noise covariance of Qx ) 10-8In, a disturbance

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

Figure 7. Closed-loop simulation results (prediction errors) obtained with sampling time Ts ) 20 min for four predictive controllers: s, DMC 1; - -, DMC 2; - ‚ -, MPC 1; ‚‚‚, MPC 2. Table 2. Closed-Loop Objective Function in Different Cases case 1 case 2 case 2 (noisy)

DMC 1

DMC 2

MPC 1

MPC 2

2485 2992 -

420 383 403

215 118 -

245 137 169

noise covariance of Qd ) Ip, and an output noise covariance of Rv ) 10-7Ip (where Is is the identity matrix of dimension s). The simulation results for the rejection of (filtered) step disturbances in the feed flow rate and feed composition, as depicted in Figure 4, are reported in Figures 5-7 (outputs, inputs, and prediction errors, respectively). To have a quantitative measure of the closedloop performance, in Table 2, the value of the closedloop objective function, i.e. ∞

Φ)

(yj - yk)T Q(yj - yk) + ∑ k)0 (uk - uk-1)T S(uk - uk-1) (23)

is reported for each controller. To clarify the effect of the rotation factor on the closedloop performance, the simulation results obtained with DMC 2 based on different rotation factor values (as reported in the figure caption) are presented in Figure 8 (outputs only). 3.4. Case 2: Comparison of Controllers with Short Sampling Times. In this subsection, a shorter sampling time of Ts ) 2 min, close to what one would

Figure 8. Closed-loop simulation results (outputs) obtained with sampling time Ts ) 20 min for DMC 2 based on different rotation factor values: R ) s, 0.3; - -, 0.2; - ‚ -, 0.1; ‚‚‚, 0.01.

actually use in practice, is chosen. As in the previous case, four predictive controllers are compared: (1) DMC 1: same as in case 1, but with horizons of L ) 15 and N ) 90. (2) DMC 2: same as in case 1, but with horizons of L ) 15 and N ) 90. (3) MPC 1: same as in case 1, but with horizons of L ) N ) 15. (4) MPC 2: same as in case 1, but with horizons of L ) N ) 15. The same regulator weight matrices as in case 1 are chosen, but the move suppression penalty, S, is increased by a factor of 100 (because the sampling time is shorter). Notice that, in this case, the disturbance and, consequently, the steady-state targets are updated only every 10 sampling times, i.e., only every 20 min when a new measurement is available. Strictly speaking, this means that for the computation of the filter gain matrices, one must consider

[ ] A Bd 0 I

10

as the augmented dynamic matrix. This also means that, to have the same model correction, a rotation factor of 0.2 for case 1 corresponds to 0.02 in this case. The simulation results for the rejection of the same disturbances in the feed flow rate and feed composition as in the previous case are reported in Figures 9 and 10 (outputs and inputs, respectively). The corresponding closed-loop objective function is also reported in Table 2. Notice that, for a clear comparison with case 1, the move suppression penalty S used in eq 23 is the same

2694

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

Figure 9. Closed-loop simulation results (outputs) obtained with sampling time Ts ) 2 min for four predictive controllers: s, DMC 1; - -, DMC 2; - ‚ -, MPC 1; ‚‚‚, MPC 2.

Figure 10. Closed-loop simulation results (inputs) obtained with sampling time Ts ) 2 min for four predictive controllers: s, DMC 1; - -, DMC 2; - ‚ -, MPC 1; ‚‚‚, MPC 2.

for all cases (i.e., not multiplied by 100 as that actually used in the optimization module), and consequently, only inputs every 20 min are used in computing the difference uk - uk-1. A potential objection to the use of a rotation factor or an input disturbance model with respect to the output disturbance model could be made with regard of the output noise sensitivity of these disturbance models. To clarify this point, the simulation results in the presence of random output noise (normal distribution with zero mean and standard deviations of 0.0001 and 0.001 for xD and xB, respectively) are shown in Figures 11 and 12 (outputs and inputs, respectively), and the corresponding closed-loop objective functions are reported in Table 2. Notice that the same random output noise sequence has been used for all controllers, and for clarity of presentation, only results obtained with DMC 2 and MPC 2 are shown.

cantly. The slow correction of the output disturbance model is also emphasized in Figure 7 by the prediction errors that, after a disturbance affects the process, do not go to zero quickly. It is also interesting to see that MPC 1, which is based on the full model like DMC 1 but with an input disturbance model, shows good performance (actually, it shows the lowest closed-loop objective function; see Table 2), and this is more evidence that it is the output disturbance model that makes DMC 1 behave so poorly. It is also important to notice that MPC 1 and MPC 2 both perform very well and that the performance of MPC 2 is only slightly worse than that of MPC 1, which confirms that the approximation of the process dynamics as pseudo-ramps is valid. It also shows (see, in particular, the prediction errors in Figure 7) that the input disturbance model is an effective choice for dealing with unmeasured disturbances and plant-model mismatch. Moreover, by comparing DMC 2 and MPC 2 (see, in particular, Table 2) it can be seen that the use of the “pure” input disturbance model is a better choice than the rotation factor model (even though the improvement in performance is not as large as it was between DMC 1 and DMC 2). The results in the presence of output noise show a minor performance degradation. For the rotation factor model, the sensitivity to noise is increased when a larger value of R is used, which, on the other hand, would improve the load disturbance rejection performance as shown in Figure 8. From these results, it is clear that, when the rotation factor becomes too small, sluggish disturbance rejection similar to that of the output

4. Discussion The results presented in the previous section merit appropriate comments. First, it is important to emphasize that the use of a ramp model improved the performance of the DMC controller significantly, as shown in Figures 5, 6 and 9, 10 and also in Table 2. Indeed the closed-loop performance of DMC 1 is poor, and the source of this behavior can be found in the output disturbance model. Given the slow (process and disturbance) dynamics, the output disturbance model correction is too weak, so that the controlled variables deviate from the set point signifi-

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

Figure 11. Closed-loop simulation results (outputs) obtained with sampling time Ts ) 2 min in the presence of output noise for two predictive controllers: - -, DMC 2; ‚‚‚, MPC 2.

Figure 12. Closed-loop simulation results (inputs) obtained with sampling time Ts ) 2 min in the presence of output noise for two predictive controllers: - -, DMC 2; ‚‚‚, MPC 2.

disturbance model is obtained (indeed, if a zero rotation factor is used, the corresponding augmented system is not detectable, and the controller would not remove steady-state offset). For the input disturbance model instead, resiliency to output noise and disturbance rejection performance can be more effectively traded off by adjusting the estimator tuning (even though this goes beyond the scope of the present work). One last important consideration concerns the effect of the sampling time. For controllers based on the ramp model (DMC 2 and MPC 2), the use of a shorter sampling time improves the performance, and this is in agreement with what one would normally expect. For controllers based on the full model (DMC 1 and MPC 1) instead, the effect of a shorter sampling is different: the performance of DMC 1 deteriorates, whereas the performance of MPC 1 improves. These results can be explained by the fact that MPC 1 uses a terminal penalty that “approximates” the optimal infinite-horizon cost (see refs 28 and 29), and hence, the control sequence computed in open-loop control is close to closed-loop optimality. In DMC 1, when a shorter sampling time is used, the prediction horizon is too short to capture the steady-state process characteristics, and therefore, the computed control sequence is far from closed-loop optimality. It is important to notice that the prediction horizon (N) used by MPC 1 and MPC 2 is much smaller than that used by DMC 1 and DMC 2, and this shows how a properly chosen terminal penalty can make the predictive controller easier to implement (with a smaller prediction horizon) and more effective at the same time. Furthermore, the use of FSR models (as in several

industrial algorithms) would also introduce a model truncation error when a short sampling time is used. 5. Conclusions The product quality control of superfractionators with model predictive control algorithms has been addressed in this paper. The main difficulties that arise when designing an MPC regulator for a superfractionator column are due to the large time constants that characterize the product quality response to changes in the manipulated variables. Large time constants make the identification of a full model for the controller very difficult and time-consuming, and the choice of important controller parameters (such the sampling time and the prediction horizon) can be complicated by the presence of other controlled variables with faster dynamics. The use of a “reduced” model that is easier to identify, represented by a pseudo-ramp (an integrating system with delay), has been tested for different types of model predictive controllers. The results show that the use of a reduced model can give performance comparable to that obtained with the full model when it is implemented with an effective disturbance model and that the reduced model can be used with a reasonably small sampling time. The “rotation factor” disturbance model has been considered with the reduced model, and results show that, although the use of the input disturbance model represents a better choice, the rotation factor disturbance model can represent a simple and relatively effective practical alternative.

2696

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

Acknowledgment We thank Harmon Ray for his contributions to our research community. Literature Cited (1) Skogestad, S. Dynamics and Control of Distillation Columns: A Tutorial Introduction. Trans. Inst. Chem. Eng. 1997, 75A, 539-562. (2) Hurowitz, S.; Anderson, J.; Duvall, M.; Riggs, J. B. Distillation Control Configuration Selection. J. Process Control 2003, 13, 357-362. (3) Finco, M. V.; Luyben, W. L.; Polleck, R. E. Control of Distillation Columns with Low Relative Volatilities. Ind. Eng. Chem. Res. 1989, 28, 75-83. (4) Skogestad, S.; Jacobsen, E. W.; Morari, M. Inadequacy of Steady-State Analysis for Feedback Control: Distillate-Bottom Control of Distillation Columns. Ind. Eng. Chem. Res. 1990, 29, 2339-2346. (5) Qin, S. J.; Badgwell, T. A. A Survey of Industrial Model Predictive Control Technology. Control Eng. Pract. 2003, 11, 733764. (6) Huang, H.; Riggs, J. B. Including Levels in MPC to Improve Distillation Control. Ind. Eng. Chem. Res. 2002, 41, 4048-4053. (7) McDonald, K. A.; McAvoy, T. J. Application of Dynamic Matrix Control to Moderate and High-Purity Distillation Towers. Ind. Eng. Chem. Res. 1987, 26, 1011-1018. (8) Georgiou, A.; Georgakis, C.; Luyben, W. L. Nonlinear Dynamic Matrix Control for High-Purity Distillation Columns. AIChE J. 1988, 34, 1287-1298. (9) Yu, Z. H.; Li, W.; Lee, J. H.; Morari, M. State Estimation Based Model Predictive Control Applied to Shell Control Problem: A Case Study. Chem. Eng. Sci. 1994, 49, 285-301. (10) Gokhale, V.; Hurowitz, S.; Riggs, J. B. A Comparison of Advanced Distillation Control Techniques for a Propylene/Propane Splitter. Ind. Eng. Chem. Res. 1995, 34, 4413-4419. (11) Hovd, M.; Michaelsen, R.; Montin, T. Model Predictive Control of a Crude Distillation Column. Model Identif. Control 1999, 20, 75-81. (12) Norquay, S. J.; Palazoglu, A.; Romagnoli, J. A. Application of Wiener Model Predictive Control (WMPC) to an Industrial C2Splitter. J. Process Control 1999, 9, 461-473. (13) Mathur, U.; Conroy, R. J. Successful Multivariable Control without Plant Tests. Hydrocarbon Process. 2003, 82, 55-65. (14) Shinskey, F. G. Feedback Controllers for the Process Industries; McGraw-Hill: New York, 1994. (15) Lundstro¨m, P.; Lee, J. H.; Morari, M.; Skogestad, S. Limitations of Dynamic Matrix Control. Comput. Chem. Eng. 1995, 19, 409-421.

(16) Morari, M.; Lee, J. H. Model Predictive Control: Past, Present and Future. Comput. Chem. Eng. 1999, 23, 667-682. (17) Camacho, E. F.; Bordons, C. Model Predictive Control; Advanced Textbooks in Control and Signal Processing; Springer: London, 1999. (18) Mayne, D. Q.; Rawlings, J. B.; Rao, C. V.; Scokaert, P. O. M. Constrained Model Predictive Control: Stability and Optimality. Automatica 2000, 36, 789-814. (19) Maciejowski, J. M. Predictive Control with Constraints; Prentice Hall: Upper Saddle River, NJ, 2002. (20) Smith, H. W.; Davison, E. J. Design of industrial regulators. Integral feedback and feedforward control. Proc. IEE 1972, 119, 1210-1216. (21) Pannocchia, G.; Rawlings, J. B. Disturbance Models for Offset-free Model Predictive Control. AIChE J. 2003, 49, 426437. (22) Muske, K. R.; Badgwell, T. A. Disturbance Modeling for Offset-Free Linear Model Predictive Control. J. Process Control 2002, 12, 617-632. (23) Pannocchia, G., Robust Disturbance Modeling for Model Predictive Control with Application to Multivariable Ill-Conditioned Processes. J. Process Control 2003, 13, 693-701. (24) Rajamani, M. R.; Rawlings, J. B.; Qin, J. Optimal Estimation in the Presence of Incorrect Disturbance Model; Technical Report 2004-09 (Draft); TWMCC, Department of Chemical Engineering, University of Wisconsin-Madison: Madison, WI, 2004. (25) Muske, K. R.; Rawlings, J. B. Model Predictive Control with Linear Models. AIChE J. 1993, 39, 262-287. (26) Rawlings, J. B.; Muske, K. R. Stability of Constrained Receding Horizon Control. IEEE Trans. Autom. Control 1993, 38, 1512-1516. (27) Chmielewski, D.; Manousiouthakis, V. On Constrained Infinite-Time Linear Quadratic Optimal Control. Syst. Control Lett. 1996, 29, 121-129. (28) Scokaert, P. O. M.; Rawlings, J. B. Constrained Linear Quadratic Regulation. IEEE Trans. Autom. Control 1998, 43, 1163-1169. (29) Rawlings, J. B. Tutorial: Model Predictive Control Technology. In American Control Conference; IEEE Press: Piscataway, NJ, 1999; pp 662-676. (30) Ljung, L. System Identification: Theory for the User, 2nd ed.; Prentice Hall Inc.: Upper Saddle River, NJ, 1999.

Received for review May 17, 2004 Revised manuscript received September 29, 2004 Accepted September 29, 2004 IE0495832