Relay with Additional Delay for Identification and Autotuning of

The paper presents a relay method for the identification of completely unknown processes for autotune purposes. It is an extension of a previous techn...
1 downloads 0 Views 127KB Size
Download PDF
Ind. Eng. Chem. Res. 1999, 38, 1987-1997

1987

Relay with Additional Delay for Identification and Autotuning of Completely Unknown Processes Claudio Scali,* Gabriele Marchetti, and Daniele Semino Department of Chemical Engineering, University of Pisa, Via Diotisalvi 2, 56126 Pisa, Italy

The paper presents a relay method for the identification of completely unknown processes for autotune purposes. It is an extension of a previous technique (ATV; Li, W.; et al. Ind. Eng. Chem. Res. 1991, 30, 1530), which assumed the delay of the process to be known. By means of a maximum of three relay tests, with additional delay, models with up to five parameters can be built. The proposed procedure does not present any convergence problem and is equivalent to the original one in terms of relative accuracy and duration of tests. The application of the identification procedure to a much wider set than in the original paper shows that the model obtained from the identification presents good accuracy in the high-frequency region, while some discrepancies of different nature may be present in the low frequency. When a suitable tuning method, which exploits this characteristic, is adopted, reasonably good closed-loop performance can be achieved for proportional-integral control in all cases. The low sensitivity to experimental errors and the simple implementation make the method interesting for application in industrial autotuners. 1. Introduction A process model may be unknown for several reasons including the lack of basic knowledge to build it from first principles, the complexity of the process, and its variability with the operating conditions. In these cases relay methods are a viable alternative for process identification; in fact, they allow one to build a low-order model, in a relatively simple and fast way; in its turn, the model can be used for the design of a conventional [proportional-integral (PI) or proportional-integralderivative (PID)] controller with the objective of giving acceptable performance. Being fast and easy, relay tests can be frequently repeated to detect variations in process parameters and to redesign the controller on the basis of the new model. This is why relay identification is generally incorporated in autotuning packages. The basic idea of the method (A° stro¨m and Ha¨gglund, 1984) is to bring the system into conditions of permanent oscillation by means of a relay which replaces the controller in the closed loop (Figure 1a). Under mild conditions required for application of the describing function analysis (Cook, 1986), including low-pass characteristics and the small influence of all harmonics higher than the first one, the ultimate gain (ku) and frequency (ωu) can be computed as a function of relay height (h), amplitude (a), and period (Pu) of the output oscillation as

ku )

4h 2π ; ωu ) πa Pu

Opposite to these advantages, relay methods can have two main drawbacks. The first one concerns the accuracy of the identification; the relay being a nonlinear element, the values of the critical parameters determined experimentally can be different with respect to the real values of the process. Li et al. (1991) report that in many cases errors up to 25% must be expected. Precision can be improved by adding other nonlinear elements (e.g., saturations, as illustrated by Friman and

Waller, 1995), thus increasing somewhat the complexity of the apparatus. The second limitation concerns the fact that from a limited knowledge of the process, as given by the parameters ku and ωu, the Ziegler and Nichols technique (ZN; Ziegler and Nichols, 1942) becomes an obliged choice for the tuning of the controller. However, this does not guarantee the stability of the closed-loop response under PI control, even for simple classes of processes. Improved tuning techniques have been introduced to enlarge the stability region. Among them are the following: Hang et al. (1991) propose a refinement of the ZN technique, which requires knowledge of the process gain; Shen and Yu (1994) adopt larger detuning factors, without any other additional knowledge on the process. Anyway, these techniques are not able to guarantee stability for completely unknown processes when the only available information is the critical point (ku, ωu) (Semino et al., 1996). Certainly, this fact is a limitation to the success of these techniques for autotuning application and for the design of decentralized controllers for multivariable processes by means of sequential relay tests, as illustrated by Loh et al. (1993) and by Shen and Yu (1994). In fact, in these cases the response of the process with some loops closed may become unstable with ZN tuning, even starting from individual processes with very simple dynamics. Also the approach via simultaneous multirelay tests (Friman and Waller, 1994; Palmor et al., 1995), which is based on the same two parameters, can be affected by the problem of not guaranteeing stability. This inconvenience can be eliminated by acquiring some more knowledge about the process, for example, by performing an identification in the third quadrant of the Nyquist plane. One appropriate point in the third quadrant, able to guarantee some gain and phase margin, can be obtained by means of a two-channel relay (Friman and Waller, 1997) or by means of a relay with hysteresis (Scali and Costaggiu, 1998). More points

10.1021/ie980616l CCC: $18.00 © 1999 American Chemical Society Published on Web 04/15/1999

1988 Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999

Figure 1. (a) Pure relay: scheme and relationships among parameters. (b) Relay with additional delay: scheme and representative points in a Nyquist plane.

