Identification of integrating and unstable processes from relay feedback

Computers and Chemical Engineering 32 (2008) 3038–3056

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

Identification of integrating and unstable processes from relay feedback Tao Liu, Furong Gao ∗ Department of Chemical Engineering, Hong Kong University of Science & Technology, Clear Water Bay, Kowloon, Hong Kong

a r t i c l e

i n f o

Article history: Received 14 February 2007 Received in revised form 18 April 2008 Accepted 22 April 2008 Available online 29 April 2008 Keywords: Integrating and unstable processes Identification Relay feedback First-order-plus-dead-time (FOPDT) Second-order-plus-dead-time (SOPDT) Limit cycle

a b s t r a c t In this paper, the exact relay response expressions are derived for integrating and unstable 1st or 2nd order processes. Identification algorithms are subsequently developed in terms of the measurable parameters under relay feedback and fitting conditions of the limit cycle. It is demonstrated that the limit cycle can be undoubtedly formed for an integrating process under a biased/unbiased relay feedback test. A tighter limiting condition is given for identification of unstable processes with the standard relay feedback structure. Denoising methods are presented to cope with measurement noise. Illustrative examples are given to demonstrate the effectiveness and merits of the proposed identification algorithms. © 2008 Elsevier Ltd. All rights reserved.

1. Introduction Many chemical processes, such as heating boilers and continuous-stirred-tank-reactors (CSTRs), have integrating or unstable dynamics. Model-based control strategies have been developed for such processes with effective set-point tracking and disturbance rejection (e.g., Liu, Cai, Gu, & Zhang, 2005a; Liu, Zhang, & Gu, 2005b; Marchetti, Scali, & Lewin, 2001). For the purpose of modeling, relay feedback identification has attracted increasing attention because it is such a closed-loop exercise that the process response does not drift far away from the set-point during the test. Based on a standard and a modified relay feedback scheme, Kaya and Atherton (2001) and Kaya (2006) proposed identification methods for integrating and unstable processes in terms of the so-called A-locus analysis. Compared to Luyben (2003) using the step response data, Gu, Ou, Wang, and Zhang (2006) reported an improved relay identification method using the time domain response analysis for a class of integrating processes with inverse response. Majhi and Atherton (2000) and Majhi (2007) developed relay identification algorithms for integrating and unstable processes from a state-space approach. For integrating processes, Ho, Feng, and Gan (1996) proposed a modified relay identification structure by inserting a differentiator in front of the relay to enhance identification accuracy, and by comparison, Tsang, Lo, and Rad (2000) suggested the use of a derivative filter in the relay feedback path. For unstable processes, Tan, Wang, and Lee (1998) proposed a FOPDT modeling method according to the gain and phase conditions for steady oscillation; Shiu, Hwang, and Li (1998) used two relay response points with the maximal and minimal gains to estimate the first-order-plus-dead-time (FOPDT) and second-order-plus-dead-time (SOPDT) models; Vivek and Chidambaram (2005) and Ramakrishnan and Chidambaram (2003) recently reported two improved FOPDT and SOPDT identification algorithms using the fitting conditions established from the Fourier or Laplace transform for the relay response. Few papers, however, reported the limiting conditions to form steady oscillation from relay feedback for integrating and unstable processes, with the exceptional works of Tan et al. (1998), Majhi and Atherton (2000), Thyagarajan and Yu (2003), and Lin, Wang, and Lee (2004) that were specifically for unstable processes. To ensure identification stability, closed-loop identification methods with a proportional (P), proportional-integral (PI) or proportionalintegral-derivative (PID) controller have been simultaneously developed (to enumerate a few, Kwak, Sung, & Lee, 1997; Lee & Sung, 1993; Paraskevopoulos, Pasgianos, & Arvanitis, 2004; Sree & Chidambaram, 2006; Sung & Lee, 2006). The key idea behind these identification methods is the use of a P/PI/PID controller to stabilize an integrating or unstable process first and then to perform the set-point step tests for ‘open-loop’ identification. A modified relay feedback structure with a proportional-derivative (PD) controller to stabilize an unsta-

∗ Corresponding author. Tel.: +852 2358 7139; fax: +852 2358 0054. E-mail addresses: [email protected] (T. Liu), [email protected] (F. Gao). 0098-1354/$ – see front matter © 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2008.04.006

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

3039

Fig. 1. Block diagram of relay feedback identification.

ble process was recently presented by Padhy and Majhi (2006) for identification of a FOPDT model based on the describing function analysis. In this paper, a systematic identification method using a single relay feedback test is proposed to obtain a FOPDT or SOPDT model for an integrating or unstable process. The corresponding relay response expressions are derived analytically for evaluation of the relay response and choice of the model structure. It is proved herein that the limit cycle can be undoubtedly formed for an integrating process under a biased/unbiased relay test. Besides, a tighter limiting condition to form steady oscillation is also derived for unstable processes. The proposed identification algorithms are developed using the measurable parameters in steady oscillation and fitting conditions of the limit cycle. To cope with measurement noise in practice, a statistical averaging method for measurement is given to ensure identification robustness in the presence of low noise level, and a low-pass Butterworth filter is suggested to deal with high noise level. The paper is organized as follows: Section 2 briefly introduces the implementation of relay feedback identification, and then, the exact relay response expressions and corresponding identification algorithms for obtaining FOPDT and SOPDT models of integrating and unstable processes are detailed, respectively, in four subsections. Section 3 presents two denoising methods to cope with measurement noise. Several examples are given in Section 4 to demonstrate the effectiveness of the proposed identification algorithms. Finally, conclusions are drawn in Section 5. 2. Relay feedback identification A general relay feedback structure for identification is shown in Fig. 1, where the relay function is usually specified as



u(t) =

u+ u−

for{e(t) > ε+ } or {e(t) ≥ ε− and u(t− ) = u+ } for{e(t) < ε− } or {e(t) ≤ ε+ and u(t− ) = u− }

(1)

where u+ =  + 0 and u− =  − 0 denote, respectively, the positive and negative relay magnitudes; ε+ and ε− denote, respectively, the positive and negative relay switch hystereses. Note that u+ = −u− and ε+ = −ε− leads to an unbiased (or named symmetrical) relay function. To avoid measurement noises causing incorrect relay switches, a short ‘listening period’ (e.g., 10–50 samples) may be adopted for reference to set ε+ and ε− properly. It is usually suggested to set ε+ and ε− at least twice over the noise band, together with an upper limit almost equal to 0.95 times of the absolute minimum of u+ and u− (Wang, Lee, & Lin, 2003). Without loss of generality, a biased relay feedback test is used here to derive the exact relay response expressions for integrating and unstable processes, according to the most commonly used FOPDT and SOPDT models which are respectively in the form of GI−1 =

kp e−s s

(2)

GI−2 =

kp e−s s(p s + 1)

(3)

GU−1 =

kp e−s p s − 1

(4)

GU−2 =

kp e−s (1 s − 1)(2 s + 1)

(5)

where kp denotes the process gain,  the process time delay and  p ( 1 and  2 ) the process time constant(s). All of these parameters need to be identified for desirable representation of the process response characteristics, in particular for the low frequency range with a phase change from 0 to −, which covers the most important response characteristics for control system design (Wang et al., 2003; Yu, 2006). Note that the relay response expressions to be derived can be easily generalized to unbiased relay test, for which the limit cycle is symmetrical across the time axis and P+ = P− = Pu /2, where P+ and P− are respectively the positive and negative half periods of the limit cycle. For clarity, the corresponding identification algorithms, along with these relay response expressions, are presented in the following four subsections, respectively. It should be noted that the proposed identification algorithms for obtaining a FOPDT or SOPDT integrating process model may be implemented as well under an unbiased relay feedback test. Besides, if the process gain is known a priori, the proposed identification algorithms for obtaining a FOPDT or SOPDT unstable process model can also be implemented for unbiased relay test. 2.1. FOPDT integrating process model For a first-order integrating process modeled by Eq. (2), the following proposition gives the exact relay response expression for steady oscillation:

3040

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

Fig. 2. Relay response depiction for a FOPDT integrating process.

Proposition 1. Under a biased relay feedback test as shown in Fig. 2, the output response of a first-order integrating process converges to a limit cycle characterized by y+ (t) = kp ( + 0 )t + kp ( − 0 )t0 ,

t ∈ [0, P+ ]

