Ind. Eng. Chem. Process Des. Dev., Vol. 18, No. 1, 1979
Langmuir, I., O.S.R.D. Report No. 865, Part I V (1942). Pearson, E., M.A.Sc. Thesis, University of Waterloo, Waterloo, Canada, 1977. Stechkina, I. B., Kirsh, A. A,, Fuchs. N. A,, (1969) KolloU. Zh.,31, No. 1, 121 (1969). Strauss, W., "Industrial Gas Cleaning", Pergamon Press, Australia, 1975 Sullivan, R. R., Hertel, K. L., J . Appl. Phys., 11, 761 (1940).
177
Whitby, K. T., et a / . , A . P . C . A . J . , 11, 11 (1961). Yeh, H. C., Liu, B. Y . H., Aerosol Sci., 5 , 191 (1974).
Received f o r reuieul March 30, 1978 Accepted August I, 1978
Multivariable Control System Design for an Industrial Distillation Column Bjorn D. Tyr6us AB Bofors, Nobel Kemi, S-690 20 Bofors, Sweden
A multivariable control system has been designed and tested on a 30-tray, 2-m diameter industrial distillation column. The design was performed with the Inverse Nyquist Array method which provided a rationale for variable pairing and compensator selection. The compensated system showed considerable improvement in dynamic behavior compared to a set of interacting controllers in a multi-loop configuration.
Introduction As material and energy prices increase, the incentive for improved process control also increases. For example, benefits from better distillation control can come from the possibility to operate closer to a quality constraint and thereby increase the yield of a more valuable product. Improvements in profitability can also come from decreased variability in the controlled qualities of a column, since nonlinear separation effects tend to increase the average energy consumption when the qualities are subject to variations. However, design of advanced process control systems is not always straightforward since the process often turns out to be multivariable and interacting in nature. A problem therefore arises in selecting a useful design method. Some of the more well known techniques have been listed in Table I. The choice between these methods is simplified when the required process models are considered. The last five methods work with state space representations while the first four are designed for transfer function models. It is usually quite difficult to obtain a manageable state space representation for a large distillation column while a small matrix of transfer functions can be derived from relatively simple plant tests. The following discussion is therefore confined to the first four techniques in Table I. Dominant interaction methods have been proposed by Bristol (1966) and Niederlinski (1971). The methods provide a rationale for variable pairing and in the latter case also controller tuning in a system of several interacting loops. The idea is to utilize the dominant interactions in the process while other interactions are ignored. In contrast to this the noninteraction methods seek to eliminate all interactions between the loops. This is accomplished with various decoupling techniques. Decouplers for distillation systems have been discussed by Luyben (1970) and Waller (1974). Even if specific examples have been given where noninteracting systems are superior to dominant interaction configurations, the indications are not clear as to when decouplers should be used in a general case. Furthermore there are several ways to decouple a system (i.e,, ideal decoupling, simplified decoupling, and partial decoupling) and there seems to be no general rule to guide in the selection of the appropriate variant. 0019-7882/79/1118-0177$01 .OO/O
Table I. List of Well Known Design Methods for Multivariable Systems Dominant Interaction Noninterac ting Inverse and Direct Nyquist Array Characteristic Loci Optimal Control Pole Placement Lyapunov Functions Incomplete State Feedback Multivariable Adaptive Control
The above difficulties are resolved with the other methods considered here. For example, both the Inverse Nyquist Array method and the Characteristic Loci technique give indications about the need for compensators. Furthermore the methods outline a procedure for selecting the appropriate compensator structure if compensators are needed. The choice between the two methods is probably quite arbitrary. In some situations as described by Hughes and Mallouppa (19761, both techniques can be used in combination to facilitate design. In the present study the Inverse Nyquist Array method is employed in the design of a control system for an industrial distillation column. The purpose of the investigation was to analyze an existing control system and to propose ways of improving it. The Inverse Nyquist technique turned out to be ideal for this type of exercise. Inverse Nyquist A r r a y Method The Inverse Nyquist Array method is outlined in full mathematical rigor elsewhere (e.g., Rosenbrock, 1974) and only a shorter description is given here in order to explain the main features of the design procedure. Consider a system in which G ( s ) is the open loop transfer function matrix and K ( s ) represents a combined controller and compensator matrix. Let Q ( s ) denote the product Q(s) = G(s)K(s)and let H ( s ) be the matrix of closed loop transfer functions. Straightforward manipulations then give
H
= (I
+ QF)-'Q
where I is the identity matrix and F is a diagonal matrix of feedback gains. The Laplace variable "s" has been dropped for clarity. By taking the inverse of both sides GZ 1978 American Chemical Society
178
Ind. Eng. Chem. Process Des Dev , Vol 18, No 1, 1979
of eq 1 and denoting H-‘ by H and Q by Q the following relations are obtained directly
Ei
+ QF)
= Q(Z
Ei=Q+F 6 1 1
=
(3)
+ f,
Qll
(2)
(4)
I? these equations h, Q, and f are elements of the matrices H , Q, and F , respectively. The closed loop characteristic equation of a multivariable system is defined as det(1 + QF). The existence of zeroes with positive real values in this equation determines the stability of the closed loops. These zeroes can be detected by counting the number of clockwise encirclements, N , that det(Z + QF)makes around the origin as the variable “s” encircles the right-hand side of the complex plane. For a stable system we have
N = -p 0 (5) where p o is the number of open loop poles in the right half plan?. From eq 2 it follows that det(Z + QF) = det(H)/ det(Q), which together with eq 5 gives (6) N H - NQ = -Po N H and are the number of clockwise_encirclements around the origin made by det(H) and det(Q), respectively. Evaluation of eq 6 is further simplified when the following theorem is used
Figure 1. Original control system.
depend on the controller settings in these other loops. This complicates the tuning of a multivariable interactive or partly interactive control system. Since the Inverse Nyquist Array method gives rise to an interactive or partly interactive system it is essential that the closed loop dynamic properties can be determined or estimated relatively easily. T o that end use is made of the following theorem which holds for a dominant system
k
LN, = Nz
(7)
1=1
If and only if k
1 4 > CIz1,l J=1
(8)
/#1
The theorem says that the sum of clockwise encirclements made by the diagonal elements in a complex matrix 2 equals the number of clockwise encirclements made by de@) provided the absolute values of the diagonal elements are greater than the sum of absolute values of the nondiagonal elements in the same row. Such a matrix is said to be diagonally row dominant. Equation 6 now becomes k
k
ENq1- 1=l LNh,
= Po
(9)
1=1
For a dominant system eq 9 makes it possible to assess stability of a closed loop multivariable system by inspection of the Nyquist plots of the diagonal elements in the inverse of Q. The number of encirclements around the origin made by Q,,forms the first term in eq 9, and the number of encirclements around -fl made by Q,,forms the second term in eq 9. This latter statement follows from eq 4. Diagonal dominance is assured for both Q and H if circles with radii k
along the paths of Q,,do not encompass neither the origin nor the point -f,. It is usually desirable to impose further criteria on a control system in addition to stability. However, in order to evaluate the dynamic properties of a particular loop in a multivariable control system the open loop transfer function between the controlled and manipulated variables in that loop is required. Since the other loops have to be closed it follows that the required transfer function, h,, will
Equation 10 says that the inverse of the open loop transfer function h, is located within a band around Q,,defined by the circles of radii r,. R e can thus use the same inverse Nyquist plots to assess stability as well as closed loop dynamic properties. The dynamic properties can be estimated with single input, single output criteria such as gain and phase margins. To summarize, we have indicated that without enforcing complete noninteraction the Inverse Nyquist Array method provides a means of treating a multivariable interacting system as a set of independent, single input, single output loops. The loop gains characterized by f,can be selected from a single set of plots to give total system stability and desirable closed loop behavior. Process Description The column studied in this work is a 30-tray, 2-m diameter industrial distillation column separating a partly vaporized multicomponent feed into two streams. The distillate which leaves the reflux drum as a vapor makes up roughly 25 mol YO of the feed. The bottom stream leaves the column as a liquid and is later separated into two salable products. Since the bottom product is more valuable than the distillate, the operational target is to maximize the light key component in the base. The permissible maximum depends on a concentration constraint for the lighter of the two salable products. Normally the maximum concentration for the light key in the base is less than 1 mol YO. Prior to this study the column was equipped with standard instrumentation. See Figure 1. A pressure controller sensing the pressure in the reflux drum manipulated the distillate flow. A temperature controller sensing a tray temperature close to the base manipulated the stream flow to the reboiler. The reflux flow was on flow control. The rate of condensation in the overhead could be altered via the set point of a pressure controller sensing the pressure in the condenser which contained a boiling refrigerant. The selection of the control tray near the base of the column was a compromise when two op-
Ind. Eng. Chem. Process Des. Dev., Vol. 18, No. 1, 1979
posing effects were considered. On the one hand, it was desirable to have the control tray a couple of plates up in the column since this increases the magnitude of the temperature response with respect to key component changes. On the other hand, it was desirable to keep the sensor in the base since the temperature is meant to reflect composition changes in the bottom stream. One of the reasons for studying the column was some dynamic problems following intermittent changes of, for example, reflux flow in order to affect the level in the reflux drum. The level could be changed either with the condensor load or with reflux flow. However, both these variables affected the column pressure, which in turn changed the tray temperatures. The temperature controller responded to this change and affected both the reflux drum level and the column pressure. The interactive oscillations that resulted made tight temperature control difficult and the column had to be operated with low concentrations of light key component in the base. This safety margin in composition prevented the salable products to go off specification after a disturbance but inevitably meant lower recovery of one of the valuable products. A Process Model For reasons stated previously, a transfer function model was considered for design purposes. The size of the process matrix was determined by the number of feedback loops included. The column already had two interacting loops but it seemed important to add a third loop, namely one including automatic control of the reflux drum level. A third manipulative variable was therefore needed. Reflux flow was selected because steady-state simulations had shown that the condensor should operate at the lowest permissible pressure at all times. This excluded the condensor load from being a candidate for control. In order to obtain a model of the column, pulses in the three manipulative variables were injected. The resulting responses on the output variables were analyzed with techniques described by Luyben (1973). By fitting loworder transfer functions to the frequency response data obtained from the pulse test analysis a matrix of process transfer functions G(s) was obtained
The values of the G matrix are given in Table 11. As seen from this table, several of the transfer functions contained integrators which partly explained some of the control difficulties experienced in practice. Design Procedure T h e first step in the design procedure is to perform variable pairing and to check the need for compensators. This is done by plotting the Nyquist curves of the diagonal elements- in the inverted process matrix. Therefore, the inverse Q = G was evaluated numerically for a range of frequencies wmin < w < cornax.At each frequency the sum k
ri =
C IBijl
j=1 j#i
was computed. Some of the circles with radius rLwere drawn with centers along the plots of Qii. Figure 2 shows these plots for the process matrix described by eq 11. From Figure 2 it is evident that the system was not diagonally dominant since some of the r L circles include the origin. The easiest way to attempt to improve on dom-
179
Table 11. Open Loop Transfer Functions (Time Constants in Minutes) e-2.5s
g,,= 1.318 (20s e- 4s g,, = 0.333
+
1)
S
e- 11s
-
g,, = -0.11
S
g2,= 0.038
(182s + 1) ( 2 7 s + 1)(1Os t 1)(6.5s + 1)
0.36 g12 = S
e- 38
g,, = -0.12
s(14s + 1)(8s + 1) 0.313
g31 =
+
s(29s
1)(17s+ 1) 2.222
g32 = -
s(l0s t 1)(6.7s + 1) 1.563
g33 =
s(17s + 1)
Figure 2. Inverse Nyquist array for Q = G.
inance is to check the pairing of controlled and manipulated variables. A different combination of these variables corresponding to an interchange of the columns in G(s) might lead to a dominant system. For the distillation tower considered here this was not the case and the variable pairing defined by eq 11 gave a system which was closest to dominant. The next step in the design procedure is to achieve dominance. In the method outlined by Rosenbrock (1974) this is done with a computerized design package and ipteractive computing. The computer displays the whole Q matrix and from inspection of this matrix various compensator matrices are tried. When a design package is not available a simplified procedure can be adopted. This procedure is based on a study of the open loop matrix G(s). For example, inspection of eq 11 showed the elements g,, and g,, could be reduced in relation to g,, if a
180
Ind. Eng. Chem. Process Des. Dev., Vol. 18, No. 1, 1979
Id
m'
Figure 4. Block diagram of closed loop, multivariable control system with interaction compensators K1 and K2.
is the process transfer function matrix, the K matrices are compensators, and the B matrix represents the feedback controllers. From Figure 4 it follows that x* = Kl(x d ) (14) x = Gm (15)
=&
+
m = K2m* x* = K1GK2m* Kld
(16)
+
Figure 3. Inverse Nyquist array for Q = K1GK2.
multiple of row 2 was subtracted from row 1. In matrix notation this is described as 1 -0.92e'4s Q=KIG=[i
0 :]G
;: )
(12)
A T - 0.92e-"AP
x*=(iE)'-..=(
.
The effect of the premultiplier K 1 is to improve dominance in the temperature loop. Similarly, dominance for the level loop was improved by making the elements g,, and g32 small in relation to g33. A postmultiplier K2 was needed in order to add multiples of column 3 to columns 1 and 2 in the G matrix. Q = K,GK, = K , G
(17) Comparing eq 17 with eq 13 it is evident that Figure 4 represents the compensator structure for the distillation column. The physical intepretations of the compensators K 1 and K 2 are given by eq 14 and 16. For example, eq 14 gives
c-
"I
1 0 0 1 -0.2 1 . 4 2 ( 1 7 ~+ 1) 1 29s + 1 (10s + 1)(6.7s + 1)
(I3)
The new matrix Q = KlGK2 was inverted numerically to give the Inverse Nyquist Array shown in Figure 3. Diagonal dominance had been achieved a t this stage. The last step in the design procedure is to determine suitable controller gains. One of the attractive features of the design method is that plots of the type shown in Figure 3 can be used to that end. Suitable values for the feedback gains f, are selected and it is then realized that the loop gains will remain the same if the feedback gains are set to unity and the controller gains are given the values f,. The following points must, however, be considered in selecting the f,. First, the f, must be located outside the r, circles to assure dominance for H. Secondly, for an open loop stable system the f, must be located to the right of the intersection between the Nyquist plot and the negative real axis. This will assure closed-loop stability according to eq 9. Finally, the f,should be located such that desirable closed-loop dynamic properties are obtained. These properties can be estimatedfrom the plots since the inverse of the open-loop transfer function is located within the band swept out by the r, circles according to eq 10. The selected values for the f, are shown in Figure 3. These gains provided good starting points for later controller tuning in the plant. Figure 3 also indicated that integral action could be used in the temperature and pressure loops but not in the level loop. Physical Realization Consider the closed loop multivariable control system shown as a block diagram in Figure 4. In this system G(s)
(18)
which is nothing more than pressure compensation of the control tray temperature. Similarly, eq 16 gives A QB
K,m* =
L+ i -0.2m,*
m1* m2* 1.42(17s + l ) m , *
(10s + 1)(6.7s + 1)
+ m,*
(19)
which shows how the outputs from the three feedback controllers should be connected to the control valves. In order to facilitate the implementation of eq 18 and 19 in standard analog hardware the following approximations were made 0.92 0.92e-4s = 5s 1
+
1.42(17s + 1) 1.42 z (10s 1)(6.7s + 1) 5s + 1
+ 0.2 --- 0.2 29s + 1 20s + 1
(21) (22)
Since these first-order lags are fairly rough approximations of the designed compensators, it was decided to simulate the closed loop system. The process transfer functions, the first-order lags, and the set of feedback controllers were therefore translated into difference-differential equations which were numerically integrated in time. Several different disturbances ( d in Figure 4) were injected and the output results were plotted. An example of a computer-simulated response of the controlled variables after a load change in pressure is shown in Figures 5 through 7. Also shown in these figures are the responses from a well-tuned multiloop interactive control system subject to the same pressure disturbance. The beneficial effects of the simple compensators are clearly revealed.
Ind. Eng. Chem. Process Des. Dev., Vol. 18, No. 1, 1979
181
Figure 5. Simulated control tray temperature response after load change in column pressure: -, compensated control system; ---, uncompensated (interactive) control system.
--_ n
_____
'
I
8.
- W E (br)
Figure 6. Simulated pressure response after load change in column pressure.
Figure 8. Physical realization of final control system. -
I:-L-
I
,, 1 br -
~
I.
T ME (hr)
Figure 7. Simulated change in reflux drum level after load change in column pressure.
T h e control system shown in Figure 1 was augmented with standard analog hardware (commercial adderssubtractors and lead-lags) to give the final control system shown in Figure 8. The square root extractor shown in Figure 8 was inserted in order to make the reflux flow respond linearly to the flow set point changes. Results After initial instrument and wiring checkout the system was commissioned in three steps. First the level and pressure loops were closed with the temperature controller in manual and pressure compensation disconnected. When stable operation was obtained the temperature controller and the pressure compensation were brought in. Finally the temperature controller was tuned to give the best possible control of the light key component in the base. Improvements over the old system were mainly achieved through a sixfold increase in the temperature controller gain made possible with the new system. T h e overall improvements in system dynamics are shown in Figure 9 reflecting actual plant data. Figure 9 shows the flows of distillate ( D ) and salable products before and after the new system was implemented. The salable products are obtained from a downstream distillation unit which is fed directly with the bottom flow of our column. Variations in the bottom flow rate are therefore reflected in the flow rates of these products. T h e new control system permits closer control of the base temperature. This allows the operators t o run the unit a t high key component concentrations and still not overshoot the target and cause off specification products. In addition, less heavy components are boiled over top
~
1-
I
Figure 9. Strip charts with product flow rates before (top) and after (bottom) the new control system was implemented.
which lowers the top temperature. This gives less condensation in the condenser when the refrigerant temperature is constant (constant pressure). The net result is a lower reflux flow and reboiler heat duty with a consequent saving in energy. The total benefits of the control system have been estimated to give a pay back period of less than 3 months. Conclusions The Inverse Nyquist Array method has proved useful in a multivariable control system study for an industrial distillation column. The method gave information about variable pairing and the need for interaction compensators. Since compensators were needed, the method also provided a design objective which could be met with standard analog hardware. Considerable dynamic improvements were achieved over a multiloop interacting system. Nomenclature B , matrix of controllers AD, change in distillate flow rate d , vector of disturbances FC, flow controller FL,,feed rate of liquid Fv, feed rate of vapor
182
Ind. Eng.
Chem. Process Des. Dev.,
Vol. 18, No. 1, 1979
f,, feedback gain
po, number of open loop poles in the right half side of the
G, inverse of G g element of G matrix of closed-loop transfer functions H, inverse of H h,, open-loop transfer function in loop '7" when the other loops are closed h,,? element of A I , identity matrix K , matrix of dynamic compensators k , number of feedback loops in the system defined by G AI,, change in reflux drum level LC, level controller LI, level indicator m, manipulative variable m*, vector of controller outputs N, number of encirclements around the origin N,, number of encirclements around the origin made by an element in a matrix NZ, number of encirclements around the origin made by the determinant of the complex matrix 2 .QH, number of encirclements around the origin made by det(H) NQ,number of encirclements around the origin made by det(Q) N h , ,number of encirclements around the origin made by h,, N number of encirclements around the origin made by GLL Ap,' change in column pressure PC, pressure controller P T , pressure transmitter
AQB, change in reboiler load
C, matrix of process transfer functions
4,'
complex plane
Q, product between process matrix and compensator matrices Q, inverse of Q ,
4..element of Q A%, change in reflux flow r i , suqi of absolute values of nondiagonal elements in a row
of
Q
s, Laplace variable
A T , change in control tray temperature T C , temperature controller T T , temperature transmitter x, vector of controlled variables x*, vector of compensated controlled variables x*set,set points for x * 2, complex matrix z,,, elements in 2 L i t e r a t u r e Cited Bristol, E. H., I€€€ Trans. Autom. Control, AC-11, 133 (1966). Hughes, F. M., Mallouppa, A,, Automatica, 12, 201-210 (1976). Luyben, W. L., AIChE J., 16 (2),198 (1970). Luyben, W. L., "Process Modeling, Simuhtion and Control for Chemical Engineers", McGraw-Hill, New York, N.Y., 1973. Niederlinski, A,, Automatica, 7, 691 (1971). Rosenbrock, H. H., "Computer-Aided Control System Design", Academic Press, London, 1974. Waller (Toijala), K. V. T., AIChE J . , 20 (3), 592 (1974).
Received f o r review May 19, 1978 Accepted September 8, 1978
COMMUNICATIONS The Number of Roots in the NRTL and LEMF Equations and the Effect on Their Performance The number of roots in the NRTL and LEMF equations that can be obtained from a set of binary vapor liquid equilibrium data is investigated. When the two infinite dilution activity coefficients are used, it is shown that for positive deviations from Raoult's law, one pair of roots exists for the NRTL equation and up to three for the LEMF equation; for negative deviations, the reverse occurs. When all (T-P-X- V ) data are regressed, a multiplicity of roots was observed for the LEMF equation for y's both larger and smaller than 1. For the NRTL equation and for y > 1.0, one set of roots was observed at a = 0.2 and 0.3 and two at LY = 0.47; for y < 1.0, three sets of roots were found. For both equations, starting values of (0.0 and 0.0) in the regression subroutine lead to the smallest in absolute value roots and gave the best fit of the data. Roots obtained by using the 7,-and y2- values correlate the binary data of this study with accuracy comparable with that obtained from regressing all data; for the LEMF equation, however, and for positive deviations from Raoult's law, the other two pairs predict partial miscibility.
Introduction The Wilson, NRTL, and LEMF equations have found extensive use in the correlation and prediction of vapor liquid equilibrium data (Wilson, 1964; Orye and Prausnitz, 1965; Renon and Prausnitz, 1968; Marina and Tassios, 1973). For this purpose available binary experimental data are regressed through a nonlinear regression subroutine t o obtain the parameters in these equations that best fit the data. When the Wilson equation is used, a parametric substitution shows that one pair of roots is obtained for positive deviations from Raoult's law and up to three for negative deviations (Aristovich e t al., 1969; Silverman and Tassios, 1977). Multiplicity of roots in the correlations of mutual solubility data has been discussed for the NRTL
equation by Mattelin and Verhoeye (1975) and by Heidemann and Mandhane (1973). I t is the purpose of this paper to investigate the number of pairs of roots that can be obtained from a set of binary VLE data for the NRTL and LEMF equations and the effect of the multiple roots, whenever they exist, on the accuracy of correlating these binary data. References for the data used in this paper are presented in Table I. Activity coefficients were calculated from the experimental X, Y , P, T data according to the method of Prausnitz e t al. (1967). The NRTL a n d L E M F Equations Renon and Prausnitz (1968), using Wilson's (1964) concept of local mole fractions and Scott's (1956) two1978 American Chemical Society
