A Methodology for Autotuning of Multivariable Systems - Industrial

Jul 20, 2002 - (Sequential Identification and Autotuning by Relay Techniques of Decentralised Controllers for MIMO Processes, ADCHEM 2000, Pisa, Italy...
0 downloads 0 Views 242KB Size
Download PDF
Ind. Eng. Chem. Res. 2002, 41, 4605-4615

4605

A Methodology for Autotuning of Multivariable Systems W. K. Toh and G. P. Rangaiah* Department of Chemical and Environmental Engineering, National University of Singapore, Singapore 119260, Singapore

A methodology for autotuning decentralized proportional-integral-derivative (PID) controllers for multivariable systems is proposed. It combines the autotuning procedure involving sequential loop closing and relay tests with time-domain curve fitting via least squares to identify suitable transfer function models that account for interactions among the loops, using limited response data. Simulation results reveal that the curve fitting is able to identify models using only the first few oscillations of the relay-test response. Thus, there is no need to wait for the process to reach steady oscillations. A controller tuning method developed by Semino and Scali (J. Process Control 1998, 8 (3), 219-227) is adopted to tune the PI controllers. Detailed results on several 2 × 2, 3 × 3, and 4 × 4 multivariable processes show that the proposed methodology gives PI controllers having performances comparable with those tuned by the ATV+ method of Parabita et al. (Sequential Identification and Autotuning by Relay Techniques of Decentralised Controllers for MIMO Processes, ADCHEM 2000, Pisa, Italy, 2000) or Luyben’s BLT technique (Ind. Eng. Chem. Process Des. Dev. 1986, 25 (3), 654-660). The major advantage of the proposed autotuning with least-squares is the significant reduction in the duration of relay tests. Computations involved in the methodology can be carried out on a personal computer in reasonable time. 1. Introduction Decentralized proportional-integral-derivative (PID) controllers are used to control multivariable processes in the chemical and process industry because of their simplicity and the availability of tuning rules. However, many of these tuning rules require accurate process models or critical data. Astrom and Hagglund4 first developed an autotuning procedure for single-input single-output (SISO) processes to estimate the critical gain and frequency from a relay test. They observed that a process could be made to oscillate with period Pu if the system has a phase lag of at least π at high frequencies. The critical gain and frequency can then be estimated as

Ku ) 4h/πa

(1)

ωu ) 2π/Pu

(2)

and

respectively, where h is the relay height and a is the amplitude of the oscillations. With the critical data, one may then use suitable tuning rules to find the PID controller settings for the process. Major advantages of the autotuning method are that no knowledge of the system is required a priori and it involves only a closedloop test. One of the earliest applications of the autotuning procedure to chemical processes is by Luyben,5 who proposed a method (known as the ATV method) to obtain transfer function models for highly nonlinear distillation columns. Subsequently, other researchers improved on Luyben’s ATV procedure. The improved ATV method of Li et al.6 removes the need for the * Corresponding author. E-mail: [email protected]. Tel: (65) 874 2187. Fax: (65) 779 1936.

steady-state gain to be known a priori, by using one or more additional relay tests. Chang et al.7 derived analytical expressions for the amplitude ratio and period of oscillation, while Shen et al.8 used a biased relay to obtain the steady-state gain and critical data. These data were combined with the ATV method to obtain better process models. Other procedures by Sung et al.,9 Huang et al.,10 Wang et al.,11 and Sung et al.12 incorporated relay tests to identify parametric models. Model-based controller tuning techniques may then be used. A standard autotuning procedure identifies only one point on the negative real axis of the Nyquist plot. To identify a point in the third quadrant of the Nyquist curve, Friman and Waller13 developed a twochannel relay (TCR) technique. Scali et al.14 extended the improved ATV method to completely unknown processes. In this method (known as the ATV+ method), three relay tests are performed, with an additional time delay added to two of the relay tests to determine two other points in the third quadrant of the Nyquist plot. The tuning procedure developed by Semino and Scali1 is then used to obtain suitable PI controller parameters. All of the above methods are targeted at SISO systems. The use of autotuning for multivariable processes to obtain proportional-integral (PI) controller parameters was independently proposed by Loh et al.15 and Shen and Yu.16 The former combined sequential loop closing with relay autotuning to obtain the critical gains and frequencies of individual SISO loops while accounting for process interactions. Ziegler-Nichols (Z-N) tuning rules were then used to determine suitable PI controller parameters. The sequential identification and autotuning method for multiinput multioutput (MIMO) systems by Shen and Yu16 is similar, except for the controller tuning method adopted. Noting the potential stability problem in tuning underdamped systems by the Z-N method, which may arise in sequential loop closing,

10.1021/ie010712d CCC: $22.00 © 2002 American Chemical Society Published on Web 07/20/2002

4606

Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002