y− (t) = kp ( − 0 )(t + t0 ) + kp ( + 0 )P+ ,

t ∈ [0, P− ]

(1.1) (1.2)

where y+ (t) denotes the monotonically ascending part for t ∈ [0, P− ], corresponding to t ∈ [P+ ,Pu ] in the limit cycle while y− (t) indicates the monotonically descending part for t ∈ [P+ ,Pu ], and Pu = P+ + P− is the oscillation period. Proof.

See Appendix A.



It follows from Proposition 1 that the process time delay can be directly measured as the time to reach the peak of the process output response from the initial relay switch point in a half period of the relay, which is denoted as tp∗ in Fig. 2. Correspondingly, it can be obtained from (1.1) and (1.2) that y+ (0) = kp ( − 0 )t0 = A−

(1.3)

y− (0) = kp ( − 0 )t0 + kp ( + 0 )P+ = A+

(1.4)

Substituting (1.3) into (1.4) yields kp =

A+ − A− ( + 0 )P+

(1.5)

It should be noted that the process gain cannot be derived from a biased relay test as

 pu

kp = G(0) =

0pu 0

y(t) dt

(1.6)

u(t) dt

which is commonly used for identification of a stable process (Wang et al., 2003). The reason lies in G(0) → ∞ for an integrating process, as can be seen from Eq. (2) or Eq. (3). Therefore, relay feedback identification of a FOPDT integrating process model can be summarized in the following algorithm named Algorithm-I-1a: (i) Measure P+ , P− , A+ and A− from the limit cycle; (ii) Measure the process time delay as tp∗ , which is the time to reach the positive peak (A+ ) of the process output response from the initial relay switch point in a negative half period of the relay (P− ); (iii) Compute the process gain, kp , from (1.5). Note that the process response at the oscillation frequency may be formulated as

 Pu

G(jωu ) =

0Pu 0

y(t) e−jωu t dt u(t) e−jωu t

dt

= Au ejϕu ,

ωu =

2 Pu

(1.7)

Substituting the FOPDT integrating process model into (1.7), we obtain kp = Au ωu −ωu −

 = ϕu , 2

(1.8) ϕu ∈



−,

− 2

 (1.9)

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

3041

Fig. 3. Relay response depiction for a SOPDT integrating process.

It should be mentioned that ϕu → − at the oscillation frequency according to the describing function analysis (Yu, 2006), but it is not equal to − due to the phase lag caused by the relay function. Consequently, the process gain and time delay can be derived from (1.8) and (1.9) as kp = Au ωu =−

(1.10)

2ϕu +  2ωu

(1.11)

Therefore, an alternative identification of a FOPDT integrating process model can be summarized in the following algorithm named Algorithm-I-1b: (i) (ii) (iii) (iv)

Measure Pu from the limit cycle; Compute G(jωu ) from (1.7); Compute the process gain, kp , from (1.10); Compute the process time delay, , from (1.11).

Remark 1. Owing to the process response at the oscillation frequency is exactly represented by the FOPDT model derived from AlgorithmI-1b, it may be utilized to obtain better identification accuracy for a high order integrating process compared with Algorithm-I-1a, but at the cost of more computation effort due to the numerical integral involved. It should be noted that the periodic integral for computing G(jωu ) is independent of the initial integral moment owing to the integral property of a periodic integrand. This integral can be practically computed from

 tos +Pu

G(jωu ) =

t

 tosos +Pu tos

y(t) e−jωu t dt u(t) e−jωu t dt

where tos may be taken as any relay switch point in steady oscillation, and the trapezoidal rule can be utilized to implement the numerical integral as done in the examples given in Section 4. 2.2. SOPDT integrating process model For a second-order integrating process modeled by Eq. (3), the following proposition gives the exact relay response expression for steady oscillation: Proposition 2. Under a biased relay feedback test as shown in Fig. 3, the output response of a second-order integrating process converges to a limit cycle characterized by y+ (t) = kp ( + 0 )(t − p ) + kp ( − 0 )t0 + 2kp 0 p E e−t/p , y− (t) = kp ( − 0 )(t + t0 − p ) + kp ( + 0 )P+ + 2kp 0 p F e

t ∈ [0, P+ ]

−t/p

,

t ∈ (0, P− ]

(2.1) (2.2)

where E=

1 − e−p− /p 1 − e−pu /p

F =− Proof.

(2.3)

1 − e−p+ /p 1 − e−pu /p

See Appendix B.

(2.4) 

3042

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

It can be derived from (2.1) and (2.2) that dy+ (t) = kp ( + 0 ) − 2kp 0 E e−t/p dt

(2.5)

dy− (t) = kp ( − 0 ) − 2kp 0 F e−t/p dt

(2.6)

By letting dy+ (t)/dt = 0 and dy− (t)/dt = 0, the time to reach the single extreme value of y+ (t) and y− (t) can be derived, respectively, as tp+ = p ln

20 E  + 0

(2.7)

tp− = p ln

20 F  − 0

(2.8)

It can be easily verified that d2 y+ (t)/dt2 > 0 and d2 y− (t)/dt2 < 0. Therefore, tp+ is the time to reach the minimum of y+ (t), while tp− is the time to reach the maximum of y− (t). Note that the time, tp∗− , to reach the maximum of y− (t) from the initial relay switch point in a negative half period of the relay can be measured, as shown in Fig. 3. It follows that tp− = tp∗− − 

(2.9)

Likewise, there exists tp+ = tp∗+ − 

(2.10)

tp∗+

where is the time to reach the minimum of y+ (t) from the initial relay switch point in a positive half period of the relay, as shown in Fig. 3. Subtracting (2.9) from (2.10) yields tp+ − tp− = tp∗+ − tp∗−

(2.11)

By substituting (2.7) and (2.8) into (2.11), we obtain ln

tp∗+ − tp∗− 0 −  1 − e−p− /p + ln = 1 − e−p+ /p p 0 + 

(2.12)

Here define z = x/P− , a = P− /P+ , and f (x) =

1 − e−P− /x , 1 − e−P+ /x

x ∈ (0, ∞)

(2.13)

It follows that f (z) =

1 − e−1/z 1 − e−1/az

(2.14)

df (z) e−(a+1)/az (e1/z − a e1/az + a − 1) = 2 dz az 2 (1 − e−1/az )

(2.15)

g(z) = e1/z − a e1/az + a − 1

(2.16)

Let

It can be verified that dg(z) 1 = 2 (e1/az − e1/z ) < 0 dz z

(2.17)

lim g(z) = 0

(2.18)

z→∞

Therefore, we can conclude that g(z) > 0 for z ∈ (0,∞) and df (z) >0 dz

(2.19)

Hence, f(z) increases monotonically with respect to z and so does for f(x) with respect to x. Consequently, it can be concluded that the left-hand side of (2.12) increases monotonically with respect to  p , and meanwhile, the right-hand side of (2.12) obviously decreases monotonically with respect to  p . Note that for  p ∈ (0,+∞), there exists 0 < ln

1 − e−p− /p p− < ln 1 − e−p+ /p p+

(2.20)

We can thus confirm that there exist finite solutions of  p for (2.12), which can be solved using any numerical algorithm such as the Newton–Raphson method.

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

3043

By substituting the SOPDT model of Eq. (2) into the process response fitting condition at the oscillation frequency as shown in (1.7), we obtain ωu

kp



p2 ωu2 + 1

−ωu −

= Au

(2.21)

 − arctan(p ωu ) = ϕu , 2

ϕu ∈



−,

− 2

 (2.22)

It can be derived from (2.21) and (2.22) that



p2 ωu2 + 1

kp = Au ωu =−

(2.23)





 1 ϕu + + arctan(p ωu ) ωu 2

(2.24)

It should be noted that for the case that / p > 1, y+ (t) may decrease monotonically for t ∈ [0, P+ ] while y− (t) may increase monotonically for t ∈ [0, P− ]. Correspondingly, there exists tp∗+ = tp∗− = 

(2.25)

In such a case, the process time constant can be inversely derived from (2.24) as p =



1  tan −ωu − − ϕu ωu 2



(2.26)

Correspondingly, the process gain can be derived from (2.23). However, the process time constant should not be derived from (2.26) if



tan −ωu −

 − ϕu 2



<0

(2.27)

