Steady State Optimal Test Signal Design for Multivariable Model

In current model predictive control (MPC) applications, the accuracy of the models utilized is usually the main factor for the final successes. These ...
0 downloads 0 Views 203KB Size
Download PDF
8514

Ind. Eng. Chem. Res. 2006, 45, 8514-8527

Steady State Optimal Test Signal Design for Multivariable Model Based Control Qiang Zhan and Tong Li Chemical Process Modeling and Control Center, Department of Chemical Engineering, Lehigh UniVersity, Bethlehem, PennsylVania 18015

Christos Georgakis* Department of Chemical and Biological Engineering and Systems Research Institute, Tufts UniVersity, Medford, Massachusetts 02155

In current model predictive control (MPC) applications, the accuracy of the models utilized is usually the main factor for the final successes. These models, especially for complex multivariable processes, are often obtained by system identification. Then, the design of the identification input signals is the key step for a successful MPC scheme. In this paper, a new design method is proposed based on a preliminary steady state gain model of a multivariable system and process operability. This is a general design method, and the existing methods can be viewed as its special cases. Previous publications have only considered an orthogonally rotated input signal whose aspect ratios are fixed by the singular values of the gain matrix. The one presented here also allows oblique projections, and the aspect ratios are calculated by the use of explicit considerations on the outputs. As a result, the identified model from this method is expected to be more accurate. From the parameter covariance analysis, it is shown that this method leads to the D-optimal design for the gain matrix parameter estimation. Because the preliminary model is usually an approximate one, model uncertainty is explicitly accounted for in the proposed method. Two examples are given to prove the effectiveness of the proposed method. This paper is a more detailed account of an earlier communication [Zhan, Q.; Georgakis, C. In Proceedings of the 12th IFAC Symposium on System Identication, Santa Barbara, CA, 2000]. 1. Introduction In the last two decades, model predictive control (MPC) has emerged as a powerful practical control technique. Various MPC techniques such as dynamic matrix control (DMC),2 model algorithmic control (MAC),3 internal model control (IMC),4 and generalized predictive control (GPC)5,6 were proposed and used in both theory and practice. Significant profit was claimed by implementing these kinds of technology. However, it is also very well-known that the success of MPC relies heavily on the model it utilizes. The modeling of a system is by far the most time and money consuming part in the MPC project. Interaction between methods of identification and control is a recent topic of interest. It is well-accepted that the control objective should be taken into account at the identification step to yield a “control-relevant model”. Once data have been collected, there is little one can do but try to minimize bias and covariance errors through the proper selection of the data prefilter and noise model. On the other hand, if data have not yet been collected, it is the design of experiments that will have the greatest impact on identified models. A good identification test plays a dominant role in system identification. No identification algorithm can give a good model from data containing poor information. It is necessary to know which characteristics of the process are important to be identified accurately in order to perform an efficient control-relevant experiment design. For single-input single-output (SISO) systems, it has been shown that the crossover frequency part of the model is of utmost importance for robust controller design.7 For multi-input multi-output (MIMO) systems, it has been pointed out that precise estimation * To whom correspondence should be addressed. 4 Colby Street, Medford, MA 02155. Phone: (617) 627-2573. E-mail: [email protected].

of gain directionality is a necessary condition for the design of high performance multivariable controllers.8,9 However, as pointed out by Zhu,10 although the multivariable control concept is well accepted by the MPC industrial community, the identification practice is mostly based on single-input thinking. The conventional way for obtaining a MIMO system model is to perturb the inputs separately. The SISO models are identified first and then combined into a MIMO one.11 The current practice in the MPC industry is to use a series of PRBStype single variable tests for model identification. The biggest advantage of performing a single test for one input at a time is that it is conservative and relatively safe. Since in most cases an operator knows from process knowledge how and to what extent the step or pseudo-random-binary-sequence (PRBS) change of a certain input will affect the outputs, the plant safety and product quality can be relatively easily controlled by the operator. However, this one input at a time approach makes the identification experiments last weeks to months. Furthermore, it also contains poor information for a MIMO system because gain directionality is not considered. When SISO models are used together for model-based control, small errors in the individual SISO models may cause significant performance degradation and even instability in the MIMO system. The reason is that a certain combination of errors in the SISO models has a much more profound effect on the closed-loop controller behavior than other combinations of errors of a similar amplitude. In other words, the gain directionality of a system should be considered. In fact, the accuracy of individual transfer functions is not a good measure for assessing the quality of a model to be used for multivariable model based control. Due to the inefficiency of the one input at a time approach, it is argued that several inputs should be tested simultaneously. It is shown8,9 that independent multivariable perturbation gives better gain directionality information. The independent multi-