Shen and Yu advocated using a modified Z-N method with a conservative detuning factor. Using these sequential loop tuning methods, only n sets of critical data are required to obtain good PI controller settings for the multivariable system without a priori knowledge of the process. Semino and Scali17 examined the sequential autotuning procedure of Loh et al.15 and extended it to MIMO systems with negative RGA, unstable systems and systems with cascade. They also noted the limitation of Z-N settings for processes with complex dynamics and very small delay. Subsequently, Semino and Scali1 extended the improved ATV method of Li et al.6 to MIMO systems by following the sequential autotuning procedure. The authors also developed a model-based controller tuning procedure that gave stable and acceptable responses. A comparison of the various autotune identification techniques for multivariable processes, namely, Shen and Yu’s16 (SY) method, the TCR method of Friman and Waller,13 and the ATV+ technique of Scali et al.,14 was carried out by Parabita et al.2 The TCR and ATV+ techniques were originally presented in SISO context but can be combined with a sequential loop closing approach to tune multivariable processes. From the results on three 2 × 2 examples, Parabita et al.2 concluded that the SY method does not always guarantee a stable response, while the TCR and ATV+ methods yielded comparable responses. However, the ATV+ technique had the longest total duration of tests because of the additional relay tests required. Autotuning methods are often computationally simple, but they involve a longer duration of relay tests because of the need to wait for steady oscillations, and some methods require more than one relay test for each loop tuning. This is particularly disadvantageous for chemical processes because of their slow dynamics. As was shown by Viswanathan and Rangaiah,18 the first two cycles of process response in a relay test along with relay output are sufficient to accurately identify second order plus dead time (SOPDT) models for SISO processes. They have successfully used time-domain curve fitting along with global optimization for model identification. The approach was later extended to simultaneous identification of 2 × 2 processes (Viswanathan et al.19). These findings provide the motivation to propose the use of time-domain curve fitting via least squares in the autotuning of multivariable processes for identifying a suitable transfer function model from the initial relay response and thereby reduce the duration of relay tests. The success and reliability of this approach needs to be investigated because multivariable processes will have interactions among loops, which will affect model identification from the initial relay response, and the suitability of the identified model for controller tuning is not known. Further, in all of the literature surveyed above on sequential relay tuning of multivariable processes, only Yu20 studied a system greater than 2 × 2. He successfully applied Shen and Yu’s16 sequential relay tuning technique to a 3 × 3 system. Hence, the main objectives of this study are to evaluate the efficacy of curve fitting the initial response for autotuning of multivariable processes, to apply and evaluate ATV+ for 3 × 3 and 4 × 4 processes, and to compare the efficacy of the proposed approach with that of ATV+, TCR, or Luyben’s3 BLT techniques. Several 2 × 2, 3 × 3, and 4 × 4 processes and measurement noise are

Figure 1. Simplified SISO closed-loop relay feedback system.

considered in the study. The results show that the proposed methodology for autotuning with least squares (ATLS) is effective for tuning of MIMO systems. 2. ATLS The autotuning method of Loh et al.15 and Shen and Yu16 treats the multivariable process as a series of SISO feedback loops. Relay autotuning is done for individual loops while taking into account the interactions between loops. This is achieved by the following steps for an nthorder MIMO system: 1. Rank the loops according to their speed, and the fastest loop is to be tuned first. 2. Perform a closed-loop relay test on the first loop, keeping all other loops open. 3. From the critical gain and frequency information obtained from the steady oscillations, determine suitable PID controller parameters. 4. Close the first loop with the controller, and conduct a relay test on the second loop while keeping other loops open. Select the controller parameters for the second loop. 5. Close the second loop as well, and carry out a relay test and controller design for the next loop until all n loops are covered. This completes one round. 6. The above process is repeated starting from the first loop (step 2) until the controller parameters converge. Ranking of loops in step 1 is possible based on the nature of the process and earlier experiences on similar processes. However, for some situations, the relative speed of the loops may be unknown, and hence the loops cannot be ranked according to speed. In these cases, the loops may be ranked arbitrarily because the method will still work albeit with an increase in the number of iterations required for convergence (Loh et al.15). The proposed ATLS method follows the same steps as above except step 3 is modified. Instead of waiting for steady oscillations, a suitable transfer function model is identified by curve fitting the initial response of the relay test and then the identified model is used to tune the PID controller using an appropriate technique. The methodology for this is developed starting from the relay test in a SISO system (Figure 1), which is the case in step 2 of the above autotuning procedure in the first iteration. Assume that relay output, m(t), and process output, y(t), from a relay test on this feedback loop are available. The goal of identification by time-domain curve fitting is to find the best model for the process that minimizes the mean sum of squares of errors (MSE) between the actual and calculated closed-loop responses:

MSE )

1

N

∑

N k)1

[

]

ya(tk) - yc(tk) ypeak

2

(3)

Here, ya(tk) and yc(tk) refer to the actual and calculated

Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002 4607

Figure 2. A 2 × 2 system with a controller in the first loop and a relay in the second loop.

closed-loop responses at time tk (an integer multiple of the sampling period), respectively, N is the number of sampling periods in the test duration Tf, and ypeak is the peak value of ya(tk) used. In the examples presented later, only the first two oscillations of the relay response were used for model identification. The predicted response, yc, can be calculated assuming the process Gp is approximated by a model such as a SOPDT model:

Kme-θms

Gm(s) )

τm2s2 + 2τmζms + 1

(4)

where Km is the model gain, τm is the time constant, ζm is the damping factor, and θm is the time delay. (As will be discussed later, a first-order integrating model may have to be used in some situations if the curve fitting using a SOPDT is unsatisfactory.) Calculations for yc can be conveniently done using a software such as SIMULINK. The actual relay output, M(t), is used to excite the simulation model in order to improve the accuracy of the SOPDT model, particularly in the presence of measurement noise.18 For a 2 × 2 process, the ATLS procedure begins by closing loop 1 with a relay and placing loop 2 in manual mode. A SOPDT model is identified for loop 1 (Gm1) by curve fitting/least squares as described above. Suitable PI controller parameters are then obtained using Semino and Scali’s1 tuning technique. According to this, the integral time is

∑i τi) +

τI ) (

θm 2

(5)

where the summation is for all time constants in the process model (eq 4). The controller gain Kc is found such that the complementary sensitivity function has a maximum of 1.26:

max ω

|

|

Gm C ) 1.26 1 + GmC

(6)

where C is the decentralized PI controller. After the first loop is tuned and closed (considered as one iteration), the second loop is closed with a relay (Figure 2) and another test is conducted. To identify a SOPDT model for the second loop (Gm2) from the relay output (m2) and the process output (y2), time-domain curve fitting is used and the objective function to be minimized is the same as eq 3. Note that the model Gm2 includes interaction between loops 1 and 2 as well. Once the model is