Thereby (2.25) and (2.27) may be used to check whether the process time delay should be derived from (2.26). The identification of a SOPDT integrating process model can be summarized in the following algorithm named Algorithm-I-2: Measure P+ , P− , tp∗+ and tp∗− from the limit cycle; Compute G(jωu ) from (1.7); Compute the process time delay from (2.25) as  = tp∗+ , and then check whether (2.27) is satisfied. If yes, go to step (vi); Compute the process time constant,  p , from (2.26), and then check whether / p < 1 is satisfied. If yes, go to step (vi); Compute the process gain, kp , from (2.23). This is the end of the algorithm if both (2.27) and / p < 1 are not satisfied. Compute the process time constant,  p , from (2.12) using the Newton–Raphson iteration method. The initial estimation of  p for iteration may be taken as tp∗+ (or tp∗− ); (vii) Compute the process gain, kp , from (2.23); (viii) Compute the process time delay, , from (2.24); (i) (ii) (iii) (iv) (v) (vi)

(ix) End the algorithm if the fitting constraint of relay response,

N



2



[y(kTs + tos ) − y(kTs + tos )] /N < ε, is satisfied, where y(kTs + tos ) and

k=1

y(kTs + tos ) denotes respectively the model and process responses in the limit cycle, Ts the sampling period with N = Pu /Ts , and ε is a user-specified fitting threshold that may be set between 0.01–1%. Otherwise, change the initial estimation of  p and then go to step (vi).

Fig. 4. Relay response depiction for a FOPDT unstable process.

3044

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

2.3. FOPDT unstable process model For a first-order unstable process modeled by Eq. (4), the following proposition gives a limiting condition for relay feedback identification and the exact relay response expression for steady oscillation: Proposition 3. process is  < min p

Under a biased relay feedback test as shown in Fig. 4, a limiting condition to form steady oscillation for a first-order unstable

 ln

2kp 0 2kp 0 , ln kp (0 − ) + ε+ kp ( + 0 ) − ε−

 ,

(3.1)

and the resulting limit cycle is characterized by y+ (t) = −kp ( + 0 ) + y− (t) = kp (0 − ) +

2kp 0 (1 − eP− /p ) et/p , 1 − ePu /p

2kp 0 (eP+ /p − 1)et/p , 1 − ePu /p

t ∈ [0, P+ ] t ∈ [0, P− ]

(3.2) (3.3)

where y+ (t) increases monotonically for t ∈ [0, P+ ] while y− (t) decreases monotonically for t ∈ [0, P− ]. Proof.

See Appendix C. 

Remark 2. If an unbiased relay feedback test is used, we can obtain the limiting condition for steady oscillation by substituting  = 0 and ε+ = ε− = 0 into (3.1), i.e., / p < ln 2, which is just the limiting condition reported in the recent literature (Majhi & Atherton, 2000; Tan et al., 1998; Thyagarajan & Yu, 2003). It follows from Proposition 3 that the process time delay can be directly measured as the time to reach the peak of the process output response from the initial relay switch point in a half period of the relay, which is denoted as tp∗ in Fig. 4. Correspondingly, it can be obtained from (3.2) and (3.3) that y+ (P+ ) = −kp ( + 0 ) + y− (P− ) = kp (0 − ) +

2kp 0 eP+ /p (1 − eP− /p ) = A+ 1 − ePu /p

2kp 0 eP− /p (eP+ /p − 1) = A− 1 − ePu /p

y+ (P+ − ) = −kp ( + 0 ) + y− (P− − ) = kp (0 − ) +

2kp 0 e(P+ −)/p (1 − eP− /p ) = −ε− 1 − ePu /p

2kp 0 e(P− −)/p (eP+ /p − 1) = −ε+ 1 − ePu /p

(3.4) (3.5) (3.6) (3.7)

Substituting (3.4) into (3.6), we obtain p =

 ln{[kp (0 + ) + A+ ]/[kp (0 + ) − ε− ]}

(3.8)

Alternatively, substituting (3.5) into (3.7) yields p =

 ln{[kp (0 − ) − A− ]/[kp (0 − ) + ε+ ]}

(3.9)

Note that the process gain, kp , can be derived as

 tos +Pu t

kp = −G(0) = −  tos +P os

tos

u

y(t) dt u(t) dt

(3.10)

The process time constant can then be derived from (3.8) or (3.9). It is preferred in practice to use (3.8) for better accuracy, in view of the fact that the positive fitting part in a half period of the limit cycle occupies a larger percentage compared to the negative fitting part in the other half period. Hence, the identification of a FOPDT unstable process model can be summarized in the following algorithm named Algorithm-U-1a: (i) Measure P+ , P− and A+ from the limit cycle; (ii) Measure the process time delay as tp∗ , which is the time taken to reach the positive peak (A+ ) of the process output response from the initial relay switch point in a negative half period of the relay (P− ); (iii) Compute the process gain, kp , from (3.10); (iv) Compute the process time constant,  p , from (3.8).

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

3045

Fig. 5. Relay response depiction for a SOPDT unstable process.

Note that by substituting the FOPDT unstable process model into the process response fitting condition at the oscillation frequency as shown in (1.7), we can obtain



1 p = ωu =−

kp2 A2u

−1

(3.11)

1 [ϕu +  − arctan(p ωu )] ωu

(3.12)

Therefore, an alternative identification of a FOPDT unstable process model can be summarized in the following algorithm named Algorithm-U-1b: (i) (ii) (iii) (iv) (v)

Measure Pu from the limit cycle; Compute G(jωu ) from (1.7); Compute the process gain, kp , from (3.10); Compute the process time constant,  p , from (3.11); Compute the process time delay, , from (3.12).

2.4. SOPDT unstable process model For a second-order unstable process modeled by Eq. (5), the following proposition gives the exact relay response expression for steady oscillation: Proposition 4. Under a biased relay feedback test as shown in Fig. 5, the resulting limit cycle for a second-order unstable process is characterized by y+ (t) = −kp ( + 0 ) + y− (t) = kp (0 − ) +

2kp 0 1 E1 t/ 2kp 0 2 E2 −t/ 2, e 1+ e 1 + 2 1 + 2

2kp 0 1 F1 t/ 2kp 0 2 F2 −t/ 2, e 1+ e 1 + 2 1 + 2

t ∈ [0, P+ ] t ∈ [0, P− ]

(4.1) (4.2)

where 0 < E1 < E2 , F2 < F1 < 0, and E1 =

1 − ep− /1 1 − epu /1

(4.3)

E2 =

1 − e−p− /2 1 − e−pu /2

(4.4)

F1 = −

1 − ep+ /1 1 − epu /1

(4.5)

F2 = −

1 − e−p+ /2 1 − e−pu /2

(4.6)

Proof.

See Appendix D.



It can be derived from (4.1) and (4.2) that 2kp 0 dy+ (t) (E1 et/1 − E2 e−t/2 ) = 1 + 2 dt

(4.7)

3046

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

2kp 0 dy− (t) (F1 et/1 − F2 e−t/2 ) = 1 + 2 dt

(4.8)

By letting dy+ (t)/dt = 0 and dy− (t)/dt = 0, the time to reach the single extreme value of y+ (t) and y− (t) can be derived, respectively, as tp+ =

1 2 E2 ln 1 + 2 E1

(4.9)

tp− =

1 2 F2 ln 1 + 2 F1

(4.10)

It can be easily verified that d2 y+ (t)/dt2 > 0 and d2 y− (t)/dt2 < 0. Therefore, tp+ is the time to reach the minimum of y+ (t), while tp− is the time to reach the maximum of y− (t). Note that the time, tp∗− , to reach the maximum of y− (t) from the initial relay switch point in a negative half period of the relay can be measured as shown in Fig. 5, and so does for tp∗+ to reach the minimum of y+ (t) from the initial relay switch point in a positive half period of the relay. Therefore, it follows that tp− = tp∗− − 

(4.11)

tp+ = tp∗+ − 

(4.12)

Subtracting (4.11) from (4.12) yields tp+ − tp− = tp∗+ − tp∗−

(4.13)

By substituting the SOPDT model into the process response fitting condition at the oscillation frequency as shown in (1.7), we obtain



kp (12 ωu2

+ 1)(22 ωu2 + 1)

= Au

−ωu −  + arctan(1 ωu ) − arctan(2 ωu ) = ϕu ,

(4.14)

ϕu ∈



−,

− 2

 (4.15)