10.1021/ie0601201 CCC: $33.50 © 2006 American Chemical Society Published on Web 08/23/2006

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8515

Figure 1. Process steady state operability.

variable PRBS perturbation has also been reported for a real industrial MPC project.10 A significant benefit was claimed both in the testing time and the resulting model accuracy. The open-loop excitation for MIMO system identification was studied by Koung and MacGregor.12,13 Singular value decomposition (SVD) is performed on the steady state gain model of the system, and the input signal is designed so that all the output directions will be excited equally well at steady state. For illconditioned systems, the amplitudes of the input signal in the weak directions will be much higher than the amplitudes in the strong directions. It has been proven that this design method is optimal for the identification of the input rotation matrix and the minimum singular value, which are very important for the robust multivariable model based controller design. However, some limitations can be found. First, the accurate steady state gain model is needed in order to perform the optimal input design. When the steady state gain model is not accurate, there is no systematic way of considering the model uncertainty. Second, there is no systematic way of incorporating the input and output constraints so as to design the signal with the maximal signal-to-noise ratio. An identification method to fit both a process model and its inverse for an ill-conditioned multivariable process was proposed by Li and Lee.14 Achieving a small error in the gain matrix guarantees the accuracy of the large singular values and vectors, and a small error in the inverse gain matrix guarantees the accuracy of the small singular values and vectors. The model found through this method will have accurate gain directionality in both the strong and weak directions. Both open-loop and closed-loop experiments are needed. However, the formidable task of controller design and a large number of experiments for high dimensional problems make it impractical. The closedloop experiment is also recommended8,9,12,13 to guarantee the gain directionality. Although the closed-loop experiment is desired, there are also certain disadvantages. First, for a highdimensional ill-conditioned system, the design of the stabilizing controller itself is not a trivial task. Second, in closed-loop system identification, any error in the noise model will be brought into the process model and cause additional bias error compared to the open-loop system identification. In this paper, only open-loop experiment design will be considered. A new multivariable steady state test signal design method is proposed. The physical meaning of this design method

will be given by using the process operability concept proposed by Vinson and Georgakis.15,16 It will be shown that the existing open-loop testing methods can be unified under this design framework. 2. Preliminaries 2.1. Process Operability. Vinson and Georgakis16 proposed a new measure of process operability. Several important new definitions were introduced. The available input space (AIS) is defined as the set of values that input variables can take based on the design of the process. The desired output space (DOS) represents the set of desired values of the output variables. In addition to the desired values of the output variables, one is also interested in the output values that can be achieved based on the available input values and the process characteristics described by the steady state process model. The set of all such achievable output values is denoted by the achievable output space (AOS) and can be calculated by AOS ) G(AIS), where G represents the steady state model. The desired input space (DIS) refers to the values of the input variables that enable us to achieve the desired output values. The DIS can be calculated by DIS ) G-1(DOS). All the above definitions can be illustrated by a simple example. Consider a 2 × 2 system with steady state gain matrix G as follows:

G)

[

0.55 -0.45 0.45 -0.55

]

(1)

