Integrals of Relay Feedback Responses for Extracting Process Information

Jietae Lee and Su Whan Sung*

Department of Chemical Engineering

Kyungpook National University

Daegu, 702-701, KOREA

*Corresponding author. E-mail: suwhansung@knu.ac.kr, Tel: +82-53-950-6838, Fax: +82-53-950-6615

and

Thomas F. Edgar,
Department of Chemical Engineering (C0400),

The University of Texas at Austin,
Austin, TX 78712-1062, USA
E-mail: edgar@che.utexas.edu

Keywords : Relay feedback, Integrals, Identification, Ultimate information, PID, autotuning

Abstract

Relay feedback methods have been widely used in tuning proportional-integral-derivative controllers automatically. Previous approaches usually use specific point data such as the oscillation amplitude and period of the relay response. In this paper, we propose new identification methods which use integrals of the relay response instead of the point data. The proposed methods guarantee better accuracy and advantages in obtaining the ultimate information of process as well as parametric models compared with previous approaches because effects of the high order harmonic terms are suppressed significantly by using the integrals of the relay responses.

Introduction

Since Astrom and Hagglund (1984) introduced the autotuning method, which used the relay feedback test, many variations have been proposed for autotuning of PID controllers (Astrom and Hagglund, 1995; Hang et al., 2002; Yu, 2006). Several methods such as a saturation relay (Yu, 2006), relay with a P control preload (Tan et al., 2006) and a two level relay (Sung et al., 1995) were introduced to obtain more accurate ultimate information of the process by suppressing the effects of the high order harmonic terms. To obtain a Nyquist point other than the critical point, a relay with hysteresis or a dynamic element such as time delay has been used (Astrom and Hagglund, 1995; Kim, 1995; Tan et al., 1996; Chiang and Yu, 1993). Recently, a two channel relay has been proposed to obtain a Nyquist point information corresponding to a given phase angle (Friman and Waller, 1997; Sung et al., 2006). Methods to reject unknown load disturbances and restore symmetric relay oscillations have been available (Hang et al., 1993; Shen et al, 1996b; Sung and Lee, 2006). A biased relay has been used to obtain the process steady state gain as well as the ultimate information (Shen et al., 1996a) from only one relay test. Huang et al. (2005) used the integral of the relay transient to obtain the steady state gain of the process.

Many Nyquist points of the process dynamics can be extracted from only one relay experiment by applying the FFT (fast Fourier transformation) technique to the whole transient responses from the start to the final cyclic steady-state part of the relay responses (Wang et al., 1997a). However, the computations are somewhat complex and the complete transient responses must be stored. For the same purpose, Laplace transformation of a periodic function has been used to obtain many frequency responses from one relay test (Ma and Zhu, 2006).

The shape factor (Luyben, 2001) has been used to extract a three-parameter model from the cyclic steady state part of the relay response. Several authors (Kaya and Atherton, 2001; Panda and Yu, 2003) derived exact expressions relating the parameters of the FOPTD process to the measured data of the relay response. They used the analytic expressions to extract parameters of the FOPTD model. However, the methods based only on the cyclic steady state data cannot provide acceptable robustness for uncertainty such as process/model mismatches and nonlinearity. They may provide poor model parameter estimates such as negative gain when the model structure is different from that of the process (Panda and Yu, 2005). The second order plus time delay (SOPTD) model can be also extracted from the cyclic steady state part of the relay response using analytic equations. However, as in the FOPTD model case, the method is also not robust. A relay experiment with a subsequent P control experiment or another relay feedback test can be used to obtain an SOPTD model robustly (Sung et al., 1996; Sung and Lee, 1997).

In this research, identification methods to extract more accurate frequency response information and parametric models from a single conventional relay feedback test are proposed. We use various integrals of the original relay feedback responses to enhance identification performances without modifying the relay feedback system. As in the step response method (Astrom and Hagglund, 1995), areas (equivalently, integrals) will have merits over point data and they are investigated here. Since the conventional relay feedback method is used, the proposed method shares its practical and theoretical merits. Integrals of the relay responses let the fundamental frequency term to be dominant compared to the high harmonic terms, resulting in better accuracy in estimating frequency response information and model parameters. Because it is not required to store the whole trajectories and computations are simple, the proposed methods can be incorporated easily in commercial PID controllers.