It can be derived from (4.14) and (4.15) that



1 2 = ωu =−

kp2

−1

(4.16)

1 [ϕu +  − arctan(1 ωu ) + arctan(2 ωu )] ωu

(4.17)

A2u (12 ωu2 + 1)

Note that the process gain, kp , can be derived using (3.10). By substituting (4.9), (4.10) and (4.16) into (4.13), we can obtain a transcendental equation with respect to  1 . This equation can be numerically solved with the Newton–Raphson iteration method. The initial estimation of  1 for iteration may be taken as



1 1 = ωu

kp2 A2u

−1

(4.18)

which has been given in Algorithm-U-1b for computing the single time constant of a FOPDT unstable process model. Therefore, the identification of a SOPDT unstable process model can be summarized in the following algorithm named Algorithm-U-2: Measure P+ , P− , tp∗+ and tp∗− from the limit cycle; Compute G(jωu ) from (1.7); Compute the process gain, kp , from (3.10); Compute the process time constant,  1 , from the equation resulting from substituting (4.9), (4.10) and (4.16) into (4.13) by using the Newton–Raphson iteration method. The initial estimation of  1 for iteration may be taken from (4.18); (v) Compute the process time constant,  2 , from (4.16); (vi) Compute the process time delay, , from (4.17);

(i) (ii) (iii) (iv)

(vii) End the algorithm if the fitting constraint of relay response, initial estimation of  1 and then go to step (iv).

N



2

[y(kTs + tos ) − y(kTs + tos )] /N < ε, is satisfied. Otherwise, change the

k=1

3. Measurement noise issue Measurement noises are often unavoidable in process operation. In the context of process identification, noise-to-signal ratio (NSR) is commonly evaluated as NSR =

mean(abs(noise)) mean(abs(signal))

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

3047

To reduce the influence from measurement noise, it is suggested that the measured values of P+ , P− , A+ and A− for 10–20 periods in steady oscillation may be respectively averaged as P¯ + , P¯ − , A¯ + and A¯ − for computation in the proposed identification algorithms. The oscillation period can thus be obtained as P¯ u = P¯ + + P¯ − . Accordingly, the process response at the oscillation frequency can be computed as

 tos +N P¯ u tos

G(jω ¯ u) =

 tos +N P¯ u tos

y(t) e−jω¯ u t dt

= A¯ u ejϕ¯ u

u(t) e−jω¯ u t dt

where ω ¯ u = 2/P¯ u , N is the number of periods used for averaging that may be taken in the range of 10–20, and tos may be chosen as any relay switch point in steady oscillation. Note that the relay output, u(t), remains as a constant for each half periods, so it may be used for measurement of the oscillation period. For FOPDT and SOPDT unstable processes, the process gain can be derived as

 tos +N P¯ u

kp = − 

tos tos +N P¯ u tos

y(t) dt u(t) dt

Note that the above usage of the measured parameters is based on the statistical principle for eliminating the random measurement errors, and therefore, is capable of identification robustness for low noise level (e.g., NSR < 10%). Besides, in view of that



Pu

e

−jωu t

Pu

dt = 

0

e−jωu t dt = 0

0

where  ∈ R, we can ensure that

 tos +Pu

G(jωu ) =

tos

 tos +Pu tos

y(t) e−jωu t dt

 tos +Pu =

u(t) e−jωu t dt

tos

(y(t) + ) e−jωu t dt

 tos +Pu tos

u(t) e−jωu t dt

Hence, the computation of the process response at the oscillation frequency is indeed not affected by static load disturbances. The influence of nonstatic load disturbances can be intuitively excluded by comparing the identity between the sequential periods of the limit cycle used for the above computation. To cope with measurement noises causing high noise level, a low-pass digital Butterworth filter is suggested to recover the corrupted limit cycle data. It may be determined by specifying the filter order, nf , and the cutoff frequency, fc , i.e., Butter (nf , fc ) =

b0 + b1 z −1 + b2 z −2 + · · · + bnf z −nf 1 + a1 z −1 + a2 z −2 + · · · + anf z −nf

where Butter (nf , fc ) denotes the filter function with two input parameters of nf and fc . In view of that measurement noise is mainly of high-frequency, the guideline for choosing the cutoff frequency is fc ≥ 5fu =

5 Pu

Thereby, only the signal components within the frequency band around ωu (i.e., 2fu ) can be passed through. Note that the phase lag caused by the low-pass filter almost does not affect measurement of the oscillation amplitude and period of the limit cycle, because the relay output has a similar phase lag under the filtered output feedback. 4. Illustration Six examples studied in the recent literature are performed here to illustrate the effectiveness and accuracy of the proposed identification algorithms. Examples 1–4 are given to demonstrate the achievable accuracy of the proposed algorithms for identifying first- and secondorder integrating and unstable processes, respectively. Measurement noise tests are introduced in Example 4 to illustrate identification robustness of the proposed algorithms. Examples 5 and 6 are given to show the effectiveness of the proposed algorithms for identification of higher order integrating and unstable processes. Table 1 Identified models for several integrating and unstable processes Process

Limit cycle P+

e−5s s e−10s s(20s+1) −0.4s e s−1 e−0.5s (2s−1)(0.5s+1) (−s+1) e−5s s(s+1)5 e−0.5s (5s−1)(2s+1)(0.5s+1)

Integrating/unstable process model P−

A+

A−

10.8

10.8

2.7

−2.7

42.9

60.06

9.1041

−6.4432

ϕu

FOPDT

3.4377

−3.0267

1.0000 e−5.0000s s

10.3853

−3.0657

Au

1.0001 e−0.4s

1.1

3.17

0.7432

−0.5532

0.5621

−2.7639

4.24

6.35

0.8569

−0.6441

0.6177

−2.8587

17.99

29.99

4.9569

−2.8459

7.3798

−3.0074

0.9954s−1 1.0001e−1.0486s 2.1459s−1 0.9664 e−10.9703s s

13.76

18.88

0.8698

−0.7199

0.6694

−2.9359

1.0001 e−3.2821s 5.7663s−1

SOPDT 0.9983 e−10.027s s(19.9443s+1) 1.0001 e−0.5051s (1.9975s−1)(0.4988s+1) 1.018 e−8.5278s s(2.5293s+1) 1.0001 e−0.9422s (4.9996s−1)(2.07s+1)

3048

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

Fig. 6. Nyquist plots of FOPDT models for Example 4.

Example 1. G=

Consider the FOPDT integrating process studied by Majhi and Atherton (2000)

e−5s s

Majhi and Atherton (2000) gave a FOPDT model with each parameter error within 0.1%. By performing an unbiased relay test with u+ = −u− = 0.5 and ε+ = −ε− = 0.2, the proposed Algorithm-I-1a results in the exact process model as listed in Table 1. It should be noted that the proposed Algorithm-I-1b can also obtain the exact process model, which may be verified using the measured parameters of the limit cycle for computation as listed in Table 1. Example 2. G=

Consider the SOPDT integrating process studied by Kaya (2006) and Ho et al. (1996)

e−10s s(20s + 1)

Kaya (2006) derived an almost exact SOPDT model, showing quite improved accuracy in comparison with Ho et al. (1996). By performing a biased relay test with u+ = 0.7, u− = 0.5 and ε+ = −ε− = 0.1 as in Kaya (2006), the proposed Algorithm-I-2 yields the process model listed in Table 1, indicating good identification accuracy. Example 3. G=

Consider the FOPDT unstable process widely studied in the literature (e.g., Marchetti et al., 2001; Shen, Wu, & Yu, 1996)

e−0.4s s−1

Based on two different relay tests, Marchetti et al. (2001) derived a FOPDT model, Gm = 0.928 e−0.392s /(0.757s − 1), showing further improved accuracy in comparison with Shen et al. (1996). By performing a biased relay test with u+ = 1.2, u− = −0.8 and ε+ = −ε− = 0.1, the proposed Algorithm-U-1a yields a model listed in Table 1, indicating high accuracy. It should be noted that the proposed Algorithm-U-1b can also obtain the exact process model using the measured parameters of the limit cycle shown in Table 1. Example 4. Consider the SOPDT unstable process studied by Vivek and Chidambaram (2005) and Ramakrishnan and Chidambaram (2003)

G=

e−0.5s (2s − 1)(0.5s + 1)