Suppose that the input variables u1 and u2 can both be changed within the range (4, and both the output variables y1 and y2 are not expected to be out of the range (1. Then the AIS, AOS, DIS, DOS are shown in Figure 1. The operability concepts introduced above are very useful in interpreting the input test signal design. The purpose of introducing the test signal is to excite the system enough to obtain experimental data containing rich information about the system. On the other hand, it is desirable not to excite the process too much since the effect on the plant operation should be minimized. Some tradeoff has to be made on the extent to which the system can be excited. Generally speaking, the input signal should be designed to be within both the AIS and the DIS. It needs to be inside the AIS because the AIS represents

8516

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

Figure 2. Process steady state operability subject to output directionality change.

Figure 3. Process steady state operability subject to input directionality change.

physical limitations of the input moves. And it also needs to be within the DIS because the output is expected to be within the DOS when the input excitation signal is injected. 2.2. Gain Directionality. SVD is an important tool for the analysis of MIMO systems, as the singular values represent a multivariable generalization of the classical concept of gain for SISO systems. A system is called “ill-conditioned” when σ1 . σn, where σ1‚‚‚σn are its singular values. The ith direction defined by singular vectors wi and Vi is called a “strong direction” if the corresponding σi is large and is called a “weak direction” if the corresponding σi is small. For an ill-conditioned process, the output will exhibit strong directionality. For simplicity, consider a 2 × 2 system with a steady state gain matrix G. We have the following input and output relationship:

y ) Gu ) (WΣVT)u ) WΣ(VTu) ) WΣz ) σ 0 z [w1 w2 ] 1 σ z1 ) (σ1z1)w1 + (σ2z2)w2 (2) 0 2 2

[ ][ ]

It can be seen that the output vector y is the combination of two directions w1 and w2 weighted by σ1z1 and σ2z2, respectively. Thus, if the system is ill-conditioned and hence σ1 . σ2, the

output vector will be completely dominated by w1 and exhibit strong directionality in w1 if z1 and z2 have the same magnitude. Taking the same 2 × 2 steady state gain matrix G as in eq 1 and performing SVD on G, we will get the following:

G)

[

]

0.55 -0.45 ) 0.45 -0.55 0.7071 -0.7071 1 0 0.7071 -0.7071 0.7071 0.7071 0 0.1 -0.7071 -0.7071

[

][ ][

]

(3)

The ratio of the two singular values is 10, so the output will exhibit the strong direction in w1 ) [0.7071 0.7071]T, or [1 1]T, and the weak direction in w2 ) [-0.7071 0.7071]T, or [-1 1]T. From Figure 1, it can be seen that the magnitude of the AOS in the w1 direction is about 10 times bigger than that in the w2 direction. For ill-conditioned systems, the output directions are mainly affected by the output rotation matrix W in most situations. The input rotation matrix V has a lesser effect. From the operability point of view, the directions of the AOS are mainly decided by the W matrix for ill-conditioned systems, while the V matrix has little effect on the directions of the AOS. On the other hand, the directions of the DIS are mainly decided by the V matrix while W plays a less important role on it. To

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8517

illustrate this phenomenon, define GR as the following:

GR )

[

1

xR

2

+1

][ ][

R -1 1 0 0.7071 -0.7071 1 R 0 0.1 -0.7071 -0.7071

]

(4)

The value of R will be changed from -10 to 0 to 10 at the step change of 1. That is, we allow the output matrix to change in all the directions while the singular values and the input rotation matrix are fixed. The operability region is shown in Figure 2 with the same AIS and DOS as defined before. It can be seen that the directions of the DIS are almost the same for all the values of R, while the direction of the AOS is changing with R correspondingly. Similarly, the Gβ matrix can be defined as follows:

Gβ )

1

xβ

2

+1

[

][ ][

0.707 -0.707 1 0 β -1 0.707 0.707 0 0.1 -1 -β

]

(5)

The value of β will be changed from -10 to 0 to 10 at the step change of 1. That is, we allow the input rotation matrix to change in all directions while the singular values and the output rotation matrix are fixed. The operability region is shown in Figure 3 with the same AIS and DOS as used before. It can be seen that the directions of the DIS are changing with β correspondingly, while there is almost no change in the AOS. In most situations, the input rotation matrix V has little effect on the AOS directionality, and the output rotation matrix W has little effect on the DIS directionality. 2.3. Robust Models for Designing Control Systems. Classical SISO identification does not necessarily provide adequate models to the design of control systems for ill-conditioned MIMO processes. Consider the multivariable control system in Figure 4, the transfer function matrix between the disturbances and the outputs is

