# Integration of design and control: a robust control approach using MPC.

INTRODUCTIONIn the traditional plant design approach the design of the process and the controller are performed in two separate sequential steps. Once the process has been designed based on capital and operating costs calculated at steady state, a controller is then synthesized to satisfy a robust performance criterion that is related to closed-loop stability and performance.

On the other hand, it has been widely recognized that better designs can be obtained by optimizing the system based on an objective function that takes into account simultaneously both the process steady state and dynamic behaviours (Seider et al., 1999). Due to this fact, researchers have resorted to integrate process design and process control into one single step. Several methodologies have been proposed to tackle this problem. These techniques can be classified into two main groups as follows. The first group consists of methodologies based on multi-objective optimization problems where the cost of control was quantified by using controllability indicators such as a disturbance cost (Lewin, 1996), the maximal eigenvalue (Blanco and Bandoni, 2003), the gap metric (Lee et al., 2000) or the condition number (Luyben and Floudas, 1994). Since these controllability indices are based on linear theory, linear models are often used by these methods to approximate the process transient behaviour. The key disadvantage in using controllability indices is that it is no trivial task to relate these indices to an actual economic cost and that they ignore the non-linear nature of the process. The second group of methodologies propose a single performance objective function applying a set of weighting functions and use the complete process dynamic model together with feedback regulatory structures to assess the integration of design and control (Mohideen et al., 1996; Bansal et al., 2000; Georgiadis et al., 2002). Complex and time consuming non-linear constrained dynamic optimization problems have resulted when using these approaches. Consequently, these approaches have solved the problem over pre-specified finite time horizons only. A more comprehensive review of the methods that perform the integration of design and control can be found in Sakizlis et al. (2004) and Seferlis and Georgiadis (2004).

This paper, therefore, addresses the following issues:

1. The integration of the cost of variability in the controlled variable during dynamic operation into one objective function together with the design cost given by the sum of capital and operating costs.

2. The development of a robust model used to represent the nonlinear model by a family of linear models. This family of models is represented by a nominal linear model supplemented with model parameter uncertainty. This representation has been previously used in robust control studies (e.g. Doyle et al., 1989). Thus, the uncertainty accounts for model nonlinearity with respect to an identified nominal linear model that serves as the internal model in the predictive controller. The worst-case variability for this family of models is then analytically quantified and its associated economic cost is calculated. This term will be referred to as the robust variability cost. This approach offers significant computational advantages since, by calculating analytical bounds for the variability, lengthy non-linear dynamic simulations are avoided. Also, by using the proposed analytical bounds, the variability can be estimated over an infinite time horizon.

3. The integration problem is solved using a Model Predictive Control algorithm. In previous studies dealing with integration of control and design, only multi-loop PI controllers were considered (e.g. Lee et al., 2000), whereas control strategies such as Model Predictive Control (MPC), especially suited for multivariable problems, have not been addressed. The calculation of closed-loop variability resulting from the application of MPC to a non-linear process represented by a nominal linear model with uncertainty is novel and it has not been previously reported in the literature.

In the following sections, the proposed integrated design approach is introduced and applied to the design of a distillation column. Finally, the traditional approach where design and control are performed sequentially is then compared to the proposed integrated control and design approach.

INTEGRATED APPROACH TO DESIGN AND CONTROL

Objective Function

The objective can be represented as follows:

min CC([d.sub.e],[u.sub.op]) + OC([d.sub.e], [u.sub.op]) + VC([d.sub.e], [u.sub.op], c) (1)

[d.sub.e], [u.sub.op], c s.t. h([d.sub.e], [u.sub.op], c)=0 g([d.sub.e], [u.sub.op], c) [less than or equal to] 0

where CC is the capital cost at steady state conditions, OC is the operating cost at steady state conditions, VC is the variability cost due to imperfect control, h(x) are equality constraints, e.g. mathematical relations obtained from energy and mass balances, g(x) are inequality constraints, e.g. constraints on controlled or manipulated variables or stability conditions, [d.sub.e] are design variables, e.g. geometric characteristics of the process to be designed, [u.sub.op] are operating variables, e.g. steady state values of flow rates or heat fluxes, and c are controller tuning parameters. In this work the only controller tuning parameter considered in optimization is, for simplicity, the input weight magnitude [lambda], which appears in the MPC controller gain matrix (see Appendix A, Equation (A-6)). Although, the control and prediction horizons m and p could also be used for tuning the controller, they were assumed constant in the case study for simplicity.

The calculation of the different terms in the above optimization problem is described in the following subsections. The capital and operating costs are process specific and they will be discussed later in the manuscript. The variability can be calculated in a systematic fashion by using a robust control approach as explained in the following subsection.

Output Variability under MPC Control

The process variability cost, VC([d.sub.e], [u.sub.op], c), can be calculated in general form based on the representation of the non-linear model as a nominal model with uncertainty and the MPC algorithm as follows.

The non-linearity of the process is captured by the so-called model uncertainty. The model uncertainty, [l.sub.j,k,i], can be described as shown in Figure 1 and is defined later in the section.

Figure 1 shows qualitatively the step response to steps of different magnitudes around a specific operating condition. Also, for non-linear systems, different responses will be obtained for step changes of equal magnitude but different signs. Given a nominal step response model, a robust model can be defined with respect to that nominal model by using the concept of uncertainties as follows. For each step response coefficient, the non-linear nature of the process, resulting in these different step responses, is quantified by the multiplicative model uncertainty, [l.sub.j,k,i], as follows:

[l.sub.j,k,i] = [S.sub.j,k,i] - [S.sup.nom.sub.j,k,i]/[S.sup.nom.sub.j,k,i] (2)

where [S.sub.j,k,i] and [S.sup.nom.sub.j,k,i] are the actual and nominal step response coefficients for each combination of input k and output j and at each time interval i, respectively. The largest model uncertainty to be referred to as [l.sup.max.sub.j,k,i], is the largest deviation between the nominal and upper or lower bound models shown in Figure 1 defined as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (3)

and the nominal model is defined as follows:

[S.sup.nom.sub.j,k,i] = [S.sup.up.sub.j,k,i] + [S.sup.lw.sub.j,k,i]/2 (4)

[FIGURE 1 OMITTED]

where [S.sup.up.sub.j,k,i] and [S.sup.lw.sub.j,k,i] are the upper and lower bounds on each one of the coefficients [S.sub.j,k,i]. Based on above description, a robust model, [S.sub.n], can be defined to describe the process non-linear dynamic behaviour as a function of a nominal linear model [S.sup.nom.sub.n] and the uncertainty related to each of the nominal model parameters. The robust model [S.sub.n] has elements [S.sub.j,k,i] such that:

