Ind. Eng. Chem. Res. 1994,33, 1867-1871
1867
RESEARCH NOTES Estimating Second-Order plus Dead Time Model Parameters Gade
P.Rangaiah a n d Peruvemba R. Krishnaswamy’
Department of Chemical Engineering, National University of Singapore, Singapore 0511
A simple method for determination of second-order plus dead time model parameters from process transients is presented. The method relies on three points of the process step response curve to minimize the error integral between the process and the model responses. Illustrative examples show that the method will allow a rapid and reliable estimate of parameters. Introduction In a recent paper, Huang and Huang (1993) discussed a procedure for estimating a second-order plus dead time (SODT) model from process transients. The method was developed for a specific model transfer function of the form G(s) =
exp(-8s) 5 > 0.707 T 2s2 + 25Ts + 1’
(1)
It may be noted that in eq 1 the gain is normalized to unity. For overdamped responses,the above equation can be written equivalently as G(s) =
exp(-8s) 5>1 (T,s + 1)(T# + 1)’
A major advantage of Huang and Huang’s method is that, unlike other previous methods (e.g., Meyer et al., 1967; Sten, 1970;Sundaresan et d,1978;Huang and Clements, 1982), it does not rely on graphical procedures, nor does it involve a computer-based parameter optimization. The method is so simple that its application requires only four data points read from the process step response curve. However the method lacks a theoretical basis for parameter estimation and the four data points selected are rather arbitrary. The objectivehere is to propose an alternative technique, which, while retaining the simplicity of the earlier Huang and Huang (1993) technique, would at the same time provide a sufficient theoretical basis or criterion to yield reliable parameter estimates. Application of the method will be illustrated by systematicallyvarying [in the SODT process transfer function (with and without modeling error). In addition, a dispersion model will also be used to illustrate the method developed. Parameter Estimation Let us consider a typical higher-order overdamped process reaction curve y’(t) and a SODT model approximation y ( t ) . The two curves, as shown in Figure 1, are allowed to meet at three points, namely, ( t ~ , y (tz,y2), ~), and (t3,y3). The three intersections are sufficient to estimate the three parameters involved (i.e., T I ,T2, and 8 or alternately T , ,$,and 8). Increasing the number of
* Author to whom correspondence should be addressed.
Time
Figure 1. Process step response and SODT approximation.
intersecting points beyond three could force the y’(t)and y ( t ) curves to intertwine, with the result that a smooth curve fit may not be assured. The shaded area in Figure 1,representingintegral of absolute error (IAE), is a measure of closeness of the model fit. The IAE can be obtained by suitable addition and subtraction of areas as IAE = 8 + TI+ T2- m,’+ 2 1 y [ f ‘ ( t )- f(t)l dt +
Here ml’represents first moment of the actual response and f’(t) and f ( t ) are the complements of y’(t) and y ( t ) , respectively, with respect to unity. Thus
f’(t) = 1-y’(t)
(3b)
fW = 1- Y W
(34
Derivation of eq 3 is indicated in the Appendix. In order to minimize IAE, eq 3 may be differentiated with respect to 8, T I ,and T2 and the derivatives set equal to zero. Using the analytical time domain solution of eq 2 and the derivative of IAE with respect to 8, we get, as shown in the Appendix, Y 1 + Y3
= Y2 + 0.5
(4)
The remaining two derivatives,Le., d(IAE)/dTIand a(IAE)/ dT2, resulted in complex expressions that were not amenable to algebraic manipulation. Hence the three
Q888-5885/94/2633-1867$Q4.5Q/Q0 1994 American Chemical Society
1868 Ind. Eng. Chem. Res., Vol. 33, No. 7, 1994 lllC-,.
O
t
.
I
"
.
.
.
.
-
.
.
1.0
0.7
t
.
.
1.5
Yl .
.
.
2.0
25
2.0
2.5
10 t
0 0.7
--~
"
"
1.0
"
"
"
"
"
"
1.5
'
5 Figure 2.
versus pointa of intersection (y&),
(y&),
and (y&.
derivative expressions were solved numerically covering a set of values for T, 8, and 5. For each set of T, 8, and 5, the derivative equations were solved first for y1, y2, and y3 and then for the corresponding tl, t2, and t3. As the effect of 8 is limited to shifting a given response in the time axis, it has no influence,on y1, y2, and y3. Thus the intersecting valuesof y1, y2, and y3 are found to be functions of [ only. These results are shown as solid lines in Figure 2. The corresponding values of tl, t2, and t3, with dead time removed, will also be functions of 5 only provided they are first normalized with respect to T. -These normalized times are designated as tl = (tl - OM', t2 = (t2 - 8)/T,and t 3 = (t3 - 8)/Tand are also included in Figure 2. As 5 increases (Le., as the system tends toward a firstorder plus dead time (FODT) dynamics), it can be seen from Figure 2 that the first intersecting point (t1,yl)moves closer to the base line of the actual response. Thus in the limit, for a FODT approximation, two intersecting points would suffice to evaluate the two parameters (namely the dead time and the time constant) involved. It is interesting to note that the y2 and y3 values (for 5 > 3) approach the two values (namely 0.353 and 0.853) reported earlier (Sundaresan and Krishnaswamy, 1978) for FODT approximation. Generally, 5 is not known a priori and hence Figure 2 cannot be readily used to read off values of the intersecting points that represent the optimum for a specific system. However variations in y1, y2, and y3 for the range of 5 considered (about 1.0-3.0) are limitedand hence an average value for each of the three points can be obtained as an approximation. Use of this average (estimated as y1 = 0.14,y2 = 0.55, andy3 = 0.911,although it may not produce a minimum M E for all systems, may yield integral error values closer to the minimum. Thus, based on the above development, it is recommended that for an overdamped SODT type of process the model curve be forced to pass through the points corresponding to 14%, 55%, and 91% of the actual response. The corresponding time values (i.e., tl, t2, and t3, respectively) could then be read off from the actual response and used to estimate 8,TI,and TZ(or equivalently 8, T,and €1. However, before developing the required correlations for parameter estimation, it may be pointed out that similar graphical relationships could also be obtained for the €
range 0.707 < 5 < 1.0. It is necessary to consider this range because it has been pointed out (Huang and Clements, 1982) that this region of 5 could also be applied to nonoscillatory processes. The resulting graphs, based on eq 3 coupled with the analytical solution to eq 1 for 5 < 1, are shown as dashed lines in Figure 2. It appears from this portion of Figure 2 that location of the (t.9~~3) point in this narrow region of t (Le., 0.707-1.0) could in practice be error-prone. This is because the (t3,y3) point is now located at or very close to the final steady state of the actual response. Hence the following discussion utilizes the three pointa already chosen for 5 > 1, namely, the values corresponding to 14%,55 % ,and 91 % of the actual response. It is presumed that the equations to be developed on this basis may also be applicable for the range 0.707 < 5 < 1. Necessary correlations for parameter estimation were developed by making use of the observation that tl, t 2 , and t3 depend only on 5. Consequently, the ratios (t3 t2)/(t2 - tl), (tz - tl)/T, and (t2 - 8)/T,upon which the correlations are based, also depend on 5 only. Data on tl, t2, and t3 were obtained by noting the time needed for a set of normalized SODT model responses (generated by assuming several discrete values for 5 in the range 0.7073.0 with T = 1.0 and 8 = 0.4) to attain 0.14,0.55, and 0.91 of the final response. The correlations were then obtained by fitting the resulting data to suitable polynomials in by the least-squares method. Let a = (t3 - t2)/(t2 - tl)
(5)
p = ln[a/(2.485 - a)]
(5a)
and
The inverse correlation for 5 in terms of a is found to be
5'
= 0.50906
+ 0.517438 - 0.076284P2+ 0.041363@ 0.0049224@
+ 0.00021234@ (6)
for 1.2323 < a < 2.4850 which corresponds to 0.707 < 5 < 3.000. The corresponding multiple correlation coefficient R is 1.oooO. The other correlations developed are (t2 - tl)/T = 0.85818 - 0.629075 + 1.289712- 0.36859e
+
0.038891e (7)
and
+
(t2- e)/T = 1.3920 - 0.52536t + i.299it2- 0.3601453 0.037605e (8)
The above two correlations are also for 0.707 < [ < 3.0 and R = 1.0000. Thus, the parameter estimation method developed involves the following steps: (1) calculate the process steady-state gain from the magnitude of the step input and that of the corresponding response; (2) find tl, t2, and t3 from the response data corresponding to 14%, 55%, and 91 % of the actual response; (3) calculate a,8, and 5 using eqs 5, 5a, and 6, respectively; (4) compute T and 8 by eqs 7 and 8, respectively. Illustrative Examples Example 1. An error-free second-order transfer function, given by eq 9, is considered to be the actual process:
G'(s) =
ea' T2s225'7"s
+
+1
In other words, in this example, the process transfer function G'(s) and model transfer function G(s) are meant to be identical. TI is fixed at a nominal value of unity and 8' is taken as 0.4. Several discrete values for r ,viz., 0.75,0.9,1.0,1.5, 2.0, 2.7, and 2.95, are considered. The purpose of this study basically is to see if the method proposed is (i) able to recover the original process parameter values over a range of process dynamics and (ii) applicable to the extended range 0.707 < 5 < 1.0. Some representative results are given in Tables 1 and 2. Table 1 summarizes the tl, t 2 , and t 3 values (as read from simulated process response curves) that are required for the purpose of parameter estimation. In Table 2, the parameter estimates resulting from the application of the present method are provided. For the purpose of comparison the results due to the method of Huang and Huang (1993) are included in Table 2 along with those due to an optimization procedure that minimizes IAE. The results clearly show the reliability of the present three-point optimization procedure and its applicability for [ in the range 0.707 < 5 < 3.0. Further simulation study based on eq 9 was carried out for additional values of 8' ranging from 0.1 to 1.0. The results, as expected, provided a confirmation that the choice of 8' value does not influence the relative trend observed. Example 2. A transfer function, given by eq 10, is assumed to represent the actual process: G'(s) =
(T2s2 + 25'T's +
ea' 1)(0.15T's
Ind. Eng. Chem. Res., Vol. 33, No. 7,1994 1869 Table 1. Time for the Response in Example 1 (with T' = 1.0 and 0' = 0.4) To Reach 14%, SS%, and 91% of Its Final Value .? tl tz ts 0.75 0.90 1.00 2.00 2.95
1.020 1.040 1.054 1.226 1.440
1.988 2.135 2.244 3.658 5.149
3.254 3.903 4.422 9.665 14.364
Table 2. Parameter Estimates for Example 1 (with T' = 1.0 and ?L = 0.4)
5' 0.75
model parameters f
T e M E X 102 0.90
f
T
e
IAEx102 1.00
6 T
e
IAE x 102
2.00
(
T B
IAE x 102 2.95
f
T
e
I A E X 102
parameter estimation method Hung and IAE present Huang (1993) mind 0.7500 0.9996 0.4003 0.0540 0.8998 1.0013 0.3992 0.1173 1.0004
1.oooO 0.3999 0.0759 2.0005 0.9998 0.3998 0.0256 2.9487 1.0004 0.3991 0.0823
0.7558 1.0002 0.3996 1.1821 0.9165 0.9579 0.4349 2.3075 0.9886 1.0049 0.3986 1.4504 1.9596 1.0214 0.3939 0.5618 2.8845 1.0242 0.3917 0.1909
0.7500 0.9996 0.4003 0.0540 0.8994 10009 0.3990 0.0732 1.0002 0.9998 0.3998 0.0345 1.9979 1.0012 0.3991 0.0183 2.9489 1.0003 0.4001 0.0191
Table 3. Parameter Estimates for Example 2 (with P = 1.0
+ l)(O.lTs+ 1)
and 8' = 0.4)
(10)
Equation 10 could be viewed as the same process transfer function (eq 9) considered earlier with additional time constants (viz.,O.l57" and0.17") to provide somemodeling error or deviation from second-order dynamics. As before 5' is varied between 0.707 and 3.0 with 7" and 8' taken as 1.0 and 0.4, respectively. For each of the resulting cases the corresponding tl, tz, and t 3 values are obtained. Substitution of these values in the proposed equations for parameter estimation yielded SODT approximations to the dynamics governed by eq 10. Illustrative results, obtained by various methods, are summarized in Table 3. Example 3. This example is the same as that considered by Huang and Huang (1993). It involves SODT approximation of the transient responses generated by Brenner's (1962) solution for the axially dispersed plug flow model coupled with Danckwert's boundary conditions. In contrast to the previous examples, this one is based on distributed parameter dynamics. Many chemical process systems (e.g., tubular flow systems such as packed, fluidized, and trickle beds, tubular heat exchangers, absorption towers, extraction columns) exhibit distributed parameter characteristics. Thus this example signifies wider applicability of the estimation method. The shape of the response curvesfor the axiallydispersed plug flow model varied with Pe (i.e., Peclet number describing axial dispersion), and this enabled a range of t to be covered. The results obtained by the methods under consideration are shown in Table 4. Variation in Pe from 0.2 to 6.0 is considered as this provided a usable range o f t (0.8 C 5 < 1.97). It is evident that the present method is applicable to this system as well.
6' 0.75
model parameters f
T
e
IAEx102 0.90
(
1.00
B IAEx102 f
T
T 0
IAE x 102 2.00
f
T
e 2.95
IAEx102 f
T
e
IAEx102
parameter estimation method Huang and IAE present Huang (1993) mind 0.7407 1.0292 0.6167 1.8252 0.8897 1.0290 0.6188 1.2041 0.9861 1.0320 0.6129 0.8541 1.9114 1.0544 0.6217 0.2887 2.7038 1.0980 0.6148 0.2849
0.7293 1.0692 0.5864 4.4945 0.9023
0.9908 0.6500 2.5603 0.9706 1.0418 0.6076 2.3133 1.8610 1.0857 0.6090 0.5677 2.6979 1.1021 0.6084 0.7452
0.7524 1.0073 0.6390 1.0546 0.8978 1.0126 0.6353 0.9255 0.9928 1.0191 0.6269 0.7366 1.9108 1.0545 0.6214 0.2602 2.7035 1.0979 0.6147 0.2247
Discussion Referring to the results in Tables 2-4 based on examples 1-3, it can be seen that IAE is always less than 0.05, which indicates a good match between the actual and the model responses. The present method generally gives an IAE which is less than half of that due to the Huang and Huang (1993) method. Further reduction in IAE achieved by the rigorous IAE minimization procedure is often marginal. Hence the present method, using only three points from step response data, provides a SODT model prediction
1870 Ind. Eng. Chem. Res., Vol. 33, No. 7,1994 Table 4. Parameter Estimates for Esamde 3 Pe T € e Present Method 0.2 0.241 1.967 0.054 0.6 0.343 1.265 0.133 1.0 0.362 1.103 0.202 2.0 0.355 0.954 0.322 6.0 0.283 0.802 0.534 Huang and Huang (1993) 0.2 0.273 1.762 0.037 0.6 0.368 1.208 0.111 1.0 0.388 1.047 0.184 2.0 0.350 0.946 0.332 6.0 0.306 0.754 0.519 ~~~~
IAE x 102 ~
0.157 0.198 0.357 1.006 2.357 0.181 0.579 0.938 1.169 3.652
IAE Minimization 0.2 0.6 1.0 2.0 6.0
0.242 0.343 0.368 0.337 0.241
1.952 1.265 1.112 0.981 0.894
0.054 0.133 0.210 0.339 0.573
0.071 0.198 0.310 0.769 1.445
that is comparable to that by IAE minimization and is more accurate than that by the Huang and Huang (1993) method. The wide majority of chemical process systems exhibit an overdamped higher-order type of response, and the illustrative examples chosen show that the method developed can be readily applied to such systems as well as to lightly underdamped systems. However a smaller number of systems may show a pure second-order (with zero dead time) or a FODT behavior. Application of the present method to the pure secondorder case could result in a small value for 8 which sometimes could be negative (due to inevitable numerical errors). It is worth pointing out that this will not have any adverse influence in the estimation of the parameters that are relevant (namely T and 5). This is because 8 here is calculated at the end (unlike in many other previous methods) and the resulting small 8 (positive or negative) could be ignored. As for the FODT process, application of the method developed here could yield an a 2 2.485 (a being an intermediate variable, the upper limit for which is set at 2.485). This implies that 5 for the process under consideration is large P3.0) and hence the process may be better represented by FODT dynamics. In fact the proposed method was applied to both these cases (pure second order and FODT) to demonstrate that, even when applied to situations outside its scope, the method could provide valuable clues for selecting a more appropriate model.
G(s), G’(s) = model and process transfer function,respectively IAE = integral of absolute error defined by eq 3 ml= first moment of y ( t ) ml’ = first moment of y’(t) Pe = Peclet number R = multiple correlation coefficient s = Laplace transform variable SODT = second-order plus dead time t = tine tl, tz, t 3 = time corresponding to the three intersections tetween y ( t ) and y’(t) (see Figure 1) ?I, t 2 , t 3 = (tl - 8)/T,( t 2 - f?)/T, and (t3 - 8)/T,respectively T, T I ,T2 = time constants of the SODT model 7“ = time constant of the process y ( t ) , y’(t) = model and process step response, respectively y1, yz, y3 = value of the response correspondingto tl, t 2 , and t 3 , respectively a = defined in eq 5 6 = defined in eq 5a 5, = damping coefficient of model and process, respectively 8, 8’ = dead time of model and process, respectively
Appendix Derivation of Eq 3. It may be recognized (Gibilaro and Waldram, 1972) that the first moment of a step response is simply the area enclosed by the response and the projection of its eventual steady-state value (see eq 3a). It may be further noted that the first moment of a model response is related to the model transfer function by (Shemilt and Mehta, 1971) (A.1) which for eq 2 is m, = 8 + T , + T2
(A.2)
Noting that IAE (the shaded area in Figure 1)is divided into four parts (A, B, C,and D) and making w e of the definitions of first moment as above, we get by suitable addition and subtraction of areas
Conclusions
A simpler technique for estimating SODT process parameters has been developed using only three points (corresponding to 14%, 55%, and 91%) of the actual response. The points chosen seek to minimize the discrepancy between the actual and the model responses in terms of IAE,thereby providing some theoretical consideration for parameter estimation. As the method involves use of a set of new correlations for parameter estimates, the need for a graphical technique or computer search is eliminated. The method is applicable for 5 in the range 0.707 < 5 < 3.0. Several illustrative examples, showing a comparison with an earlier published method (Huang and Huang, 1993), reveal that the proposed technique is capable of yielding good parameter estimates.
which is the same as eq 3. Derivation of Eq 4. This involves setting the partial derivative of IAE with respect to 0 to zero. For a given process, as ml’and f ( t ) remain invariant with respect to the model parameter 8, we get
Since the model is assumed to be an overdamped SODT model (eq 2),
Nomenclature
-
f ( t ) , f ( t ) = 1 y ( t ) and 1 - y’(t), respectively
FODT = first-order plus dead time
and hence
Ind. Eng. Chem. Res., Vol. 33, No. 7, 1994 1871 Gibilaro, L. G.; Waldram, S. P. The Evaluation of System Momenta from Step and Ramp Response Experiments. Chern.Eng. J. 1972, 4,197-198.
Huang, C. T.; Clementa, W. C., Jr. Parameter Estimation for the Second-Order-Plus-Dead-TimeModel. Ind. Eng. Chem. Process Des. Dev. 1982,21,601-603. Huang, C. T.; Huang, M. F. Estimation of Second-Order Parameters from the Process Transient by Simple Calculation. Ind. Eng. Chem. Res. 1993,32,228-230. Meyer, J. R.; Whitehouse, G. D.; Smith, C. L.; Murrill, P. W. Simplifying Process Response Approximation. Instrum. Control
= f ( t J- f ( t 2 )
Syst. 1967, 40 (12), 76-79.
Similarly it can be shown that
Shemilt, L. W.;Mehta, S. C. Evaluating Parameters in Modela for Flow Systems. Chem. Eng. Sci. 1971,26,1777-1779. Sten, J. W. Evaluating Second-OrderParameters. Instrum. Technol. 1970,17 (9), 39-41.
= f(t3) (A.7) Substituting eqs A.6 and A.7 in eq A.4 and rearranging,
f(t1)+ f(t3) = f ( t 2 )
+ 0.5
which in terms of y ( t ) is Y 1 + Y3
= Y2 + 0.5
Sundaresan, K. R.; Krishnaewamy, P. R. Estimation of Time Delay Time Constant Parameters in Time, Frequency, and Laplace Domains. Can. J. Chem. Eng. 1978,56, 257-262. Sundaresan, K. R.; Praead, C. C.; Krishnaawamy, P. R. Evaluating Parameters from Process Transients. Ind. Eng. Chem. Process Des. Dev. 1978,17,237-241.
Received for review November 16, 1993 Revised manuscript received March 21, 1994 Accepted April 1, 1994.
(A.8)
Literature Cited Brenner, H. The Diffusion Model of Longitudinal Mixing in Beds of Finite Length. Chern. Eng. Sci. 1962,17, 229-243.
0
Abstract published in Advance ACS Abstracts, May 1,1994.
