Estimation of the Dynamic Matrix and Noise Model ... - ACS Publications

842

Ind. Eng. Chem. Res. 2002, 41, 842-852

Estimation of the Dynamic Matrix and Noise Model for Model Predictive Control Using Closed-Loop Data Ramesh Kadali and Biao Huang* Department of Chemical & Materials Engineering, University of Alberta, Edmonton, Alberta, Canada T6G 2G6

A dynamic matrix is a lower triangular matrix containing the step response coefficients of the deterministic input used in the model predictive control schemes such as the dynamic matrix controller. Subspace matrices (defined in subspace state-space identification methods) corresponding to the deterministic input and the stochastic input contain the impulse response coefficients of the deterministic and stochastic models, respectively. This paper proposes a new subspace identification based method for the estimation of the dynamic matrix of the deterministic input(s) directly from the closed-loop data. The noise model is simultaneously obtained from the closed-loop data in the impulse response form. The method is extendable to the case of measured disturbances. All of the results presented in this paper are applicable to the multivariate systems. Guidelines for the practical implementation of the algorithm are also presented in this paper. The proposed method is illustrated through MATLAB simulations and an application on a pilot-scale plant. 1. Introduction Model predictive controllers (MPC) have found many successful applications in process industries for more than 2 decades. One of the key aspects of MPC is the prediction of the future process response and minimization of the output deviation from the setpoint by manipulating the inputs. A model for the process is required to make these predictions based on the past data. Hence, an MPC design starts with first identifying a nominal model for the process. One of the industrially successful predictive control schemes is the dynamic matrix controller or DMC, which explicitly uses a lower triangular matrix called the “dynamic matrix” containing the step response coefficients corresponding to the deterministic input(s) to the process.3,4 Many other MPC formulations also use the dynamic matrix in one way or another.2,22,25 For constructing the dynamic matrix, in the case of DMC, a step response model for the process is first obtained from the open-loop data. The step response coefficients are arranged in a specific lower triangular form known as the dynamic matrix. However, for safety reasons and other practical limitations, open-loop operation of the process may not always be possible or in some cases there may be a hidden feedback in the system. Estimation of the dynamic matrix from closed-loop data is desired in such cases. It has been shown14 that if the model is used for modelbased control design, then the favorable experimental conditions are actually under a closed loop. Closed-loop identification refers to the identification of the process model and noise model using the data sampled under feedback control. Correlation between the disturbances entering the process and the input offers the fundamental limitation1,11,13,23,24,36 for utilizing the standard open-loop identification methods with * To whom all correspondence should be addressed. E-mail: [email protected]. Tel: (780) 492-9016. Fax: (780) 492-2881.

Figure 1. Experimental setup.

Figure 2. Closed-loop system.

closed-loop data (see Figure 2). Several closed-loop parametric model identification methods have been suggested in the literature which require either certain

10.1021/ie000909q CCC: $22.00 © 2002 American Chemical Society Published on Web 01/18/2002

Ind. Eng. Chem. Res., Vol. 41, No. 4, 2002 843

assumptions about the model structure or knowledge of the controller model. The closed-loop identification methods found in the literature are broadly classified11 into direct,35 indirect,6,15,38 and joint input/output16 identification methods. See refs 5, 11, and 26 for a review of the features and limitations of different closedloop identification methods. The subspace identification method is a relatively new approach used for the state-space model identification. In this approach, certain subspace matrices of the process are first identified, by regression of the data Hankel matrices (refer to the standard book on subspace identification such as that in ref 29 and references therein), from which the state-space matrices are extracted. Three subspace matrices are obtained as a first step,27-29 which correspond to the states, the deterministic input(s), and the stochastic input(s) to the system. These subspace matrices are directly calculated from the input/output data matrices in a single iteration compared to the iterative schemes used in the prediction error methods. Moreover, the subspace identification methods minimize the summation of the multistep ahead prediction errors during the estimation of the subspace matrices. This advantage is, however, lost when lower order state-space system matrices are estimated from the subspace matrices. Hence, directly using the subspace matrices is a very appealing idea for designing the predictive controllers.7,9,10,18,32-34 The subspace matrix corresponding to the deterministic input contains the impulse response coefficients (Markov parameters for multivariate processes) of the deterministic input(s) in a lower triangular form. Similarly, the subspace matrix corresponding to the stochastic input contains the impulse response coefficients/Markov parameters of the noise model. This allows the alternative approach for direct estimation of the dynamic matrix and noise model from the open-loop input/output data matrices. However, it has been shown11,23 that the open-loop subspace identification methods cannot be directly applied to the closed-loop data. Identification of the subspace matrices from closedloop data has recently received attention by a number of researchers.8-10,24,30 Van Overschee and De Moor30 proposed an N4SID (numerical subspace state-space identification) based method for closed-loop subspace identification which requires the knowledge of the first N impulse response coefficients of the controller, where N is the maximum order of the state-space model to be identified. Knowledge of the impulse response (IR) coefficients of the controller is required if one wants to identify all three subspace matrices and subsequently a state-space model for the system. Ljung and McKelvey24 presented a method for the identification of subspace matrices from closed-loop data using estimated predictors and stated that their algorithm is merely an illustration of a “feasible” method rather than the “best way” of identifying systems operating in a closed loop. MOESP (MIMO output error state-space model identification) and CVA (canonical variate analysis) approaches are also proposed for the identification of a state-space model using closed-loop data.21,37,40 In addition to the setpoint excitation, MOESP/CVA approaches use an external white noise signal addition to the controller output to make it independent of the noise. A closed-loop state-space model is first identified from the closed-loop data from which the open-loop state-space matrices are retrieved. The principal goal