Using an unbiased relay test, Vivek and Chidambaram (2005) derived a FOPDT model, Gm−1 = 0.7534 e−1.0412s /(2.1642s − 1), and Ramakrishnan and Chidambaram (2003) gave a SOPDT model, Gm−2 = e−0.52s /(1.9999s − 1)(0.4837s + 1), from a biased relay test. By performing an input biased relay test with u+ = u− = 1.0, ε+ = 0.1 and ε− = −0.2, the proposed Algorithm-U-1a yields a FOPDT model, Gm = 1.0001 e−1.82s /(4.1688s − 1), and the proposed Algorithm-U-1b and Algorithm-U-2 result in FOPDT and SOPDT models listed in Table 1. For comparison, the Nyquist plots of the FOPDT models obtained from the proposed algorithms and Vivek and Chidambaram (2005) are shown in Fig. 6. It is seen that enhanced fitting is captured by the proposed algorithms, while the proposed Algorithm-U-1b procures better fitting over the low frequency range owing to that the model response coincides with the process at the oscillation frequency, i.e. (−0.5932, −j0.1724). Note that Vivek and Chidambaram (2005) gave a FOPDT model, Gm−1 = 0.5332 e−0.8807s /(1.3524s − 1), for the case that a random noise with zero mean and standard deviation of 1% is added to the process output measurement and then used for relay feedback. It can be verified that the measurement noise causes NSR = 1.6%. By using the measured parameters in 10 periods of steady oscillation, the proposed Algorithm-U-1a and Algorithm-U-1b yield FOPDT models, Gm = 0.9992 e−1.806s /(4.0701s − 1) and Gm = 0.9992 e−1.047s /(2.1449s − 1), respectively. The identification errors resulting from the proposed Algorithm-U-2 are listed in Table 2, indicating good identification robustness. To further demonstrate identification robustness, assume that the noise level is increased to NSR = 10% for measurement of the limit cycle, it can be seen from Table 2 that notable errors occur for measuring the positive and negative magnitudes of the limit cycle, but the computation of the process response at the oscillation frequency is less affected owing to the use of all the measured data in the 10 oscillation

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

3049

Fig. 7. Recovered limit cycle for Example 4 with NSR=20%.

Fig. 8. Nyquist plots of identified models for Example 5.

periods. Accordingly, only the stable time constant of the process is apparently underestimated. To enhance the identification accuracy, a third-order digital Butterworth filter with a cutoff frequency fc = 0.5 (HZ) according to the guideline given in Section 3 is therefore used to recover the corrupted limit cycle data. The corresponding results listed in Table 2 show that obviously improved accuracy for the measured parameters and model identification can thus be obtained. For illustration, the recovery effect by using the above filter in both the forward and reverse directions is shown in Fig. 7 for the case that a random noise causing NSR = 20% is added to the process output measurement except for relay feedback. It can be seen that the corrupted limit cycle has been well recovered against the severe noise level, leading to the small identification errors listed in Table 2. Example 5. G=

¨ Consider the high order integrating process studied by Ingimundarson and Hagglund (2001)

(−s + 1) e−5s s(s + 1)5

¨ Based on the process step response data, Ingimundarson and Hagglund (2001) derived a FOPDT model, Gm = 1.0000 e−11.0000s /s. By performing a biased relay test with u+ = 0.5, u− = −0.3 and ε+ = −ε− = 0.1, the proposed Algorithm-I-1a yields a FOPDT model, Gm = 0.8675 e−11.25s /s, and the proposed Algorithm-I-1b and Algorithm-I-2 result in the FOPDT and SOPDT models listed in Table 1. For comparison, the Nyquist plots of these models are shown in Fig. 8. It can be seen that the FOPDT model obtained from the proposed Algorithm-I-1b procures better fitting compared to the other FOPDT models, owing to that the model response coincides with the process Table 2 Identification of Example 4 against measurement noise NSR

Denoising

1.6% 10% 10% 20%

Averaging Averaging Filtering Filtering

% error in the measured values

% error in the model parameters

A+

A−

Au

ϕu

kp

1

2



1.49 20.35 1.96 2.21

−2.6 −26.19 −0.97 −1.86

−0.96 1.39 0.58 1.34

−0.05 0.47 0.41 0.58

−0.08 −0.02 −0.01 0.02

0.21 −0.82 0.76 0.87

−2.44 −11.8 −0.78 −6.04

3.52 5.94 1.94 4.36

3050

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

at the oscillation frequency, i.e. (−7.3134, −j0.9873). In contrast with these FOPDT models, the proposed SOPDT model shows apparently improved fitting effect. Example 6. G=

Consider the high order unstable process studied by Majhi (2007)

e−0.5s (5s − 1)(2s + 1)(0.5s + 1)

Majhi (2007) derived a SOPDT model, Gm = 1.001 e−0.939s /(10.354s2 + 2.932s − 1) from a state-space analysis on the process response under relay feedback. By performing a biased relay test with u+ = −u− = 1.0, ε+ = 0.1 and ε− = −0.15, the proposed Algorithm-U-1a yields a FOPDT model, Gm = 1.0001 e−6.55s /(13.4749s − 1), and the proposed Algorithm-U-1b and Algorithm-U-2 result in FOPDT and SOPDT models listed in Table 1. For comparison, the Nyquist plots of these models are shown in Fig. 9. It is seen that both of the SOPDT models obtained from the proposed Algorithm-U-2 and Majhi (2007) obtain very close fitting. Compared to Algorithm-U-1a, Algorithm-U-1b again shows better fitting over the low frequency range owing to that the model response coincides with the process at the oscillation frequency, i.e. (−0.6553, −j0.1367). 5. Conclusions In view of that FOPDT and SOPDT models are widely adopted in engineering practice for model-based control design for integrating and unstable processes, the exact relay response expressions for these low-order process models under a relay feedback test have been derived in this paper for evaluation of the relay response and choice of a suitable model structure. It has been therefore demonstrated that the limit cycle can be definitely formed for an integrating process, and a tighter limiting condition to form steady oscillation for unstable processes has also been cultivated. The corresponding identification algorithms have been developed transparently using the measurable parameters in steady oscillation and fitting conditions of the limit cycle. To guarantee identification robustness, a statistical averaging method has been presented for measurement against low noise level, and a low-pass Butterworth filter has been suggested to recover the corrupted limit cycle data from high noise level. The applications to several examples studied in the recent literature have verified the effectiveness and accuracy of the proposed identification algorithms, and also demonstrated that a SOPDT model thus obtained evidently outperforms the FOPDT model for frequency response fitting for a real high order integrating or unstable process. Acknowledgement This work is supported in part by the Hong Kong Research Grants Council under Project No. 612906. Appendix A. Proof of Proposition 1 The initial step response of a FOPDT integrating process arising from the relay output  − 0 can be obtained as y0 (t) = kp ( − 0 )(t − )

(I.1)

When it comes to the first relay switch point denoted as t0 in Fig. 2, the relay output changes to be  + 0 , indicating that a step change of 2 0 is added to the process input. Using the linear superposition principle, the process output response can be derived as y1 (t) = y0 (t + t0 ) + 2kp 0 (t − )

(I.2)

By using a time shift of t0 + , (I.2) can be rewritten as y1 |shift (t) = y0 (t + t0 + ) + 2kp 0 t

(I.3)

Fig. 9. Nyquist plots of identified models for Example 6.

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

3051

When it comes to the second relay switch point, the relay output changes to be  − 0 , indicating that a step change of −20 is added to the process input. Again using the linear superposition principle, the process output response can be obtained as y2 (t) = y1 (t + P+ ) − 2kp 0 (t − )

(I.4)

Using a time shift of t0 +  + P+ (I.4) can be rewritten as y2 |shift (t) = y0 (t + t0 +  + P+ ) + 2kp 0 t(1 − 1) + 2kp 0 P+

(I.5)

At the third relay switch point, the relay output changes again to be  + 0 , indicating that a step change of 20 is once again added to the process input. The process output response can thus be obtained as y3 (t) = y2 (t + P− ) + 2kp 0 (t − )

(I.6)

Using a time shift of t0 +  + Pu (I.6) can be rewritten as y3 |shift (t) = y0 (t + t0 +  + Pu ) + 2kp 0 t(1 − 1 + 1) + 2kp 0 P+

(I.7)

