Relay Feedback Autotuning Method for Integrating Processes with

Mar 30, 2006 - Hence, the best tradeoff between robust stability and nominal performance of the closed-loop system can be achieved by adjusting the co...
16 downloads 13 Views 460KB Size
Ind. Eng. Chem. Res. 2006, 45, 3119-3132

3119

Relay Feedback Autotuning Method for Integrating Processes with Inverse Response and Time Delay Danying Gu,* Linlin Ou, Ping Wang, and Weidong Zhang Department of Automation, Shanghai Jiaotong UniVersity, Shanghai 200240, People’s Republic of China

Integrating processes with inverse response and time delay are frequently encountered in the level control of a boiler steam drum and are difficult to identify and control. In this paper, a systematic approach for autotuning of a proportional-integral (PI) or proportional-integral-derivative (PID) controller is proposed for this type of process. A single run of a biased relay feedback test is carried out to obtain four parameters of the model without any prior information about the time delay or the steady-state gain. Since the exact limit cycle expressions are derived, accurate parameter estimations are achieved. Then, according to the estimated model, an analytical procedure for PI/PID controller design is developed based on H∞ optimization and internal model control (IMC) theory. The internal stability of the closed-loop feedback system is guaranteed. Moreover, robust stability analysis for the proposed control scheme is provided in the presence of model mismatches. Hence, the best tradeoff between robust stability and nominal performance of the closed-loop system can be achieved by adjusting the controller parameter λ conveniently. For illustration, simulation examples are given to compare the proposed method with those previously used. 1. Introduction The processes with inverse response initially show a step response in the opposite direction to where it eventually ends up.1 The paradoxical and problematic characteristics of an inverse-response system result from two (or more) competing dynamic effects with different time constants from the same manipulated variable. This type of response can be observed in some practicle control loops in the chemical industry, such as level controls in the drum boiler and the reboiler section of a distillation column and also the exit temperature control of a tubular catalytic reactor via the reactor feed temperature.2 Inverse-response processes including integrator and time delay are also frequently encountered in the level control of a boiler steam drum. Actually, almost every manufacturing complex has a steam utility system. The heart of the steam system is the boiler. The steam drum level controller, which is one of the important control loops on a steam boiler, is illustrated in Figure 1.3 If there is an increase in the boiler feedwater (BFW) and no corresponding increase in outlet steam flow, the level will necessarily increase continually. Hence, the level control of a boiler steam drum possesses an integral character. The inverse response occurs because the cold BFW causes some of the bubbles in the steam drum to collapse and results in a slight drop of the liquid level. However, since the flow rate of BFW is higher than the outgoing flow rate, the level will eventually be driven up. In practice, the constraints on the steam drum level are important. If the liquid level grows too high, liquid will enter the superheater section above the steam drum, expanding rapidly and causing the “riser” pipe to rupture. If the liquid level gets too low, then water will be no longer in the “downcomer” pipes in the radiation section below the steam drum, causing the pipes to get too hot and break down.3 On the other hand, because of the combination of integration and inverse response, the control * To whom correspondence should be addressed. E-mail: [email protected]. Tel.: +86-21-34202021. Fax: +86-2154260762.

Figure 1. Control instrumentation diagram for boiler level.

of the steam drum level is tougher than most other level-control problems. Many papers1,2,4,5 have discussed proportional-integral (PI) and proportional-integral-derivative (PID) controller tuning methods for the inverse-response systems. Scali and Rachid4 provided a design procedure of a PID controller for an inverseresponse process without time delay based on internal model control (IMC) theory. This control scheme is able to give acceptable performance both for set-point tracking and for loaddisturbance rejection. A PI controller tuning method for processes with both inverse response and time delay was proposed by Luyben.5 Zhang et al.1 developed a modified Smith predictor for the inverse-response system based on modern H∞ control theory. Chien et al.2 derived PID and PI-only tuning methods from a direct-synthesis controller design method. However, there have been very few discussions on the integrating processes with inverse response and time delay. Luyben6 presented an identification method for this type of system from step-response data. On the basis of this, the PI and PID tuning methods were discussed in this excellent paper. For autotuning of the above PI/PID controllers, the process model transfer functions first have to be identified. The relay feedback method has become a widely accepted approach to system identification in process control. Many research works on modifying the relay feedback identification and autotuning methods for open-loop stable or unstable processes have been reported in recent years.7-9 However, only a few investigations have been carried out on relay feedback identification of inverse-

10.1021/ie050739n CCC: $33.50 © 2006 American Chemical Society Published on Web 03/30/2006

3120

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006

Figure 2. Block diagram of a relay feedback system.

response systems. It is rather difficult to identify the integrating processes with inverse response and time delay. On the other hand, since experiments for identification are usually expensive and subjected to unknown disturbances, it is most desirable to complete the identification within a few constant cycles in a single run. In the work presented herein, the biased relay feedback identification method is proposed under guidance of this requirement. Without any prior information about the time delay or the steady-state gain, four parameters of the integrating model with inverse response and time delay are obtained effectively by means of only a single relay feedback test. Since the exact analytical expressions of limit cycle are derived in terms of biased relay feedback measurements, more accurate parameter estimations are achieved. Once the process model is found, the PI/PID controllers are designed analytically based on H∞ optimization and IMC theory,10 so that the internal stability is guaranteed and the adjustable robustness is achieved.

Figure 3. Biased relay.

2. Identification with Biased Relay Feedback 2.1. Proposed Identification Procedure. Consider the following integrating process with inverse response and time delay

G(s) )

Kp(-τ1s + 1) (τ2s + 1)s

e-θs

[

]

2µ(1 - e-Pu2/τ2) Kp(µ-t1 + µ+t*) - Kp(τ1 + τ2) µ+ - e-t*/τ2 1 - e-(Pu1+Pu2)/τ2 (2)

[

y2n+1 ) Kp(µ-t1 + µ+Pu1 + µ-t*) - Kp(τ1 + τ2) µ- + e

-t*/τ2 2µ(1

the value of the process output y at t* ) t/min is the local minimum. This yields

(1)

where Kp is the steady-state gain, τ1 and τ2 are the time constants, and θ is the time delay. In this case, Kp, τ1, τ2, and θ are all positive. To obtain these four parameters from a single relay experiment and to reduce the influence of the measurement noise, a relay feedback system of Figure 2 is set up and the biased relay as shown in Figure 3 is introduced. The outputs of the biased relay are µ0 + µ and µ0 - µ, respectively, and the hysteresis of relay is . The resulting input u and output y of the process are shown in Figure 4. The following theorem gives the analytical expressions of this limit cycle oscillation. Theorem 1: For the process of eq 1 under the biased relay feedback of Figure 3, the process output y converges to the stationary oscillation in one period (Pu1 + Pu2), and the limit cycle oscillation is characterized by

