Approximate Dynamic Modeling of Large Staged Systems Changse Kim and John C. Friedly* Department of Chemical .Engineering. University of Rochester. Rochester, New York 14627
The general applicability of the commonly used approximate transfer function d(s) = exp(- TdS)/( 1 Tls)(l T2s) for large staged Systems is investigated. It is shown that idealized staged systems, no matter how large, always have properties such that the approximation is applicable. Only the spectrum of eigenvalues (poles and zeroes) of the exact system model determines whether the approximation is valid in general. Methods are suggested by which the phenomenological parameters TI, TP, and Td are directly related to the system design parameters and can be calculated a priori. These methods are based on the exact system eigenvalues if known and on readily computed moments otherwise, giving comparable results. It is shown that the effective delay time thus obtained can be optimal in a precise sense. Systems with numerator dynamics can be treated in the same way. Examples of an ideal cascade with an arbitrary number of stages andI a previously modeled experimental distillation column are considered.
+
Introduction It is well recognized that the linear dynamics of many real chemical processes can be adequately approximated by transfer functions of the form
Although approximate, such transfer functions have proved to be quite adequate in conventional control system design (e.g., Koppel, 1968). Using the graphical procedures of Oldenbourg and Sartorius (1948) or Smith (1959), transfer functions such as eq 1 can be readily obtained from experimental step response data. This is most often the origin of such approximate models. However, the simplicity of eq 1 and its apparent wide applicability have led other workers (Moczek, et al., 1965; Bhat and Williams, 1969a,b; Biery and Boylan, 1963, for example) to fit such an approximation to the computed step response from much more complicated and presumably more fundamentally correct mathematical models. Then eq 1can be used in more extensive simulation work. The same idea is inherent in the use of the method of moments as suggested by Paynter and Takahashi (1956) and Gibilaro and Lees (1969). The advantage of such an approach for large systems should be quite clear. Instead of using a minimum of 100 differential equations for a 100-stage distillation column, two differential equations plus delay logic yields eq 1. Plants containing such processes can be economically simulated using such approximations but perhaps cannot be otherwise. Although models such as eq 1 are naturally the first to be tried in approximating staged processes, there is little guidance available for when or why it works. The phenomenological parameters of eq 1 cannot be directly related a priori to all design parameters of a staged system. Work by Moczek, et al. (1965), Bhat and Williams (1969a,b), and more recent improvements by Weigand, et al. (1972), suggest that the delay time and the time constants of eq 1 are directly related to the inventory time of a distillation column. The inventory time is related to the total liquid holdup divided by the liquid flow rate from stage to stage. These papers extend the concept empirically to cover nonlinear column responses also. The simulation work of Biery and Boylan (1963) suggests that the time delay term in eq 1 is directly related to the sum of the smaller time constants of the system response, a result essentially equivalent to the inventory time concept.
In this paper we investigate when and why eq 1is a reasonable approximation to the dynamics of large staged systems. We give results indicating that eq 1 will always be adequate for large, but idealized staged systems. Explicit relations are derived between the phenomenological parameters of eq 1 and the process design parameters. Although these results cannot be guaranteed for all real staged processes, suggestions are made as to what conditions must be checked before they are likely to work.
Properties of Staged Systems To obtain general theoretical results for large staged systems an idealized model is necessary. We will use the linear (or linearized) model for a binary separation in a general countercurrent cascade of equilibrium stages (Marshall and Pigford, 1947) dx,/dt
j = 1,..., N
= X,.,-(~+CY)~,+CYX,+,
(2)
where x , is the mole fraction of one component in the j t h stage, t is the dimensionless time, and CY is the stripping factor, assumed independent of stage, Inlet compositions at either end x g or X N + ~will be considered candidates for disturbances. Equation 2 can be written in vector form dxldt = Ax+Bu
(3)
where
The system matrix A is the familiar tridiagonal form
0
,
.l
Although idealized this model is convenient and informative because it can be readily manipulated to show the system properties in the limit as the number of stages gets very large. In order to compare the results with more realistic systems the simplest model used by Levy, et al. (1969), will be considered as well. Their simulation work of the experInd. Eng. Chem.,
Process Des. Develop., Vol. 13,
No. 2, 1974 177
Table I. Elements of System Matrix A for a 10-Plate
Distillation Column a,.,--1
-ai,
ai.,+ i
8.4316 4.065 3.4758 3.1179 2.9346 3.3749 3.3008 3.2211 3.1473 3.0868 0,13716
0.52832 16.448 11.222 10.947 10.858 10.851 6,9273 6.8813 6,8546 6,8287 6.8057 0.12278
54266 10.733 7.7784 7.8831 7,9729 3.0522 3.662 3.6997 3.7425 3.7619 3.0753
j
1 2 3 4 5 6 7 8 9 10 11 12
,-. , -, The behavior of large systems can be investigated by looking a t the limit as N m . The first few eigenvalues then approach
-
'
The case of a = 1, which is known to be optimal for a staged system (Voetter, 1957), is the most important for wefi-designed large staged systems. In this case the following is true for the time constants ( - l / A j ) -11x1
imental column of Luyben, et al. (1964), showed that the concentration dynamics described by the linearized eq 3 gave the essence of the column response. The state vector x is the set of tray compositions starting from the reboiler up to the condenser. Table I gives the values of the elements of A obtained by linearizing about the nominal experimental operating conditions. It should be noted that A is again of the tridiagonal form as eq 4 with negative diagonal elements and positive off-diagonal elements. The elements are not independent of stage number however. Rosenbrock (1963) has shown that the same form of A always results even when the stages are not at equilibrium and the binary equilibrium curve is arbitrarily shaped. The properties of the tridiagonal A matrix are well known (see, for example, Marcus and Ming, 1964). Provided all diagonal elements are negative while the rest are positive, the eigenvalues of A are all real and distinct. Furthermore, Rosenbrock (1963) has shown that all the eigenvalues for binary distillation columns must be negative as well. (It is of course expected that eq 1 would be only used to approximate those systems known to be stable.) Levy, et al. (1969), found that some of the lesser important eigenvalues are complex when holdup and energy dynamics of the column were included. The spectrum of eigenvalues of the tridiagonal matrix eq 4 is interesting as well. Equation 4 is convenient because all eigenvalues can be obtained analytically as shown by Lapidus and Amundson (1950)
-
hIA,
5
A,/X,
1=1
Nfl
(7)
-
114
(7)
+
+I6
Therefore the time constants are well separated in magnitude. The first approaches four times the second, while the sum of all the rest amounts to only the same magnitude. Table I1 shows the computed time constants for the matrix eq 4 for various values of the system size, N . In these computations the parameter a was permitted to vary in such a way that the total separation (%I %N+&. the change in the steady-state value of concentration in the light phase from one end to the other of the set of stages, was constant. The time constants for the distillation column matrix in Table I are also included. Clearly there are only a couple of dominant time constants in all cases, the sums of the rest being smaller. Levy, e t al. (1969), found domination of the first eigenvalues even when a more complete model including holdup and energy dynamics was used. It can be expected then that staged systems are characterized by eigenvalues which are real, distinct, stable, and well separated. They tend to be dominated by the first couple of time constants. Furthermore, this tends to be more true the larger the system is, although it can be true for a modest number of stages as well. These properties
Table 11. Comparison of Exact and Approximate Time Constants for Idealized Stage Systems and Distillation Column Exact system N
u
- l/A,
Moment approximation
.v
c
--1/Xz
(-l/Aj)
TP
Ti
Optimal delay time (TI = -l/Ai, T2 =
-1 / A 2 )
Td
Td*
j=3
Staged System (El - E . v + ~ ) / ~=o 0 . 8 4 1,11 1.27 3.74
5 10 20 40
0.984 1.138 1.184 1,190
3.76 11.00 31.18 68.22
2.91 9.57 29.56
5 10 20 40
0.454 0,715 0.88 0.959
3.48 10.83 40.34 161.92
1.28 3.42 11.47 42.77
1.01
4.21 15.64 57.69
Staged System (2, *
1.54 5.21 18.27 68.08
10.95 30.95 65.62
3.58 11.83 40.28
0.87 3.59 13.62 49.58
1.11 4.24 15.72 58.00
1.17 4.40 15.98 60.19
1,54 5.25 18.36 68.34
f , ~ + ~ ) / Z= o 0.99
3.41 10.72 40.13 161.32
1.73 4.33 13.98 51.28
Staged System Asymptotically Large
N-..
j=3
1
(1 - 0,003)
Distillation Column 11.5 178
Ind. Eng. Chem., Process Des.
9.1
2.36
Develop., Vol. 13, No. 2, 1974
11.25
9.48
2.24
2.36
tend to be sufficient to justify the validity of eq 1 as an approximation to the system dynamics.
Development of Approximate Models .We will consider three ways of developing an approximate model from the system equations. For simplicity start with a normalized transfer function of the form
This can be obtained from eq 2 by transforming and solving for a concentration response at one end of the cascade to a disturbance a t the opposite end. More general transfer functions with numerator dynamics will be considered below. The usual techniques of the method of moments (Paynter and Takahashi, 1956; Gibilaro and Lees, 1969) can be used to relate the parameters in the approximation eq 1to those in eq 8. By matching the first, second, and third moments the three parameters T I , Tz, and Td can be obtained from the three simultaneous equations
+
N
Ti2 Tz2 =
+
T I 3 T23=
(-1/AJi2
(9)
J-1
c N
(-l/A1I3
J-1
It is clear that if the eigenvalues are quite well separated, i.e., -1/X1>> - l / h z >> - l / X 3 , . . . , the solution to eq 9 yields simply ' T
T,=
- l / A l + o[
(2)$1
Then again T I = - l / X 1 , T2 = -1/Xz, and Td = ZI=3N. (- l/X,). The error in the frequency domain becomes important only when the frequency is of order w = l / [ 2 1 = 3 N ( l / Xj)'I1/'. Using the asymptotic form of X,'s, that is -W/XI = 7.1. The error is then important only after the frequency response is damped by over a n order of magnitude by the time constants T1 = - 1 / X 1 and T Z= - 1 / X 2 . Both the moment method and the low-frequency approximation suggest that the time constants in the approximation should be equal to or closely related to the dominant time constants of the system. The delay time is equal to the sum of the remainder of the time constants. It is intuitively appealing to have the dominant time constant in the approximation agree with the true one because that governs the rate of approach to a new steady state. It is not so clear though that the delay time really should be related physically to the sum of the smaller time constants. The third method used to develop an approximate transfer function involves the determination of an optimal delay time Td*, given that TI = XI and TZ = -1/Xz. Td* was selected to minimize the square of the difference between the approximate and exact model impulse response
where
p 1 -where O[.] represents small correction terms of order [.I. However, using the asymptotic forms of the 1 / X j obtained in eq 7 permits us to obtain the exact numerical solutions of eq 9
(-1)j-l
2 sin? N+1
~
(&)
and
This single variable minimization can readily be carried out to yield the equation N "lA2 i = 1 (A1
PIA1
f A])(A, +A,)
exp(h,Td*) +[G(Td*)I2 = 0
(14) with their degree of approximation shown. Clearly the approximate time constants are closely related to the actual dominant time constants and the delay time is closely related to the sum of the smaller time constants. This latter result is that of Biery and Boylan (1963). Equation 9a, which represents the matching of the first moments, is essentially the result of Bhat and Williams (1969a,b), Weigand, et al. (1972), and Moczek, et al. (1965). The second way of developing the approximate model is to consider a low-frequency approximation in the Laplace transformed domain. Since,
The approximation can also be derived by expanding all of the smaller time constants in a similar manner
which Td* must satisfy. Again by using the asymptotic form of the X,'s the exact numerical solution to eq 14 is found to be ,1.
Td* =
(-l/Aj)(l-O.O3)
(15)
I =3
The last two sets of columns in Table I1 give the cornputed results for the parameters of the approximation eq 1 by the moment method and the optimal delay time method. The suggested low-frequency approximation uses the exact model time constants. Either method approximates the time constants closely. This is true for both the large and small idealized staged systems and for the 10-plate distillation column. Figures 1 and 2 show the exact and approximate step responses computed for a relatively small idealized staged system ( N = 5, (zl- z ~ + ~ ) =/ 0.84) z ~ and the distillation Ind. Eng. Chem., Process Des. Develop., Vol. 13, No. 2 , 1974
179
-1
I
I
I
/
“I
I
I
I
I
0
5
10
IS 20 25 DIMENSIONLESS TIME
I
I
I
I
1
08
O6
ti -J
f
-EXACT MODEL
I
x MOMENT MODEL 4
04
LOW FREOUENCY OR OPTlM4L DEL4Y TIME MOOEL
1 I
N :5
02
0
2
4 OIMENSIONLESS TIME
-
10
12
14
Figure 1. Approximate model responses for a five-stage idealized system.
I
I
I
I
I
30
35
40
-
I 45
Figure 3. Response of bottoms concentration of ten-stage distillation column to a feed concentration change.
Table 111. Zeroes of Distillation Column Models and Their Moment Approximations
/
Exact bottom product transfer function
-I
d
04
R
-EXACT MODEL MOMENT MODEL OLOW FREOUENCY OR OPTlM4L MODEL DELAY TIME MODEL
-l/vI -1/~z -~ 0.82 0.50
20 DIMENSIONLESS TIME
-
30
(-1/vy)
j=2
1.026
Exact top product transfer function
1
p/
c
Moment approximation T3 0.996
~
-T d ,
0.846
Moment approximation
5
40
Figure 2. Approximate model responses for ten-stage distillation column without numerator dynamics. column. The continuous approximation responses have been represented by discrete points on these figures to avoid the confusion of overlapping curves. The optimal delay time model response is virtually coincident with that for the low-frequency approximation. Clearly the agreement is quite good for any of the three methods used to obtain the approximation. The agreement tends to be even better when the staged system is larger and when the steady-state separation is better. Numerator Dynamics The problems treated above had no dynamics in the numerator of the transfer function. This situation results in staged systems when considering the response in a variable a t one end of the cascade to a disturbance a t the opposite end. When other transfer functions are considered numerator dynamics result. For example, for the ideal cascade of eq 2 the normalized transfer function relating the response a t one end X I to a disturbance at the same end xo can be written
Ming, 1964) guarantees that the value of each v j falls between X j and X j + l . The second example to be considered is the distillation column with system matrix given in Table I. The normalized transfer function relating the bottoms product concentration to a feed disturbance on the fifth tray is of the form
where the zeroes are given in Table III. The eigenvalues are of course the same as those given in Table 11. The response of the top product G T ( s ) is of the same form as eq 18 except that there are only five zeroes. For this column it is not possible to prove that the spectrum of zeroes is as nicely behaved as the A’s. However, the computed values indicate that they are. If then the spectrum of roots of the numerator is the same as that for the denominator, the approximation can be treated as before. The most convenient method, either the method of moments or the low-frequency expansion, can be used to approximate (1 - s / v d . . . (1 - s/v,) by (1 T3s)exp(-Tdl,S). Only a single important zero is included in the numerator approximation because physically it must be of lower order than the denominator. Table 111 includes the moment approximations to the numerator of the distillation column transfer function eq 18. Combining these approximations with that for the denominator shown in Table I1 yields the overall approximate transfer function of the form
+
where the V k are the zeroes given by the formula
cos(xk/N) (17) It is important to note here that the numerator has exactly one less zero than the denominator and that its roots are the eigenvalues of an ( N - 1)st order square matrix of exactly the same form as A in eq 4. Therefore the character of the zeroes of eq 16 is exactly the same as for the eigenvalues discussed above. They are real, distinct, stable, and well separated. Of course the Aj’s in eq 16 are exactly those obtained above for the nonnumerator dynamics problem. In addition, Sturm’s chain theorem (Marcus and Uk =
180
-(l+a)+2G
Ind. Eng. Chem., Process Des. Develop., Vol. 13, No. 2, 1974
If the resulting T3 is significantly lower than Tz indicating weak numerator dynamics, as is the case for the bottoms response in Table 111, it is more consistent to set Ts = 0 and increase - T d v accordingly. Figure 3 shows that the
step response for the bottoms response of the distillation column is quite well approximated by eq 19 (or equivalently in this case, eq 1).
Discussion and Conclusions It is concluded from the results presented here and others obtained by Kim (1973) that the well-known transfer function eq 1 has a great deal of theoretical justification for its use on large staged systems. It is not only convenient but also predictable a priori. Its success depends on the existence of a set of eigenvalues of the system matrix which are well separated and consequently dominated by a couple of real ones. Properties of the matrices involved for simple binary stages suggest that this is satisfied for general staged systems. It is shown to be increasingly true as the ideal system becomes larger. From a knowledge of the spectrum of eigenvalues a priori accuracy estimates can be readily made on the higher moments of the impulse response and the frequency response. It is suggested that the best method for predicting the parameters in eq 1 is to use eq 10, setting the time constants TIand Tz equal to the dominant time constants of the system and the delay time equal to the sum of the remainder. Then the rate of approach to a final steady state will be precisely correct, and also there is evidence that the effective delay time will be nearly optimal in a precise sense as well. This procedure, of course, suffers from the fact that all eigenvalues must be known. This is inconvenient if not impractical for large systems. If the eigenvalues are not known but are expected to be as above, then the method of moments (eq 9) should give equivalent results. As shown by Friedly (1972) the moments can be obtained from an exact model, even though large, with much greater ease than the eigenvalues. These suggested procedures are entirely consistent with previous results obtained by Biery and Boylan (1963), Moczek, et al. (1965), Bhat and Williams (1969a,b), and Weigand, et al. (1972). More complicated problems with numerator dynamics can be treated in exactly the same way. It is suggested that the numerator of the transfer function be treated separately from the denominator, combining the results into a final transfer function. If the numerator dynamics are strong, eq 19 is a better form of approximation than eq 1.
The approximation is of course good only for relatively low-frequency disturbances. The time responses shown suggest that it would be good for most practical process applications for responses to step change and pulses. It of course cannot be used to predict highly damped responses a t high frequencies.
Nomenclature A = system matrix in eq 3 a j k = elementsof A B = matrix i n e q 3 G = transfer function I = performanceindex N = order of staged system s = Laplace transform parameter Ta Tdv= time delays Tj = time constants t = dimensionless time u = disturbance vector x j = dimensionless concentration in j t h stage Greek Letters a = dimensionless separation factor in eq 2 P j = coefficients defined for eq 12 X j = eigenvalues or singularities of transfer functions v j = zeroes of transfer function w
= frequency
Superscripts (:) = approximate model ( - ) * = optimalvalue Subscripts B = bottoms product response T = top product response Literature Cited Bhat, P. V.. Williams, T. J., R o c . JACC, 8th, 7969,8, 135 (1969a). Bhat, P. V., Williams, T. J., Inst. Chem. Eng. Syrnp. Ser.. 32, 22 (1969b). Biery, J. C.. Boyian, D. R., Ind. Eng. Chem., Fundam., 2, 44 (1963). Friedly, J. C., "Dynamic Behavior of Processes," Prentice-Hail. Englewood Cliffs, N. J., 1972, Chapters 6 and 10. Gibilaro, L. G., Lees, F. P., Chern. Eng. Sci., 2 4 , 85 (1969). Kim, C., Ph.D. Thesis, University of Rochester, 1973. Koppel, L. B., "Introduction to Control Theory," Prentice-Hall, Englewwd Cliffs, N. J., 1968, Chapter 9. Lapidus, L., Amundson, N. R., Ind. Eng. Chem., 42, 1071 (1950). Levy, R. E., Foss, A. S.,Grens, E. A,, 11, Ind. Eng. Chem., Fundam., 8, 765 (1969). Luyben, W. L., Verneuil, V. S.,Jr., Gerster, J. A,, AIChE J . . 10, 357 (1964). Marcus, M., Ming, H., "A Survey of Matrix Theory and Matrix Inequaiities," Allyn and Bacon, Boston, Mass., 1964, p 166. Marshall, W. R., Pigford, R. L., "Application of Differential Equations to Chemical Engineering Problems," University of Delaware Press, Newark, Del., 1947. Moczek. J. S., Otto, R. E., Williams, T. J., Chem. Eng. Progr. Symp. Ser., 61 (55), 136 (1965). Oldenbourg, R. C., Sartorius, H., "Dynamics of Automatic Control," American Society of Mechanical Engineers, New York. N. Y . . 1948, pp 8889. Paynter, H. M., Takahashi, Y . , Trans. ASME, 78, 749 (1956). Rosenbrock, H. H., Automatica, 1, 31, 97 (1963). Smith, 0. J. M., /SA J , 6 (2), 28 (1959). Voetter, H., in "Plant and Process Dynamics Characteristics," Butterworths. London, 1957, p 73. Weigand, W. A.. Jhawar, A. K . , Williams, T. J., AIChE J.. 18, 1243 (1972).
Received for reuieu: September 10, 1973 Accepted January 9, 1974
Ind. Eng. Chem., Process Des. Develop., Vol. 13, No. 2, 1974
181