obtained, controller for loop 2 is tuned and the loop is closed. The second round then begins with a relay on the first loop, except that this time the second loop is closed. The above procedure is repeated until the controller parameters for loop 1 or 2 deviate by less than 10% from the previous round, which is the criterion used by Parabita et al.2 The procedure for 2 × 2 processes can be generalized to an n × n system, in which case each round consists of relay testing and tuning n loops sequentially. Optimization techniques used in this study for minimizing MSE (eq 3) are genetic algorithms (GAs) and the Nelder and Mead (NM) method. The use of local optimization alone, as was observed during the initial trials, does not guarantee that the global minimum is found unless very good initial estimates are available. GA is a global optimization technique suitable for cases where the objective function is nonconvex and multiple minima are expected. The five main steps in a typical GA are initialization, evaluation, selection, recombination, and termination. More details on GAs are available in Goldberg,21 Haupt and Haupt,22 and Viswanathan et al.19 A GA with floating-point coding, population size ) 100, maximum number of generations ) 50, tournament selection with a tournament size of 5, shape parameter ) 3, probability of crossover ) 0.6, and probability of mutation ) 0.2 was used in this study. These parameters were chosen to achieve both efficiency and a high probability of obtaining the global minimum. Generally, the GA does not find the global optimum accurately although it reliably reaches the vicinity of the global minimum. Hence, the solution obtained using a GA needs to be further refined using a local optimization technique such as the NM method. The NM method is a direct search algorithm not requiring derivatives of the objective function and has been used in controller tuning9 and several other applications.23 For model identification in the first round of relay tests, the solution from the GA is used as the initial estimate for the NM method to obtain the global minimum reliably and accurately. This combined strategy also obviates the need for good initial estimates for the success of local optimization methods. However, it is not necessary to use a GA in the second and subsequent rounds because good model parameter estimates for the NM method are available from the first round. This reduces the total computation time required. Hence, in the present study, a GA was used in the first round only. A MATLAB program for the GA used earlier19 and the NM method from the optimization toolbox of MATLAB were employed. Appropriate models for simulating MIMO systems with some or all loops closed were developed using SIMULINK. To simulate a typical industrial process, samplers were used in the feedback loop with a sampling time of about 1/20 of the largest process time constant. Integration of the governing differential equations was carried out using a fourthorder Runge-Kutta method with a reasonable step size for each example. Reported CPU times in the tables refer to the computational time in minutes on a Pentium II/300 MHz computer. 3. Results and Discussion For studying ATLS, the eight examples given in Chart 1 were considered:

4608

Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002

Chart 1

Example WB:

Example TY:

(

12.8e-s -18.9e-3s 16.7s + 1 21s + 1 Gp1(s) ) 6.6e-7s -19.4e-3s 10.9s + 1 14.4s + 1

(

)

(7)

Gp3(s) ) -0.1153(10s + 1)e

-0.1s

-2s

0.2429e (33s + 1)2 0.2429e-0.17s (44s + 1)(20s + 1)

3

(4s + 1) -0.087e-12.6s (43s + 1)(22s + 1)

(

Example T4: Gp5(s) ) -0.71s

-60s

-2.24s

-2.986e -5.24e -5.984e 66.67s + 1 400s + 1 14.29s + 1 0.0204e-0.59s -0.33e-0.68s 2.38e-0.42s (7.14s + 1)2 (2.38s + 1)2 (1.43s + 1)2 -11.3e-3.79s -9.811e-1.59s 0.374e-7.75s 22.22s + 1 (21.74s + 1)2 11.36s + 1

(

)

)

Example WW:

(

Example OR:

(9)

Example A2:

(

)

(8)

Gp4(s) )

0.66e-2.6s 6.7s + 1 1.11e-6.5s 3.25s + 1 -34.68e-9.2s 8.15s + 1

-0.61e-3.5s -0.0049e-s 8.64s + 1 9.06s + 1 -2.36e-3s -0.01e-1.2s 5s + 1 7.09s + 1 0.87(11.61s + 1)e-s 46.2e-9.4s 10.9s + 1 (3.89s + 1)(18.8s + 1)

Example NI: Gp6(s) )

(

0.5 -0.6 0.1 1 0.2 0.8 0.3 (s + 1)(2s + 1)2(0.5s + 1) -1.0 0.1 1.0 (11)

-9.811e-1.59s 0.374e-7.75s 11.36s + 1 22.22s + 1

Example DL:

(

-0.101e-12s 0.126e-6s 60s + 1 (48s + 1)(45s + 1) Gp2(s) ) 0.094e-8s -0.12e-8s 38s + 1 35s + 1

5.984e-2.24s 14.29s + 1 Gp7(s) ) 2.38e-0.42s (1.43s + 1)2 -11.3e-3.79s (21.74s + 1)2

-3.79s -2.368e-27.33s -11.3e 33.3s + 1 (21.74s + 1)2 -8.72s 0.422e 5.24e-60s 2 400s + 1 (250s + 1) -0.33e-0.68s 0.513e-s s+1 (2.38s + 1)2

-1.986e-0.71s 66.67s + 1 0.0204e-0.59s (7.14s + 1)2 -0.176e-0.48s 15.54e-s s+1 (6.9s + 1)2

4.48e-0.52s 11.11s + 1

4.09e-1.3s -6.36e-0.2s -0.25e-0.4s (33s + 1)(8.3s + 1) (31.6s + 1)(20s + 1) 21s + 1 -5s -4s -1.01s -0.05e -4.17e 6.93e 45s + 1 44.6s + 1 (34.5s + 1)2 Gp8(s) ) -17s -11s 5.11e -1.73e 4.61e-1.02s 18.5s + 1 (13s + 1)2 (13.3s + 1)2 -11.18e-2.6s 14.04e-0.02s -0.1e-0.05s (43s + 1)(6.5s + 1) (45s + 1)(10s + 1) (31.6s + 1)(5s + 1)

