Use of Inverse Repeat Sequence (IRS) for Identification in Chemical

response models of chemical process systems has been quite popular. ... identification studies conducted on three nonlinear chemical .... τ ) 1, 3, 5...
0 downloads 0 Views 132KB Size
3420

Ind. Eng. Chem. Res. 1999, 38, 3420-3429

Use of Inverse Repeat Sequence (IRS) for Identification in Chemical Process Systems Ranganathan Srinivasan† and Raghunathan Rengaswamy* Department of Chemical Engineering, Indian Institute of Technology, Bombay, Powai, Mumbai 400 076, India

Use of a pseudo random binary sequence as an input signal for identification of finite impulse response models of chemical process systems has been quite popular. In this paper, we show the utility of another binary signal, the inverse repeat sequence (IRS) signal in identification of chemical processes. The attractiveness of IRS lies in the fact that, by its design, even-order kernels in the process response are canceled out. This reduces the effect of nonlinearity in identification. The improvements that can be derived by using IRS are shown through the aid of continuous stirred tank reactor and fluidized catalytic cracking unit case studies. Further, the utility of the identified models in the model predictive control framework is also demonstrated. Introduction With increasing competitiveness, the sophistication and flexibility of the chemical processes have increased considerably. This has significantly increased the complexity of the encountered control problem, thus leading to the development and implementation of advanced control strategies such as model predictive control (MPC), generalized predictive control (GPC), etc. Practical experience with advanced control has demonstrated that process identification is an important factor in the resulting control performance. In the design of advanced controllers for industrial processes, a considerable (8090%) amount of time is spent in process identification.1 Hence, robust and efficient methods for identification can greatly enhance the acceptability and deployment of advanced control techniques for industrial processes. Several system identification methods are available to estimate process models from experimental data.2-4 Identifying parsimonious transfer function models using maximum likelihood estimation (MLE) or prediction error methods is a well-established field.3,4 Even though the open-loop behavior of most individual unit operations in chemical processes can be modeled sufficiently well with lower-order transfer functions, these models fail when the dynamic behavior of these units are coupled with the whole process consisting of interconnected units with recycles. Another problem with parsimonious parametric models is the estimation of the order of the model. Direct identification of nonparsimonious finite impulse response (FIR) models overcomes some of these problems. Because of sufficiently high order model structure, the FIR models have increased flexibility to model a large class of linear systems. Other than the number of terms to include for capturing full dynamics until steady state, no structural decisions need to be made. Furthermore, FIR models are widely used in MPC in industry because these models are intuitive and easy to understand.5 Additionally, nonparsimonious FIR models are also used in the context of subsequent development of parsimonious transfer function models. * Author to whom all correspondence should be addressed. † Part of this work was performed while the candidate was doing his M.Tech thesis in the Systems & Control group, IIT Bombay, Mumbai, India. Current address: Honeywell India Software Operations, Bangalore, India.

This has prompted several researchers to work on the problem of direct identification of nonparsimonious FIR models from input-output data.6,7 For direct identification of FIR models, many researchers have used a pseudo random binary sequence (PRBS) signal in identification, in both chemical engineering and other fields.8-14 The reasons for the wide use of PRBS are as follows: (i) PRBS is an easily implementable signal; (ii) because PRBS is a predetermined sequence, the identification experiments can be repeated; (iii) PRBS can be proved to be persistently exciting; and (iv) the duration of experimentation can be reduced considerably. A PRBS signal is usually parametrized by three parameters, viz., amplitude (R), switching time (∆t), and time period (T). An assumption made while designing the PRBS signal is that the process is linear. Most processes are to some extent nonlinear. As yet, methods identifying the nonlinearities are not well developed. As a result, advanced control approaches like MPC use linear models to characterize the process. Hence, for the purpose of control, an accurate characterization of the linear dynamics of the system will be desirable. Hence, in the identification of a linear model to characterize a nonlinear system, it is important to minimize the effects of nonlinearities in the process of identification. In this paper, through theoretical arguments and simulation studies, it is shown that this could be achieved by using a binary signal called an inverse repeat sequence (IRS). The IRS signal17-20 is generated by doubling the PRBS sequence and toggling every other digit of the doubled sequence. The outline of the paper is as follows. We first discuss the theoretical arguments which show that even-order kernels in process response are canceled out through the use of IRS to explain how the effects of nonlinearities are minimized in the identification process. We then demonstrate the practical utility of this approach through identification studies conducted on three nonlinear chemical engineering examples. We also show that better performance can be achieved in the MPC framework through the use of FIR models identified using the IRS signal due to the reduction of effects of nonlinearities. An important point to note is that the models are estimated using the least-squares technique. The theo-

10.1021/ie980761z CCC: $18.00 © 1999 American Chemical Society Published on Web 08/21/1999

Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999 3421 Table 1. Generation of Inverse Repeat SequencesMethodology9 no. of clock pulses

m sequence

PRBS u1(t)

inverse repeat sequence

inverse repeat signal u2(t)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

1 0 0 0 1 1 1 1 0 1 0 1 1 0 0

