I
R. S. SCHECHTER and E. H. WISSLER Department of Chemical Engineering, The University of Texas, Austin, Tex.
Frequency Response from Step Input Response This is a systematic and accurate method of computing the frequency response of a system when the step response is given
THE
quality of automatic control which can be obtained with a given controller is determined by the process characteristics. Techniques have been developed to enable the control engineer to compute the closed loop response to load changes if the open loop frequency response of the process system is known (2). The frequency response of the plant also yields an immediate qualitative indication of the controllability and stability of the system. Unfortunately, the frequency response of an industrial unit is difficult to obtain directly, as a sinusoidal signal must be applied to the final control element, and if the system has a large time constant, a long period is required before the transient portion of the response damps out and the phase lag and amplitude ratio can be measured. Furthermore, phase lag and amplitude ratio must be determined for a spectrum of input frequencies to define completely the phase and gain curves of the system. This may require many days of testing, during which an off-specification product is produced. A far more desirable procedure involves putting a unit step change into the final control element and measuring the output of the controlled variable if frequency response can be derived accurately and easily from the measured output. This article derives and demonstrates a technique which permits calculation of the frequency response characteristics of the system from the response to a unit step change. Although strictly valid only for linear systems, the techniques apply to many nonlinear systems, if the magnitude of the step change is small. Bilous, Block, and Piret (7) applied perturbation theory to a system which is described by nonlinear equations and analyzed the resulting linear equations. They demonstrated that, in one example, an input
signal amplitude of 10% of the steadystate input resulted in an output amplitude only 0.16% different from that computed from the linearized equations.
where A is the amplitude and 4 is the phase lag. The problem is simply to evaluate A ( w ) and +(w) for a spectrum of values of w, given U ( t ) .
Formulation of Problem
Short-Cut Techniques
T h e response, E,(t), of a linear system to an arbitrary bounded input, E,(t), with at most a finite number of discontinuities can be expressed (6) as
This problem has received attention from investigators in the field of servomechanisms. A great deal of work has been devoted to short-cut, approximate graphical procedures in which accuracy is sacrified for simplicity and speed (8, 9, 73, 74). Special computers have been constructed to take the Fourier transform of transient data ( 4 , 70).
where ~ ( 7 ) is called the weighting function. Once the weighting function has been established, the output from the particular system can be computed for any input by application of Equation 1. If E, is a unit step function, defined as E,(t) = 0 E,(t) = 1
to
(2)
Numerical Technique
T h e output response of a system to a sinusoidal input is given by Equation 5 , which can also be written as
Equation 1 becomes cos and the value of the weighting function at any time can be determined by measuring the slope of the response curve. However, this procedure is unsatisfactory because of the inherent inaccuracy of numerical differentiation. The response, &(ut),of the plant to a sinusoidal input of unit amplitude,
E,(t) = sin
(4)
ut
can be expressed as Sa(wt) =
Jtt.(T)
sin w ( t
- T ) dT
(5)
After the sinusoidal input has been applied for a long time, the output can be represented as Lim SO( t u ) = A ( @ )sin t- m
[ut
wt
sin
W T } ~ T
(7)
Integrating Equation 7 by parts yields
w
COS
wt
1
U ( T )cos W
T ~ T
(8)
T h e sine input must be applied for a long time to damp out the transients and obtain a sinusoidal output from which the phase lag and amplitude ratio can be determined. If 7 becomes very large, Equation 8 can be expressed approximately as follows : &(ut) g w sin w t ( ~ ’ U ( Tsin ) wrdr
+
- +(w)l (6) VOL. 51, NO. 8
AUGUST 1959
945
The problem is now resolved into evaluating the two integrals, (9)
Lb
U ( T )sin
t' is defined such that
.
j 1 - U ( t ' )1 p, and are tabulated by Filon for 1 ' increments in e up to 20' and for So increments between 20' and 45'. This table is also given by Sneddon (72). The values of these parameters can be computed from the defining equations; however, for values of e less than 20°, the series approximations given below are recommended because of the extremely accurate values of the sines and cosines necessary to yield accurate values of p: and 7 . For small values of 8 :
second-order
2
1
where h is the length of the subinterval. Then over the range (7, - h, 7 , -I- h) the function U ( r ) is approximated by a quadratic expression. The integral expressed by Equation 1 2 is evaluated by integrating over the limits T, - h to T , + h and summing the resulting expressions. This technique yields
where
sin W T d T
f U ( b )sin wb
His method entails dividing the interval ( a b ) into 2n subintervals, so that
S,(wt) =
= wJO":l(i)
1 6 = 2 [
wu
The authors, employing the procedure outlined by Filon, have developed the expression :
1
d m s i n [ut - arc tan-Fl/F*]
$ [U ( u ) sin
U ( T )cos W T ~ T
-
or
Fz
Lb
When w is large (but not infinite), it is exceedingly difficult to evaluate this type of integral, because of the rapid oscillation of the sine and cosine terms, and in general the ordinary quadrature formulas such as Simpson's rule, require dividing the range of integration into extremely small steps. Filon has devised a procedure ( 3 ) to evaluate integrals of the form
U ( t ' ) ( can be made arbitrarily small by choosing a large enough value of t ' . Integrating and simplifying Equation 9.
so(wt) E
and
W T ~ T
where
system
60
-~ 4
+
6
22,215
Largest deviations occur a t w = 6
-8
40
88
997,920
20
0 0
946
I
2
3 w
4
INDUSTRIAL AND ENGINEERiNG CHEMISTRY
5
6
'
'
'
The computations can be easily performed on a digital computer. Second-Order System. The accuracy of the method was tested for a first-order system and found to give excellent agreement with theorelical values. A more severe test is provided by applying
FREOUENCY RESPONSE the method to an underdamped secondorder system. T h e system chosen for analysis can be represented by the following differential equation :
where T is the time constant and 5 is the damping coefficient. T h e solution to Equation 17 when Ed = 1 at t 2 0 (a unit step change) is E0
=
1
- e--t/T [:.- sin wt - + cos G t ]
(18)
where G = wn
4 - 2
w n = 1/TE
For E, = sin ut, the phase lag and amplitude ratio for the second-order system can be calculated from the following equations:
9
=
arc tan
[GI 1 - -
wn2
Figure 2. Calculated values of the a m p l i t u d e ratio for a second-order system agree with theoretical values
0.6 POINTS REPRESENT CALCULATED VALUES-THE CURVE IS THE THEORETICAL R ES P 0 N S E 0.4 -
0.2 0
I
A =
2
3
4
5
6
w
and 1
p )+z ($7 [I
-
The values of U ( t ) were computed at intervals of 0.25 time unit for T = 1 and E = 1 / 4 1 + For the selected values of the parameters, G = T . Hence, the transient response oscillates with a period of 2 time units. At the end of 8 time units, U ( t ) = 0.9997. Because the system is oscillating, it is not obvious that the deviation will not become greater at some future time. However, at t = 8 the response function is at the minimum value for that cycle, and because the oscillations are damped, the absolute value of the difference between U(t) and 1 will be less than 0.0003 for all times greater than 8. Hence, t’ is taken as 8 time units. The phase lag and amplitude ratio were calculated at values of w = 1, 2, 4,and 6. The results are compared with the theoretical values in Figures 1 and 2. The largest deviations occur at w = 6 At this frequency, the phase lag is 1.7% in error and the amplitude ratio differed from the theoretical by 4.1yc. Application to Process Data. Woods (75) has described a situation in which the field instrument men were unable to stabilize the side temperature on a fractionation column which separates monoethylene glycol from water. A step change in the reflux ratio was imposed upon the system by Woods and the transient response curve obtained. Although the data obtained leave much
to be desired, the step response curve as shown by Woods was used to compute the frequency response of the distillation column by applying the methods described (Figures 3 and 4). These curves were generated in approximately 10 minutes on an IBM 650 computer using 60 time increments. Analyzing these results is complicated because the column temperature did not level off within the range of the temperature recorder used by Woods. However, by using the assumptions of column behavior and static gain proposed by Woods, the authors have determined the proper controller settings using the information given in Figures 3 and 4. These settings correspond closely to the values used successfully by Woods. Comparison to Other Techniques
Ascertaining the frequency response of a system from the transient response has been of some interest to engineers. Samulon ( 77) has described a technique which is precise for systems having a definite cutoff frequency. However, a definite cutoff frequency cannot be assigned to most process equipment, and consequently Samulon’s method offers little utility to process computations. Samulon’s technique was tested on a simple first-order system using a time increment comparable to similar computations based on the methods described in this article. T h e deviations
between the calculated and theoretical values were excessive. The error in the phase lag of a particular frequency was approximately 50% when compared to the actual value, and the amplitude ratio was in error by 18.5%. An excellent technique is described by Levenstein ( 7 ) ) in which the transient response of a linear system is expressed in the form of a difference equation. The roots of characteristic equations of the difference equation and the differential equation are shown to be related and the coefficients of the solution of the equations are identical. Hence, by computing the coefficient8 of the difference equation, the Laplace transform of the differential equation can be synthesized. T h e method given by Levenstein is exact, if the order of the differential equation which describes the system is known. Unfortunately, this is seldom the case iE practical control problems; hence the method is of limited usefulness to engineers. Huss and Donegan (5) have considered the .problem of obtaining Fourier transforms numerically. The method entails fitting the time function by a staircase function, using intervals of time. The amplitude of the staircase function can be ) calculated for the interval (7 - T ~ by the equation
VOL. 51, NO. 8
AUGUST 1959
947
I
I
I
I
I I l l /
I
I
I
I
I I I I I
1
300 -
Conclusion
200
a Q
J
150 v)
-8
50
0.I
0.2
0.3 0.4
0.6 0.8 I
w IN RADIAN / M I N U T E S
Figure 3. Phase lag of distillation column was calculated using step response data presented by Woods
-
0.6 -
0.4 -
-I02
0.30.2-
w
0
2 0.08& 0.06-Q
(1) Bilous, O., Block, H. D., Piret, E. L., A . I.CI1.E. Journal 3, 248 (1957). (2) Brown, G. S., Campbell, D. P.,
0.03002 0 I I 1 I lllll 0.01 0.02 0.04 0.1 ~LI
I
I
I l l l l ~ l
0.5 0.8 I O IN RADIAN / M I N U T E S 0.2 0.3
2 0 3.0 5.0 8.0 10.0
Figure 4. Amplitude ratio of distillation column was calculated using step response data presented by Woods Controller settings corkespond closely to those used b y Woods
where amplitude of the step change
T h e Fourier cosine integral can be evaluated numerically as follows : U(T ) COS W T dT 4T
2
=
r , sinr cos (2n Z
n =O
I&
(21)
where AT
x=--W
2
T h e Fourier sine integral can be evaluated by a similar series,
Som
U(7)sin cord7
948
T h e frequency response of a linear system can be computed accurately using the methods presented here. Other proposed schemes have been tested by using numerical examples. The methods of Samulon and Levenstein incorporate assumptions which drastically limit their usefulness to engineers. T h e analysis shown by Samulon requires that the system have a cutoff frequency, and to apply Levenstein’s method, one must know the order of the difl‘erential equation which describes the system. T h e method proposed in this article necessitates numerical evaluation of two integrals. Equations for performing the quadratures are based on the technique proposed by Filon. The numerical procedure of Huss and Donegan has been considered as an alternate technique and found to be faster for desk computation if their tables are available. However, their quadrature formulas arr less accurate than those of Filon, and for machine computation, where their tables would not be used, the two require nearly the same labor. literature Cited
0.04 -
rn =
digital computer, the method presented in the numerical section is recommended because it is no more laborious and is more accurate.
=
Huss and Donegan - -present tables to facilitate the numerical computation. This numerical procedure was used to evaluate F1 and FZ in the two examples previously considered. T h e results did not compare with the theoretical values as well as did the methods used by the authors. For the case of the secondorder system a t w = 6, the calculated value of the amplitude ratio is in error by 22% compared with 4.170 error using Equations 14 and 15, and the phase angle differs from the theoretical value by 3.470, which is comparable to the 1.7% error noted in the example. Although the method of Huss and Donegan is not as accurate as the method described in the numerical section, the computations are easier to perform on a desk calculator, if the two tables presented in their publication are available. If the tables are not available or if the calculations are to be performed on a
INDUSTRIAL AND ENGINEERING CHEMISTRY
“Principles of Servomechanisms,” Chap. 11, Wiley, New York, 1948. (3) Filon, L. N. G., Proc. Roy. Soc. Edznburgh 49, 38 (1928). (4) Furth, R., Pringle, R. W., Phil. Mag. 37, Series 7, 1 (1946). (5) Huss, C. R., Donegan, J. J., “Tables for the Numerical Determination of the Fourier Transform of a Function of Time and the Inverse Fourier Transform of a Function of Frequency with Some Applications to Operational Calculus Methods,” Natl. Advisory Comm. Aeronaut., TN 4073 (October 1957). (6) James, H. M., Kichols, N. B., Phillips, R. S., “Theory of Servomechanisms,” pp. 33-5, McGraw-Hill, New York, 1947. (7) Levenstein, H., Control Eng. 4, 90 f A ~ r i 1957). l (Sf Ludbrook, L. C., Electronic Eng. 26, No. 311, 27 (1954). (9) Oldenburg, R. C., Sartorious, H., “Dynamics of Automatic Controls,” ASME, New York, 1948. (10) Reynolds, J. B., Control Eng. 2 , 60 (October 1955). (11) Samulon, H. A,, Pruc. Z.R.E. 39, 175 (1951). (12) Sn:$don, I. N., “Fourier Transforms, p. 521, McGraw-Hill, New York. 1951. (13) Teasdale, A. R., Control Eng. 2 , 55 (October 1935). (14) Thal-Larsen, H., T r a m . A m . Znsl. Elec. Engrs. (Applications and Industry), 74,109 (2955) \----,. F. A., Control Eng. 5 , 91 (15) r.Woods, -. %OCQ\ (May
LJJO,J.
RECEIVEDfor review June 20, 1958 ACCEPTED March 30, 1959