can also be obtained by repeating experiments with additional delay elements in series to the relay (Figure lb), according to the autotune variations (ATV) method (Luyben, 1987; Li et al., 1991). In this case, a model of the process is built at the expense of a longer duration of tests and under the assumption of the process gain or time delay to be known; this model is then used to improve the closed-loop response. The scope of this paper is to extend the ATV method to processes for which no a priori knowledge is available. Related questions which will be addressed regard the accuracy of identification via ATV techniques for process modeling and the control design technique to be coupled with the identification method to improve closed-loop performance. In particular, the sections which follow deal with the new procedure for the extension of ATV techniques, identification results and model accuracy, closed-loop performance for PI control, and sensitivity to measurement errors and perturbations. 2. Extension of the ATV Identification Method The identification techniques known as ATV have been developed by Luyben and co-workers (1987, 1991) to identify stable processes by means of low-order models. For the case of highly nonlinear distillation columns, Luyben (1987) first proposed a method (ATV0) to identify the two time constants of a model up to third order by means of only one relay test; the gain and time delay must be known and are evaluated before beginning the relay tests. A modification of this method, for processes with known delay was presented by Li et al. (ATV; 1991): it requires two or three tests, according to the number of parameters to be estimated in the model. An application of this technique for identification and autotuning of multivariable processes is described in Semino and Scali (1998).

The knowledge of the gain and of the delay of the process can be obtained by steady-state data or by performance of a preliminary step test on the plant; nevertheless, in many cases these data are not promptly available or can vary during the operation; in addition, the necessity of these data goes against the principles of autotuning procedures which should be able to detect automatically changes in all process parameters. Motivated by these considerations, a further extension (ATV+) of the ATV method for a completely unknown process has been studied and is described in detail in this section. In particular, the differences with respect to the previous ATV technique, as proposed by Li et al. (1991), will be pointed out. The steps of the procedure are the following: 1. Get ultimate gain and frequency of the process from an autotune test. 2. Perform one or two more autotune tests with additional delay in order to identify gains and phases at some additional frequencies. 3. Elaborate the experimental data to identify the unknown parameters of different possible models. 4. Discriminate among the candidate models. Steps 1 and 2 of the procedure follow the guidelines of Li et al. (1991). The first test is the usual relay feedback test proposed by A° stro¨m and Ha¨gglund (1984). With the process being completely unknown, also the sign of the gain cannot be known; in this case, if the wrong sign of the relay is used in the first run, the system will not oscillate and a second run is required, after having changed the sign of the relay. The additional delays are chosen as ∆1 ) θ2 - θ1 ) 5π/12ωu, and ∆2 ) θ3 - θ1 ) 5π/4.8ωu as suggested by Li et al. (1991), who state that for most cases these delays correspond to two additional phase angles of ≈45° and ≈75°. In order to avoid misleading results caused by the transient oscillating behavior, amplitude and frequency are identified after two stabilized cycles have taken place. This situation has been identified as correspond-

Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999 1989

ing to two successive cycles, the latter of which has an amplitude and period which are no more than 5% different from the previous ones. Average values of the amplitude and period on the two cycles are computed. This somewhat arbitrary criterion has turned out not to put limits on the precision of the identification procedure. To illustrate the procedure and the involved computational aspects, six candidate models have been considered:

G1(s) ) G2(s) )

b0e-θ1s (3 par.) a1s + 1 b0e-θ1s

a2s2 + a1s + 1

(4 par.)

Therefore, a system of six equations coming from the three tests in the three unknowns b0, a1, and θ1 is obtained; this problem can be posed as a nonlinear leastsquares problem; i.e., find the values of the parameters which minimize the sum of squares of the differences between the left and right terms of the equations. However, the individual equations are nonlinear only in the delay; therefore, for a particular choice of the delay, the problem can be easily solved through a pseudoinversion. An estimate of the delay is found through the determination of minimum and maximum bounds. The minimum value for the delay is obviously 0. An estimate of the maximum value can be found by separating the contribution of the polynomial parts of the transfer function (φ(ωi)) to the one of the delay in the phase angles as follows:

φ(ω1) - θ1ω1 ) -π

-θ1s

G3(s) )

b0e 3

a3s + a2s2 + a1s + 1

Subtracting eq 17 from eq 16 and solving for θ1, one finds

-θ1s

G4(s) ) G5(s) ) G6(s) )

(b1s + b0)e a1s + 1 (b1s + b0)e

φ(ω2) - (θ1 + ∆1)ω2 ) -π

(5 par.)

(4 par.) θ1 )

-θ1s

a2s2 + a1s + 1

(5 par.)

(b1s + b0)e-θ1s a3s3 + a2s2 + a1s + 1

(6 par.)

This choice enlarges somewhat the one made by Li et al. (1991) by including the possibility of a process zero in the models of order 1-3. The reason why the transfer functions are written as a ratio of polynomials without explicitly showing gains and time constants will be apparent in the development that follows. For each model the number of parameters is explicitly written beside the model. The number of tests is directly related to the number of parameters: because each test supplies two data, the maximum number of identifiable parameters is twice the number of tests. Therefore, in the case of no initial guess about the order of the process, as is the rule for the case of lack of any knowledge on the process, three tests are required. It is necessary to state the problem of finding the process parameters from the available data appropriately. The procedure is explicitly written for model G1 (the extension to more complex models is obvious). The three tests supply the critical frequencies (ωi) and the real and imaginary coordinates of the critical points (Xi ) 1/kui, Yi ) 0) for the original process G(s) and for two modified processes G′(s) ) G(s) e-∆1s and G′′(s) ) G(s) e-∆2s. For each point, an equation can be written as