+R -R -R -R +R +R +R +R -R +R -R +R +R -R -R

1 1 0 1 1 0 1 0 0 0 0 0 1 1 0

+R +R -R +R +R -R +R -R -R -R -R -R +R +R -R

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

1 0 0 0 1 1 1 1 0 1 0 1 1 0 0

+R -R -R -R +R +R +R +R -R +R -R +R +R -R -R

0 0 1 0 0 1 0 1 1 1 1 1 0 0 1

-R -R +R -R -R +R -R +R +R +R +R +R -R -R +R

Figure 1. ACF of the inverse repeat sequence.

Figure 2. Effects of switching interval. System I with PRBS signal parameters as T ) 2.5Ts, Tsa ) ∆t/2, and amplitude R ) (1.5.

retical analysis that is presented in the paper, though, is based on correlation analysis. It can be quite easily shown that correlation analysis and least-squares solution are equivalent.3,15 Hence, most of the analysis based on the correlation method of finding FIR models extend directly to the least-squares method. PRBS, IRS, and the Effect of Nonlinearity PRBS has been a well-studied signal in the identification of chemical processes. The most important reason for the use of PRBS in identification is that its autocorrelation and cross-correlation properties closely resemble that of the white noise. Another reason for the interest in PRBS is that FIR models could be estimated with short duration experiments. For a review of PRBS and its properties, the reader is referred to Davies.16 An IRS binary sequence is obtained by doubling the maximum length PRBS and toggling every other digit of the doubled sequence.21 This is illustrated in Table 1. Over a common period of 2N∆t, the autocorrelation function (ACF), φuu(τ), is given by

φuu(τ) ) R2

))

R2 N

R2 N

(2) Applying a correlation relationship16 to the above system, we get

∫0∞h1(λ1)φuu(τ,λ1) dλl + ∫0∞∫0∞h2(λ1,λ2) φuu(τ,λ1,λ2) dλ1 dλ2 + ... + ∫0∞...∫0∞hp(λ1,...,λp) φuu(τ,λ1,...,λp) dλ1 ... dλp

φuy(τ) )

(3)

where φuu(τ,λ1,...,λp) is the pth-order ACF of the input signal u and is defined as given below:



τ)N τ ) 2, 4, 6, ... τ ) 1, 3, 5, ...

∫0∞hl(λ1) u(t-λl) dλ1 + ∫0∞∫0∞h2(λ1,λ2) u(t-λl) u(t-λ2) dλl dλ2 + ∞ ∞ ... + ∫0 ...∫0 hp(λl,...,λp) u(t-λ1) ... u(t-λp) dλl ... dλp

y(t) )

φuu(τ,λ1,...,λp) ) 1 K u(t) u(t+τ-λ1) u(t+τ-λ2) ... u(t+τ-λp) dt lim Kf∞ K 0 (4)

τ)0

) -R2

as given below (only kernels up to order p are shown):

Because the IRS is a repeating signal with period T and antisymmetric, that is (as can be seen from Table 1)

(1)

as shown in Figure 1. For a nonlinear system, the output signal y(t) can be represented as a Volterra functional series expansion

u(t) ) -u(t + T/2)

(5)

we have all of the even-order kernels (2, 4, 6, ...) of the ACF of u(t) going to zero. This could be explained by considering the 2nd-order kernel. From eq 4 and the

3422 Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999

fact that the IRS signal is repeating, 2nd-order ACF kernel could be written as

φuu(λ1,λ2) )

1 T

∫0Tu(t) u(t+τ-λ1) u(t+τ-λ2) dt

)

1 T

∫0T/2u(t) u(t+τ-λ1) u(t+τ-λ2) dt + 1 T ∫ u(t) u(t+τ-λ1) u(t+τ-λ2) dt T T/2

(6)

The two components on the right-hand side of eq 6 are exactly equal and opposite. This can be easily seen from the antisymmetric property of IRS; refer to eq 5. For the odd-order ACF of u(t), the two components become exactly equal but positive. Because of this, the even-order kernels vanish and the odd-order kernels remain. The IRS signal suppresses the effects of nonlinear terms, canceling out the contributions from even-order kernels (2, 4, 6, ...) and leaving only the linear term and terms involving kernels of order 3, 5, 7, and so on. In many practical situations, the effects of higher order kernels on the output diminish quite rapidly with increasing order, so that with the effects of the 2ndorder kernel removed, a better estimate of the linear kernel can be obtained.18,21

Figure 3. Power spectral density of PRBS with varying switching time ∆t.