of all of the above approaches is the identification of a state-space model for the system using closed-loop data. Even though the subspace identification method is used as a vehicle, the goal of the identification method from closed-loop data proposed in this paper is not the estimation of state-space system matrices {A, B, C, D, and K} but the estimation of the dynamic matrix of the process and the noise model in impulse response form. It is shown in this paper that if we want to estimate only two of the subspace matrices, i.e., only those corresponding to the deterministic and stochastic inputs, from closed-loop data, then knowledge of the controller impulse response coefficients can be avoided. We can then obtain the process dynamic matrix from the deterministic input subspace matrix and the noise model in impulse response form from the stochastic input subspace matrix. The method proposed in this paper can be considered as a nonparametric approach for closed-loop identification. Nonparametric model identification methods, although known to give less bias error because of less model structure and order limitations, could result in higher variance (due to the larger number of parameters) compared to parametric model identification methods. This is a tradeoff between bias error and variance error in process identification. The actual process is typically high-order and nonlinear, and it is difficult to represent it by a single linear parametric model. Consequently, bias error is inevitable in practice. On the other hand, it is known that the variance error can be reduced with the increased sample size.23 Therefore, depending on the application, for example, depending on the data sample size, one can choose to use the parametric or nonparametric identification method. However, nonparametric model identification methods do have some practical advantages. Consider the case when we want to identify a process model for designing and implementing a model-based predictive controller. Even if we first identify a parametric model for the process (using closed- or open-loop process data), the process model has to be eventually transferred to a nonparametric (impulse or step response models) form for designing the MPC. The question is whether this additional intermediate step will introduce an additional error and inconvenience. When it comes to industrial implementation, nonparametric MPCs have shown a considerable success rate. The idea here is to directly identify a nonparametric model and avoid choosing a “model structure”, which is unknown and a prerequisite for parametric model identification methods. The remainder of this paper is organized as follows. Section 2 gives an overview of the existing closed-loop subspace identification methods. Section 3 is the main section where the estimation of the process dynamic matrix and noise model from closed-loop data is described. Remarks on the different steps of the closedloop identification along with some guidelines for the practical implementation of the algorithm are provided in section 4. The application of the proposed method is explained using MATLAB simulations in section 5 and implementation on a pilot-scale process in section 6. Conclusions from the above work are given in section 7. 2. Subspace Identification Methods A linear time-invariant system can be described by an innovation form of state-space representation as

844

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

xk+1 ) Axk + Buk + Kek

(1)

yk ) Cxk + Duk + ek

(2)

where A, B, C, D, and K are the state-space system matrices. K is the Kalman filter, and ek is a white noise sequence with the following covariance matrix:

E[epeTq ]

) Sδpq

where S is the covariance matrix and δpq is the Dirac delta operator. Consider a system with l inputs and m outputs. Thus, A, B, C, D, and K are n × n, n × l, m × n, m × l, and n × m matrices, respectively. We have the measurements of the inputs and outputs uk and yk for k ∈ {0, 1, ..., 2i + j - 2)}. The data block Hankel matrices with i block rows and j block columns are arranged as

[

u0 u Up ) 1 ‚‚‚ ui-1

u1 u2 ‚‚‚ ui

‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚

]

uj-1 uj ; ‚‚‚ ui+j-2

[

ui u Uf ) i+1 ‚‚‚ u2i-1

ui+1 ui+2 ‚‚‚ u2i

‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚

ui+j-1 ui+j ‚‚‚ u2i+j-2

]

(3)

Xp ) [x0, x1, ..., xj-1]; Xf ) [xi, xi+1, ..., xi+j-1] The input-output matrix equations relating the data Hankel matrices are expressed29 as

Xf ) AiXp + ∆iUp + ∆si Ep

(4)

Yp ) ΓiXp + HiUp + Hsi Ep

(5)

Yf ) ΓiXf + HiUf + Hsi Ef

(6)

where Γi is the extended observability matrix, Hi and Hsi are the lower triangular matrices containing the impulse response coefficients/Markov parameters of the system due to the deterministic inputs uk and the unknown stochastic inputs ek, respectively. p and f denote the past and future. The subscript i follows from the number of row blocks in the block Hankel matrices as shown above. The system-related matrices are expressed as By substituting eq 5 in eq 4, we can write

[]

∆i ) [Ai-1B, Ai-2B, ..., B];

[

]

∆si ) [Ai-1K, Ai-2K, ..., K]; D 0 ‚‚‚ 0 CB D ‚‚‚ 0 Hi ) ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ; CAi-1B CAi-3B ‚‚‚ D Im 0 Im CK Hsi ) ‚‚‚ ‚‚‚ i-2 CA K CAi-3K

[

‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚

0 0 ‚‚‚ Im

]

where the superscript dagger represents the MoorePenrose pseudo-inverse. Therefore, we can write

Y ˆ f ) ΓiXf + HiUf ) LwWp + HiUf

(8)

where Y ˆ f represents the predictions for Yf. Lw, Hi, and Hsi are the subspace matrices corresponding to the states, deterministic input(s), and stochastic input(s), Y respectively, and Wp ) Up . The differences between p the different open-loop subspace-based state-space identification methods arise from (1) the way the extended observability matrix, Γi, and the future states sequence, Xf, are estimated from (LwWp) and (2) the determination of the state order. For example, N4SID and MOESP take the number of dominant singular values of the matrix (LwWp) as the optimal state order,9,27,29 while CVA uses the AIC-based method21 to determine the optimal state order. User-defined pre- and postmultiplicative weighting matrices can be used on (LwWp) for deriving Γi and Xf. It has been shown28,29 that the difference between the subspace identification methods N4SID, MOESP, and CVA is in terms of the choice of the weighting matrices used in eq 8 in deriving the states and extended observability matrix from (LwWp). However, as will be shown in the next section, we are interested in identifying only the subspace matrices. If (i) the deterministic input uk is uncorrelated with ek and ek is not identically zero, (ii) uk is persistently exciting of the order 2i, and (iii) the measurements go to infinity, j f ∞, the open-loop subspace identification9,27,29 involves finding the prediction of the future outputs Yf using a linear predictor:

[ ]

Each block in the data Hankel matrices contains a column vector of the respective variables. Similarly, the block Hankel matrices Yp, Yf, Ep, and Ef are defined for yk and ek, respectively. The past and future state sequences are defined as

Yp Xf ) [AiΓ†i , (∆i - AiΓ†i Hi), (∆si - AiΓ†i Hsi )] Up Ep

Γi ) [CT, (CA)T, ..., (CAi-1)T]T;

(7)

Y ˆ f ) LwWp + HiUf The least-squares prediction Y ˆ f can be found by solving the least-squares problem:

( )

W min ||Yf - (Lw, Hi) U p ||F2 Lw,Hi f

(9)

Y ˆ f is found by the orthogonal projection of the row space of Yf into the row space spanned by Wp and Uf defined as29

[ ]

W [Lw, Hi] ) Yf U p f

†

[ ] ([ ]

W Y ˆ f ) Y f/ U p f

)

W ) Yf[WTp , UTf ] U p [WTp , UTf ] f

-1

(10)


This projection can be implemented in a numerically robust way with a QR decomposition.27-29,31,39 However, if closed-loop data are used, the above method of estimation is not suitable because of the correlation between Uf and Ef in eq 8. Several approaches have been used to overcome this constraint. The various closedloop subspace identification approaches reported in the literature can be summarized as follows. N4SID Approach.30 The key feature of this approach is extracting open-loop subspace matrices from the closed-loop subspace matrices using the knowledge of the first “N” impulse response coefficients/Markov parameters of the controller, where N is the number of rows of the data Hankel matrices. Closed-loop data are obtained from the process by persistent setpoint excitation. The closed-loop subspace matrices are estimated using the least-squares estimation method. Open-loop subspace matrices are then extracted from the closedloop subspace matrices. A state-space model for the process is subsequently identified from the open-loop subspace matrices. The limitation of this method is the requirement of the knowledge of the controller. In industrial systems, the accurate knowledge of the IR coefficients/Markov parameters of the controller may not always be available. MOESP Approach.40 Consider a controller, C, acting on a process, P. Apart from setpoint excitation, this approach relies on excitation in the controller output with external signals to make the process input independent of noise and does not require the knowledge of the controller. r1(k) is the white noise sequence added to the controller output, u(k). r2(k) is the setpoint, which is a white noise sequence. w(k) is the white noise (disturbance) added to the process output, yt, through a noise model, Fn. Then the measurable input vector is [r1(k), r2(k)] and measurable output vector is [e1(k), e2(k), u(k), y(k)], where e1 ) Ce2 ) u - r1 and e2 ) r2 - y. Using this information, a global state-space model is first identified using the MOESP technique. The global state-space model is denoted as

[ ]

r (k) x(k+1) ) Ax(k) + [B1, B2] 1 r2(k)

[ ] [ ] [ ][ e1(k) C1 C e2(k) ) C2 x(k) + u(k) 3 C4 y(k)

D11 D21 D31 D41

D12 D22 r1(k) D32 r2(k) D42

(11)

]

(12)

P ) [A, B1, C4, D41][A, B1, C3, D31]-1 ) [A, B1, C4, D41][A - B1D31-1C3, B1D31-1,

)

[(

)(

)

[] [

][ ]

r1 (I + CP)-1 C(I + CP)-1 u ) y P(I + CP)-1 CP(I + CP)-1 r2

-D31-1C3, D31-1]

]

(C4 -D41D31-1C3), D41D31-1 (13)

(14)

Note that all of the transfer function blocks have the same denominator. Hence, to preserve the common denominator term, the CVA subspace identification method is used to first identify a state-space model, which is then converted into the transfer function form. The process transfer function, P, is then extracted by any of the following inversions:

P ) CP(I + CP)-1[C(I + CP)-1]-1 -1

-1 -1

) P(I + CP) [(I + CP) ]

(15) (16)

Ljung and McKelvey’s Approach via Estimated Output Predictors.24 First an (na, nb)-order ARX model is identified from the closed-loop input and output data. The ARX model is used to calculate the j-step ahead predictors yˆ (t + j|t) from data by replacing “future” u(t) in this prediction with zero. The vector of j-step ahead predictions is formed as

[

yˆ (t + 1|t) Ym(t) ) ‚‚‚ yˆ (t + m|t)

]

(17)

This approach first estimates {Ym(t)} using the closedloop ARX model. The states are obtained using the projection

x(t) ) LYm(t)

x(k) has an order np + nc, where np is the order of P and nc is the order of C. The state-space representations for the deterministic relations between r1 and y and between r1 and u are given by the system matrices [A, B1, C4, D41] and [A, B1, C3, D31], respectively. Using the rules for concatenating and inverting of state-space models,12 the state-space model for the process, P, is obtained as

B1D31-1 A -B1D31-1C3 , , B1D31-1 0 A - B1D31-1D3

Note that P can be identified in multiple ways using the same rules of concatenation and inversion, for example, using the state-space models between {r2&y} and {r2&u}. The overall deterministic state-space model is first identified using the MOESP algorithm. The order selection in this step needs the specification of the sum of the process model order and controller order. The individual plant model order is determined/selected in the subsequent model reduction step. CVA Approach.20,37 The CVA approach first identifies the following closed-loop transfer function relationships:

(18)