b0(cos(ωiθi) - i sin(ωiθi)) ) Xi + iYi ia1ωi + 1 Cross-multiplication and separation into real and imaginary parts give two equations for each point as

Xi ) a1ωiYi + b0 cos(ωiθi) Yi ) -a1ωiXi - b0 sin(ωiθi)

∆1ω2 + φ(ω1) - φ(ω2) ω1 - ω2

The term φ(ω1) - φ(ω2) is unknown, but being ω1 > ω2, its maximum value is equal to zero for models 1-3 and equal to π/2 for models 4-6. Therefore, a larger estimate of θ1 can be found as

θmax ) 1 ) θmax 1

∆1ω2 (models 1-3) ω1 - ω2

∆1ω2 + π/2 (models 4-6) ω1 - ω2

When the system is written in compact form as

z ) Ua where

[ ]

z ) [X1 Y1 X2 Y2 X3 Y3]T ω1Y1 -ω1X1 ω2Y2 U) -ω2X2 ω3Y3 -ω3X3

cos(ω1θ1) -sin(ω1θ1) cos(ω2θ2) -sin(ω2θ2) cos(ω3θ3) -sin(ω3θ3)

a ) [a1 b0]T

with θ1 fixed, the linear least-squares solution is

a ) [UTU]-1UTz The computation procedure is briefly illustrated in the sequel. For a given model, the sum of squares as a function is of the delay θ1 in the interval 0 < θ1 < θmax 1 computed; calculations are repeated for the different candidate models. As a second step, the discrimination among the models is accomplished according to the value of the

1990 Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999

Figure 2. Residual error as a function of the delay θ1 for different order models.

residual error. The value of the delay for which the pseudoinversion problem gives the lowest value of the sum of squares errors inside the interval is chosen. Unstable models are promptly discarded, as open-loop stability is assumed, at this time, as a base hypothesis for the application of the method; the extension to an open-loop gain stabilizable method is possible with some changes in the procedure. For model order up to 2, the condition is verified by imposing that coefficients a1 and

a2 are positive; for order 3 models, poles are explicitly computed. Models with a wrong sign of the gain (b0), as compared to the one determined by the sign of the relay, are discarded as well. When two models with the same number of parameters are compared, the one with the smaller error is preferred. However, if the number of parameters is different, the more complex model is preferred only if the reduction of the residual error is significant (1 order of magnitude or more). This is illustrated in Figure 2: in this case the model G2, indicating a value of θ1 ) 1.1, is chosen, even if a model G5 would indicate a smaller square error for a different value of θ1. It is relevant to point out that the computation time is absolutely negligible, if compared with the time required for the relay experiment. It might even be reduced by adopting a minimum search routine; this is not advised though, as in general the presence of local minima cannot be excluded. By operating this way, in all of the tests that have been accomplished no computational problems have been detected (Marchetti, 1998). The proposed method therefore avoids the convergence problems of iterative methods that do not assume the knowledge of the delay, which had been anticipated by Li et al. (1991); these problems motivated their choice of taking the delay to be known. Contrary to that, the proposed identification is accomplished with no previous knowledge. It is clear, however, that, were the delay known, the method would be exactly the same as the one by Li et al. (1991).

Table 1. Examined Processes I

n ) 1-3

G(s) )

ke-θs τs + 1

II

n ) 4-9

G(s) )

ke-θs (τ1s + 1)(τ2s + 1)

III

n ) 10-15

G(s) )

ke-θs τ s + 2τξs + 1

IV

n ) 16-18

G(s) )

ke-θs (τ1s + 1)(τ2s + 1)(τ3s + 1)

V

n ) 19-27

G(s) )

ke-θs (τ1s + 1)(τ22s2 + 2τ2ξs + 1)

VI

n ) 28-39

G(s) )

k(Rs + 1)e-θs (τ1s + 1)(τ2s + 1)

VII

n ) 40-45

G(s) )

2 2

k(Rs + 1)e-θs 2 2

τ s + 2τξs + 1

k)1 θ ) 0.1, 1, 10 τ)1 k)1 θ ) 0.1, 1, 10 τ1 ) 1 τ2 ) 0.8, 0.2 k)1 θ ) 0.1, 1, 10 τ)1 ξ ) 0.5, 0.25 k)1 θ ) 0.1, 1, 10 τ1 ) 1 τ2 ) 0.9 τ3 ) 0.8 k)1 θ ) 0.1, 1, 10 τ1 ) 0.1, 1 τ2 ) 0.1, 1 ξ ) 0.25 k)1 θ ) 0.1, 1, 10 τ1 ) 1 τ2 ) 0.8 R ) 0.3, 3, -0.8, -3 k)1 θ ) 0.1, 1, 10 τ)1 ξ ) 0.5 R ) 3, -3

Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999 1991