Conventional Relay Feedback Method

We consider a classical relay feedback system as shown in Figure 1 and Figure 2 to derive the required equations in this research. We start the relay feedback system at a steady state condition. The relay is first kept on until the process output rises up to a given level and then is set to the normal mode of switching at the instant that the process output crosses a given set point. This relay feedback system will produce a stable oscillation as shown in Figure 2. It is notable that we should set the given level in the beginning of the relay feedback to a significantly large value if we want to extract the zero-frequency information of the process. Otherwise, we cannot guarantee acceptable robustness in extracting zero frequency information from the relay responses.

Astrom and Hagglund, (1984) used this oscillation to extract approximate ultimate information and tune the proportional-integral-derivative (PID) controllers automatically. Let the input and output trajectories be  and  for the conventional relay feedback system, respectively. At the time ,  and  are assumed to be fully developed (cyclic steady state). Then it can be represented by the Fourier series as

                                         (1)

where  and  is the relay amplitude.  Let  and  are the period and the frequency of the relay feedback oscillation. Equation (1) is valid for

                                                                                                       (2)

The output corresponding to  is

      (3)

where  is the process transfer function. Neglecting the high harmonic terms and assuming , we obtain the ultimate frequency w_u, and , . We obtain the following approximate ultimate period  and ultimate gain  as

                                                                                                                 (4)

                                                                                                              (5)

where  is the measured amplitude of . It should be noted that the estimated ultimate information are approximate because we neglect the high harmonic terms. As a result, the ultimate period of equation (4) and the ultimate gain of equation (5) show relative errors up to 5% and 18%, respectively, for the first order plus time delay (FOPTD) process. In this case, the ultimate gain error may not be acceptable.

Proposed Methods for Ultimate Data Estimation

In this research, we use integrals of the process input and output instead of the point data to obtain more accurate process frequency information and parametric process models by suppressing the effects of high harmonic terms.

Proposed equation 1

Let  and  be integrals of the relay responses as

                                                                                                   (6)

                                                                                                   (7)

Then, from equation (1), the response  after  is

                                (8)

where  is the mean value of

                                                                                           (9)

From equation (3), the response  after  can be represented as

   (10)

where  is the mean value of .

                                                                                          (11)

Figure 2 shows typical plots of these responses.

By neglecting the high harmonic terms and assuming  (equivalently, the relay period is the ultimate period), we obtain , . Then, we have the following approximate ultimate gain .

                                                                                                           (12)

where  is the amplitude of . The quantity  is physically the half of the shaded area of  in Figure 2. The proposed method of equation (12) will be superior to equation (5) because the ratios of the high harmonic terms to the fundamental frequency term in  are much smaller than those in  as shown in equations (1),(3),(8) and (10).

Proposed equation 2

From the relay feedback responses in Figure 2, we can construct the following responses as

           (13)

                                                               (14)

We should remark that  is a rectangular wave,  is a triangular wave and is their combination such that the third harmonic term vanishes. The forcing functions of  and  are closer to a sinusoidal wave than . Hence,  and  are closer to the sinusoidal wave than . Figure 3 shows typical plots of these responses.

Approximating the maximum of  as , we have

                                                             (15)

Proposed equation 3

Consider the following quantity.

                                    (16)

The right hand side in equation (16) is derived by applying orthogonality of sine and cosine functions (Appendix A) to equations (3). It is remarked that, to compute the above quantity, the trajectory of  does not need to be stored. By ignoring high harmonic terms in , we obtain

                                                                                                        (17)

When  is sinusoidal,  and equations (5) and (17) provide the same results.

Proposed equation 4

Consider the following quantity.

                       (18)

As for equation (16), the right hand side in equation (18) is derived by applying orthogonality of sine and cosine functions to equation (10). By ignoring high harmonic terms in , we have

                                                                                                      (19)

When  is sinusoidal,  and equations (12) and (19) provide the same results.

Figure 4 shows relative errors in estimated ultimate period and ultimate gains. For FOPTD processes with ratios of time delays to time constants between 0.1 and 5, equations (12), (15), (17) and (19) have relative errors below about 6%. We can see that errors in equation (5) for the ultimate gain can be improved considerably.

Proposed Methods for Nyquist Point Data Estimation

