Heuristic-Based Mathematical Programming Framework for Control

Ind. Eng. Chem. Res. 2001, 40, 2079-2088

2079

Heuristic-Based Mathematical Programming Framework for Control Structure Selection Ioannis K. Kookos and John D. Perkins* Centre for Process Systems Engineering, Department of Chemical Engineering, Imperial College, London SW7 2BY, U.K.

This paper presents an optimization-based method for selecting manipulated variables for regulatory control schemes. The objective of this mathematical programming technique is the minimization of the overall interaction and sensitivity of the closed-loop system to disturbances. In addition, a general methodology for incorporating qualitative knowledge as linear constraints to the problem is demonstrated. The main advantage of the method is that the plantwide nature of the problem is preserved, because decisions related to different levels of the base regulatory control scheme are made simultaneously. The usefulness of the method is demonstrated using a double-effect evaporator case study, the hydrodealkylation of toluene case study, and the Tennessee Eastman challenge problem. Introduction The design of plantwide regulatory control systems is a creative challenge that has received a lot of attention in recent years. Pressures to reduce capital investment and working capital as well as safety and environmental concerns have led to highly integrated designs for industrial processes. As a result, the complexity involved in the design of regulatory control systems has increased. The coordination of decisions related to the control of isolated process units has proven to be a nontrivial task in a plantwide environment. Furthermore, particularly effective control schemes developed for simple process units may become totally ineffective or even incompatible in a plantwide context. In the past decade, a number of methodologies have been proposed in the chemical engineering literature for the generation of promising plantwide regulatory control structures. These methodologies range from pure mathematical programming based methods1 to heuristicbased methods.2 Other interesting approaches such as decomposition3,4 or hierarchical approaches5,6 have also been presented. A characteristic shared by most of these methodologies is the division of the plant and its control system into subsets. This is not a new idea and can be found even in the early and inspiring development of the concept of material balance control by Buckley,7 which is still dominant in most approaches. The aim of this decomposition is to divide the general and usually intractable problem into a number of serial or parallel subproblems, amenable to systematic solutions. The design of the base control scheme, defined as the set of regulatory control loops related to the energy management, inventory, production level, and quality control has been identified as a crucial step in synthesizing workable plantwide control schemes.3,4,7-9 This is due to the fact that these elements of the regulatory * To whom all correspondence should be addressed. Tel: +44(0)207 5945556. Fax: +44(0)207 5945639. E-mail: [email protected].

control scheme are related to fast disturbances, and thus the overall performance of a plantwide control scheme depends strongly on the performance of these base level control loops. Tyreus10 presented a methodology for identification of dominant variables for partial control. This methodology can be useful for the development of the base control schemes where a limited set of dominant variables are controlled close to their set points. The selection of suitable sets of manipulated variables for controlling these dominant variables is crucial for the success or failure of this base control scheme. However, as can be seen from a recent paper by McAvoy,9 the selection of this set of manipulated inputs is a challenging problem because of its combinatorial nature. When ny manipulated variables have to be chosen among nu potential manipulated variables and assigned to ny controlled variables, then the number of alternatives is given by

nu!/(nu - ny)!

(1)

For the base control system of a medium-scale industrial process such as the Tennessee Eastman (TE) process,11 this number can be as large as 4 × 107. This paper presents a systematic methodology for selecting manipulated variables for regulatory control schemes. The proposed methodology is based on the formulation of a linear mathematical programming problem for the selection of promising sets of manipulated variables. The objective of this mathematical programming technique is the minimization of the overall interaction as expressed in terms of the steady-state relative gain array12 (RGA) and the minimization of the sensitivity of the closed-loop system to simultaneously acting disturbances. In addition, a general methodology, first proposed by Raman and Grossmann,13 for incorporating qualitative knowledge (heuristics) as constraints to the process synthesis problem is demonstrated. The main advantage of the method is that the plantwide nature of the problem is preserved. This is achieved by the simultaneous consideration of all deci-

10.1021/ie000415t CCC: $20.00 © 2001 American Chemical Society Published on Web 04/10/2001

2080

Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001

sions related to different levels of the structure of the base control system. Furthermore, process insight and knowledge about the problem under consideration can easily be incorporated into the mathematical programming framework in a highly structured manner. Thus, the effects of considering alternative heuristics can be analyzed and assessed quite effectively. It should be noted that when more than one alternative sets of basic controlled variables are available, the proposed methodology can be applied for each one of them, and thus all alternatives can be evaluated and compared. The structure of this paper is as follows. First, a mixed-integer, linear programming (MILP) formulation for the calculation of the RGA is presented followed by the development of the optimization problem for simultaneous consideration of the multivariable interaction and sensitivity to disturbances. Next, a general methodology for incorporating qualitative knowledge into the MILP formulation is demonstrated through a number of simple examples. Finally, the methodology is applied to a double-effect evaporator case study, the hydrodealkylation of toluene (HDA) case study, and the TE challenge problem. 2. MILP Formulation for Control Structure Selection In this paper we consider systems described by the following transfer function matrixes

y(s) ) G(s) u(s) + Gd(s) d(s)

(2)