[S.sub.j,k,i] = [DELTA][S.sup.nom.sub.j,k,i](1 + [l.sub.j,k,i] (5)

where [S.sup.nom.sub.j,k,i] denotes the elements of nominal linear model [S.sup.nom.sub.n] and [l.sub.j,k,i] was given in (2).

With regards to the MPC algorithm, the reader is referred to Appendix A where the algorithm and corresponding notation are presented. In this work, it is assumed that the measured disturbances and set-point vectors are zero, i.e., [S.sub.d][DELTA]d(i-1)=0 and R(i + 1) = 0. If only rejection of unmeasured disturbances is to be considered within the MPC control algorithm, the error at time i will be as follows:

[epsilon](i) = -[M.sub.p]Y(i - 1) - w(i) (6)

where Y(i) is the output prediction vector, [M.sub.p] is the truncated shift matrix within the MPC algorithm and the symbol w denotes the unmeasured disturbance vector at time step i. Substituting the unmeasured disturbance vector (A-4) and (6) into the updated model (A-2), Y(i) is expressed as follows:

Y(i) = ([M.sub.n] - [S.sub.n][K.sub.MPC][M.sub.p])Y(i - 1) - [S.sub.n][K.sub.MPC]w(i) (7)

Rearranging Equation (7) and after applying Z-transform, the transfer matrix between prediction to disturbances can be expressed as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (8)

The unmeasured disturbance is assumed to be sinusoidal for the calculations. Then, the discrete frequency response can be used to compute a bound over the maximal value of the variability in the prediction output vector Y(i) after neglecting the initial transients by setting Z = [e.sup.lwT] where [omega] is the frequency of the unmeasured disturbance w and T is the sampling period. Accordingly, the following transfer matrix can be defined:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (9)

Each of the elements of the transfer matrix G([e.sup.twT]) is a transfer matrix [G.sup.Y.sub.j,k,i] that describes the prediction step i of the response of output k to a change in the unmeasured disturbance k. Since in this work the disturbance is assumed to be known a priori, e.g. a sine of a specific frequency, it is possible to obtain an upper bound over the "actual" expected output. Since the disturbance is also assumed to be periodic, any row matrix corresponding to a specific i in Equation (9) can be used to calculate the maximum output variability of the frequency response. For example, setting i=1 in Equation (9) and using the fact that Y(i) = w + Y(i) and Y = Gw, the robust variability or the largest variability for the family of linear plants representing the non-linear system is calculated from the frequency response using Equations (8) and (9) as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (10)

where [Y.sup.max.sub.k,l] is a scalar value that specifies a bound over the largest variability of the output k at the time step i=1. The symbol [G.sup.Y.sub.j,k,i] denotes the transfer matrix in the first ny rows of the transfer matrix in Equation (9) and W denotes the amplitude of the disturbance. The sum in the right-hand side of Equation (10) is due to the fact that the prediction at each step i is given by the multiplication of a vector of transfer functions by the disturbance vector w of length p, i.e., the prediction horizon as defined in Appendix A. The solution to Equation (10) was found with the function fminsearch in Matlab and different initial guesses were tried to ensure that a global optimum is approached.

If sinusoidal disturbances of different frequencies are considered, the maximization in (10) has to be also performed with respect to frequency since the maximum variability over time is bounded by the worst variability as a function of frequency (Morari and Zafiriou, 1989, p. 232) as per the following equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (11)

Then, the above relation can be used to calculate a bound over the standard deviation of the output. Clearly, for the single sinusoidal case the ratio of the norms of W and Y is equal to the ratio between the corresponding amplitudes resulting in Equation (10). It should also be noted that the accuracy of this calculation depends on the length of the control horizon p, which also determines the dimension of the disturbance w (see Equation (A-4)). Then, for better accuracy p has to be equal or larger than the settling time of the closed-loop. It should be noted that the variability could have been calculated from dynamic simulations of the non-linear process. However, the key advantage of calculating analytically a bound on the variability based on a nominal model with uncertainty as per Equation (10) is that long dynamic simulations are avoided and this is essential in order to reduce the duration of the optimization. One obvious disadvantage of using uncertain models is that some conservatism may result. One possibility to further reduce the variability is to use a non-linear MPC approach. However, this is beyond the scope of the present study.

Manipulated Variable Constraint

The same idea used to compute the maximum output magnitude of the frequency response can be used to compute the maximum variability in the manipulated variables during closed-loop operation. Based on the MPC algorithm, the following transfer matrix in the z-domain between the manipulated variable moves and the unmeasured disturbance can be obtained:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (12)

Using as before, frequency response concepts, the maximum value of the manipulated variable move, [DELTA][u.sup.max.sub.k,i] can be calculated from the transfer matrix in Equation (12) as follows (i=1):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (13)

where [DELTA][u.sup.max.sub.k,l] is a scalar value that specifies a bound over the largest variability of the manipulated variable k at the time step i=1. The symbol [G.sup.U.sub.d,k,1] represents the transfer matrix in the first [n.sub.u] rows of the transfer matrix in Equation (12). Then, the largest manipulated variable change is bounded as follows (i=1):

[u.sub.k,nom] + [DELTA][u.sup.max.sub.k,l] [less than or equal to] [u.sub.k,max] (14)

Equation (14) was incorporated within the proposed optimization problem as a constraint.

Robust Stability Constraint

The robust stability condition, developed by Zanovello and Budman (1999), in terms of a Structured Singular Value condition referred to as [[micro].sub.[DELTA]] analysis (Morari and Zafiriou, 1989) is expressed as:

[[micro].sub.[DELTA]](M([omega])) < 1 (15)

where M([omega]) is the interconnection matrix defined as follows:

[u (i + 1)/[y.sub.[DELTA]]] = M [u(i) [u.sub.[DELTA]] (15a)

[u.sub.[DELTA]] = [[DELTA].sub.u][y.sub.[[DELTA]] (15b)

and is a function of the MPC controller and the nominal linear step response model, [S.sup.nom.sub.n]. In Equation (15a) u(i) represents the vector of manipulated variables as defined by Equation (A-19). The variables [u.sub.[DELTA]] and [y.sub.[DELTA]] are dummy variables that can be eliminated by combining Equations (15a) and (15b) to result in the closed-loop system equation:

u (i + 1) = F[u(i), [[DELTA].sub.u]] (15c)

where the operator F can be obtained by combining the controller equations given in Appendix A together with the robust model equation presented above.

The quantity [[micro].sub.[[DELTA] is defined as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (16)

The interconnection matrix is generated applying linear fractional transformations (Zhou and Doyle, 1997); for details regarding the construction of the matrix M the reader is urged to refer to the paper by Zanovello and Budman (1999). The uncertainty block [[DELTA].sub.u] is a diagonal matrix with scalar elements as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (17)

The elements of the above matrix are related to the model parameter uncertainties [l.sub.j,k,i] given by Equation (2) as follows:

[l.sub.j,k,i] = [l.sup.max.sub.j,k,i][[DELTA].sub.l]

Each [[DELTA].sub.l] in (17) represents a perturbation that is bounded by the following condition (Zhou and Doyle, 1997):

[sigma]([DELTA]) [less than or equal to] 1 (18)

where [sigma] denotes the maximum singular values of the perturbation matrix [DELTA]. The Structured Singular Value problem can be easily calculated using the Matlab Robust Control Toolbox.

CASE STUDY

To test the methodology proposed in previous section, a depropanizer process (Lee, 1993) is used as a case study where the mole fractions of propane in distillate and isobutane in bottom product are controlled by simultaneously adjusting the reflux ratio (RR) and boil-up rate (Q) in the face of sinusoidal disturbances in feed composition. Thus, the controlled variables are the mole fraction of propane and isobutane in distillate and bottom products, respectively, and the manipulated variables are the reflux ratio RR and boil-up rate Q. The mole fractions of propane in distillate and isobutane in bottom product were specified to be 0.783 and 0.1, respectively, whereas the manipulated variables, i.e., the reflux ratio (RR) and boil-up rate (Q), were varied accordingly to meet these product targets. The distillate flow rate and bottom flow rate were allowed to vary to leave enough degrees of freedom for the steady state calculations. The process model was simulated with various distillation column designs (19-59 stages) using the rigorous equilibrium stage model, Radfrac, in Aspen Plus. In order to carry out the optimization problem described in Steady-State Correlations section, two types of correlations are needed:

1. Steady-state model correlations introduced as equality constraints in the optimization problem; that is, h([d.sub.e], [u.sub.op], c)=0 in Equation (1).

2. Dynamic model correlations to calculate the output variability in the cost function and the manipulated variable and robust stability inequality constraints described above; that is, g([d.sub.e], [u.sub.op], c) [less than or equal to] 0 in Equation (1).

These two types of correlations are discussed in the following subsections. Then, the model uncertainty can be calculated and the variability in the output variables can be computed from the nominal models and the uncertainty around different operating conditions.

Steady-State Correlations

The steady-state solution was obtained from the simulation of the case study at different conditions using Aspen Plus. In principle the steady state values of manipulated and controlled variables could be solved along the optimization problem for any set of given conditions. Instead, to simplify the optimization procedure, results were obtained for a wide range of conditions and empirical correlations were obtained between the variables in the form of polynomials calculated from simple least squares fitting, the structure of the polynomials used in this work to represent the transients in Q and RR were obtained using the software TableCurve2D[R] (TableCurve2D, 2002) based on steady state results obtained with ASPEN. To test the accuracy of these correlations they were verified by using additional values that were not used for model calibration. For example, for the target top and bottom composition given above, the manipulated variables, i.e. the reflux (RR) and the boil-up rate (Q) are effectively correlated according to the following empirical relation calculated from regression of ASPEN results of a large set of values:

Q = 209.58 x [RR.sup.2] - 178.01 x RR + 3548.40 (19)

The values of Q and RR that meet this relation provide the required mole fraction of propane in distillate of 0.783 and the mole fraction of isobutane in bottom product of 0.1 at steady-state conditions. For the dynamic situations these two variables are not directly correlated to each other. Also, due to the steady-state correlation between Q and RR, the design variables at steady state can be obtained as a function of RR only. Thus, the number of stages N, column diameter D (m), heat duty in reboiler HD (GJ/h), distillate flow rate DF (lbmol/h) and bottom flow rate BF (lbmol/h) as functions of the reflux ratio RR were all determined using Radfrac. After applying least squares regression to the results obtained from ASPEN at different RR values the following empirical correlations were obtained:

N = -156.06 + 21.58/ln(RR) + 161.35/[(ln(RR)).sup.2] - 242.13/[(ln(RR)).sup.3] + 72.91/[(ln(RR)).sup.4] (20)

D = 0.19 x [RR.sup.2] - 0.56 x RR + 3.76 (21)

HD = 1.09 x [RR.sup.2] + 1.38 x RR + 25.49 (22)

DF = 42.92 x [RR.sup.2] - 203.96 x RR + 1210.35 (23)

BF = -44.21 x [RR.sup.2] + 209.27 x RR + 2013.48 (24)

The steady-state correlations given by Equations (19) to (24) were incorporated into the optimization problem as equality constraints. Although the two manipulated variables in this case study, i.e. RR and Q are correlated at steady state, for the closed-loop process under MPC control the values of RR and Q are independent during dynamic operation.

Dynamic Modelling Correlations

The dynamic models are formulated to represent the dynamic behaviour of the input and output variables in terms of deviation values defined with respect to different steady state values that are calculated from the correlations shown in the previous section. Although step response models are needed for MPC control, in order to simplify the identification and the optimization it was decided to parameterize the models by using a first order plus dead time (FOPDT) approximation as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (25)

In that way the dynamic models could be simply represented by three variables, i.e., a gain [K.sub.p], a time constant [[tau].sub.p] and a time delay [theta]. The gains were identified from Aspen Plus calculations whereas the time constants and time delays were obtained from empirical design correlations proposed by Seider et al. (1999). Thus, an additional advantage in using the FOPDT approximation is that empirical design correlations can be used. Although one could identify the time constants and time delays directly from Aspen Plus dynamic simulations, Seider's design correlations were preferred to simplify the identification process. These empirical correlations, summarized in Appendix B relate the dead time and the time constants to the liquid flow rate in the trays, the number of trays and the column diameter. Since all these three variables are functions of RR as per Equations (20) and (21), both the dead times and the time constants for all the dynamic models can be expressed as a function of RR as shown later in this section.

The FOPDT assumption was also justified from the typical delayed exponential response to steps obtained for all variables with the Aspen Plus dynamic simulator of the distillation column. For example, Figures 2 and 3 show a series of responses of the top and bottom compositions to step changes in Q and RR of different magnitudes. It is also clear from these responses, that the process is highly non-linear since the gains varied as a function of the sign and the magnitude of the step change in Q.

The process gains were calculated from the Aspen Plus responses to step changes of [+ or -] 1% and [+ or -] 35% of the nominal RR value and [+ or -] 1% and [+ or -] 40% of nominal Q value used to ensure that a wide range of conditions is covered around the nominal values.

Both the responses to positive and negative changes are taken into account for the calculation of the nominal gain, [K.sub.p,RR]. For example, the nominal gain in terms of the average gains obtained for [+ or -] 1% and [+ or -] 35% changes is as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (26)

where [DELTA]y represents the steady state changes to the step changes in RR of the top or bottom compositions. It was found that all other step changes in between [DELTA]RR = [+ or -] 1% to -RR = [+ or -] 35% result in nominal gains values in between the extreme values. Thus, the model uncertainty in gain in combination with the nominal gain in Equation (26) can be used to describe all the possible [K.sub.p,RR] values. The same calculation was used to calculate the process gains of the compositions with respect to the boil-up Q.

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

For the MIMO case with two controlled variables ([x.sub.D], [x.sub.B]) and two manipulated variables (RR,Q), four gains are obtained for the four possible pairings between controlled and manipulated variables: [x.sub.D]-RR, [x.sub.D]-Q, [x.sub.B]-RR and [x.sub.B]-Q. The experiments shown in Figures 2 and 3 were performed around different steady state values of RR and their corresponding Q values according to Equation (19). Then the set of gain values obtained around the different steady states were regressed with respect to the values of RR and Q resulting in the following equations:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (27)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (28)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (29)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (30)

where a normalized boil-up rate is defined as Q' = Q/2000.

Figure 4 shows the plots of the process gain relating changes in RR to the top composition showing a non-linear behaviour at low RR values and a more linear dependency for higher RR values. Although not shown in the paper for brevity, the other gains show similar behaviour to the one in Figure 4.

For the calculation of the process time constants ([[tau].sub.p]), the formula presented by Seider et al. (1999) was used to obtain the process time constant. As shown in Appendix B, the time constants are a function of the liquid flow rate in each tray. Steady state calculations for step changes in RR and Q can be performed using Aspen Plus, and from these, the liquid flow rate in the trays can be obtained around different steady state conditions. Then, the process time constants based on the Seider's formula involving the liquid flow rate on each tray in the distillation column can be calculated. For each step change in RR ([+ or -] 1% and [+ or -] 35%) and Q ([+ or -] 1% and [+ or -] 40), the set of the time constants were calculated and the corresponding average time constant was calculated for every different column design. For example, the following expression was used to calculate the process time constant, op for each column design:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (31)

Following Seider's formula where the time constants are solely functions of the values of RR and Q, there are only two expressions of the time constants for the four possible pairings. Thus, the time constant as function of RR can be used for both parings [x.sub.D]-RR and [x.sub.B]-RR whereas the time constant as function of Q was used for pairings [x.sub.D]-Q, [x.sub.B]-Q. Then, the process time constants were represented by functions of the manipulated variables RR and Q by using least squares regression as follows:

[[tau].sub.p1] = 4.91/[(log(RR)).sup.4] - 60.59/[(log(RR)).sup.3] + 276.40/[(log(RR)).sup.2] - 548.73/(log(RR)) + 404.97 (32)

[[tau].sub.p2] = 0.11/[(log(Q')).sup.4] - 2.23/[(log(Q')).sup.3] + 13.95/[(log(Q')).sup.2] - 23.82/(log(Q')) - 14.56 (33)

[FIGURE 4 OMITTED]

The empirical correlation shown in Equation (32) was plotted in Figure 5 showing similar behaviour to the process gain described above, i.e., non-linear dependency at low values of RR versus linear dependency at high values of RR.

The process dead time, also determined using the formula presented in Seider et al. (1999), mainly depends, as shown in Appendix B, on the liquid flow rate on each tray. For the distillation column, the liquid flow rate on each tray is directly related to the reflux ratio RR. Then, the process dead time of different conditions can be formulated as function of RR as follows:

[theta] = 1.23 x RR + 2.61 (34)

The process model as represented by Equation (25) can be now formulated in terms of the step response model, [S.sub.n] for use in the MPC algorithm. Each one of the step response coefficients is a function of the reflux ratio RR and the boil-up rate Q because each one of the FOPDT model parameters, i.e., [K.sub.p], [[tau].sub.p], [theta] are functions of RR and Q.

Model Uncertainty

The uncertainty was estimated by considering positive and negative steps in the manipulated variables, RR and Q, around the nominal operating point with the design variables fixed at a specific value. The differences in the responses obtained for negative and positive changes are due to the system non-linearity. The uncertainty to be implemented into the integrated design and control problem is considered as the variation in the step response model with respect to the nominal values identified in the previous section. The bound on the uncertainty was defined as the maximum uncertainty, [l.sub.max], that can take either positive or negative values and is defined later in this section. The negative value and the positive value represent the lower and upper bounds for uncertainty, respectively.

Since MPC is based on step response models, it was a natural choice to represent the model parameter uncertainty as the maximum variation of the actual step response from the nominal one calculated using the nominal parameters from the previous section. Uncertainty in the process dead time was not considered because the dead time takes effect immediately after the input change and therefore, the dead time is mostly determined by the nominal steady state around which the step change occurs. Thus, the model parameter uncertainty is in principle related to the variation in process gain and time constant when different step sizes of manipulated variable were applied. However, based on extensive Aspen Plus dynamic simulations, it was observed that the effects of changes in the time constant due to the shift from one steady state to the next were found negligible as compared to the changes in the process gain. Thus, the uncertainty in process gain was found to produce the maximal deviation of the step response with respect to the nominal response. Figure 6 shows qualitatively the type of responses obtained in simulations and how the uncertainty in the process gain can be used to capture the responses to step changes. This figure also shows that the maximal deviation with respect to the nominal model is expected to occur at steady state.

[FIGURE 5 OMITTED]

Following the arguments above and using Equations (2) to (4), the maximum uncertainty in terms of the actual and nominal step response coefficients in the process gain was obtained as follows:

[l.sub.max] = [max.sub.i][l.sup.max.sub.j,k,i] (35)

where the elements [S.sup.pu.j,k,i], [S.sup.nom.sub.n] and [S.sup.lw.sub.j,k,i] required to estimate [l.sub.max] in (35) were obtained from the discretization of the dynamic response of the models [y.sub.upper], [y.sub.nom] and [y.sub.lower] shown in Figure 6, respectively.

In order to simplify the optimization procedure, [l.sub.j,k,i] was formulated as explicit functions of the nominal values of RR and Q. Then the expressions for the model uncertainty of each pairing of the controlled and manipulated variables were found to be as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (36)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (37)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (38)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (39)

For example, Equation (36) was plotted in Figure 7 showing that at low RR values, the uncertainty in the parameters of the model [x.sub.D]-RR is high, that is, the dynamic behaviour [x.sub.D]-RR is non-linear whereas at high RR values, the uncertain values are small, which means that [x.sub.D]-RR behaves in a linear fashion.

The model parameter uncertainties that describes the process non-linearities between the manipulated and controlled variables are lumped in an uncertainty set, L, defined as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (40)

[FIGURE 6 OMITTED]

The uncertainties in the set L are varied in order to obtain the maximum variability cost, the manipulated variable constraint and the robust stability limit as explained in the Integrated Approach to Design and Control section.

Cost of Process Variability

A sinusoidal disturbance with frequency [[omega].sub.d], was introduced as the disturbance in the feed composition resulting in sinusoidal responses in the distillate and bottom compositions. Thus, these concentrations oscillate, due to imperfect control, around their pre-specified set points when under closed-loop control. To quantify the cost of this variability in the controlled variables, it was assumed, following common industrial practice, that fixed volumes of products from the distillation column will be held downstream from the distillation column in two intermediate storage tanks, one for the top product and one for the bottom product, before delivery to the customer. Each tank will also have the effect of attenuating the variation in product composition. A large volume of material in each storage tank will naturally result in a lower variability in the exit of the storage tank but, at the same time, it will increase the inventory cost. The variability cost was then assumed to be equal to the inventory cost that is directly related to the volume of product stored in these intermediate storage tanks.

The volume of intermediate storage was calculated at the specific frequency of the disturbance, [[omega].sub.d]. From a simple mass balance of each one of the intermediate storage tank,

V = d[C.sub.out]/dt = [Q.sub.in][C.sub.in] - [Q.sub.out][C.sub.out] (41)

Assuming constant flow rates Q = [Q.sub.in] = [Q.sub.out] and defining the residence time [tau] = V/Q where V is the product volume, from Laplace transformation of Equation (40):

[[delta]C.sub.out]/[[delta]C.sub.in] = 1/[tau]S + 1 (42)

where [[delta]C.sub.in] is equal to the maximal variability of the products exiting the distillation column and entering the storage tanks calculated from Equation (10) and to be referred to as [x.sub.D,max] or [x.sub.B,max] for the top and bottom concentrations, respectively. Then converting to the frequency domain at [omega] = [[omega].sub.d]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (43)

where Equation (43) represents the volume of each product in its intermediate storage tank, V, as a function of [[delta]C.sub.in], the maximal deviations of the product concentrations [x.sub.D,max] or [x.sub.B,max], [[delta]C.sub.out], the maximum allowed product concentration variability to be delivered to the customer, Q, the product flow rate and [[omega].sub.d], the frequency of the disturbance. Therefore, the volume of product, V, in the intermediate storage tank depends only on [[delta]C.sub.in], i.e., the top or bottom maximal variability in concentration of the products exiting the distillation column, because [[delta]C.sub.out] and Q are specified by the customer and [[omega].sub.d] is the frequency of the unmeasured disturbance. In this work, [[delta]C.sub.out] was set equal to 1% of the set points for the top and bottom compositions, respectively. Then, the variability cost (VC) is calculated by considering the price of the product held in the intermediate storage tank before retailing as follows:

VC = P x V x (A /P, i, N) (44)

[FIGURE 7 OMITTED]

where V is the volume calculated from Equation (43), P is the price of the product and (A/P,i,N) is the capital recovery factor with interest rate, i, and amortization period N. The volume of the distillate product to be held in the holding tank was calculated using Equation (43). The volume V1 depends on the maximal variability in mole fraction of propane in distillate ([x.sub.D,max]) and the distillate flow rate (DF).

V1 [congruent to] [[delta]C1.sub.in]/[[delta]C1.sub.out] x DF/[[omega].sub.d] (45)

The volume V2 depends on the variability in isobutane mole fraction in bottom product ([x.sub.B,max]) and the bottom product flow rate (BF).

V2 [congruent to] [[delta]C2.sub.in]/[[delta]C2.sub.out] x BF/[[omega].sub.d] (46)

Similarly, the variability cost can be defined as the one shown in Equation (44). The following expression is the variability cost of the distillate product.

VC1 = W1 x P1 x V1 x (A/P,i,N) (47)

The following equation states the variability cost of the bottom product.

VC2 = W2 x P2 x V2 x (A/P,i,N) (48)

where W1 and W2 are the price correction factors to present variations in prices of propane and isobutane, respectively, and P1 and P2 are the nominal price of propane and isobutane. (A/P, i, N) is the capital recovery factor to convert the principal cost of the products to an annualized cost, and in this case it equals to 0.115. Both costs, i.e., VC1 and VC2 were added into the general cost function together with the capital and operating costs of the distillation column under study.

Manipulated Variable Constraint

The limits on the manipulated variable moves calculated at the disturbance's frequency were considered by bounding the elements in the first and second rows of the transfer matrix [G.sup.U.sub.d,k,l] in Equation (13). These rows represent the dynamic moves in RR and Q, respectively. Then for the disturbance rejection problem, the dynamic control moves of RR and Q can be calculated as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (49)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (50)

where T is the sampling period, [omega] is the frequency of the disturbance and W is the amplitude of the disturbance. In the optimization problem, the values of [RR.sub.max] and [Q.sub.max] were fixed, at [RR.sub.max] = 2.63 and [Q.sub.max] = 3980 lbmole/h. For example, the operating range of RR is between 1.69 and 2.63; the values of RR from 1.69 to 2.63 correspond to the nominal designs of the column 59 to 19 stages, respectively. Therefore, the constraints for RR and Q during the optimization can be as follows:

[RR.sub.nom] + [DELTA]RR [less than or equal to] 2.63 (51)

[Q.sub.nom] + [DELTA]Q [less than or equal to] 3980 (52)

Equations (50) and (51) were incorporated in the optimization problem as the manipulated constraints.

OPTIMIZATION PROBLEMS USING THE INTEGRATED METHODOLOGY VERSUS THE TRADITIONAL METHODOLOGY

This section presents the optimization problems for the two design methodologies considered here, i.e., the proposed integrated design approach and the traditional two-step design approach described in the introduction.

Integrated Design Methodology

This method consists of only one optimization problem as follows :

[min.sub.RR,Q,[lambda]] CC(RR) + OC(RR) + [max.sub.L] VC(RR,Q,L,[lambda] (53) s.t.

h(RR,Q) = 0 g(RR,Q,[lambda]) [less than or equal to] 0

where [lambda] is the controller tuning parameter that appears in the MPC gain matrix calculation shown in Appendix A. The equality constraints correspond to correlations shown in the Steady-State Correlations section and the inequality constraints correspond to the manipulated variable and robust stability constraints as stated in Equations (50), (51) and (15), respectively. The capital cost, CC, and operating costs, OC, were calculated using the following correlations proposed by Luyben and Floudas (1994) for a distillation column:

CC = 12.3(615 + 324[D.sup.2] + 486(6 + 0.76N)D) + 245N(0.7 + 1.5[D.sup.2]) (54)

OC = [[beta].sub.tax] ([Q.sub.r] x OP x UC (55)

where [Q.sub.r], OP and UC are reboiler heat duty, operating period and utility cost. All these variables can be related to the reflux ratio RR as per equations presented in the Steady-State Correlations section, thus simplifying considerably the optimization process.

Traditional Design Methodology

This method consists of two optimization steps performed sequentially as follows:

Step 1. Steady state optimal design

min CC(RR) + OC(RR) (56) RR s.t.

h(RR) = 0 g(RR) [less than or equal to] 0

The equality constraints are the correlations for number of stages, diameter and heat duty calculations given by Equations (20) to (24). The inequality constraint represents the constraint on the manipulated variable, RR, as represented by the following expression:

[RR.sub.min] = RR = [RR.sub.max] (57)

Step 2. Optimal controller design

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (58)

where [beta] is a scalar and the robust performance condition, [[micro].sub.[DELTA]]([M'.sub.perf]([beta]) < 1) was obtained from Morari and Zafiriou (1989) and it follows the [micro]-synthesis design methodology. The [[micro].sub.[DELTA]](M'.sub.perf]([beta])) problem is schematically presented in Figure 8. The perturbation matrix [DELTA] is defined as [DELTA] = diag{[[DELTA].sub.p], [[DELTA].sub.u]} whereas [[micro].sub.[DELTA]](M'.sub.perf]([beta])) is the interconnection matrix that is defined as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (59)

where [M.sub.perf] is obtained from the following relationships:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (60)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (61)

The symbol [[DELTA].sub.u] is a function of the model uncertainty set L whereas [[DELTA].sub.p] relates the feedback errors in controlled variables to the input disturbances, i.e. disturbances in feed composition. [[DELTA].sub.p] is a full matrix of perturbations with number of rows as the number of errors and number of columns as the number of disturbances. A detailed description of the construction of the matrix [M.sub.perf] is given by Chawankul et al. (2005).

RESULTS AND DISCUSSION

All the results in the following case study were obtained for a sampling period of T=5 min. For simplicity it was assumed that the prediction step n=10, the prediction horizon p=3 and the control horizon m=1. Thus, for the current study, the only controller tuning parameter to be optimized is the input weight scalar [lambda]. A sinusoidal disturbance in feed composition of amplitude of 0.0215 and frequency of 0.0004 rad/s was considered. The amplitude of the disturbance was chosen to be of 2% propane in the feed, which is equivalent to a change of 0.0215 in mole fraction of propane in the feed. In principle, the frequency at which the variability is the largest could be found from the proposed optimization problem in combination with a one-dimensional search with respect to frequency; however, this has not been done at this point to avoid complexity in the calculations. For both the traditional method and the integrated method, the optimization was carried out for different values of the price correction factors W1 and W2 in order to assess the sensitivity of the results to product prices changes.

[FIGURE 8 OMITTED]

Traditional Design Method

From Step 1, i.e., the steady-state optimization, it was found that the optimal RR was 1.918. From Step 2, i.e., the control design step, the optimum tuning variable [lambda] was obtained at 0.505 and the performance weight [beta] was found to be 0.587.

The optimum design obtained from the first step in terms of RR and the tuning parameter [lambda] obtained from the second step in the traditional method were used to calculate the variability and the total costs. In Table 1, the sensitivity of the total cost with respect to the product price was conducted.

As expected, RR and [lambda] do not depend on the product cost because they are obtained from a steady-state optimization of capital and operating costs that are independent of the product cost. On the other hand, the variability cost depends on both RR and [lambda]. It is clear from the table that the variability cost is more sensitive to the cost changes in isobutane than to the cost changes in propane.

Integrated Method

The optimization results obtained by using the integrated method are shown in Table 2. This table shows that as the product prices, W1 and W2, increase the optimum RR is smaller. This results, following Equations (20) and (21), in a larger number of stages and a smaller column diameter. Consequently, the capital cost increases whereas the operating cost decreases.

The variability costs with respect to the increases in price of propane (W1) were found larger than the ones with respect to the price of isobutane (W2) increases. This can be explained in terms of the effect of the tuning parameter [lambda]; the values of [lambda] are larger (0.2332-0.2350) in the case that the price of propane is increased whereas for the isobutane price increases the values of [lambda] were found to be smaller (0.1830-0.1886). Smaller [lambda] results in more aggressive control and consequently, smaller variability.

To investigate the optimal design and cost without any automatic control system and only large intermediate storage tanks to attenuate the variability in the product concentrations the optimization was run with a large fixed [lambda]. By fixing [lambda] to a large value, 1000 in this case, we are essentially eliminating any control action by the controller and mimicking the case of operation without any controller. Table 3, shows the results for the case when the tuning parameter [lambda] was set at 1000.

By fixing [lambda], the objective function, the summation of capital, operating and variability costs, became only a function of RR, rather than RR and [lambda] as in the previous optimizations. The results show an RR of 1.7014 that corresponds to 51 stages and a 3.365 metre diameter column. Thus, the optimal design is no longer a function of the prices of the products. The design without control is higher in capital cost because the tower contains more stages (51 stages compared to about 30 for the integrated approach), lower in operating costs because the reflux ratio is lower (1.7014 compared to about 1.9 for the integrated approach with variable [lambda]) but much higher in variability costs because the design without control requires very large tanks to attenuate the product variability.

The variability costs, VC, in the integrated approach are typically less than $100/day whereas in the case without control they are all greater than $600/day and depending on the price of the products they can be as high as $7372/day. One may use the difference in variability cost with and without control as an estimate of the benefit of control or the cost of not implementing process control. Also, the difference in the values of ??between the two approaches (Tables 1 and 2) clarifies the necessity to take into account the process dynamics and control at the design stage. In fact, the cost savings described above come from the reduction in the process variability, which is due to simultaneous calculation of the process design variables and the MPC tuning parameter ?. As shown in Table 4, for this case, the difference in the variability cost with and without control ranges from $598.60/day to $7276.01/day clearly demonstrating the benefit of feedback control for this problem.

VC = VC1 + VC2

where VC1 and VC2 are the variability costs of the distillate and isobutene, respectively.

The process interaction is also studied and the results are shown in Table 5. The RGA element [[lambda].sub.11] representing the system interaction was calculated for all the designs. It was found that as the product price increases the optimal RR decreases towards the value of one which corresponds to zero interaction. The process interaction is directly related to the design of the process and it has a direct effect on the robust stability, i.e. the larger interaction the more difficult to satisfy the robust stability condition, (Morari and Zafiriou, 1989). Thus, in order to satisfy the robust stability condition, the optimization searches for smaller RR values that correspond to smaller interaction as shown in Table 5. At the optimal points the robust stability constraint is found to be always active. If RR at the optimum will be larger, the corresponding [[lambda].sub.11] will be larger and the robust stability constraint would be violated.

Thus, the optimum designs that are presented in Table 5 were selected with less interaction as the prices of the products increased, i.e., the less interaction, the easier to satisfy the robust stability. As the prices were increased the optimum tuning parameter [lambda] slightly decreases so the robust stability condition will be satisfied. Although the smaller [lambda] can cause system instability, a design with a relatively small degree of interaction can be still stabilized with a relatively small value of [lambda], i.e., more aggressive control. It should be noticed that the traditional approach, since it does not consider controllability issues at the design stage, resulted in higher RR values that correspond to larger interaction according to RGA calculations. When the interactions are larger the controller has to be tuned less aggressively in order to satisfy the robust stability test. This explains why the traditional approach results in higher variability costs and therefore, higher total costs as compared to the integrated method.

Finally, to illustrate the points above, in Table 6 the relative savings obtained by using the design obtained with integrated method versus the traditional method were compared. The total costs obtained from both traditional and integrated methods were compared. As shown, the total cost savings were significant when the product prices were increased.

From Table 6, as the propane price factor W1 was increased by 5, 10, 15 and 20 it was found that the savings were in the range of 12% to 28%. On the other hand, when the isobutane price factor W2 was increased by the same factors the savings were found to be about 26% to 94%. Thus, the sensitivity of the savings to the price changes of isobutane is more significant than the sensitivity to the price changes of propane.

CONCLUSIONS

A new design methodology that integrates the cost of variability in the controlled variables with the capital and operating costs into one objective function is proposed. The non-linear process was modelled by a family of linear models represented by a nominal FOPDT model and model uncertainty to account for non-linearity with respect to the nominal models around different values of reflux ratio. Empirical correlations were used within the optimization to represent steady state relations identified from ASPEN simulations. An MPC controller tuned for robustness in the presence of model uncertainty was used. The proposed integrated design methodology was compared with the traditional design methodology where the process design and control design are performed sequentially.

The proposed integrated design approach resulted in lower costs because the design using the integrated method allowed for a more aggressive tuning of the controller without violating the robust stability constraint. This method allows one to estimate the benefit of using process controllers versus intermediate storage tanks to attenuate variability in the output; depending on the price of the product the savings from using process controllers varied from 25% to 75%. The differences between the integrated method and the traditional one become more significant as the price of the product increases and were found to be as great as 95%.

APPENDIX A:

MPC ALGORITHM AND DEFINITIONS

For a better comprehension, the matrices and vectors used in the algorithm are defined at the end of the appendix. The MPC algorithm consists of the following steps:

1. Initialize the model prediction vector, Y at time i = 0:

Y(0) = [[y[(0).sup.T] , ... , [y(0).sup.T].sup.T](A-1)

It is assumed that the system is at steady state, that is, y(0)=0, where y(i) is the measured output vector at the current interval i.

2. Update the model

Y(1) = [M.sub.n]Y(i - 1) + [S.sup.nom.sub.n][DELTA]u(i - 1) + [S.sup.nom.sub.d][DELTA]d(i - 1) (A-2)

where [M.sub.n] is the shift matrix vector, [S.sup.nom.sub.n] denotes the nominal step response model vector, [S.sup.nom.sub.d] is the nominal step response model vector due to the measured disturbances entering the process, [DELTA]u and [DELTA]d are the manipulated and the measured disturbances moves at time step i-1, respectively. The first element of Y(i) represents the prediction of the output at time step i. The description of the above matrices is given below.

3. Compute the error vector

[epsilon](i + 1/i) = R(i + 1) - [M.sub.p]Y(i) - [S.sup.nom.sub.n][DELTA]d(i) - w(i + 1/i) (A-3)

where R is the vector of set points, [M.sub.p] is the truncated shift matrix vector, made of the first p*[n.sub.y] rows of [M.sub.n] and w is the unmeasured disturbance vector represented as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A-4)

where p represents the prediction horizon and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] is an identity matrix of dimensions nyxny.

4. Compute the manipulated variable move vector and implement the first element given by:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A-5)

where the MPC gain matrix is obtained from a least squares solution and is expressed as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A-6)

where, the input weight [LAMBDA] = [lambda] x [I.sub.mxm], [lambda] is a user selected scalar, the output weight matrix is set to identity ([GAMMA]=I) for simplicity and [S.sub.p] is the dynamic matrix made of the first p rows of the matrix [S.sub.u] defined below.

5. Return to step 2

Matrices used in the above algorithm:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A-7)

[S.sup.nom.sub.n] = [[[S.sup.nom.sub.u,1], [S.sup.nom.sub.u,2].... , [S.sup.nom.sub.u,n]].sup.T] (A-8)

[S.sup.nom.sub.d] = [[[S.sup.nom.sub.d,1], [S.sup.nom.sub.d,2].... , [S.sup.nom.sub.d,n]].sup.T] (A-9)

[DELTA]U(I) = [[DELTA]u(i), [DELTA]u(i + 1), .... , [DELTA]u(i + m - 1]].sup.T] (A-10)

Y(i)[DELTA] = [[Y(i), y(i + 1), ... , y(i + n - 1)].sup.T] (A-11)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A-12)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A-13)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A-14)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A-15)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A-16)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A-17)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A-18)