In the previous section, we estimate the ultimate gain on the assumption that the period of relay feedback oscillation is the ultimate period. We omit that assumption in this section

Our goal is to find the following amplitude ratio and the phase lag of the process at the frequency of relay oscillation.

                                                                  (20)

where,  and . From equations (3), (10), (16) and (18), we can obtain the approximate amplitude ratio at the frequency  by neglecting high harmonic terms. Among them, we consider

                                                                                                     (21)

which is the most accurate equation for an amplitude ratio at the frequency . When , its approximation error is

                         (22)

Now, consider the following quantity

      (23)

Ignoring high harmonic terms, we have

                                                                                              (24)

Then, the approximate phase angle of process is .

Figure 5 shows relative errors of the amplitude ratio estimates of equation (21) and phase lag estimates of equation (24) for the FOPTD process. We can see that the estimates of equations (21) and (24) have relative errors below 0.5%.

Steady State Data Estimation

Process dynamic information contained in the cyclic steady state part of the relay responses corresponds to the frequency of  and its multiples. Therefore, to estimate the process information at the frequency zero, we need the transient relay response from the start to the cyclic steady state. Here a method to extract the process steady state information without storing the relay transient and complex computations is introduced. In Ma and Zhu (2006), it is shown that

                              (25)

Let

                                                                                 (26)

Since , we have (Appendix B)

                                                                    (27)

where

The above equation (27) for the process steady state gain is equivalent to that in Huang et al. (2005).

As in the moment analysis of Astrom and Hagglund (1995), we can obtain  (Appendix B)

  (28)

where

This process information is very useful in obtaining higher order process models.

Proposed Methods for First Order Plus Time Delay Model Estimation

From the relay feedback responses, we can obtain the following first order plus time delay (FOPTD) model.

                                                                                                          (29)

Since it has three unknowns, three experimental quantities are needed. From the steady state gain  of equation (27), the peak value of  and the relay oscillation period , the FOPTD model can be determined analytically by the following previous approach (Wang et al., 1997; Kaya and Atherton, 2001; Panda and Yu, 2005).

                                                                                                               (30)

                                                                                                 (31)

                                                                                     (32)

These equations are exact for the FOPTD process.

Proposed equation 5

Equation (31) (consequently, equation (32)) does not provide acceptable accuracy for the case of a large time delay with measurement errors. We overcome this problem by using the integrals of the relay feedback responses. From the oscillation amplitude  of  instead of the peak value of , we can obtain estimate of the time constant as

                                                                                          (33)

This equation can be obtained by rearranging the analytic equations for  and  in Table 1. Here, the proposed last five equations in Table 1 are derived from the first three equations (Panda and Yu, 2005; Kaya and Atherton, 2001). Equation (33) is exact for the FOPTD process, but requires solving a nonlinear algebraic equation for . For a simpler application, we use the following approximate solution:

                                                                       (34)

This equation can be obtained by applying perturbation analysis and numerical technique to the exact nonlinear equation.

Proposed equation 6

From the analytic equation for in Table 1, we can obtain

                                                                                             (35)

Instead of solving the nonlinear equation for , we use an approximate solution:

                                                                                     (36)

This approximation is obtained by utilizing (Abramowitz and Stegun, 1972)

                                                                                           (37)

Proposed equation 7

From the analytic equation for  in Table 1, we can obtain

                                                                                  (38)

Instead of solving the nonlinear equation for , we use an approximate solution:

                                                                              (39)

This approximation is obtained by applying numerical technique together with the approximation of equation (37).

Proposed equation 8

From equation (21) and , we can also obtain

                                                                                          (40)

In summary, we estimate the time constant of the FOPTD model using one of the proposed methods of equations (34), (36), (39) and (40). The static gain and time delay can be estimated by equations (30) and (32), respectively. Figure 6 shows relative approximation errors of equations in estimating the time constant . It is seen that errors are all within 1% except for approximation (40). Figure 7 shows relative error changes when the measurements of , ,  and  are not accurate. We can see that the estimated time constant based on  (previous approach of equation (31)) is very sensitive to the measurement error when the time delay is large. This sensitivity is because the time constant change does not affect the amplitude of  much when the time delay is large. One should be cautious in using the equation (31) when the time delay is expected to be large. Figure 7 shows that the proposed methods can relieve this disadvantage.