where y(s) is the Laplace transform of the ny vector of controlled variables y(t) selected for the base control scheme, u(s) is the Laplace transform of the nu vector of potential manipulated variables u(t), and d(s) is the Laplace transform of the nd vector of disturbances d(t). This mathematical description can be obtained either by identification experiments, such as step or pulse response, or from the corresponding linear state space description. The analysis that follows makes use of the steady-state description of the process only, and it is assumed that the number of available manipulated variables is at least equal to the number of controlled variables. 2.1. MILP Formulation for Calculation of the RGA. Define the indices and sets i, l ) indices over the set I of measurements that must be controlled, with I ) {1, 2, ..., ny}, j ) index over the set J of potential manipulated variables, with J ) {1, 2, ..., nu}, and m ) index over the set M of disturbances, M ) {1, 2, ..., nd}, and assume that the ny × nu transfer function matrix that relates the potential manipulated with the controlled variables is given by

[

g11 g12 g g G ) 21 22 l l gny1 gny2

... g1nu ... g2nu l l g ... nynu

]

(3)

where nu g ny. Because the transfer function matrix is not square, it is also not invertible. However, ny out of the nu potential manipulated variables are to be used in order to form the base regulatory control system. For any feasible selection of ny manipulated inputs, we can ˆ that consists of the ny define an ny × ny matrix G

columns of the G matrix that correspond to these inputs. We will require this matrix to be invertible in order for the selection to be feasible. In this work, the formation of this matrix is done in the following way. First, the following matrix is defined:

[

X11 X12 X21 X22 X) l l Xny1 Xny2

... X1nu ... X2nu l l ... Xnynu

]

(4)

The elements of the X matrix are binary variables associated with the following logic statements:

Xij )

[

1, if manipulated variable j is used to control controlled variable i 0, otherwise

]

Using the definition above, we can deduce the following: ny

Xij ) ∑ i)1

[

1, then manipulated variable j is used to control one of the measurements 0, then manipulated variable j is not used to control any of the measurements

]

Then we define the following ny × nu matrix which includes the columns of the transpose of the inverse of the G ˆ matrix plus nu - ny zero columns

[

g˜ 11 g˜ 12 g˜ g˜ G ˜ ) 21 22 l l g˜ ny1 g˜ ny2

... g˜ 1nu ... g˜ 2nu l l g ˜ ... nynu

]

(5)

The elements of this matrix are defined as follows: nu

gijg˜ lj - δil ) 0, ∑ j)1

∀ i, l

(6)

where δil is the Kronecker delta. It should be noted that G ˜ is not a square matrix. However, it consists of ny nonzero columns plus nu - ny zero columns. The ny nonzero columns correspond to the manipulated variables selected in the structure, while the nu - ny zero columns correspond to the manipulated variables not selected in the structure. The nonzero columns of G ˜ constitute the transpose of the inverse of the ny - ny matrix G ˆ . Using this approach, the inverse of any ny ny submatrix of G may be effectively generated. To enforce the fact that the columns for which the ny Xlj is zero have zero elements, we add the sum Σl)1 bounds ny

-Ω(

ny

Xlj) e g˜ lj e Ω(∑Xlj), ∑ l)1 l)1

∀ i, j

(7)

where Ω is a sufficiently large number. Having calculated the elements of the transpose of the inverse of the

Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2081

G ˆ matrix, then the elements λij of the RGA matrix can be calculated as

λij - gijg˜ ij ) 0, ∀ i, j

ny nu

∑ ∑|λij - Xij| i)1 j)1

||Λ - I||sum )

(8)

(15)

Furthermore, define the matrix

In addition, the following constraints denote the fact that one and only one manipulated variable has to be assigned to each controlled output:

ηij ) λij - Xij, ∀ i, j

(16)

as well as the absolute values of its elements µij, determined so that

ny

Xij - 1 e 0 ∑ i)1

(9)

nu

-µij e ηij e µij, ∀ i, j µij g 0

(17)

Next, observe that eq 15 can be written as follows:

Xij - 1 ) 0 ∑ j)1

(10) ny nu

Thus, we have presented a MILP formulation for the effective calculation of the RGA matrix. In the following section, this formulation will be used in order to assist us in solving the control structure selection problem based on the criteria of minimum interaction, sensitivity to input uncertainty, and disturbance rejection. 2.2. MILP Formulation for Minimization of the Overall Interaction. A multivariable process is said to have interaction when manipulated variables affect more than one controlled variable and RGA is commonly used as a measure of interaction. Specifically, input and output variables should be paired so that the diagonal elements of the RGA are as close to 1 as possible. For multivariable systems, this intuitive measure of interaction can be expressed as follows:

||Λ - I||sum )

∑ ∑µij i)1 j)1

(18)

Thus, to find the structure that minimizes the overall interaction, the following MILP should be solved: ny nu

∑∑

min µij Xij,µij i)1j)1

(19)

ny

Xij - 1 e 0, ∑ i)1

∀j

nu

Xij - 1 ) 0, ∑ j)1

∀i

ny

|λii - 1| ∑ i)1

(11)

For small sensitivity to input uncertainty, RGA must have a small induced ∞-norm (maximum row sum) defined as14

||Λ||i∞ max ( i

∑j |λij|)

(12)

Furthermore, there is a close relationship between the minimized condition number, a measure of the illconditioning of a plant and the magnitude of the RGA elements15,16

2 max {||Λ||i1, ||Λ||i∞} -

1 γmin

e γmin e ||Λ||sum + k (13)

where ||Λ||i1 ) maxj (∑i|λij|), ||Λ||sum ) ∑i,j|λij|, and k equals the size of the matrix Λ minus 2. It is particularly noteworthy that all of these measures are closely related to the absolute values of the elements of the RGA matrix. Because of this fact, these objectives are not conflicting. Thus, in this work the selection of promising sets of manipulated variables and pairings is based on the minimization of the following index:14 ny

