1530
Znd. Eng. Chem. Res. 1991,30, 1530-1541
Park,S.; Ramirez, W. F. Optimal Production of Secreted Protein in Fed-Batch Reactors. AZChE J. 1988,34, 1550-1558. Luus, R. Sensitivity of Optimal Control to Final State Rosen, 0.; Specification by a Combined Continuation and Nonlinear Programming Approach. Chem. Eng. Sci. 1989,44, 2527-2534. Rosenbrock, H. H.; Storey, C. Computational Techniques for Chemical Engineers; Pergamon Press: New York, 1966; pp 241-277.
Sage, A. P.; White, C. C. Optimum S y s t e m Control; Prentice-Hall: Englewood Cliffs, NJ, 1977. Van Dooren, R. A Chebyshev Technique Applied to a Controlled Nuclear Reactor System. Optim. Control Appl. Methods 1989, 10, 285-291.
Received for review September 24, 1990 Accepted March 19,1991
An Improved Autotune Identification Method Wei Li, Esref Eskinat, and William L. Luyben* Department of Chemical Engineering, Lehigh University, 1 1 1 Research Drive, Bethlehem, Pennsylvania 18015
The method proposed by Luyben (1987) for obtaining simple transfer function models from one autotune test requires that the steady-stategain be known. This limits its usefulness for identification when only plant data are available. This paper presents a modified procedure that does not require knowledge of the steady-state gain. The method uses two autotune tests. The first is a normal one; the second is run with an additional known dead time added so that the phase angle is shifted about 45’ and a smaller ultimate frequency is obtained. Then a least-square method is used to determine the unknown parameters: two time constants and the steady-state gain. The method was demonstrated on several simulated processes and successfully tested on an experimental heat-exchanger process. Autotune testing of both continuous and sampled-data systems were considered.
Introduction To design a control system for a process, one usually needs a good dynamic process model. But quite often the process is poorly understood or some parameters are unknown that makes it difficult to develop a model. In some cases, even if knowledge of the system is available, a rigorous model would be too complex and nonlinear. Because almost all control analysis and design techniques require linear process models, simplified linear models are usually required. One solution to these problems is the use of identification techniques to obtain a linear transfer function model suitable for controller design purposes. Several identification techniques exist, including direct sine-wave testing, step testing, pulse testing, and pseudorandom binary sequence (PRBS)testing. All of these have their limitations when applied to nonlinear processes. Another alternative identification method is autotune testing (Astrom and Hagglund, 1984). Autotuning testing uses a relay in the feedback loop to produce a limit cycle. The period of the cycle and the amplitudes of the output and input can be used to quickly determine the ultimate gain and ultimate frequency of the process. It is a closed loop test, so the process is kept in the linear region. I t provides accurate information at a frequency that is important for feedback controller design. A method for deriving transfer functions of highly nonlinear processes using information obtained from the autotune method was developed by Luyben (1987a). The proposed procedure is as follows: (1)Obtain the ultimate gain and ultimate frequency by using Astrom’s autotune method. (2) Obtain the dead time from the initial response of the system to the autotune test. (3) Obtain the steady-state gain from a steady-state model of the process. (4) Calculate the time constants (one or two) of simple transfer functions that fit the data at the zero frequency and ultimate frequency points. The principal reason for wanting to develop a parametric (transfer function) model is that models of this type are used in most controller design methodologies. This is
* To whom all correspondence should be addressed. 0888-5885/91/2630-1530$02.60/0
particularly true for model-based controllers where the process model is used explicitly in the controller structure. Transfer function models (or their equivalent step response models) are vital for the analysis of multivariable processes. The method was tried on several simulated processes, and the transfer functions derived matched quite well with those obtained from analytically derived linearized models,. It has been successfully applied to many simulated and real industrial processes. However, this method requires that the steady-state gain be known. This means that the method works well in an engineering analysis environment when rigorous nonlinear steady-state and dynamic models are available, and steady-state gains are easily obtained. However, when we are trying to determine a transfer function model from only plant data, the steady-state gain is usually not known without performing some additional testa like a step test. And obtaining accurate linear steady-state gains experimentally from plant data is very difficult because of noise and nonlinearity (Luyben, 1987b). This requirement of knowing the steady-state gain limits the usefulness of Luyben’s original method in an experimental environment. The purpose of this paper is to present a modified method that does not require knowledge of the steady-state gain.
Accuracy of the Autotune Method The autotune test places a relay in the feedback loop. If the process has a phase lag of at least -T radians, the relay feedback will cause the system to oscillate. Measurement of the oscillation will give important information about the process dynamics that can be used to calculate the ultimate gain and the ultimate frequency: K, = 4h/(a7f) (1) w , = 27f/P,
where h = height of the relay, a = amplitude of principal harmonic of the output, and P, = period of limit cycle. Since the relay is a nonlinear device, the ultimate gain and ultimate frequency obtained from an autotune test will 1991 American Chemical Society
Ind. Eng. Chem. Res., Vol. 30, No. 7, 1991 1531 be somewhat different than the real ultimate gain and ultimate frequency of the linear process. To quantify this, several processes were explored and comparisons were made between the real linear ultimate gains and ultimate frequencies and those obtained from the nonlinear autotune method. The comparisons were made for the following transfer functions: (3)
taken from the magnitude of the transfer function, the other from the argument: model 1 = (I/~,)[(K/M)~ 110.5 (11) 7 = -(l/wJ tan (A + w,D) (12) model 2 T
model 3
T
(5)
=
( 1 / w u ) [ ( ~ / ~ 2 /3 110.5
(13) (14) (15) (16)
tan [(A + w,,D)/3] K model 4 M = (17) [I + ( o , T ~ ) ~ ] + ~ ~(w,T~)~]O.~ ~[~ A = -w,D + arctan (-w,TJ + arctan ( - 0 ~ 7 2 ) (18) 7
(4)
= (I/~,)[(K/M)- 110.5 = -(l/w,,) tan [(A + w$)/2] = -(l/w,,)
77
All the parameters of the transfer function of the process (time constants, steady-state gain, and dead time) were specified. We simulated each process with relay feedback using standard numerical integration methods. From the simulation we obtained the autotune predictions of the ultimate gain and the ultimate frequency. For comparison, knowing all parameters of the transfer functions, we are able to rigorously calculate the ultimate gain and ultimate frequency by determining the frequency where the phase angle equals -T and the magnitude at that frequency. For each of the transfer functions we set the steady-state gain equal to 1and changed the time constant 7 or dead time D one at a time to see the effect of each parameter on the errors in K, and w, obtained by the autotune method. The dead time was changed from 0.1 to 4, and the time constant was changed from 0.1 to 20. Results are given in Table I. It can be seen that (1)the autotune results usually give errors in the 5 % range but can give errors as large as 20%, (2) the error usually increases as the dead time to time constant ratio increases, and (3) the error becomes less as the order of the system increases. These results indicate that autotune testing should give reasonable predictions of ultimate gains and frequencies for most processes, particularly since chemical processes are usually nonlinear anyway. However, this inherent inaccuracy of the nonlinear autotune method may lead to difficulties if we try to back-calculate parameter values from an autotune test for a linear system. Review of the Original Method (ATV1) Luyben's original method involves one autotune test to determine the ultimate gain and ultimate frequency. The steady-state gain and dead time must to be known. Five possible models were assumed: Ke-D8 (6) model 1 Gl(d = 7s + 1 model 2
%a)
=
nv -e- D- s(7s + 1 ) 2
(7)
model 3 model 4
A = -w,D + arctan (-W,71) + 2 arctan (-W,72) (20) where M = magnitude of transfer function (M= l/K,), and A = argument of transfer function ( A = - T ) , and w, = ultimate frequency. These equations are solved for each model, and the simplest model that fits the data is chosen. For models 1-3 the appropriate pairs of equations are solved for two values of T , and their values are compared. If they are approximately the same for any of the models, that model should be used. If not, other models must be tried. Models 4 and 5 have two time constants, so eqs 17 and 18 and eqs 19 and 20 must be solved simultaneously. One technique for solving these equations is to guess a value for 71 and then solve the magnitude equation for 7 2 . The argument, AdC, is obtained from the second equation. If the actual argument A (-r) is equal to Adc, the correct values of T~ and 72 have been found. Sometimes there is no physical solution; sometimes there are multiple solutions. If more than one model satisfies the data, the lowest order model should be preferred. There is, of course, no guarantee that any of the models will fit the data. In this method the unknown parameters are the time constants: 7 or 71 and T* The dead time and steady-state gain are assumed to be known.
Initial Attempts To Modify the Method Muhrer (1989) did some preliminary exploration of the possibility of deriving all four parameters of the transfer function (K, D,71, and 72). The idea was to use the information from two autotune tests, at two different frequencies, to obtain four equations in four unknowns. The first test was a conventional one. In the second test, the phase angle and frequency were shifted by adding either hysteresis or additional dead time in the feedback loop with the relay. For example, consider the second-order process of model 4 (eq 9). Equations 17 and 18 can be written for the two autotune tests. Let Dz be the additional known dead time that is added to the loop in the second test: Al = wlD + arctan
M2 =
where K = known steady-state gain, D = known process , 72 = unknown time constants. dead time, and 7 , T ~ and From each of these models, two equations can be derived for use in determining the unknowns. One equation is
+ arctan (-W172)
(-w~T~)
K
[1
+ ( ~ d 2 7 1 ) ~ ] ~ " [+1 (02T2)2]o's
A2 = -w2D - w2D2 + arctan
(W271)
(22)
(23)
+ arctan ( - 0 ~ 7 ~ ) (24)
where M1= l/KUl,M2 = l/Kd, K,, and w1 are the ultimate
1532 Ind. Eng. Chem. Res., Vol. 30, No. 7, 1991 Table I. Errors in Autotune Results W,
D
r
w
wAW
0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
0.1 0.5 1.0 3.0 5.0 10.0 20.0 50.0 0.1 0.5 1.0 3.0 5.0 10.0 20.0 50.0 0.1 0.5 1.0 3.0 5.0 10.0 20.0 50.0
20.29 16.89 16.32 15.92 15.83 15.77 15.74 15.72 5.31 4.06 3.67 3.34 3.26 3.20 3.17 3.15 2.86 2.29 2.03 1.76 1.69 1.63 1.60 1.58
21.37 17.36 16.62 16.19 16.03 15.95 15.95 15.95 5.55 4.23 3.79 3.40 3.31 3.23 3.20 3.17 2.95 2.40 2.11 1.80 1.72 1.65 1.61 1.59
0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.5 0.5 0.5 0.5 0.5 0.5 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0
0.1 0.5 1.0 3.0 5.0 10.0 20.0 0.1 0.5 1.0 3.0 5.0 10.0 20.0 0.1 0.5 1.0 3.0 5.0 10.0 20.0
13.07 6.22 4.44 2.57 2.00 1.41 1.00 4.57 2.61 1.92 1.14 0.89 0.63 0.45 2.63 1.72 1.31 0.79 0.62 0.44 0.31
13.26 6.05 4.25 2.42 1.87 1.31 0.92 4.73 2.64 1.90 1.10 0.85 0.59 0.42 2.70 1.76 1.32 0.78 0.60 0.42 0.30
0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.5 0.5 0.5 0.5 0.5 0.5 0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0
0.1 0.5 1.0 3.0 5.0 10.0 20.0 0.1 0.5 1.0 3.0 5.0 10.0 20.0 0.1 0.5 1.0 3.0 5.0 10.0 20.0
9.16 2.82 1.54 0.55 0.34 0.17 0.09 4.00 1.83 1.15 0.48 0.31 0.16 0.08 2.43 1.36 0.92 0.42 0.28 0.15 0.08
9.24 2.79 1.52 0.55 0.33 0.17 0.09 4.11 1.84 1.14 0.48 0.30 0.16 0.09 2.48 1.38 0.92 0.42 0.28 0.15 0.08
% error
5.34 2.78 1.85 1.74 1.23 1.11 1.32 1.44 4.58 4.35 3.29 1.78 1.32 0.88 0.72 0.50 2.95 4.85 4.14 2.40 1.72 1.05 0.67 0.41
Ku
KU,
D A. First-Order Model
K,
% error
2.26 8.50 16.35 47.76 79.18 157.72 314.80 786.03 1.13 2.26 3.81 10.07 16.35 32.06 63.47 157.72 1.04 1.52 2.26 5.37 8.50 16.35 23.06 79.18
2.03 7.10 13.51 39.39 64.96 129.26 258.52 646.44 1.28 2.02 3.25 8.33 13.44 26.22 51.82 128.48 1.27 1.47 2.02 4.50 7.04 13.41 26.16 64.44
1.46 -2.70 -4.15 -5.86 -6.46 -7.05 -7.56 3.55 0.94 -0.85 -3.42 -4.43 -5.64 -6.49 2.62 2.40 0.82 -1.95 -3.13 -4.54 -5.69
2.71 10.68 20.67 60.67 100.67 200.67 400.67 1.21 2.71 4.69 12.67 20.67 40.67 80.67 1.07 1.74 2.71 6.68 10.68 20.67 40.67
2.56 10.10 19.29 55.20 90.63 178.59 352.91 1.30 2.56 4.45 11.90 19.19 37.10 72.53 1.27 1.67 2.56 6.32 10.04 19.15 37.05
0.84 -1.11 -1.35 -1.43 -1.43 -1.02 4.80 2.80 0.48 -0.51 -1.31 -1.41 -1.12 5.73 2.32 1.59 0.37 -0.98 -1.27 -1.35 1.36
2.50 5.15 6.22 7.28 7.55 7.77 7.88 1.25 2.50 3.54 5.46 6.22 6.98 7.45 1.09 1.76 2.50 4.25 5.15 6.22 6.98
2.42 5.01 6.05 7.09 7.36 7.59 9.29 1.32 2.42 3.44 5.30 6.04 6.80 9.04 1.27 1.72 2.42 4.13 5.00 6.06 7.43
-10.41 -16.45 -17.39 -17.54 -17.95 -18.04 -17.88 -17.76 13.23 -10.76 -14.69 -17.27 -17.79 -18.21 -18.36 -18.54 22.41 -3.08 -10.86 -16.18 -17.22 -18.00 -18.38 -18.61
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0
W,
K", KuAw 70 error
T
w
wAw
% error
K,
0.1 0.5 1.0 3.0 5.0 10.0 20.0 50.0 0.1 0.5 1.0 3.0 5.0 10.0 20.0 50.0 0.1 0.5 1.0 3.0 5.0 10.0 20.0 50.0
1.50 1.29 1.14 0.95 0.90 0.84 0.82 0.80 1.01 0.91 0.82 0.68 0.63 0.58 0.55 0.54 0.77 0.70 0.64 0.53 0.49 0.45 0.42 0.41
1.52 1.34 1.20 0.99 0.92 0.86 0.82 0.80 1.02 0.94 0.86 0.70 0.65 0.59 0.56 0.54 0.77 0.72 0.67 0.56 0.51 0.46 0.43 0.41
1.61 4.51 4.77 3.41 2.57 1.60 0.94 0.46 1.10 3.78 4.72 4.00 3.19 2.10 1.26 0.60 0.83 3.15 4.44 4.38 3.65 2.51 1.55 0.75
1.01 1.19 1.52 3.03 4.59 8.50 16.35 39.91 1.01 1.10 1.29 2.26 3.29 5.89 11.12 26.82 1.00 1.06 1.19 1.88 2.64 4.59 8.50 20.28
1.27 1.30 1.47 2.62 3.87 7.03 13.39 32.51 1.27 1.28 1.34 2.01 2.82 4.92 9.15 21.88 1.27 1.27 1.30 1.73 2.31 3.86 7.03 16.57
25.92 9.12 -3.09 -13.59 -15.72 -17.31 -18.09 -18.55 26.67 16.28 3.69 -10.92 -14.15 -16.55 -17.72 -18.42 26.95 20.20 9.11 -8.24 -12.54 -15.77 -17.35 -18.28
1.43 1.08 0.86 0.55 0.43 0.31 0.22 0.98 0.79 0.66 0.44 0.35 0.25 0.18 0.75 0.63 0.54 0.37 0.30 0.22 0.16
1.45 1.11 0.88 0.55 0.43 0.30 0.21 0.99 0.82 0.68 0.44 0.35 0.25 0.17 0.75 0.65 0.56 0.37 0.30 0.21 0.15
1.55 3.16 2.34 -0.27 -1.56 -3.19 -4.58 1.09 3.05 2.92 0.75 -0.56 -2.28 -3.79 0.83 2.74 3.12 1.44 0.17 -1.59 -3.09
1.02 1.29 1.74 3.69 5.68 10.68 20.67 1.01 1.16 1.43 2.71 4.03 7.35 14.01 1.01 1.10 1.29 2.22 3.20 5.68 10.68
1.27 1.33 1.67 3.49 5.38 10.02 19.14 1.27 1.28 1.42 2.56 3.81 6.94 13.09 1.27 1.27 1.33 2.10 3.02 5.38 10.02
24.78 3.43 -4.23 -5.42 -5.36 -6.10 -7.42 26.11 10.84 -0.74 -5.61 -5.40 -5.59 -6.57 26.62 15.91 3.43 -5.49 -5.54 -5.41 -6.11
1.37 0.92 0.68 0.35 0.24 0.14 0.08 0.95 0.71 0.55 0.31 0.22 0.13 0.07 0.73 0.58 0.46 0.27 0.20 0.12 0.07
1.39 0.95 0.69 0.35 0.24 0.14 0.08 0.96 0.73 0.56 0.31 0.22 0.13 0.08 0.74 0.59 0.47 0.27 0.20 0.12 0.07
1.46 2.45 1.53 -0.28 -0.84 -1.08 7.81 1.04 2.53 2.12 0.32 -0.42 -1.01 3.45 0.81 2.37 2.39 0.79 -0.04 -0.83 1.58
1.03 1.34 1.76 3.08 3.93 5.15 6.22 1.01 1.19 1.48 2.50 3.24 4.44 5.63 1.01 1.13 1.34 2.15 2.80 3.93 5.15
1.27 1.37 1.72 2.98 3.81 5.01 8.10 1.27 1.29 1.47 2.42 3.15 4.32 6.43 1.27 1.28 1.37 2.09 2.72 3.82 5.53
23.84 2.26 -2.46 -2.95 -2.91 -2.85 30.34 25.61 8.35 -0.53 -3.02 -2.93 -2.67 14.24 26.31 13.36 2.26 -3.00 -2.99 -2.82 7.41
B. Second-Order Model -5.34 -5.37 -6.70 -9.01 -9.97 -11.00 -11.92 7.40 -5.47 -5.16 -6.12 -7.16 -8.78 -10.09 19.13 -4.23 -5.56 -5.39 -5.98 -7.35 -8.90
2.0 2.0 2.0 2.0 2.0 2.0 2.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0
0.1 0.5 1.0 3.0 5.0 10.0 20.0 0.1 0.5 1.0 3.0 5.0 10.0 20.0 0.1 0.5 1.0 3.0 5.0 10.0 20.0
C. Third-Order Model -2.97 -2.74 -2.73 -2.59 -2.51 -2.29 17.89 5.32 -2.90 -2.78 -2.84 -2.78 -2.55 21.41 16.91 -2.45 -3.00 -2.89 -2.88 -2.52 6.53
gain and frequency from first ATV test, Ku and o2are the ultimate gain and frequency from second A T V test, D = process dead time, and D2= additional known dead time. Muhrer (1989) attempted to solve these four equations D, T ~ and , TJ. However, he was for the four unknowns (K,
2.0 2.0 2.0 2.0 2.0 2.0 2.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 4.0 4.0 4.0 4.0 4.0 4.0 4.0
0.1 0.5 1.0 3.0 5.0 10.0 20.0 0.1 0.5 1.0 3.0 5.0 10.0 20.0 0.1 0.5 1.0 3.0 5.0 10.0 20.0
able to find solutions in only a few cases. Most of the time, the iterative solution methods that he used did not converge. It is believed that these problems were due primarily to the inherent inaccuracy in the K, and w, obtained from the nonlinear autotune procedure. This sen-
Ind. Eng. Chem. Res., Vol. 30, No. 7,1991 1533 sitivity means that applying this approach to real plant data, with ita inherent noise and inaccuracy, would be difficult. The extreme nonlinearity of the equations with respect to the parameters also makes them difficult to solve. A Modified Autotune Method (ATVB) In this paper we present a method that does not require knowledge of the steady-state gain K . It is assumed that we do know the dead time D. Dead times are usually easy to determine by simply looking at the initial response to a step test. A. Models. Most chemical processes can be adequately modeled using simple first-, second-, and third-order transfer functions with dead time and steady-state gain. Therefore we assume that the transfer function should be one of the forms given below. The method can be easily extended to other three- and four-parameter models. model 1 (first-order model)
model 2 (second-order models)
(27)
model 3 (third-order models)
GS(d =
Ke-D8 (71s + l)(r22s2 + 2l7@
+ 1)
(29)
The models can also be written in polynomial form: model 1 model 2
G(d =
boe-D8 als2 + a$ + 1
where the parameters bo,ao,a1 and a2 are functions of the unknowns: model 1: K and T model 2: K , q , and r2 or K , 7 , and l model 3: K , T ~ 72, , and 7 3 or K , T ~ T, ~ and , l B. Outline of Procedure. For each model, one autotune test yields two equations from the magnitude and the argument of the transfer function. If we do another autotune test, while adding additional known dead time D2,we can obtain another two equations for each model. We assume that the process dead time D is known. The procedure is divided into three steps: (1) Get ultimate gains and ultimate frequencies from two autotune tests. The second test includes an additional dead time. The additional dead time added is determined by the method described below. (2) Use ultimate gains and ultimate frequencies together with the dead times to calculate the unknown parameters
Figure 1. Nyquist plot of process showing two frequencies of autotune tests.
for all models: ao, bo for model 1; al, ao,and bo for model 2; a2,al, ao, and bo for model 3. As we will show below, the problem can be formulated as a linear least-squares problem. For example, model 1 has two parameters, but we have four equations. Therefore we find the values of the two parameters that do the best job in satisfying the four equations in the least-squares sense. Note that we could perform three or four autotune tests if we wished, giving us six or eight equations, and still use the same least-squares procedure to determine the two "best" parameters. ( 3 ) Select the best model by examining both the calculated parameters for each model and the residuals of the least-squares solution. For an open-loop stable process, the correct model should parameters that are all positive. In some cases, more than one model satisfies this condition, and we need to look at the residuals of the least-squares solution. After the correct model is chosen, the time constants, damping coefficient and steady-state gain can be calculated from the relationships between the parameters ao,al, a2,bo and time constants 7 , T ~q, , 73or damping coefficient l and steady-state gain K . C. Discussion of Step 1. To explore several cases, we choose different known transfer functions for simulation studies. The transfer functions we used for simulation ranged from first order to fifth order. The autotune testa were done by simulating the known transfer functions with relay feedback. The steady-state gains used for all transfer functions were equal to 1. The relay heights were 1. The additional dead time D2for the second test was chosen to make the additional phase angle roughly r / 4 . This is because we are interested in the range of the Nyquist plot from phase angle - 3 r / 4 to -r. A Nyquist plot of a typical process is shown in Figure 1. The frequency points w1 and w2 correspond to the frequencies obtained from the two autotune tests:
ei = w i ~ i
(33)
where Di = additional dead time (for the first autotune test
D 1= 0), Bi = additional phase angle, and i = the number
of autotune test (usually i = 1, 2). The additional dead time is determined in the following manner. Since we know that we want to make the additional phase angle roughly equal to r / 4 , we let B2 = 1r/4. But in eq 33 both w2 and D 2are unknown, so we must specify one and calculate the other. From the first autotune test, we known w1 and also we known w2 is going to be smaller than wl. So we choose w2 equal to a fraction of w1 and then calculate D2.For most cases, when you
1534 Ind. Eng. Chem. Res., Vol. 30, No. 7, 1991
choose w2 = 3u1/5 and calculate Dz using following equation, it will give an actual phase angle of 40-50’: D2 = 5 ~ / 1 2 w l (34)
If a third autotune test is used, the additional dead time D3 is determined in a similar way. We want an additional phase angle O3 = w3D3roughly equal to ~ / 2 . 4(75O), and we know o3is going to be smaller than 02. Therefore we choose w3 to be 2w1/5: 0 3 = 5~/4.8~1 (35) D. Discussion of Step 2. From step 1 we obtain two or more sets of ultimate gains and frequencies and dead times. To calculate the parameters ao, a,, 02, and bo by least-squares analysis, we must use an equation that is linear in the parameters. The derivation of the equations is given below. We can express the transfer function of any model in the following form: Gomi,= X i + jYi = MieiAi
(36)
where X i = real part of the transfer function at wi, Yi = imaginary part of the transfer function at ai,Mi = magnitude of the transfer function of oi,and Ai = argument of the transfer function at w.i From the ith autotune test, the frequency w i and ultimate gain Kui are calculated. Then Mi = l/& (37) Ai = -A
model 3
Gumi)=
bo[cos (wiD) - j sin (oiD)]
+ jaoq + 1
bo[cos (wiD)- j sin (wiD)] -ja2w? - ala? + jaowi + 1
Equating eq 36 to eqs 41-43 yields bO[cosCwiD) - j sin CwiD)] model 1 = xi model 2 model 3
where u = column vector with the elements x, and a = column vector with the elements a,. If we have r sets of measurements of the z’s and xm’s, then
z1 = ulTa z2 = uZTa
z3 = uBTa
(46)
Cross multiplying and separating into real and imaginary parts give the following equations relating the parameters of the transfer functions: model 1 Xi = aowiYi + bo COS ( w ~ D ) (47) Yi 3 -a@iXi - bo sin (wjD) (48) model 2 X i = a l o t X i + aooiYi+ bo cos ( q D ) Yi = alwtYi - aowiXi - bo sin ( q D )
2,
= u,Ta
The matrix U is defined
. . * .
. . . . . .
Equation 55 may be written in matrix form:
bo[cos (wiD) - j sin (wiD)] = Xi +jYi -ja2w? - alw? + jaowi + 1
+ jaoq + 1
(54)
. .
(45)
-a1w?
z = uTa
(43)
= xi + j Y ,
bO[cos( q D ) - j sin (wiD)]
(53)
or in vector form
u =
(44)
+1
+ ... + a,x,
(42)
+ jYi
jaowj
z = a1x1 + 02x2
(40)
Substituting s = Jwi and e-jDmi= cos (wiD)- j sin (wiD) into eqs 30-32 gives bO[cosCqD) - j sin (wiD)] (41) model 1 Gciw,,= 1 + jaowi -alw?
Notice that the unknown parameters (bo,ao,al, and a2) appear linearly in these equations and that all of the other quantities in these equations are known (Xi, yi, mi, and D). This means that a linear least-squares method can be used to estimate the unknown parameters in each model. The dead time D used in the equations above is the total dead time in the loop. For the first autotune test, it is the dead time of the process. For the second autotune test, it is the process dead time plus the additional dead time that has been added. E. Review of the Least-Squares Method. Suppose we have a linear equation such as that below and we want to find the values of parameters al, a2, ..., a,:
(39)
Yi = Mi sin Ai
Gciw,,=
Yi = a2wi3Xi + alw?Yi - aowiXi - bo sin (oiD) (52)
(38)
So X i and Yi can be determined from Xi = Mi COS Ai
model 2
+
model 3 Xi= -a20i3Yi + alwfXi ao~iYi+ bo COS (wiD) (51)
z = Ua
Solving the least-squares problem gives the unknown parameters (al, a2, ..., a,):
a = [UW]WTz
(58)
Now let us consider our problem. If two autotune testa are made, we obtain four measurements. The z vector for all models is
z =
(59)
(49) (50)
The U matrices and the a vectors for each model are
Ind. Eng. Chem. Res., Vol. 30, No. 7, 1991 1535 Expressing T~ in terms of a. and al gives
model 1
- a071 + a1 = 0
71'
a =
(60)
Solving for
T~
(67)
yields a0 f [(-ao)2 - 4a1]1/2
=
(68) 2 If the term inside the square root is negative, the transfer function is not G2 and G3 should be tried. (b) transfer function G3: 71
model 2
Using eqs 31 and 69, we can obtain two equations relating the time constant 7 and damping coefficient (with a. and model 3
a
=[i][ u
=
a1:
-qsy,
0 1 2 x 1
-02%
wz"2
q3x1 o:yl
dx2
0 2 ~ ~
o,Y,
cos(o1D) -olx1 -sin(olD) yY2 cos(@) 42 2 -sin(@)
1
(62) The parameters can then be determined by using eq 58. The residual for each model can also be calculated. For example, the residual for model 3 is R = lRil + + + lR4l (63) where R1 = X i + a2wl3Y1 - alw12X1 - ao.1Y1- bo COS (WID)
Rz = Yl - a2w13X1- alwl2Yl+ aowlX1+ bo sin (wlD) R3 = X 2 + a ~ 2 3 Y- 2alw22X2- aow2Y2- bo cos (w&) R4 = Y2- a2w23X2- alu22Y2+ aow2X2+ bo sin (w2D) (64) F. Discussion of Step 3. The last step is to determine the model from the relationships between the time constants or damping coefficient and parameters ao, al, and a?. First, we should decide which models fit the data best by looking at the parameters of each model. For an open-loop stable process, the parameters for the process should be greater than zero. If any of the parameters in the model is negative, the model should not be used. For all the remaining models, we can select the best by comparing residuals. If the residuals of two models are about the same, the lower order model should be selected. Finally, the time constants of the transfer function can be calculated from a0,al, and a2. The steady-state gain is always equal to bo: 1. model 1 (eqs 25 and 30)
K=bo 7=00 2. model 2 (eqs 31, 26, and 27) (a) transfer function G2: Ke-D8 -G2(s) = (71s + 1)(Tg + 1) rlrgz
(65)
Ke-D8
+ + (T1
T2)S
+1
(66) First we assume the transfer function has the form G2,so the unknown time constants are rl and T ~ From . eqs 31 and 66 we can relate time constants T~ and T~ with parameters a. and al: (Io = 71 + 7 2 a1 = 7172
a, =
72
a. = 2(7
(70)
f = Uo/27
(71)
Therefore 7
= ai1/'
3. model 3
Ke-D8 717273S3
+
(7172
+ 7173 +
72T3)S2
+ (TI+ 72 + 73)s + 1 (72)
Using eqs 32 and 72, we get
+ + T3 + 7173 +
a0 = 71 a1 = 7172
a2 =
72
7273
(73)
717273
which yield
- Q 0 7 1 +~ a171 - a2 = 0
(74) Solving for r1 should give three roots. Choose the one with a positive real part and zero imaginary part. If no root satisfies this condition, the model is not third order. After 71 has been determined, 7 2 can be found by using eq 76 A E a2/71 = 7 2 7 3 T13
- 71 = 7 2 + 73
B Solving for
72
(75)
gives
B f [(-B)2- 4A]'/' (76) 2 If the term inside the square root is positive, r2 can be found. If it is not, G6should be used instead of G,. Then 73 can be calculated from the first line in eq 73. Ke-D* (b) G6(8) = (71s + 1)(T22S2+ 2fTg + 1) Ke-Da (77) 717g3 + (27172( + S2)S2 + (71 + 2 t T z ) S + 1 If the G6transfer function is used, the relationships be, coefficient tween the time constants rl and T ~ damping ( and parameters ao, al, and a2 are 72
=
a2 =
T
~
T
+
~
= 2217zf 72' 00 = 71 + 2721: a1
~
(78) Expressing T~ in terms of a0,al, and a2 again gives eq 74,
1536 Ind. Eng. Chem. Res., Vol. 30, No. 7, 1991 Table 11. Calculated Parameters model 71 72 or 3 2 1
0.056 -0.001
8.069 8.056 8.020
or t a2 a1 A. First-Order Process (ATV2) t = 0.016 0.025 -0.011 -0.011
bo
R
8.067 8.054 8.020
0.993 0.992 0.988
O.oo00
8.423 8.566 -5.030
0.865 0.853 -0.501
O.oo00 0.0298 0.3977
25.20 14.38 -4.988
0.851 0.581 -0.202
O.oo00 0.1324 0.3133
34.13 6.336 -6.090
0.972 0.305 -0.293
O.oo00 0.2877 0.4855
a0
73
0.0017 0.0023
B. Second-Order Process (ATV2) -1.129
8.588 8.528
3 2 1
7.416 -5.030
1.150
3 2
13.314 { = 0.586
t = 0.909
1
0.990 12.282 -4.988
3 2 1
18.89 13.313 -6.090
8.359 = 0.238
{ = 0.912
3 2
58.82 23.53 -9.801
E. Fifth-Order Process, Two Autotune Tests (ATV2) 80.11 16.40 { = 0.649 15820 1521 5.972 { = 0.127 553.4 -9.801
1.223 0.199 -0.326
O.oo00 0.3585 0.5489
0.936 13.72 -6.764
F. Fifth-Order Process, Three Autotune Tests (ATV3) 25.50 13.50 { = 0.910 170.6 205.3 20.11 [ = 0.733 188.4 -6.764
0.873 0.767 -0.258
0.0181 0.1974 0.8262
C. Third-Order Process (ATV2) 175.4 201.2 150.8
D. Fourth-Order Process (ATV2)
1
3 2 1
r
1320
357.8 177.2
which is solved for T~ as before. Then T~ and ( are given by 7 2 = (a2/T1)1/2 ( = (ao- T 1 ) / 2 T 2 (79) The procedure outlined above is easily solved on a computer. Ultimate gains, ultimate frequencies, and dead times are obtained from experiments, and the program calculates the parameters ao,al, a2,bo and time constants or damping coefficient as well as the residuals for each model. From the results we can determine which model best fits the data. Some simulation examples will be given in the next section, and results of experimental tests will be given in a later section.
The Nyquist plots of the real process and the model obtained from the ATV2 method are given in Figure 2. B. Example 2 (Second Order). The process transfer function used was
Simulation Examples The modified autotune method (ATV2) was tested on a wide variety of simulated processes ranging from first order to fifth order. The additional dead times D2for the second autotune test were determined by using eq 34. Only a few illustrative examples will be given below. A. Example (First Order). The transfer function simulated was
The calculated parameters of the models are shown in Table IIB. It appears that only the second-order model fits the data well since all the parameters calculated for the second-order model are positive in value: 0.853e-b (83) Gmodel = (1.15s + 1)(7.416s + 1) Nyquist plots of the process and model are given in Figure 2b. C. Examples 3 (Third Order).
Gprocm
=
(80)
The ultimate gains K,, and ultimate frequencies wi obtained from two autotune tests (using an additional dead time D2 = 1.526 in the second test) were ~1 0.858 K,, = 7.031 w2 = 0.513 K,, = 4.287 The results calculated for different models are given in Table IIa. By looking at parameters ao,al, a2,and bo for the different models, we can tell that the first-order model fits the data well because no negative values of the parameters appear. So the estimated model is first-order with time constant equal to 8.02 (instead of the actual time constant of the process equal to 10) and steady-state gain equal to 0.988 (instead of 1.0): 0.988e-b Gmodel = 8.02s + 1
,-28 C
Gprocess
= (s
+ 1)(10s + 1)
The ultimate gains and ultimate frequencies obtained from two autotune tests (using an additional dead time D2 = 2.188) were 01 = 0.5983 K,, = 6.557 ~2
Gprocem
= 0.3644
(s
K,, = 3.614
+ 1)(1Os + 1)(20s + 1)
(84)
The ultimate gains and frequencies (D2 = 6.154) were ~1 0.2127 K,, 10.45 02
= 0.1190
K,, 3.848
Table IIc shows that two models satisfy the condition that all the parameters are greater than zero. So we need to compare the residuals of the two models. The residual for the third-order model is smaller than the residual of the second-order model, indicating that the third-order model fits the data better: 0.851e-% Gmodel= (0.99s + 1)[(13.3)2s2+ 2(13.3)(0.909)s + 11 (85)
Ind. Eng. Chem. Rea., Vol. 30,No. 7, 1991 1637
a
0.1
0 -0.1
g B
w23 51300
-0.2
-0.3 -04
-05
-06
1
-0.7
-0.2
0.2
0
0.6
04
0.8
1
-0.8 -0.4
-0.2
0.2
0
0.4
0.6
0.8
Re[G(iw)l
RelG(iw)l
0 0.1 0 -0 1
wl=.o65% W23.04197
-0 2
-0 3 -0.4
-0 5
-0.6
-0.7 -0.8 -0.7
-0.2
0.2
0
0.6
0.4
0.8
1
-0.9
-0.2
-0.4
0
0.2
0.4
0.6
0.8
1
1.2
RelGOw)l
c
0.1 wl
-
ReIOtiWl
Figure 2. Nyquist plots of process and of model obtained from A m 2 (solid line process; dashed line = model): (a) fiit-order procsasa; (b) eecond-order process; (c) third-order process; (d) fourth-order process, (e) fifth-order process.
Nyquist plots of the process and the model are shown in Figure 2c. D. Example 4 (Fourth Order). CprOcm
=
e-%
(8
+ 1)(5s + 1)(1Os + 1)(20s + 1)
(86)
The ultimate gains and frequencies (D2= 9.860) were ~1 0.1328 K,, 5.660 WZ
= 0.08207
Kuo 2.580
Results are given in Table 2D and Figure 2d. The thirdorder model fits best: 0.972e-% G'model = (18.89s+ 1)[(8.36)2s2+ 2(8.36)(0.912)s + 11 (87)
E. Example 5 (Fifth Order).
1538 Ind. Eng. Chem. Res., Vol. 30, No. 7, 1991
The ultimate gains and frequencies (D2 = 19.85) were K,,= 4.595 01 = 0.06596
K,, = 2.240 (89) Results are given in Table 2E and Figure 2e. The thirdorder model fits best: 1.233e-% Gmodel = (58.8s + 1)[(16.4)2s2+ 2(16.4)(0.649)s + 11 02
= 0.04197
(90)
The results from the five examples illustrate that the frequency responses of the estimated transfer function models are quite close to real process transfer functions near the crossover frequency, especially for the higher order transfer functions. F. Example 6 (Third Order Using Three Autotune Tests). To illustrate the use of more than two autotune testa, the third-order process of example 3 was used (eq 841, and three testa were run. The ultimate gains, ultimate frequencies and additional dead time were ~1 = 0.2127 K,,, = 10.45 02
A. First-Order Sampled-Data Prooess. The transfer function for a first-order process in the Laplace domain is
= 0.1190
Ku2= 3.848
D2 = 6.154
= 0.08056 K,,, = 2.273 D, = 15.39 Results are given in Table IIF. The estimated transfer function from three tests is given in eq 91 and is essentially the same as that obtained from two tests (eq 85): 0.873e-% Gmodel = (0.936s + 1)[(13.5)2s2+ 2(13.5)(0.91)s + 11 (91) The sensitivity of the proposed method to errors in the assumed values of dead time, in measurement errors of the ultimate gains and frequencies was fairly extensively studied. Details are given by Li (1990). Results indicate that the method can tolerate errors of the order of &lo% and still give reasonable models. wg
Autotune Identification in Sampled-Data Systems The experimental work on the proposed method was performed on a heat exchanger. Details of the experimental equipment and results will be given in the next section. The data acquisition and control system used on the experiment was a digital sampled-data system. When the autotune tests were performed on the experiment, they were done with a zero-order hold in the loop and a 3-s sampling period. This sampling period was not very short compared to the dead time and the time constants of the system. Therefore, we could not take the results of the sampled-data autotune tests and use the equations that apply to continuous analogue systems. The autotune test in a sampled-data environment includes both the hold and the process. So we will obtain ultimate gains and ultimate frequencies that are for the discrete system whose transfer function is HGG). Once the discrete transfer function has been identified, we can determine the continuous transfer function of the process G(,).The transformation from HG,,) to G(8jis, of course, not unique, but we use the lowest order continuous system that corresponds to the given HG,,). The discussion below illustrates how the proposed method can be used when the autotune testa are performed on a sampled-data system. The appropriate equations for first- and second-order processes will be derived. These correspond to eqs 60 and 61 that would be used when the autotune test were run on a continuous system.
For simplicity, we assume that the dead time D is an integer multiple of the sampling period T, (D= kT,). Using a zero-order hold gives an open-loop transfer function in the z domain: K ( l - b) HG,,) = (93) zk(z - b) where and k = integer multiple of sample time. The ultimate gain and frequency of the closed-loop sampled-data system can be found by using a proportional sampled-data controller with gain K,. The closed-loop characteristic equation for the system is K a ( 1 - b) 1+ (95) zk(z - b) Defining new variables and rearranging, we obtain z ( ' + ~ )- b 1zk
+ b2K c = O
(96)
where bl = b and bz = K ( l - b) (97) A t the ultimate gain (K, = K,) the roots of closed-loop characteristic equation lie on the unit circle in the z plane ( z = de). Substituting z = de= cos 6 j sin 0 (where 8 = UT,)into eq 96 and equating the real and imaginary parts yield the following two equations: COS [(I + k)eJ - bl COS (kBJ + b&, = 0 (98)
+
sin [(l + k)&] - bl sin (he,) = 0 (99) where 8, = olT,, w1 = ultimate frequency of first autotune test of sampled-data system, and K,, = ultimate gain from first autotune test of sampled-data system. These would be the two equations that apply to the first autotune test. In the second autotune test equations would be similar except an additional dead time is added (D2 = k2T,): cos [(l+ k + k2)e2]- b, cos [(k + k2)e2]+ b&,, = 0 (100) sin [(l+ k + k2)f12]- bl sin [(k + k2)e2]= 0 (101) where O2 = 02T,,w2 = ultimate frequency of second autotune test of sampled-data system, and K = ultimate gain from second autotune test of sampled-zata system. Equations 98-101 have the same structure as eqs 47 and 48, so we can use linear least-squares to solve for the two unknown parameters bl and bl. The least-squares estimate of the parameter vector is given by eq 58 with
Once the bl and b2 parameters have been identified, the values of steady-state gain K and time constant 7 can be
Ind. Eng. Chem. Res., Vol. 30, No. 7, 1991 1639
.................................................................................................
Relay
I .............................................................. FC
7 I
Steam
Steam valve
+
7 ,!: 1 , $ I.................................
I
TT T2
.
I
Water Cold
I Exchanger
,
r3
Figure 3. Experimental heat-exchanger process.
calculated from eqs 94 and 97. B. Second-Order Sampled-Data Process. For a second-order system the transfer function in the Laplace domain is
Substituting z = de= cos 0 + j sin 8 and equating the real and imaginary parts yield COS [(2 + ,well = bl COS [(I + milbz COS (k8J - b&, COS 81 + bdK,, (111) sin [(2 + k)O,] = bl sin [(l + k)8,] - b2 sin (he,) - b&, sin 8, (112)
Using a zero-order hold gives a z-domain transfer function: (104) where the time constants r1and 72 in eq 103 are related to the poles p 1 and p 2 in eq 104
p 1 I ~ - T I / T I p2 = e-TJ? (105) The constants a. and zlcan be calculated from p l , p2, 71, and r2: ao=1+ 21
=
- Pi71 71 - 72
P272
~ 1 ~ 2 (7 271)
71
+~
- pi72
2 7 1
- 72 + P272 - Pi71
(107)
The closed-loopcharacteristic equation with a proportional sampled-data controller is
Defining new variables and rearranging give d2+')- bldl+') b& b3Kg - b4K, = 0 (109) where bl = P1+ P2 bz = PtP2 b3 = a& (110) b4 = ad(z1
+
+
where O1 = wlTs,w1 = ultimate frequency of first autotune test of sampled-data system, and K,, = ultimate gain from first autotune test of sampled-data system. These would be the two equations that apply to the first autotune test. In the second autotune test equations would be similar except an additional dead time is added (D2 = k2Ts): COS [(2 + k + k2)e21 = bl COS [(I + k + k2)e21 b2 cos [(k + k2)e21 - b&,, cos o2 + b , ~ , , (113) sin [(2 + k
+ k2)8,] = b, sin [(l+ k + kz)d2] -
b2 sin [(k + k2)8,] - b3K,, sin O2 (114)
where d2 = w2T,,w2 = ultimate frequency of second autotune test of sampled-data system, and Ku, = ultimate gain from second autotune test of sampled-data system. Equations 111-114 are linear in the parameters bl, b2, b3, and bl. Therefore, the parameters can be calculated by least-squares analysis using
1540 Ind. Eng. Chem. Res., Vol. 30, No. 7, 1991 Steam flow rate (1st ATV tMt)
1
6
2
P
5 4 31
0
20
40
60
80 100 120 sampling time
140
160
180
1
200
Water exit tempilst ATV test)
62 61
u 6 0 59
I
-10 - 4 - 2
I
u,
40
Qo
80 loo 120 sampling timc
140
160
180
Steam flow rate (2nd ATV tat)
79
0
2
zbo
6
4
8
1 0 1 2 1 4
Re[Ww)l
Figure 5. Nyquist plot comparisonof models obtained by ATV and PRBS with auxiliary information: solid line = model obtained by PRBS; dashed line = model obtained by fvat and second ATV teats; dotted line = model obtained by first and third ATV testa.
I
6
I
B
5 4
d0
220
240
c
71
A0 300 320 sampling time
340
360
380
Water exit temp.(Znd ATV test)
62 I
5$0
260
2b
240
260
300 320 sampling time
280
3bo
M ! 4
I
360
Steam flow ratd3rd AT” test)
380
4ik
I
Sampling time 62 r
Water exit temp.(3rd ATV test)
I
A schematic drawing of the experimental process is shown in Figure 3. Two heat exchangers and a surge tank are involved. The autotune tests were done on the first heat exchanger by inserting a relay in the feedback loop. The controlled variable T2was the exit temperature of the process water leaving the first heat exchanger, and the manipulated variable was the flow rate of steam. The inlet water temperature TIwas kept constant at 30 OC by means of a temperature controller that regulated the cooling water flow rate of the second heat exchanger. The exit water temperature T2of the steam exchanger was controlled by a cascade steam flow controller whose setpoint was changed by the relay signal. The steam flow rate was measured by an orifice plate and differential pressure transmitter whose output signal was proportional to the square of the flow. All measurements were sent to a PC equipped with a MICROMAC 6OOo data acquisition and control system. Control algorithms were programmed in the form of a digital PID controller. Flow rates and temperatures were measured periodically with sampling time T,= 3 8. The results of the three experiments are given in Figure 3. The first autotune test was a conventional one. The second was run with additional dead time D2= 3 s. The third was run with an additional dead time D3= 9 s. The experimental ultimate gains and frequencies were K,, = 0.382 01 0.150
K,, = 0.312 K,,= 0.250 Sampling time
Figure 4. Experimental resulta of three ATV tests.
An example of using the discrete second-order model for the heabexchanger experiment is given in the next section.
Experimental Tests The ATV2 procedure was tested experimentally on a tube-in-shell heat exchanger in our laboratory. The exchanger used low-pressure steam on the shell side to heat cold water flowing through the tube side. The transfer function identified was that relating steam flow and water exit temperature.
02
= 0.123
03
= 0.095
The dead time of the process was determined from step response. This was achieved by suddenly changing the position of the steam valve. Results are shown in Figure 4. Dead time was about 6 8. A. Models Obtained Assuming Continuous System. If we assumed that the system was a continuous analog one, which it was not, the resulting estimated models were as follows: (1) Model derived from first and second tests: 16.31e* (41.26s + 1)[(3.83)2s2+ 2(3.83)(0.639)s + 11 (117) (2) Model derived from first and third tests:
Ind. Eng. Chem. Res. 1991,30,1541-1547
10.16e* G(n)= (25.33s+ 1)[(4.26)’s2 + 2(4.26)(0.619)s + 11 (118) These results are incorrect because the system is really a sampled-data one, and the ultimate gains and frequencies depend on sampling period. B. Models Obtained by Using Correct SampledData Relationships. Using the equations derived in the last section for sampled-data autotuning, we obtained the following discrete models for the experimental heat exchanger. (1) Model derived from first and second ATV tests: Parameter values found: bl = 1.2319,b2 = 0.3002,b3 = -0.01420,b4 = -0.8504 and p1 = 0.8974 r1 = 27.72 K = 11.84 72 = 2.739 p2 = 0.3345 ~-~(-0.04199z + 0.8504) (119) HG(z)= (z - 0.8974)(2 - 0.3345) 11.84e* (120) ()’ = (27.72s+ 1)(2.739s + 1) (2)Model derived from first and third ATV tests: ~-~(-0.1052 + 0.8505) (121) HG(z) (z - 0.8536)(z - 0.4023) 8.52e*
+
+
(122)
Gb) = (18.95s 1)(3.29s 1) The discrete ATV models can be compared to another discrete model obtained by Ekkinat (1990)using PRBS test with auxiliary information. The transfer function obtained from PRBS was
~-~(-0.0126z + 0.8423) (123) HG(z)= (z - 0.9027)(z - 0.3390) Nyquist plots comparing these three models are given in Figure 5. The two models obtained from the autotune testing match well with the model obtained from PRBS in the area near the crossover point. There are some differences between the models in the lower frequency region. This is not surprising because the information from the three autotune tests is concentrated in the high-frequency region.
1541
Conclusions A new method has been developed to determine transfer function models for open-loop stable processes using information obtained from two autotune tests: ultimate gains a t two different ultimate frequencies. Then a least-squares method can be used to calculate the parameters of several simple transfer functions. The best model can be selected by examining parameter values and residuals. The method worked very well on several known simulated systems and was tested successfully on a heat-exchanger experiment. The transfer function derived from the proposed method matched well with that obtained from PRBS tests with auxiliary information. The method is suitable for use with plant data since it provides estimates of the steady-state gain in addition to the time constants. The method has the potential to be very useful for practical identification problems. Acknowledgment Some of the preliminary work on the inaccuracy of the autotune method for getting ultimate gains and frequencies was done by Jose Marti while he was a visiting engineer at the Chemical Process Modeling and Control Center at Lehigh University in 1988. Literature Cited Astrom, K. J.; Hagglund, T. Automatic Tuning of Simple Regulators with Specifications on Phase and Amplitude Margins. Automatics 1984,20, 645-651. Eskinat, E. Use of Auxiliary Information and Nonlinear Methods in Process Identification. Ph.D. Thesis, Lehigh University, Bethlehem, PA, 1990. Li, W. An Extension of the ATV Identification Technique. M.S. Thesis, Lehigh University, Bethlehem, PA, 1990. Luyben, W. L. Derivation of Transfer Functions for Highly Nonlinear Distillation Columns. Ind. Eng. Chem. Res. 1987a, 26, 24S2495. Luyben, W. L. Sensitivity of Distillation Relative Gain Arrays to Steady-State Gains. Ind. Eng. Chem. Res. 1987b,26,2076-2078. Luyben, W. L. Process Modeling, Simulation and Control for Chemical Engineers, 2nd ed.; McGraw-Hill: New York, 1990,pp 662-668. Muhrer, C. A. Private communication, 1989.
Received for review September 4,1990 Revised manuscript received December 17, 1990 Accepted February 7, 1991
Multivariable Global Input/Output Linearized Internal Model Control of a Semibatch Reactor Jaydeva Bhat,
K.P.Madhavan,* and M. Chidambaram
Department of Chemical Engineering, Indian Institute of Technology, Bombay, Powai, Bombay 400 076, India
The control of a semibatch reactor is treated as a multivariable setpoint trajectory tracking problem. Global input/output linearization is used to obtain a decoupled linear system that is placed under internal model control. A coordination strategy is proposed for the operation of the two loops. Adaptive features are incorporated in the strategy to account for uncertain kinetics. The control system has a number of tuning parameters that can be used to obtain performance tradeoff. Introduction A semibatch mode of operation is often used for strongly exothermic reactions to balance the heat generated due to reaction with the available cooling capacity of the re-
actor. There are three phases in the temperature control problem in a semibatch reactor-startup phase, maintenance phase, and termination phase. In the startup phase, the reagent addition and flow of heating/cooling agent are
0888-5885/91/2630-1541$02.50/00 1991 American Chemical Society