y2n )

Figure 4. Input and output of a process under a biased relay feedback.

- e-Pu1/τ2)

1 - e-(Pu1+Pu2)/τ2

]

dy2n ) | dt* t*)tmin*

( )

K pµ + - K p

τ1 2µ(1 - e-Pu2/τ2) / + 1 e-tmin/τ2 ) 0 (4) τ2 1 - e-(Pu1+Pu2)/τ2

So t/min is given by

t/min ) -τ2 ln

(

µ+(1 - e-(Pu1+Pu2)/τ2)

2µ(τ1/τ2 + 1)(1 - e-Pu2/τ2)

)

) ˜t min - θ

(5)

and the corresponding minimum value of y is

ymin ) y2n(t*min) ) Kpµ-t1 Kpµ+τ2 ln

(

µ+(1 - e-(Pu1+Pu2)/τ2)

2µ(τ1/τ2 + 1)(1 - e-Pu2/τ2)

)

- Kpτ1µ+ ) Ad (6)

Similarly, the value of the process output y at t* ) t/max is the local maximum. This yields

dy2n+1 | ) K pµ - + dt* t* ) tmax* τ1 2µ(1 - e-Pu1/τ2) / Kp + 1 e-tmax/τ2 ) 0 (7) τ2 1 - e-(Pu1+Pu2)/τ2

( )

(3)

Proof: See the Appendix. Equations 2 and 3 are the exact expressions of the limit cycle oscillation for the process in eq 1. It is clear from Figure 4 that

So t/max is given by

t/max )

-τ2 ln

(

-µ-(1 - e-(Pu1+Pu2)/τ2)

2µ(τ1/τ2 + 1)(1 - e-Pu1/τ2)

)

) ˜t max - θ

(8)

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006 3121

and the corresponding maximum value of y is

or

ymax ) y2n+1(t*max) ) Kp(µ-t1 + µ+Pu1) Kpµ-τ2 ln

(

-µ-(1 - e-(Pu1+Pu2)/τ2)

2µ(τ1/τ2 + 1)(1 - e-Pu1/τ2)

)

θ ) ˜t max + τ2 ln

- Kpτ1µ- ) Ap (9)

Accordingly, by subtracting eq 8 from eq 5 and simplifying it, we get

-τ2 ln

(

µ+(1 - e-Pu1/τ2)

-µ-(1 - e-Pu2/τ2)

)

) ˜t min - ˜t max

(10)

Kp(µ+Pu1 + 2µτ1) -

[ ( (

µ+ ln

-µ-(1 - e-(Pu1+Pu2)/τ2)

2µ(τ1/τ2 + 1)(1 - e-Pu1/τ2)

µ+(1 - e-(Pu1+Pu2)/τ2)

)]

2µ(τ1/τ2 + 1)(1 - e-Pu2/τ2)

)

-

) Ap - Ad (11)

Figure 4 shows that when t* ) Pu1 - θ, y2n ) . Applying this boundary condition in eq 2, we have

y2n|t*)Pu1-θ ) Kp(µ-t1 + µ+(Pu1 - θ)) -

[

]

2µ(1 - e-Pu2/τ2) Kp(τ1 + τ2) µ+ - e-(Pu1-θ)/τ2 )  (12) 1 - e-(Pu1+Pu2)/τ2 Similarly, by substituting the boundary conditions t* ) Pu2 θ and y2n+1 ) - in eq 3, we have

y2n+1|t*)Pu2-θ ) Kp(µ-t1 + µ+Pu1 + µ-(Pu2 - θ)) -

[

]

2µ(1 - e-Pu1/τ2) Kp(τ1 + τ2) µ- + e-(Pu2-θ)/τ2 ) - (13) 1 - e-(Pu1+Pu2)/τ2 Accordingly, by subtracting eq 13 from eq 12 and simplifying it, we get

[

-Kp(µ-Pu2 + 2µθ) - 2Kpµ(τ1 + τ2) 1 e

θ/τ2

-(Pu1+Pu2)/τ2

1-e

]

(e-Pu1/τ2 + e-Pu2/τ2 - 2 e-(Pu1+Pu2)/τ2) ) 2 (14)

In the biased relay feedback test, µ, µ0, and  are known and Pu1, Pu2, Ap, Ad, ˜tmin, and ˜tmax can be directly tested from the limit cycle oscillation. So eq 10 contains only one unknown parameter τ2, and it can be solved with a simple iterative method such as Newton interpolation. Furthermore, the Matlab function “Fsolve” can easily be used to numerically find the solution of τ2. Thus, with τ2 known, the time delay θ is determined only by τ1 from eq 5 or 8 as

θ ) ˜t min + τ2 ln

[

]

2µ(τ1/τ2 + 1)(1 - e-Pu1/τ2)

(16)

Similarly, the steady-state gain Kp is determined only by τ1 from eq 11 as Kp )

[ (

(µ+ Pu1 + 2µτ1) - τ2 µ- ln

A p - Ad -µ-(1 - e-(Pu1+Pu2)/τ2)

2µ(τ1/τ2 + 1)(1 - e-Pu1/τ2)

) ( - µ+ ln

µ+(1 - e-(Pu1+Pu2)/τ2)

2µ(τ1/τ2 + 1)(1 - e-Pu2/τ2)

)]

(17)

Similarly, subtracting eq 6 from eq 9 yields that

Kpτ2 µ- ln

[

-µ-(1 - e-(Pu1+Pu2)/τ2)

µ+(1 - e-(Pu1+Pu2)/τ2)

2µ(τ1/τ2 + 1)(1 - e-Pu2/τ2)

]

(15)

