An Evaluation of MIMO Input Designs for Process Identification

Ind. Eng. Chem. Res. 2004, 43, 3847-3854

3847

An Evaluation of MIMO Input Designs for Process Identification Jeremy S. Conner and Dale E. Seborg* Department of Chemical Engineering, University of California, Santa Barbara, California 93106

The identification of a process model is a critical step in the implementation of model-based control schemes. Process identification tests for multiple-input-multiple-output (MIMO) control problems can be quite time-consuming and disruptive to the process. Thus, there is significant motivation to reduce the duration of the identification tests. In this paper, we compare alternative methods for conducting the plant tests. A detailed simulation case study is used to evaluate the advantages and disadvantages of various input excitation strategies in terms of the accuracy of the estimated parameters, prediction errors, process gains, and relative gain arrays. 1. Introduction System identification experiments are often expensive and time-consuming. This is especially true in the process industries, where the process dynamics can be very slow. During plant tests, production is disrupted, and the product quality can be off-specification. Therefore, it is important to find the most expedient means of identifying the process, to minimize disruption to production. Re-identification might be necessary if the process or operating conditions change significantly. Identification experiments involving multiple inputs require special attention to avoid cross-correlation among the inputs. Rivera and Jun1 suggest using delayed versions of a single pseudo-random binary signal (PRBS) to reduce cross-correlation in a multiinput PRBS. Godfrey’s book on input signal design2 contains a wealth of information pertaining to specialized identification signals. Additionally, model predictive control applications usually allow specification of input and output constraints. Process identification when the estimated parameters are constrained has been considered by Murakami and Seborg.3 The objective of this research is to provide an in-depth evaluation of alternative methods of input excitation for multiple-input-multiple-output (MIMO) plant tests. Apparently, this interesting theoretical and practical problem has received little attention. In sequential input testing, each input is perturbed separately while the other inputs are kept at their nominal values. This approach simplifies the test administration and data interpretation by plant personnel, because the effect of each input on all of the outputs is readily apparent. Sequential input testing is the current standard for model predictive control (MPC) applications.4 An alternative method, simultaneous input testing, excites more than one input at a time. As a result, more input excitation occurs for a given time period, as compared with a sequential input test. However, the test implementation and data interpretation are more difficult. It is possible to design specialized input sequences called rotated inputs using the concept of gain directionality. Gain directionality is a term used to describe the combinations of inputs that lead to either very strong changes in the outputs or very weak changes in * To whom correspondence should be addressed. Tel.: +1-805-893-3352. Fax: +1-805-893-4731. E-mail: seborg@ engineering.ucsb.edu.

the outputs. For example, it might happen that very little change is seen in a process output if some inputs are increased and some are decreased. This combination of input moves is an example of a weak gain direction. It is important to note that gain directionality does not directly refer to whether an individual output is increasing or decreasing. Gain directions are determined by a singular value decomposition (SVD) of the steadystate gain matrix. Knowledge of gain directionality can be used to calculate a set of rotated inputs that make only very small input moves corresponding to the strong gain direction and very large input moves corresponding to the weak gain direction. The advantage is that, on average, larger input changes can be made than with any other standard excitation sequence. Andersen and Ku¨mmel5 have demonstrated that simultaneous input excitation improves the accuracy of gain directionality estimates. Koung and MacGregor6 designed input sequences that provide equal excitation in strong and weak gain directions. Zhan and Georgakis7 describe the formulation of optimal input signals for a constrained, MIMO process using rotated input signals. These signals are constructed from preliminary knowledge of the steady-state gain matrix. The rotated inputs allow more excitation to be applied in the weak direction of the process and less excitation in the strong gain direction. In other words, combinations of input moves that do not have a strong effect on the outputs are highly emphasized, and vice versa. As a result, the output variability is considerably reduced. Therefore, it is possible to apply larger input signals without violating output constraints. This type of input design can be particularly useful for identification of ill-conditioned processes. Rotated inputs are described in more detail in section 2.3. This paper provides a detailed comparison of three types of MIMO excitation: (i) sequential testing, (ii) simultaneous testing, and (iii) rotated input testing. Two main concerns are the accuracy of the identified models and the test duration required to identify a model with a prescribed accuracy. These issues are analyzed via a simulation case study. The major contribution of this research is to provide a quantitative evaluation of the relative merits of different forms of input excitation during MIMO process identification.