The first three examples are 2 × 2 systems, which were employed by Parabita et al.2 for comparing three autotuning techniques. Example WB is Wood and Berry’s model of a pilot-scale distillation column, example WW is from Wardle and Wood, and example TY is a model of stabilizer reported by Tyreus; model parameters in these three examples were taken from Luyben.3 Examples OR and T4 are 3 × 3 systems reported by Ogunnaike and Ray and Tyreus (case 4), with model parameters taken from Luyben3 and Yu,20 respectively. Example NI is a 3 × 3 system with higher order transfer functions from Niederlinski.24 Examples DL and A2 are 4 × 4 systems by Doukas and Luyben and Alatiqi (case 2), with model parameters taken from Luyben.3 3.1. Two-Input Two-Output Systems. In all of the three 2 × 2 examples, loop 1 is the faster loop and is thus tuned first. The results from the final round of model identification by the proposed ATLS for example WB are summarized in Table 1. Two rounds consisting of four iterations were required for autotuning examples WB and WW while example TY required three iterations only (Table 2). For all the 2 × 2 examples, a near perfect model was obtained by ATLS for the SISO loop in the first iteration, and the calculated responses in

)

)

)

(10)

(12)

(13)

-0.49e-5s (22s + 1)2 1.53e-2.8s 48s + 1 -5.48e-0.5s 15s + 1 4.49e-0.6s (48s + 1)(6.3s + 1)

)

(14)

Table 1. Final Models Identified in Autotuning of Some Examples by ATLS example WB WB with noise OR OR with noise

model

Km

Gm1 Gm2 Gm1 Gm2 Gm1 Gm2 Gm3 Gm1 Gm2 Gm3

10.8 -6.35 9.14 -7.85 0.612 -2.45 0.630 0.693 -2.56 0.569

τm

ζm

1.31 5.24 2.53 0.560 1.56 3.71 2.92 0.559 0.278 12.9 0.906 3.27 0.181 12.7 0.235 15.4 2.63 1.01 0.537 3.80

θm

MSE

0.911 2.46 0.880 2.05 2.59 2.96 1.000 2.73 2.16 0.912

5.90 × 10-4 3.51 × 10-3 3.71 × 10-2 1.23 × 10-2 3.05 × 10-2 2.99 × 10-2 1.50 × 10-5 6.42 × 10-2 0.399 3.16 × 10-2

subsequent iterations match the actual process responses well (see MSE values for example WB in Table 1). Semino and Scali’s1 tuning method was used for finding the PI controller parameters in each loop, and the final settings are shown in Table 2. The ATV+ method2 was also applied to the three TITO examples. This method requires three relay tests (until steady oscillations are established) in each iteration. As can be seen from Table 2, the number of iterations required and/or controller settings obtained by ATV+ may differ from those by ATLS (due to different models).

Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002 4609

Figure 3. Performance of controllers tuned by ATLS (solid line), the ATV+ method (dashed line), and the BLT method (dots and crosses) for a set-point change in (a) y1 and (b) y2 of example WB. Table 2. Comparison of ATLS with ATV+ and BLT: 2 × 2 Examples IAE for step change in example

method

Kc1

τI1

Kc2

τI2

WB

ATLS ATV+ BLT(+4) ATLS ATV+ BLT(+4) ATLS ATV+ BLT(+4)

0.87 0.68 0.375 56.6 42.1 27.4 -16.9 -13.1 -16.6

14.2 5.98 8.29 45.4 48.1 41.4 5.22 4.62 20.6

-0.078 -0.077 -0.075 -23.1 -18.0 -13.3 16.7 11.0 70.6

4.06 3.93 23.6 29.4 26.5 52.9 64.1 47.9 80.1

WW TY

The performance of the controllers tuned by ATLS and the ATV+ method for the three TITO examples was studied for a step change in the set point of each loop. A unit step change was first applied to the set point of the first loop, and the process responses were recorded; then another step change was applied to the second loop. A comparison of the results (Figure 3 as well as IAE values in Table 2) reveals that the performance of controllers tuned by ATLS is comparable to that tuned by ATV+. The IAE value is computed as the sum of the absolute difference between the closed-loop process response and the final steady-state value. (For examples WB, WW, TY, OR, and NI, IAE was computed for time t ) 0-100 units, whereas for example T4, the IAE calculation was up to 400 units to allow the process responses to reach steady state. IAE was computed to 200 units for examples DL and A2.) The response curves

total Tf

total CPU time

no. of rounds (iterations)

89.2 638

30

2 (4) 3 (6)

267 934

76

2 (4) 2 (3)

15.6 1381

10

2 (3) 2 (4)

loop 1

loop 2

y1

y2

y1

y2

3.6 4.4 4.5 17.8 18.2 19.4 5.0 5.3 14.1

6.6 8.0 14.6 30.8 24.7 30.6 9.0 11.4 2.5

2.5 1.7 3.0 6.2 7.9 8.2 3.4 3.8 13.1

8.4 9.0 28.4 22.8 22.9 38.2 26.8 27.6 23.6

for the ATV+ method for example WB (Figure 3) are similar to those shown by Parabita et al.,2 thus confirming the validity of the present study. Parabita et al.2 reported total Tf of 591, 1146, and 721 for tuning examples WB, WW, and TY, respectively, by the ATV+ method. These are somewhat different from those obtained in this study (Table 2) probably because of differences in the criterion for determining steady oscillations. Parabita et al.2 also reported total Tf for the TCR method as 124, 423, and 635 for examples WB, WW, and TY, respectively. Compared to these Tf for the ATV+ and TCR methods, results in Table 2 show that ATLS requires significantly shorter test duration for all of the three TITO examples tested. 3.2. Three-Input Three-Output Systems. For 3 × 3 and 4 × 4 systems, we found that it was necessary to repeat the identification and tuning procedure until the

