1538
Ind. Eng. Chem. Res. 2007, 46, 1538-1546
Modified Relay Feedback Identification Based on Describing Function Analysis Ping Wang, Danying Gu, and Weidong Zhang* Department of Automation, Shanghai Jiaotong UniVersity, Shanghai 200240, People’s Republic of China
Describing function analysis has been widely adopted in many existing identification methods that are based on relay feedback tests. In these methods, the critical point of the Nyquist curve is usually approximated to the oscillation point. Nevertheless, such an approximation may result in poor accuracy for identification, in particular, for those processes with large dead time. In this paper, we introduce the concept of phase deviation to compensate for the span between the critical point and the oscillation point, so that the ultimate gain and ultimate frequency can be accurately obtained using only a single relay feedback test. The phase deviation is then analytically derived, with which all the model parameters can be identified without any prior information of the dead time or static gain. Numerical examples are given to illustrate the proposed method. 1. Introduction Model identification has an important role in industrial process control, because the performance of the controller is ultimately dependent on how the relationship between the input and output of the process is efficiently represented.1 Based on the estimated models and various performance specifications, most controllers in industrial and chemical practice are tuned using the trial-and-error method. Although the manual control capabilities of control engineers have been proven adequate for a large number of control loops, the manual tuning procedures are time-consuming, in particular, for those plants with slow responses. To overcome this disadvantage, Astrom2 proposed an automatic tuning controller that was based on the relay feedback technique. This method soon became a superior alternative to the conventional tuning. Consequently, many modifications for improved accuracy and performance have been reported in the literature. Luyben3 summarized the applications and extensions of Astrom’s auto-tuning method and proposed the Auto-Tuning Variation (ATV) method. In the ATV method, the ultimate frequency wu and the ultimate gain ku (these represent the most important process information) can be directly extracted using a describing function linear approximation. However, some researchers1,4 have reported that the approximation value may deviate significantly from the true value (for example, an approximation error in ku of 5%-20%). The approximation error comes from the assumption that all the periodic signals are truncated as the uniform sinusoid waveform of their primary harmonic in a Fourier series, or, equivalently, the oscillation point is used as the critical point of the Nyquist curve in the describing function analysis. This observation implies that the observed ultimate period is the oscillation period Tosc, rather than the true ultimate period Tu. Much research has been devoted to overcoming this problem and developing higher-accuracy identification strategies for autotuning controllers. Sung5 used a six-step signal instead of the square signal to better fit the sinusoid curve. Lee6 introduced a mapping function to update the periodic signals iteratively. Noteworthy improvement was obtained on the identification accuracy. Given that the two methods require much additional time on repetitive relay tests, Yu1 suggested using an additional saturated relay to reshape the square signals, to make them more * To whom all correspondence should be addressed. Tel: +86-2134204021,Fax: +86-21-34204019.E-mailaddress:
[email protected].
Figure 1. Block diagram of a relay feedback system.
similar to being sinusoidal. A modified saturated relay method proposed by Sivakumar7 considered more high-order harmonics. These saturated relay methods have attained satisfactory results with relay tests that involve two times. However, they cannot be directly used for various processes, because of a trade-off between accuracy and the existence of the sustained limit cycle. Luyben8 noted that even if wu and ku are obtained with high accuracy, they may be insufficient to identify the model parameters for some processes. To realize auto-tuning, the ATV method applies the prior information of the dead time or the static gain. To avoid the use of prior information, many modified relay methods have been developed. Luyben8 derived an analytical expression of the shape factor to determine model parameters. Friman9 devised a two-channel experiment to obtain more information with two continuous relay tests. Wang10 designed an optimal controller based on the two-channel experiment and achieved satisfactory results. Wang11 and Gu12 derived analytical expressions without any prior information in the time domain, using a single biased relay test. In this paper, the cause of the approximation error is analyzed. A new concept, called phase derivation, is then introduced, to compensate for the loss of those truncated high-order harmonics. The phase deviation is then analytically derived, with which the ultimate gain and the ultimate frequency can be accurately obtained using only a single relay feedback test. In the proposed method, all the model parameters can be identified without any prior information of the dead time or static gain. This paper is organized as follows: In section 2, the phase derivation is defined and analytically derived. In section 3, a new identification algorithm is developed. The proposed identification method is illustrated by several numerical examples in section 4. Finally, some conclusions are given in section 5.
10.1021/ie061141y CCC: $37.00 © 2007 American Chemical Society Published on Web 02/01/2007
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007 1539
Figure 2. Plot of the Nyquist curve. Table 1. Critical Information Estimation for the Equation Gp1(s) ) [1/(s + 1)]e-θs
2. The Phase Deviation It has been noticed that most existing identification methods have restricted attention on processes in the form of a firstorder process with dead time (FOPDT), which, in fact, cannot represent a variety of industrial and chemical processes sufficiently well.13 Without a loss of generality, three typical models are studied here: Model 1, which is a first-order process with dead time (FOPDT):11
G(s) )
(τs +k 1)e
-θs
(
)
[
]
k(-τ1s + 1) (τ2s + 1)s
e-θs
θ
ku
wu
ku
2.0 5.0 10.0
1.52 1.13 1.04
1.14 0.53 0.29
1.45 1.28 1.27
4.2 13.2 22.4
a
Proposed Method
wu
PEa
ku
PEa
wu
PEa
1.16 0.55 0.29
1.7 3.7 2.3
1.52 1.13 1.04
0.0 0.0 0.0
1.14 0.53 0.29
0.0 0.0 0.0
Percentage error.
Table 2. Critical Information Estimation for the Equation Gp2(s) ) [1/(s2 + 1.5s + 1)]e-θs Real Process θ
ku
wu
2.0 5.0 10.0
1.33 1.05 1.01
0.87 0.48 0.27
a
ATV Method ku
PEa
wu
Proposed Method PEa
ku
1.30 2.26 0.89 2.30 1.33 1.20 14.3 0.49 2.08 1.05 1.20 18.8 0.27 0.0 1.04
PEa
wu
PEa
0.0 0.0 0.0
0.87 0.48 0.27
0.0 0.0 0.0
Percentage error.
(2)
where a, b, c, and θ are all positive. Model 3, which is an integrating process with inverse response and dead time:12
G(s) )
ATV Method PEa
(1)
where k, τ, and θ are all positive. Model 2, which is a second-order process with dead time (SOPDT):10
1 e-θs G(s) ) as2 + bs + c
Real Process
(3)
where k, τ1, τ2, and θ are all positive. It is well-known that the oscillation period and the ultimate period are different from each other. Correspondingly, on the Nyquist curve, the oscillation point and the critical point are not superposed. In a review of relay feedback, Astrom2 provided a theorem in discrete time state space to calculate the oscillation period for model 1. Further extensions to high-order models were given by Chang14 and Yu.1 However, the calculation of the oscillation period in these methods is complicated, compared to the calculation of the ultimate period using the describing function linear approximation. Suppose that we can determine the relationship between the oscillation period and the ultimate
period; the oscillation period can then be directly obtained from the ultimate period. In the following discussion, a new concept, called phase deviation, is introduced to measure the gap between the two periods. Figure 1 shows the structure of relay feedback experiments. With the help of the relay, it is possible to identify various models as the sustained limit cycle is established. In the conventional describing function analysis, the negative reciprocal curve (noted by -1/N(a)) of the ideal relay is used as the negative real axis, or, equivalently, the oscillation point is used as the cross point (which is called the critical point) between the negative real axis and the Nyquist curve. If the situation is true, there will be no errors in the estimation of the ultimate information. Unfortunately, -1/N(a) is not exactly the same as the negative real axis. The oscillation point then drifts away from the critical point in the complex plane, as indicated in Figure 2. As a result, the parameters derived from this theory can deviate significantly from the real system parameters. Our idea is, first, to derive the span between the oscillation period and the ultimate period analytically and then use it to compensate for the approximation error.
1540
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007
Table 3. Models for Typical Processes case
real process
ATV method
1)]e-0.2s
1 2 3 4
proposed method
1)]e-0.2s
+ 1.6088s + 1.4279)]e-0.2205s [1/(4.0003s2 + 4.0005s + 1.0012)]e-2.006s [1/(1.831s2 + 1.277s + 1.095)]e-8.723s [1/(2.779s2 + 2.282s + 1.246)]e-6.706s
[1/(2.0743s + [1/(3.5686s + 1)]e-2s [1/(6.6843s + 1)]e-6s [1/(1.665s + 1)]e-8s
+ 1.5s + [1/((2s + 1)2)]e-2s [1/((s + 3)(s2 + s + 1))]e-6s [(-s + 1)/((s + 1)4)]e-4s
[1/(0.5s2
[1/(0.5010s2
Table 4. Parameter Estimation from New Relay Feedback Tests case 1 2
real process
proposed method without noise 1)]e-2s
[(-3s + 3)/(2.8s + [(-3s + 1)/(1.2s + 1)]e-6s
2.9586[(-0.9777s + 1)/(2.7148s + 1)]e-2.0490s 0.9851[(-3.1990s + 1)/(1.3970s + 1)]e-5.6241s
(-3.000[(1.000s + 1))/(2.801s + [(-3.000s + 1)/(1.200s + 1)]e-6.000s
Let us now determine how to introduce a concept to compensate for the approximation error. As shown in Figure 2, the negative reciprocal of the ideal relay (-1/N(a)) is represented by a radial from the origin point to the oscillation point in the second and third quadrants. The angle between the oscillation point and the critical point is defined as the phase deviation (δp1). It will be used to compensate for the approximation error later. The exact definition of the phase deviation δp1 is given in Definition 1. Furthermore, by extending the analysis, phase derivation δpn (for n ) 1, 2, 3, ...) for the nth-order harmonics can also be obtained, which is given in Definition 2. Definition 1. The negative reciprocal curve -1/N(a) intersects G(jw) at the oscillation point in the complex plane. Let wosc be the oscillation frequency. The first-order harmonic phase derivation is δp1 ) -∠G(jwosc) - π + 2Rπ, where R ) 1, 2, 3, ... keeps δp1 ⊂ (-π, π]. Definition 2. The negative reciprocal curve parameter -1/N(a) intersects G(jw) at the oscillation point in the complex plane. Let wosc be the oscillation frequency. The nth-order harmonic phase derivation is δpn ) -∠G(jnwosc) - π + 2Rπ (for n ) 1, 2, 3, ...), where R ) 1, 2, 3, ... keeps δpn ⊂ (-π ,π]. When model 1 is adopted in relay feedback identification, the first-order phase derivation can be derived from the theorem in Astrom,2 which expresses the oscillation period in an analytical way. The detailed calculation procedure is given in the following (Theorem 1). Moreover, in regard to general situations of various models, the nth-order phase derivation (n ) 1, 2, 3, ...), including the first-order phase derivation, can be calculated based on Theorem 2. Theorem 1. Assume that the sustained limit cycle has been established for Model 1. The phase deriVation is giVen as
{
proposed method with noise
1)]e-1.999s
σθ (for σ g 0) + tan-1 σ - π τ δp1 ) σθ + tan-1 σ (for σ < 0) τ
the process. Let Fsn ) (2/Tosc)∫Tosc y(t) sin(nwosct) dt and Fcn ) (2/Tosc)∫Tosc y(t) cos(nwosct) dt. The phase deriVation is
{
( ) ( ) ( ) ( )
-tan-1
δpn )
-tan
-1
Fcn Fcn for g0 Fsn Fsn Fcn Fcn +π for 1
(A4)
According to the definition of δp1, the phase at the oscillation point on the Nyquist curve can be expressed as
∠G(jwosc) ) -π - δp1
(A5)
For Model 1, we have
-θwosc - tan-1(τwosc) ) -π - δp1
(A6)
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007 1545
Substituting wosc in eq A3 into eq A6 yields
δp1 )
[
or, equivalently,
]
πθ/τ π + mπ - π + tan-1 θ/τ θ/τ ln(2e - 1) ln(2e - 1) (for m ) 0, 1) (A7)
yn(t) )
4h |G(jnwosc)|[sin(nwosct) cos(-π - δpn) + nπ cos(nwosct) sin(-π - δpn)]
)
4h |G(jnwosc)|[-sin(nwosct) cos(-δpn) nπ cos(nwosct) sin(-δpn)] (A18)
Some algebraic manipulation gives
{
σθ + tan-1 σ - π (for σ g 0) τ δp1 ) σθ (for σ < 0) + tan-1 σ τ
(A8)
where
σ)
π ln(2eθ/τ - 1)
(A9)
Let
Fsn )
4h |G(jnwosc)| cos(-δpn) nπ
(for n ) 1, 3, 5, ...) (A19)
Fcn )
4h |G(jnwosc)| sin(-δpn) nπ
(for n ) 1, 3, 5, ...) (A20)
Appendix A-2: Proof of Theorem 2 In the relay feedback test, the output of the relay (u(t)) is expanded into a Fourier series:
y(t) can be written as ∞
y(t) )
∞
u(t) )
F2l-1 sin[(2l-1)wosct] ∑ l)1
(l ) 1, 2, 3, ...) (A10)
where
Because
(A11)
Let un(t) denote the nth-order components. The fundamental component then can be written as
u1(t) )
4h sin(wosct) π
{ {
0 (for m* n) cos(mwosct)cos(nwosct) dt ) Tosc osc (for m ) n) 2 (A23)
∫T
∫T
(A12)
When u(t) is used as the input of the process in the identification, the corresponding output of process, y(t), can be expanded into the Fourier series, in the form of
sin(mwosct) cos(nwosct) dt ) 0
osc
π
∞
A*
∑ l)1
∫T
)
2l - 1
osc
A* ) |G[j(2l - 1)wosc]|
(A14)
φ* ) ∠G[j(2l - 1)wosc]
(A15)
∠G(jnwosc) ) -π - δpn
(A17)
osc
Fsn sin2(nwosct) dt
FsnTosc 2
(A25)
∫T
osc
Fsn cos2(nwosct) dt
FsnTosc 2
(A26)
Therefore,
Fsn )
According to Definition 2, we have
y(t) cos(nwosct) dt ) )
The nth-order component of y(t) is
4h |G(jnwosc)| sin[nwosct + ∠G(jnwosc)] nπ (n ) 1, 3, 5, ...) (A16)
∫T
and
∫T
where
yn(t) )
y(t) sin(nwosct) dt )
sin[(2l - 1)wosct + φ*] (for l ) 1, 2, 3, ...) (A13)
(for all n, m) (A24)
we get osc
4h
(A21)
0 (for m* n) T sin(mw t) sin(nw t) dt ) osc osc osc Tosc (for m ) n) 2 (A22)
∫ 4h F2l-1 ) π(2l - 1)
y(t) )
∑ [-Fsn sin(nwosct) - Fcn cos(nwosct)] n)1
Fcn )
2 Tosc
∫T
2 Tosc
∫T
osc
y(t) sin(nwosct) dt
(n ) 1, 3, 5, ...) (A27)
osc
y(t) cos(nwosct) dt
(n ) 1, 3, 5, ...) (A28)
In the relay feedback test, the values of y(t) have been obtained at sampling points. With this information, both eqs
1546
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007
A27 and A28 can be transferred into Riemann integrations. Then, δpn can be obtained by the division of eqs A27 and A28:
δpn ) -tan-1
( )
Fcn + mπ Fsn (for m ) 0, 1 and n ) 1, 3, 5, ...) (A29)
or, equivalently,
{
( ) ( ) ( ) ( )
-tan-1
δpn )
-tan-1
Fcn Fcn for g0 Fsn Fsn Fcn Fcn +π for