Substituting eqs 17 and 15 (or 16) into eq 14, we can find the solution of τ1. The above derivation can be summarized as the following identification procedure. 2.1.a. Identification Procedure. The biased relay feedback experiment is performed. The relay output µ0 + µ and µ0 - µ and the relay hysteresis  are known. The process input u and output y are recorded. The periods (i.e., Pu1 and Pu2), the amplitudes (i.e., Ap and Ad), and the time intervals (i.e., ˜tmin and ˜tmax) are measured from the limit cycle oscillation. Step 1. Compute the time constant τ2 from eq 10. Step 2. Substitute eqs 17 and 15 (or 16) into eq 14. Then compute the time constant τ1. Step 3. Compute the time delay θ from eq 15 (or 16). Step 4. Compute the steady-state gain Kp from eq 17. The proposed identification method makes it possible to identify four unknown parameters of the integrating model with inverse response and time delay with only a single relay test. It does not require any prior information about the time delay or the steady-state gain. The exact analytical expressions are derived for the desired parameters in terms of simple limit cycle measurements. 2.2. Examples. For assessment of the accuracy of the proposed method, numerical examples are carried out for different integrating processes with inverse response and time delay, where the outputs of the biased relay are 1.2 and -0.8, respectively, and the hysteresis of relay is 0.2. The measured results of limit cycle and the parameter identification results are presented in Table 1. The results show that very small errors exist in the proposed relay feedback identification method. The reliability of the relay feedback identification approach is compared with that of the identification method based on step response proposed by Luyben6 and the least-squares curve fitting method provided in Matlab toolbox (see Table 2). In Luyben’s method, the time delay is directly determined from the step-response curve. Then, the other three parameters are obtained by solving three nonlinear equations with the help of Matlab. As we all know, the Matlab curve fitting function “Polyfit” can find the coefficients of a specified-order polynomial that best fits a given set of data in the least-squares sense. Hence, the curve fitting method in the relay feedback test contains the following procedure: First, the process output y in the interval [θ, t1 + θ] is best fitted as a fourth-order polynomial, which corresponds to eq 3. Second, the Matlab function “Taylor” returns the fourth-order Maclaurin polynomial approximation to the process output y. Finally, since the corresponding coefficients of the above two fourth-order polynomials should have the same value, four parameters of the model can be calculated in theory. To simplify the computations, it is assumed that the time delay θ is read exactly

3122

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006

Table 1. Parameter Estimation from Biased Relay Feedback Test measured results

identified results

Pu1

Pu2

Ap

Ad

˜tmax

˜tmin

Kp

θ

τ1

τ2

-0.4s + 1 -0.1s e (s + 1)s

2.8804

4.3206

0.7952

-0.5838

1.2959

0.9348

0.6000

0.1002

0.3999

1.0000

-s + 1 -s e (2s + 1)s

7.8116

11.7174

1.0298

-0.7450

3.6031

2.8271

0.2500

1.0002

0.9999

2.0000

case

process

1

0.6

2

0.25

Table 2. Comparison with Other Identification Methods case

process

new method

-0.4s + 1 (s + 1)s

1

0.6

2

0.25

-s + 1 -s e (2s + 1)s

Luyben’s method

Curve fitting method

0.6000

-0.3999s + 1 -0.1002s e (1.0000s + 1)s

0.8674

-0.4629s + 1 -0.1s e (0.9633s + 1)s

0.6

0.2500

-0.9999s + 1 -1.0002s e (2.0000s + 1)s

0.4308

-1.2708s + 1 -s e (1.8137s + 1)s

0.25

-0.5508s + 1 -0.1s e (1.3740s + 1)s -1.519s + 1 -s e (2.6330s + 1)s

Table 3. Model Parameter Estimation from Biased Relay Experiments for High-Order Processes case

process

model

3

-10s + 1 -s e (s + 1)5 s

1.0095

-9.1679s + 1 -4.9844s e (1.8497s + 1)s

4

-5s + 1 e-0.5s (s + 1)(s2 + s + 1)s

1.1079

-3.9072s + 1 -3.6144s e (0.5179s + 1)s

(-s + 1)s

1.0550

-0.5408s + 1 -3.4874s e (3.1416s + 1)s

5

(2s + 1) s 2

e-s

from the limit cycle of relay feedback test. Moreover, Theorem 1 indicates that the waveforms of the process input u and output y are periodic with (Pu1 + Pu2). They can be expanded into a Fourier series. The direct-current components of these periodic waves are extracted, so the steady-state gain Kp can be exactly calculated as the ratio of the integral of process y and u via the following formula:11

∫0P +P Kp ) G(0) ) P +P ∫0 u1

u2

y(t) dt

u1

u2

u(t) dt

(18)

Table 2 shows that the proposed method can give nearly the exact estimation of the process parameters. Furthermore, the primary advantage of the new method is that it allows the identification of four unknown parameters with only a single

relay test and requires no prior information. The Nyquist curves of the real processes and the corresponding models are shown in Figure 5. 2.2.a. Model-Mismatch Issue. Up to now, we have examined the proposed identification method for low-order integrating processes with inverse response and time delay with the same model structure. Here, the proposed identification procedure is tested against model mismatches. Three typical model structures are employed, including second-order systems and high-order systems. The identification results are listed in Table 3. The Nyquist curves of the real processes and the low-order models are shown in Figure 6, and they are very close to each other over a phase range of 0 to π. With these model-mismatch tests, the proposed method can also be used to approximate the highorder integrating processes with inverse response and time delay by low-order models.

Figure 5. Nyquist curves of real processes and corresponding identified models.

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006 3123

Figure 6. Nyquist curves of high-order processes and corresponding low-order models.

Figure 7. Input-output responses of biased relay experiments with noise.

2.2.b. Noise Issue. In a realistic environment, the major concerns for any identification method are measurement noise. Different identification methods are more or less sensitive to noise. As for the measurement noise in the relay test, Astro¨m and Ha¨gglund17 pointed out that a hysteresis in the relay is a simple way to reduce the influence of the measurement noise. For assessment of the accuracy of the proposed method, the biased relay tests are performed in Cases 1 and 2 when the

measured process output is corrupted by normally distributed random noise with a variance of 0.006. Here, hysteresis of 0.2 is used to prevent a severe fluctuation of the relay output. The input-output responses for the biased relay tests with measurement noise are illustrated in Figure 7. It is noted that 10 cycles of periodical stationary oscillations are used in the identification procedure, and Pu1, Pu2, Ap, Ad, ˜tmin, and ˜tmax are calculated as the mean value of those in the same data, respectively. The

3124

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006

Figure 8. Nyquist curves of real processes and corresponding models in biased relay tests with noise.

Under the nominal condition, the sensitivity and complementary sensitivity functions of the closed-loop system can be expressed respectively as

S(s) )

1 ) 1 - G(s)Q(s) 1 - G(s)C(s)

(20)

G(s)C(s) ) G(s)Q(s) 1 - G(s)C(s)

(21)

Figure 9. Unit feedback control system.

T(s) ) Table 4. Model Parameter Estimation from Biased Relay Experiments with Noise case

process

-0.4s + 1 -s e (s + 1)s

1

0.6

2

-s + 1 -s e 0.25 (2s + 1)s