||Λ - I||sum )

|λii - 1| + ∑|λij| ∑ i)1 i,j

(14)

i*j

From eq 14 and the definition of the matrix X, we obtain

[

nu

gijg˜ lj - δil ) 0, ∑ j)1 ny

ny

Xlj) e g˜ lj e Ω(∑Xlj) ∑ l)1 l)1

-Ω(

λij - gijg˜ ij ) 0 ηij - λij + Xij ) 0 -µij - ηij e 0 ηij - µij e 0

Xij ) {0, 1}, ∀ i, j,

]

∀ i, l

∀ i, j

µij g 0, ∀ i, j

2.3. MILP Formulation for Minimization of the Sensitivity to Disturbances. A further objective that can be posed in the control structure selection problem is the rejection of disturbances. A measure of the ability of a specific control structure to reject disturbances is the following:14

||G ˆ -1 Gd||i∞

(20)

That is, the induced ∞-norm of the product of the inverse of the plant transfer function and the disturbance transfer function is a measure of the sensitivity of a specific structure to disturbances. The smaller the induced ∞-norm, the more capable the system is to reject disturbances. To form the objective function in this case, we need to construct the inverse of the ny × ny matrix that corresponds to the selected columns (or manipulated

2082

Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001

inputs) of the original G matrix. However, this inverse matrix has been calculated as part of the formulation for the calculation of the RGA matrix. Define the following ny × nd matrix:

S)G ˆ -1Gd

formulation can be used: ny nu

∑∑

min µij} {Fφ + (1 - F) Xij,µij,j,φ i)1 j)1

(21)

(25)

ny

Xij - 1 e 0, ∑ i)1

Define also the nu × nd matrix Σ as a matrix defined by the following equation:

∀j

nu

Xij - 1 ) 0, ∑ j)1

ny

g˜ ijgd,im - σjm ) 0, ∑ i)1

∀ j, m

(22)

It should be observed that the Σ matrix consists of the ny rows of the S matrix that correspond to the inputs selected in the structure as well as nu - ny additional zero rows that correspond to the inputs not selected in the structure. Because of this property, it follows that

||S||i∞ ) ||Σ||i∞

(23)

In summary, the following MILP formulation can be used in order to find the set of manipulated inputs that is least sensitive to disturbances

min φ Xij,j,φ

(24)

∀j

Xij - 1 ) 0, ∀ i

-Ω(

∀ i, j

ny

g˜ ijgd,im - σjm ) 0, ∑ i)1

[

nd

-j e

∑ σjm e j

m)1

j - φ e 0

Xij ) {0, 1}, ∀ i, j,

∀ j, m

]

ny

Xlj) ∑ l)1

Xlj) e g˜ ij e Ω(

λij - gijg˜ ij ) 0 ηij - λij + Xij ) 0 -µij - ηij e 0 ηij - µij e 0 ny

g˜ ijgd,im - σjm ) 0, ∑ i)1 nd

∑ σjm e j

m)1

j - φ e 0

]

∀ i, j

∀ j, m

]

∀j

µij g 0, ∀ i, j,

j g 0, ∀ j

3. Incorporation of Qualitative Knowledge

ny

Xlj) e g˜ ij e Ω(∑Xlj), ∑ l)1 l)1

∑ l)1

∀ i, l

The weighting coefficient F ∈ [0, 1) is used in order to assign different weights to the contribution of interaction or disturbance sensitivity to the objective function.

gijg˜ lj - δil ) 0, ∀ i, l

ny

ny

-Ω(

Xij ) {0, 1}, ∀ i, j,

nu

∑ j)1

gijg˜ lj - δil ) 0, ∑ j)1

-j e

nu

∑ j)1

nu

[

ny

Xij - 1 e 0, ∑ i)1