The process output response following the fourth relay switch point, is the result of four interlaced step changes respectively with a magnitude of 20 . Using a time shift of t0 +  + Pu + P+ , the process output response can be written as y4 |shift (t) = y0 (t + t0 +  + Pu + P+ ) + 2kp 0 t(1 − 1 + 1 − 1) + 2kp 0 · 2P+

(I.8)

Hence, the time shifted output response after each relay switch point can be generally expressed as



y2n+1 

 y2n+2 

shift

(t) = y0 (t + t0 +  + nPu ) + 2kp 0 t + 2kp 0 · nP+

(I.9)

shift

(t) = y0 (t + t0 +  + nPu + P+ ) + 2kp 0 · (n + 1)P+

(I.10)

where n = 0,1,2,. . .. By substituting (I.1) into (I.9) and (I.10), respectively, we obtain



y2n+1 

 y2n+2 

shift

(t) = kp ( + 0 )t + kp ( − 0 )t0 + kp n[( − 0 )Pu + 20 P+ ]

(I.11)

shift

(t) = kp ( − 0 )(t + t0 ) + kp ( + 0 )P+ + kp n[( − 0 )Pu + 20 P+ ]

(I.12)

Note that n → ∞ as t → ∞. The condition for the process output response to move into steady oscillation is therefore required as ( − 0 )Pu + 20 P+ = 0

(I.13)

Hence, in the limit cycle, it follows that y+ (t) = kp ( + 0 )t + kp ( − 0 )t0 ,

t ∈ [0, P+ ]

y− (t) = kp ( − 0 )(t + t0 ) + kp ( + 0 )P+ ,

t ∈ [0, P− ]

(I.14) (I.15)

where y+ (t) denotes the ascending part in the half period P+ as shown in Fig. 2, while y− (t) indicates the descending part in the other half period P− . It can be easily seen from (I.14) and (I.15) that y+ (t) increases monotonically for t ∈ [0, P+ ] while y− (t) decreases monotonically for t ∈ [0, P− ], corresponding to t ∈ [P+ , Pu ] in the limit cycle. Note that there exists an inherent connection between y+ (t) and y− (t), y+ (0) = y− (P− )

(I.16)

which can be reformulated using (I.14) and (I.15) as ( − 0 )P− + ( + 0 )P+ = 0

(I.17)

It can thus be verified that (I.16) coincides with (I.13). Accordingly, we can conclude that under a conventional relay feedback test, the limit cycle can be definitely formed for a FOPDT integrating process. In particular, there exist  = 0 and Pu = 2P+ when an unbiased relay feedback test is used, so that (I.13) is obviously satisfied for ∀0 . It should be noted that before the process output response moves into the limit cycle, the time intervals between each relay switch points may not equal the corresponding half periods in the limit cycle. Nevertheless, we can equalize them as done in (I.4)–(I.8) to derive the exact expression of the limit cycle, in view of that the limit cycle does not include these initial response characteristics. In other words, the process output response in steady oscillation has the same limit cycle with that of the corresponding ideal oscillation which has identical time intervals between the sequential relay switch points from beginning to end. Appendix B. Proof of Proposition 2 The initial step response of a SOPDT integrating process arising from the relay output  − 0 can be derived as y0 (t) = kp ( − 0 )[t −  + p (e−(t−)/p − 1)]

(II.1)

Following a similar analysis as in Appendix A from the initial to the fourth relay switch point, the time shifted output responses after each relay switch point can be respectively derived as y1 |shift (t) = y0 (t + t0 + ) + 2kp 0 (t − p ) + 2kp 0 p e−t/p

(II.2)

3052

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

y2 |shift (t) = y0 (t + t0 +  + P+ ) + 2kp 0 (t − p ) · (1 − 1) + 2kp 0 P+ + 2kp 0 p e−t/p (e−P+ /p − 1)

(II.3)

y3 |shift (t) = y0 (t + t0 +  + Pu ) + 2kp 0 (t − p ) · (1 − 1 + 1) + 2kp 0 P+ + 2kp 0 p e−t/p (e−Pu /p − e−P− /p + 1)

(II.4)

y4 |shift (t) = y0 (t + t0 +  + Pu + P+ ) + 2kp 0 (t − p ) · (1 − 1 + 1 − 1) +2kp 0 · 2P+ + 2kp 0 p e−t/p (e−(Pu +P+ )/p − e−Pu /p + e−P+ /p − 1)

(II.5)

The general relay response can be therefore summarized as



y2n+1 



y2n+2 

shift

(t) = y0 (t + t0 +  + nPu ) + 2kp 0 (t − p ) + 2kp 0 · nP+ + 2kp 0 p E e−t/p

(II.6)

shift

(t) = y0 (t + t0 +  + nPu + P+ ) + 2kp 0 · (n + 1)P+ + 2kp 0 p F e−t/p

(II.7)

where n = 0,1,2,. . ., and E =1+

n 

(e−kPu /p − e−[(k−1)Pu +P− ]/p )

(II.8)

k=1

F=

n 

(e−(kPu +P+ )/p − e−kPu /p )

(II.9)

k=0

Note that 0 < e−Pu /p < 1. The infinite series in (II.8) and (II.9) are uniformly convergent as n → ∞, that is, n 

e−kPu /p =

k=0

1 1 − e−Pu /p

(II.10)

Substituting (II.10) into (II.8) and (II.9), respectively, we obtain E=

1 e−p− /p 1 − e−p− /p − = 1 − e−pu /p 1 − e−pu /p 1 − e−pu /p

(II.11)

F=

e−p+ /p − 1 1 − e−p+ /p =− −p / u p 1−e 1 − e−pu /p

(II.12)

Note that lim e−(t+t0 +nPu )/p = lim e−(t+t0 +nPu +P+ )/p = 0

n→∞

(II.13)

n→∞

By substituting (II.1) and (II.13) into (II.6) and (II.7), respectively, we obtain



y2n+1 

 y2n+2 

shift

(t) = kp ( + 0 )(t − p ) + kp ( − 0 )t0 + kp n[( − 0 )Pu + 20 P+ ] + 2kp 0 p E e−t/p

(II.14)

shift

(t) = kp ( − 0 )(t + t0 − p ) + kp ( + 0 )P+ + kp n[( − 0 )Pu + 20 P+ ] + 2kp 0 p F e−t/p

(II.15)

In view of that n → ∞ as t → ∞, the condition for the output response to move into steady oscillation is therefore required as ( − 0 )Pu + 20 P+ = 0

(II.16)

Note that (II.16) is the same as (I.13) for a FOPDT integrating process. Hence, in the limit cycle, it follows that y+ (t) = kp ( + 0 )(t − p ) + kp ( − 0 )t0 + 2kp 0 p E e−t/p ,

t ∈ [0, P+ ]

y− (t) = kp ( − 0 )(t + t0 − p ) + kp ( + 0 )P+ + 2kp 0 p F e−t/p

t ∈ [0, P− ]

(II.17) (II.18)

where y+ (t) denotes the ascending output response in the half period P+ , due to a positive step change of the relay output as shown in Fig. 3, while y− (t) denotes the descending output response in the other half period P− , due to a negative step change of the relay output. According to the inherent connection between y+ (t) and y− (t), there exists y+ (0) = y− (P− )

(II.19)

By substituting (II.17) and (II.18) into (II.19), we obtain ( − 0 )P− + ( + 0 )P+ = 0

(II.20)

Thereby, it is seen that (II.20) coincides with (II.16), indicating that the limit cycle can be definitely formed for a SOPDT integrating process under biased/unbiased relay feedback. Following a similar analysis, the above conclusion can be transparently generalized to a high order integrating process.

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

3053

Appendix C. Proof of Proposition 3 A unity step response of a FOPDT unstable process can be derived as y(t) = kp (e(t−)/p − 1)

(III.1)

It can be seen from Fig. 4 that according to the initial relay output  − 0 , there exist y0 (t) = kp ( − 0 )(e(t−)/p − 1) t ∈ [, t0 + ]

(III.2)

y0 (t0 ) = −ε+

(III.3)

By using (III.3), we can obtain y0 (t0 + ) = −kp ( − 0 ) + [kp ( − 0 ) − ε+ ] e/p

(III.4)

After the time of t0 + , the process output begins to respond to the switched relay output  + 0 in the form of y(t) = y0 (t0 + ) e(t−t0 −)/p + kp ( + 0 )(e(t−t0 −)/p − 1),

t ∈ (t0 + , t1 + ]