If one denotes by H(s) the transfer function matrix from r and d to y and u, it follows that

model

0.7077

-0.4583s + 1 -0.0515s e (1.0666s + 1)

[ ]

[ ]

y(s) r(s) ) H(s) u(s) d(s)

-0.4296s + 1 -1.4251s e 0.3044 (2.1826s + 1)s

(22)

where resultant model parameters are presented in Table 4. The Nyquist curves of the real processes and the models are shown in Figure 8. With these antinoise tests, the proposed method can reject noise quite effectively.

3. PI/PID Controller Design The essential characteristic of the integrating process with inverse response and time delay is that the process transfer function has both right-half-plane (RHP) zero and integrator. Generally, these processes are particularly difficult to control. In this section, H∞ optimization and IMC theory10 are used to design the PI/PID controller that guarantees the stability and performance of the closed-loop system. 3.1. Internal Stability of the Closed-Loop System. Consider the unit feedback control loop shown in Figure 9, where G(s) is the process given by eq 1 and C(s) is the PI/PID controller. The unit feedback loop can be equivalent to an IMC structure through10

C(s) )

Q(s) 1 - G(s)Q(s)

Here Q(s) is actually an IMC controller.

(19)

H(s) )

[

] [

]

G(s)Q(s) G(s)(1 - G(s)Q(s)) T(s) G(s)S(s) ) Q(s) -G(s)Q(s) Q(s) -T(s) (23)

The closed-loop system is regarded to be internally stable if all transfer functions in H(s) are stable. Because G(s) has an integral term 1/s, the sufficient and necessary conditions of internal stability are (1) Q(s) is stable and (2) G(s)(1 - G(s)Q(s)) is stable. The above conditions are equivalent to that Q(s) is stable and satisfies the following constraints:

lim(1 - T(s)) ) 0

(24)

d lim (1 - T(s)) ) 0 sf0 ds

(25)

sf0

3.2. PI/PID Controller Design Procedure. The primary objective of any feedback control scheme is to keep the controlled output y close to the desired set point r or, equivalently, to keep the complementary sensitivity function T close to 1. What is meant by “close” can be defined in terms of various performance specifications for the closed-loop system. Here, the H∞ optimal performance objective min |W(s)(1 G(s)Q(s))|∞ is utilized to derive the PI/PID controller, where W(s) is the weighting function and can be chosen as 1/s for the

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006 3125

step change of the set-point inputs. Then the optimal IMC controller Qim(s) can be designed as10

Qim(s) )

(τ2s + 1)s Kp

(b1s + b0)(τ2s + 1)s

(b1s + b0)(-τ1s + 1) (λs + 1)

3

e-θs

KC ) f′(0), TI )

(33)

f′(0) f(0)

(34)

If the process signal is not too noisy, the application of derivative action can improve the closed-loop dynamics in this system because the derivative action permits the use of smaller integral times. Therefore, using the first three terms of eq 32, the standard PID controller is derived by

(

C(s) ) KC 1 +

1 + TDs TIs

)

(35)

where

KC ) f′(0), TI )

(28)

The constraints, eqs 24 and 25, mean that the desired sensitivity function possesses two zeros at s ) 0 so that the step load disturbances injected into the integrating process input can be counteracted. Substituting eq 28 into eqs 24 and 25 yields b0 ) 1 and b1 ) 3λ + τ1 + θ. Therefore, the desired closed-loop transfer function is given by

)

1 TIs

where

(27)

(λs + 1)3Kp

The corresponding complementary sensitivity function of the closed-loop system is

T(s) )

(

C(s) ) KC 1 +

(26)

Obviously Qim(s) is not proper and cannot be physically realized in practice. Hence, the low-pass filter has to be introduced to roll Qim(s) off at high frequency. To simplify the design task, typically we use a single parameter filter of the form J(s) ) (b1s + b0)/(λs + 1)3, which makes Q(s) proper and satisfies the internal stability constraints, eqs 24 and 25. Here, λ > 0 is an adjustable parameter and represents the tradeoff between the disturbance rejection and robustness. Decreasing λ will result in better disturbance rejection and worse robustness, while increasing λ will result in worse disturbance rejection and better robustness. Then we have

Q(s) )

Obviously, the first two terms of the above expansion are exactly a standard PI controller in the form of

f′(0) f′′(0) , and TD ) f(0) 2f′(0)

(36)

To simply the calculation, let

N(s) M(s)

(37)

(τ2s + 1)[(3λ + τ1 + θ)s + 1] Kp

(38)

f(s) ) where

T(s) )

[(3λ + τ1 + θ)s + 1](-τ1s + 1) (λs + 1)3

e-θs

(29) N(s) )

Correspondingly, employing eq 21, the controller C(s) in the following form is obtained.

[(3λ + τ1 + θ)s + 1](τ2s + 1)s 1 C(s) ) 3 Kp (λs + 1) - [(3λ + τ + θ)s + 1](-τ s + 1) e-θs 1 1 (30) Since the internal stability constraints result in

it is clear that zero-pole canceling occurs at eq 30, which cannot be removed directly and will cause the control system to work unstably in practice. Therefore, the mathematical Maclaurin expansion formula is employed to reproduce the controller in a simple way.12,13 Accordingly, C(s) ) f(s)/(s), then C(s) is expanded in Maclaurin series as

[

f′′(0) 2 f′′′(0) 3 1 s + s ... f(0) + f′(0)s + s 2! 3!

{

]

(32)

In the level-control system of a boiler steam drum, the use of a PI controller is more appropriate than the use of a PID controller because the level signal is typically quite noisy.6

(λs + 1)3 - [(3λ + τ1 + θ)s + 1](-τ1s + 1) e-θs s2

(39)

Then we have f(0) )

lim{(λs + 1)3 - [(3λ + τ1 + θ)s + 1](-τ1s + 1) e-θs} ) 0 sf0 (31)

C(s) )

M(s) )

f′(0) )

N(0) M(0) N′(0)M(0) - M′(0)N(0)

f′′(0) )

M(0)2 N′′(0)M(0)2 - M′′(0)N(0)M(0) - 2M′(0)N′(0)M(0) + 2M′(0)2N(0) M(0)3

(40)

Since all the parameters are known, its computation is, in fact, very simple. To implement the PID controller physically, a lag 1/(1 + sTd/10) is usually introduced to the derivative term. 3.3. Stability and Performance Analysis. There always exist model mismatches between the “true” process transfer function and the “fitted” model transfer function. An important merit of the proposed controller design method is that the robust analysis is very simple. Suppose the process represented by eq 1 has uncertainty on four parameters: K ˜ p ) Kp(1 + ∆Kp), τ˜ 1 ) τ1(1 + ∆τ1), τ˜ 2 ) τ2(1 + ∆τ2), and θ˜ ) θ(1 + ∆θ), where ∆i