10.1021/ie034068+ CCC: $27.50 © 2004 American Chemical Society Published on Web 04/21/2004

3848 Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004

If the autoregressive terms are omitted, A(q) ) 1, and a finite impulse response (FIR) model results

2. Process Identification for Model Predictive Control Model predictive control is a well-known advanced control strategy that has been widely used in the petrochemical industry and elsewhere.4 However, a reasonably accurate dynamic model is required for a model predictive controller to function properly. Because of process changes such as new operating conditions, new (or malfunctioning) equipment, or even daily weather variations, a plant-model mismatch can arise. If the plant-model mismatch becomes significant, MPC system performance can degrade significantly. If allowed to continue unchecked, the loss of economic benefits can become significant, or the closed-loop system can even become unstable. Thus, process reidentification might be required to maintain optimal MPC performance. A series of preliminary experiments typically precedes the actual MPC process identification test.8 During the pretest experiments, estimates of the process gains, time constants, and time delays are obtained. The pretest is also important to ensure that plant instrumentation is working properly. The full process identification usually involves open-loop tests. Because most processes have both input and output constraints, it is necessary to ensure that the inequality constraints are satisfied during the experimental tests. Identification of large MIMO process models for MPC applications is typically done by manipulating one input at a time. MPC applications as large as 603 × 283 have been reported.4 A recommendation given by some MPC vendors is that the total process test duration should be Ttest ) 6(nu + nd) × Tmax, where nu and nd are the numbers of inputs and measured disturbance variables, respectively, and Tmax is the maximum settling time of the process.8 Identification of large models can require weeks of around-theclock process testing. If the duration of the plant test can be significantly reduced by simultaneously adjusting some or all of the inputs, significant time and cost savings could be realized. 2.1. Model Structures. The model structures used in most MPC applications are empirical linear dynamic models such as difference equations, finite step response (FSR) models, and state-space models. Each particular structure has its own advantages and disadvantages. Often, one of the major reasons for choosing a particular MPC package is the available model structures. In this paper, we assume that the MIMO model can be represented as a set of multiple-input-single-output (MISO) models. Each submodel is an autoregressive with exogenous input (ARX) model

A(q) y(k) ) B1(q) u1(k) + ‚‚‚ + Bnu(q) unu(k) + e(k) (1) where y(k) is the process output at time k, {u1(k), ..., unu(k)} are the nu process inputs, e(k) is the process noise, and q-1 is the backward shift operator. Polynomials A(q) and B(q) are defined as

A(q) ) 1 + a1q-1 + ‚‚‚ + anaq-na

(2)

Bi(q) ) bi,0 + bi,1q-1 + ‚‚‚ + bi,nbq-nb

(3)

The autoregressive order is na; the order with respect to each exogenous input is nb.

y(k) ) B1(q) u1(k) + ‚‚‚ + Bnu(q) unu(k) + e(k) (4) 2.2. Parameter Estimation. Standard least-squares estimation is commonly used for MPC system identification. It is used here to estimate ARX and FIR model parameters from process data, as briefly described below. First, the MISO model in eq 1 is converted to the standard theta format9

y(k) ) φ(k)Tθ + e(k)

(5)

where θ is the vector of unknown model parameters

θ } [a1, ..., ana, b1,0, b1,1, ..., b1,nb, b2,0, b2,1, ..., b2,nb, ..., bnu,0, bnu,1, ..., bnu,nb]T (6) and φ(k) is the regressor matrix at time k

φ } [-y(k - 1), ..., -y(k - na), u1(k), u1(k -1), ..., u1(k - nb), u2(k), ..., u2(k - 1), u2(k - nb), ..., unu(k), unu(k - 1), ..., unu(k - nb)]T (7) The output can be predicted using the system model in eq 8

yˆ (k) ) φ(k)Tθˆ

(8)

The well-known least-squares estimate of θ for N data points is given by

θˆ ) [ΦNTΦN]-1ΦNTYN

(9)

YN } [y(1) ‚‚‚ y(N)]T

(10)

ΦN } [φ(1) ‚‚‚ φ(N)]T

(11)

where

and