[

∀i

∀j

j g 0, ∀ i, j

It should be noted that formulation (24) does not give an answer to the input-output pairing problem but only the selection of a set of manipulated variables. This is due to the fact that the disturbance sensitivity measure used and defined in eq 20 does not depend on the pairing of the input-output variables but only on the sets of controlled and manipulated variables selected. Finally, to simultaneously consider the interactions as well as the disturbance sensitivity, the following

In the previous sections an MILP formulation was presented for inventing promising control structures based on the ideas of minimum overall interaction and disturbance sensitivity. Process control engineers have developed a number of heuristics for the quick generation of promising control structures. It is particularly desirable to be able to incorporate a number of these heuristic rules or qualitative knowledge into the proposed formulation for control structure selection. The objective is to use qualitative knowledge in order (a) to expedite the solution of the resulting MILP problem by reducing the size of the feasible region and (b) to eliminate a number of promising (from the mathematical programming point of view) but undesirable (from the practical point of view) control structures. Raman and Grossmann13 were the first to propose the introduction of qualitative knowledge into the solution of process synthesis problems using mathematical programming methods. They presented a systematic method, first discussed in Cavalier et al.,17 for converting complex logical expressions into equivalent mathematical representations. The interested reader is referred to Raman and Grossmann13 for a comprehensive presentation of the general methodology. In the following,

Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2083

tank. The corresponding Boolean expressions are as follows:

P31 w P12 ∧ P23

(32)

(¬P31 ∨ P12) ∧ (¬P31 ∨ P23)

(33)

or Figure 1. Two-tank process example.

the application of this methodology in transforming a number of well-known heuristics into mathematical expressions (linear constraints) is presented. A well-established heuristic in designing plantwide control structures is related to the construction of “selfconsistent” inventory control structures.3,4,7 An inventory control structure is said to be “self-consistent” if it is able to propagate a production rate change throughout the process so that such a change produces changes in the flow rates of all major feed and product streams. Self-consistency requires that the chain of level controls be constructed to radiate outward from the throughput manipulator. This implies that level controllers between the feed and the throughput manipulator are in the direction opposite to flow, while those between the throughput manipulator and the main product are in the direction of flow. Furthermore, throughput manipulators take two forms. Process streamflows can be throttled directly to serve as “explicit” manipulators (such as feed and product streams), or other process variables (often utilities) can be used as “implicit” throughput manipulators. The mathematical formulation of the requirement for self-consistent inventory control schemes is derived next for a simple example. Consider the simple process shown in Figure 1 and consisting of two tanks in series. The objective is to control the liquid levels in the tanks (y1 and y2) as well as the product flow rate (y3) by using the openings of the available valves (u1, u2, u3, and u4). The literal Pij is used to represent the following logical statement: manipulated variable j is used to control controlled variable i, and we assign to this literal the binary variable Xij. First, observe that either u1, u2, or u4 can be used to control y1, either u2, u3, or u4 can be used to control y2, and either u1, u2, or u3 can be used to control y3. These statements are equivalent to the Boolean expressions

P11 ∨ P12 ∨ P14

(26)

P22 ∨ P23 ∨ P24

(27)

P31 ∨ P32 ∨ P33

(28)

where ∨ is the EXCLUSIVE OR operator. Translating each expression into the equivalent mathematical form, we obtain (note that x ∨ y is equivalent to x + y ) 1)

X11 + X12 + X14 ) 1

(29)

X22 + X23 + X24 ) 1

(30)

X31 + X32 + X33 ) 1

(31)

The heuristic for self-consistent inventory control structures requires that when u1 is used as a throughput manipulator, then u2 should be used to control the level in the first tank and u3 to control the level in the second

where ∧, ¬, and ∨ are the AND, NOT, and OR operators. Similar arguments apply when u2 or u3 is used as a throughput manipulator. The equivalent mathematical expressions for all possible throughput manipulators are summarized below (note that x ∨ y is equivalent to x + y g 1):

-X12 + X31 e 0, -X23 + X31 e 0

(34)

-X11 + X32 e 0, -X23 + X32 e 0

(35)

-X11 + X33 e 0, -X22 + X33 e 0

(36)

These constraints are hard constraints because no feasible solution can violate them. However, most of the heuristics are only recommendations (soft constraints) and not proven facts (hard constraints). Thus, it would be desirable to allow violation of these heuristics. This can be done by modifying these constraints and introducing the binary slack variables ν13

-(X12 + ν31) + X31 e 0, -(X23 + ν31) + X31 e 0 (37) -(X11 + ν32) + X32 e 0, -(X23 + ν32) + X32 e 0 (38) -(X11 + ν33) + X33 e 0, -(X22 + ν33) + X33 e 0 (39) It should be obvious that when the heuristic is violated, the corresponding slack variable is forced to have the value of 1. The objective function is also modified to include terms that penalize the weighted violation of the heuristics

w(ν31 + ν32 + ν33)

(40)

where w represents weights that reflect the relative importance of violating the corresponding heuristic. Consider the more general case of n tanks in series where yi is the level of the ith tank and yn+1 is the product flow rate. To control the n levels and the product flow rate, we define for each level the set of available manipulated variables in the direction opposite to flow (OFi) as well as the set of manipulated variables in the direction of flow (DFi). Then for each selection of a total production manipulator, the requirement of a selfconsistent level control structure can be expressed using the following Boolean expressions:

Pn+1,j w

Plk ∧ lj,k∈DFi

(41)

i.e., level controllers between the feed and the throughput manipulator are in the direction opposite to flow, while those between the throughput manipulator and the main product are in the direction of flow. It is noteworthy that eq 41 is applicable to processes having only one major feed stream and one major product stream. In cases where more than one feed or product

2084

Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001

stream (or both) are involved, eq 41 applies to each path connecting any feed to any product stream. However, in this case, it might not be possible to construct selfconsistent inventory control structures for all product streams. Following a similar procedure, the well-known heuristic according to which input-output pairings that correspond to negative RGA elements should be avoided has the following equivalent mathematical description:

-Ω(1 - Xij) e λij

(42)

where Ω is a sufficiently large number. When the heuristic is violated (λij < 0), then the binary variable Xij is forced to have the value of zero. Furthermore, to allow violation of this constraint, a binary slack variable can be used

-Ω(1 + νij - Xij) e λij

(43)

together with an appropriate weighting term included in the objective function. Concluding, assume that the specific logical relationships among the controlled variables and the potential manipulated variables are expressed in the form of logical propositions, namely, by a set of conjunctions of clauses. Then, the logical propositions are converted into conjunctive normal form (CNF) or disjunctive normal form (DNF) as illustrated in the examples presented above (see Grossmann and Daichendt18 and the reference therein for details). Next, based on the one-to-one correspondence with the binary variables that are used to denote variable pairing, the heuristic rules are transformed into linear constraints. Finally, formulation (25) together with the integer constraints presented in this section has the general form T T T min c x + e X + w v x,X,v

(44)