The curvature factor by Luyben (2001) can be used to check how large  is. Without measuring an additional data for the Luyben’s curvature factor, as an alternative, a quantity

                                                                                                           (41)

is used here. Figure 8 shows the curvature factor by Luyben (2001) and equation (41) for FOPTD processes. The curvature factor is useful to select an appropriate modeling method before we estimate the model parameters. For example, the method of equation (31) is not appropriate if the proposed curvature factor is larger than about 0.7 (which corresponds to =1.0), as shown in Figure 6.

A FOPTD model can also be obtained without process steady state gain information (Luyben, 2001; Panda and Yu, 2005). However, because information in the relay oscillation responses is concentrated around ultimate frequencies, model parameters can be very inaccurate. So such methods are not recommended in general.

Figure 9 shows the integral of absolute error (IAE) in the frequency domain for high order processes, critically-damped second order plus time delay (SOPTD) processes, under-damped SOPTD processes and processes with inverse responses.

                                                                         (42)

It is seen that the proposed estimations of equations (34), (36) and (39) improve the application of the previous approach of equation (31).

Proposed Method for Critically Damped Plus Time Delay Model Estimation

Panda and Yu (2005) have shown that a three-parameter model of

                                                                                                          (43)

is applicable to a wide range of processes. Model parameters for this critically damped plus time delay (CDPTD) model are estimated. For the steady state gain , equation (30) is used. From equation (21) and , we have

                                                                                             (44)

An analytic equation for the oscillation period (Panda and Yu, 2005; Yu, 2006) is

                                                                                     (45)

where

Equation (45) can be solved for d/t by the Newton-Raphson method with an initial estimate

                                                                                 (46)

Proposed Method for Second Order Plus Time Delay Model Estimation

A general second order plus time delay (SOPTD) model,

                                                                                          (47)

has four unknowns, thus four experimental quantities are needed. We use the process steady state gain information of equation (27), its derivative of equation (28),  of equation (18) and the oscillation period p. We solve

                                                                                                   (48)

                                                                                              (49)

                                             (50)

                                                   (51)

The computational procedure is as follows:

Step 0: From relay feedback response, obtain , ,  and the oscillation period . Let .

Step 1: Calculate  from equation (49).

Step 2: Calculate  from equation (50).

Step 3: Adjust  and repeat Steps 1 and 2 so that equation (51) is satisfied.

Figure 10 shows IAE values in the frequency domain for various processes. For most processes tested, the CDPTD model by equations (43)-(46) shows better results than the FOPTD model by equations (30), (32) and (39). We can see that the SOPTD model by equations (47)-(51) is the best. However, the SOPTD model needs more measurements and iterative computations.

Simulations

Simulations are performed to investigate the performances of methods under sampling, discretization errors and noisy environments. The following process is considered.

                                                                                     (52)

Uniformly distributed noise with mean value and magnitude of 0 and 0.2 is introduced in the output. For noisy responses,  should be adjusted (Astrom and Hagglund, 1995) as

                                                  (53)

The noise characteristics, , can be estimated before starting the relay feedback test. Similarly,  needs to be adjusted as in  when noise is involved.

Simulation results are shown in Figure 11. Models obtained are in Table 2 and their Nyquist plots are in Figure 12. For the FOPTD model, estimation equations of (30), (32) with (39) are used. Because the process (52) is an overdamped process with a moderate time delay, all models estimated have excellent agreement with the true model, as shown by the Nyquist plot of Figure 12.

Conclusions

Without modifying the original relay feedback system, simple and accurate estimates of process ultimate information are easily obtained. We use various integrals of the relay responses instead of point data to reduce the effects of the high harmonic terms. For FOPTD processes, errors over 15% of the previous approaches in estimating the ultimate gain can be reduced to below 5%. We derived the equations to extract a FOPTD, CDPTD and SOPTD models from the process information at the steady state and near the ultimate frequency. They are very simple and can be applied easily to commercial PID controllers.

Appendix A (Orthogonality of trigonometric functions (Kreyzig, 1999))

For a set of trigonometric functions, {1, , , … ,, , …}, we have

                                                           (A1)

for . Applying these relationships, we can derive equations (16) and (17). For example, for a function of ,

                                                 (A2)