Case Studies The utility of IRS in identification is demonstrated through the use of three chemical engineering case studies available from the literature. These case studies are also used in explaining some issues in designing the IRS. The systems that have been used in the simulation studies are as listed below: 1. Continuous stirred tank reactor (CSTR) as described in Uppal et al.22sSystem I. In this CSTR, process product conversion is the controlled variable (CV) and the manipulated variable (MV) is coolant temperature. 2. CSTR as described in Seborg et al.23sSystem II. This CSTR considered here is a well-known example used extensively for simulation purposes in the process control literature. Product concentration is taken as the CV, and the MV is feed flow rate. 3. Fluidized catalytic cracking unit (FCCU) as described in Hovd and Skogestad24sSystem III. A simplistic model combining the three-lump kinetic scheme of Weekman and Nace,25 the riser model of Shah et al.,26 and the regenerator model of Errazu et al.27 has been used by Hovd and Skogestad.24 The values of different parameters for the Hovd and Skogestad24 model and the spatial equations of the riser temperature are taken from Balchen et al.28 For the single-input single-output (SISO) system considered here, the manipulated variable is the air flow rate (Fa), the disturbance variable is the feed oil flow rate (Foil), and the CV is the regenerator temperature (Trg). Results of Case Studies In this section we show the improvements that could be obtained by using IRS as an input signal. As mentioned before, PRBS/IRS signals are parametrized by three parameters, viz., amplitude (R), switching time (∆t), and time period (T). In our experiments, we found that choosing good values for these parameters is crucial

Figure 4. Input amplitude effect on system I.

for good identification results. In the following discussion, we briefly discuss how these parameters were chosen for the PRBS signal. The IRS signal design is performed using the same parameters except that the time period is doubled. The period T of the PRBS is chosen equal to 2.5Ts, where Ts is the settling time of the process. This length gives a reasonable amount of data for the robust estimation of the impulse response coefficients. This correlates well with the observation that is made by Box and Jenkins.15 The switching time was found to be a crucial parameter for good identification. There has already been some work done on designing ∆t for identification.8,12,13,21 To see the effect of switching time on identification, consider Figure 2, which shows the effect of switching time on system I. In this case, no noise was added to the process output. The square error that is shown in Figure 2 is the sum of errors in prediction for both the signal that was used to identify the process and a test signal that was used to test the goodness of the model. There were many simulations that were performed, and here we show only a representative example. The switching time is important because the power spectrum of PRBS/ IRS varies considerably with respect to switching time. The power spectra of the PRBS at various ∆t values are

Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999 3423

Figure 5. Effect of noise on system I identification with IRS as the input signal.

Figure 7. (a) System I identification data/model performance. Input signal: IRS. (b) System I validation data/model performance. Input signal: random. Table 2. Comparison of PRBS and IRS in Terms of MSE MSE train data

MSE test data

no. of coefficients

PRBS IRS

system I 0.0043 0.0014

0.0071 0.0022

60 60

PRBS IRS

system II 0.0006 0.0004

0.0012 0.0004

60 60

0.3980 0.0917

0.9670 0.1077

80 80

0.5054 0.1497

1.6681 0.6544

80 80

system III manipulated variable PRBS IRS disturbance variable PRBS IRS

Figure 6. (a) System I identification data/model performance. Input signal: PRBS. (b) System I validation data/model performance. Input signal: random.

given in Figure 3. The power spectrum broadens with decreasing ∆t. From eq 8, which shows the mean-square error (MSE) in estimation,16 it is clear that input power acts as a weighting function in the least-squares estimation and hence is an important factor in identification. In this paper, because we wanted to show the improvements possible using IRS, we chose the best possible ∆t based on experimentation. In general, from simulation results, it is seen that a good value for switching time is between τp/5 and τp/10, where τp is the time constant of the process. In all of the studies, an approximate value of τp was taken from the settling time Ts with the thumb rule τp ≈ Ts/4. Ts (settling to within 98% of the steady-state gain) was experimentally obtained by giving a step test.

e2(t) )

∫-∞+∞||H(jω) - Hˆ (jω)||2Φuu(ω) dω

1 2π

(7)

Amplitude R is chosen as the largest input signal amplitude that will not cause undue disturbance to the normal plant operation. The largest input amplitude that will not cause undue disturbance to the plant can

Figure 8. (a) System II identification data/model performance. Input signal: PRBS. (b) System II validation data/model performance. Input signal: random.

be found by giving a step test which forces the output to move away from the present operating region into a different operating region and exhibit strong nonlinearities. Figure 4 shows the effect of amplitude on system I. It can be clearly seen that an input amplitude R greater than 2.0 moves the system into regions of

3424 Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999

Figure 9. (a) System II identification data/model performance. Input signal: IRS. (b) System II validation data/model performance. Input signal: random.

Figure 10. (a) System III manipulated variable identification data/model performance. Input signal: PRBS. (b) System III manipulated variable validation data/model performance. Input signal: random.

nonlinearity. In practical situations, though, such multiple tests might not be feasible, and the maximum amplitude that can be chosen would be restricted by other considerations. In systems like CSTR (system I), one could perform analytical calculations to infer the region of linearity, but this would not be feasible for more complex systems. One could also use other approaches for identifying the region of linearity based on input-output data,29 and this can be used for selecting the maximum input signal amplitude. The system noise levels are taken as 10% of the output amplitude for simulation purposes. The amplitude of the input signal chosen should be sufficient to perturb the process over the noise level. To illustrate this, consider Figure 5. With very low amplitude of the input signal, the MSE in identification increases because the noise swamps out the process response. The number of FIR coefficients should approximately equal Ts/∆t. Because Ts is an approximate number, one could still optimize the number of FIR coefficients using some criteria like AIC; the improvements, though, may