s.t. Ax + BX ) b Cx + DX + Ev e f x∈X νij, Xij ) {0, 1} This is a MILP problem for which effective methods and software for solving to global optimality are available. The reader is referred to Williams19 and Wolsey20 for comprehensive treatment of the related theory and useful guidelines for efficient use of the available software. 4. Case Studies 4.1. Double-Effect Evaporator Case Study. This case study considers the concentration operation shown in Figure 2. The concentration of a solution of A is to be raised from 3.5% to 10% in a double evaporation process. The liquor from an upstream process is fed into the first evaporator through the feed pump. A utility steam is supplied to the first effect, and steam to the second effect is supplied by vapor from the first effect. The product liquor is fed through the product pump to a downstream storage facility. The details of the mathematical model used to describe the dynamics of the double-effect evaporator process can be found in Narr-

Figure 2. Double-effect evaporator process.

away.21 The manipulated variables for this case study are the positions of the six valves, while the measured variables include the composition, level, temperature, pressure, steam chest temperature, and pressure in each evaporator. The disturbances are variations in the feed pressure, temperature and composition as well as the utility steam temperature. Previous studies on this process1 have shown that effective control can be achieved by controlling the two levels and the pressure and temperature in the second effect. It is noteworthy that this structure does not rely upon any composition analyzer, a characteristic particularly desirable from the cost and maintenance point of view. The steady-state transfer functions are as follows:

[

G(0) ) 4.7226 -6.3975 -3.7862 -1.9313 -1.9972 -2.6689 2.5032 10.5727 -9.2682 -2.0181 -3.0346 -3.1967 -0.0300 -0.1232 -0.0750 0.4702 0.6887 -4.0700 0.0262 0.1178 0.0669 0.7827 1.1365 -8.4526

[

0.3837 0.2034 Gd(0) ) -0.0024 0.0021 where

-0.1670 -0.1745 0.0407 0.0677

0.0219 0.0282 0.0176 -0.0120

-0.2403 -0.2511 0.0585 0.0974

[][

y1 liquid level in the first effect y2 liquid level in the second effect y) y ) temperature in the second effect 3 y4 pressure in the second effect

[][

u1 opening of liquid valve 0 u2 opening of liquid valve 1 u3 opening of liquid valve 2 u) u ) opening of vapor valve 0 4 u5 opening of vapor valve 1 u6 opening of vapor valve 2

[][

d1 feed pressure d2 feed temperature d) d ) feed composition 3 d4 steam temperature

]

]

]

]

]

(45)

(46)

(47)

(48)

(49)

These matrixes have been scaled so that

||u||∞ max |uj| e 1, ||d||∞ max |dm| ) 1 j m

(50)

Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2085

For the system under consideration, 360 alternative control structures can be obtained. Thus, even for this relatively small example, the evaluation of all alternatives (complete enumeration) requires considerable time and effort. The MILP formulation given by formulation (25) was applied to this case study without using any qualitative knowledge. The value used for the large number Ω used in the formulation is 100 (the same value was found to be appropriate for all case studies). For F close to 1, the following solution is obtained using the GAMS interface with CPLEX22 in 0.42 s (Ultra Sparc 360 MHz)

y1 y X ) y2 3 y4

[

u1 1 0 0 0

u2 0 0 1 0

u3 0 1 0 0

u4 0 0 0 1

u5 0 0 0 0

u6 0 0 0 0

]

(51)

Table 1. HDA Process Potential Manipulated Variables

i.e., the pairings y1-u1, y2-u3, y3-u2, and y4-u4 are obtained. The measure of sensitivity to disturbances has the value ||G ˜ -1Gd||i∞ ) 0.222, and the corresponding RGA is given by

y1 y Λ ) y2 3 y4

[

u1 0.7583 0.0002 0.1531 0.0884

u3 -0.0002 0.5959 0.2699 0.1344

u2 0.2424 0.4034 0.2216 0.1326

u4 -0.0004 0.0005 0.3553 0.6446

]

(52)

For F close to zero, the following solution is obtained:

y1 y X ) y2 3 y4 y1 y Λ ) y2 3 y4

[

[

u1 1 0 0 0

u1 0.7570 0.0007 0.1715 0.0708

u2 0 0 1 0

u3 0 1 0 0

u3 0.0000 0.5954 0.2819 0.1227

u4 0 0 0 0

u5 0 0 0 0

u6 0 0 0 1

u2 0.2428 0.4040 0.2420 0.1111

Figure 3. HDA process flow sheet.

]

u6 0.0000 0.0000 0.3046 0.6954

(53)

]

variable

variable name

u1 u2 u3 u4 u5 u6 u7 u8 u9 u10 u11 u12 u13

product column reflux ratio product column reboiler vapor ratio compressor power cooling water flow gas feed flow toluene feed flow separator liquid valve opening separator vapor valve opening separator inlet valve opening purge valve opening stabilizer column reboiler vapor ratio recycling column reflux ratio recycling column reboiler vapor ratio

Table 2. HDA Process Selected Controlled Variables variable

variable name

y1 y2 y3 y4 y5

separator inlet temperature production rate product purity reactor H2/aromatics ratio separator outlet pressure

Table 3. HDA Process Disturbances

(54)