3126

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006

Figure 11. Block diagram for autotuning of a PI/PID controller.

Furthermore, recall the nominal performance condition of the feedback system which is internally stable. Then the sufficient and necessary condition that guarantees the nominal performance and the robust stability of the closed-loop system is10,14

||T(jω)∆m(ω)| + |W(jω)S(jω)||∞ < 1, ∀ω Figure 10. Multiplicative uncertainty profile ∆m(ω).

is the multiplicative uncertainty on each parameter, with the possibility of assuming negative or positive values. Define

(

(

) )

Gp(jω) - Gm(jω) Kp(1 + ∆Kp) |)| Kp Gm(jω) τ2ωj + 1 -τ1(1 + ∆τ1)ωj + 1 -j∆θω e - 1| -τ1ωj + 1 τ2(1 + ∆τ2)ωj + 1 (41)

∆m(ω) g |

)(

where ∆m(ω) is the multiplicative uncertainty profile, Gp(jω) is the real process, and Gm(jω) is the process model. For a quantitative evaluation of ∆m(ω), let us assume equal maximum uncertainty on all the system parameters: ∆Kp ) ∆τ1 ) ∆τ2 ) ∆θ ) δ. A typical trend of ∆m(ω) is illustrated in Figure 10, for increasing values of parametric uncertainty δ. Especially when only the steady-state gain Kp is uncertain, the expression is simplified to

∆Kp | ∆m(ω) ) | Kp

(42)

When the time constants τ1 and τ2 are uncertain, the expression is simplified to

∆m(ω) g |

(

τ2ωj + 1

)(

)

-τ1(1 + ∆τ1)ωj + 1 - 1| -τ1ωj + 1 τ2(1 + ∆τ2)ωj + 1 (43)

When only the time delay θ is uncertain, the expression is simplified to

{

-j∆θω - 1| ω max ∆θ < π ∆m(ω) ) |e 2 ω max ∆θ g π

(44)

The sufficient and necessary condition that guarantees the robust stability of the closed-loop system is10,14

|T(jω)∆m(ω)|∞ < 1, ∀ω

(45)

Substituting eq 29 into eq 45, we obtain the closed-loop robust stability constraint as

(λ2ω2 + 1)x(λ2ω2 + 1)

x[(3λ + τ1 + θ)τ1ω2 + 1]2 + (3λ + θ)2ω2

< ∆m-1(ω), ∀ω (46)

(47)

where W(s) is a weight function of the sensitivity function S(s) ) 1 - T(s) of the closed-loop system, which usually can be selected as 1/s. As we can see, the adjustable parameter λ aims at the tradeoff between the nominal performance and the robust stability of the closed-loop system. When λ increases, S(s) increases and T(s) decreases. The system can endure larger uncertainty with a worse nominal performance. When λ decreases, the system tends to have optimal nominal performance with a worse robust stability. The relationship is monotonic. Since there is a single adjustable parameter λ in both C(s) and its Maclaurin approximation, the practically proposed PI/PID controller can be conveniently tuned by λ to obtain the required robustness. Although none of the above constraints of robust performance can be solved analytically, the quantitative relationship between λ and the performance of the closed-loop system can be obtained by numerical simulation tests based on the MATLAB toolbox. For the proposed PID controller, the quantitative relationship is similar to that of conventional H∞ or H2 PID controller. Further research and illustration will be given in the simulation examples. 3.4. Autotuning Procedure. The block diagram for autotuning of the PI/PID controller is shown in Figure 11. The autotuning procedure to find controller parameters can be performed as follows: Step 1. Switch from the controller mode to relay mode, so that the biased relay feedback configuration is obtained. Step 2. Carry out the identification procedure proposed in Section 2.1 to estimate parameters for the model of the integrating process with inverse response and time delay. Step 3. Fix λ. Generally, it is recommended to fix λ around the process time delay in the first place. Step 4. Calculate tuning parameters using eqs 33-40. Step 5. Switch from relay mode to controller mode with calculated parameters for the PI/PID controller. Step 6. If the resultant closed-loop performance or robustness is not acceptable, monotonically increase or decrease λ on-line until the desirable tradeoff between the closed-loop nominal performance and robust stability is achieved. 4. Simulation Examples 4.1. Example 1. Consider the following integrating process with inverse response and time delay from Luyben:6

G(s) ) 0.547

-0.418s + 1 -0.1s e (1.06s + 1)s

The proposed PI/PID method is compared to some other tuning

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006 3127 Table 5. Various Tuning Rules of PI/PID Controller for Example 1 PI

PID

tuning method

KC

TI

ISEa

ZN RZN AH Luyben proposed

b)1

1.3978 0.9142

4.5441 7.7682

17.9244 14.7404

λ ) 2.0

0.854 0.8438

23 7.0759

16.1512 14.9358

a

b ) 0.7965 R ) 4, Φm ) π/4 λ ) 0.8

KC

TI

TD

1.8637 1.8637 2.1964 1.61 2.0883

2.6730 2.6730 1.0271 5.75 3.8664

0.6683 0.6683 4.1082 1.15 0.6879

ISE evaluated on the PI control system with process noise.

Figure 12. Responses of nominal system with PI control for Example 1.

Figure 14. Responses of PI control system with process noise for Example 1.

Figure 13. Responses of nominal system with PID control for Example 1.

methods: Ziegler-Nichols (ZN) method,15 refined ZieglerNichols (RZN) method,16 A° stro¨m-Ha¨gglund (AH) method,17 and Luyben’s method.6 The process under ideal relay with unit height control provides a limit cycle frequency of 1.1753 with an amplitude of 0.4099. As a result, the ultimate gain is 3.1062 and the ultimate period is 5.3460, which are necessary for ZN, RZN, and AH methods. In Luyben’s PID tuning method, the integral time is 25% of the corresponding PI integral time so that the load response is improved. In our proposed method, take λ ) 20θ, i.e., λ ) 2.0 for PI controller and λ ) 0.8 for PID controller in order to obtain the similar set-point rising speed as those of the other controllers for comparison. By employing the design formulas eqs 33-40, we can obtain the parameters of the PI/PID controller. All the tuning rules of the PI/PID controller are illustrated in Table 5. In the closed-loop control system, a unit step input is added at t ) 0s and an inverse step load disturbance with magnitude