u(i) = u(i - 1) + [DELTA]u(i) (A-19)

APPENDIX B: EMPIRICAL CORRELATIONS TO CALCULATE THE TIME CONSTANT AND THE DEAD TIME

The process time constant was calculated using the formulae proposed by Seider et al. (1999). The dominant time constant, [[tau].sub.p], is estimated as:

[[tau].sub.p] = [[tau].sub.col] + [[tau].sub.con] + [[tau].sub.reb] (B-1)

where [[tau].sub.con] and [[tau].sub.reb] are the time constants (minute) associated with the condenser and reboiler, respectively, and [[tau].sub.col] is the time constant (minute) for the column.

The column time constant is estimated according to

[[tau].sub.col] = [N.summation over (i=1)][M.sub.i]/[L.sub.i] (B-2)

where [M.sub.i] is the volumetric holdup ([m.sup.3]) on tray i, [L.sub.i] is the liquid flow rate ([m.sup.3]/min) from tray i, and N is the number of trays. The liquid holdup is expressed as

[M.sub.i] = [A.sub.c] ([h.sub.w] + [h.sub.ow]) = [pi][D.sub.c.sup.2]/4 ([h.sub.w] + [h.sub.ow]) (B-3)

where [D.sub.c] is the column diameter (m), and [h.sub.w] and [h.sub.ow] are the weir height and fluid height above the weir (m), respectively. The latter can be expressed in terms of the weir length, [l.sub.w](m), using the Francis weir equation:

[h.sub.ow] = [([L.sub.i]/111 x [l.sub.w]).sup.2/3] (B-4)

Also the following heuristics are used:

[[tau].sub.con] = [[tau].sub.reb] = 0.5[[tau].sub.col]

[l.sub.w] = 0.65 [D.sub.c]

[h.sub.w] = 2 in

The process dead time, [theta], was calculated using a formulae presented by Seider et al. (1999). A process dead time, [theta], is related to the change in the fluid holdup above the weir. For a single tray, the process dead time each tray, [theta]', can be estimated by considering the time taken for [h.sub.ow] to stabilize after a change in the liquid flow rate.

The calculation of the dead time of each single tray, [theta]', is represented as show in the following equation.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (B-5)

The process dead time for the whole column, [theta], will be the summation of the dead time of each single tray:

[theta] = [N.summation over (1)][theta]' (B-6)

NOMENCLATURE c controller tuning parameters C concentration CC capital costs [d.sub.e] design variables [G.sup.Y] transfer function between disturbance to controlled variable [G.sup.U] transfer function between disturbance to manipulated variable [K.sub.MPC] MPC controller gain matrix [K.sub.p] process gain l model parameter uncertainty M interconnection matrix to test robust stability using [micro] [M.sub.n] shift matrix with dimension nxn [M.sub.p] truncated shift matrix with dimension pxn [M.sub.perf] interconnection matrix to test robust performance using [micro] OC operating costs P price of the product Q volumetric flow RR reflux ratio R(i) set point vector at time interval i [S.sub.n] step response model vector [S.sub.nom] nominal linear step response model u(i) vector of future manipulated variables [u.sub.op] process operating variables [u.sub.k,nom] value of the manipulated variable at the nominal operating point V volume VC variability costs W amplitude of the disturbances w(k) vector of unmeasured perturbations at time step k Y(i) output prediction vector z z-transform Greek Symbols [DELTA] structured perturbation matrix used in the SSV analysis [[DELTA]u.sub.i,k] moves of the ith manipulated variable at time step k calculated by the MPC controller [DELTA]d change in the measured disturbance [epsilon] prediction error vector used for MPC calculations [omega] frequency [LAMBDA] controlled variable weight in MPC objective function [GAMMA] manipulated variable weight in MPC objective function [micro] structured singular value [theta] dead time between manipulated to controlled variables [sigma] maximum singular value [[tau].sub.p] process time constant [delta]C deviation in concentration with respect to the nominal value [beta] scalar weight multiplying performance weight [lambda] diagonal elements of the weight matrix [GAMMA] (assumed equal) Indexes i time interval n number of prediction steps nd number of measured disturbances nu number of manipulated variables or inputs ny number of controller variables or outputs m control horizon, number of future input moves p prediction horizon l complex number T sampling period