i.e., the pairings y1-u1, y2-u3, y3-u2, and y4-u6 are obtained with ||G ˆ -1Gd||i∞ ) 0.3968. These results are in close agreement with previous results presented by Heath23 for the same case study. The liquid level and temperature control is assigned to the liquid valves and the pressure control to one of the vapor valves. In summary, for this case study, the proposed methodology gives results consistent with those of previous studies in negligible computational time. 4.2. HDA Process Case Study. The process for HDA to produce benzene is shown in Figure 3. The detailed design of this process is presented in work by Douglas.24 The process consists of nine units. The fresh feeds of toluene and hydrogen are mixed with the two recycle streams and fed into the reactor. Two vapor-phase reactions generate benzene (product) and methane and diphenyl (byproducts). The reactor effluent is quenched. Then, the reactor products are fed into the separation section consisting of a vapor/liquid separator, a stabilizer column where the methane is removed, a product column where the benzene is removed, and a recycle

variable

variable name

d1 d2

cooling water temperature downstream pressure

column for recycling the unreacted toluene and removing the undesirable byproduct (diphenyl). Cao and Rossiter25,26 have considered the selection of five manipulated variables, among the 13 potential manipulated variables given in Table 1, for controlling the measurements given in Table 2. The steady-state gain matrices relating the measured variables with the manipulated variables and disturbances can be found in Cao and Rossiter.26 The application of formulation (25) is considered for this case study. The number of alternative control structures for this case study is 1.5444 × 105. Formulation (25) was solved for F ) 0.5 in 1.28 s (Ultra Sparc 360 MHz) using the GAMS interface with CPLEX,22 and the results are summarized in Table 4. From this table one can observe that the pairings y1-u4, y2-u5, y3-u11 or y3-u1, and y5-u10 are dominant in all promising control structures. This is in close agreement with the results and simulations previously presented by Cao and Rossiter25,26 on the same case study. 4.3. TE Process Case Study. In Figure 4, a schematic representation of the TE process is given. A detailed description of this process can be found in

2086

Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001

Table 4. Best Control Structures for the HDA Process 1 2 3 4 5 6 7 8 9 10

Table 5. TE Process Potential Manipulated Variables

y1

y2

y3

y4

y5

objective

variable

variable name

u4 u4 u4 u4 u4 u4 u4 u4 u4 u4

u5 u5 u5 u5 u5 u5 u5 u5 u5 u5

u11 u11 u11 u11 u11 u11 u11 u1 u1 u1

u9 u8 u3 u1 u2 u6 u12 u13 u9 u8

u10 u10 u10 u10 u10 u10 u10 u10 u10 u10

1.3053 1.3244 1.3810 1.4095 1.4095 1.4598 1.5126 1.5245 1.7250 1.7464

u1 u2 u3 u4 u5 u6 u7 u8 u9 u10 u11 u12

D feed flow E feed flow A feed flow A and C feed flow compressor recycle valve purge valve separator pot liquid flow stripper liquid product flow stripper steam valve reactor cooling water temperature set point condenser cooling water flow agitator speed

Table 6. TE Process Selected Controlled Variables variable

variable name

y7 y8 y9 y12 y15 y17 y30 y42

reactor pressure reactor level reactor temperature product separator level stripper level stripper underflow (production rate) B composition in the purge ratio G/H in the product (product mix)

Table 7. TE Process Potential Manipulated Variables for Each Controlled Variable Figure 4. TE process.

Vogel.11

Downs and The process has four gaseous feeds, a product stream, and a purge stream to remove an inert component. It consists of four main process units: a twophase chemical reactor, a condenser, a vapor/liquid separator, and a stripper. A catalyst, which resides in the liquid phase, promotes the following exothermic gasphase reactions:

A(g) + C(g) + D(g) f G(l)

product 1

A(g) + C(g) + E(g) f H(l)

product 2

A(g) + E(g) f F(l) 2D(g) f F(l)

byproduct byproduct

(55)

The heat of reaction causes the liquid fed to the reactor to vaporize. The vapor leaving the reactor is fed into a water-cooled condenser. The partially condensed stream is fed into a vapor/liquid separator. The vapor stream coming out of the separator, consisting mostly of reactants, is returned to the reactor after increasing its pressure. The liquid stream, mostly products, passes to a stripper for the final split of the products, reactants, inert, and byproduct. The difficulties in controlling the TE process can be attributed to the severe interactions among the temperature, pressure, and liquid level in the reactor. The liquid level strongly depends on the reaction rate because reaction products are liquids. The reaction rate, in turn, depends on the temperature (Arrhenius exponential dependence) and pressure (roughly third-order dependence). The pressure is also affected by the reaction rate because the reaction is accompanied by a significant reduction in the volume of the mixture (6 mol of reactants gives roughly 2 mol of products). Furthermore, the available manipulators for controlling the reactor are limited. Fresh flow rates of D and E must be held in a constant ratio for a given product mix, the fresh feed flow rate of A is negligible (both in mass and molar flow rate terms) when compared to the total

controlled variable

potential manipulated variables

y7 y8 y9 y12 y15 y17 y30 y42

u4, u6, u10, u11 u1, u2, u10, u11 u10, u12 u7, u11 u7, u8, u9 u1, u2, u7, u8, u11 u4, u6 u1, u2

reactor feed, and the fresh feed of C is not directly fed into the reactor. The available manipulated variables are summarized in Table 5. Note that u10 is the set point to the slave controller that adjusts the temperature of the reactor cooling water by manipulating the opening of the reactor cooling water valve.8 These manipulated variables are also indicated in Figure 4. The set of controlled variables selected among the 42 available measured variables is summarized in Table 6. Variable 42 is defined as the ratio G/H in the product stream. Similar sets of base controlled variables have been proposed in the literature.2,6,9,27-36 The linear time-invariant, state space description of the TE process was obtained by applying numerical differentiation to the routines supplied by Downs and Vogel11

