Ind. Eng. Chem. Res. 2010, 49, 7807–7813
7807
Area Methods for Relay Feedback Tests Jietae Lee* and Su Whan Sung Department of Chemical Engineering, Kyungpook National UniVersity, Taegu 702-701, Korea
Thomas F. Edgar Department of Chemical Engineering, UniVersity of Texas, Austin, TX78712
The amplitude and period of the relay feedback oscillation are used to obtain the ultimate parameter of a dynamic process and to tune a proportional-integral-derivative (PID) controller automatically. Equations for the ultimate gain and period are based on ignoring higher harmonic terms in the relay feedback response. The integral of the relay feedback response can suppress the higher harmonic terms, and the amplitude of this integral has been found to be better for estimating the ultimate gain of a process. As extensions, simpler methods that use several area calculations of the relay feedback response are proposed and shown to provide more accurate frequency responses, ultimate gain and period, and parametric models. Introduction A feedback system where a relay is put in the feedback loop will oscillate. The oscillation period is close to the ultimate period when the higher harmonic terms are attenuated sufficiently by the process dynamics. The ultimate gain can also be obtained by measuring the amplitude of oscillation. The ultimate gain and period are used to tune proportional-integralderivative (PID) controllers.1 This relay feedback system can be easily implemented in field with minimal efforts and is one of the standard methods to tune PID controllers automatically.2,3 When attenuation of the higher harmonic terms by the process dynamics is not sufficient, there exist considerable errors in estimation of the ultimate gain and period. Several methods such as a saturation relay,3 relay with a P control preload,4 and a two level relay5 were introduced to obtain more accurate ultimate gain and period of the process by suppressing the high order harmonic terms in the relay output. Asymmetric relay oscillations due to unknown load disturbances cause additional errors in the ultimate parameter estimations. Methods to reject unknown load disturbances and restore symmetric relay oscillations have been proposed.6-8 Parametric models from several measurements of the relay feedback oscillation including the oscillation amplitude can be obtained.9-11 For example, Luyben9 used the shape factor to extract a three-parameter model from the cyclic steady state part of the relay response. However, methods that use only the cyclic steady state part cannot provide acceptable robustness for uncertainty such as process/model mismatch and nonlinearity. They may provide poor model parameters such as negative gain when the model structure is different from that of the process.11 A relay experiment with a subsequent P control experiment or another relay feedback test can be used to obtain a parametric model robustly.12,13 Huang et al.14 used the integral of the relay transients to obtain the steady state gain of the process. Recently, Lee et al.15 proposed methods to extract more accurate frequency response data and parametric models from the response of a single conventional relay feedback test and its integral. Integrals of the relay responses make the fundamental frequency term more dominant compared to the higher harmonic terms, resulting in * To whom correspondence should be addressed. Tel.: +82-53-9505620. Fax: +82-53-950-6615. E-mail:
[email protected].
better accuracy in estimating frequency responses and model parameters. This paper extends the previous method of integral of the relay feedback response and proposes methods using integrals of the conventional relay feedback responses. These modifications enhance identification performances without modifying the relay feedback system. Accurate estimates of frequency responses of the first and third harmonics, ultimate gain/period and parametric models are obtained. Since the conventional relay feedback method is used, the proposed methods share the practical and theoretical merits of the conventional relay feedback method. Because it is not necessary to store whole trajectories and the computations are simple, the proposed methods can be implemented easily as an enhancement of standard commercial PID controllers. Conventional Relay Feedback Method. A classical relay feedback system shown in Figure 1a is considered. This relay feedback system will produce a stable oscillation. Astrom and Hagglund1 used this oscillation to extract approximate ultimate data and tune PID controllers automatically. The approximate ultimate period Pu and ultimate gain Kcu are Pu ) p 4h Kcu ) πay
Figure 1. Relay feedback system and responses.
10.1021/ie901546j 2010 American Chemical Society Published on Web 11/20/2009
(1)
7808
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
where p is the period of oscillation and ay is the measured amplitude of process output. The ultimate period and gain in eq 1 are approximate due to ignoring the high harmonic terms in the process output. As a result, the ultimate period and gain of eq 1 show relative errors up to 5% and 18%, respectively, for first-order plus time delay (FOPTD) processes, which may not be acceptable. Let the input and output trajectories which are fully developed (cyclic steady state) be u(t) and y(t) for the conventional relay feedback system, respectively (Figure 1b). The relay output (a square wave) can be represented by a Fourier series17 1 4h 1 sin(ωt) + sin(3ωt) + sin(5ωt) + · · · , π 3 5 ω ) 2π/p(2) where h is the relay amplitude and p is the period of relay feedback oscillation. The output corresponding to u(t) of eq 2 is16
(
u(t) )
)
4h |G(jω)| sin(ωt + ∠G(jω)) + π 1 |G(j3ω)| sin(3ωt + ∠G(j3ω)) + · · · (3) 3 where G(s) is the process transfer function. From eq 3,
(
y(t) )
)
4h y(p/4) ) |G(jω)| cos(∠G(jω)) π 1 |G(j3ω)| cos(∠G(j3ω)) + · · · ≈ ay 3 4h y(0) ) |G(jω)| sin(∠G(jω)) + π 1 |G(j3ω)| sin(∠G(j3ω)) + · · · ) 0(4) 3 By neglecting the high harmonic terms, eq 4 leads to eq 1, the estimates of ultimate gain/period of G(s). Methods such as a saturation relay,3 relay with a P control preload,4 and a two level relay5 use u(t) with smaller higher harmonic terms to reduce errors in eq 1. The integral of y(t) reduces the higher harmonic terms in y(t), and methods by Lee et al.15 use this. In fact, we can find exact G(jω) from the two integrals of ∫p0y(t) cos(ωt) dt and ∫p0y(t) sin(ωt) dt.16 Here, simple methods that approximate these integrals are proposed. Areas of Relay Oscillation. First the following two areas as shown in Figure 1b are considered
4h |G(jω)| sin(ωt + ∠G(jω)) + π 1 |G(j3ω)| sin(3ωt + ∠G(j3ω)) + · · · 3 4h ) |G(jω)| cos(∠G(jω)) sin(ωt) + π |G(jω)| sin(∠G(jω))cos(ωt) + 1 |G(j3ω)| cos(∠G(j3ω)) sin(3ωt) + 3 1 |G(j3ω)| sin(∠G(j3ω)) cos(3ωt) + · · · )(8) 3 we have (Appendix)17
(
y(t) )
)
(
8hp 1 |G(jω)| cos(∠G(jω)) + |G(j3ω)| cos(∠G(j3ω)) + · · · 9 π2 8hp 1 q2 ) 2 |G(jω)| sin(∠G(jω)) - |G(j3ω)| sin(∠G(j3ω)) + · · · 9 π
( (
q1 )
(9) In addition to the above areas, an additional two areas as shown in Figure 1b are considered for better estimates
∫ D) ∫
∫ ∫
0 p/2
B)
p/4
∫
p
0
q(t)y(t) dt ) 2(
) 2(A + B) q2 )
∫
p
0
(
qt-
) 2(A - B)
∫
p/4
0
p y(t) dt ) 2( 4
)
{
q3 )
y(t) dt
(5)
p/4
∫
p/2
p/4
y(t) dt -
y(t) dt)
r(t)y(t) dt ) 2(
∫
∫
p/2
p/4
q4 )
p/4
0
ty(t) dt +
(p/2 - t)y(t) dt) ) 2(C - D) + pB
∫ r(t - 4p )y(t) dt ) 2( ∫ (p/4 - t)y(t) dt ∫ (t - p/4)y(t) dt) ) -2(C + D) + 2p (A + B)(12) p
p/4
0
0
p/2
The Fourier series of r(t) is17
∫
p/2
p/4
2p 1 sin(ωt) - sin(3ωt) + · · · 9 π2
(
)
(13)
Hence we have
y(t) dt)
4hp2 1 (|G(jω)| cos(∠G(jω)) |G(j3ω)| cos(∠G(j3ω) + · · · ) 27 π3 4hp2 1 q4 ) 3 (|G(jω)| sin(∠G(jω)) + |G(j3ω)| sin(∠G(j3ω) + · · · ) 27 π
q3 )
where (Table 1)
{
p
0
y(t) dt
(6)
1, t ∈ (0, p/2) -1, t ∈ (p/2, p) (7) q(t ( p) ) q(t) The Fourier series of the square wave q(t) is q(t) ) (4/π)(sin(ωt) + (1/3) sin(3ωt) + ...).17 Since q(t - p/4) ) (4/π)(cos(ωt) - (1/3) cos(3ωt) + ...) and q(t) )
∫
p/4
y(t) dt +
0
(11)
we have
r(t) )
∫
(10) ty(t) dt
t, t ∈ [0, p/4) r(t) ) -t + p/2, t ∈ [p/4, 3p/4) t - p, t ∈ [3p/4, p) r(t ( p) ) r(t)
These can describe the integral q1 )
ty(t) dt
For the triangular wave (Table 1)
)
A)
0 p/2 p/4
(
p/4
p/4
C)
(
)
) )
(14) These quantities are summarized in Table 2. Two-Area (2A) Method. Utilizing the above quantities and ignoring higher harmonic terms, we can obtain the amplitude ratio and phase angle of G(jω). One of the simplest equations is
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010 2
2
π π √(A + B)2 + (A - B)2 q 2 + q22 ) 8hp √ 1 4hp q2 A-B φ ) ∠G(jω) + π ) arctan ) arctan q1 A+B (15) |G(jω)| )
()
(
)
Equation 15 uses only two areas of A and B. Because |A B| , |A + B|, ∠G(jw) ≈ π and |G(jw)| ≈ (π2/(4hp))(A + B). Consequently we have Kcu ≈ 4hp/(π2(A + B)), which is the same as the expression given in the work of Lee et al.15 The performance of eq 15 is very similar to that in the work of Lee et al.15 However, eq 15 provides the phase angle information. Errors in eq 15 are due to the third and higher harmonic terms.
3 1 2p 8hp 2 q5 ) q1 + a ) 2 (|G(jω)| cos(∠G(jω)) + 0 + 4 4π y 25 π |G(j5ω)| cos(∠G(j5ω) + · · · ) 3 3 1 2p 8hp × 0 ) q2 ) 2 (|G(jω)| sin(∠G(jω)) + q6 ) q2 + 4 4π 4 π 2 0 + |G(j5ω)| sin(∠G(j5ω) + · · · )(16) 25 Then we have |G(jω)| ) π2 8hp
Two-Area with Amplitude (2AA) Method. To reduce errors due to the third harmonic term of G(j3ω), we combine eqs 4 and 9. That is, let
π2 q 2 + q62 ) 8hp √ 5
( 43 qq + 41 2pπ a ) + (43q ) 2
1
φ ) arctan
() 6
q5
y
2
( (
(17)
2
) arctan q2 / q1 +
2pay 3π
))
Table 1. Weighting Functions for Integrations of the Output
Table 2. Equations for Some Quantitiesa quantity
a
series expression
y(0) ) 0
4h 1 |G(jω)| sin(∠G(jω)) + |G(j3ω)| sin(∠G(j3ω)) + · · · π 3
y(p/4) ) ay
1 4h |G(jω)| cos(∠G(jω)) - |G(j3ω)| cos(∠G(j3ω)) + · · · π 3
q1 ) ∫p0 q(t)y(t) dt ) 2(A + B)
1 8hp |G(jω)| cos(∠G(jω)) + |G(j3ω)| cos(∠G(j3ω)) + · · · 9 π2
q2 ) ∫p0 q(t -(p/4))y(t) dt ) 2(A - B)
1 8hp |G(jω)| sin(∠G(jω)) - |G(j3ω)| sin(∠G(j3ω)) + · · · 9 π2
q3 ) ∫p0 r(t)y(t) dt ) 2(C - D) + pB
1 4hp2 |G(j3ω)| cos(∠G(j3ω)) + · · · |G(jω)| cos(∠G(jω)) 27 π3
q4 ) ∫p0 r(t -(p/4))y(t) dt ) -2(C + D) + 0.5p(A + B)
1 4hp2 |G(j3ω)| sin(∠G(j3ω)) + · · · |G(jω)| sin(∠G(jω)) + 27 π3
p/2 p/2 A ) ∫p/4 y(t) dt, B ) ∫p/4 y(t) dt, C ) ∫p/4 ty(t)dt, D ) ∫p/4 ty(t)dt. 0 0
7809
(
)
(
(
(
) )
)
(
(
)
)
7810
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
Figure 2. Relative errors in |G(jw)| and the angle φ for processes, G(s) ) exp(-θs)/(s + 1)n.
Figure 4. Estimates (circles) of G(jω) and G(j3ω).
Figure 3. Relative errors in Kcu and Pu for processes, G(s) ) exp(-θs)/ (s + 1).
Figure 2 shows errors in eq 17 for FOPTD processes and critically damped second-order plus time delay (SOPTD) processes. We can see that errors in the amplitude ratios are within 5%. Errors in eq 17 are due to approximation of y(p/4) ) ay and the fifth and higher harmonic terms. Four-Area (4A) Method. Consider the following quantities 6π 32hp q ) |G(jω)| cos(∠G(jω)) + 0 + p 3 π2 2 |G(j5ω)| cos(∠G(j5ω)) + · · · 125 6π 32hp ( q ) q8 ) q2 + |G(jω)| sin(∠G(jω)) + 0 + p 4 π2 2 |G(j5ω)| sin(∠G(j5ω)) + · · · (18) 125 Consequently, q7 ) q1 +
(
)
π2 q 2 + q87 32hp √ 7 q8 φ ) arctan q7 |G(jω)| )
()
(19)
Figure 5. Estimated models for the process, G(s) ) exp(-θs)/(s + 1)(2s + 1)(4s + 1).
will have estimation errors smaller than eq 17 because eq 18 has smaller fifth harmonic terms than eq 16. The weighting function corresponding to eq 18 is (Table 1) qr6(t) ) q(t) +
16 6π 2 r(t) ) sin(5ωt) + · · · sin(ωt) + 0 + p π 25
(
)
(20)
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
7811
This has no third harmonic term and eliminates the third harmonic terms in eq 18. Figure 2 shows errors in eq 19 for FOPTD processes and critically damped SOPTD processes. We can see that errors in the amplitude ratios are within 1%. Errors in eq 19 are due to the fifth and higher harmonic terms. Ultimate Gain and Period. Previously, Kcu ) 1/|G(jω)| was used. In this case, the accurate estimation of |G(jω)| does not guarantee the accuracy of Kcu and there can be a significant error in Kcu because the oscillation period can differ from the ultimate period. Previous methods without changing the shape of relay output have no way to reduce errors in estimating the ultimate gain because the information of ∠G(jω) is not estimated. On the other hand, we have information about ∠G(jω) and the estimation errors in the ultimate gain and period can be reduced with correlations. For this, we use a linear correction as Kcu ) (1 + 0.15φ)/|G(jω)| Pu ) (1 - 0.3φ)p
(21)
These correlations are obtained from results of the low-order processes. Figure 3 shows errors of estimates of the ultimate gain and period for FOPTD processes. For a higher-order process, the phase angle φ is very small and both estimates with and without correlations will yield accurate ultimate gain and period. Figure 6. Flowchart of the 4A method to compute the ultimate gain and the ultimate period of the process.
Estimation of G(j3ω). The frequency response at the frequency of 3ω can be obtained with combining quantities of q1-q4 differently. Consider
Table 3. Simulation Results ultimate data process e-θs (s + 1)(2s + 1)(4s + 1) (-0.5s + 1)e-θs (s + 1)3 a
0.5 (Kcu 3 (Kcu 0.1 (Kcu 3 (Kcu
) 6.26 Pu ) 8.96) ) 2.44 Pu ) 16.54) ) 2.94 Pu ) 5.63) ) 1.357 Pu ) 12.54)
parametric model
eq 1
proposed
Kcu ) 6.05 (-3.5%) Pu ) 9.10 (1.5%) Kcu ) 2.35 (-3.8%) Pu ) 16.4 (-0.9%) Kcu ) 2.79 (-4.9%) Pu ) 5.74 (2.0%) Kcu ) 1.32 (-2.5%) Pu ) 12.2 (-2.6%)
Kcu ) 6.23 (-0.6%) Pu ) 9.00 (0.4%) Kcu ) 2.46 (0.7%) Pu ) 16.5 (-0.2%) Kcu ) 3.06 (4.1%) Pu ) 5.54 (-1.5%) Kcu ) 1.36 (0.3%) Pu ) 12.5 (-0.2%)
θ
IAE (integral of absolute error) for the unit step response.
Table 4. Correction Factors for Relay Feedback Systems
FOPTD15 exp(-2.48s)/(8.70s IAEa ) 3.66 exp(-5.24s)/(5.86s IAE ) 1.10 exp(-1.73s)/(2.45s IAE ) 0.61 exp(-4.79s)/(1.85s IAE ) 0.22
proposed + 1) + 1) + 1) + 1)
exp(-1.18s)/(3.27s IAE ) 0.22 exp(-3.67s)/(3.15s IAE ) 0.06 exp(-1.12s)/(1.25s IAE ) 0.07 exp(-4.12s)/(1.20s IAE ) 0.06
+ 1)2 + 1)2 + 1)2 + 1)2
7812
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010
q9 ) q1 -
q10
2π 32hp q ) 0 + |G(j3ω)| cos(∠G(j3ω)) + p 3 27π2 27 |G(j5ω)| cos(∠G(j5ω)) + · · · 125
channel relay.18 With slight modifications as shown in Table 4, we can apply the proposed methods to them.
(
)
2π 32hp 27 q ) ) q2 0 - |G(j3ω)| sin(∠G(j3ω)) + p 4 125 27π2
(
)
|G(j5ω)| sin(∠G(j5ω)) + · · · (22) Then we have 27π2 2 q + q102 32hp √ 9 q10 ∠G(j3ω) ) -arctan q9 |G(j3ω)| )
( )
(23)
2π 9 16 qr2(t) ) q(t) 0 + sin(3ωt) + r(t) ) sin(5ωt) + · · · p 9π 25
)
(24) This has no first harmonic term and eliminates the first harmonic terms in eq 21, enabling the estimation of G(j3ω). Figure 4 shows estimation performances for G(jω) and G(j3ω). Parametric Model Identification. Parametric models of processes are very useful. Information in the relay oscillation is usually not enough. Relay oscillation with the process steady state gain is considered. The process steady state gain can be obtained easily in the course of relay feedback test14 or after the relay feedback test with a step set point change. We consider a process model G(s) )
ke-θs (τs + 1)n
Methods to obtain accurate frequency responses of the first and third harmonics and the ultimate gain/period of a process are proposed, along with parametric models from oscillations of the conventional relay feedback system. These methods use two or four areas of the relay feedback response that are simple to compute. The first and third harmonic frequency responses can be used to adjust tuning rules based on the ultimate gain and period of process and to determine the order of parametric models. Appendix (Orthogonality of Trigonometric Functions17)
The weighting function corresponding to eq 21 is (Table 1)
(
Conclusions
(25)
For a set of trigonometric functions, {1, cos(ωt), cos(2ωt), ..., sin(ωt), sin(2ωt), ...}, we have
∫ sin(nωt) sin(mωt) dt ) {0,p/2, nn *) mm ∫ sin(nωt) cos(mωt) dt ) 0 0, n*m ∫ cos(nωt) cos(mωt) dt ) p/2, n ) m * 0 p
0 p
{
0
p
0
p,
(A1)
n)m)0
for p ) 2π/ω. Applying these, we can obtain simple equations for integrals such as
∫
p
0
p p f(t)g(t) dt ) Rmξa + βmξb 2 2
(A2)
for functions of f(t) ) R0 + R1 cos(ωt) + β1 sin(ωt) + R2 cos(2ωt) + β2 sin(2ωt) + ... and g(t) ) ξa cos(mωt) + ξb sin(mωt).
The process order n is determined as n)
{
1, |G(j3ω)|/|G(jω)| g 0.3 2, |G(j3ω)|/|G(jω)| < 0.3
Literature Cited
(26)
Then, from G(jω) ) ke-jθω/(jτω + 1)n, we have3 τ)
1 ω
k ) ( |G(jω)|
2/n
-1
(27)
-∠G(jω) - n arctan(τω) (28) ω For FOPTD processes, |G(j3ω)|/|G(jω)| g 1/3, and the threshold value in eq 26 is chosen from this. Panda and Yu11 showed that the critically damped SOPTD model (n ) 2) is effective for wide range of processes. Figure 5 shows Nyquist plots and step responses of parametric models obtained. We can see that the proposed SOPTD models are very accurate for processes simulated. Applications. Figure 6 shows the flowchart to compute the ultimate gain and period by the 4A method. The proposed methods are applied to various processes. Two of them are listed in Table 3. Better results are obtained with the proposed 4A method. Both processes are high order and attenuate higher harmonic terms; the conventional method of eq 1 also provides accurate estimates of ultimate gain and period. On the other hand, in estimating the parametric model, the proposed models are more accurate than models of the conventional method. The proposed area methods can be applied to other types of relay feedback systems whenever they show symmetric responses. These include the relay with hysteresis2 and the twoθ)
(1) Astrom, K. J.; Hagglund, T. Automatic Tuning of Simple Regulators with Specifications on Phase and Amplitude margins. Automatica 1984, 20, 645. (2) Astrom, K. J.; Hagglund, T. PID Controllers; Instrument Society of America: North Carolina, 1995. (3) Yu, C. C. Autotuning of PID Controllers: A Relay Feedback Approach; Springer: London, 2006. (4) Tan, K. K.; Lee, T. H.; Huang, S.; Chua, K. Y.; Ferdous, R. Improved Critical Point Estimation Using a Preload Relay. J. Process Control 2005, 16, 445. (5) Sung, S. W.; Park, J. H.; Lee, I. Modified Relay Feedback Method. Ind. Eng. Chem. Res. 1995, 34, 4133. (6) Hang, C. C.; Lee, T. H.; Ho, W. K. AdaptiVe Control; Instrument Society of America: North Carolina, 1993. (7) Shen, S. H.; Wu, J. S.; Yu, C. C. Autotune Identification under Load Disturbance. Ind. Eng. Chem. Res. 1996, 35, 1642. (8) Sung, S. W.; Lee, J. Relay Feedback Method under Large Static Disturbances. Automatica 2006, 42, 353. (9) Luyben, W. L. Getting More Information from Relay-feedback Tests. Ind. Eng. Chem. Res. 2001, 40, 4391–4402. (10) Kaya, I.; Atherton, D. P. Parameter Estimation from Relay Autotuning with Asymmetric Limit Cycle Data. J. Process Control 2001, 11, 429. (11) Panda, R. C.; Yu, C. C. Analytic Expressions for Relay Feedback Responses. J. Process Control 2003, 13, 48. (12) Sung, S. W.; O, J.; Lee, I.; Lee, J.; Yi, S. H. Automatic Tuning of PID Controller Using Second-Order Plus Time Delay Model. J. Chem. Eng. Jpn. 1996, 29, 990. (13) Sung, S. W.; Lee, I. Enhanced Relay Feedback Method. Ind. Eng. Chem. Res. 1997, 36, 5526. (14) Huang, H. P.; Jeng, J. C.; Luo, K. Y. Auto-tune System Using Single-run Relay Feedback Test and Model-based Controller Design. J. Process Control 2005, 15, 713.
Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010 (15) Lee, J.; Sung, S. W.; Edgar, T. F. Integrals of Relay Feedback Responses for Extracting Process Information. AIChE J. 2007, 53, 2329– 2338. (16) Seborg, D. E.; Edgar, T. F.; Mellichamp, D. A. Process Dynamics and Control, 2nd ed.; Wiley: New York, 2004. (17) Kreyszig, E. AdVanced Engineering Mathematics; Wiley: New York, 1999.
7813
(18) Friman, M.; Waller, K. V. A Two Channel Relay for Autotuning. Ind. Eng. Chem. Res. 1997, 36, 2662.
ReceiVed for reView October 1, 2009 ReVised manuscript receiVed November 2, 2009 Accepted November 3, 2009 IE901546J