Figure 15. Complementary sensitivity function T(jω) satisfying robustness bound.

0.5 is added to the process input at t ) 40s. The responses of the nominal system with PI/PID control are shown in Figures 12 and 13. Figure 12 shows that the proposed PI method gives a similar result to that of the RZN method and is apparently superior to the ZN method. Moreover, the new method is better than the Luyben method in terms of load-disturbance rejection, but with a slightly larger overshoot. Figure 13 shows that the proposed PID controller is better than those of the Luyben and RZN methods for load-disturbance rejection and is obviously superior to that of the ZN method for set-point response. Compared with the AH and Luyben methods in terms of set-point response, the new controller provides a shorter settling time, but with a slightly larger overshoot.

3128

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006

Figure 16. Responses of perturbed system with PI control for Example 1 (∆ ) 10%).

Figure 17. Responses of perturbed system with PID control for Example 1 (∆ ) 10%).

Figure 19. Relationship between overshoot and λ.

Figure 20. Relationship between perturbance peak and λ. Table 6. Parameters of PID Controller with Adjusting λ adjustable parameter

parameters of PID controller

λ

KC

TI

TD

0.5 0.8 1.5 2

3.1986 2.0883 1.1281 0.8438

3.0592 3.8664 5.7402 7.0759

0.6796 0.6879 0.5886 0.4742

gives a response closer to that of the RZN method and is more robust for process noise than the ZN and Luyben methods. The robustness is more important than the nominal performance in practice. One key in designing the controller is to achieve the tradeoff between the robust stability and the nominal performance. Now consider the model-mismatch issues. Suppose that the process has equal maximum uncertainty on gain, time delay, and time constants: ∆Kp ) ∆τ1 ) ∆τ2 ) ∆θ ) 10%. In this case, the multiplicative uncertainty profile can be calculated from eq 41 to be Figure 18. Responses of perturbed system with increasing λ.

With the existence of process noise, the use of a PI controller is more popular than the use of a PID controller. Here, the random noise which is normally distributed with a variance of 0.005 is introduced to the PI control system to illustrate the effectiveness of the proposed PI tuning method and some other PI methods. Figure 14 shows the closed-loop responses for the PI control system with a noisy process signal. The PI controller settings and the closed-loop integral-squared-error (ISE) values are given in Table 5. This indicates that the proposed PI method

∆m(ω) g |1.1

1.06ωj + 1 -0.46ωj + 1 e (1.166ωj + 1)(-0.418ωj + 1)

-j0.01ω

- 1|

According to eq 46, the closed-loop system with the proposed PI/PID controller has robust stability if and only if

(λ2ω2 + 1)x(λ2ω2 + 1)

x[(3λ + 0.518)0.418ω

2

+ 1] + (3λ + 0.1) ω 2

2

2

< ∆m-1(ω), ∀ω

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006 3129

stability of the closed-loop system. For a more detailed study, the closed-loop responses of the perturbed system with PI/PID control are shown in Figures 16 and 17. Figures 16 and 17 show that the proposed method maintains robust stability in terms of set-point response and loaddisturbance rejection. It should be noted that, by increasing λ in the proposed PID controller, the robust stability of the closedloop system could become better at the cost of degraded nominal performance. To illustrate this, we carried out simulation tests with the perturbed process, adjusting λ for four cases listed with the corresponding PID controller forms in Table 6. The simulation results of the perturbed system responses are shown in Figure 18.

Figure 21. Relationship between recovery time and λ.

It can be clearly seen from Figure 15 that the magnitude of the complementary sensitivity function |T(jω)| always satisfies the robustness bound ∆m-1(ω) for increasing the controller parameter λ. This means that the proposed control method can guarantee the robust stability of the closed-loop system. Furthermore, larger λ values decrease the magnitude of the complementary sensitivity function and, thus, improve the robust

Figure 22. Responses of nominal system with PI control for Example 2.

Figure 18 shows that adjusting λ can achieve the best tradeoff between robust stability and nominal performance of the closedloop system. Further research shows the quantitative relationship between λ and the closed-loop performance. The perturbance peak is defined as the peak of system output when a step disturbance is added to the process input, and recovery time is defined as the perturbed system output coming to its final value. The quantitative performance is shown in Figures 19-21. Hence, the time-domain performance of the closed-loop system, such as overshoot, perturbance peak, and recovery time, can be estimated quantitatively. This is very attractive in a practical autotuning procedure.

3130

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006

Figure 23. Responses of nominal system with PID control for Example 2. Table 7. Various Tuning Rules of PI/PID Controller for Example 2 proposed PI

ZN PI

proposed PID

ZN PID

case

τ1

θ

λ

KC

TI

KC

TI

λ

KC

TI

TD

KC

TI

TD

A B C D

0.2 1.6 0.2 1.6

0.2 0.2 1.6 1.6

1.20 2.80 2.88 4.16

0.810 0.259 0.267 0.168

4.745 10.754 11.247 16.532

1.058 0.210 0.275 0.138

3.788 9.438 9.357 14.263

0.80 1.60 1.76 2.88

1.236 0.390 0.407 0.223

3.683 7.502 8.327 13.160

0.642 0.780 1.056 1.180

1.411 0.280 0.367 0.184

2.229 5.552 5.504 8.390

0.557 1.388 1.376 2.098

4.2. Example 2. For further assessment of the closed-loop performance of the proposed PI/PID controller, the following process with different parameters is considered:

G(s) )

(-τ1s + 1) e-θs (s + 1)s

Two values of the time delay (θ ) 0.2 and 1.6) and also two values of the τ1 parameter (τ1 ) 0.2 and 1.6) are tested, so there are four subexamples. Because the method by Luyben6 is based on an empirical relationship, this process with different values of θ and τ1 requires a different empirical tuning method. The advantage of the proposed method is that it is derived from the model parameters; thus, it is more direct and convenient than the Luyben method. Here, the proposed method is compared to the ZN PI and PID tuning rules. Table 7 gives the PI/PID tuning rules used in the simulation. The resultant responses of the nominal system with PI/PID control are shown in Figures 22 and 23, respectively. Notice that the ZN PI/PID settings give