The state-space system matrices are then identified from the vectors of estimated states and the outputs. The authors state that this method is only a “feasible” method rather than the “best way” of identifying systems operating in a closed loop. Comparison of the N4SID approach and MOESP/CVA approach with the new approach, which will be explored in the next section, is presented in Figure 3. 3. Estimation of the Process Dynamic Matrix and the Noise Model from Closed-Loop Data Consider the case when the system (1)-(2) is operating under a closed loop with a linear time-invariant feedback-only controller Q, expressed in transfer function form as

uk ) Q(rk - yk)

(19)

where rk is the setpoint for the process output at the sampling instant k and rk - yk is the output deviation from the setpoint. Assume that the controller does not cancel any plant dynamics. The control system can be

846


Figure 3. Comparing the existing closed-loop subspace state-space identification methods and the new approach.

expressed in a state-space representation as

Using eqs 6 and 24, we can derive

c x(k+1) ) Acxck + Bc(rk - yk)

uk )

Ccxck

(20)

+ Dc(rk - yk)

(21)

By recursively using the above state-space equations, we can write the input-output equations for the control system as

Up ) Γci Xcp + Hci (Rp - Yp)

(22)

Xcf ) Aci Xcp + ∆ci (Rp - Yp)

(23)

Uf ) Γci Xcf + Hci (Rf - Yf)

(24)

[] [ [

][ ] ]

Yf (I + HiHci )-1Γi (I + HiHci )-1HiΓci Xf ) + -1 c c Uf Xcf -(I + Hi Hi) Hi Γi (I + Hci Hi)-1Γci (I + (I +

HiHci )-1HiHci Hci Hi)-1Hci

] [ Rf +

(I + HiHci )-1Hsi -(I + Hci Hi)-1Hci

Hsi

Ef (26)

Theorem 1. The input-output equations for the closed-loop system in eq 26 can be equivalently expressed as follows:

[]

[ ] [ ]

Yf LO LO c e ) L D + R + E p p f Uf LIc LIe f

(27)

where

Dp ) [YTp , UTp , RTp , ETp ]T; LIc ) (I + Hci Hi)-1Hci ; c -1 c LO c ) (I + HiHi ) HiHi (28)

where

c -1 s LIe ) -(I + Hci Hi)-1Hci Hsi ; LO e ) (I + HiHi ) Hi (29)

T Γci ) [CTc , (CcAc)T, ..., (CcAi-1 c ) ]; i-2 c c c c ∆ci ) [Ai-1 c Bc, Ac Bc, ..., Bc]; Xp ) [x0, x1, ..., xj-1];

[

Dc CcBc Hci ) ‚‚‚ CcAi-2 c Bc

c , ..., xi+j-1] Xcf ) [xci , xj+1

0 Dc ‚‚‚ CcAi-3 c Bc

‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚

0 0 ‚‚‚ Dc

]

[] [ ] [ ]

Y ˆf LO LO c ) M + p I Rf I U ˆf L L c

The matrices Rp and Rf are constructed in the same way as that shown in eq 3. Using eq 22 in eq 23, we can derive

[

U Xcf ) [Aic(Γci )†, (∆ci - Aic(Γci )†Hci )] R p - Y p p

Proof. See the appendix. With the above theorem, the estimation of the closedloop subspace matrices using the closed-loop data is essentially an open-loop identification problem. We can define

]

(25)

where Mp ) [YTp , UTp , RTp ]T. From the above equation and because Rf can be chosen as white noise and is uncorrelated with Mp and Ef, the closed-loop subspace matrices {LI, LIc, LO, LO c } can be obtained as a solution ˆ f} are of the least-squares estimation problem. {U ˆ f and Y found by the orthogonal projection of the row space of


{Uf and Yf} into the row space spanned by {Mp and Rf}.

[ ]

M [LI, LIc] ) Uf R p f

†

[ ]

Mp [LO, LO c ] ) Yf R f

†

([ ]

)

([ ]

)

M ) Uf[MTp , RTf ] R p [MTp , RTf ] f

M ) Yf[MTp , RTf ] R p [MTp , RTf ] f

-1

(30) -1

ing the system step response coefficients, Si, can be obtained as

[

s0 s Si ) 1 ‚‚‚ si-1

0 s0 ‚‚‚ ‚‚‚

‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚

][

g0 0 g 0 ) 1 0 ‚‚‚ gi-1 0

(31)

This projection can be implemented in a numerically robust way with a QR decomposition. The first row of Y ˆ f represents the one-step ahead predictions for the input. Therefore, the white noise disturbance sequence entering the process can be estimated as T

ˆ f(1,:) ef ) [ei, ei+1, ..., ei+j-1] ) Yf(1,:) - Y

(32)

where ‚(1,:) represents a vector containing the elements from the first row and all columns of the matrix. Let us define

ˆ f ) LIeEf Ξ f ) Uf - U

(33)

The block Hankel matrix, Ef, for the noise can be constructed using the estimated noise, ef. Therefore, LIe is estimated as

LIe ) Ξf/Ef ) ΞfE†f

(34)

[ ]

I I ) Hi ‚‚‚ I

0 I ‚‚‚ I

‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚

(AB)

-1

-1

-1

-1

-1 -1

) B A ; A B ) A (B )

-1

) (B A)

-1

) {(HiHci )-1 + I}-1

(35)

[ ] ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚

0 0 0 I

(39)

(40)

Hsi ) -(LIc)-1LIe

(41)

Hsi contains the impulse response coefficients corresponding to the stochastic input.

Hsi

[

I l ) 1 ‚‚‚ li-1

0 I ‚‚‚ ‚‚‚

‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚

0 0 0 0

]

(42)

where lk represents the kth impulse response coefficient of the stochastic input. Thus, the first column of Hsi represents the noise model N(z-1) in IR form.

N(z-1) ) I + l1z-1 + l2z-2 + ... + lkz-k + ... +

LIc ) (I + Hci Hi)-1Hci ) {(Hci )-1(I + Hci Hi)}-1 ) {(Hci )-1 + Hi}-1

]