(III.5)

It follows from (III.5) that dy(t) 1 = [y0 (t0 + ) + kp ( + 0 )] e(t−t0 −)/p p dt

(III.6)

To form steady oscillation, y(t) should have an ascending tendency in the time interval (t0 + , t1 + ), thus requiring that y(t0 + ) + kp ( + 0 ) > 0

(III.7)

By substituting (III.4) into (III.7), we can obtain 2kp 0  < ln p kp (0 − ) + ε+

(III.8)

Analogously, it follows from (III.5) and y(t1 ) = −ε+ that y(t1 + ) = −kp ( + 0 )[kp ( + 0 ) − ε− ] e/p

(III.9)

It can be seen from Fig. 4 that after the time of t1 + , the process output begins to respond to the switched relay output  − 0 in the form of y(t) = y(t1 + ) e(t−t1 −)/p + kp ( − 0 )(e(t−t1 −)/p − 1),

t ∈ (t1 + , t2 + ]

(III.10)

It follows from (III.10) that dy(t) 1 = [y(t1 + ) + kp ( − 0 )] e(t−t1 −)/p p dt

(III.11)

To form steady oscillation, y(t) should have an descending tendency in the time interval (t1 + , t2 + ], thus requiring that y(t1 + ) + kp ( − 0 ) < 0

(III.12)

By substituting (III.9) into (III.12), we can obtain 2kp 0  < ln p kp (0 + ) − ε−

(III.13)

Combining (III.8) with (III.13), we can obtain the limiting condition to form steady oscillation as shown in Proposition 3. In the sequel, following a similar analysis as in Appendix A from the initial to the fourth relay switch point, the time shifted process output responses after each relay switch point can be respectively derived as y1 |shift (t) = y0 (t + t0 + ) − 2kp 0 + 2kp 0 et/p

(III.14)

y2 |shift (t) = y0 (t + t0 +  + P+ ) − 2kp 0 (1 − 1) + 2kp 0 et/p (eP+ /p − 1)

(III.15)

y3 |shift (t) = y0 (t + t0 +  + Pu ) − 2kp 0 (1 − 1 + 1) + 2kp 0 et/p (ePu /p − eP− /p + 1) t/p

y4 |shift (t) = y0 (t + t0 +  + Pu + P+ ) − 2kp 0 (1 − 1 + 1 − 1) + 2kp 0 e

(Pu +P+ )/p

(e

(III.16) Pu /p

−e

P+ /p

+e

− 1)

(III.17)

The general relay response can be therefore summarized as



y2n+1 

 y2n+2 

shift

(t) = y0 (t + t0 +  + nPu ) − 2kp 0 + 2kp 0 E et/p

(III.18)

shift

(t) = y0 (t + t0 +  + nPu + P+ ) + 2kp 0 F et/p

(III.19)

where n = 0,1,2,. . ., and E =1+

n  k=1

(ekPu /p − e[(k−1)Pu +P− ]/p ) = 1 + (ePu /p − eP− /p ) ·

1 − e(n−1)Pu /p 1 − ePu /p

(III.20)

3054

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

F=

n 

(e(kPu +P+ )/p − ekPu /p ) = (eP+ /p − 1) ·

k=0

1 − enPu /p 1 − ePu /p

(III.21)

It follows from (III.2) that y0 (t + t0 +  + nPu ) = kp ( − 0 )(e(t+t0 +nPu )/p − 1)

(III.22)

(t+t0 +nPu +P+ )/p

y0 (t + t0 +  + nPu + P+ ) = kp ( − 0 )(e

− 1)

(III.23)

Substituting (III.20)–(III.23) into (III.18) and (III.19), respectively, we obtain



y2n+1 

(III.24)

kp et/p {enPu /p [( − 0 )(1 − ePu /p ) e(t0 +P+ )/p − 20 (eP+ /p − 1)] + 20 (eP+ /p − 1)} 1 − ePu /p

(III.25)

(t) = −kp ( + 0 ) +

shift

(t) = kp (0 − ) +



y2n+2 

kp et/p {enPu /p [( − 0 )(1 − ePu /p ) et0 /p − 20 (1 − e−P+ /p )] + 20 (1 − eP− /p )} 1 − ePu /p

shift

Note that n → ∞ as t → ∞. It can be seen from (III.24) and (III.25) that there exists an inherent constraint for steady oscillation, ( − 0 )(1 − ePu /p ) et0 /p − 20 (1 − e−P+ /p ) = 0

(III.26)

Hence, in the limit cycle, it follows that y+ (t) = −kp ( + 0 ) + y− (t) = kp (0 − ) +

2kp 0 (1 − eP− /p ) et/p , 1 − ePu /p

2kp 0 (eP+ /p − 1) et/p , 1 − ePu /p

t ∈ [0, P+ ]

(III.27)

t ∈ [0, P− ]

(III.28)

where y+ (t) denotes the monotonically ascending part in the half period P+ , while y− (t) the monotonically descending part in the other half period P− . It can be verified from (III.27) and (III.28) that the inherent connection between y+ (t) and y− (t) is satisfied, i.e., y+ (0) = y− (P− )

(III.29)

y+ (P+ ) = y− (0)

(III.30)

Appendix D. Proof of Proposition 4 The initial step response of a SOPDT unstable process arising from the relay output  − 0 can be derived as y0 (t) = kp ( − 0 )





2 1 e(t−)/1 + e−(t−)/2 − 1 1 + 2 1 + 2

(IV.1)

Following a similar analysis as in Appendix A from the initial to the fourth relay switch point, the time shifted process output responses after each relay switch point can be respectively derived as y1 |shift (t) = y0 (t + t0 + ) + 2kp 0





2 1 et/1 + e−t/2 − 1 1 + 2 1 + 2

y2 |shift (t) = y0 (t + t0 +  + P+ ) − 2kp 0 · (1 − 1) +