Of course, many alternative parameter estimation algorithms are available. For example, a more advanced algorithm can avoid the possibly ill-conditioned matrix inversion in eq 9. For more information, the reader is referred to the textbooks by Ljung9 and van Overschee and De Moor.10 2.3. Input Excitation. For MPC applications, the process identification tests are typically implemented using sequential input excitation.4 For large control problems with many inputs, sequential input excitation can lead to extremely long test durations. A few previous papers have suggested that simultaneous excitation of the inputs will lead to more efficient use of the plant testing time.6,11-13 For MIMO processes, some or all of the inputs can potentially be adjusted simultaneously, but only a few commercial MPC identification packages currently allow this option.4 However, some problems can be encountered when implementing simultaneous input excitation. First, the inputs must be uncorrelated or the parameter estimation problem can be ill-conditioned. For FIR model identification, this condition is easily detected by generating the input sequences a

Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3849

priori and checking for the degree of correlation before performing the actual experimental test. Freudenberg and Middleton14 discuss issues involved with illconditioned two-input-two-output plants. Second, when several inputs are simultaneously moved in a strong gain direction of the process, there is a danger that output constraints will be violated. This problem can be avoided in several ways, such as performing a simulation using estimated process gains from a pretest or using rotated input designs where the excitation signals are smaller in the strong gain directions and larger in the weak gain directions. Andersen and Ku¨mmel have studied gain directionality estimation for MIMO processes.5,15 Rotated input signals for process identification were proposed by Koung and MacGregor.6 To illustrate the motivation for using rotated inputs, consider the singular value decomposition (SVD) of the steady-state gain matrix K ) MΣVT of a MIMO process that relates nu inputs to ny outputs. In the following, an overbar indicates steadystate values of a variable if inputs are held constant from the current time forward.

Y h T ) KUT ) MΣVTUT

(12)

j ny ] j 2 ‚‚‚ y j1 y Y h } [y

(13)

U } [u1 u2 ‚‚‚ unu ]

(14)

y j i } [yji(1) yji(2) ‚‚‚ yji(N) ]T

(15)

uj } [uj(1) uj(2) ‚‚‚ uj(N) ]T

(16)

where

and

The scalar yji(k) values are the steady-state outputs that would result if the inputs were held constant at {uj(k)}. In eq 12, the matrix M is called the output rotation matrix, V is called the input rotation matrix, and the singular values (σ1 g σ2 g ‚‚‚ g σn) are the diagonal elements of Σ

[ ]

σ1 · · · 0 · · · Σ) · · · · · · 0 · · · σn

(17)

[

rotation matrix V and scaled by a factor of R so that the outputs do not exceed prespecified limits.

Ξ } RU ˜ VT

(19)

If the new set of rotated inputs, Ξ, is used in eq 12 instead of the original inputs, U, the modes of the steady-state gain matrix are individually excited by scaled, uncorrelated PRBS signals. Because V is orthonormal (VTV ) I), it follows that

Y h T ) MΣVTΞT ) RMΣVT(U ˜ VT)T ) RMΣU ˜ T (20)

It is advantageous to design rotated inputs so that these uncorrelated PRBS signals individually excite the modes of the steady-state gain matrix. Furthermore, the rotated input signals should be appropriately scaled so that the modes corresponding to small singular values are excited more strongly. Thus, the original inputs, {ui}, are first scaled by the singular values {σi} to give U ˜

σ1 σ1 σ1 u3 ‚‚‚ unu U ˜ } u1 u2 σ nu σ2 σ3

Figure 1. Sequential and simultaneous input sequences.

]

(18)

where the u1, ..., unu are uncorrelated unit-amplitude PRBS vector signals of length N as shown in eq 16. These new inputs, U ˜ , are then multiplied by the input

For purposes of comparison, examples of sequential and simultaneous input sequences are shown in Figure 1. A rotated input sequence is shown in Figure 2. Note the high degree of input correlation present when a rotated input design is used. These three types of input excitation are evaluated in the case study in section 4. Parametric accuracy and prediction error accuracy are compared on the basis of models identified using these input signal designs. The use of rotated inputs allows some guarantees to be placed on process outputs for simple cases. We now present the following theorem to demonstrate this concept. Theorem 1. Consider a rotated input signal, designed using eqs 12-19, with each element of U in eq 14 bounded by (1 and with R equal to the reciprocal of


It indicates the amount of variability explained by the model.17 The normalized sum squared of errors (SSE) for a scalar output y can be expressed as N

SSEN )

[y(k) - yˆ (k)]2 ∑ k)1 N

(25)

∑ [y(k) - yj]2

k)1

Then R2 is given by

R2 ) (1 - SSEN) × 100

and R2(n) is defined as the R2 value based on n-stepahead prediction errors