4610

Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002

Table 3. Comparison of ATLS with ATV+ and BLT: 3 × 3 Examples

example

method

Kc1

τI1

Kc2

τI2

Kc3

τI3

OR

ATLS ATV+ BLT(+6) ATLS ATV+ BLT(+6) ATLS ATV+

1.68 1.09 1.51 -11.1 -15.7 -11.3 0.955 0.590

16.9 4.46 16.4 136.4 13.6 7.09 6.18 4.22

-0.300 -0.260 -0.295 -2.38 -4.72 -3.52 0.560 1.80

14.8 3.69 18.0 6.12 4.11 14.5 5.60 3.59

2.61 4.07 2.63 -0.19 -0.03 -0.182 0.425 1.41

10.2 3.63 6.61 26.6 13.6 15.1 5.27 3.46

T4 NI

variations in all three (or four) sets of controller parameters compared to those in the previous round were less than 10%. (This is in contrast to the 2 × 2 systems studied, for which it was sufficient to stop the iterations when controller settings for either loop 1 or 2 are within 10% of those in the previous round.) If two of the three controllers in a 3 × 3 or 4 × 4 system converge, it does not necessarily mean that the third controller will converge too. For 3 × 3 systems, the sampling time used is 0.1, lower bounds are Km ) -10, τm ) 0.1, ζm ) 0.1, and θm ) 0, and upper bounds are Km ) 10, τm ) 50, ζm ) 25, and θm ) 10. The integration step size is 0.02 for examples OR and T4 and 0.1 for example NI. For examples T4 and NI, initial tests showed that accurate models were recovered using ATLS with small MSE values. However, servo tests done using the PI controller settings obtained from ATLS yielded unstable responses. Semino and Scali1 proposed controller tuning by eqs 5 and 6 for SISO and used it successfully for 2 × 2 processes only. It was suspected that they may be too aggressive for higher order systems, and a simple modification was made to the tuning equations by including a detuning factor. After a few trials using various values, a detuning factor of 2 was selected. Therefore, the modified tuning formulas used in ATLS for all 3 × 3 examples are Kc′ ) Kc/2 and τI′ ) 2τI, where Kc and τI were obtained from eqs 5 and 6. However, the original Kc and τI are used in the ATV+ method. For example OR, loop 3 is the fastest, followed by loop 1 and then loop 2. Thus, the tuning sequence is loop 3, loop 1, and loop 2. The small MSE values obtained (Table 1) show that the final models simulate the dynamics of the actual process well. In the SISO identification of Gp33 in the first iteration, the ATV+ method identified a less accurate FOPDT model with parameters Km3 ) 0.53, τm3 ) 3.16, and θm3 ) 1.12 compared to Km3 ) 0.66, τm3 ) 0.17, ζm3 ) 13.9 (corresponding to two time constants ) 4.72 and 0.0061), and θm3 ) 0.99 obtained by ATLS (see the 3,3 element in the transfer function matrix, eq 10). Six iterations were required for autotuning example OR by ATLS compared to 11 iterations for the ATV+ method (Table 3). The total test duration for ATLS was much smaller (109 versus 1544 for ATV+). To evaluate the controllers designed, a step change of 0.1 units was introduced in turn in loops 1 and 2 (whose controlled variables are mole fractions), and then a step change of 5 units was introduced in loop 3 (whose controlled variable is temperature). A comparison between the closed-loop responses of controllers tuned by ATLS and the ATV+ method is shown in Figure 4, which indicates that the controller settings obtained by ATLS are inferior to those by ATV+. This is confirmed by the sum of IAEs

total Tf

total CPU time

no. of rounds (iterations)

109 1544

29

2 (6) 4 (11)

83.1 1851

16

3 (8) 2 (6)

86.1 3034

8

3 (7) 5 (15)

sum of IAE for all loops for step change in loop 1 loop 2 loop 3 29.9 14.1 22.6 19.2 22.8 17.2 32.3 16.3

10.1 4.3 8.0 112 437 64.0 25.3 19.7

34.5 18.0 22.9 38.4 43.6 28.8 14.7 14.1

of the individual loops for step change in a particular loop, summarized in Table 3. In example T4, the identification sequence is loop 1, loop 2, and finally loop 3. Difficulties were encountered in this example when obtaining the model for loop 3 in the first round. It was found that ATLS failed to find a reasonable SOPDT model because the model gain (Km) was unreasonably large. Further inspection of the results showed that a first-order integrating model would yield a better fit for the relay response in loop 3. The PI controller settings for a first-order integrating process were calculated using the method of Tyreus and Luyben,25 where

Kc ) 0.487/KPθ

(15)

τI ) 8.75θ

(16)

When a first-order integrating model is used for loop 3 and the detuning factor of 2 is applied to eqs 15 and 16, ATLS converged after eight iterations. Servo control tests were performed using the final PI controller settings, and the resulting IAE values are summarized in Table 3. Final PI settings obtained using the ATV+ method on example T4 are shown in Table 3. Difficulties were encountered in this case because there was a change in the magnitude of the relay response in loop 3. Attempts to use only the initial response failed to yield a valid model; all six models obtained for loop 3 using the ATV+ method had a pole in the right half-plane and were thus open-loop unstable. Therefore, all six candidate models were invalid because the ATV methods assume stable processes (Luyben5). The ATV+ method succeeded only if the subsequent steady oscillations after the change in the process magnitude were used in the model identification. In the SISO identification of Gp11 in the first iteration, the ATV+ method identified a FOPDT model with parameters Km1 ) -0.72, τm1 ) 13.28, and θm1 ) 0.82. This is a very poor approximation of the actual model compared to the ATLS model with Km1 ) -3.04, τm1 ) 1.42, ζm1 ) 23.89 (corresponding to two time constants ) 68.04 and 0.03), and θm1 ) 0.68. Furthermore, even though only six iterations were needed in the ATV+ method, the total test duration was more than 20 times that required for ATLS (Table 3). IAE values for responses to unit step changes in the set point (Table 3) reveal that the controller settings obtained by ATLS for example T4 are somewhat better than those by the ATV+ method. For example NI, only seven iterations were required for ATLS and the final controller settings are shown in Table 3. Using the ATV+ method, 15 iterations were required, and the total Tf needed is substantially more than that required by ATLS. There was also a drastic

Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002 4611