2kp 0 1 t/ P+ / 2kp 0 2 −t/ −P+ / 1 − 1) + 2 (e 2 − 1) e 1 (e e 1 + 2 1 + 2

y3 |shift (t) = y0 (t + t0 +  + Pu ) − 2kp 0 · (1 − 1 + 1) +

(IV.3)

2kp 0 1 t/ Pu / 2kp 0 2 −t/ −Pu / 1 − eP− /1 + 1) + 2 (e 2 − e−P− /2 + 1) e 1 (e e 1 + 2 1 + 2 (IV.4)

y4 |shift (t) = y0 (t + t0 +  + Pu + P+ ) − 2kp 0 · (1 − 1 + 1 − 1) + +

(IV.2)

2kp 0 1 t/ (Pu +P+ )/ 1 − ePu /1 + eP+ /1 − 1) e 1 (e 1 + 2

2kp 0 2 −t/ −(Pu +P+ )/ 2 (e 2 − e−Pu /2 + e−P+ /2 − 1) e 1 + 2

(IV.5)

The general relay response can be therefore summarized as



y2n+1 

shift

(t) = y0 (t + t0 +  + nPu ) − 2kp 0 +

shift

(t) = y0 (t + t0 +  + nPu + P+ ) +



y2n+2 

2kp 0 E1∗ 1 1 + 2

2kp 0 F1∗ 1 1 + 2

et/1 +

et/1 +

2kp 0 E2 2 −t/ 2 e 1 + 2

2kp 0 F2 2 −t/ 2 e 1 + 2

(IV.6) (IV.7)

where n = 0,1,2,. . ., and E1∗ = 1 +

n  k=1

(ekPu /1 − e[(k−1)Pu +P− ]/1 ) = 1 +

(ePu /1 − eP− /1 )(1 − e(n−1)Pu /1 ) 1 − ePu /1

(IV.8)

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056 n 

E2 = 1 +

(e−kPu /2 − e−[(k−1)Pu +P− ]/2 ) =

k=1

F1∗ =

n 

(e(kPu +P+ )/1 − ekPu /1 ) = −

k=0

F2 =

n  k=0

y0 (t + t0 +  + nPu ) = kp ( − 0 )

1 − e−P− /2 1 − e−Pu /2

(IV.9)

(1 − eP+ /1 )(1 − enPu /1 ) 1 − ePu /1

(e−(kPu +P+ )/2 − e−kPu /2 ) = −

It follows from (IV.1) that

3055



(IV.10)

1 − e−P+ /2 1 − e−Pu /2

(IV.11)



2 1 e(t+t0 +nPu )/1 + e−(t+t0 +nPu )/2 − 1 1 + 2 1 + 2

y0 (t + t0 +  + nPu + P+ ) = kp ( − 0 )



(IV.12)



1 2 e(t+t0 +nPu +P+ )/1 + e−(t+t0 +nPu +P+ )/2 − 1 1 + 2 1 + 2

(IV.13)

Substituting (IV.8)–(IV.13) into (IV.6) and (IV.7), respectively, we obtain



y2n+1 

shift

(t) = −kp ( + 0 ) + +



y2n+2 

shift

2kp 0 1 E1 t/ e 1 1 + 2



2kp 0 2 E2 −t/ kp 1 20 (ePu /1 − eP− /1 ) 2 + e e[t+(n−1)Pu ]/1 ( − 0 ) e(t0 +Pu )/1 − 1 + 2 1 + 2 1 − ePu /1

(t) = kp (0 − ) +

(IV.14)





2kp 0 1 F1 t/ 2kp 0 2 F2 −t/ kp 1 20 (eP+ /1 − 1) 2 + e 1+ e e(t+nPu )/1 ( − 0 ) e(t0 +P+ )/1 − 1 + 2 1 + 2 1 + 2 1 − ePu /1 (IV.15)

where E1 =

1 − ep− /1 1 − epu /1

F1 = −

(IV.16)

1 − ep+ /1 1 − epu /1

(IV.17)

In view of that n → ∞ as t → ∞, we can see from (IV.14) and (IV.15) that there exists an inherent constraint for steady oscillation, ( − 0 )(1 − ePu /1 ) et0 /1 − 20 (1 − e−P+ /1 ) = 0

(IV.18)

which is similar to (III.26). Hence, in the limit cycle, it follows that y+ (t) = −kp ( + 0 ) + y− (t) = kp (0 − ) +

2kp 0 1 E1 t/ 2kp 0 2 E2 −t/ 2, e 1+ e 1 + 2 1 + 2

2kp 0 1 F1 t/ 2kp 0 2 F2 −t/ 2, e 1+ e 1 + 2 1 + 2

t ∈ [0, P+ ] t ∈ [0, P− ]

(IV.19) (IV.20)

where y+ (t) denotes the ascending response in the half period P+ , due to a positive step change of the relay output as shown in Fig. 5, while y− (t) denotes the descending response in the other half period P− , due to a negative step change of the relay output. It can be verified from (IV.19) and (IV.20) that the inherent connection between y+ (t) and y− (t) is satisfied, i.e., y+ (0) = y− (P− )

(IV.21)

y+ (P+ ) = y− (0)

(IV.22)

Here define z = x/ 1 , a =  1 / 2 , and f (x) =

1 − ex/1 , 1 − e−x/2

x ∈ (0, ∞)

(IV.23)

It follows that f (z) =

1 − ez 1 − e−az

(IV.24)

ez [(a + 1) e−az − a e−(a+1)z − 1] df (z) = dz (1 − e−az )2

(IV.25)

g(z) = (a + 1)e−az − a e−(a+1)z − 1

(IV.26)

Let

3056

T. Liu, F. Gao / Computers and Chemical Engineering 32 (2008) 3038–3056

It can be seen that dg(z) = a(a + 1) e−az (e−z − 1) < 0 dz

(IV.27)

lim g(z) = 0

(IV.28)

z→0

Thus, it can be concluded that g(z) < 0 for z ∈ (0,∞) and correspondingly, df (z) <0 dz

(IV.29)

Hence, f(z) decreases monotonically with respect to z and so does for f(x) with respect to x. We can then conclude from f (p− ) > f (pu ) that 0 < E1 < E2

(IV.30)

Following a similar analysis, it can be concluded that F2 < F1 < 0

(IV.31)

References Gu, D. Y., Ou, L. L., Wang, P., & Zhang, W. D. (2006). Relay feedback autotuning method for integrating processes with inverse response and time delay. Industrial and Engineering Chemistry Research, 45, 3119–3132. Ho, W. K., Feng, E. B., & Gan, O. P. (1996). A novel relay autotuning technique for processes with integration. Control Engineering Practice, 4(7), 923–928. ¨ Ingimundarson, A., & Hagglund, T. (2001). Robust tuning procedures of dead-time compensating controllers. Control Engineering Practice, 9, 1195–1208. Kaya, I. (2006). Parameter estimation for integrating processes using relay feedback control under static load disturbances. Industrial and Engineering Chemistry Research, 45, 4726–4731. Kaya, I., & Atherton, D. P. (2001). Parameter estimation from relay autotuning with asymmetric limit cycle data. Journal of Process Control, 11, 429–439. Kwak, H. J., Sung, S. W., & Lee, I. B. (1997). On-line process identification and autotuning for integrating processes. Industrial and Engineering Chemistry Research, 36, 5329–5338. Lee, J., & Sung, S. W. (1993). Comparison of two identification methods for PID controller tuning. AIChE Journal, 39(4), 695–697. Lin, C., Wang, Q. G., & Lee, T. H. (2004). Relay feedback: A complete analysis for first-order systems. Industrial and Engineering Chemistry Research, 43, 8400–8402. Liu, T., Cai, Y. Z., Gu, D. Y., & Zhang, W. D. (2005). New modified Smith predictor scheme for controlling integrating and unstable processes. IET Control Theory & Applications, 152(2), 238–246. Liu, T., Zhang, W. D., & Gu, D. Y. (2005). Analytical two-degree-of-freedom tuning design for open-loop unstable processes with time delay. Journal of Process Control, 15(5), 559–572. Luyben, W. L. (2003). Identification and tuning of integrating processes with deadtime and inverse response. Industrial and Engineering Chemistry Research, 42, 3030–3035. Majhi, S. (2007). Relay based identification of processes with time delay. Journal of Process Control, 17, 93–101. Majhi, M., & Atherton, D. P. (2000). Obtaining controller parameters for a new Smith predictor using autotuning. Automatica, 36, 1651–1658. Marchetti, G., Scali, C., & Lewin, D. R. (2001). Identification and control of open-loop unstable processes by relay methods. Automatica, 37, 2049–2055. Padhy, P. K., & Majhi, S. (2006). Relay based PI-PD design for stable and unstable FOPDT processes. Computers and Chemical Engineering, 30, 790–796. Paraskevopoulos, P. N., Pasgianos, G. D., & Arvanitis, K. G. (2004). New tuning and identification methods for unstable first order plus dead-time. IEEE Transactions on Control Systems Technology, 12(3), 455–464. Ramakrishnan, V., & Chidambaram, M. (2003). Estimation of a SOPTD transfer function model using a single asymmetrical relay feedback test. Computers and Chemical Engineering, 27, 1779–1784. Shen, S. H., Wu, J. S., & Yu, C. C. (1996). Use of biased-relay feedback for system identification. AIChE Journal, 42, 1174–1180. Shiu, S. J., Hwang, S. H., & Li, M. L. (1998). Automatic tuning of systems with one or two unstable poles. Chemical Engineering Communication, 167, 51–72. Sree, R. P., & Chidambaram, M. (2006). Improved closed loop identification of transfer function model for unstable systems. Journal of the Franklin Institute, 343, 152–160. Sung, S. W., & Lee, J. (2006). Relay feedback method under large static disturbances. Automatica, 42, 353–356. Tan, K. K., Wang, Q. G., & Lee, T. H. (1998). Finite spectrum assignment control of unstable time delay processes with relay tuning. Industrial and Engineering Chemistry Research, 37, 1351–1357. Thyagarajan, T., & Yu, C. C. (2003). Improved autotuning using the shape factor from feedback. Industrial and Engineering Chemistry Research, 42, 4425–4440. Tsang, K. M., Lo, W. L., & Rad, A. B. (2000). Autotuning of phase-lead controller for integrating systems. IEEE Transactions on Industrial Electronics, 47(1), 203–210. Vivek, S., & Chidambaram, M. (2005). An improved relay auto tuning of PID controllers for unstable FOPTD systems. Computers and Chemical Engineering, 29, 2060–2068. Wang, Q. G., Lee, T. H., & Lin, C. (2003). Relay feedback: Analysis, identification and control. London: Springer-Verlag. Yu, C. C. (2006). Autotuning of PID controllers: A relay feedback approach (2nd ed.). London: Springer-Verlag.