x3 (t) ) Ax(t) + Bu(t) + Fd(t) y(t) ) Cx(t) + Du(t) + Ed(t)

(56)

Then the procedure described by McAvoy8 and based on the work by Arkun and Downs37 was followed in order to obtain the steady-state, transfer function matrices. Next, allowable pairings are selected (see Table 7) and transformed into integer constraints as shown in section 3 (referred to as heuristic 1). Further integer constraints were added to impose the requirements for self-consistent inventory control structures (heuristic 2) and pairing on positive RGA elements (heuristic 3). The weighting coefficients for violation of these constraints were set equal to 0.5, 10, and 0.5,

Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001 2087 Table 8. Best Control Structures for the TE Case Study y12 y15 y8

y17

y9

y7 y42 y30

1 2 3 4 5

u11 u11 u11 u11 u7

u7 u7 u7 u7 u8

u2 u1 u2 u2 u2

u8 u8 u8 u8 u11

u10 u10 u10 u10 u10

u4 u4 u4 u6 u4

u1 u2 u1 u1 u1

u6 u6 u3 u4 u6

6 7 8 9

u11 u11 u7 u7

u7 u7 u8 u8

u2 u1 u1 u2

u8 u8 u11 u11

u10 u10 u10 u10

u4 u6 u4 u6

u3 u2 u2 u1

u1 u4 u6 u4

violated heuristic

similar structures

3 (2) Luyben 3 (3) 1 (1), 3 (2) 3 (3) 3 (3) Scali and Cortonesi 1 (2), 3 (2) 3 (4) 3 (4) 3 (4)

respectively. Finally the optimization problem given by formulation (44) was solved (F ) 0.25) in 31.78 s (Ultra Sparc 360 MHz) using the GAMS interface with CPLEX22 and the index given by eq 11 as a measure of interaction. The results are summarized in Table 8. In Table 8, the nine best control structures obtained are given together with the heuristic(s) violated in each structure and the number of violations (in parentheses) plus a reference where a similar base control scheme is proposed. As can be seen from this table, the heuristic to prefer pairings on positive RGA elements is violated by all structures. The best structure uses product flow to directly control the production rate, inventories are controlled in the direction opposite to the flow, and C fresh feed flow rate is used to control the reactor pressure, purge stream flow rate to control the inert composition in the purge, and D fresh feed flow rate to control the product mix. This structure was previously derived by Luyben34 based solely on process insight. Structures 3, 4, 6, and 8 have minor differences, while structure 5 is a structure that has previously been proposed by Scali and Cortonesi.29 This structure uses the condenser cooling water flow rate to set the production rate. 5. Conclusions This paper presents a systematic methodology for selecting manipulated variables for regulatory control schemes. The selection of promising sets of manipulated variables is based on minimization of the overall interaction as expressed in terms of the steady-state RGA and the minimization of the sensitivity of the closedloop system to disturbances. A general methodology for incorporating heuristics as constraints to the problem is demonstrated. The main advantage of the method is that the plantwide nature of the problem is preserved because decisions related to different levels of the structure of the base control system are obtained simultaneously. The usefulness of the method is demonstrated using a double-effect evaporation process, the HDA process, and the TE challenge problem. Nomenclature g˜ ) element of the G ˜ matrix G ˆ ) submatrix of G corresponding to the selected manipulated and controlled variables G ˜ ) transpose of the inverse of the G ˆ matrix d ) vector of disturbances G ) transfer function matrix that relates the controlled variables with the potential manipulated variables g ) element of the G matrix Gd ) transfer function matrix that relates the controlled variables with the disturbance variables

i ) index over the set I I ) set of controlled variables j ) index over the set J J ) set of potential manipulated variables l ) index over the set I m ) index over the set M M ) set of disturbance variables n ) integer number P ) literal used to represent the logical statement s ) complex variable u ) vector of potential manipulated variables v ) integer slack variables w ) weight coefficient X ) matrix with integer elements y ) vector of measured variables Greek Letters δ ) Kronecker delta  ) positive variable used in formulation (24) η ) defined in eq 16 φ ) positive variable used in formulation (24) λ ) element of the RGA array Λ ) relative gain array µ ) positive variable used in formulation (24) σ ) element of the Σ matrix Σ ) matrix defined by eq 22 Ω ) large positive constant Subscripts d ) disturbance variables u ) control variables y ) measured variables