Figure 11. (a) System III manipulated variable identification data/model performance. Input signal: IRS. (b) System III manipulated variable validation data/model performance. Input signal: random.

Figure 12. (a) System III disturbance variable identification data/ model performance. Input signal: PRBS. (b) System III disturbance variable validation data/model performance. Input signal: random.

not be substantial. In cases where because of computational considerations in MPC one would like to have a smaller number of coefficients, one could optimize the number of coefficients between two objectives, viz., minimizing the number of coefficients and capturing enough dynamics for a reasonable representation of the process. System I: CSTR. Following the above discussion, the approximate values of settling time Ts and the time constant τ were found by performing a step test on system I and are as follows: (a) settling time (Ts) ≈ 6 s (b) time constant (τp) ≈ 1.5 s The PRBS parameters were fixed based on the above discussion and are given by the following: (a) amplitude 1.5 (b) time period T ≈ 2.5Ts ) 15 s (c) switching time ∆t ) 0.2 s (see Figure 2)

Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999 3425

Figure 16. Model predictive control of system III (disturbance rejection with measured disturbance). Figure 13. (a) System III disturbance variable identification data/ model performance. Input signal: IRS. (b) System III disturbance variable validation data/model performance. Input signal: random.

experiments, the sampling time is taken as half of the switching time. (a) sampling time Tsa (∆t/2) ) 0.1 s (b) n ) 7 (c) no. of coefficients M ) Ts/Tsa ) 6/0.1 ) 60

Figure 14. Model predictive control of system III (set-point tracking).

A PRBS of length N ) 75 with the specified R and ∆t was generated using the MATLAB function MLBS. The PRBS generated is applied to the system, and the FIR coefficients are then estimated. Figure 6 gives the response of the system, for both the identification data and the validation data. The CSTR taken here is highly amplitude-dependent, and the model obtained here is a good fit for the process in the amplitude region between 0 and 1.5. For the same process considered, i.e., system I, an IRS with parameters fixed to the same values as above was generated through the procedure given in Table 1 and applied. The FIR model identified was tested with validation data, where the test signal was a random signal. Figure 7 gives the response of the system, for both identification data and validation data. It could be seen from Figures 6 and 7 that the model identified using IRS does better because of the removal of even-order kernels. The performance improvement in terms of MSE is given in Table 2. It can be noticed that there is a slight persistent bias that is present in validation sequences. This is because of the inherent nonlinearity present in the system identified. In this work, because a good linear fit to the process is identified, the characterization is better than models identified using other signals, but a small bias is still present because of the characterization of an inherently nonlinear process using a linear model. System II: CSTR. Similar to system I, the settling time and time constant of system II were estimated from a step test and are as given below. (a) settling time (Ts) ≈ 1.2 h (b) time constant (τp) ≈ 0.3 h The PRBS was generated with parameters fixed as follows:

Figure 15. Model predictive control of system III (disturbance rejection with unmeasured disturbance).

The sampling time has to be less than or equal to ∆t. This is typically based on instrumentation limitations and/or the noise level in the system. In all of our

(a) amplitude 10 (b) time period T ≈ 2.5Ts ) 3.0 h (c) switching time ∆t ) 0.04 (best value based on experimentation) (d) sampling time Tsa (∆t/2) ) 0.02 h (e) n ) 7

3426 Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999 Table 3. Step Response Model for the Manipulated Variable of System III Using the PRBS Signal time (s)

step response coefficient

time (s)

step response coefficient

0.000 000 150.000 000 300.000 000 450.000 000 600.000 000 750.000 000 900.000 000 1050.000 000 1200.000 000 1350.000 000 1500.000 000 1650.000 000 1800.000 000 1950.000 000 2100.000 000 2250.000 000 2400.000 000 2550.000 000 2700.000 000 2850.000 000 3000.000 000 3150.000 000 3300.000 000 3450.000 000 3600.000 000 3750.000 000 3900.000 000

0.277 268 1.646 474 2.848 962 3.931 628 4.881 079 5.744 635 6.500 262 7.193 120 7.799 156 8.359 937 8.849 397 9.307 945 9.710 006 10.085 329 10.409 956 10.724 670 11.000 603 11.264 303 11.494 642 11.716 477 11.910 132 12.097 933 12.271 720 12.431 991 12.581 538 12.715 668 12.853 898

4050.000 000 4200.000 000 4350.000 000 4500.000 000 4650.000 000 4800.000 000 4950.000 000 5100.000 000 5250.000 000 5400.000 000 5550.000 000 5700.000 000 5850.000 000 6000.000 000 6150.000 000 6300.000 000 6450.000 000 6600.000 000 6750.000 000 6900.000 000 7050.000 000 7200.000 000 7350.000 000 7500.000 000 7650.000 000 7800.000 000 7950.000 000

12.968 820 13.081 225 13.180 543 13.286 642 13.373 285 13.465 391 13.540 934 13.612 510 13.676 208 13.747 611 13.805 779 13.860 126 13.907 344 13.954 555 13.998 375 14.032 435 14.067 697 14.091 520 14.123 952 14.152 413 14.181 475 14.205 694 14.229 102 14.247 135 14.268 637 14.277 875 14.294 915