Y(s) ) [I + G(s)C(s)]-1D(s)

(6)

When the process model G ˆ (s) is used for designing the control system, it becomes17

Y(s) ) [I + G ˆ (s)C(s)]-1[I + E(s)H(s)]D(s)

(7)

E(s) ) G ˆ (s)G-1(s) - I

(8)

H(s) ) [I + G(s)C(s)]-1G(s)C(s)

(9)

where

At steady state, the complimentary sensitivity function H(s) is usually identity. Hence, if the term E(s) is small, then the performance of the controller on the real system G(s) will not be much different from that on the model G ˆ (s). For robust performance, a model with small E(s) ) ||G ˆ (s)G-1(s) - I|| is desired. In this paper, we will only consider the mismatch in the steady state gain model. 3. Steady State Optimal Test Signal Design 3.1. Insight from the Existing Methods. There are several ways in the literature and practice to perform open-loop input excitation for MIMO systems. Some of the methods used can be summarized as follows. • Method 1: PRBS-type signals in only one input at a time.

Figure 4. Control system.

• Method 2: Multiple independent PRBS signals simultaneously applied in the input variables. • Method 3: Multiple independent PRBS signals simultaneously applied along the orthogonally rotated directions of the original input variables. The way of rotating the input signals defines a specific subcase. In Method 3, one can include as a subcase the method proposed by Koung and MacGregor12,13 and further examined by Conner and Seborg.18 However, both of these publications examined only orthogonal rotations with the aspect ratio determined by the relative magnitudes of the singular values of the gain matrix. This method is based on the SVD of the steady state gain matrix. The steady state relationship between inputs and outputs can be expressed as

y ) Gu ) WΣVTu ) WΣz ) (σ1z1)w1 + (σ2z2)w2 + ‚‚‚ (σnzn)wn (10) They argued that the optimal input test signal should be able to excite all the output directions wi equally. Since the different output directions are weighted by different terms σizi, the rotated input z ) VTu is designed as independent PRBS signals. The amplitude of zi is chosen to be inversely proportional to the singular values, (zi/zj) ) (σj/σi), so that all the output directions are equally excited. Assume there are n independent PRBS signals (si, i ) 1, 2, ‚‚‚, n) with normalized amplitude (1. The test signal for method 1 can be expressed as the following:

[ ] [ ][ ]

u1 a1 u2 0 u) · ) · ·· ·· un 0

0 a2 ·· · 0

··· ··· ··· ···

0 s1 0 ··· 0 0 s2 ··· ·· ·· ··· ·· · · · an 0 0 ···

0 0 ·· · sn

(11)

Where, ai is the amplitude of the corresponding input signal. The testing signal for method 2 can be expressed as the following:

[ ] [ ][ ]

u1 a1 u2 0 u) · ) · ·· ·· un 0

0 a2 ·· · 0

··· ··· ··· ···

0 s1 0 s2 ·· ·· · · a n sn

(12)

Here, ai is also the amplitude of the corresponding input signal. We are interested in selecting the maximal possible values of the amplitudes ai when both the input and output constraints are satisfied. Since all the inputs are injected simultaneously, the amplitude here of each input will be smaller than the amplitude of the previous method under the same input and output constraints. For method 3, the Koung-MacGregor (K-M) subcase that corresponds to the approach proposed by Koung and

8518

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

[ ]

MacGregor12,13 and uses the VT matrix for the rotation of the inputs. The rotated input z can be expressed as the following:

[]

z1 z z ) ·2 ·· zn

[]

k σ1 0 · · · k 0 σ ··· ) 2 ·· ·· ··· · ·

0

0 ·· · k 0 0 ··· σ n

[]

max abs(det(T))