This is known as the Parseval theorem.

Appendix B (Moment analyses (Astrom and Hagglund, 1995))

Since

                              (B1)

we have

   (B2)

Then applying integration by parts:

                   (B3)

where  and  and  are periodic after .

By applying the same procedure to

                             (B4)

we can obtain equations for  and :

                                                              (B5)

Literature Cited

Abramowitz, M., and I. A. Stegun, Handbook of mathematical functions, Dover Publications, NY (1972).

Astrom, K. J., and T. Hagglund, “Automatic tuning of simple regulators with specifications on phase and amplitude margins,” Automatica, 20, 645 (1984).

Astrom, K. J. and T. Hagglund, PID controllers., Instrument Society of America, NC (1995).

Bi, Q., Q. G. Wang, and C. C. Hang, “Relay-based Estimation of multiple points on process frequency response,” Automatica, 33, 1753 (1997).

Chiang, R. and C. Yu, “Monitoring Procedure for Intelligent Control: On-Line Identification of Maximum Closed-Loop Log Modulus,” Ind. Eng. Chem. Res., 32, 90 (1993).

Friman, M., and K. V. Waller, “A two channel relay for autotuning,” Ind. Eng. Chem. Res., 36, 2662 (1997).

Hang, C. C., K. J. Astrom, and Q. G. Wang, “Relay feedback auto-tuning of process controllers – a tutorial review,” J. Process Control, 12, 143 (2002).

Hang, C. C., K. J. Astrom, and W. K. Ho, “Relay auto-tuning in the presence of static load disturbance,” Automatica, 29, 563 (1993).

Huang, H. P., J. C. Jeng and K. Y. Luo, “Auto-tune system using single-run relay feedback test and model-based controller design,” J. Process Control, 15, 713 (2005).

Kaya, I., and D. P. Atherton, “Parameter estimation from relay autotuning with asymmetric limit cycle data,” J. Process Control, 11, 429 (2001).

Kim, Y. H., “PI Controller Tuning Using Modified Relay Feedback Method,” J. Chem. Eng. Japan, 28, 118 (1995).

Kreyszig, E., Advanced engineering mathematics, Wiley, New York (1999).

Luyben, W. L., “Getting more information from relay feedback tests,” Ind. Eng. Chem. Res., 40, 4391 (2001).

Ma, M. D. and X. J. Zhu, “A simple auto-tuner in frequency domain,” Comp. Chem. Eng., 30, 581 (2006).

Panda, R. C. and C. C. Yu, “Analytic expressions for relay feedback responses,” J. Process Control, 13, 48 (2003).

Panda, R. C. and C. C. Yu, “Shape factor of relay response curves and its use in autotuning,” J. Process Control, 15, 893 (2005).

Shen, S. H., J. S. Wu, and C. C. Yu, “Use of biased relay feedback method for system identification,” AIChE J., 42, 1174 (1996a).

Shen, S. H., J. S. Wu, and C. C. Yu, “Autotune identification under load disturbance,” Ind. Eng. Chem. Res., 35, 1642 (1996b).

Sung, S. W., J. H. Park and I. Lee, “Modified relay feedback method,” Ind. Eng. Chem. Res., 34, 4133 (1995).

Sung, S. W., and I. Lee, “Enhanced relay feedback method,” Ind. Eng. Chem. Res., 36, 5526 (1997).

Sung, S. W., and J. Lee, “A Two-Channel Relay Feedback Method under Static Disturbances,” Ind. Eng. Chem. Res., 45, 4071 (2006).

Sung, S. W., and J. Lee, “Relay feedback method under large static disturbances,” Automatica, 42, 353 (2006).

Sung, S. W., J. O, I. Lee, J. Lee and S. H. Yi, “Automatic Tuning of PID Controller using Second-Order Plus Time Delay Model,” J. Chem. Eng. Japan, 29, 990 (1996).

Tan, K. K., T. H. Lee, S. Huang, K. Y. Chua and R. Ferdous, “Improved critical point estimation using a preload relay,” J. Process Control, 16, 445 (2005).

Tan, K. K., T. H. Lee and Q. G. Wang, “Enhanced automatic tuning procedure for process control of PI/PID controllers,” AIChE J., 42, 2555 (1996).