time (s)

step response coefficient

8100.000 000 8250.000 000 8400.000 000 8550.000 000 8700.000 000 8850.000 000 9000.000 000 9150.000 000 9300.000 000 9450.000 000 9600.000 000 9750.000 000 9900.000 000 10050.000 000 10200.000 000 10350.000 000 10500.000 000 10650.000 000 10800.000 000 10950.000 000 11100.000 000 11250.000 000 11400.000 000 11550.000 000 11700.000 000 11850.000 000

14.298 732 14.313 489 14.324 289 14.335 776 14.349 436 14.357 081 14.365 753 14.371 120 14.387 893 14.394 066 14.404 476 14.404 817 14.425 047 14.426 041 14.435 387 14.427 919 14.436 526 14.429 442 14.438 054 14.429 039 14.440 088 14.425 289 14.430 668 14.418 275 14.410 645 14.394 511

Table 4. Step Response Model for the Manipulated Variable of System III Using the IRS Signal time (s)

step response coefficient

time (s)

step response coefficient

0.000 000 150.000 000 300.000 000 450.000 000 600.000 000 750.000 000 900.000 000 1050.000 000 1200.000 000 1350.000 000 1500.000 000 1650.000 000 1800.000 000 1950.000 000 2100.000 000 2250.000 000 2400.000 000 2550.000 000 2700.000 000 2850.000 000 3000.000 000 3150.000 000 3300.000 000 3450.000 000 3600.000 000 3750.000 000 3900.000 000

3.717 484 5.174 953 6.472 912 7.623 319 8.651 978 9.566 722 10.388 732 11.121 093 11.783 220 12.374 145 12.911 638 13.392 506 13.830 270 14.225 601 14.585 358 14.913 987 15.212 081 15.487 646 15.736 788 15.968 247 16.178 138 16.374 993 16.554 698 16.722 927 16.876 124 17.022 484 17.158 243

4050.000 000 4200.000 000 4350.000 000 4500.000 000 4650.000 000 4800.000 000 4950.000 000 5100.000 000 5250.000 000 5400.000 000 5550.000 000 5700.000 000 5850.000 000 6000.000 000 6150.000 000 6300.000 000 6450.000 000 6600.000 000 6750.000 000 6900.000 000 7050.000 000 7200.000 000 7350.000 000 7500.000 000 7650.000 000 7800.000 000 7950.000 000

17.283 895 17.400 133 17.508 915 17.608 124 17.703 892 17.793 018 17.877 393 17.952 088 18.027 585 18.093 757 18.159 391 18.218 660 18.277 948 18.332 597 18.383 566 18.430 490 18.475 813 18.519 810 18.560 333 18.601 279 18.638 487 18.675 592 18.710 663 18.744 206 18.776 818 18.807 738 18.838 967

(f) no. of coefficients M ) Ts/Tsa ) 1.2/0.02 ) 60 The simulation is done on a noise-free system, and the output was recorded. The identified FIR model response for identification and random test inputs as compared with the process output is shown in Figure 8. It could be observed that the identified model behaves well. An IRS signal generated was applied, and the product concentration was recorded. The FIR model identified was tested, and the model, plant responses are shown in Figure 9. It could be clearly seen that IRS indeed identifies a better model compared to PRBS. This can also be observed from the MSE given in Table 2.

time (s)

step response coefficient

8100.000 000 8250.000 000 8400.000 000 8550.000 000 8700.000 000 8850.000 000 9000.000 000 9150.000 000 9300.000 000 9450.000 000 9600.000 000 9750.000 000 9900.000 000 10050.000 000 10200.000 000 10350.000 000 10500.000 000 10650.000 000 10800.000 000 10950.000 000 11100.000 000 11250.000 000 11400.000 000 11550.000 000 11700.000 000 11850.000 000

18.867 071 18.895 132 18.917 661 18.943 362 18.965 542 18.988 500 19.005 221 19.026 760 19.043 819 19.063 320 19.076 895 19.095 746 19.109 007 19.126 167 19.140 198 19.156 004 19.167 194 19.181 289 19.188 506 19.200 346 19.203 005 19.214 532 19.219 785 19.230 015 19.232 709 19.241 381

System III: Nonlinear FCCU System. In the FCCU case study, we discuss the utility of the identified models in MPC. As mentioned before, the manipulated variable is the air flow rate (Fa), the disturbance variable is Foil, the feed oil flow rate, and the CV is Trg, the regenerator temperature. From the open-loop response of the system, the time constant and settling time of Trg with respect to both Foil and Fa are nearly the same and are as given below: (a) settling time (Ts) ≈ 12 000 s (b) time constant (τp) ≈ 3000 s A PRBS sequence was generated with the following:

Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999 3427 Table 5. Step Response Model for the Disturbance Variable of System III Using the PRBS Signal time (s)

step response coefficient

time (s)

step response coefficient