0 0 0 0 0 I ‚‚‚ I

where sk represents the kth step response coefficient of the deterministic input. 3.1. Closed-Loop Estimation of the Noise Model. The noise model can be estimated from the input data residuals. Using the definitions for LIc and LIe provided in theorem 1, the noise dynamic matrix is obtained as

we can write c -1 c c -1 c -1 LO c ) (I + HiHi ) HiHi ) {(HiHi ) (I + HiHi )}

‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ I I ‚‚‚ I

0 0 0 I

Now, using the identities -1

0 g0 ‚‚‚ ‚‚‚

li-1z-i+1 (43)

(36) 4. The Algorithm

Therefore, I -1 ) {(HiHci )-1 + I}-1{(Hci )-1 + Hi} LO c (Lc)

) {[(Hci )-1 + Hi]-1[(HiHci )-1 + I]}-1 ) {[(Hci )-1 + Hi]-1[(Hci )-1 + Hi]Hi-1}-1 ) (Hi-1)-1 ) Hi

(37)

Hence, Hi, which contains the impulse response coefficients corresponding to the deterministic input, can be identified as

[

g0 g I -1 ) 1 Hi ) LO c (Lc) ‚‚‚ gi-1

0 g0 ‚‚‚ ‚‚‚

‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚

0 0 0 0

]

(38)

where gi represents the ith impulse response coefficient of the deterministic input. The dynamic matrix contain-

The following are the steps in the closed-loop identification: Step I. Construction of the data Hankel matrices {Up, Uf, Yp, Yf, Rp, Rf} using the closed-loop data. By linear regression, the deterministic closed-loop subspace matrices are identified. Remarks: The guidelines presented in section 4.1 can be used in the selection of the number of rows and columns. By adding a persistent exciting signal, which is uncorrelated with the process noise, in the setpoint, we ensure unbiased estimation of the closed-loop subspace matrices. This step is an open-loop identification problem, with the setpoint change being the external inputs and the closed-loop subspace matrices being the models to be identified. Step II. Estimation of the vector of noise data from the “output data Hankel matrix” and “residual data Hankel matrix” corresponding to “input data Hankel matrix”. Estimation of the stochastic closed-loop subspace matrices. Remarks: The first row of the residual ˆ f) represents one-step ahead prediction matrix (Yf - Y errors and an unbiased estimate of the noise entering the process because the feedback does not effect the

848


current noise. The matrix Ξf ) (Uf - U ˆ f) is the residual data Hankel matrix corresponding to the input. The noise data Hankel matrix, Ef, is constructed using the vector of estimated process noise. By linear regression, the stochastic closed-loop subspace matrix, LIe, is estimated. Step III. Retrieving the open-loop deterministic subspace matrix from the closed-loop subspace matrices. Remarks: Closed-loop subspace matrices are just the open-loop subspace matrices weighted by the inverse of the subspace matrices corresponding to the sensitivity function. The analogies between the process/noise transfer functions and the open-loop deterministic/stochastic subspace matrices are obvious. The method presented in this paper is parallel to the approach used in the “joint input-output closed-loop identification method”, which is well-known in the transfer function domain. However, with the “joint input-output closed-loop identification method”, inverting the transfer function (or transfer function matrices for the multivariate case) can give problems such as the resultant transfer function (matrix) maybe being improper or of high order. No such problems are encountered in the subspace matrices based approach proposed in this paper because we are dealing with matrices instead of transfer functions, provided the closed-loop subspace matrices are of full rank. See the guidelines in section 4.1 for keeping the closed-loop subspace matrices from becoming rank deficient. Step IV. Retrieving the open-loop stochastic subspace matrix from the closed-loop subspace matrices. 4.1. Some Guidelines for the Practical Implementation of the Algorithm. Building the data Hankel matrices is the first step in all of the subspace-based identification methods. If we want to identify a statespace model for the system using the subspace identification methods, then the number of rows, i, is chosen to be higher than the order of the state-space model to be identified.29 The number of columns, j, should tend to infinity. Because we have finite data in real situations, we can only choose a finite j, the maximum number of columns that can be constructed with the available data. As far as the closed-loop identification method presented in this paper is concerned, we are identifying only the subspace matrices and not the statespace system matrices. Here are some guidelines that can help in deciding the number of rows and columns of the data Hankel matrices: (a) To obtain the complete process model, the number of rows, i, should be chosen such that the last impulse response coefficient (last element of the first column of Hi) approaches zero. (b) The choice of the number of columns really depends on the excitation signal used for identification. The richer the excitation signal is, the fewer the number of columns required. The number of columns is chosen in such a way that the corresponding impulse response coefficients in the columns of the subspace matrices are very close. The higher the number of columns taken in the data Hankel matrices is, the closer the corresponding coefficients in the columns will be. (c) It may be a good idea to check the rank of the closed-loop subspace matrices before retrieving the open-loop subspace matrices. If the closed-loop subspace matrices are rank deficient, then either decrease the number of rows or increase the number of columns of the data Hankel matrices.