Figure 2. Rotated input sequences.

the largest singular value of the steady-state gain matrix K. Then, the steady-state output vector is bounded by the largest row sum of the absolute values of the elements of M in eq 12. Proof. T

T T

Y h ) MΣV Ξ

(21)

From eq 20

˜T Y h T ) RMΣU

(22)

Keeping nu significant singular values and substituting eqs 17 and 18 with R ) σ1-1 yields

[ ][ σ1 · · · 0

· Y h T ) σ1-1M · · 0

]

σ1 σ1 σ1 T · u1 u2 u3 ‚‚‚ unu · · · σ2 σ3 σnu σ · · · nu (23)

·

·

(26)

Simplifying gives

Y h T ) σ1-1M[u1σ1 u2σ1 ‚‚‚ unuσ1]T ) M[u1 u2 ‚‚‚ unu]T (24) which is, of course, bounded by the largest row sum of M, given the assumption that the elements of the column vectors, u1,..., unu, are bounded by (1. 9 3. Model Validation After a model has been identified, validation is required to demonstrate that the model is adequate to represent the true process. One of the most widely used validation techniques is to test the model by using it to predict future outputs from an independent set of validation data. For example, the prediction error sum of squares (PRESS) statistic can be used.16 Another commonly used measure of model accuracy is the coefficient of determination, R2, which is described below. Further, measures of parametric or spectral accuracy can also be defined and used for model validation. The metrics used in this paper are described in sections 3.1 and 3.2. 3.1. Prediction Residuals. The coefficient of determination, R2, is a classical measure of model accuracy.

{

N

R2(n) ) 1 -

}

[y(k) - yˆ (k| k - n)]2 ∑ k)1 N

∑ [y(k) - yj]

2

k)1

× 100%

(27)

Caution must be taken when using R2 to determine the “best” model. It is possible to make the R2 value arbitrarily close to 100% by increasing the number of model parameters. Thus, only parsimonious models should be considered. Additionally, the standard deviations of the prediction errors can be compared. Identification methods that yield model prediction errors with smaller standard deviations are preferred. 3.2. Parametric Accuracy. In simulation case studies such as the one presented in section 4, the true model of the process is known. Thus, it is possible to compare the accuracy of the estimated model parameters to the actual values. Three measures of parametric accuracy are presented below. First, the improvement in parametric accuracy, Ip, is used to characterize identification using either simultaneous input excitation or rotated inputs. For example, the improvement in parametric accuracy that occurs for simultaneous excitation relative to sequential excitation is given by

Ip }

hˆ sim - θ0|| hˆ seq - θ0|| - ||θ ||θ hˆ seq - θ0|| ||θ

× 100%

(28)

hˆ seq and θ hˆ sim where θ0 is the true parameter vector and θ are parameter estimates averaged over M repeated simulation runs

hˆ x } θ

1

M

∑θˆ x

M i)1

(29)

where the subscript x denotes either “sim” or “seq”. A positive value of Ip indicates that a better model was identified than the one obtained using standard sequential input excitation. Next, a normalized comparison of steady-state gain estimates is considered. Index JK in eq 30 is the 2-norm

Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3851 Table 1. Properties of PRBS Signals switching time (min) period (min)

PRBS 1

PRBS 2

PRBS 3

2 2050

4 4090

10 10 230

variables. Several papers have evaluated multiloop control system performance for the Wood-Berry model.20-22 4.1. Input Sequences. Several types of input sequences were used in the case study to obtain conclusions that are independent of the characteristics of any particular input sequence. Gaikwad and Rivera23 presented the following guideline for the frequency band over which excitation should be provided

Rs 1 eωe L H βsτdom τdom

Figure 3. Schematic diagram of the Wood-Berry distillation column.

of the matrix of normalized steady-state gain estimation errors, averaged over M repeated simulation runs.

JK }

1

M

∑

x ∑ ∑(

M k)1

ny

ny

i)1 j)1

)

K ˆ ij - Kij Kij

2

(30)

The third parametric accuracy measure that is considered is the relative gain array (RGA) accuracy. RGA accuracy is particularly important for control system design.18 In the special case of a 2 × 2 system, only one element of the RGA need be calculated, λ11. The average RGA estimation error over M repeated simulation runs is given by

Jλ }

1

M

∑ | λˆ 11 - λ11|

(31)

M k)1