0.000 000 150.000 000 300.000 000 450.000 000 600.000 000 750.000 000 900.000 000 1050.000 000 1200.000 000 1350.000 000 1500.000 000 1650.000 000 1800.000 000 1950.000 000 2100.000 000 2250.000 000 2400.000 000 2550.000 000 2700.000 000 2850.000 000 3000.000 000 3150.000 000 3300.000 000 3450.000 000 3600.000 000 3750.000 000 3900.000 000

0.0 -0.651 376 -1.328 497 -1.694 355 -2.314 441 -2.668 648 -3.227 387 -3.499 103 -4.009 157 -4.193 514 -4.558 276 -4.790 934 -5.189 955 -5.334 822 -5.623 838 -5.762 493 -6.072 355 -6.209 748 -6.530 826 -6.624 356 -6.921 869 -6.925 681 -7.348 787 -7.341 093 -7.739 396 -7.575 966 ‚-8.048 660

4050.000 000 4200.000 000 4350.000 000 4500.000 000 4650.000 000 4800.000 000 4950.000 000 5100.000 000 5250.000 000 5400.000 000 5550.000 000 5700.000 000 5850.000 000 6000.000 000 6150.000 000 6300.000 000 6450.000 000 6600.000 000 6750.000 000 6900.000 000 7050.000 000 7200.000 000 7350.000 000 7500.000 000 7650.000 000 7800.000 000 7950.000 000

-7.882 161 -8.243 758 -8.037 224 -8.471 482 -8.211 541 -8.512 223 -8.294 011 -8.578 652 -8.314 955 -8.619 044 -8.432 193 -8.572 971 -8.365 292 -8.557 128 -8.464 907 -8.640 948 -8.425 960 -8.610 625 -8.419 680 -8.667 252 -8.524 709 -8.795 276 -8.555 787 -8.738 389 -8.570 071 -8.800 634 -8.563 703

time (s)

step response coefficient

8100.000 000 8250.000 000 8400.000 000 8550.000 000 8700.000 000 8850.000 000 9000.000 000 9150.000 000 9300.000 000 9450.000 000 9600.000 000 9750.000 000 9900.000 000 10050.000 000 10200.000 000 10350.000 000 10500.000 000 10650.000 000 10800.000 000 10950.000 000 11100.000 000 11250.000 000 11400.000 000 11550.000 000 11700.000 000 11850.000 000

-8.663 872 -8.415 507 -8.613 078 -8.423 259 -8.501 419 -8.338 467 -8.395 370 -8.275 436 -8.317 582 -8.288 992 -8.288 746 -8.161 303 -8.352 320 -8.281 934 -8.324 469 -8.139 275 -8.272 182 -8.170 420 -8.273 832 -8.119 572 -8.203 992 -7.979 578 -8.078 985 -7.9371 07 -7.9710 17 -7.703 465

Table 6. Step Response Model for the Disturbance Variable of System III Using the IRS Signal time (s)

step response coefficient

time (s)

step response coefficient

0.000 000 150.000 000 300.000 000 450.000 000 600.000 000 750.000 000 900.000 000 1050.000 000 1200.000 000 1350.000 000 1500.000 000 1650.000 000 1800.000 000 1950.000 000 2100.000 000 2250.000 000 2400.000 000 2550.000 000 2700.000 000 2850.000 000 3000.000 000 3150.000 000 3300.000 000 3450.000 000 3600.000 000 3750.000 000 3900.000 000

0.000 00 -1.497 720 -2.054 555 -2.558 116 -3.043 459 -3.482 220 -3.906 951 -4.293 264 -4.665 132 -5.004 380 -5.333 760 -5.635 873 -5.928 187 -6.196 945 -6.457 083 -6.700 528 -6.932 606 -7.151 437 -7.357 362 -7.551 573 -7.738 595 -7.919 038 -8.085 725 -8.246 548 -8.398 154 -8.548 041 -8.685 337

4050.000 000 4200.000 000 4350.000 000 4500.000 000 4650.000 000 4800.000 000 4950.000 000 5100.000 000 5250.000 000 5400.000 000 5550.000 000 5700.000 000 5850.000 000 6000.000 000 6150.000 000 6300.000 000 6450.000 000 6600.000 000 6750.000 000 6900.000 000 7050.000 000 7200.000 000 7350.000 000 7500.000 000 7650.000 000 7800.000 000 7950.000 000

-8.822 092 -8.946 095 -9.067 955 -9.179 096 -9.292 366 -9.391 411 -9.494 053 -9.582 477 -9.677 428 -9.754 501 -9.838 115 -9.909 370 -9.985 518 -10.051 412 -10.119 267 -10.177 201 -10.237 926 -10.292 769 -10.351 357 -10.402 203 -10.456 061 -10.503 622 -10.554 870 -10.597 786 -10.646 222 -10.687 406 -10.731 837

Table 7. Tuning Constants of MPC for System III MPC constant

values

control horizon (M) output horizon (P) output weight (ywt) input weight (uwt)

10 1 1 0