Figure 4. Performance of controllers tuned by ATLS (solid line), the ATV+ method (dashed line), and the BLT method (dots and crosses) for a set-point change in (a) y1, (b) y2, and (c) y3 of example OR.

Figure 5. Relay response for loop 1, round 4 using the ATV+ method for example NI.

increase in the magnitude of the response in loop 1 after the third round by the ATV+ method, as shown in Figure 5 for round 4. An important observation from Figure 5 is that the first few oscillations are similar in magnitude, which may mislead one into thinking that the oscillations are steady, particularly in the presence of noise. Critical gain and frequency obtained from these

initial oscillations may not be appropriate. Final controller parameters obtained using the ATV+ method are shown in Table 3. Servo tests using unit step changes in the set point showed that the ATLS settings were sluggish and less oscillatory compared to ATV+, and this is reflected in the somewhat larger IAE values for ATLS compared to ATV+ (Table 3). 3.3. Four-Input Four-Output Systems. For examples DL and A2, the lower bounds are Km ) -20, τm ) 0.1, ζm ) 0.1, and θm ) 0 and upper bounds are Km ) 20, τm ) 80, ζm ) 25, and θm ) 10. In example DL, the sampling time and integration step size are 0.1, and for example A2, the sampling time and integration step size are 1.0. Initial attempts in using ATLS on the two 4 × 4 examples failed after a few iterations because the PI controller settings obtained caused the process response in the next iteration to become unstable, even though the recovered models had very small MSE values, which indicate that they accurately represent process dynamics. Hence, as before, a detuning factor of 2 was applied to the original tuning relations (eqs 5 and 6). For example DL, the tuning sequence is loop 4, loop 2, loop 3, and loop 1. The small MSE values obtained show that the final models simulate the dynamics of the actual process very well, and the final controller settings

4612

Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002

Figure 6. Performance of controllers tuned by ATLS (solid line), and the BLT method (dots and crosses) for a set-point change in (a) y1, (b) y2, (c) y3, and (d) y4 of example DL. Table 4. Comparison of ATLS with BLT: 4 × 4 Examples

example DL A2

method

Kc1

τI1

Kc2

τI2

Kc3

τI3

ATLS -0.245 32.6 -15.8 134.1 0.09 2.88 BLT (+8) -0.118 23.5 -7.26 11.0 0.429 12.1 ATLS 0.675 14.9 1.49 26.0 0.960 28.5 BLT (+8) 0.923 61.7 1.16 13.2 0.727 13.2

(Table 4) yield good control. Attempts to use the ATV+ method on this example failed when the PI controller settings obtained in previous iterations caused the relay response in the fifth iteration to become unstable. Hence, the controller parameters obtained by ATLS were compared with Luyben’s3 settings obtained using the BLT (+8) technique (Table 4). Unit step set-point changes were introduced in each loop separately and the servo responses compared. Results in Figure 6 and

Kc4

τI4

1.70 104.6 0.743 7.94 1.22 18.1 2.17 40.0

sum of IAE for all loops total no. of for step change in rounds total CPU Tf time (iterations) loop 1 loop 2 loop 3 loop 4 122

14

3 (11)

316

9

4 (15)

55.2 51.6 210.1 116.9

33.0 42.3 276.2 219.3

155.7 205.9 17.4 11.1

27.2 27.6 142.7 97.6

Table 4 show that the controller settings by ATLS provide slightly better control than those by BLT. The tuning sequence for example A2 is loop 3, loop 2, loop 4, and loop 1. The recovered models in this example also simulate the actual process dynamics very well, with MSE values less than 10-4. Final controller parameters are shown in Table 4. Application of the ATV+ method to this example failed at the fifth iteration, where the process response in the relay test with

Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002 4613