4. Case Study: Wood-Berry Distillation Column Model A detailed simulation case study was performed to evaluate the three input excitation methods described above. The motivation for this case study was to determine the relative advantages and disadvantages of simultaneous input excitation during plant tests. The system chosen for the case study was the Wood-Berry distillation column model, a dynamic model of an eighttray pilot-scale methanol-water column at the University of Alberta. The two-input-two-output transfer function model has been widely used in the process control literature.8 Figure 3 is a schematic diagram of the distillation column, and eq 32 is the transferfunction model reported by Wood and Berry.19

[

][

12.8 e-s -18.9 e-3s XD(s) R(s) + 1 21s + 1 ) 16.7s -7s -3s XB(s) S(s) -19.4 e 6.6 e 10.9s + 1 14.4s + 1

[ ]

]

(32)

The model consists of four first-order plus time delay (FOPTD) transfer functions. The reflux and steam flow rates, R and S, are the input variables, and the distillate and bottoms compositions, XD and XB, are the output

(33)

where Rs ≈ 2 is the fractional closed-loop speed of response of the process, i.e., the ratio of the open-loop settling time to the closed-loop settling time. βs ≈ 4 is an integer representing the number of time constants that correspond to the 95% process settling time. τH dom ) 21 min and τLdom ) 10.9 min are the high and low estimates, respectively, of the dominant process time constant. Using the guideline in eq 33, the frequency range over which excitation should be provided is 0.012 rad/min e ω e 0.18 rad/min. To ensure excitation in the desired frequency range, Gaikwad and Rivera23,24 suggest that the switching time, Tsw, of a PRBS signal satisfy the following inequality

2.8τLdom Tsw e Rs For the Wood-Berry model, Tsw should thus be less than 15 min. Four types of input excitation with different spectral characteristics were considered: three PRBS sequences with the characteristics described in Table 1 and a typical MPC test sequence consisting of a series of steps of different durations. The MPC sequence consisted of a sequence of binary steps in the input with durations of 0.25, 0.5, 0.75, and 1-1.25 times the settling time of the process (approximately 90 min). The spectra of the four input sequences are shown in Figure 4. The MPC step test sequence is the strongest signal in the lowfrequency range, whereas the PRBS sequence with Tsw ) 2 min is strongest in the high-frequency range. Clearly, there is sufficient excitation over the frequency range indicated by eq 33. The critical frequencies for the two individual control loops of the Wood-Berry model are also shown in Figure 4. Input rotation was evaluated as a means to provide larger input signals while still constraining the output variables. Input rotation allows larger variations in the weak gain direction of the process, while only very small changes are made in the strong gain direction. Using the methods described in section 2.3, rotated versions of the four types of input sequences described above were evaluated for process identification. These sequences were rescaled to provide standard deviations in the process output variables that were within the 99% confidence limits for the standard deviations of the process output variables from simulations using unrotated simultanous excitation sequences. The spectra of

3852 Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 Table 2. R2 (%) Values for Identified FIR Models and Validation Data y1

Figure 4. Standard input sequence spectra.

Figure 5. Rotated input sequence spectra.