(a) amplitude 0.5-1.0 (b) time period T ≈ 2.5Ts ) 30 000 s (c) switching time ∆t ) 300 s (d) sampling time Tsa (∆t/2) ) 150 s (e) n ) 7 (f) no. of coefficients M ) Ts/Tsa ) 12 000/150 ) 80

time (s)

step response coefficient

8100.000 000 8250.000 000 8400.000 000 8550.000 000 8700.000 000 8850.000 000 9000.000 000 9150.000 000 9300.000 000 9450.000 000 9600.000 000 9750.000 000 9900.000 000 10050.000 000 10200.000 000 10350.000 000 10500.000 000 10650.000 000 10800.000 000 10950.000 000 11100.000 000 11250.000 000 11400.000 000 11550.000 000 11700.000 000 11850.000 000

-10.768 341 -10.812 748 -10.841 161 -10.880 879 -10.909 303 -10.948 506 -10.970 317 -11.004 422 -11.024 309 -11.058 502 -11.076 967 -11.109 024 -11.124 039 -11.150 736 -11.163 531 -11.193 896 -11.203 343 -11.231 377 -11.238 465 -11.266 808 -11.273 108 -11.300 131 -11.305 386 -11.328 235 -11.332 145 -11.356 781

A PRBS sequence was applied. Models were built for both the manipulated and disturbance variables. Figures 10 and 11 show the model performance in the case of the manipulated variable using PRBS and IRS signals, respectively. Figures 12 and 13 show the model performance in the case of the disturbance variable using PRBS and IRS signals, respectively. The improved performance in using IRS over PRBS can also be gathered from the MSE shown in Table 2. The actual performance of the identified models in the MPC framework was also examined. The MPC algorithm, the linear models, and the MPC parameters

3428 Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999

are summarized in Appendix A. The same MPC tuning parameters were used for comparison. Many simulations were performed, but only some representative results are shown here. Three different simulations are presented here: 1. A set-point change in the regenerator temperature from 971.95 to 975 K (Figure 14). 2. Rejection of disturbance in Foil assuming that the disturbance is unmeasured (Figure 15). 3. Rejection of disturbance in Foil assuming that the disturbance is measured (Figure 16). From the set-point tracking (Figure 14), it can be clearly seen that the use of the IRS model leads to very quick and smooth set-point tracking with very little overshoot. For the purpose of comparison, all of the MPC parameters are kept the same in both cases. Of course, it might be possible to improve the performance of the PRBS model by tuning the MPC parameters, but this can be achieved through an increased control effort only. Figure 15 shows an example of disturbance rejection assuming that the disturbance is unmeasured. Again smooth disturbance rejection with very little overshoot is achieved through the use of the IRS model. Similar results are seen in the case where the disturbance is assumed measured also (Figure 16). Figure 15 presents the rejection of disturbance in Foil assuming that the disturbance is unmeasured; i.e., disturbance model is not incorporated in MPC. Figure 16 rejects the disturbance faster than the previous case because of the incorporation of the disturbance model in the MPC. In both cases, though, the model generated using the IRS signal performs better. Conclusion and Future Work The accuracy of estimated FIR coefficients obtained from the regression method (or any method) is mainly influenced by the quality of the data. PRBS has been a popular input signal for identification of FIR models in chemical process systems. In this paper we have shown the utility of IRS as a perturbation signal for identification. The improvements in identification result from the fact that even-order kernels are canceled out when IRS is used as an identification signal. This cancellation of even-order kernels reduces the effect of nonlinearities in the identification experiment. The improvements that can be attained by using IRS are demonstrated through three realistic chemical engineering examples. It is also shown that the models identified using IRS perform much better in the MPC framework. There is tremendous scope for future work in this area. This work has to be extended for multiple-input multiple output (MIMO) identification. In MIMO identification, when there are different outputs with different time constants, the choice of ∆t for proper identification will become a crucial issue. More work needs to be done to address this issue. Further, the experimental duration will also become long in MIMO problems. Methods for reducing the experimental duration can also be addressed. Use of delayed and multilevel sequences could be investigated in this context. Nomenclature e2(t) ) mean-square error in estimation hi ) impulse response coefficient H(jω) ) frequency response of the actual system H ˆ (jω) ) frequency response of the estimated system

M ) no. of FIR coefficients n ) number of shift registers p ) order of the volterra kernel T ) time period of PRBS Ts ) settling time of the system Tsa ) sampling time of the system Abbreviations ACF ) auto correlation function AIC ) Akaike information criteria CCF ) cross-correlation function CSTR ) continuous stirred tank reactor FCCU ) fluidized catalytic cracking unit FIR ) finite impulse response GPC ) generalized predictive control IRS ) inverse repeat sequence MIMO ) multiple-input multiple-output MPC ) model predictive control PRBS ) pseudo random binary sequence SISO ) single-input single-output Greek Symbols φuu(τ) ) autocorrelation of u(t) at τth lag φuy(τ) ) cross-correlation between u(t) and y(t) at τth lag Φuu(ω) ) power spectral density of the signal u(t) ∆t ) switching interval of PRBS R ) amplitude of the PRBS signal σ ) standard deviation of white noise τ ) lag at which ACF and FIR are estimated τp ) approximate time constant of the process