quite oscillatory responses for all four different cases. On the other hand, the proposed method gives very satisfactory closedloop responses despite the changes in θ and τ1. To test our method in a realistic environment, the normally distributed random noise with a variance of 0.005 is also introduced to the PI control system in Case A (θ ) 0.2 and τ1 ) 0.2). Figure 24 shows the closed-loop responses for the PI control system with process noise. The closed-loop ISE values of the proposed method and the ZN method are 10.7857 and 15.2725, respectively. This clearly indicates that the proposed method is more effective to overcome process noise. A simulation run is also made to demonstrate the closedloop performance of the proposed tuning method under the model-mismatch condition. Assuming the PID tuning constants are calculated based on the model parameters in Case D (θ ) 1.6 and τ1 ) 1.6), the actual process has 10% model mismatches in gain, time delay, and time constants, i.e., ∆Kp ) ∆τ1 ) ∆τ2 ) ∆θ ) 10%. Closed-loop responses in Figure 25 verify that

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006 3131

Phosphor Program of Shanghai (04QMH1405), National Key Fundamental Research Program (2002cb312200), and NCET (04-0383). Appendix: Proof of Theorem 1 We derive the analytical expressions of the limit cycle by using the time-domain analysis.18 The transfer function of a typical integrating process with inverse response and time delay is given in eq 1. The process input in the relay feedback test consists of a series of step changes with down-amplitude µ0 µ and up-amplitude µ0 + µ. Suppose the relay switches to u ) µ0 - µ at t ) 0. After a time delay θ, the relay feedback response is shown in Figure 4. Let µ- ) µ0 - µ, µ+ ) µ0 + µ. Thus, the system output is

Figure 24. Responses of PI control system in Case A with process noise.

Y(s) ) U(s)Kp

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

e-θs )

µ- (-τ1s + 1) -θs K (A1) e s p (τ2s + 1)s

This equation can be rearranged as

[

Y(s) ) Kpµ-

]

1/τ2 1 - (τ1 + τ2) e-θs 2 s(s + 1/τ2) s

(A2)

For convenience, we shift the coordinate origin to t ) θ. With an inverse Laplace transform, the time-domain expressions of y for different time intervals can be written as follows:

y1(t) ) Kpµ-[t - (τ1 + τ2)(1 - e-t/τ2)] 0 e t < t1 (A3) y2(t) ) y1(t) + 2Kpµ[(t - t1) - (τ1 + τ2)(1 - e-(t-t1)/τ2)] ) Kpµ-[t - (τ1 + τ2)(1 - e-t/τ2)] + 2Kpµ[(t - t1) Figure 25. Responses of perturbed system with PID control in Case D (∆ ) 10%).

the proposed method maintains better robustness in comparison to the ZN method.

(τ1 + τ2)(1 - e-(t-t1)/τ2)] t1 e t < t1 + Pu1 (A4) y3(t) ) y2(t) - 2Kpµ[(t - t1 - Pu1) - (τ1 + τ2)(1 e-(t-t1-Pu1)/τ2)] t1 + Pu1 e t < t1 + Pu1 + Pu2 (A5)

5. Conclusions This paper considers autotuning of PI/PID controller for integrating processes with inverse response and time delay. The identification method based on a biased relay feedback test makes it possible to identify four parameters of the model with only a single relay test. It does not require any prior information about the time delay or the steady-state gain. The existence of both right-half-plane zero and integrator makes it difficult to obtain a higher control quality. This paper not only improves the identification method but also gives an analytical design procedure for a PI/PID controller based on modern H∞ control theory and IMC theory. The most important feature of the proposed method is that the robust stability and nominal performance of the closed-loop system can be conveniently adjusted by one controller parameter λ. Moreover, the quantitative relationship between λ and the closed-loop performance is researched. For a known nominal process, the time-domain performance of the closed-loop system can be estimated quantitatively, which is very attractive in a practical autotuning procedure. Acknowledgment This work was supported by the National Natural Science Foundation of China (60474031), Science and Technology

y4(t) ) y3(t) + 2Kpµ[(t - t1 - Pu1 - Pu2) (τ1 + τ2)(1 - e-(t-t1-Pu1-Pu2)/τ2)] t1 + Pu1 + Pu2 e t < t1 + 2Pu1 + Pu2 (A6) y5(t) ) y4(t) - 2Kpµ[(t - t1 - 2Pu1 - Pu2) - (τ1 + τ2)(1 e-(t-t1-2Pu1-Pu2)/τ2)] t1 + 2Pu1 + Pu2 e t < t1 + 2Pu1 + 2Pu2 (A7) Let t* ) t - t1. Equation A4 can be rewritten as

y2(t* + t1) ) Kpµ- [(t* + t1) - (τ1 + τ2)(1 - e-(t*+t1)/τ2)] + 2Kpµ[t* - (τ1 + τ2)(1 - e-t*/τ2)] (A8) Equation A8 can be further simplified as follows:

y2(t* + t1) ) Kp(µ-t1 + µ+t*) - Kp(τ1 + τ2)[µ+ e-t*/τ2(2µ + µ- e- t1/τ2)] (A9) Similarly, when t* ) t - t1 - Pu1, eq A5, after simplification, can be rewritten as

3132

Ind. Eng. Chem. Res., Vol. 45, No. 9, 2006

y3(t* + t1 + Pu1) ) Kp(µ-t1 + µ+(t* + Pu1)) -

lim ((1 - x1)(1 + x1x2 + x12x22 + ‚‚‚ + x1n-1x2n-1) -

nf∞

Kp(τ1 + τ2)[µ+ - e-(t*+Pu1)/τ2 (2µ + µ- e-t1/τ2)] 2Kpµ[t* - (τ1 + τ2)(1 - e

x1nx2n-1(µ - µ0) e-t1/τ2) ) (1 - x1)

-t*/τ2

)] ) Kp(µ-t1 + µ+Pu1 +

µ-t*) - Kp(τ1 + τ2)[µ- + e-t*/τ2 (2µ - e-Pu1/τ2 (2µ + -t1/τ2

µ- e

lim x1nx2n-1(µ - µ0) e-t1/τ2 )

))] (A10)

nf∞



x1n-1x2n-1 ∑ n)1

1 - x1

1 - x1x2

When t* ) t - t1 - Pu1 - Pu2, eq A6, after simplification, can be rewritten as

)

1 - e-Pu1/τ2 1 - e-(Pu1+Pu2)/τ2 (A19)

Thus, it is clear that

y4(t* + t1 + Pu1 + Pu2) ) Kp(µ-t1 + µ+Pu1 + µ-(t* + Pu2)) - Kp(τ1 + τ2)[µ- + e- (t*+Pu2)/τ2 (2µ - e-Pu1/τ2 (2µ + µ- e-t1/τ2))] + 2Kpµ[t* - (τ1 + τ2)(1 - e-t*/τ2)] ) Kp(µ-t1 + µ+Pu1 + µ-Pu2 + µ+t*) - Kp(τ1 + τ2)[µ+ e-t*/τ2(2µ - e-Pu2/τ2 (2µ - e-Pu1/τ2 (2µ + µ- e-t1/τ2)))] (A11)