As a last remark, for sensitivity reasons it is advised to keep the number of equations always larger at least by one than the number of parameters. For this reason models with up to five parameters have been considered here (models G1-G5) for ATV+, while the model G6 has been used only to compare results with the case of known delay (ATV). The accuracy of the model for identification and control purposes and the sensitivity to errors are addressed in the next sections. 3. Identification Results and Model Accuracy A large number of processes have been examined in order to cover a wide range of dynamic behaviors: first, second, and third order, with time delay; overdamped and underdamped; with one negative or positive zero. The seven categories of processes (I-VII) with the parameter variations listed in Table 1 originate a total number of 45 meaningful cases, referred by their number in the table of results. The cases examined here enlarge significantly those investigated in Li et al. (1991), where only six overdamped cases (up to 5th order for the process and up to 3rd order for the model) have been studied. A first comparison of results limited to these six cases has been presented by Marchetti et al. (1998); comments and conclusions presented here generalize those first results. In order to test the ATV+ identification method, each of the 45 cases has been investigated by performing three relay tests; the resulting experimental data supply six equations for the five candidate models G1-G5 (model G6 is a candidate only for the ATV technique). It is worthwhile mentioning that two relay tests could be sufficient to identify low-order models with both techniques. However, because it is desirable to consider the more general case of no previous knowledge whatsoever, three tests are performed in all cases to include the complete spectrum of possible models in the identification procedure. Identification results are summarized in Table 2; some observations follow. When columns 3-5 of Table 2 are compared, it can be noted that the order of the model and the number of parameters identified by means of the ATV or ATV+ procedures can be different from the process. Generally, the model is lower order (25/45 cases for ATV+ and 16/ 45 for ATV); this happens when the process time constants have different orders of magnitude. In a few cases it is higher order (4/45 for ATV+ and 2/45 for ATV); this is caused by errors in the experimental data obtained via relay, which would indicate a nonmonotonous trend in G(ω) requiring higher order to be matched. In some cases the type of model is different from the process, even though the number of parameters is the same; for instance, an error on the time delay can be compensated by introducing some inverse response element. ATV+ in general, does not identify the same value of delay as ATV (for which it is known); only in the hypothetical case of perfect data would the two methods give exactly the same results. However, given that some identification error is inherent in the relay test, the possibility of changing the delay in order to fit the data may be beneficial. More precise data could be identified with sinusoidal set-point variations (Marchetti, 1998) but at the expenses of longer times for the tests and higher sensitivity to perturbations.

Table 2. Summary of Identification Results: Number of Parameters and Values of the IME Index for ATV and ATV+ for the 45 Examined Cases no. of params

IME

cat.

no.

G

G ˜ ATV

G ˜ ATV+

G ˜ ATV

G ˜ ATV+

I I I II II II II II II III III III III III III IV IV IV V V V V V V V V V VI VI VI VI VI VI VI VI VI VI VI VI VII VII VII VII VII VII

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45

3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5

3 3 3 4 4 4 4 4 4 3 4 5 3 4 4 4 5 5 4 3 5 4 5 5 4 5 5 4 5 5 3 3 3 4 5 6 4 5 5 3 3 3 5 5 5

3 5 3 3 4 4 3 5 4 3 4 4 3 5 5 3 4 5 3 3 5 5 5 5 3 5 4 4 4 5 3 3 3 3 3 5 3 3 5 3 3 3 3 4 3

0.204 0.141 0.210 0.196 0.085 0.053 0.204 0.116 0.067 0.071 0.039 0.044 0.456 0.022 0.022 0.190 0.074 0.041 0.206 0.187 0.127 0.217 0.043 0.020 0.456 0.039 0.046 0.197 0.124 0.202 0.521 0.221 0.214 0.101 0.188 0.340 1.217 0.598 0.500 1.276 0.592 0.280 1.476 0.488 0.296

0.205 0.161 0.211 0.196 0.086 0.054 0.204 0.125 0.060 0.057 0.057 0.044 0.457 0.061 0.107 0.190 0.068 0.039 0.206 0.187 0.149 0.210 0.080 0.021 0.458 0.063 0.053 0.198 0.122 0.191 0.521 0.210 0.214 0.104 0.213 0.403 1.225 0.720 0.913 1.276 0.531 0.279 1.484 0.629 0.700

The comparison of the number of parameters is not in any case a reliable criterion for fitness of identification. However, it is difficult to define a synthetic (numeric) index to evaluate at a glance the accuracy of the identification, in an absolute sense (referring to the true process G) and in a relative sense (comparing the models G ˜ from ATV and ATV+). The error on single parameters is not meaningful, as both identification techniques change the values of parameters in order to compensate for the different orders between the models and the process. In columns 6 and 7 of the table, values of the index IME are reported, where the index is defined as the integral of the multiplicative error in the frequency domain:

IME )

G ˜ (ω) - G(ω) dω G(ω)

∫0ω*

The higher frequency (ω*) in the range of integration is defined such that φ(ω*) ) -3π/2, that is, the frequency at which the Nyquist plot of the identified process reaches the positive imaginary axis. This index

1992 Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999

Figure 3. Four typical frequency behaviors: G, s; G ˜ ATV, - -; G ˜ ATV+, - - -.

takes into account the error caused by the approximation and therefore allows one to compare the relative merits of the two methods; the low-frequency error has a smaller weight, as the two techniques identify in the high-frequency range; in addition, low-frequency errors can be considered less important in the case of an identification for control purposes, as will be shown later. When the values of the index for the two techniques (columns 6 and 7) are compared, it can be seen that the difference is contained within 10% for 31/45 cases and is larger than 30% only for 6/45. On this basis, it can be concluded that the two methods give practically the same error with respect to the true process G. A more complete comparison can be obtained through the analysis of the behavior of the three systems G, ˜ ATV+ in the frequency domain. The appearance G ˜ ATV, G of some typical behaviors allows one to integrate previous considerations based on the (subjective) index IME. In particular, four classes of processes can be identified. The typical behavior of a representative of each class is reported in Figure 3. 1. The first class of processes (14/45 cases) shows negligible differences among the three functions, thus indicating a very good identification. This happens both for underdamped and for overdamped systems having comparable values of delay and time constants.