Manuscript received January 29, 2007; revised manuscript received April 23, 2007; accepted for publication May 30, 2007.

REFERENCES

Bansal, V., J. D. Perkins, E. N. Pistikopoulos, R. Ross and J. M. G. Schijndel, "Simultaneous Design and Control Optimization under Uncertainty," Comput. Chem. Eng. 24, 261-266 (2000).

Blanco, A. M. and J. A. Bandoni, "Interaction between Process Design and Operability of Chemical Processes: An Eigenvalue Optimization Approach," Comput. Chem. Eng. 27, 1291-1301 (2003).

Chawankul, N., P. L. Douglas and H. Budman, "The Integration of Design and Control: IMC and Robustness," Comput. Chem. Eng. 29(2), 261-271 (2005).

Doyle, F. J., A. K. Packard and M. Morari, "Robust Controller Design for a Nonlinear CSTR," Chem. Eng. Sci. 44(9), 1929-1947 (1989).

Georgiadis, M. C., M. Schenk and E. N. Pistikopoulos, "The Interactions of Design, Control and Operability in Reactive Distillation Systems," Comput. Chem. Eng. 26, 735-746 (2002).

Lee, L. P., "Nonlinear Process Control," Springer-Verlag, London (1993).

Lee, L. P., H. Li and T. Cameron, "Decentralized Control Design for Nonlinear Multi-Input Plants: A Gap Metric Approach," Chem. Eng. Sci. 55(18), 3743-3758 (2000).