(A12)

µ+Pu1 + µ-Pu2 ) 0

(A13)

Hence, eq A11 can be further simplified as

y4(t* + t1 + Pu1 + Pu2) ) Kp(µ-t1 + µ+t*) Kp(τ1 + τ2)[µ+ - e-t*/τ2(2µ - e-Pu2/τ2(2µ - e-Pu1/τ2(2µ + µ- e-t1/τ2)))] (A14) Similarly, when t* ) t - t1 - 2Pu1 - Pu2, eq A6, after simplification, can be rewritten as

y5(t* + t1 + 2Pu1 + Pu2) ) Kp(µ-t1 + µ+Pu1 + µ-t*) Kp(τ1 + τ2)[µ- + e-t*/τ2(2µ - e-Pu1/τ2(2µ - e-Pu2/τ2 (2µ e-Pu1/τ2(2µ + µ- e-t1/τ2))))] (A15) Let x1 ) e-Pu1/τ2 and x2 ) e-Pu2/τ2. Thus, the generalized expressions can be written as

y2n ) Kp(µ-t1 + µ+t*) - Kp(τ1 + τ2)[µ+ - e-t*/τ2(2µ((1 + x1x2 + x12x22 + ‚‚‚ + x1nx2n) -x2(1 + x1x2 + x12x22 + ‚‚‚ + x1n-1x2n-1)) + x1nx2nµ- e-t1/τ2)] (A16) y2n+1 ) Kp(µ-t1 + µ+Pu1 + µ-t*) - Kp(τ1 + τ2)[µ- + e-t*/τ2(2µ(1 - x1)(1 + x1x2 + x12x22 + ‚‚‚ + x1n-1x2n-1) + x1nx2n-1µ- e-t1/τ2)] (A17) As n f ∞, the infinite series in eqs A16 and A17 become

lim ((1 + x1x2 + x12x22 + ‚‚‚ + x1nx2n) - x2(1 + x1x2 +

nf∞

x12x22 + ‚‚‚ + x1n-1x2n-1) - x1nx2n(µ - µ0) e-t1/τ2) )

∑ n)0

]

1 - e-(Pu1+Pu2)/τ2 2µ(1 - e-Pu1/τ2) y2n+1 ) Kp(µ-t1 + µ+Pu1 + µ-t*) - Kp(τ1 + τ2) µ- + e-t*/τ2 1 - e-(Pu1+Pu2)/τ2

[

]

(A20)

Literature Cited

Pu1 µ| |) µ+ Pu2 i.e.,



x1nx2n - x2

2µ(1 - e-Pu2/τ2)

The proof is completed.

For the limit cycle, it is satisfied that



{

[

y2n ) Kp(µ-t1 + µ+t*) - Kp(τ1 + τ2) µ+ - e-t*/τ2

x1n-1x2n-1 - lim x1nx2n(µ - µ0) e-t /τ ) ∑ nf∞ n)1

(1) Zhang, W. D.; Xu, X. M.; Sun, Y. X. Quantitative Performance Design for Inverse-Response Processes. Ind. Eng. Chem. Res. 2000, 39, 2056-2061. (2) Chien, I. L.; Chung, Y. C.; Chen, B. S.; Chuang, C. Y. Simple PID Controller Tuning Method for Processes with Inverse Response Plus Dead Time or Large Overshoot Response Plus Dead Time. Ind. Eng. Chem. Res. 2003, 42, 4461-4477. (3) Bequette, B. W. Process Control: Modeling, Design and Simulation; Prentice Hall: Upper Saddle River, NJ, 2003. (4) Scali, C.; Rachid, A. Analytical Design of Proportional-IntegralDerivation Controllers for Inverse Response Processes. Ind. Eng. Chem. Res. 1998, 37, 1372-1379. (5) Luyben, W. L. Tuning of Proportional-Integral Controllers for Processes with Both Inverse Response and Deadtime. Ind. Eng. Chem. Res. 2000, 39, 973-976. (6) Luyben, W. L. Identification and Tuning of Integrating Processes with Deadtime and Inverse Response. Ind. Eng. Chem. Res. 2003, 42, 30303035. (7) Luyben, W. L. Processing Modeling, Simulation, and Control for Chemical Engineers, 2nd ed.; McGraw-Hill: New York, 1990. (8) Yu, C. C. Autotuning of PID Controllers: Relay Feedback Approach; Springer-Verlag: London, 1999. (9) Wang, Q. G.; Lee, T. H.; Lin, C. Relay Feedback: Analysis, Identification and Control; Springer-Verlag: London, 2003. (10) Morari, M.; Zafiriou, E. Robust Process Control; Prentice Hall: Englewood Cliffs, NJ, 1989. (11) Ramirez, R. W. The FFT Fundamentals and Concepts; Prentice Hall: Englewood Cliffs, NJ, 1985. (12) Lee, Y.; Lee, M.; Park, S.; Brisilow, C. PID Controller Tuning for Desired Closed Loop Response for SISO System. AIChE J. 1998, 44, 106155. (13) Zhang, W. D.; Xu, X. M.; Sun, Y. X. Quantitative Performance Design for Integrating Processes with Time Delay. Automatica 1999, 35, 719-723. (14) Doyle, J. C.; Francis, B. A.; Tannenbaum, A. R. Feedback Control Theory; Macmillan Publishing Company: New York, 1992. (15) Ziegler, J. G.; Nichols, N. B. Optimum Settings for Automatic Controllers. Trans. ASME 1942, 64, 759-768. (16) Hang, C. C.; A° stro¨m, K. J.; Ho, W. K. Refinements of the ZieglerNichols Tuning Formula. IEE Proc. Control Theory Appl. 1991, 138 (2), 111-118. (17) A° stro¨m, K. J.; Ha¨gglund, T. Automatic Tuning of Simple Regulators with Specifications on Phase and Amplitude Margins. Automatica 1984, 20 (5), 645-651. (18) Thyagarajan, T.; Yu, C. C.; Huang, H. P. Assessment of Control Performance: A Relay Feedback Approach. Chem. Eng. Sci. 2003, 58, 497512.

1 2

1 - x2

1 - x1x2

)

1 - e-Pu2/τ2 1 - e-(Pu1+Pu2)/τ2

(A18)

ReceiVed for reView June 21, 2005 ReVised manuscript receiVed February 23, 2006 Accepted February 24, 2006 IE050739N