Vivek, S., and M. Chidambaram, “Identification using single symmetrical relay feedback test,” Comp. Chem. Eng., 29, 1625 (2005).

Wang, Q. G.., C. C. Hang and Q. Bi, “Process frequency response estimation from relay feedback,” Control Eng. Practice, 5, 1293 (1997).

Yu, CC., Autotuning of PID controllers: A relay feedback approach, Springer, London (2006).
Captions

Table 1. Analytic equations for relay feedback responses.

Table 2. Identification results of the proposed method.

Figure 1. A conventional relay feedback system with integrals of responses.

Figure 2. Typical relay feedback responses and their integrals.

Figure 3. Normalized wave forms of fully-developed relay feedback responses.

Figure 4. Relative error in estimation of the ultimate period and gain for FOPTD processes (G(s)=exp(-ds)/(ts+1)).

Figure 5. Relative error in estimation of the amplitude ratio and phase angle for FOPTD processes (G(s)=exp(-ds)/(ts+1)).

Figure 6. Relative error in estimations of the time constant t for FOPTD processes (G(s)=exp(-ds)/(ts+1)).

Figure 7. Sensitivities of the estimated time constant t for FOPTD processes (G(s)=exp(-ds)/(ts+1)).

Figure 8. Curvature factor plot for FOPTD processes (G(s)=exp(-ds)/(ts+1)).

Figure 9. IAE plots of FOPTD models for several higher order processes (Solid line: Eq. (31), dotted line: Eq. (34), dashed line: Eq. (36), dash-dotted line: Eq. (39)).

Figure 10. IAE plots for several higher order processes (Solid line: FOPTD model (Eqs. (30), (32) and (39)), dotted line: CDPTD model (Eqs. (43)-(46)), dashed line: SOPTD model  (Eqs. (47)-(51))).

Figure 11. Relay feedback responses subject to measurement noise.

Figure 12. Nyquist plots of the process of equation (52) and its estimates (Solid line: FOPTD model (Eqs. (30), (32) and (39)), dotted line: CDPTD model (Eqs. (43)-(46)), dashed line: SOPTD model  (Eqs. (47)-(51))).
Table 1. Analytic equations for relay feedback responses.

Process	Analytic equations

Table 2. Identification results of the proposed method.

Process	Run No.	Output Noise Size	(1.00)*	(-3.70)*	Identified Models	IAE
	1	0	1.009	-3.780		0.093
						0.019
						0.013
	2	0.1	1.025	-4.114		0.062
						0.039
						0.085
	3	0.1	1.002	-3.432		0.076
						0.040
						0.053

* Exact value

Figure 1. A conventional relay feedback system with integrals of responses.

Figure 2. Typical relay feedback responses and their integrals.

Figure 3. Normalized wave forms of fully-developed relay feedback responses.

Relative Errors

Figure 4. Relative error in estimation of the ultimate period and gain for FOPTD processes (G(s)=exp(-ds)/(ts+1)).

Relative Errors

Figure 5. Relative error in estimation of the amplitude ratio and phase angle for FOPTD processes (G(s)=exp(-ds)/(ts+1)).

Relative Errors

Figure 6. Relative error in estimations of the time constant t for FOPTD processes (G(s)=exp(-ds)/(ts+1)).

Relative Errors

Figure 7. Sensitivities of the estimated time constant t for FOPTD processes (G(s)=exp(-ds)/(ts+1)).

Curvature Factor

Figure 8. Curvature factor plot for FOPTD processes (G(s)=exp(-ds)/(ts+1)).

IAE

Figure 9. IAE plots of FOPTD models for several higher order processes (Solid line: Eq. (31), dotted line: Eq. (34), dashed line: Eq. (36), dash-dotted line: Eq. (39)).

Figure 10. IAE plots for several higher order processes (Solid line: FOPTD model (Eqs. (30), (32) and (39)), dotted line: CDPTD model (Eqs. (43)-(46)), dashed line: SOPTD model  (Eqs. (47)-(51))).

Figure 11. Relay feedback responses subject to measurement noise.

Figure 12. Nyquist plots of the process of equation (52) and its estimates (Solid line: FOPTD model (Eqs. (30), (32) and (39)), dotted line: CDPTD model (Eqs. (43)-(46)), dashed line: SOPTD model  (Eqs. (47)-(51))).