Lewin, D. R., "A Simple Tool for Disturbance Resiliency Diagnosis and Feedforward Control Design," Comput. Chem. Eng. 20(1), 13-25 (1996).

Luyben, M. L. and C. A. Floudas, "Analyzing the Interaction of Design and Control-1: A Multi Objective Framework and Application to Binary Distillation Synthesis," Comput. Chem. Eng. 18, 933-969 (1994).

Mohideen, M. J., J. D. Perkins and E. N. Pistikopoulos, "Optimal Design of Dynamic Systems under Uncertainty," AIChE J. 42, 2251-2272 (1996).

Morari, M. and E. Zafiriou, "Robust Process Control," Prentice Hall, New Jersey (1989).

Sakizlis, V., J. D. Perkins and E. N. Pistikopoulos, "Recent Advances in Optimization-Based Simultaneous Process and Control Design," Comput. Chem. Eng. 28, 2069-2086 (2004).

Seferlis, P. and M. C. Georgiadis, "The Integration of Process Design and Control," 1st ed., Elsevier, Amsterdam (2004).

Seider, W. D., J. D. Seader and D. R. Lewin, "Process Design Principle," Wiley, New York (1999).

TableCurve 2D, Systat Software Inc.[R], San Jose, CA, U.S.A. (2002).