Literature Cited (1) Narraway, L. T.; Perkins, J. D. Selection of Process Control Structure Based on Linear Dynamic Economics. Ind. Eng. Chem. Res. 1993, 32, 2681. (2) Luyben, M. L.; Tyreus, B. D.; Luyben, W. L. Plantwide Control Design Procedure. AIChE J. 1997, 43, 3161. (3) Price, R. M.; Georgakis, C. Plantwide Regulatory Control Design Procedure Using a Tiered Framework. Ind. Eng. Chem. Res. 1993, 32, 2693. (4) Price, R. M.; Lyman, P. R.; Georgakis, C. Throughput Manipulation in Plantwide Control Structures. Ind. Eng. Chem. Res. 1994, 33, 1197. (5) Ponton, J. W.; Laing, D. M. A Hierarchical Approach to the Design of Process Control Systems. Trans. Inst. Chem. Eng. 1993, 71, 182. (6) Ng, C.; Stephanopoulos, G. Plant-wide control structures and strategies. Preprints of the 5th IFAC Symposium on Dynamics and Control of Process Systems (DYCOPS-5), Corfu, Greece, June 1998; Georgakis, C., Ed.; Elsevier. (7) Buckley, P. Techniques of Process Control; John Wiley & Sons: New York, 1964. (8) McAvoy, T. J. A methodology for screening level control structures in plantwide control systems. Comput. Chem. Eng. 1998, 22, 1543. (9) McAvoy, T. J. Synthesis of Plantwide Control Systems Using Optimization. Ind. Eng. Chem. Res. 1999, 38, 2984. (10) Tyreus, B. D. Dominant Variables for Partial Control. 1. Thermodynamic Method for their Identification. Ind. Eng. Chem. Res. 1999, 38, 1432. (11) Downs, J. J.; Vogel, E. F. A Plant-wide Industrial Process Control Problem. Comput. Chem. Eng. 1993, 17, 245. (12) Bristol, E. H. On a new measure of interaction for multivariable process control. IEEE Trans. Autom. Control 1966, 11, 133. (13) Raman, R.; Grossmann, I. E. Relation between MILP modelling and logical inference for chemical process synthesis. Comput. Chem. Eng. 1991, 15, 73. (14) Skogestad, S.; Postlethwaite, I. Multivariable Feedback Control; John Wiley & Sons: Chichester, U.K., 1996.

2088

Ind. Eng. Chem. Res., Vol. 40, No. 9, 2001

(15) Nett, C. N.; Manousiouthakis, V. Euclidean condition and block relative gainsconnections conjectures and clarifications. IEEE Trans. Autom. Control 1987, 32, 405. (16) Grosdidier, P.; Morari, M.; Holt, B. R. Closed-loop properties from steady-state gain information. Ind. Eng. Chem. Res. 1985, 24, 221. (17) Cavalier, T. M.; Pardalos, P. M.; Soyster, A. L. Modeling and integer programming techniques applied to propositional calculus. Computers & Operations Research 1990, 17, 561. (18) Grossmann, I.; Daichendt, M. M. New trends in optimization-based approaches to process synthesis. Comput. Chem. Eng. 1996, 20, 665. (19) Williams, H. P. Model Building in Mathematical Programming; John Wiley & Sons: Chichester, U.K., 1993. (20) Wolsey, L. A. Integer Programming; John Wiley & Sons: New York, 1998. (21) Narraway, L. T. Selection of Process Control Structures based on Economics. Ph.D. Thesis, Imperial College, University of London, London, 1992. (22) Brooke, A.; Kendrick, D.; Meeraus, A. GAMSsA User’s Guide; Scientific Press: Redwood City, CA, 1988. (23) Heath, J. A. Process Control Structure Selection Based on Economics. Ph.D. Thesis, Imperial College, University of London, London, 1997. (24) Douglas, J. M. Conceptual Design of Chemical Processes; McGraw-Hill: New York, 1988. (25) Cao, Y.; Rossiter, D. An input pre-screening technique for control structure selection. Comput. Chem. Eng. 1997, 21, 563. (26) Cao, Y.; Rossiter, D. Input selection and localized disturbance rejection. J. Process Control 1998, 8, 175. (27) McAvoy, T. J.; Ye, N. Base Control for the Tennessee Eastman Problem. Comput. Chem. Eng. 1994, 18, 383. (28) Lyman, P. R.; Georgakis, C. Plant-Wide Control of the Tennessee Eastman Problem. Comput. Chem. Eng. 1995, 19, 321.

(29) Scali, C.; Cortonesi, C. Control of the Tennessee Eastman benchmark: Performance versus integrity tradeoff. Proceedings of the 3rd European Control Conference, Rome, Italy, 1995; pp 3913-3918. (30) Walsh, S.; Heath, J. A.; Perkins, J. D. The Tennessee Eastman problem in teaching and research at Imperial College. Proceedings of the 3rd European Control Conference, Rome, Italy, 1995; European Union Control Association, pp 3925-3930. (31) Banerjee, A.; Arkun, Y. Control Configuration Design Applied to the Tennessee Eastman Plant-Wide Control Problem. Comput. Chem. Eng. 1995, 19, 453. (32) Ye, N.; McAvoy, T. J.; Kosanovich, K. A.; Piovoso, M. J. Optimal Averaging Level Control for the Tennessee Eastman Problem. Can. J. Chem. Eng. 1995, 73, 234. (33) Ricker, L. N. Decentralized Control of the Tennessee Eastman Challenge Process. J. Process Control 1996, 6, 205. (34) Luyben, W. L. Simple Regulatory Control of the Eastman Process. Ind. Eng. Chem. Res. 1996, 35, 3280. (35) McAvoy, T. J.; Ye, N.; Gang, C. Nonlinear Inferential Parallel Cascade Control. Ind. Eng. Chem. Res. 1996, 35, 130. (36) Tyreus, B. D. Dominant Variables for Partial Control. 2. Application to the Tennessee Eastman Challenge Process. Ind. Eng. Chem. Res. 1999, 38, 1444. (37) Arkun, Y.; Downs, J. A General Method to Calculate Input-Output Gains and the Relative Gain Array for Integrating Processes. Comput. Chem. Eng. 1990, 14, 1101.

Received for review April 20, 2000 Revised manuscript received August 7, 2000 Accepted February 5, 2001 IE000415T