s1 s2 ·· · sn

(13)

[ ] [

]

a1 0 1 0 0 -1 ) 0 a2 0 1 -1 0 0.5a1 0.5a1 1 1 -1 -1 0.5a2 -0.5a2 1 -1 1 -1

(14)

][

[ ][

a1 0 1 1 -1 -1 0 a2 1 -1 1 -1

]

]

(15)

(16)

• Method 3 (K-M subcase):

[

][

k k 1 1 -1 -1 u ) σ V1 σ V2 1 -1 1 -1 1 2

]

(17)

From the above analysis, it is clear now that all the existing methods can be put into one unified format for the system with two inputs:

[ ] [

(1 1 1 -1 -1 )T u)T (1 1 -1 1 -1

]

] [ ]

(20)

(18)

The matrix T is the linear transformation matrix to be designed. Different prior methods correspond to the specific structures of the T matrix. It is clear here that the signal design problem can be generalized to the one where the structure of the T matrix is not prespecified. Then, one needs to find the transformation matrix T that can maximize the information contained in the input-output data. In fact, the objective used in this paper is to find the T with the maximal determinant under the input/output constraints. It will be proven later that the T with the maximal determinant will lead to the D-optimal design.11

(21)

All three previous methods can be readily transformed into the same optimization formulation. Since different methods are based on different geometric assumptions, additional constraints need to be imposed: • Additional constraints for method 1 (one input at a time)

t11 - t12 ) 0 and t21 + t22 ) 0

(22)

The additional constraints are given because the geometric shape is a rhombus and the vortices must be on the axes. • Additional constraints for method 2 (independent multivariable PRBS)

t12 ) 0 and t21 ) 0

(23)

The additional constraints are imposed due to the fact that the test input space is a rectangle with the side parallel to the axes. T is a diagonal matrix. • Additional constraints for method 3

t11/t21 ) V11/V21 and t12/t22 ) V12/V22

• Method 2:

u)

] [ ]

yl1 yh1 1 1 -1 -1 yl2 < GT 1 -1 1 -1 < yh2

As stated before, this method was originally proposed for the unconstrained system. For the constrained system, k should be given the biggest possible value which can satisfy both the input and output constraints. Conner and Seborg18 provided an initial suggestion for such a correction by the use of the largest singular value of G(0). 3.2. Steady State Test Signal Design. In this section, we first examine the case of 2 × 2 systems. However, the result is easily expandable to the general m × n systems. It can be found that all the existing three methods can be put into one unified design format. Since every unit-amplitude PRBS signal has only two levels of 1 and -1, the possible input signals of the existing methods can be expressed as follows. • Method 1:

[

[ ] [ DOS constraints

[]

[

]

AIS constraints

ul1 uh1 1 1 -1 -1 ul2 < T 1 -1 1 -1 < uh2

u1 s1 k k k u2 s u ) · ) Vz ) σ V1 σ V2 · · · σ Vn ·2 ·· ·· 1 2 n un sn

[ ][

(19)

tij

and the resulting input signal can be expressed as

u)

The input signal design problem can be formulated as the following optimization problem: Objective

(24)

Where, Vij is the element of the input rotation matrix V of the nominal model G ˆ . The additional constraints are imposed due to the fact that the test input space is a rectangle with the sides parallel to the V matrix directions. In the original formulation of this method,13 the amplitude of Vij needs to be adjusted according to the ratio of the singular values. 3.3. Geometric Interpretation. The determinant of a transformation matrix is related to the change of the volume of a region in the initial space (x) when it is transformed into the corresponding region in the new space (y). A linear transformation T expands (or shrinks) the volume of the region by |det(T)|. So, the volume of the image y ) Tx has a simple formula:

volume(y) ) |det(T)|volume(x)

(25)