(d) Numerical tools such as QR decomposition can be used to avoid numerical problems associated with the inversion of big matrices, especially in step I. (e) Multivariate systems: Although, in principle, all of the derivations in this paper are applicable to multivariate systems, the numerical problem is a potential concern. As the number of variables increases, the size of the data Hankel matrices can be prohibitively high, especially for systems with a long settling time. MISO identification instead of MIMO identification can reduce the size of the data Hankel matrices. Slower sampling can be used for processes with slow dynamics. Numerical techniques such as QR decomposition, principal component analysis, and CVA can be used to deal with the inversion of large matrices for multivariate systems. (f) Studies on the derivation of statistical properties for subspace-based identification methods are an area of active research and have been considered in refs 17, 19, 29, and 39 and references therein, where it has been shown that under open-loop conditions subspace identification can yield a consistent estimation of the parameters. Because the proposed method in this paper is equivalent to an open-loop subspace identification problem, the same conclusion can be applied. Other statistical properties remain as open problems that need to be further pursued. 5. Closed-Loop Simulations Certain comparative simulations were carried out in MATLAB for two cases, univariate and multivariate systems, between the nonparametric approaches presented in this paper, MOESP40 and CVA.20,37 The purpose of this exercise is to check the validity of the proposed nonparametric approach of closed-loop identification for both univariate and multivariate cases and also to see how it performs compared to the existing subspace identification methods, the MOESP/CVA approaches, which do not require controller knowledge for closed-loop identification. [It has been shown by Van Overschee and De Moor28 that the difference between the three subspace identification algorithms N4SID/ MOESP/CVA is the difference in the way the weighting matrices are used in the subspace identification algorithm. In fact, MATLAB-6 offers a feature called N4weight for the “N4SID” command wherein the user can specify MOESP or CVA and the respective weighting matrices will be used, so that it is equivalent to using MOESP/CVA subspace identification algorithms. Hence, in the simulations presented in this section, the MOESP/CVA weighting matrices are used in the “N4SID” algorithm in MATLAB, instead of writing separate algorithms for MOESP/CVA.] 5.1. Univariate System. Consider the following system:29

[

] [ ]

0.6 0.6 0 1.6161 xk+1 ) -0.6 0.6 0 xk + -0.3481 uk + 0 0 0.7 2.6319 -1.1472 -1.5204 ek -3.1993

[ ]

yk ) [-0.4373, -0.5046, 0.0936]xk + [-0.7759]uk + ek


Figure 4. Univariate system: closed-loop system data.

Figure 6. Univariate system: comparing the true (solid line) IR coefficients with those identified by MOESP/CVA approaches.

same set of closed-loop data. MOESP/CVA identified impulse response coefficients are plotted against the true impulse response coefficients in Figure 6. We can see that both MOESP and CVA identified impulse response coefficients are reasonably well matching with the true impulse response coefficients in this univariate case. Therefore, all three methods yield similar results for this univariate case. 5.2. Multivariate System. Consider the following system taken from the MATLAB/MPC toolbox help manual.

Figure 5. Univariate system: comparing the true (solid line) IR coefficients with those obtained in the subspace matrices (dotted line).

where xk, yk, uk, and ek represent the system state, output, input, and unmeasured random noise, respectively, at time k. A PID controller, 0.1 + 0.08/s + 0.08s, is tuned online for the above system for good setpoint tracking and disturbance rejection performance. We assume that the controller knowledge is unknown for the closed-loop identification. “Closed-loop input-output/ setpoint” data are obtained by exciting the system using a designed “RBS” signal of magnitude 1 for the system output setpoint and random white noise of standard deviation 0.1 in MATLAB-Simulink. The closed-loop data plotted in Figure 4. Using the closed-loop subspace identification method presented in section 3, with rows (i) ) 30 and columns (j) ) 2000 in the data Hankel matrices, the subspace matrices Hi and Hsi are identified. Because of the presence of noise, the upper nondiagonal elements in Hi and Hsi will not be exactly zero but very small numbers (they approach zero as j f ∞). The true impulse response coefficients of the system can be calculated from the state-space system matrices provided above. The identified impulse response coefficients are plotted against the true impulse response coefficients in Figure 5. It is illustrated that the identified impulse response coefficients match very well with the true coefficients. Closed-loop MOESP/CVA is used to identify the deterministic part of the system using the

[

12.8e-s y1(s) +1 ) 16.7s -7s y2(s) 6.6e 10.9s + 1

[ ]

][

-18.9e-3s 21.0s + 1 u1(s) + -19.4e-3s u2(s) 14.4s + 1 3.8e-8s 14.9s + 1 w(s) 4.9e-3s 13.2s + 1

]

[ ]

where {y1(s), y2(s)}, {u1(s), u2(s)}, and w(s) represent the system outputs, inputs, and random noise disturbance, respectively. A state-space-based MPC controller is designed in MATLAB for the system. A sampling period of T ) 2 time units is used in the simulations. Closedloop input-output data are obtained by exciting the system using a designed “RBS” signal of magnitude 1 for the setpoint (rt) and random white noise (wt) of standard deviation 0.1 in MATLAB-Simulink. We assume that the controller and wt are not known for the closed-loop identification. The closed-loop data are plotted in Figure 7. The proposed nonparametric closed-loop identification algorithm is used with rows (i) ) 50 and columns (j) ) 2500 in the data Hankel matrices to identify the subspace matrices Hi and Hsi . The identified impulse response coefficients are plotted against the true impulse response coefficients in Figure 8. It can be seen from the plot that the identified impulse response coefficients match very closely with the true coefficients. Next, the MOESP approach is used to identify the deterministic part of the system using closed-loop data. The impulse response coefficients identified using the

850


more). The order of the end model could probably be reduced using a standard model reduction method but with a compromise in terms of bias error and complexity. The CVA approach20 involves the inversion of a transfer function matrix, which may not always be possible in MATLAB because of the time delays or because of the fact that the resultant matrix could contain improper transfer functions.

6. Identification of the Dynamic Matrix: Practical Application on a Pilot-Scale Plant

Figure 7. Multivariate system: closed-loop system data.