an additional dead time became unstable. As in the previous example, the performance of controller settings obtained by ATLS is compared with that obtained by the BLT (+8) technique.3 The total IAE values from servo control tests (in Table 4) show that controller settings by ATLS yield poorer control compared to those by the BLT technique. 3.4. Measurement Noise. The effect of measurement noise is studied on examples WB and OR. Uniformly distributed noise with a zero mean was introduced at the process output, as was done by Choi et al.26 The amplitude of the noise used in this study was 10% of the peak relay response. The introduction of noise may cause the relay output to switch arbitrarily between the two states. Astrom and Hagglund4 suggested using a relay with hysteresis, such that the relay switch-on and switch-off points are larger than the absolute value of the noise. Introduction of noise in example WB did not cause the relay output to oscillate arbitrarily, and thus an ordinary relay without hysteresis was used. Models identified from a noisy response differ somewhat from those obtained in the absence of noise (Table 1), and the MSE values are also larger because of noise. However, the calculated responses match the actual responses well. Final controller settings obtained using a noisy response are Kc1 ) 0.87, τI1 ) 14.21, Kc2 ) -0.067, and τI2 ) 4.29, which are comparable to those obtained in the absence of noise in Table 2. A simulation study using these two sets of controller parameters gave very similar closed-loop responses and IAE values. A relay with hysteresis had to be used in example OR because the relay output oscillated wildly in the presence of noise. The magnitude of the relay with hysteresis was 1.1 times the magnitude of the noise, which was 30% of the peak relay response. The models identified in the presence of noise are shown in Table 1, and the final controller parameters are Kc1 ) 1.44, τI1 ) 17.2, Kc2 ) -0.22, τI2 ) 12.8, Kc3 ) 2.60, and τI3 ) 9.08. All of these results as well as closed-loop responses to set-point changes are comparable to those obtained in the absence of noise (Table 3). A comparison between the calculated and actual responses (Figure 7) shows that the SOPDT models obtained represent the process and interactions very well in the presence of noise also. Thus, the results on examples WB and OR in the presence of noise show that autotuning of MIMO systems by ATLS is largely unaffected by noise. 3.5. Computational Time and Controller Tuning. The total CPU time for modeling relay responses in all iterations and rounds on a Pentium II/300 MHz personal computer is in the range of 10-30 min for the processes considered in this study, except for example WW, which required a CPU time of 76 min. This is partly because of a large Tf. In general, the CPU time depends on other factors such as the complexity of the example and the number of rounds. With the increasing availability of faster computers, these CPU times can be substantially reduced. On the other hand, dynamics of processes in the chemical and process industry are unlikely to change. This study, perhaps for the first time, reports extensive results on autotuning of 3 × 3 and 4 × 4 systems. These results show that controller tuning relations are critical for the success of sequential relay tuning. Parabita et al.2 noted that Shen and Yu’s16 technique gave unstable responses in a 2 × 2 example. Use of Semino and Scali’s1 relations for SISO, although successful for 2 × 2 examples, were aggressive and

Figure 7. Comparison of actual (dashed line) and calculated (solid line) relay responses with noise in round 2 for example OR.

unsuccessful for larger MIMO systems. There is thus a need to develop sound controller tuning relations for use in sequential relay tuning. The results of the present study further show that the magnitude of the relay response may change after several initial oscillations, which, if not realized, may lead to improper controller tuning. 3.6. Comparison with Other Methods. In addition to autotuning, there are other methods for tuning MIMO systems involving either open- or closed-loop tests. Luyben5 proposed the use of n relay tests (with only one loop closed in each test) for n × n systems. Responses of all n variables in each relay test along with steady-state gains obtained separately are employed to find n × n transfer functions. Knowing these transfer functions, decentralized PID controllers can be tuned using BLT or other techniques, which may involve a detuning factor. Loh et al.15 remarked that the detuned loop may not reject disturbances well, leading to poor transients. Semino and Scali17 noted that the basic advantage of sequential loop closing over detuning techniques is that the former takes into account the relationship between the controlled variable and the corresponding manipulated variable when all other loops are closed. As can be seen from Shen and Yu16 and the present study, some detuning of controllers based on the SISO tuning formula is necessary for successful autotuning of a variety of MIMO processes.

4614

Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002

Table 5. Experimental Time for ATLS, n Relay Tests, and n Step Tests experimental time (Tf) for example

ATLS

n relay tests

n step tests

WB WW TY OR T4 NI DL A2

89 267 16 109 83 86 122 316

35 302 87 52 45 84 35 164

54 240 374 114 280 122 620 685

For comparison, the performances of PI controllers tuned by BLT, reported by Luyben,3 are also tested. Results on 2 × 2 examples (Table 2) show that PI controller settings by BLT are generally inferior to both ATLS and ATV+ methods. For example OR, PI controllers by both ATLS and BLT(+6) methods yield similar closed-loop responses (Figure 4), and ATV+ tuning is better than both of these (Table 3). However, as can be seen by IAE values in Table 3, BLT is the best method, followed by ATLS and ATV+, for example T4. For the two 4 × 4 examples studied (Table 4), controller settings by ATLS are better for one example while those by BLT are better for another example. The experimental time for n relay tests (until steady oscillations) to find n × n transfer functions in the procedure of Luyben5 is compared with that for ATLS in Table 5. Except for WW, TY, and NI examples, n relay tests require substantially less experimental time than ATLS. In principle, time-domain curve fitting via least squares has a good potential for finding n × n transfer functions from the initial responses of n relay tests, which will further reduce the experimental time. Hence, a detailed study is required to compare autotuning of MIMO systems with other methods. Viswanathan et al.19 studied the simultaneous identification of 2 × 2 processes from closed-loop responses to step changes in the set point. Although this method is yet to be tested for larger processes, the experimental time required for it is compared with that for ATLS in Table 5. For each example, step tests were carried out using the PI controller settings found in the previous sections and the total time was calculated as the sum of the time required for the response to reach 95% of the final steady-state value when a unit step change is made in the set point of each loop in turn. This is assumed to be adequate for closed-loop identification. Results in Table 5 show that the total time required for tuning 2 × 2 processes is somewhat more by ATLS compared to that for closed-loop step tests. However, a clear advantage for using ATLS can be seen in the experimental times for 3 × 3 and 4 × 4 processes. Because of the large optimization problem involved in the process identification from closed-loop responses, the CPU time for it can be expected to be much more compared to that for ATLS. Thus, it is more advantageous to use ATLS than the closed-loop step tests. 4. Conclusions In the standard autotuning of MIMO systems, one must wait for at least two steady oscillations before measuring the amplitude and period of the limit cycle. The results of this study demonstrate that, by using ATLS, suitable models can be obtained from the first two oscillations of the relay response. There is thus no

