Ind. Eng. Chem. Process Des. Dev. 1985, 2 4 , 1155-1160
1155
Controller Tuning of a Third-Order Process under Proportional-Integral Control Ll-Yang Tan” and Thomas W. Weber Department of Chemical Engineering, State University of New York at Buffalo, Buffalo, New York 14260
The control of a thlrd-order process under proportional-integralfeedback control was studied. A previous paper of Weber has been extended to include processes that are inherently underdamped. Stabilii criterla are developed, and it is shown that the commonly used continuous cycling method (CCM) of Ziegler and Nichols sometimes leads to unstable behavior for proportional-Integral control. The region of instabillty is defified in terms of the reciprocal of the maximum gain and critical frequency under proportional control. The paper by Weber is further extended to explore how the unstable region varies when the “standard”CCM parameters (0.45 for gain and 1.2 for reset time) are changed. A procedure Is then suggested for determining satisfactory controller settings through the use of a gain margin. The method is illustrated for an actual process.
Introduction A number of models have been proposed to characterize the behavior of processes for the purpose of estimating suitable controller parameters. One of the simplest is a first-order element plus a pure time delay. The model has been used in a number of studies (Cohen and Coon, 1953; Murrill and Smith, 1966; Lopez et al., 1967; Miller et al., 1967; Smith and Murrill, 1966; and Yuwana and Seborg, 1982). These studies basically differed in the criterion of optimization used in determining the best controller settings. Other investigators have used a second-order process plus a pure time delay (Latour et al., 1967; McAvoy and Johnson, 1967; and Weigand and Kegerreis, 1972). This paper is basically an extension, and to some extent a generalization, of an earlier paper by Weber and concerns a process that can be described by a third-order model. The earlier paper, Weber and Bhalodia (1979), was confined to overdampled cases, whereas the present paper includes the situation in which two of the elements combine to exhibit underdamped behavior. As explained in the earlier paper, a third-order model is reasonable because some processes consist of a series of essentially first-order elements. Moreover, if the process exhibits very little pure delay, it can often be approximated by a third-order model. An example of this is provided at the end of this paper where a process having the characteristics of a pure delay plus a first-order element is very well modeled by a third-order model used in this paper. The third-order model has several computational advantages over models which include a pure delay function, most notably in the computation of the integral of the square of the error; however, the present work did not involve that computation. The focus of this study concerns the application of the continuous cycling method (CCM) for tuning controllers which was presented by Ziegler and Nichols (1942). The paper begins with the development of equations for determining the stability of the process model under proportional-integral control. As was shown in the 1979 study, some processes will be unstable in feedback control if the Ziegler-Nichols controller parameters are used. The conditions under which this occurs are explored in some detail. Some modifications of the Ziegler-Nichols param-
* On leave from the Department of Chemical Engineering, Polytechnic University of Beijing, East Suburb, Beijing, China (P.R.C.).
eters are suggested, by which a set of proper controller settings can be easily found in practical application. Finally, the modified scheme is applied to an actual process.
Development of the Model and Stability Criteria The open-looptransfer function of a third-order process under proportional-integral controller has the form
wnL
wn
(1) where Kc and TRare the proportional-integral controller settings and T3is one of the process time constants. The remaining process elements are described in terms of the customary underdampled second-order format, where w, is the undamped natural frequency and l is the damping factor. For frequency-domain analysis, s is replaced by j w
where K = KCK2K3 = KcKp
(3)
The frequency function, GGw), is then broken down into its real and imaginary parts: GGw) = Re (a)+ j Im (w) (4) It is convenient to define the following three dimensionless parameters:
The first of these, namely, E,, renders the time constant, T3, and the natural frequency, w,, as a single dimensionless parameter. Similarly, is the dimensionless reset time, 0 1985 American Chemical Society
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 4, 1985
1156
I
-10
-08
,
I
-06
,
I
E7 , : hTRw, w , = =1 i00 i
1 ,
-04
m
-02
,
0
I
,
02
I
,
I
04
06
I
,
108
,
~_Re(w) 10
, :-Asymptote
~
, ,
Referring to eq 9b for the dimensionless critical frequency, only the positive value of the square root is meaningful because the frequencies wc and w, are positive; furthermore, only one value of XpI results, confirming that there is only one intersection of the curve shown in Figure 1 with the real axis. Special Case of Proportional Control. In the practical implementation of the Ziegler-Nichols CCM method, the system is placed under proportional action only, and the maximum gain and critical frequency are determined. These two parameters are sufficient to characterize the process. The case of proportional action only is a limiting one of that developed above as the reset time, TR, and correspondingly, q, approach infinity. Equation 9 reduces to
-i4
Figure 1. Plot of the frequency function, C(ju), in the complex plane.
and X is the dimensionless frequency. Three frequency-dependent functions are now defined: u1 = t n ( X 2 - 1) - 2{ (64 62
= XZ(1
+ 2tn{) - 1
Re
(w)
(Xu,)2
=
u
+ ;)1'2
It should be noted in passing that if the system is overdamped where { 2 1, one obtains the equations found in the textbooks of Harriott (1964) and Weber (1973)
(6b)
+ a22
(64 Equation 5a-c is now substituted into eq 2; rearrangement and the substitution of eq 6 give u =
xc = (1
T 3 ~ c=
1 + Kmax
Feasible Regions and Unstable Regions of the fmax-wC Plane In the earlier paper of Weber and Bhalodia (1979), a third-order process was characterized in terms of parameters that are related to the CCM tuning method and to a Bode diagram. The maximum gain K, is the reciprocal of the normalized amplitude ratio of a Bode diagram at the frequency where the phase lag is 180'. The value of the reciprocal at K,, is denoted f,,,; thus, 1 (12) fmax = Kmax
The dimensionless critical frequency, tc,is given by E c = tnhc = T3wc (13) Both of these parameters are positive. As was pointed out in the paper cited above, only a portion of the fm,-wc plane is feasible for real systems. That paper only dealt with the case of an overdamped third-order process; however, the concern here is to extend the analysis to underdamped processes. Turning to eq 10a and b, the dimensionless critical frequency, A,, is eliminated between the two of them; furthermore, eq 13 is used to substitute tCfor En. The result is
where
x(
&
- 1)
+ i-(;-
(llb)
where R1 = Tl/T3 and Rz = T2/T3. Furthermore, for the overdamped situation,
- u2)
Figure 1 illustrates a typical plot of Im ( w ) vs. Re (0) as w increases from zero to infinity. For the purposes of this example, the process parameters chosen were { = 0.75, w, = 10, T3 = 0.1, TR = 1, and K = 1; thus [, = T3w, = 1.0 and q = TRw, = 10. For the process chosen for this study, the curve has only one intersection with the real axis. This is evident from the shape that a Bode diagram would take, and it will be proven mathematically below. The components of this system are minimum-phase elements and the Bode stability criterion can be applied, namely, at the point where Im ( w ) = 0, (84 the system will be stable only if Re (0) 2 -1 (8b) The maximum gain occurs under the equality condition. Upon substituting eq 7 into 8, one obtains the maximum gain, KpI,and the dimensionless critical frequency, XpI, of the closed-loop system
e=
)'"
R1 + R2 + RIRz
$)
(9c)
The function 8 depends only on the system parameters, whereas the three functions u, ul, and uz are evaluated at the critical frequency here.
It can be shown then that Kmax
2(1
+ EC2Y
P tc2
I-
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 4, 1985 1157 Table I. Unstable Regions in f--[,
5, fmax
E, fmax
Plane under CCM Settings High Values of & 9.00 0.02620 3.60 0.02439
10.00 0.02466 3.80 0.02618
8.00 0.02780 3.50 0.02331
7.00 0.02934 3.30 0.02071
6.00 0.03053 3.00 0.01540
5.00 0.03060 2.50 0.00095
4.00 0.02756 2.00 0
Middle Values of 6, EC
fmax
1.50 0
0.80 0
0.75 0
0.70 0
0.45 0
0.40 0.00516
0.35 0.09131
0.65 0
0.60 0
0.55 0
0.50 0
0.25 0.56836
0.20 1.65297
0.15 0
0.10 0
Low Values o f f , EC
fmax
0.30 0.24391
The equality holds if l = 0. The three cases for the damping factor, 5; must now be considered. t = 1. In this case, eq 14 is quadratic in K,,, with the result that
Feasible but unstable under CCM settings I
1
,
S : i O ! 5.09
io 8
6
f
< 1. Here, eq 14 and 15 lead to the conclusion that
4 2
0
oi
02
03
04
fmax
f
> 1. When eq 14 is used, the result is
In terms off,,,
Figure 2. Various regions of the fmax-& plane.
and t,, these relations become
Figure 3. Three-dimensional plot of Figure 2.
found by eq 13 and 10. From these, the CCM controller settings are given by the equations (19c)
KCCM= 0.45-
1
(204 fmax These constraints are illustrated in Figure 2 which presents contour plots for various values of f in the fmm-lc 2a f C C M = 1.2h, plane. It is evident from these as well as from eq 19c that the upper bound or size of the feasible region is dependent Finally, in order for the system to be stable, KCCM must upon the value of f, and the bound vanishes as f 0. It be less than KpIas calculated from eq 9a. may be noted in passing that the contour for c = 1 corRather than just randomly selecting points from the responds to the bound used in the previous paper. Figure plane, specific values of 5, were chosen. For each one, a 3 is a three-dimensional picture of Figure 2. Regions of Instability for CCM Controller Settings. computer search was carried out to find the corresponding value off,,, if any, at which instability occurred. The In the earlier paper, regions of instability were found for results are given in Table I and also are shown in Figure the overdamped case. As shown above, any point of the 2. f,-E, plane is feasible if underdamped processes are included. The procedure for determining if a given point The results substantiate the earlier findings, namely, that for overdamped processes, if the maximum gain under is stable under CCM settings is as follows. Having selected proportional control alone is large-say greater than 30a point, the corresponding values of A,; c, and tn can be
-
1158
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 4, 1985
while the dimensionless frequency parameter, F,, is greater than about 3, the system is likely to be unstable a t the CCM settings. The present study indicates that for all processes having a very low value of [, say less than 0.5, there is a possibility that instability will occur. However, for all values of [, between about 0.5 and 2.5, instability will never occur. Modification of the CCM Parameters It is evident from the preceding section that the CCM parameters, namely 0.45 for gain and 1.2 for reset time, are not appropriate for all processes; indeed, even if the fmU-tcvalues for a particular process fell within the stable region, if they were close to the unstable region, the response of the process would be very oscillatory and, needless to say, unsatisfactory. Clearly, the CCM parameters should be a function of the process parameters, 5, and fi moreover, the gain and reset parameters, should, in principle, depend upon each other-in other words, there should be a trade-off between the controller gain and the reset time. In view of these points, a study was carried out to find some relationship between the process parameters, &, and {, on the one hand and the two controller tuning parameters on the other. Call the gain parameter a and reset time parameter p. Thus, the generalized CCM equations take the form KCCM = aKmax (214 P’WC
or
v = - 2ir P A C
K,,
and A, can be evaluated by eq 10; thus, substituting eq 10b into 10a and into 21c, one obtains
The process parameters, f , and {, could, in principle, be estimated for a process by some identification method and then be used to calculate controller settings by eq 22,21a, and 21d. It follows from eq 9a that the calculated settings will be stable if the following inequality holds:
4 Stable line at pmin= 1 2 under CCM settings 2 Unstable line a1 Bmln=4 2 under CCM settings 3 Stable line at amax: 0 45 under CCM settings 4 Unstable line at amax= 0 45 under iCCM settings
3.06 46ci\
2 0
I
1
I
I
I
02
04
06
08
10
I 12
1 44
1 46
&I,
Figure 4. amsrvs. &,, for various sets of process parameters.
expect to find a, to be less than 0.45. In particular, the point where f, = 0.0156 and 4, = 0.3775 was chosen. cy,, was found to be 0.2132. Following along these lines, ten hypothetical processes were assumed with damping factors of 0.6,0.8, 1.0, 1.5, and 2.0; furthermore, for each of these, 5, was given a value of 0.1 or 0.5. Finally, p was assigned various values between 0.2 and 1.4 in increments of 0.2. For example, one case would be a process with { = 0.6,tn = 0.1, and p = 0.2. The corresponding value of amax was calculated. The results of all these computations are presented in Figure 4 where a,, is plotted against Omin for each of the ten processes. It should be pointed out here that these are minimum values of the reset parameter, p, because each Pmin-cy,,, pair represents a system on the verge of instability; the system could be stabilized by either decreasing the controller gain or increasing the reset time. In practice, if a,, was known for a process, then some gain margin, g,, could be used to calculate a reasonable or conservative value of a; thus, (y=-
ffmax gm
A gain margin of 1.7, as mentioned in the texts by Now we shall turn the problem around. Assume that the reset parameter, p, has been set. We seek to find the maximum value of a,namely a., This maximum can be calculated from the equation (24) The value off,,, which is the reciprocal of KmaX, can be obtained from eq 22; XpI can be calculated from eq 9b, with 0 coming from eq 9c. The values of u1 and u2 can be calculated from eq 6a and b. Equation 21d is used to obtain ‘ICCM. Thus, am= is only a function of the process parameters tnand f and of the reset tuning parameter p: a m a x = J’(En,S;P) (25) For example, if p = 1.2 and we select a feasible fmax-& pair in the plane, but within the unstable region, we should
Coughanowr and Koppel (1965) and Weber (1973), is a fairly common specification. Before going on with an example, the relationship between Figures 2 and 4 should be pointed out. The unstable line in Figure 2 corresponds to the point am, = 0.45 and @- = 1.2 in Figure 4. Note the two dashed lines in Figure 4 which cross throught that point. If &in = 1.2 has been selected, then any point in the plane of Figure 2 that lies within the stable region will fall on the portion of the vertial line in Figure 4 above the value of 0.45, while all points in the unstable region will fall below that value. Conversely, if ,a has been selected as 0.45, then all points in Figure 2 within the stable region will lie to the left of 1.2 on the horizontal line, while unstable points will lie to the right of it. It is apparent then that for each point on Figure 4,there is a corresponding unstable region boundary on Figure 2.
Ind. Eng. Chem. Process Des. Dev., Vol. 24, No. 4, 1985 A The r e d process response 0 The third order mods1 response
io
Table 11. Modified Settings and CCM Settings
I case
08
c
Pmin
06
amax
04
a b m = 1.7)
02
KC TR
0
i0
20
30
40
1159
50
60
Time (sec)
ln = 0.2440
{ = 1.5633
1 0.8 0.4989 or 0.5 0.2911 7.3650 4.3942
2 1.2 0.4096 or 0.4 0.2409 6.0949 2.9294
3 1.2* 0.45" 11.3852 2.9294
Ziegler-Nichols settings.
Figure 5. Step-input response of a flow process.
Response I Use settings in Run 4 2 Use settings in Run 2 3 Use settinqs in Run 3
Application of the Controller Tuning Method The process is assumed to be described by the model used in this study, namely,
The values of the process parameters T3,a, and b would have to be identified by some method. The parameters a and b would then be recast into En and {by the equations
w, = 0.4809
llil
I
0
4
I 2
I
I
I
3
4
5
Time (sec)
Figure 6. Closed-loop response of the flow process under proportional-integral control for several controller settings.
and En
= T3wn
The process parameters a and b are both positive because the process is assumed to be inherently stable. The next step is to assume a value of the reset time parameter, P. The commonly used value of 1.2 would be a good starting value. From this, sCCM can be calculated from eq 21d, and the reset time, TR, follows directly from eq 5b. Knowing the values oft,, {, and p, one can evaluate a, from Figure 4. When a gain margin of say 1.7 is used, a can be determined from eq 26. The maximum gain, K-, can then be determined from eq 22, and the recommended controller gain, KCCM, follows from eq 21a. Finally, since the process gain, K p ,is known from the process identification, the controller gain to be used can be calculated from KC
= KCCM/KP
Example. Some tests were carried out on a liquid-flow process and the response was found to be modeled quite closely by a first-order model plus a pure delay: 0.923e-$ Gp = (30) 6.19s + 1 A typical response curve is shown in Figure 5. However, it was found that the response could be fit about equally well by the third-order model of this study; in particular, 0.923 Gp = (0.5074s + 1)(0.7520s + 1)(5.7494s + 1) The closeness of the fit is indicated in Figure 5. Anyone of the three time constants may be designated as T3. The first was so chosen, and the remaining two were combined to form the second-order form 0.923 Gp = (0.5074s + 1)(4.3235s2+ 6.5014s+ 1) The resulting values of the model parameters were calculated as outlined above, with the results that w , = 0.4809, C;, = 0.2440, and { = 1.5633. When these values and the
graphs in Figure 4 are used, the results in Table I1 were obtained. Closed-loop response tests were then carried out for the controller settings in each of the three cases. These are shown in Figure 6. In case 3 where the Ziegler-Nichols settings were used, the theory would predict that the system would be unstable. The response in Figure 6 for this case basically confirms this prediction. Comparing the responses for cases 1 and 2, about all one can say is that case 2 results in a better response than case 1. Comparing the values of Kc and TR for the two cases, it is simply a trade-off between controller gain and reset time, and apparently in this situation, the shorter reset time in case 2 is advantageous. Summary and Conclusions The previous paper of Weber and Bhalodia (1979) was extended to include third-order processes that are inherently underdamped. That paper assumed the standard CCM parameters of 0.45 for the controller gain and 1.2 for the reset time. Stability criteria were developed, and the unstable region under these CCM settings was defined in terms of the reciprocal of the maximum gain and the critical frequency under proportional control. This region is only a function of these tuning parameters. The unstable region falls within a rather irregular and not easily predictable boundary in Figure 2. However, it was found that by, in a sense, recasting Figure 2 into the a,,-bmin plane of Figure 4, the prediction of the location of the stable and unstable regions is considerably simplified. Furthermore, in the format of Figure 4, one is not constrained to the standard CCM parameters. As pointed out in the previous paper of Weber and Bhalodia (1979), p values between 0.1 and 0.6 are often better than the value of 1.2. Having chosen a value for P, one can then use Figure 4 to determine the stable region boundary and, hence, an appropriate value of a through the use of a gain margin such as 1.7. Although Figure 4 was only developed for selected values of { and tn,the calculation for other values would be very easy. Thus, the method of selecting suitable controller parameters suggested here could be easily implemented in practice. One apparent limitation of the method outlined here is the use of a third-order model. As mentioned in the Introduction section and illustrated in the example, the
1160
Ind. Eng. Chem. Process Des. Dev. 1985, 2 4 , 1160- 1 165
third-order model can be used to approximate other models in some cases. Of course, in using the ZieglerNichols method, one does need to know the process model, but as noted in the paper, the Ziegler-Nichols may fail to give satisfactory performance, whereas the method proposed here will not fail. One might well question what the advantages are of the method in comparison with such other classical methods as root-locus and frequency response (Le., Bode plots) since these latter methods are not limited to third-order models. The root-locus method is a rather complicated method and becomes especially so if a pure delay is involved. In using the Bode method, one would need to construct several frequency response plots to determine reasonable controller settings in terms of gain and phase margin. Some trial and error would be involved. On the other hand, the method proposed here is computationally simple and straightforward. The initial controller settings provided by the method will be better than those given by the method of Ziegler and Nichols. If desired, some trial and error could be used to find the "best" value of p, but this additional step may not always be necessary. Nomenclature
g, = gain margin G = transfer function K = product of controller and process gains Kc = controller proportional gain
K,, K , = process gains Kp = product of process gains KpI= maximum loop gain under proportional-integralcontrol K,,, = maximum loop gain s = Laplace operator R1 = TIIT3 R2 = T2IT3 T I ,T,, T3 = process time constants
TR = reset time Greek Symbols a = tuning parameter for gain in eq 21a
P = tuning parameter for reset time in eq 21b { = damping factor 7 = dimensionless reset time defined by eq 5b 0 = system parameter defined by eq 9c X = dimensionless frequency defined by eq 5c A, = dimensionless critical frequency E, = dimensionless process parameter defined by eq 5a Xp, = dimensionless critical frequency defined by eq 9b u, (TI, u2 = frequency-dependentparameters defined in eq 6a-c w = frequency, radltime w, = natural frequency of second-order component of process L i t e r a t u r e Cited Cohen, G. H.; Coon, G. A. Trans. ASME 1953, 75, 827. Coughanowr, D. R.; Koppel, L. B. "Process Systems Analysis and Control"; McGraw-Hill: New York, 1965. Harriott, P. "Process Control"; McGraw-Hill: New York, 1964. Latour, P. R.; Koppel, L. B.; Coughanowr, D. R . Ind. Eng. Chem. Process Des. Dev. 1967, 6 , 452. Lopez, A. M.; Miller, J. A.; Smith, C. L.; Murrill, P. W. Instrum. Techno/. 1967, 14, 57. McAvoy, T. J.; Johnson, E. F. Ind. Eng. Chem. Process Des. Dev. 1967, 6 , 440. Miller, J. A.; Lopez, A. M.; Smith, C. L.; Murrill, P. W. ControlEng. 1967, 14, 72. Murrill, P. W.; Smith, C. L. Hydrocarbon Process. 1966, 4 5 , 105 Smith, C. L.; Murrill, P. W. ISA J. 1966, 13, 50. Weber, T. W. "An Introduction to Process Dynamics and Control"; Wiley-Interscience: New York, 1973. Weber, T. W.; Bhalodia, M. Ind. Eng. Chem. Process Des. Dev. 1979, 18, 217. Weigand, W. A.; Kegerreis, J. E. Ind. Eng. Chem. Process Des. Dev. 1972, 11, 86. Yuwana, M.; Seborg, D. E. AIChE J. 1982, 28, 434 Ziegler, J.; Nichols, N. Trans. ASME 1942, 6 4 , 759.
Received for review June 8, 1984 Revised manuscript received December 5 , 1984 Accepted December 12, 1984
Infinite Dilution Activity Coefficients by Inert Gas Stripping Method: Extension to the Study of Viscous and Foaming Mixtures Domlnique Rlchon Ecole Nationale Sup6rieure des Mines de Paris, Centre R6acteur.s et Processus, Equipe de recherche associ6e au CNRS n O 768, Laboratoire de Thermodynamique, 77305 Fontainebleau, France
FranGols Sorrentlno and Andrde Vollley €cole Nationale Sup6rieure de Biologie Appliquh a la Nutrition et B L 'alimentation, Centre Universitaire Montmuzard, Laboratoire de Biologie Physico-Chimique, 2 1 100 DJon, France
The inert gas stripping method allows rapid determination of accurate infinite dilution activity coefficients. Unfortunately, previous works were limited to the study of nonfoaming mixtures with low viscosity (less than 50 cP) due to the design of the equilibrium cell. In this paper, a new cell design is proposed to extend the validity of the method up to 1000 cP. In addition, a special device was developed to break foams without disturbing phase equilibrium inside the equilibrium cell. Preliminary tests prove that the new apparatus is well suited to investigations of aqueous mixtures of polyols, glucides, and proteins.
Rapid determination of vapor-liquid equilibria is of great interest for industry: for example, in the selection of solvents to be used in extractive distillation, liquid extraction, or optimal conditions in food processing. In the
last case, one important question is how to provide consumers with the food having the best sensory qualities at the lowest prices. To answer this question, it is necessary to known the retention of aroma in food. A first approach
0196-4305/85/1124-1160$01.50/00 1985 American Chemical Society