The proposed method for the estimation of the dynamic matrix from closed-loop data is tested on a pilotscale system. The system considered is shown in Figure 1. The input (u) is the inlet water flow rate, and the process variable to be controlled (y) is the level of water in the tank. The tank outlet flow valve is kept at a constant position. The head of the water in the inlet pipe can be considered as (an unmeasured) disturbance. The tank level is controlled by a PID controller, 2.5 + 0.05/s + s. An “RBS” signal of series of setpoint changes to the level is designed in MATLAB. Closed-loop data of the process input, setpoint, and output are collected and plotted in Figure 10. Data Hankel matrices of dimensions rows (i) ) 200 and columns (j) ) 1500 are constructed for the closedloop data, and the subspace matrices Hi and Hsi are identified using the closed-loop identification method presented in the previous sections. The columns of the subspace matrices are plotted in Figure 11. It is illustrated in the figure that the impulse response coefficients in the columns of Hi match each other.

Figure 8. Multivariate system: comparing the true (solid line) IR coefficients with those obtained in the subspace matrices (dotted line).

The accuracy of the impulse response coefficients in the matrix Hi is checked by doing an open-loop identification. Open-loop data are collected by exciting the process with an input “RBS” signal of magnitude 1. The impulse response coefficients identified using the openloop subspace identification method are plotted with the coefficients identified using closed-loop data in Figure 12. We can see that there is some mismatch in the impulse response models in the subspace matrices identified using closed-loop data and those identified using open-loop data. The mismatch may be due to the different operating regions excited between closed-loop and open-loop identifications and the effect of feedback control. The noise model mismatch may be traced to the time-varying nature of the disturbances entering the process.

Figure 9. Multivariate system: comparing the true (solid line) IR coefficients with those identified by MOESP approached.

MOESP approach are compared with the true coefficients in Figure 9, and the result is not quite comparable with the one obtained using the proposed method. It should be noted that a better match between the identified and true coefficients using MOESP could be achieved only with a very high order model (the resultant deterministic model being at least of 20th order or

7. Conclusions This paper provides a subspace identification based method for the identification of process dynamic matrix and the noise model from closed-loop data. The closedloop subspace matrices are first obtained by persistent setpoint excitation of the closed-loop system. The openloop subspace matrices are then retrieved from the closed-loop subspace matrices. The process dynamic matrix is obtained from the deterministic subspace


Appendix: Proof of Theorem 1 Consider eq 26:

[] [ [

][ ]

Yf (I + HiHci )-1Γi (I + HiHci )-1HiΓci Xf + Uf ) -(I + HcH )-1HcΓ (I + HcH )-1Γc Xcf i i i i i i i (I + (I +

HiHci )-1HiHci Hci Hi)-1Hci

] [

(I + HiHci )-1Hsi -(I + Hci Hi)-1Hci

Rf +

Hsi

]

Ef

Using eqs 7 and 25 in the above expression, we can express

[

][ ]

(I + HiHci )-1Γi (I + HiHci )-1HiΓci Xf ) Xcf -(I + Hci Hi)-1Hci Γi (I + Hci Hi)-1Γci L LI

Figure 10. Pilot-scale process: closed-loop system data.

[]

[ ] O

Yp Up Rp Ep

where I O O O I I I I LO ) [LO 1 , L2 , L3 , L4 ]; L ) [L1, L2, L3, L4] c -1 i † LO 1 ) (I + HiHi ) ΓiA Γi -

(I + HiHci )-1HiΓci (∆ci - Aic(Γci )†Hci ) i † c -1 LO 2 ) (I + HiHi ) Γi(∆i - A Γi Hi) +

(I + HiHci )-1HiΓci Aic(Γci )† c -1 c c i c † c LO 3 ) (I + HiHi ) HiΓi (∆i - Ac(Γi ) Hi )

Figure 11. Pilot-scale process: IR coefficients from the columns of the identified subspace matrices.

i † s c -1 s LO 4 ) (I + HiHi ) Γi(∆i - A Γi Hi )

LI1 ) -(I + Hci Hi)-1Hci ΓiAiΓ†i (I + Hci Hi)-1Γci [∆ci - Aic(Γci )†Hci ] LI2 ) -(I + Hci Hi)-1Hci Γi(∆i - AiΓ†i Hi) + (I + Hci Hi)-1Γci Aic(Γci )† LI3 ) (I + Hci Hi)-1Γci [∆ci - Aic(Γci )†Hci ] LI4 ) -(I + Hci Hi)-1Hci Γi(∆si - AiΓ†i Hsi ) Therefore, we can write

[]

[ ] [ ]

Yf LO LO c e ) L D + R + E p p f I Uf Lc LIe f

(44)

Literature Cited Figure 12. Pilot-scale process: comparing the IR coefficients from subspace matrices identified using the open-loop data (dotted line) and the closed-loop data (solid line).

matrix, and the noise model in the impulse response form is obtained from the stochastic subspace matrix. The method can be easily extended to the case of measured disturbances. The results from computer simulations and practical application on a pilot-scale plant are provided to illustrate the proposed closed-loop identification method.

(1) Anderson, B. D. O.; Gevers, M. R. Identifiability of linear stochastic systems operating under linear feedback. Automatica 1982, 18, 195-213. (2) Camacho, C. F.; Bordons, C. Model Predictive Control; Springer-Verlag Ltd.: London, 1999. (3) Cutler, C. R.; Ramaker, B. L. Dynamic matrix controlsa computer control algorithm. AIChE National Meeting, 1979. (4) Cutler, C. R.; Ramaker, B. L. Dynamic matrix controlsa computer control algorithm. Proceedings of the Joint Automatic Control Conference, 1980. (5) Van den Hof, P. M. J.; Schrama, R. J. P. Identification and controlsclosed-loop issues. Automatica 1995, 31 (12), 1751-1770.

852