2. A second class of processes (9/45 cases) shows similar behavior in the identification in the higher ˜ ATV(ω) = G ˜ ATV+(ω)), frequency range (ω > ω1: G(ω) = G while at low frequency the true process differs from the two models which behave similarly (ω f 0: G(ω) * ˜ ATV+(ω)). This is observed mainly for overG ˜ ATV(ω) = G damped systems having time constants of different magnitude: in the approximation due to identification, fast dynamics are filtered and the process gain is changed accordingly in order to compensate for the different dynamics, thus originating the low-frequency error. 3. For a third class of processes (16/45 cases), the two identification techniques give very similar results in the entire frequency range, but there is a bias with respect ˜ ATV+(ω). to the true process (∀ ω: G(ω) * G ˜ ATV(ω) = G This happens for the case of process dynamics dominated by a large delay or a positive zero. It can be explained by considering that both techniques are affected by the errors which are inherent in the relay method; in particular, the dominant element of the dynamics imposes a frequency behavior such that neglect of the high harmonics at the identification frequencies brings systematic errors in the amplitude, which are compensated by a variation of the gain. 4. For a fourth class of processes (6/45 cases), the two identified models are similar only in the investigated

Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999 1993

Figure 4. Set-point responses for the four classes of processes with PI control based on three models: G, s; G ˜ ATV, - -; G ˜ ATV+, - - -.

range (ω3-ω1) and differ from the process in most of the frequency range. This is observed, in general, with processes having complex dynamics (categories VI and VII only). When the previous classification, which must be seen as a simplification of all of the possible situations (some processes show intermediate behaviors), is referred to, two partial conclusions can be drawn at this point: (1) The two methods ATV and ATV+ give very similar models for 39/45 cases, despite the fact that the delay is unknown for ATV+. (2) The identification by means of relay techniques is less accurate for delay dominant processes and for inverse response (large zeroes) processes. As illustrated below, inconveniences due to identification errors can be in large part compensated by adopting an appropriate design of the controller.

should be able to outperform the lower order PI. It has been pointed out that the model obtained by ATV methods in many cases is affected by error in the lowfrequency region: this implies that an advanced controller must be detuned to guarantee robustness. As a consequence, performance becomes comparable with respect to the PI control. The PI algorithm can be considered as a reduction of a higher order structure, derived from the model by approximating the multiplicative delay by means of Pade` polynomial and neglecting terms of order superior to 1. This fact is clearly pointed out by adopting the internal model control structure (IMC; Morari and Zafiriou, 1989) as a design framework. For a process having a number n of poles pi (real or complex) and a multiplicative delay θ, the IMC design for step inputs leads to a PI controller with the integral time constant

4. Closed-Loop Results In the perspective of applying the proposed identification method in an autotuning package, a PI controller has been taken into consideration. This choice might not seem in agreement with the efforts to obtain a model of the process because, as a matter of principle, a model-based advanced controller

n

τI )

∑1 (-1/pi) + θ/2

The proportional constant kc is correlated with the speed of the closed-loop response or, equivalently, with the robustness of the control system to unavoidable

1994 Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999

Figure 5. Disturbance suppression for the four classes of processes with PI control based on three models: G, s; G ˜ ATV, - -; G ˜ ATV+, - - -.

errors on the model. Following considerations regarding the optimization of performance of conventional controllers, as given by Lee et al. (1990), the proportional gain (kc) can be set as the result of the maximization problem:

kc: max [|η(ω)|] ) 1.26 where η ) GC/(1 + GC) is the closed-loop sensitivity function of the system. When this condition is imposed, the closed-loop response generally will be slightly underdamped with an overshoot of about 20% and a relatively fast initial response; this corresponds to excitation of the system in the high-frequency range where the derived model has been identified more accurately. The closed-loop response has been analyzed for the 45 cases of Table 2; for each case a PI regulator designed, as illustrated above, on the basis of the two ˜ ATV+ obtained by the two identificamodels G ˜ ATV and G tion techniques has been compared with a design based on the true process G. Results for the four typical processes, which have been analyzed by the frequency behavior, are reported in Figure 4 for set-point tracking and in Figure 5 for the suppression of a step disturbance acting at the process input.

It can be seen that for class 1 and 2 processes, the three responses are almost coincident; on the contrary, some differences can appear for class 3 and 4 processes. ˜ ATV+ are The responses obtained with models G ˜ ATV and G very much similar for all of the cases. These can be considered as general results and are confirmed by analyzing columns 2 vs 3 and 6 vs 7 of Table 3, where performances for the two cases of setpoint tracking and disturbance suppression are evaluated by means of the difference of the integral of square errors (DISE), defined as

DISE ) 100

∫0∞eG˜ 2 dt - ∫0∞eG2 dt ∫0∞eG2 dt

where eG2 indicates the closed-loop square error when the control system has been designed on the true process G, while eG˜ 2 refers to the same error when the control system has been designed either on model G ˜ ATV or on model G ˜ ATV+. The DISE index can be considered appropriate for a comparison of absolute and relative merits of the two techniques. As far as the relative merits are concerned, the two columns of Table 3 give values which are very close to each other to indicate equivalent performances

Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999 1995 Table 3. Summary of Closed-Loop Results: Values of the DISE Index for Set-Point Tracking and Disturbance Rejection for a Control System Based on the Model Identified via ATV and ATV+ Techniques, in the Nominal Case and in the Presence of Measurement Errors set-point change nominal

disturbance suppression

20% perturbation

nominal

20% perturbation

cat.

no.

G ˜ ATV

G ˜ ATV

G ˜ ATV

G ˜ ATV

G ˜ ATV

G ˜ ATV

G ˜ ATV

G ˜ ATV+

I I I II II II II II II III III III III III III IV IV IV V V V V V V VI VI VI VI VI VI VI VI VI VI VI VI VI VI VI VII VII VII VII VII VII

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45

9.52 2.19 4.42 9.74 1.47 8.04 10.32 1.58 6.73 0.19 -1.01 4.58 6.97 -0.91 0.75 8.96 2.29 3.52 8.96 2.42 -0.88 1.90 1.96 -0.10 4.35 0.54 -1.17 9.90 1.82 1.32 10.57 -38.67 -18.72 1.47 3.11 2.06 33.73 12.99 8.40 31.19 -16.91 9.06 15.02 12.84 7.31

9.93 3.47 4.31 9.07 1.43 7.91 9.70 2.85 7.12 0.41 -1.13 3.96 5.74 4.89 -3.71 8.56 1.29 3.98 10.22 2.32 0.98 2.20 8.62 0.41 4.04 5.20 -0.47 9.72 2.40 5.67 10.50 -38.78 -18.73 1.23 2.64 1.47 33.96 12.62 7.16 30.93 -19.83 8.86 14.04 14.64 7.04

2.21 2.64 5.24 8.18 5.13 5.84 1.89 -0.04 0.55 0.24 -1.46 3.79 9.15 -2.34 -0.41 7.14 7.42 49.27 4.92 1.61 26.83 1.96 -8.65 -4.61 3.65 -8.12 10.86 1.52 -0.85 21.73 12.91 -38.66 -18.76 2.65 13.98 12.93 58.86 12.97 18.79 37.50 -18.40 25.81 26.60 25.07 9.46

3.16 2.03 8.60 9.89 4.67 8.54 1.81 0.16 0.85 0.59 0.58 2.46 10.14 -2.40 0.05 9.75 5.28 12.14 4.12 3.02 29.86 2.48 30.62 -5.61 4.71 -8.12 8.93 0.37 1.40 7.47 12.73 -38.47 -19.10 2.75 14.03 2.56 50.10 19.53 13.41 64.26 -20.22 26.89 20.40 11.93 54.05

8.91 1.79 22.93 9.14 0.84 9.90 9.95 -0.20 7.12 0.12 -0.71 9.11 5.70 0.20 1.01 8.02 0.61 6.44 8.33 -1.10 -27.58 1.57 1.73 -0.15 3.19 0.52 -0.42 9.15 2.96 51.12 2.70 -17.41 -46.30 1.69 6.42 8.55 26.19 -2.39 -7.91 -17.97 1.01 52.89 10.98 -7.08 -11.40

9.44 3.27 23.68 8.45 0.73 9.55 9.33 0.91 8.28 0.30 -1.06 7.70 4.28 1.60 -2.97 7.62 -0.26 7.03 9.70 -1.16 -33.52 1.97 6.65 0.08 3.12 2.27 1.94 8.97 3.06 34.56 2.62 -16.37 -46.28 1.43 5.75 7.71 26.40 -2.33 -7.42 -18.10 4.22 51.95 9.96 -3.94 -2.80

1.84 9.76 41.28 7.33 7.83 -2.09 1.58 -2.39 -20.16 0.16 -1.49 6.77 7.44 0.24 -0.22 5.85 6.97 87.86 4.40 5.30 -34.83 1.84 -5.98 -0.18 3.34 5.82 4.67 0.97 -6.78 60.01 4.30 -17.47 -57.31 3.06 22.82 18.11 47.77 -2.51 -3.88 -15.48 1.40 -22.71 20.35 1.36 -12.68

2.89 10.04 31.44 9.25 4.19 2.04 1.51 -2.49 -19.12 0.46 -2.24 3.90 7.86 0.32 0.12 8.77 -1.45 19.96 3.51 6.02 -35.35 2.11 25.98 -0.27 4.59 6.17 2.53 -0.27 -4.19 50.87 4.15 -18.25 -57.22 3.17 22.95 2.89 40.16 1.49 -6.20 -1.53 3.48 -22.87 14.65 -3.82 3.63

+

for the two identification techniques. In regards to absolute accuracy, it can be pointed out that, for 35/45 cases DISE e +10% to confirm equivalence also with respect to a control system based on the true process; only for 2/45 cases, DISE g +30% (to be pointed out that negative values indicate smaller errors for the controllers based on the identified models than those based on the true process, because of the suboptimality of the PI algorithm). After examination of closed-loop results, it can be concluded that equivalent performances are obtained by the two controllers designed on the basis of the identified models (ATV and ATV+), as expected from the almost overlapping Nyquist plots (Figure 3). The fact that in most cases closed-loop performance is very similar for the three controllers based on different models confirms that the design method is very effective in exploiting the accuracy of the model in the frequency range of interest. Therefore, identification techniques based on relay, even if not accurate in the entire frequency range, are a valid tool to obtain good models for control purposes,

+

+

when standard PI/PID controllers are used. Reasonably enough, a model-based control technique could point out more clearly the pitfalls of the identification; therefore, whenever a higher performance is desired to motivate the choice of an advanced control system, coupling with an appropriate (more accurate) identification technique becomes of primary importance. 5. Sensitivity to Measurement Errors and Perturbations Previous results were obtained in the hypothesis that the only source of errors were the approximation inherent to relay techniques. Hints have been given for more sophisticated methods that eliminate this mismatch. Here the presence of possible perturbations on experimental values is considered through simulations where a disturbance able to influence the amplitude and the frequency of the signal analyzed for identification is added. Random disturbances, with variable amplitude up to 30% around the average value of amplitude of the output signal and filtered by first-order dynamics, were

1996 Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999

Figure 6. Perturbation on the output signal (a) and its effect on the experimental value (b) (nominal, - - -; perturbed, s).

used for this scope. An example of such a disturbance and its effect on the output signal are reported in Figure 6. To minimize the effect of the presence of perturbations, some further attention must be paid in the transient period to detect the appearance of stable limit cycles. In this case, in order to comply with the abovestated criterion for detection of stable cycles, it may be necessary to observe more than two oscillations, thus increasing the test time. In any case, it is advised to make calculations on the basis of the average values of at (least) two stable oscillation cycles. In Table 3 values of DISE for the closed-loop response of controllers based on the models obtained by the two identification methods are reported also for the presence of a random disturbance on the output signal during identification, as described above. When the value of DISE with respect to the situation without perturbation is compared on each row, it can be seen that there are very few cases in which there is a relevant variation of the integral square error. In particular, only in 12/45 cases for y(r) and in 10/45 cases for y(d) is the change of the value of DISE larger than 10%; most of the cases are the same for the two control objectives and for the two identification methods. When we concentrate on the analysis of the set-point tracking problem, (a) the 12/45 cases refer to process numbers 18, 21, 23, 27, 30, 35, 37, 40, 42 43, 44, and 45; (b) in cases 23 and 45, ATV deteriorates less than ATV+ and vice versa for cases 18, 30, and 33; and (c) in all of the other cases, the two methods behave in the same way. Therefore, it can be concluded that by operating this way the effect of the perturbation on closed-loop performance is not very significant for both techniques. 6. Conclusions The proposed method (ATV+) shows to be effective in the identification, for autotuning purposes, of processes without any a priori knowledge on the number and values of model parameters. The analysis of a large number of process dynamics (including systems up to third order, with time delay and inverse response) allows one to conclude that the ATV+ gives equivalent accuracy with respect to the previous ATV technique, which assumed the process delay to be known. Models with up to five parameters can be identified by means of three relay experiments, thus allowing one to acquire more information about the process dynamics

with respect to the simple relay. Of course, this supplementary knowledge is paid with longer times spent for the relay tests. The identified model presents a good accuracy in the high-frequency region, while some discrepancies with the real process can appear in the low-frequency region. In any case, an excellent closed-loop performance can be obtained for PI controllers by adopting a design technique in which the integral time is correlated with the time constants and delay of the identified process, while the proportional gain is computed to give a desired maximum value of the complementary sensitivity function. This design shows to be appropriate to exploit the characteristics of the identification. The procedure to estimate the value of the delay and then the model of the process does not present any convergence problem and shows to be robust to errors deriving from perturbations on experimental measurements. Therefore, the method can be proposed for application in industrial autotuners operating on completely unknown processes. Nomenclature a ) amplitude of output oscillations a1, a2, a3 ) coefficients of the denominator a ) parameter vector b0, b1 ) coefficients of the numerator C ) controller d ) disturbance e ) error G, G′, G′′ ) transfer functions h ) amplitude of the relay k ) gain Pu ) period of the oscillation p ) poles of the system r ) set point s ) Laplace variable u ) manipulated variable U ) matrix of coefficients in the LSM solution z ) vector of coordinates in the LSM solution Xi ) real coordinate Yi ) imaginary coordinate y ) output response Greek Symbols R ) time constant ∆ ) additional delay θ ) delay φ ) phase lag

Ind. Eng. Chem. Res., Vol. 38, No. 5, 1999 1997 ω ) frequency τ ) time constant ξ ) damping coefficient Pedices ATV ) modeled via the ATV technique ATV+ ) modeled via the ATV+ technique c ) proportional 1, 2, 3 ) identified points I ) integral i ) generic element u ) ultimate Indices max ) maximum value T ) transpose

Literature Cited A° kstro¨m, K. J.; Ha¨gglund, T. Automatic Tuning of Simple Regulators with Specifications on Phase and Amplitude Margins. Automatica 1984, 20, 645. Cook, P. Non-Linear Dynamical Systems; Prentice Hall: Englewood Cliffs, NJ, 1985; pp 52-64. Friman, M.; Waller, K. V. Autotuning of Multiloop Control Systems. Ind. Eng. Chem. Res. 1994, 33, 1707. Friman, M.; Waller, K. V. Closed-Loop Identification by Use of Single-Valued Nonlinearities. Ind. Eng. Chem. Res. 1995, 34, 3052. Friman, M.; Waller, K. V. A Two-channel Relay for Autotuning. Ind. Eng. Chem. Res. 1997, 37, 2662. Hang, C. C.; A° stro¨m, K. J.; Ho, W. K. Refinements of the Ziegler Nichols Tuning Formula. IEEE Proc., Part D 1991, 138 (March), 111. Lee, J.; Cho, W.; Edgar, T. F. An Improved Technique for PID Controller Tuning from Closed-Loop Tests. AIChE J. 1990, 36, (12), 1891.

Li, W.; Eskinat, E.; Luyben, W. L. An Improved Autotune Identification Method. Ind. Eng. Chem. Res. 1991, 30, 1530. Loh, A. P.; Hang, C. C.; Quek, C. K.; Vasnani, V. U. Autotuning of Multiloop Proportional-Integral Controllers Using Relay Feedback. Ind. Eng. Chem. Res. 1993, 32, 1102. Luyben, W. L. Derivation of Transfer Functions for Highly Nonlinear Distillation Columns. Ind. Eng. Chem. Res. 1987, 26, 2490. Marchetti, G. Tesi di Laurea (in Italian); Department of Chemical Engineering, University of Pisa: Pisa, Italy, 1998. Marchetti, G.; Semino, D.; Scali, C. Identification of Unknown Processes for Autotune Purpose by Relay Tests. Proc. MCEA ‘98: International Conference on Elettronics and Applications; Marrakesh, Morocco, 1998; p 454. Morari, M.; Zafiriou, E. Robust Process Control; Prentice Hall: Englewood Cliffs, NJ, 1989; pp 114-124. Palmor, Z. J.; Halevi, Y.; Krasney, N. Automatic Tuning of Decentralized PID Controllers for TITO Processes. Automatica 1995, 31, 1001. Scali, C.; Costaggiu, R. Relay with Hysteresis for Monitoring and Controller Design. Proceedings of DYCOPS-5, IFAC Symposium: Dynamics and Control of Process Systems; Corfu, Greece, 1998; p 596. Semino, D.; Scali, C. Improved Identification and Autotuning of PI Controllers for MIMO Processes by Relay Techniques. J. Process Control 1998, 8 (3), 219. Semino, D.; Mazzanti, I.; Scali, C. Design of Decentralized Controllers by a Relay Technique: Extensions of Tuning Rules. Proc. UKACC Int. Conf. CONTROL ‘96, 1996, 1090. Shen, S. H.; Yu, C. C. Use of Relay-Feedback Test or Automatic Tuning of Multivariable Systems. AIChE J. 1994, 40, 627. Ziegler, J. G.; Nichols, N. B. Optimum Setting for Automatic Controllers. Trans. ASME 1942, 65, 433.

Received for review September 24, 1998 Revised manuscript received February 23, 1999 Accepted February 24, 1999 IE980616L