Zanovello, R. and H. M. Budman, "Model Predictive Control with Soft Constraints with Application to Lime Kiln Control," Comput. Chem. Eng. 23, 791-806 (1999).

Zhou, K. and J. C. Doyle, "Essentials of Robust Control," 1st ed., Prentice Hall, New Jersey (1997).

N. Chawankul, L. A. Ricardez Sandoval, H. Budman * and P. L. Douglas Department of Chemical Engineering, University of Waterloo, Waterloo, ON, Canada N2L 3G1

* Author to whom correspondence may be addressed. E-mail address: hbudman@uwaterloo.ca

Table 1. Traditional method: sensitivity of the total cost with respect to the product price ([theta] [not equal to] 0) W1 W2 RR * Q * (lbmole/h) [lambda] * N D (m) CC ($/day) 1 1 1.918 3978 0.50523 26 3.393 195.56 5 1 1.918 3978 0.50523 26 3.393 195.56 10 1 1.918 3978 0.50523 26 3.393 195.56 15 1 1.918 3978 0.50523 26 3.393 195.56 20 1 1.918 3978 0.50523 26 3.393 195.56 1 5 1.918 3978 0.50523 26 3.393 195.56 1 10 1.918 3978 0.50523 26 3.393 195.56 1 15 1.918 3978 0.50523 26 3.393 195.56 1 20 1.918 3978 0.50523 26 3.393 195.56 W1 W2 OC ($/day) VC ($/day) TC ($/day) 1 1 586.61 65.61 847.78 5 1 586.61 139.35 921.52 10 1 586.61 237.09 1019.26 15 1 586.61 337.74 1119.91 20 1 586.61 434.92 1217.09 1 5 586.61 245.09 1027.26 1 10 586.61 473.24 1255.41 1 15 586.61 705.24 1487.41 1 20 586.61 923.69 1705.86 Table 2. Integrated method: sensitivity of RR *, Q * and [lambda] * with respect to the product price ([theta] [not equal to] 0) W1 W2 RR * Q * (lbmole/h) [lambda] * N D (m) CC ($/day) 1 1 1.913 3975.5 0.2350 26 3.392 195.98 5 1 1.911 3973.4 0.2341 26 3.392 196.35 10 1 1.908 3971.5 0.2338 27 3.391 196.70 15 1 1.753 3880.2 0.2331 38 3.370 261.15 20 1 1.753 3880.2 0.2332 38 3.370 261.14 1 5 1.912 3974.0 0.1886 26 3.392 196.23 1 10 1.909 3972.2 0.1848 26 3.391 196.56 1 15 1.906 3970.6 0.1836 27 3.391 196.88 1 20 1.904 3969.0 0.1830 27 3.390 197.19 W1 W2 OC ($/day) VC ($/day) TC ($/day) 1 1 586.21 23.34 805.53 5 1 585.87 44.11 826.33 10 1 585.56 70.89 853.15 15 1 570.34 93.18 924.67 20 1 570.35 118.87 950.36 1 5 585.98 35.10 817.31 1 10 585.68 55.61 837.85 1 15 585.41 76.03 858.32 1 20 585.16 96.59 878.94 Table 3. Integrated method: without control ([lambda] = 1000) W1 W2 RR * Q * (lbmole/h) N D (m) CC ($/day) 1 1 1.7014 3851.5 51 3.365 331.65 5 1 1.7014 3851.5 51 3.365 331.65 10 1 1.7014 3851.5 51 3.365 331.65 15 1 1.7014 3851.5 51 3.365 331.65 20 1 1.7014 3851.5 51 3.365 331.65 1 5 1.7014 3851.5 51 3.365 331.65 1 10 1.7014 3851.5 51 3.365 331.65 1 15 1.7014 3851.5 51 3.365 331.65 1 20 1.7014 3851.5 51 3.365 331.65 W1 W2 OC ($/day) VC ($/day) TC ($/day) 1 1 565.40 621.70 1518.75 5 1 565.40 1687.30 2584.30 10 1 565.40 3019.20 3916.30 15 1 565.40 4351.20 5248.20 20 1 565.40 5683.10 6580.15 1 5 565.40 2042.90 2940.00 1 10 565.40 3819.50 4716.60 1 15 565.40 5596.00 6493.10 1 20 565.40 7372.60 8269.65 Table 4. Compare variability costs from the design with and without control V[C.sub.without MPC]/ W1 W2 V[C.sub.with MPC] 1 1 26.65 5 1 38.25 10 1 42.59 15 1 46.69 20 1 47.81 1 5 58.20 1 10 68.68 1 15 72.85 1 20 76.33 Table 5. The sensitivity of the first RGA element [[lambda].sub.11] with respect to the product price, ([theta] [not equal to] 0) W1 W2 RR * [lambda]11 1 1 1.913 3.65 5 1 1.911 3.63 10 1 1.908 3.62 15 1 1.753 3.03 20 1 1.753 3.03 1 5 1.912 3.64 1 10 1.909 3.62 1 15 1.906 3.61 1 20 1.904 3.60 Table 6. Savings from the integrated design method as compared to the traditional design methods TC--integrated TC--traditional W1 W2 method ($/day) method ($/day) % savings 1 1 805.53 847.78 5.24 5 1 826.33 921.52 11.52 10 1 853.15 1019.26 19.47 15 1 924.67 1119.91 21.11 20 1 950.36 1217.09 28.07 1 5 817.31 1027.26 25.69 1 10 837.85 1255.41 49.84 1 15 858.32 1487.41 73.29 1 20 878.94 1705.86 94.08