Both the previous and proposed new design methods can be interpreted geometrically. All the methods try to maximize the areas of different geometric shapes inside the AIS and the DIS. The geometric shapes for different methods can be summarized as follows. • Method 1: the single input signal (SIS) region can be viewed as a rhombus in the 2 × 2 case with the diagonals parallel to the axes of the ui coordinate system. • Method 2: the multiple input signal (MIS) region can be viewed as a rectangle in the 2 × 2 case and as an orthogonal parallelepiped in higher dimensions with sides parallel to the ui coordinate system. • Method 3: the orthogonally rotated multiple input signal (ORMIS), of which the K-M method is a subcase, has the input

Ind. Eng. Chem. Res., Vol. 45, No. 25, 2006 8519

G(z) ) a0G ˆ (z), where al < a0 < au

(26)

With this model uncertainty description, it will be possible to calculate the output range for the whole model uncertainty range for certain inputs. 3.4.1. Output Range for a SISO System. It is very easy to deal with the SISO system steady state gain uncertainty. Suppose the current model identified is G ˆ (z), the real system is denoted ˆ (z). The upper bound for the gain a is denoted as as G(z) ) a0G au > 1, and the lower bound, as al < 1. Then, the real system output will always be bounded by the upper bound of the model predicted output. 3.4.2. Output Range for a MIMO System. For MIMO systems, the situation is a bit more complicated. Let us first study a 1 × 2 system

y ) y1 + y2 ) g1(z)u1 + g2(z)u2 Figure 5. Process steady state operability.

signal region that can be viewed as a rotated rectangle in the 2 × 2 case and a rotated orthogonal parallelepiped in higher dimensional cases. The rotation is defined by the input rotation matrix VT. • New Method: methods 2 and 3 (including the K-M subcase) mentioned above are generalized by the proposed method where the input signal region is an oblique parallelogram or parallelepiped. The signal is named as constrained transformed multiple input signal (CTMIS). Actually, methods 2 and 3 are special cases of the new method which is proposed in this paper. For a general 2 × 2 system with the same AIS and DOS constraints, the resulting test input space areas are illustrated in Figure 5 for different methods. The larger the input space area is, the bigger the determinant of the transformation matrix T. As can be seen, method 3 and the new method give much larger input space areas than methods 1 and 2. 3.4. Model Uncertainty. From the above design procedure, the maximal test input space with parallelogram shape can be found inside the intersection of the AIS and the DIS. In such a case, both the input and output constraints can be satisfied at steady state. However, the output is obtained by using the available steady state model of the process and the model cannot be realistically assumed to be 100% correct. Since the experiment is performed in order to get a good model, the accurate model is never available at the input design stage. Although the predicted output will satisfy the constraints based on the nominal model, the actual output may not always satisfy the constraints. Then, the model uncertainty must be incorporated into the design. The dynamic model uncertainty is very difficult to be dealt with. The popular description of model uncertainty is in the frequency domain. However, there is no existing methods to obtain the time domain boundary of the output under the frequency domain model uncertainty description. To simplify the problem, the assumption made here is that the model uncertainty only exists at the steady state. The model uncertainty will be described by a upper-bound parameter au > 1 and a lower-bound parameter 0 < al < 1. Suppose the current dynamic nominal model is G ˆ (z) with model uncertainty au and al; then, the real model G(z) is assumed to belong to the model uncertainty range:

(27)

The currently identified models are denoted as gˆ 1(z) and gˆ 2(z). Define the upper and lower bounds as au > 1, bu > 1, 0 < al < 1, and 0 < bl < 1. Then, the boundary of y(k) can be shown in Table 1 for different situation of y1 and y2. To satisfy the real output constraints by using the current model, the worst case scenario should be considered. For the 1 × 2 system above, at time k, to satisfy the constraints

yl < y(k) < yh

(28)

[ ] [ ][ ] [ ]

we will need to satisfy the following

yl(k) yh(k) au bu yl(k) y (k) au bl yˆ 1(k) e a b e h (k) yl(k) y ˆ yh(k) l u 2 a b yl(k) yh(k) l l

(29)

The above formulation can be very easily generalized for m × n MIMO systems. 3.4.3. Steady State Input Signal Design under Model Uncertainty. For a 2 × 2 system, suppose the current nominal steady state model with model uncertainty description is given by the following:

G ˆ )

[

(a11l ∼ a11u)gˆ 11 (a12l ∼ a12u)gˆ 12 (a21l ∼ a21u)gˆ 21 (a22l ∼ a22u)gˆ 22

]

(30)

The optimization problem can be formulated as: Objective

max abs(det(T))

(31)

tij

AIS constraints

[ ] [

] [ ]

ul1 uh1 1 1 -1 -1 ul2 < T 1 -1 -1 1 < uh2

(32)

DOS constraints

[ ][

ai1u a yil < ai1u i1l ai1l

ai2u ai2l gˆ i1 0 1 1 -1 -1 ai2u 0 gˆ i2 T 1 -1 1 -1 < yih ai2l ∀ i ) 1, 2 (33)

][

]

8520

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

Table 1. Output Boundary for the Uncertainty Model Range yˆ 1(k)

yˆ 2(k)

y(k)-upper bound

y(k)-lower bound

>0 0 0 0. And, to eliminate the effect of column interchange, we can add such a constraint as t11 < t12 < ‚‚‚ < t1n. This is important for the high dimensional problems where there are a lot of local optima. Literature Cited (1) Zhan, Q.; Georgakis, C. Steady state optimal test signal design for multivariable model based control. In Proceeding of the 12th IFAC Symposium on System Identication, Santa Barbara, CA, 2000. (2) Cutler, C. R.; Ramaker, B. L. Dynamic Matrix Control- A Computer Control Algorithm. Proceedings of the 1980 Joint Automatic Control Conference, San Francisco, CA, August 1980; WP5-B/6. (3) Rouhani, R.; Mehra, R. K. Automatica 1982, 18, 401-414. (4) Garcia, C.; Morari, M. Ind. Eng. Chem. Process Des. DeV. 1982, 21, 308-323. (5) Clarke, D. W.; Mohtadi, C.; Tuffs, P. S. Automatica 1987, 23, 137148. (6) Clarke, D. W.; Mohtadi, C.; Tuffs, P. S. Automatica 1987, 23, 149160. (7) Astrom; Wittenmark AdaptiVe Control; Addison-Wesley: Reading, MA, 1989. (8) Andersen, H. W.; Kummel, M. J. Process Control 1992, 2, 59-66. (9) Andersen, H.; Kummel, M. J. Process Control 1992, 2, 67-86. (10) Zhu, Y. J. Process Control 1998, 2, 101-115. (11) Ljung, L. System identication-theory for the user; Prentice Hall: Englewood Cliffs, NJ, 1987. (12) Koung, C. W.; MacGregor, J. F. Ind. Eng. Chem. Res. 1993, 32, 1658-1666. (13) Koung, C. W.; MacGregor, J. F. Automatica 1994, 30, 1541-1554. (14) Li, W.; Lee, J. H. Comput. Chem. Eng. 1996, 20, 1023-1042. (15) Vinson, D. R.; Georgakis, C. A new measure of process output controllability. In Proceedings of the IFAC Symposium on Dynamics and Control of Process Systems (DYCOPS-5), Corfu, Greece, June 1998; IFAC, 1999; Vol. 2, Part 2, pp 663-672. (16) Vinson, D. R.; Georgakis, C. J. Process Control 2000, 10, 185194. (17) Lee, J.; Cho, W.; Edgar, T. F. Ind. Eng. Chem. Res. 1998, 37, 10181023. (18) Conner, J. S.; Seborg, D. E. Ind. Eng. Chem. Res. 2004, 43, 38473854. (19) Goodwin, G. C.; Payne, R. L. Dynamic system identication: experiment design and data analysis; Academic Press: New York, 1977. (20) Ogunnaike, B. A.; Lemaire, J. P.; Morari, M.; Ray, W. H. AIChE J. 1983, 29, 632-640.

ReceiVed for reView January 30, 2006 ReVised manuscript receiVed July 22, 2006 Accepted July 25, 2006 IE0601201