Appendix A DMC algorithm is used in the MPC implementation. The step response models obtained from PRBS and IRS signals for system III are given in Tables 3-6. The MPC tuning parameters are kept constant for both PRBS and IRS models for comparison of model performance. Table 7 gives the tuning constants of MPC. Literature Cited (1) Andersen, H. W.; Ramussen, K. H.; Jorgensen, S. B. Advances in Process Identification. CPC IV Proc. 1994, 237-269. (2) Wellstead, P. E. Non-parametric methods of system identification. Automatica 1981, 17, 55. (3) Soderstorm, T.; Stoica, P. System Identification; PrenticeHall: Englewood Cliffs, NJ, 1989. (4) Ljung, L. System Identification: theory for the user; PrenticeHall: Englewood Cliffs, NJ, 1987. (5) Richalet, J. A. A.; Rault, A.; Testud, J. L. Model Predictive Hueristic control. Automatica 1978, 14, 413. (6) Ricker, N. L. The use of biased least squares estimates for parameters in discrete-time impulse response models. Ind. Eng. Chem. Res. 1988, 27, 343. (7) Dayal, B. S.; Macgregor, J. F. Identification of Impulse response models: Methods and robustness issues. Ind. Eng. Chem. Res. 1996, 35, 4078. (8) Godfrey, K. R. Theory of correlation method of dynamic analysis and its application to industrial process and nuclear plant. Meas. Control 1969, 10, 65. (9) Godfrey, K. R. Dynamic analysis of an oil refinery unit under normal operating conditions. Proc. IEE 1969, 116, 879. (10) Barker, H. A.; Davy, R. W. System identification using pseudorandom signals and discrete fourier transform. Proc. IEE 1975, 305. (11) Harris, C. L.; Mellichamp, D. A. On-line identification of process dynamics: use of multifrequency binary sequences. Ind. Eng. Chem. Process Des. Dev. 1980, 19, 166. (12) Rivera, D. E.; Pollard, J. F.; Garcia, C. E. Control-relevant Prefiltering: A systematic design approach and case study. IEEE Trans. Autom. Control 1992, 37 (7), 964.

Ind. Eng. Chem. Res., Vol. 38, No. 9, 1999 3429 (13) Palavajjhala, S.; Motard, R. L.; Joseph, B. Process Identification Using Discrete Wavelet Transforms: Design of Prefilters. AIChE J. 1996, 42 (3), 777. (14) Astrom, K. J.; Wittenmark, B. Computer-controlled Systems Theory and Design; Prentice-Hall: Englewood Cliffs, NJ, 1997. (15) Box, G. E. P.; Jenkins, G. M. Time series analysiss forecasting and control; Holden Day: CA, 1976. (16) Davies, A. D. T. System identification for self-adaptive control; Wiley: New York, 1970. (17) Barker, H. A.; Davy, R. W. Measurement of second-order volterra kernels using pseudo random ternery signals. Int. J. Control 1978, 27, 277. (18) Barker, H. A.; Obidegwu, S. N.; Pradisthayan, T. Performance of anti-symmetric pseudorandom signals based on msequences. Proc. IEE 1972, 119, 353. (19) Ream, N. Non-linear identification using inverse repeat m-sequences. Proc. IEE 1970, 117, 213. (20) Simpson, A. J.; Power, H. M. Correlation techniques for identification of nonlinear systems. Meas. Control 1972, 5, 316. (21) Godfrey, K. R. Perturbation signals for system identification; Prentice-Hall: Englewood Cliffs, NJ, 1993. (22) Uppal, A.; Ray, W. H.; Poore, A. B. On the dynamic behaviour of continuous stirred tank reactor. Chem. Eng. Sci. 1974, 29, 967.

(23) Seborg, D. E.; Edgar, T. F.; Mellichamp, D. A. Process Dynamics and Control; Wiley: New York, 1989. (24) Hovd, M.; Skogestad, S. Procedure for regulatory control structure selection with application to the FCC process. AIChE J. 1993, 39, 1938. (25) Weekman, V. W., Jr.; Nace, D. M. Kinetics of catalytic cracking selectivity in fixed, moving and fluid bed reactors. AIChE J. 1970, 16, 397. (26) Shah, Y. T.; Huling, G. T.; Parakos, J. A.; Mckinney, J. D. A kinetic model for an adiabatic transfer line catalytic cracking reactor. Ind. Eng. Chem. Process Des. Dev. 1977, 16, 81. (27) Errazu, A. F.; de lasa, H. I.; Sarti, F. Fluidized bed catalytic cracking regenerator model-grid effects. Can. J. Chem. Eng. 1979, 57, 191. (28) Balchen, J. G.; Ljungquist, D.; Strand, S. State space predictive control. Chem. Eng. Sci. 1992, 47, 787. (29) Billings, S. A; Voon, W. S. F. Structure detection and model validity tests in the identification of nonlinear systems. IEEE Proc., Part D 1983, 117, 130.

Received for review December 1, 1998 Revised manuscript received May 10, 1999 Accepted May 19, 1999 IE980761Z