the rotated and rescaled input sequences are shown in Figure 5. It is clear that the rotation does not significantly affect the shape of the spectra, but the power levels in Figure 5 are much greater than in Figure 4 because stronger excitation is allowed in the low gain directions of the process. A comparison of the results for the rotated, sequential, and simultaneous input excitation is presented in section 4.3. 4.2. Simulation Data. Standard PRBS sequences were used in both the sequential and simultaneous input excitation. The binary levels of the PRBS were (1, and the duration was 1000 min. Thus, sequential excitation was provided for 500 min for each input, resulting in a total test duration of 1000 min. In the simultaneous tests, excitation was provided for the full 1000 min for both inputs. In the second series of simulations, the length of the simultaneous excitation was reduced until the identified model was found to have approximately the same degree of parametric accuracy as the model identified from the sequential excitation. A trial-and-error procedure was used to determine the best sequence length to match the degrees of parametric accuracy. The resulting sequence lengths are notably different for the various types of input signals. Finally, the third set of experiments consisted of 1000 min of excitation using rotated and

a

y2a

sequential simultaneous rotated

PRBS 1 44 52 50

78 ( 13 81 ( 16 79 ( 17


PRBS 2 55 59 54

82 ( 14 83 ( 14 80 ( 20


PRBS 3 60 63 61

82 ( 15 84 ( 14 85 ( 12


MPC Input 62 65 65

85 ( 14 86 ( 13 85 ( 14

Standard deviations are included for the y2 results.

rescaled input sequences. As mentioned in section 4.1, the rescaling was adjusted so that the process output standard deviations were equivalent to those of regular simultaneous input excitation, to within 99% confidence limits. Each identification input sequence was independently generated. Normally distributed, zero-mean measurement noise with a standard deviation σ ) 0.1 was added to each output. The Matlab System Identification Toolbox25 was used to identify ARX and FIR models for each dataset. The ARX models were identified under the assumption of known time delays. The FIR models had 90 parameters for each input-output pair, and no time delay structure was assumed. Following identification, a short set of input data (200-1000 min) for validation was applied to each identified model. Validation data were independently generated using the same frequency characteristics as for the identification sequences. The four types of input sequences used for identification were equally represented in the validation data. Identification results for each type of input excitation were averaged over 100 simulation runs. 4.3. Results. For each of the four types of input excitation, two 90th-order MISO FIR models were identified. The sampling period was 1 min, and a 90thorder model was chosen to correspond to the settling time of the process (4-5 times the largest time constant). Consequently, a total of 24 MISO FIR models were identified in each set of simulations. One hundred runs with independently generated input and noise sequences were performed to assess the variability of the identification algorithm. Table 2 shows R2 values for the identified FIR models and the validation data. The R2 values are the mean values for the 100 runs. The R2 values for the y2 outputs and their standard deviations are also included. Because the standard deviations are very large compared to the R2 values, there is considerable variability for the identified models. It appears that using simultaneous and rotated input sequences can lead to a small improvement in R2 values for PRBS input sequences, but the standard deviations indicate that this improvement is insignificant compared to the variability of the individual runs. Next, low-order MISO ARX models were identified, and corresponding R2 values were calculated. Let R(n)2 denote the R2 value for a prediction horizon of n, as

Ind. Eng. Chem. Res., Vol. 43, No. 14, 2004 3853 Table 3. R2 (%) Values for Identified ARX Models and Validation Data n ) 1 n ) 20 n ) ∞

sequential simultaneous rotated sequential simultaneous rotated sequential simultaneous rotated a

JARX K

JFIR K

JARX λ

JFIR λ


45 ( 28 34 ( 18 35 ( 21

PRBS 1 130 ( 61 90 ( 45 83 ( 48

19 ( 21 18 ( 21 6.4 ( 5.9

190 ( 700 150 ( 460 28 ( 27


35 ( 17 24 ( 11 29 ( 17

PRBS 2 77 ( 36 67 ( 29 68 ( 42

27 ( 56 18 ( 28 6.6 ( 5.0

260 ( 1300 110 ( 350 28 ( 87


26 ( 13 21 ( 10 21 ( 13

PRBS 3 59 ( 28 43 ( 29 42 ( 24

22 ( 28 15 ( 17 5.7 ( 4.4

160 ( 510 150 ( 530 12 ( 9


19 ( 9 17 ( 8 18 ( 8

MPC Input 39 ( 17 32 ( 14 29 ( 15

21 ( 31 14 ( 17 5.5 ( 4.5

120 ( 430 140 ( 550 8.8 ( 6.3

y2a

y1


Table 5. Steady-State Gain and RGA Estimation Error (%)

n)1

n ) 20

n)∞

>99 >99 >99

PRBS 1 85 64 >99 ( 0 91 ( 6 86 65 >99 ( 0 91 ( 6 86 66 >99 ( 0 89 ( 7

>99 >99 >99

PRBS 2 86 67 >99 ( 0 91 ( 6 86 67 >99 ( 0 92 ( 6 86 67 >99 ( 0 90 ( 7

>99 >99 >99

PRBS 3 86 65 >99 ( 0 91 ( 6 86 66 >99 ( 0 92 ( 6 86 67 >99 ( 0 90 ( 7

>99 >99 >99

MPC Input 85 63 >99 ( 0 87 ( 14 79 ( 26 86 64 >99 ( 0 87 ( 14 79 ( 26 84 62 >99 ( 0 89 ( 11 83 ( 20

86 ( 13 86 ( 11 82 ( 16 86 ( 11 87 ( 11 83 ( 15 86 ( 11 87 ( 11 83 ( 15

Standard deviations are included for the y2 results.

Table 4. Improvements in Parametric Accuracy and Test Duration Obtained Using Simultaneous Excitation: ARX Models input excitation

Ip (%)

reduction in test duration (%)

PRBS 1 rotated PRBS 2 rotated PRBS 3 rotated MPC input rotated

9.6 22 8.7 19 16 42 22 57

17 17 18 18 26 26 31 31

average

24

23

defined in eq 27. Values were calculated for n ) 1, 20, and ∞. Only R(1)2 was calculated for the FIR case because there are no output terms for a FIR model. The n ) 20 and n ) ∞ results are included to provide some insight into the steady-state accuracy of the identified ARX models. Table 3 summarizes the results for the ARX models. A comparison of Tables 2 and 3 indicates that the ARX models were identified much more accurately than the FIR models. This general result is expected, because the FIR models contain many more parameters. The onestep-ahead predictions for the identified ARX models are extremely accurate. As was the case for the FIR models, the ARX models identified using simultaneous input excitation have slightly larger R2 values than those identified using sequential input excitation. These improvements are very small compared to the standard deviations. Note that the standard deviations of the R2 values are much smaller for identified ARX models than for identified FIR models, because of the lower variance achieved by using a lower-order model that has the correct model structure. The improvement in parametric accuracy was defined in eq 28 using sequential excitation as the benchmark for the comparison. For ARX models, the improvements provided by simultaneous excitation or rotated input sequences are reported in Table 4. The first column shows the improvement in parametric accuracy. The second column shows the percent reduction in test duration that can be achieved using simultaneous excitation to give a model that is equivalent to the model identified using sequential input manipulation (that is,

Ip ) 0%). The trends in parametric accuracy improvement are very clear. Simultaneous input excitation leads to more accurate models. When rotated input designs are used for identification, the most accurate models are estimated. It is also clear that PRBS signals with longer switching times and the representative MPC identification test signal lead to more accurate models than PRBS signals with short switching times. Two other measures of model parameter accuracy are presented in Table 5. The first is the normalized steadystate gain estimation error, JK, defined in eq 30. Improvement in gain estimation error occurs when either simultaneous excitation or rotated inputs are employed. Although this improvement occurs for both ARX and FIR models, it is more pronounced for the FIR models because of the large number of parameters that were identified. However, the gain estimation is much more accurate for ARX models. The measure of RGA accuracy, Jλ, defined in eq 31 and reported in Table 5, show that improvements similar to the increase in steady-state gain accuracy occur for simultaneous excitation. However, the improvement in Jλ is much more dramatic when rotated input designs are used. One possible explanation is the fact that rotated input designs provide equal amounts of excitation in both the strong and weak gain directions. As a result, the weak gain direction is estimated much more accurately for rotated inputs than for the other two types of input excitation. Thus, whereas the overall gain estimation error, JK, obtained using rotated inputs is comparable to that found in the simultaneous excitation case, the RGA is estimated much more accurately because the simultaneous excitation case primarily increases the accuracy of the strong gain direction. Rotated input excitation increases accuracy in both the strong and weak gain directions. This trend is observed for both the ARX and FIR models. The improvement is even more striking for the FIR case because the steady-state gain accuracy is quite poor using standard input excitation. The rotated input designs provide a more accurate estimation of gains in both the strong and weak directions of the process. Figure 6 provides another demonstration of how simultaneous input excitation can lead to lower model variability. In Figure 6, the amplitude ratio (AR) portions of the Bode plots for the 100 identification runs are shown. The models identified using sequential excitation have a larger variability than the models


Figure 6. Amplitude ratio plots of identified models.

identified using simultaneous excitation. These results are consistent with the standard deviations of the steady-state gain estimation error reported in Table 5. 5. Conclusions Comparisons of several methods of input excitation for process identification demonstrate that a substantial improvement in parametric accuracy can be achieved for ARX models by using simultaneous and rotated input designs. Simultaneous input excitation also reduces the variability of the identified models. For highorder FIR models, the model variability is so large that meaningful comparisons are not possible. The R2 values for the prediction error results show that little, if any, improvement is gained by using simultaneous input manipulation. Improvement in the steady-state gain estimation error JK for ARX and FIR models occurs for both simultaneous excitation and rotated input excitation. Finally, it was shown that rotated input designs lead to much more accurate estimates of the relative gain array for this case study. The reader is cautioned against extending these conclusions to nonlinear case studies. The results obtained from this case study are believed to be applicable to a wide range of linear systems. Literature Cited (1) Rivera, D. E.; Jun, K. S. An Integrated Identification and Control Design Methodology for Multivariable Process System Applications. IEEE Control Syst. Mag. 2000, 20 (3), 25-36. (2) Godfrey, K. Perturbation Signals for System Identification; Prentice Hall: New York, 1993. (3) Murakami, K.; Seborg, D. E. Constrained Parameter Estimation with Applications to Blending Operations. J. Process Control 2000, 10, 195-202. (4) Qin, S. J.; Badgwell, T. A. A Survey of Industrial Model Predictive Control Technology. Control Eng. Pract. 2003, 11, 733764.

(5) Andersen, H. W.; Ku¨mmel, M. Identifying Gain Directionality of Multivariable Processes. In Proceedings of the European Control Conference; He´rmes, Paris, 1991. (6) Koung, C. W.; MacGregor, J. F. Design of Identification Experiments for Robust Control. A Geometric Approach for Bivariate Processes. Ind. Eng. Chem. Res. 1993, 32, 1658-1666. (7) Zhan, Q.; Georgakis, C. Steady State Optimal Test-Signal Design for Constrained Multivariable System. In IFAC SYSID 2000: Symposium on System Identification; Elsevier Science: New York, 2000. (8) Seborg, D. E., Edgar, T. F.; Mellichamp, D. A. Process Dynamics and Control, 2nd ed.; Wiley: New York, 2004; Chapter 20. (9) Ljung, L. System Identification Theory for the User, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, 1999. (10) van Overschee, P.; De Moor, B. Subspace Identification for Linear Systems; Kluwer Academic Publishers: Norwell, MA, 1996. (11) Chien, H.; Ogunnaike, B. A. Modeling and Control of a Temperature-Based High-Purity Distillation Column. Chem. Eng. Commun. 1997, 158, 71-105. (12) Zhu, Y. Multivariable Process Identification for MPC: the Asymptotic Method and its Applications. J. Process Control 1998, 8, 101-115. (13) Zhu, Y.; Butoyi, F. Case Studies on Closed-Loop Identification for MPC. Control Eng. Pract. 2002, 10, 403-417. (14) Freudenberg, J. S.; Middleton, R. H. Scaling and Redundancy for Ill-Conditioned Two Input, Two Output Plants. Automatica 2002, 38, 499-505. (15) Andersen, H. W.; Ku¨mmel, M. Evaluating Estimation of Gain Directionality, Part 1: Methodology. J. Process Control 1992, 2, 59-66. (16) Wang, L.; Cluett, W. R. Use of PRESS Residuals in Dynamic System Identification. Automatica 1996, 32, 781-784. (17) Montgomery, D. C.; Runger, G. C. Applied Statistics and Probability for Engineers, 3rd ed.; Wiley: New York, 2002. (18) Chen, D.; Seborg, D. E. Relative Gain Array Analysis for Uncertain Process Models. AIChE J. 2002, 48, 302-310. (19) Wood, R. K.; Berry, M. W. Terminal Composition Control of a Binary Distillation Column. Chem. Eng. Sci. 1973, 28, 17071717. (20) Chen, D.; Seborg, D. E. Design of Decentralized PI Control Systems Based on Nyquist Stability Analysis. J. Process Control 2003, 13, 27-39. (21) Åstro¨m, K. J., Johansson, K. H.; Wang, Q. G. Design of Decoupled PI Controllers for Two-by-two Systems. IEE Process Control Theory Appl. 2002, 149, 74-81. (22) McAvoy, T. J. Connection Between Relative Gain and Control Loop Stability and Design. AIChE J. 1981, 27, 613. (23) Gaikwad, S. V.; Rivera, D. E. Control-Relevant Input Signal Design for Multivariable System Identification: Application to High-Purity Distillation. In IFAC 13th World Congress; Elsevier Science: New York, 1996; pp 349-354. (24) Rivera, D. E. Monitoring Tools for PRBS Testing in ClosedLoop System Identification. In Proceedings of the 1992 AIChE Annual Meeting; American Institute of Chemical Engineers (AIChE): New York, 1992; p 131d. (25) Ljung, L. System Identification Toolbox for Use with MATLAB; The MathWorks, Inc.: Natick, MA, 2000.

Received for review August 14, 2003 Revised manuscript received January 30, 2004 Accepted February 3, 2004 IE034068+

An Evaluation of MIMO Input Designs for Process Identification

Recommend Documents