need to wait for the process to reach steady oscillations. This significantly reduces the time required for relay tests. In the presence of noise too, ATLS yields models that are just as good for controller tuning. The performance of controllers tuned based on these models is comparable to that obtained using the ATV+ or BLT technique. Some detuning of controllers based on the SISO tuning formula is necessary for successful autotuning of a variety of MIMO processes. Compared to ATV+ and TCR, ATLS requires significantly shorter total test duration, which could translate to savings of several hours of experimentation on the generally slow processes in the chemical and process industry. Nomenclature C ) transfer function of a decentralized PI controller E ) error Gm, Gp ) transfer function of the model and process Kc ) controller gain Km ) steady-state gain for the process model M ) manipulated variable N ) number of sampling periods R ) set point s ) Laplace transform variable t ) time Tf ) test duration tk ) sampling period y ) closed-loop response yi ) closed-loop response in loop i ya ) actual closed-loop response yc ) calculated closed-loop response Greek Symbols τI ) integral time τm ) time constant of the process model θm ) time delay of the process model ζm ) damping coefficient of the process model

Literature Cited (1) Semino, D.; Scali, C. Improved Identification and Autotuning of PI Controllers for MIMO Processes by Relay Techniques. J. Process Control 1998, 8 (3), 219-227. (2) Parabita, P.; Marchetti, G.; Scali, C. Sequential Identification and Autotuning by Relay Techniques of Decentralised Controllers for MIMO Processes, ADCHEM 2000, Pisa, Italy, 2000; IFAC: pp 767-772. (3) Luyben, W. L. Simple Method for Tuning SISO Controllers in Multivariable Systems. Ind. Eng. Chem. Process Des. Dev. 1986, 25 (3), 654-660. (4) Astrom, K. J.; Hagglund, T. Automatic Tuning of Simple Regulators with Specifications on Phase and Amplitude Margins. Automatica 1984, 20 (5), 645-651. (5) Luyben, W. L. Derivation of Transfer Functions for Highly Nonlinear Distillation Columns. Ind. Eng. Chem. Res. 1987, 26 (12), 2490-2495. (6) Li, W.; Eskinat, E.; Luyben, W. L. An Improved Autotune Method. Ind. Eng. Chem. Res. 1991, 30, 1530-1541. (7) Chang, C. R.; Shen, S. H.; Yu, C. C. Derivation of Transfer Function from Relay Feedback Systems. Ind. Eng. Chem. Res. 1992, 31, 855-860. (8) Shen, S. H.; Wu, J. S.; Yu, C. C. Use of Biased-Relay Feedback for System Identification. AIChE J. 1996, 42 (4), 11741180. (9) Sung, S. W.; Lee, J.; Yi, S. H. Automatic Tuning of PID Controller Using Second-Order Plus Dead Time Model. J. Chem. Eng. Jpn. 1996, 29 (6), 990-999. (10) Huang, H. P.; Chen, C. L.; Lai, C. W.; Wang, G. B. Autotuning for Model-Based PID Control. AIChE J. 1996, 42 (9), 2687-2691. (11) Wang, Q. G.; Hang, C. C.; Zou, B. Low Order Modeling from Relay Feedback. Ind. Eng. Chem. Res. 1997, 36 (2), 375381.

Ind. Eng. Chem. Res., Vol. 41, No. 18, 2002 4615 (12) Sung, S. W.; Lee, I. B.; Lee, B. K. On-Line Process Identification and Automatic Tuning Method for PID Controllers. Chem. Eng. Sci. 1998, 53 (10), 1847-1859. (13) Friman, M.; Waller, K. V. Two-Channel Relay for Autotuning. Ind. Eng. Chem. Res. 1997, 36 (7), 2662-2671. (14) Scali, C.; Marchetti, G.; Semino, D. Relay With Additional Delay for Identification and Autotuning of Completely Unknown Processes. Ind. Eng. Chem. Res. 1999, 38 (5), 1987-1997. (15) Loh, A. P.; Hang, C. C.; Quek, C. K.; Vasnani, V. U. Autotuning of Multiloop Proportional-Integral Controllers Using Relay Feedback. Ind. Eng. Chem. Res. 1993, 32 (6), 1102-1107. (16) Shen, S. H.; Yu, C. C. Use of Relay-Feedback Test for Automatic Tuning of Multivariable Systems. AIChE J. 1994, 40 (4), 627-646. (17) Semino, D.; Scali, C. Multiloop Autotuning Using Relay Feedback: Limits and Extensions. Comput. Chem. Eng. Suppl. 1996, 20, S907-S912. (18) Viswanathan, P. K.; Rangaiah, G. P. Process Identification from Closed-Loop Response Using Optimization Methods. Chem. Eng. Res. Des. 2000, 78, 528-541. (19) Viswanathan, P. K.; Toh, W. K.; Rangaiah, G. P. ClosedLoop Identification of TITO Processes Using Time-Domain Curve Fitting and Genetic Algorithms. Ind. Eng. Chem. Res. 2001, 40, 2818-2826.

(20) Yu, C. C. Autotuning of PID Controllers: Relay Feedback Approach; Springer: London, 1999. (21) Goldberg, D. E. Genetic Algorithms in Search, Optimization and Machine Learning; Addison-Wesley: Reading, MA, 1989. (22) Haupt, R. L.; Haupt, S. E. Practical Genetic Algorithms; John Wiley: New York, 1998. (23) Walters, F. H.; Parker, L. R., Jr.; Morgan, S. L.; Deming, S. N. Sequential Simplex Optimization; CRC Press: Boca Raton, FL, 1991. (24) Niederlinski, A. A Heuristic Approach to the Design of Linear Multivariable Interacting Control Systems. Automatica 1971, 7, 691-701. (25) Chidambaram, M. Applied Process Control; Allied: New Delhi, India, 1998. (26) Choi, J. Y.; Lee, J.; Jung, J. H.; Lee, M.; Han, C. Sequential Loop Closing Identification of Multivariable Process Models. Comput. Chem. Eng. 2000, 24, 809-814.

Received for review August 28, 2001 Revised manuscript received March 1, 2002 Accepted June 11, 2002 IE010712D