(6) Van den Hof, P. M. J.; Schrama, R. J. P.; Bosgra, O. H. An indirect method for transfer function estimation from closed-loop data. Proceedings of the 31st Conference on Decision and Control, Tucson, AZ, Dec 1992; pp 1702-1706. (7) Favoreel, W.; De Moor, B. SPC: Subspace predictive control; Technology Report 98-49; Katholieke Universiteit Leuven: Leuven, Belgium, 1998. (8) Favoreel, W.; De Moor, B.; Gevers, M.; van Overschee, P. Closed loop model-free subspace-based LQG-design; Technology Report 98-108; Katholieke Universiteit Leuven: Leuven, Belgium, 1998. (9) Favoreel, W.; De Moor, B.; Gevers, M.; van Overschee, P. Model-free subspace based LQG-design; Technology Report 98-34; Katholieke Universiteit Leuven: Leuven, Belgium, 1998. (10) Favoreel, W.; De Moor, B.; Van Overschee, P.; Gevers, M. Model-free subspace-based LQG-design. Proceedings of the American Control Conference, June 1999; pp 3372-3376. (11) Forssell, U.; Ljung, L. Closed-loop identification revisited. Automatica 1999, 35 (7), 1215-1241. (12) Francis, B. A. A course in H∞ control theory; Lecture Notes in Control and Information Sciences; Springer-Verlag: Berlin, 1987. (13) Gustavsson, I.; Ljung, L.; Soderstrom, T. Identification of processes in closed loopsidentifiability and accuracy aspects. Automatica 1977, 13, 59-75. (14) Hjalmarsson, H.; Gevers, M.; Bruyne, F. For model-based control design, closed-loop identification gives better performance. Automatica 1996, 32 (12), 1659-1673. (15) Van Den Hof, P. M. J.; Schrama, R. J. P. An indirect method for transfer function estimation from closed loop data. Automatica 1993, 29 (6), 1523-1527. (16) Huang, B.; Shah, S. L. Closed-loop identification: a twostep approach. J. Process Control 1997, 17 (6). (17) Jansson, M.; Wahlberg, B. On consistency of subspace methods for system identification. Automatica 1998, 34 (12), 1507-1519. (18) Kadali, R.; Huang, B. A data based subspace approach to predictive controller design; Technical Report; Department of Chemical and Materials Engineering, University of Alberta: Alberta, Canada, 2000; submitted for publication. (19) Knudsen, T. Consistency analysis of subspace identification methods based on linear regression approach. Automatica 2001, 37, 81-89. (20) Lakshminarayanan, S.; Emoto, G.; Ebara, S.; Tomida, K.; Shah, S. L. Closed loop identification and control loop reconfiguration: an industrial case study. J. Process Control 2001, 11, 587599,. (21) Larimore, W. E. Statistical optimality and canonical variate analysis system identification. Signal Process. 1996, 52, 131144. (22) Lee, J. H.; Cooley, B. Recent advances in model predictive control and other related areas. Chemical Process ControlsCPC V; CACHE: 1996; pp 201-216. (23) Ljung, L. System Identification, 2nd ed.; Prentice-Hall: New York, 1987. (24) Ljung, L.; McKelvey, T. Subspace identification from closed loop data. Signal Process. 1996, 52, 209-215.

(25) Morari, M. Chapter model predictive control: Multivariable control technique of choice in the 1990s? Advances in ModelBased Predictive Control; 1994. (26) Ng, T. S.; Goodwin, G. C.; Anderson, B. D. O. Identifiability of linear dynamic system operating in closed-loop. Automatica 1977, 13, 477-485. (27) Van Overschee, P.; De Moor, B. N4SID: Subspace algorithm for the identification of combined deterministic-stochastic systems. Automatica 1994, 30 (1), 75-93. (28) Van Overschee, P.; De Moor, B. A unifying theorem for three subspace system identification algorithms. Automatica 1995, 31 (12), 1877-1883. (29) Van Overschee, P.; De Moor, B. Subspace Identification for Linear Systems: Theory implementation applications; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996. (30) Van Overschee, P.; De Moor, B. Closed loop subspace system identification. Proceedings of the 36th Conference on Decision and Control, 1997; pp 1848-1853. (31) Ruscio, D. D. A method for identification of combined deterministic stochastic systems. In Applications of Computer Aided Time Series Modeling; Aoki, M.; Hevenner, A., Eds.; Springer-Verlag: Berlin, 1997; pp 181-235. (32) Ruscio, D. D. Model based predictive control: An extended state space approach. Proceedings of the 36th Conference on Decision and Control, 1997; pp 3210-3217. (33) Ruscio, D. D. Model predictive control and identification: A linear state space model approach. Proceedings of the 36th Conference on Decision and Control, 1997; pp 3202-3209. (34) Ruscio, D. D.; Fos, B. On state space model based predictive control. IFAC Dynamics and Control of Process Systems; 1998. (35) Schrama, R. J. P. An open-loop solution to the approximate closed-loop approximation problem. Proceedings of IFAC identification and system parameter estimation, Budapest, Hungary, 1991; pp 761-766. (36) Soderstrom, T.; Stoica, P. System Identification; PrenticeHall International: U.K., 1989. (37) Tangirala, A. K.; Lakshminarayanan, S.; Shah, S. L. Closed-loop identification using canonical variate analysis. The 47th CSChE Conference, Edmonton, Canada, Nov 1997. (38) Van den Hof, P. M. J.; Schrama, R. J. P. An indirect method for transfer function estimation from closed loop data. Automatica 1993, 29 (6), 1523-1527. (39) Verhaegen, M. Identification of the deterministic part of mimo state space models given in innovations form from inputoutput data. Automatica 1994, 30 (1), 61-74. (40) Verhaegen, M. Application of a subspace model identification technique to identify LTI systems operating on closed-loop. Automatica 1993, 29 (4), 1027-1040.

Received for review October 23, 2000 Revised manuscript received September 26, 2001 Accepted November 2, 2001 IE000909Q

Estimation of the Dynamic Matrix and Noise Model ... - ACS Publications

Recommend Documents