Ind. Eng. Chem. Res. 2002, 41, 4585-4590
4585
PROCESS DESIGN AND CONTROL Analyzing Dynamic Interaction of Control Loops in the Time Domain F. Michiel Meeuse† Department of Chemical Technology, Delft University of Technology, Julianalaan 136, 2628 BL Delft, The Netherlands
Adrie E. M. Huesman* Department of Mechanical Engineering, Delft University of Technology, Mekelweg 2, 2628 CD Delft, The Netherlands
Process control structures are still mainly based on multiple interacting single-input singleoutput loops. The well-known relative gain array (RGA) is often used as an indication of the amount of interaction in the system. The RGA is limited in the sense that it is based only on static responses. We extend the RGA to compare the entire time domain open- and closed-loop step responses as interaction measurements. The closed-loop response is based on the best achievable control performance using internal model control. The results are shown graphically by regarding the open- and closed-loop responses of each pairing as a parametric representation. For each possible pairing, a graph is produced and placed in an array. This dynamic array is completely consistent with the static RGA and a useful extension in the case that the static RGA might prefer the wrong pairing (because the dynamics have been excluded) or does not prefer any pairing (when all RGA elements have more or less the same value). 1. Introduction Almost all chemical processes are multivariable by nature. Despite the fast developments in multivariable control, most control systems are still based on singlevariable control. This generally leads to a set of interacting single-input single-output (SISO) control loops. Bristol1 introduced the well-known relative gain array (RGA) to quantify steady-state interaction. For a system with n manipulated and n controlled variables, the RGA is a n × n matrix consisting of the ratios of the corresponding open- and closed-loop gains. The RGA can be used to determine which pairing of manipulated and controlled variables should be used in order to minimize interaction. For an overview of the RGA, the reader is referred to textbooks such as those by Stephanopoulos2 and Skogestad and Postlethwaite.3 The RGA is based on static gains; in other words, the process dynamics are neglected. Especially for systems containing time delays, right-half-plane zeros, or large differences in time constants, a dynamic analysis of interaction can be very important. On the basis of the steady-state analysis, one would prefer a certain pairing, whereas the inclusion of the dynamics would lead to a different pairing. It should be noted, however, that an acceptable static RGA remains a necessary condition. * To whom correspondence should be addressed. E-mail:
[email protected]. Phone: +31 (0) 15 2788131. Fax: +31 (0)15 2789387. † E-mail:
[email protected]. Phone: +31 (0) 15 2784466. Fax: +31 (0) 15 2784452.
Therefore, a dynamic analysis is only necessary for systems with an acceptable static RGA. Another motivation for a dynamic interaction analysis is that sometimes no preference for the pairing can be determined based on steady-state analysis only (e.g., a 2 × 2 system with all RGA elements close to 0.5). Dynamic analysis might reveal which pairing should be preferred. Witcher and McAvoy4 proposed to extend the RGA to the frequency domain. Gagnepain and Seborg5 modified the dynamic approach of Witcher and McAvoy4 by averaging the result over the desired frequency range. The application of this is rather limited because now all dynamic effects are averaged. Tung and Edgar6 calculated the dynamic RGA out of a state-space description of the process. However, their method is not essentially different from the method proposed by Witcher and McAvoy.4 The methods discussed above are all based on perfect control. In the steady-state this might be realizable, but dynamic perfect control is, in general, not possible. Huang et al.7 define a generalized dynamic relative gain (GDRG) using more realistic control based on the internal model control (IMC) principle.8 However, the GDRG is also based on averaging of the dynamic response over the frequency domain so the applicability is limited. In this paper we propose a new method for the analysis of dynamic interaction. The open- and closedloop step responses are compared over a finite time horizon. The closed-loop responses are based on practical realizable controllers. The analysis is done in the time domain because of the large advantages it offers for the interpretation.
10.1021/ie000886f CCC: $22.00 © 2002 American Chemical Society Published on Web 08/02/2002
4586
Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002
This paper is organized as follows. First a review will be given of the dynamic RGA. Then our new approach will be presented, followed by a number of examples. Finally some conclusions will be drawn about the applicability of this new approach. 2. Review of Dynamic RGA Consider a square system with n manipulated variables u and n controlled variables y. The RGA element λij that is the result of pairing of the controlled variable yi with the manipulated variable uj is defined as
λij )
( ) ( ) ∂yi ∂uj
ul*j constant
∂yi ∂uj
yk*i constant
(1)
The numerator in expression (1) is the gain of the loop ij, with all other manipulated variables kept constant (open-loop gain). The denominator is the gain of the same loop, but now with control applied on all other loops keeping all other controlled variables constant (closed-loop gain). The complete RGA matrix Λ can be computed from the transfer function matrix G:
Λ(s) ) G(s) X G(s)-1
T
(2)
where X denotes element by element multiplication and s the Laplace variable. When we insert s ) 0 into eq 2, we obtain the static gains and hence the original RGA as proposed by Bristol.1 Witcher and MacAvoy4 were the first to extend this to the frequency domain. The frequency-dependent RGA is simply calculated by inserting s ) jω into eq 2. There are, however, some arguments against the use of this frequency-dependent RGA: (i) The physical interpretation is much less transparent than the steady-state interpretation. (ii) It relies on unrealizable controllers. (iii) Chemical engineers are not used to working in the frequency domain. All of these arguments will be discussed below. Calculation of the frequency-dependent RGA according to eq 2 will lead to a matrix with complex elements. These complex elements can be represented by a magnitude and a phase. The interpretation of these complex elements is not very transparent. Witcher and McAvoy4 claim that “the magnitude can be interpreted in exactly the same way as the standard relative gain array”. However, the standard RGA can have negative elements while the magnitude is always nonnegative. Moreover, the property that the sum of elements in a row or column equals 1 does not hold for the magnitudes. Seider et al.9 try to get around this by using a slightly modified expression for the dynamic RGA elements. For one element (let us suppose the 1,1 element), the following expression is used:
λ11 ) sign (λ11(0))|λ11(jω)|
(3)
Then the 1,2 element is calculated with
λ12 ) 1 - λ11
(4)
Figure 1. Structure of IMC. The plant is indicated by P, P ˆ is a plant model, and Q an approximation of the inverse of P ˆ.
However, different results will be obtained when the 1,2 element is used in eq 3 rather than the 1,1 element. So, this frequency-dependent RGA is not uniquely defined. The calculation of the dynamics of any closed-loop response relies on the selected controller. When the RGA is calculated according to eq 2, the implicit assumption is made that the controller is the inverse of the process. However, this will lead to nonproper controllers that can also be unstable when right-half-plane zeros are present or unrealizable when time delays are present.10 The frequency domain stems from the systems and control community. It has, however, never been embraced by chemical engineers, especially not by those working outside academia. It is also our personal experience that in the industrial process control practice the frequency domain is completely omitted. Especially, the RGA is intended to be used by the industrial chemical engineers and by process control engineers. So, a more accessible dynamic interaction method is required. The physical relevance of the frequency domain is also limited because disturbances are generally not sine-shaped. The same holds for set-point changes. 3. Proposed Approach The basic idea of our approach is to compare two different step responses of a specified manipulated variable on a specified controlled variable. The first step response is the open-loop response, so all other manipulated variables are kept constant. The second is the closed-loop response, so all other control loops are put on automatic. Both responses will be calculated in the Laplace domain and then transformed to the time domain. We will consider the general square system described by
y ) G‚u
(5)
with y the vector of outputs, u the vector of inputs, and G the transfer matrix. Open-Loop Response. The open-loop step response of yi to a change in uj is then given by
yij,ol(t) ) L
-1
[G 1s] ij
(6)
where L denotes the Laplace transform. Closed-Loop Response. For the closed-loop response, a controller needs to be designed for all other control loops. Obviously, this response depends on the selected control algorithm. Our approach uses IMC as presented by Garcia and Morari.11 Figure 1 shows the
Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002 4587
Figure 2. Closed-loop structure containing the internal model controller.
structure of IMC. The IMC design consists of two steps: (1) decomposition of the plant model, P ˆ , into an invertible part and a part that cannot be inverted; (2) design of a filter F such that the controller becomes proper. The plant model is first decomposed according to ref 11:
P ˆ )P ˆ +1 P ˆ +2 P ˆ-
(7)
where P ˆ +1 is selected to make P ˆ -1P ˆ + realizable and P ˆ +2 is selected to make P ˆ -1P ˆ +1P ˆ +2 stable. The following two steps are now required: (i) The transfer matrix P ˆ +1 is selected to be diagonal, as described by Morari10 and Garcia and Morari.11 (ii) For the matrix P ˆ +2, an inner/outer factorization is applied. The filter is required to make the system proper. The typical filter used in IMC is given by
F)
1 (s + 1)r
(8)
where r is selected large enough to make FP ˆ --1 proper. The parameter is a tuning parameter. This parameter can be used to trade off performance vs robustness. The transfer matrix Q now becomes FP ˆ --1. Another approach is to use the concept of perfect control. This is done by Witcher and McAvoy,4 Tung and Edgar,6 and Gagnepain and Seborg.5 Perfect control is obtained when the complementary sensitivity function equals 1.10 In terms of IMC (see Figure 1), this means that if P ˆ ) P (perfect model), Q equals P ˆ -1. In reality, this is not possible for the following reasons: (i) The inverse is not realizable for systems with time delays because the controller then contains prediction terms. (ii) The inverse contains right-half-plane poles when the plant contains right-half-plane zeros; this will lead to an unstable controller. (iii) The inverse will not be proper in the general case that the plant is strictly proper. It should be noted that while calculating λij the plant P equals G without the ith row and the jth column. Figure 2 shows the complete signal diagram for the closed-loop system. From this figure it can easily be seen
Figure 3. Typical graph. The open- and closed-loop responses are given by eqs 10 and 11. Starting at the origin, the circles indicate t ) 0-3 and 10 time units.
that the closed-loop step response of yi to a change in uj is given by
yij,cl(t) ) L
-1
[(G
ij
- Gi,l*j QGk*i,j )
1 s
]
(9)
Dynamic RGA. In the traditional RGA, the openloop is divided by a closed-loop gain. A logical extension would be to divide the open-loop time response by the closed-loop time response. However, in the case of time delays, this might lead to division by zero. We propose to present the data as a parametric representation. That is, for every time instant, a point is plotted with the closed-loop step response as an x coordinate and the open-loop step response as a y coordinate. It should be clear that if the parametric representation coincides with the line y ) x, both responses are the same. Suppose that the open- and closed-loop responses are given by
yol(t) ) L ycl(t) ) L
-1
[(
-1
[( s +1 1) 1s]
(10)
)]
1 1 1 + s + 1 s2 + 3s + 1 s
(11)
Then Figure 3 shows the graph. The graph runs from the point t ) 0 to t ) ∞. All graphs are placed in an array; so, for a 2 × 2 system the result would be four graphs placed in an array of two rows and two columns. In the rest of this paper, this array of graphs will be referred to as the dynamic RGA. Properties of the Dynamic RGA. We will now present some properties of the dynamic RGA. (i) The graph always starts (t ) 0) in the origin. The point of the graph corresponding to t ) 0 is (ycl(0), yol(0)). For smooth systems, ycl(0) ) yol(0) ) 0 and hence the starting point of the graph is the origin. (ii) The static RGA values are given by dividing the final value of the open-loop step response by the final value of the closed-loop step response. The end point of the graph, corresponding to t f ∞, is (ycl(∞), yol(∞)). These are the open- and closed-loop
4588
Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002
gains, respectively, and they do not depend on the applied control algorithm. So, the static RGA is simply the open-loop gain divided by the closed-loop gain. (iii) For stable systems, the graph is finite. This is an obvious result from the previous property. (iv) The property that the sum of the elements in a row or column equals 1 is lost. For any matrix M, given by
M ) K X K-1
T
Example 1. The plant P1 is given by
P1(s) )
(
e-s/(s + 1) 1/(s + 1) -1/(s + 1) e-2s/(s + 1)
(12)
4. Examples As explained in the Introduction, there are two reasons for extending the static RGA to a dynamic RGA: (1) In some cases the static RGA prefers the wrong pairing (because the dynamics are excluded from the analysis). (2) In other cases the static RGA does not prefer any pairing (when all elements have more or less the same value). In this chapter two examples are given. In the first example, the static RGA does not prefer any pairing, and in the second example, the static RGA prefers the wrong pairing.
(13)
The static RGA can be calculated with eq 2:
ΛP1 )
the sum of the elements in a row or column equals 1. Equation 12 is the basis of the traditional RGA and of the frequency-dependent RGA approaches. So, for those approaches the sum of the elements in a row or column equals 1. Note, however, that in the frequency-dependent approaches often the magnitude of the elements is selected. The property does not apply for the magnitude of the elements. In our dynamic RGA, we do not perform the calculation described by eq 12; hence, we lose this property. (v) For a completely decoupled system, the graph coincides with the line y ) x. For a completely decoupled system, ycl(t) ) yol(t), ∀ t ∈ R+. So, the parametric representation coincides with the line y ) x. (vi) If the graph intersects the line y ) ax, the RGA at that time instant has the value a. This results from the fact that at any time instant the x coordinate of the graph is the closed-loop response and the y coordinate the open-loop response. Interpretation of the Graphs. As was already hinted before, the most desirable situation, from an interaction point of view, is a completely decoupled system. This means that the graph coincides with the line y ) x. Deviations from this line are a measure for the amount of interaction. As a matter a fact, it is possible to “enclose” a graph by two lines through the origin. For example, in Figure 3 the graph can be enclosed by the lines y ) 0 and y ) 0.6x. So, the dynamic RGA covers the interval [0, 0.6]. This interval should be compared with the “ideal” value of 1. If the graph initially moves over the x axis, it implies that the closed- and open-loop responses have different time delays. This represents an unfavorable pairing because opening other loops might result in a stability problem. Something similar happens when the graph moves into the second or fourth quadrant. This means that the open- and closed-loop responses have different signs. So, again this represents an unfavorable pairing because of a potential stability problem.
)
(
0.5 0.5 0.5 0.5
)
(14)
So, the static RGA does not prefer any pairing. The proposed dynamic RGA can be developed as follows. Equation 13 gives the various open-loop responses; for example, the open-loop response of the 1,1 element is given by
y11,ol(t) ) L
-1
[( ) ] e-s 1 s+1 s
(15)
The closed-loop step response of the same element can be derived by substituting the correct transfer functions in eq 9:
Note that Q ) FP ˆ _-1, where the filter, F, to make Q proper is 1/(0.1s + 1). Simplifying gives
y11,cl(t) ) L
-1
[(
)]
10 e-s 1 + s + 1 s2 + 11s + 10 s
(17)
Repeating the same exercise for the other elements gives
y12,cl(t) ) L
-1
y21,cl(t) ) L
-1
y22,cl(t) ) L
-1
[( [( [(
)] )] )]
1 1 10e-3s + 2 s + 1 s + 11s + 10 s
(18)
-1 1 -10e-3s + 2 s + 1 s + 11s + 10 s
(19)
e-2s 1 10 + s + 1 s2 + 11s + 10 s
(20)
Equations 13 and 17-20 are used to construct Figure 4. Figure 4 confirms that the static RGA does not prefer any pairing. All graphs except the one of the 2,1 element end at (2, 1). So, the static RGA values of these elements is 1/2 ) 0.5. The graph of the 2,1 element ends at (-2, -1), so the static RGA value of this element also equals -1/-2 ) 0.5. However, Figure 4 reveals that the diagonal elements approach this static RGA value quite differently than the off-diagonal elements. The graphs of the diagonal elements can be enclosed by the lines y ) 0 and y ) 0.5x, so the dynamic RGA of the diagonal elements covers the interval [0, 0.5]. The dynamic RGA of the off-diagonal elements moves over the interval [0.5, 1.0]. Note that the interval of the diagonal elements is further from the ideal value 1.0 than the interval covered by the off-diagonal elements. Another thing is that the diagonal elements move away from y ) x right from the start over the line y ) 0. The off-diagonal elements also move away from y ) x, but this happens
Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002 4589
Figure 4. Dynamic RGA of plant P1. Starting at the origin, the circles indicate t ) 0-3 and 10 time units.
Figure 5. Dynamic RGA of plant P2. Starting at the origin, the circles indicate t ) 0, 5, 10, 20, and 100 time units.
only after three time units. So, the open- and closed-loop responses have different time delays. In other words, the static RGA does not prefer any pairing, but the dynamic RGA clearly prefers the off-diagonal pairing. Example 2. Plant P2 is taken from work by Seider et al.:9
The RGA value of the diagonal elements moves over the interval [0, 0.667], while the RGA value of the offdiagonal elements covers the interval [0.333, 1.0]. Figure 5 also shows that the diagonal elements move away from y ) x over the line y ) 0 during the first five time units. The off-diagonal elements only move away from y ) x after 10 time units. So, the static RGA prefers diagonal pairing, but the dynamic RGA clearly prefers off-diagonal pairing.
P2(s) )
(
2.5e-5s/(15s + 1)(2s + 1) 5/(4s + 1)
)
-4e-5s/(20s + 1) (21)
1/(3s + 1)
The static RGA can be calculated with eq 2:
ΛP2 )
(
0.667 0.333 0.333 0.667
)
(22)
So, the static RGA prefers diagonal pairing. Substitution of the correct transfer functions in eq 9 gives
y11,cl(t) ) L
-1
y12,cl(t) ) L
-1
[(
2.5e-5s + (15s + 1)(2s + 1) 1 250s + 12.5 (23) 12s3 + 127s2 + 71s + 10 s
)]
[(
5 + 4s + 1 (300s + 100)e-10s
)]
1 (24) 600s + 6370s + 3737s + 371s + 10 s 4
y21,cl(t) ) L
3
-1
[(
2
1 + 3s + 1 (80s + 20)e-10s
)]
1 (25) 600s + 6370s + 3737s + 371s + 10 s 4
y22,cl(t) ) L
-1
3
[(
-4e-5s 20s + 1
2
)]
1 6000s2 + 3400s + 200 (26) 4 3 2 12s + 247s + 1341s + 720s + 100 s Equations 21 and 23-26 are used to construct Figure 5.
5. Conclusions In this paper the RGA is seen as a structured way to compare open- and closed-loop responses for different pairings. The static RGA was extended to include the dynamics of the plant based on the following assumptions: (i) A unit step response is sufficient to characterize the dynamics. (ii) The closed-loop response should approach perfect control in a realistic way. This is achieved by IMC. Rather than taking the quotient of the open- and closed-loop responses, as was done in the static RGA, the open- and closed-loop responses are compared by means of a graph. In such a graph every time instant is a point, with the closed-loop step response being the x coordinate and the open-loop step response being the y coordinate. For each possible pairing, a graph is produced and placed in an array. This dynamic array is completely consistent with the static RGA (each value of the static RGA agrees with one point in each graph of the dynamic RGA). From the two examples in this paper, it is clear that the approach described above works. The dynamic RGA provides extra and useful information in the following cases: (i) The static RGA prefers the wrong pairing (for systems containing time delays, right-half plane zeros, or large differences in time constants). (ii) The static RGA does not prefer any pairing (when all elements have more or less the same value). For 2 × 2 systems, the approach is straightforward in the sense that the calculations (the IMC design) can be done manually. However, for larger systems it is advisable to do the calculations by dedicated software (e.g., Matlab). Because the dynamic RGA is an extension
4590
Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002
of the static RGA, an acceptable static RGA remains a necessary condition for an acceptable degree of interaction. So, in other words, it is a good idea to check the static RGA first. Literature Cited (1) Bristol, E. H. On a new measure of interaction for multivariable process control. IEEE Trans. Autom. Control 1966, 133134. (2) Stephanopoulos, G. Chemical process control; PrenticeHall: London, 1984. (3) Skogestad, S.; Postlethwaite, I. Multivariable feedback control; Wiley: Chichester, U.K., 1996. (4) Witcher, M. F.; McAvoy, T. J. Interacting control systems: steady state and dynamic measurement of interaction. ISA Trans. 1977, 16, 35-41. (5) Tung, L. S.; Edgar, T. F. Analysis of control-output interactions in dynamic systems. AIChE J. 1981, 27, 690-693.
(6) Gagnepain, J. P.; Seborg, D. E. Analysis of process interactions with applications to multiloop control system design. Ind. Eng. Chem. Process Des. Dev. 1982, 21, 5-11. (7) Huang, H. P.; Ohshima, M.; Hashimoto, I. Dynamic interaction and multiloop control system design. J. Process Control 1982, 4, 15-27. (8) Garcia, C. E.; Morari, M. Internal model controlsI. Ind. Eng. Chem. Process Des. Dev. 1982, 21, 308-323. (9) Seider, W. D.; Seader, J. D.; Lewin, D. R. Process design principles; John Wiley & Sons: New York, 1999. (10) Morari, M. Design of resilient processing plantssIII. Chem. Eng. Sci. 1983, 38, 1881-1891. (11) Garcia, C. E.; Morari, M. Internal model controlsII. Ind. Eng. Chem. Process Des. Dev. 1985, 24, 472-484.
Received for review October 10, 2000 Revised manuscript received December 14, 2001 Accepted June 12, 2002 IE000886F
