Identification on Demand Using a Blockwise Recursive Partial Least

540

Ind. Eng. Chem. Res. 2003, 42, 540-554

Identification on Demand Using a Blockwise Recursive Partial Least-Squares Technique† P. Vijaysai,‡ R. D. Gudi,*,‡ and S. Lakshminarayanan§ Department of Chemical Engineering, Indian Institute of Technology, Bombay, Powai, Mumbai 400 076, India, and Department of Chemical and Environmental Engineering, National University of Singapore, Singapore, Singapore

Regression techniques that are used for online estimation and control generally yield poor models if the data is not rich enough. Further, restrictions on the lower limit of the forgetting factor, used to prevent ill conditioning of the covariance matrix, confine applications of such techniques to slowly changing processes only. In this paper a modified blockwise dynamic recursive PLS technique that is based on selection of rich data has been proposed for online adaptation and control. Because of its ability to accommodate a wider range of forgetting factors, the technique is found to track the dynamics of slow as well as fast changes in the processes. The proposed technique has been evaluated for adaptation and control using representative chemical processes taken from the literature. 1. Introduction A majority of the chemical processes are nonlinear and time-varying in nature. Linear models are frequently employed to model and control such processes. As the operating conditions change, it is important that the model and its parameters be updated to realize the benefits of model-based advanced control. Many adaptation schemes have been proposed in the literature, with the most widely used technique being the recursive least-squares (RLS) algorithm. The RLS updates the model parameters with the availability of a new data vector and is computationally less expensive than the batch ordinary least-squares (OLS) technique. The improvement in the computational cost over OLS is achieved by avoiding the direct inversion of the entire covariance matrix of the input block. Any adaptation scheme should preferably be aimed toward updating the model parameters regularly in order to capture the time-varying behavior of the plant. Under the RLS framework, this is achieved through discounting of the past information with a forgetting factor, thus assigning a lower emphasis to the past information versus the new samples collected. There has been a significant modification in the RLS scheme to enhance its adaptive ability. The lack of information in the new data would lead to an exponential growth of the covariance matrix if the past data was continuously discounted. This problem, which is usually referred as “blowup”, has been addressed by introducing a variable forgetting factor.1,2 The improved least-squares technique of Sripada and Fisher3 uses an on-off strategy to avoid parameter drifting and bursting. Despite these modifications, the RLS algorithm encounters a major difficulty in reliable parameter estimation whenever the predictor variables are correlated. This problem is often seen under a closed-loop feedback condition where the inputs and outputs are strongly correlated. † An earlier version of the paper was presented at the 6th IFAC Symposium on Dynamics and Control of Process Systems, Jejudo Island, Korea, 2001. * To whom correspondence should be addressed. Tel.: +91.22. 5767204. Fax: +91.22.5726895. E-mail: [email protected]. ‡ Indian Institute of Technology. § National University of Singapore.

Niu et al.4 have also proposed augmented UD identification (AUDI) for the simultaneous model order and parameter estimation. This technique successfully overcomes some of the numerical drawbacks of RLS and elegantly handles issues such as estimator turnoff, covariance blowup, parameter drift, variance excitation, etc. Unlike the OLS technique and its families, the PLS algorithm has been shown to be an effective tool toward estimating a reliable model when data is correlated.5 Ricker6 investigated the use of PLS and singular value decomposition (SVD) methods to estimate finite impulse response coefficients (IRCs) for dynamic model identification. A dynamic PLS framework that essentially captures the dominant dynamic features of the process under study was proposed by Lakshminarayanan et al.7 The inherent ability of the PLS algorithm to accommodate correlated and noisy data can be used for the recursive parameter estimation under closed-loop operation. Dayal and MacGregor8 proposed an exponentially weighted recursive partial least-squares (RPLS) technique, which uses an improved kernel algorithm for updating the covariance matrix with a relatively lesser computational effort. The RPLS algorithm of Dayal and MacGregor showed a significant improvement in the performance when compared to the RLS. However, their case studies did show a few blowups even with RPLS, and they attributed it to the insufficient process information in the covariance matrix. A blockwise RPLS algorithm for online adaptation was proposed by Qin,9 in which a modified version of the basic RPLS algorithm (proposed initially by Helland et al.10) was used. Thus, blockwise RPLS uses the modified nonlinear iterative partial least-squares (NIPALS) algorithm where the scores rather than the loadings and weighing vectors are normalized. This modification was shown by Qin9 to bring a remarkable improvement in the computational efficiency when used for the adaptation. Most of the adaptation schemes proposed in the literature have been developed for processes that change gradually/slowly with time. However, in recent years, abrupt changes in the plant parameters have also drawn considerable attention in the control community.11,12 Hence, it is desirable that any updation

10.1021/ie020042r CCC: $25.00 © 2003 American Chemical Society Published on Web 01/04/2003

Ind. Eng. Chem. Res., Vol. 42, No. 3, 2003 541

technique be flexible enough to accommodate gradual as well as rapid changes in the process parameters. This paper proposes an extension of the blockwise RPLS algorithm of Qin9 for dynamic model identification under closed-loop conditions. A lack of richness in the data is usually compensated for by the introduction of dither signals. Prior to this, however, it is appropriate to check (i) whether the updation is really necessary in terms of the invalidity of the existing model and (ii) if the proposed updation of the model is possible with the existing data. Toward this end, we propose the use of a metric that is based on the condition number (CN) of an appropriate information matrix, which indicates the adequacy of the data toward model updation. In keeping with the requirement of adaptation for gradual as well as rapid change in the process parameters, a strategy based on the division of the dynamic space, to segregate old and the relatively newer information, is proposed here. It is shown here that this blockwise decomposition enables an automated choice of a wider range of forgetting factors and the subsequent pooling of these blocks greatly facilitates the task of rapid model adaptation. The identification strategy proposed here thus enables rapid model updation only on demand and is in keeping with the Just-in-Time (JIT) modeling paradigm of Cybenko;13 i.e., a model is not reestimated/updated until it is necessary and is feasible. The proposed strategy has been validated by simulations on two representative nonlinear case studies taken from the literature and has been shown to give substantially better identification and control under open-loop and closed-loop conditions. This paper is organized as follows. Section 2 encompasses primary issues related to the possible simplification of the RPLS algorithm for the univariate regression and to the assessment of the information content in the data prior to the model updation. The development of mathematical proofs required for ascertaining the validity of the estimated PLS model, in terms of the CN of the loading matrix, is also covered in the same section. A step-by-step procedure using a modified blockwise adaptation scheme under the RPLS framework is discussed in detail in section 3. Subsequently a new method of estimating the forgetting factor for the blockwise updation is presented in section 4 along with the necessary modifications required. Section 5 dwells upon the some of the practical guidelines and the criteria related to the proposed blockwise RPLSD scheme. Furthermore, the efficacy of the methodology proposed has been illustrated through the representative case studies in section 6 under both the open- and closedloop conditions. 2. RPLS on Demand In this section, we present some of the basic concepts required for the formulation of the blockwise RPLSD algorithm. For simplicity, we focus on the development of the algorithm for single input-single output (SISO) and multiple input-single output (MISO) systems. Let X ∈ Rm×n be the regressor or predictor block with m > n and the corresponding output vector be Y ∈ Rm. The X and Y are related as

assume the model structure to be parametric (ARX), then X will have the lagged values of outputs and inputs against the present output. For simplicity, we assume that the system has no dead time. Thus, the ARX in its general form is represented as

A(z-1) yk ) B(z-1) uk + ek

In eq 2, the polynomials A and B are expressed in terms of backward shift operators (z-1) as A(z-1) ) (1 + R1z-1 + R2z-2 + ... + Rnaz-na) and B(z-1) ) (β1z-1 + β2z-2 + ... + βnbz-nb). Here na and nb denote the number of past outputs and inputs, respectively, that enable the prediction of the current value of the output yk. On expanding eq 2, we get na

nb

Riyk-i + ∑βjuk-j + ek ∑ i)1 j)1

yk ) -

(1)

where E ∈ Rm is an error vector. The above representation can be easily extended for identifying the parameters of linear models. If we

(3)

When eqs 1 and 3 are compared, X, Y, and θ can be expressed as in eq 4. Similarly, for the nonparametric

[

-yk-1 -yk X) l -yk+m-1

uk-1 ... -yk-na uk ... -yk-na+1 l l l ... -yk+m-na uk+m-1

... uk-nb ... uk-nb+1 ··· l ... uk+m-nb

]

Y ) [yk yk+1 ... yk+m]T θ ) [R1 R2 ... Rna β1 β2 ... βnb]T

(4)

case (far-IR or FIR), the X block will have only the past inputs. The general form of the FIR structure is represented in eq 5. The term hi in eq 5 is called the ST

yk )

hiuk-i + ek ∑ i)1

(5)

IRC and the stochastic term (stationary disturbance) is denoted as ek. The number of IRCs (ST) is usually determined as the ratio of the settling time to the sampling time. The above equation can be expressed in a form similar to eq 4, such that the estimated IRC is given by θ ) [h1, h2, ..., hST]T. For identification of processes using either SISO or MISO structures, as was shown by Ricker,6 the calculation of each latent vector using the NIPALS algorithm is a single-step procedure and is free from the problems of convergence. To enhance the computational efficiency of the blockwise updation methodology using the NIPALS algorithm, Qin9 proposed that the scores, rather than the loadings and weighing vectors, be normalized. Let r be the rank of matrix X. The modified NIPALS algorithm for the univariate case can then be written as shown in Table 1. The PLS model between the X and Y blocks can be represented (in compact form) as PLS

Y ) Xθ + E

(2)

(X, Y) 98 {T, P, W, B}

(6)

where T, P, and W are respectively the scores matrix, loadings matrix, and matrix of the weight vectors wi. B ) [b1, b2, ..., bn]T is a n × 1 column vector of inner relation coefficients.

542

Ind. Eng. Chem. Res., Vol. 42, No. 3, 2003

Table 1. Modified NIPALS Algorithm for Univariate PLS 1. initialize E0 ) X and i ) 1 2. calculate wi ) (Ei-1TY)/(|Ei-1Ei-1TY|) 3. calculate ti ) Ei-1wi 4. determine pi ) XTti 5. determine PLS coefficients bi ) YTti 6. calculate Ei ) Ei-1 - tipiT 7. increment i by 1 and continue until the number of latent vectors calculated is equal to r 8. increment i by 1 and continue until the number of latent variables calculated is equal to r

Remark 1. It must be noted that, unlike the traditional PLS decomposition of matrices,5,14 the matrices U, Q, and F are deliberately dropped from the above algorithm because they are not necessary for the univariate case. Thus, the relation between the original and latent variables is given by

Lemma 2. Let X be a full rank matrix. If the number of latent dimensions used is less than the rank r of the matrix X, then the CN of X is always greater than the CN of P. Proof. See Appendix B. Estimating the CN. As mentioned earlier, lack of information in the data may lead to poor parameter estimation. The reliability of the estimation during the online adaptation is measured by the CN of the information matrix. Tridiagonal Property of the Loading Matrix. Various methods have been proposed in the literature for estimating the CN of the matrix X. If the conditions specified in lemma 1 are met, then P is a square matrix. Therefore, EigenDecomposition

(XTX) ≡ (PPT) ≡ (PTP) 98 λhi, 1 < i < n, r (10)

n

X)

tipiT ) TPT ∑ i)1

(7)

If all of the latent variables are not taken into account in reconstructing of the X block, then eq 7 can be expressed as

X ) TkPkT + E; 1 < k < r

(8)

When X is a full rank matrix, r ) min(m,n) ) n. The relationships between other variables, viz., T, P, W, and B can directly be deduced from Table 1. Estimation of the regression coefficient (θ) using PLS is discussed in the sections to follow. When one considers the task of PLS model building using multiple blocks of data, the above algorithm as shown by Qin9 reduces the computational cost significantly. If all of the latent variables equal to r are used in the model building, then performing PLS on any two blocks of input and output data, namely, X1, Y1 and X2, Y2, results in the same regression model as would result from a PLS on P1T, B1 and P2T, B2. Hence

[ ] [

] [

X Y P T B1 P T B1 ≡ PLS 1T PLS X1 Y1 ≡ PLS 1 X2 Y2 P 2 B2 2 2

]

(9)

For estimating good models, rich data are always imperative. Lack of excitation in the input results in a model with poorly estimated model parameters and consequently poor predictive ability. Hence, it is necessary to build or update the model only when the input is sufficiently rich. The information matrix of a rich data set is characterized by a low CN. Hence, the ability of the input data set to result in a good and reliable model is dependent on the CN of the appropriate information matrix (dependent on the choice of model structure such as ARX, FIR, etc.). The CN can, therefore, be used as a key parameter in assessing the data quality for estimation. The following lemmas unveil the relationship between the CN of the input space and the loading matrix, P. Lemma 1. Let X be a full rank matrix. If the number of latent dimensions used in building the PLS model is equal to the rank of X, then the CN of X is the same as the CN of P. Proof. See Appendix A.

which means the eigenvalues, and hence the CN, are the same for the matrices XTX, PPT, and PTP. In this paper, we show that the orthogonality in the vectors of “weight” matrix (W) imparts a special property to the loadings, so that PTP is a tridiagonal matrix (see Appendix C for the derivation). Hence, the estimation of eigenvalues of the tridiagonal matrix (PTP) becomes computationally less expensive15 than the original matrix, i.e., XTX. 3. Online Adaptation Using a Modified Blockwise RPLS Technique The RPLS approach proposed by Qin9 is based on pooling of the data blocks having new information with relatively older information. While such an adaptation mechanism could achieve successful model updating, it is also necessary to check the adequacy of the data block/ data vector in terms of its information content that could be useful for meaningful parameter adaptation. Furthermore, for faster tracking of the set-point or abrupt parameter changes in the plant, it may be necessary to incorporate a wider range of forgetting factors and yet retain the numerical stability of the algorithm. To achieve this, we propose an extension of the blockwise updating mechanism of Qin9 based on the information content in the data matrix and the requirement of model updation. The proposed adaptation strategy is divided into three steps as mentioned below. The key feature of this technique is the division of the overall data block into main and supplementary data blocks. Step 1. Let X and Y be the data sets comprising of the inputs and outputs collected for the initial model estimation. The data are divided into two parts, namely, the main block M and the supplementary data block SM. It is assumed that the model order is known a priori in constructing the input and output blocks. The regression models for the main and supplementary blocks are accordingly denoted as

Y M ) X M BM + E

(11)

YSM ) XSMBSM + E

(12)

It is necessary to ensure that plant dynamics remains invariant while the data are collected to build initial (offline) models.


Using step 1, the PLS models developed for eqs 11 and 12 can be represented as PLS

(XM, YM) 98 {TM, PM, WM, BM}

(13)

PLS

(XSM, YSM) 98 {TSM, PSM, WSM, BSM}

(14)

As an example, eq 15 shows the supplementary data block for estimating the FIR model.

[] [

... um-ST+1 ... um-ST+2 ··· l ... uk-ST

ym+1 um um-1 ym+2 um+1 um , XSM ) YSM ) l l uk-1 uk-2 yk

]

(15)

The main model is estimated using the main data block, which is constructed similarly to the supplementary data block shown in eq 15. The PLS regression on the submodel pooled with the main model yields an aggregate/combined model as shown in eq 16, and the parameters are accordingly denoted by C

[

]

[

]

X Y P T BM f PLS XM YM ≡ PLS M T PSM BSM SM SM {TC, PC, WC, BC} (16) Step 2. The estimated model is deployed online. The new input and output samples collected are augmented to the supplementary data block at the bottom. It is to be noted that the supplementary data block has a window of constant size. Hence, when the new rows of data are augmented to the supplementary block, an equal number of rows representing the older data are deleted from the top of the supplementary data block. The rows eliminated from the supplementary blocks are merged with the main block, and the main model is reestimated. Thus, the second step updates both the main model and the supplementary block. The above procedure is shown in eqs 17 and 18. Let yk+1 and uk be output and input corresponding to the instances k + 1 and k, respectively. The new samples dislodge the first row from the matrices shown in eq 15. The new block is denoted as YSM,N and XSM,N.

[] [

ym+2 um+1 um ym+2 u u , XSM,N ) m+2 m+1 YSM,N ) l l yk+1 uk uk-1

]

... um-ST+2 ... um-ST+3 ··· l ... uk-ST+1 (17)

The main model is updated by augmenting the observations removed from the supplementary data, namely, um and ym+1.

[

]

X Y PLS u M y M f {TM,N, PM,N, WM,N, BM,N} (18) m m+1 Step 3. The updated main model and the submodel estimated are merged to form an aggregate model. The new combined model is shown in eq 19.

PLS

[

]

[

]

P TB PM,NTBM,N ≡ PLS M,N T M,N f XSM,NYSM,N PSM,N BSM,N {TC,N, PC,N, WC,N, BC,N} (19)

Figure 1. Schematic representation of blockwise RPLSD identification.

This aggregate model is used for closed-loop control. Each of the above-mentioned steps are depicted in Figure 1. It can be seen in the above steps that the submodel acts as a forerunner and is the sensitive part of the modeling procedure. Any change in the process parameters is sensed through the submodel, and accordingly the main model can be weighed. The proposed methodology has two important advantages. (i) When the process information is drastically changing, the past information can readily be discounted by using a lower forgetting factor (discussed in the following section). (ii) The division of the data helps to assess the information content in the new data. The lack of excitation in the input signals makes the supplementary data more ill-conditioned. A growth in CN thus gives an indication to stop the updation so that the “blowup” can be avoided. Remark 2. While the proposed methodology presented in this paper takes only the parametric changes into account, detection of changes in the model order is also an important requirement. An overparametrization of the model could take care of the situations where the model order changes with the operating region. However, the overparametrization may require additional care in the choice of the threshold CN. Issues related to identification of the changing model order based on partial correlation analysis are currently under investigation. 4. Forgetting Factor The main feature of the proposed blockwise updation is its ability to track both gradual and sudden changes

544


Table 2. Adaptation Scheme for a Blockwise RPLSD Algorithm 1. divide the data available for the offline analysis into the main and supplementary blocks; if a highly nonlinear process is the object of study, it is advisable to use inverse repeat sequence (IRS) signals to estimate the linear model around the region of operation;21 estimate the main model, submodel, and initial aggregate model 2. put the model adaptation algorithm into online 3. augment the new samples collected to the supplementary data and estimate the submodel 4. estimate the CN of the submodel; if the submodel is well conditioned, update the main model by augmenting the latest information truncated from the supplementary data; the results proved in lemma 1 are used in achieving this 5. estimate the forgetting factor using the updated main model and submodel 6. on the contrary, if the CN is large, stop the updation and repeat step 3

in the parameters. This can be achieved because the algorithm can accommodate the value of the forgetting factor (λ) over the entire range from 0 to 1. Equation 19 with the forgetting factor can be expressed as

PLS

[

]

λPM,NTλ BM,N f {TC,N, PC,N, WC,N, BC,N} (20) XSM,N YSM,N

It should be noted here that when λ is 0, the aggregate model converges to the submodel. On the other hand, if there are no innovations (new information) in the supplementary data, λ could assume a value of 1. The proposed method of estimating λ is based on assessing the change in the information content of the main model and the submodel. In other words, the extent of information to be discounted is estimated as the change in the angle between the first vector of the weight matrices of the main model and the submodel. This is shown in eq 21. As shown in Table 1, the first

λ)

w1,MTw1,SM |w1,M||w1,SM|

(21)

vector of the weight matrix W is a measure of the correlation between the input and output data blocks. Thus, the forgetting factor estimated using eq 21 could be construed as a measure of the change in the correlation between the input and output blocks of the main and supplementary data, respectively. With no change in process parameters, the correlation between the predictor and the response blocks remains fairly constant so that λ is equal to 1. It is important to note that, for the angle criterion to be satisfied, the data blocks must be of equal variance. Otherwise, they have to be scaled suitably (e.g., unit variance). However, because the updation uses only the raw information, the raw vector must be transformed to the scaled form. Let XAS and YAS be the autoscaled blocks of the raw input X and output Y data blocks, respectively. Similarly to eq 6, the PLS on-scaled block is given by PLS

(XAS, YAS) 98 {TAS, PAS, WAS, BAS}

(22)

w1,AS is the first column vector of WAS. Let X h ∈ Rm×n be the matrix, such that all of the elements in each column represent the mean of the corresponding columns of X. Similarly, Y h ∈ Rm represents the mean of the Y block. Let ΣX ) diag(σX1,σX2,...,σXn) and σY represent the standard deviations of X and Y blocks, respectively. Then the relation between w1,AS and w1 is given by

h TY h ]σY-1 w1,AS ) CΣX-1[w1(PTPB)T1n - X

(23)

where 1n is a n × 1 vector with all of its elements equal to 1 and C is a constant. The derivation for eq 23 is given

in Appendix E. The standard deviation of the blocks (ΣX, σY) also needs to be updated with the accumulation of the new data and can be done using the method proposed by Li et al.16 Remark 3. It is to be noted that the difference in the size of the two data blocks may affect the interpretation of the forgetting factor via eq 21. However, the inclusion of the CN test for the supplementary block and provision of an appropriate lower bound on the forgetting factor can circumvent the problem resulting from the unequal sizes of the two data blocks. 5. Data Selection In this section, we discuss some practical guidelines to meet the requirements that lead to good identification under closed-loop conditions. As discussed in earlier sections, the proposed scheme estimates the submodel and updates the aggregate model only if the supplementary block is well conditioned. (It may be recalled that the CN of the supplementary block can be estimated from the loading matrix PSM.) In other words, the aggregate model is estimated only when the supplementary block is well conditioned. The algorithm for updation is summarized in Table 2. 5.1. Size and Threshold CN of the Supplementary Data Block. The supplementary block is the most sensitive part of the algorithm. The speed of adaptation depends on the size of this block; the larger the size of the supplementary block, the slower is the adaptation. Additionally, issues related to the tradeoffs between variance and bias errors, as discussed in various JIT approaches,17 need to be considered when the innovations brought by the new data result in the forgetting factor estimated from eq 21 being very small. In such cases, the model estimated via the supplementary block gets higher weighting after block pooling; therefore, the size of the supplementary block needs to be chosen, keeping the tradeoff between variance and bias errors in mind. For the case studies presented later in this paper, the minimum number of rows in the supplementary block was chosen to be at least 2 times the number of parameters to be estimated. The updation algorithm also requires the specification of a threshold CN. This can be specified reasonably well in terms of the CN of an equivalent supplementary block that results from the usage of a persistently exciting dither signal. Issues related to this specification as well as the analyses in the case when a dither signal is actually applied are currently under investigation. For the case studies considered here, the threshold CN has been fixed based on heuristics. 5.2. Selection of Latent Vectors. As already mentioned, for the recursive updation using PLS, the number of latent vectors must always be equal to the rank r. Lemmas 1 and 2 also show the need for calculation of all of the latent variables in order to


satisfy the CN criteria. On the contrary, using all of the latent variables may sometimes yield poor results6,18 when the model is used for prediction using new data. Under such circumstances, the updation is done using all of the latent variables, but the coefficients used for the predictions are estimated from the first few dimensions, which results in a minimum cross validation error. The selection of the optimal number of latent variables has been discussed in the literature.9,18 If the required number of dimensions (k) is less than the rank r, the updation scheme remains unchanged but the coefficients are estimated for only the selected number of components as follows. T +

H ) (PC,kPC,k ) PC,kBC,k; k e r

(24)

Superscript + signifies the generalized PLS inverse, and the subscript C stands for the combined model. When all of the latent variables are used for the model prediction, it can be proved that the regression coefficient

H ) (PCT)-1BC

(25)

In any case the use of eq 24 is recommended for reasons discussed below. 5.3. Mean Centering and Updation. For the dynamic model identification, the data must be corrected for the mean. Consequently, to capture the time-varying behavior of the plant, the adaptation scheme should include the updation of the mean as well. Many procedures have been proposed in the literature to correct the data for the mean,19 but those most often used are (1) mean centering of the data and (2) adding a column vector to the regressor block with all of its elements equal to 1 (unit vector, 1m ∈ Rm). The proposed method uses the raw data block for the updation. Hence, we propose a new method for correcting the estimated model parameters for the mean. From eqs 1 and 22, the regression coefficient estimated for the raw X and Y data is given by

h ) (PPT)-1PB

(26)

If MX and MY are the means of the X and Y blocks, respectively, the coefficients of the mean-centered data can be estimated as

H ) (PPT - $)-1(PB - ψ)

(28)

The development of eq 27 is shown in Appendix D. Similarly, the equation that follows shows the aggregate PLS model for the two raw input/output blocks, X1, Y1 and X2, Y2.

[ ]

[

HC ) (PCPCT - $C)-1(PCBC - ψC)

(30)

h 1TX h1 + X h 2TX h 2 and ψC ) X h 1TY h1 + X h 2TY h 2. where $C ) X For the updation with the forgetting factor, the mean is also weighed appropriately and is given by

h 1TX h1 + X h 2TX h 2 and ψC ) λX h 1TY h1 + X h 2TY h2 $C ) λX (31) Remark 4. It must be noted that addition of the unit vector to the aggregate block is equivalent to pooling the individual raw data blocks and then mean centering the data. However, the addition of the unit vector to the aggregate block can be shown to give biased parameter estimates, especially when the operating region changes but the dynamics are still invariant. Therefore, we use the above proposed method of correction for the mean, which can be shown to be equivalent to mean centering the individual data blocks and then pooling them. Remark 5. Equation 30 can be used if the regression coefficients are computed using all of the latent variables used for model building and updation. Furthermore, the proposed blockwise RPLSD scheme converges to RLS if r ) n and if all of the n latent variables are included in the model. 6. Case Studies In this section we examine the suitability of the proposed methodology for two different systems under open and closed loops. 6.1. Open-Loop Adaptation. Here, we consider the adaptive ability of the blockwise RPLSD algorithm for the continuous polymerization reactor considered by Doyle et al.20 for the sudden changes in the plant parameters. The reaction under consideration is the free-radical polymerization of methyl methacrylate with azobis(isobutyronitrile) as the initiator and toluene as the solvent. The model for the polymerization is as follows:

x˘ 1 ) 10(6 - x1) + c1x1xx2

(32)

x˘ 2 ) 80u - 10.1022x2

(33)

x˘ 3 ) c2x1xx2 + c3x2 - 10x3

(34)

x˘ 4 ) c4x1xx2 - 10x4

(35)

y ) x4/x3

(36)

(27)

where $ ) X h TX h and ψ ) X h TY h and

h ) Y - 1mMY X h ) X - 1mTMX and Y

are thus determined as

]

X Y P T B1 ) {TC, PC, WC, BC} PLS X1 Y1 ≡ PLS 1T P 2 B2 2 2 (29) The coefficients of the mean-corrected aggregate block

Here x1 is the monomer concentration and x2 is the initiator concentration. The variable y represents the number-average molecular weight and is the ratio of x4 to x3. The volumetric flow rate of initiator u is the manipulated variable. The steady-state conditions are as follows: x1,0 ) 5.506 77, x2,0 ) 0.132 906, x3,0 ) 0.001 975 2, x4,0 ) 49.3818, u0 ) 0.016 783, and y0 ) 25 000.5.

546


Figure 2. Adaptation using a blockwise RPLSD algorithm for sudden changes in plant parameters under open-loop conditions.

Figure 3. Adaptation using blockwise RPLSD (a) before the changes in parameters, (b) during the parametric changes, and (c) after the parameters have changed.

It is assumed that the sudden changes in kinetic parameters change the values of constants c1, c2, c3, and c4 in the above equations significantly at the 1000th instant. The changes are shown in Table 3. The FIR structure was chosen for the adaptation. The plant was perturbed with the IRS signal21 in order to get the best possible linear model. The output was corrupted with Gaussian noise to the extent of 5% of the nominal output values. The number of FIR coefficients required was found to be 21.

Table 3. Assumed Changes in Parameters during Polymerization parameter

before the change

after the change

c1 c2 c3 c4

-2.456 8 0.002 412 1 0.112 191 245.978

-1.8245 0.0018 0.0619 182.66

Figure 2 depicts the adaptive ability of the blockwise RPLSD algorithm. The size of the supplementary data


Figure 4. Adaptation using an exponentially weighted RPLS technique.

Figure 5. Adaptation using an exponentially weighted RPLS (a) before the change in plant parameters, (b) during the change in plant parameters, and (c) after the change in plant parameters.

block was taken as 1/10th of the main block. The changes in the plant parameter caused a shift in the mean of the output. The enlarged views of the updation before, during, and after the parametric changes are shown in

Figure 3. The performance of the blockwise updation is compared with the exponentially weighted RPLS approach shown in Figures 4 and 5. The forgetting factor used for the RPLS was estimated using the equation

548


Figure 6. Blockwise updation under closed-loop conditions using well-conditioned supplementary data blocks.

Figure 7. Closed-loop control without adaptation.

proposed by Fortescue et al.2 and is given as follows.

λ)1-

(

)

xt(XTX)-1xtT σ02N0

t

2

(37)

The forgetting factor is estimated for every new sample of collected (xt). t is the prediction error at any time t, and σ0 is the standard deviation of the expected measurement noise associated with the output. The asymptotic memory length N0 was taken as 50. Figures 3a and 5a depict that the exponentially weighted RPLS algorithm and the proposed blockwise RPLSD perform equally initially. However, as shown in parts b and c of Figure 3, when the parameter changed suddenly at the 1000th instant, the blockwise updation captured the changes much faster. The output estimated using an exponentially weighted RPLS algorithm deviated considerably from the true value, as shown in Figure 5b. The algorithm tends to converge relatively slowly over a long period of time. Table 4 shows the difference between the two techniques in terms of the mean-square prediction error (MSE). 6.2. Closed-Loop Studies. The efficacy of the adaptation scheme was further validated through control of two representative nonlinear systems in a model predic-

Table 4. Performance Evaluation of the Two Techniques Using MSE

blockwise RPLSD exponentially weighted RPLS

MSE before the change

MSE after the change

336.0547 362.2873

341.357 436.417

tive control (MPC) framework. The systems used for the validation study were as follows: (1) a polymerization reactor (discussed above in section 6.1) and (2) a wellmixed continuous fermenter presented by Henson and Seborg.22 The performance of the adaptive algorithm is tested against the locally linearized mechanistic models. 6.2.1. Polymerization Reactor. As in the previous section, the volumetric flow rate u was used to control the number-average molecular weight of the polymer. The FIR structure was used to estimate and update the step response coefficients of the MPC. The output was corrupted with 5% noise of Gaussian nature. The output prediction and the control horizon were selected to be 10 and 3, respectively. The maximum CN for turning the updation off was found to be 200 (refer to section 5.1 for the decision on this choice) and is shown as a dotted line in Figure 6b. The difference in the performance with and without adaptation is shown in Figures 6 and 7, respectively. It


Figure 8. Comparison of the set-point tracking abilities with and without adaptation on change in parameters; an enlarged view.

for the fermenter is as follows:

Table 5. Performance Evaluation of the Blockwise RPLSD in Terms of RMSE RMSE no. av. mol wt (output, y)

RMSE vol. flow rate (∆u)

4310 4989

0.0547 0.0602

with adaptation without adaptation

Table 6. Nominal Model Parameters of the Continuous Fermenter parameter

value

parameter

value

YX/S (g/g) a (g/g) b (h-1)

0.4 2.2 0.2

µm (h-1) Pm (g/L)

0.48 50

parameter km (g/L) ki (g/L)

value 1.2 22

could be seen in Figure 6 that the CN prior to the set-point change at the 50th instant is much higher than the required threshold determined for the updation. After the first set-point change, the CN drops drastically, but it is still not advisable to update the model because the CN is higher than the determined limits and any attempt to update the model at this instant would lead to “bursting” of the adaptive algorithm. The process parameters were changed at the 150th instant. This sudden change resulted in a loss of correlation between the inputs and outputs and led to the drastic drop in CN. The submodel thus estimated for the supplementary data block with a lower CN was augmented to the earlier model with an appropriate forgetting factor. The updated model tracked the changes much faster than the original model. Figure 8 depicts an enlarged view of the differences in the tracking abilities of the updated and original models. Similar changes in CN and differences in the performance can seen at the 250th instant when the set point was changed. The overall differences in the performance are shown in Table 5 in terms of the root-mean-square error (RMSE). 6.2.2. Continuous Fermenter. The continuous fermenter analyzed by Henson and Seborg22 is considered next in order to validate the proposed strategy for identification and control. The reactor uses the productfree and sterile feed. The inputs to the reactor are the dilution rate (u1) and feed substrate concentration (u2). The effluent biomass or cell concentration (x1), substrate concentration (x2), and the product concentration (x3) are the process state variables. The model equation set

dx1 ) -u1x1 + µx1 dt

(38)

µx1 dx2 ) u1(u2 - x2) dt YX/S

(39)

dx3 ) -u1x3 + (Rµ + β)x1 dt

(40)

y1 ) x1

(41)

y2 ) x3

(42)

and

where µ is the specific growth rate, Rµ + β is the specific product formation rate, and YX/S is the yield of the cell mass with respect to the limiting substrate. Using substrate and product inhibition models, the specific growth rate µ is given by

[

µm 1 µ)

]

x3 x Pm 2

x22 km + x2 + k1

(43)

The parameters of the model are listed in Table 6. In this work we restrict the analysis to the SISO case. The objective is to obtain a dynamic model that describes the effect of u2 on y1 while the dilution rate is maintained constant at 0.1636. An ARX structure of order 6 was used for capturing the local dynamic behavior of the fermenter. Figure 9a shows the closed-loop control of the continuous fermenter without adaptation. The output is corrupted with identical and independently distributed noise to the extent of 5%. The control and prediction horizon was chosen to be 3 and 7, respectively. In this problem, the difference in the performance with and without adaptation for the three different operating regimes is studied. The initial set point is 5.88, which corresponds to the first 50 samples shown in Figure 9a. The local dynamic model estimated around this operating region ensures satisfactory closed-loop performance with a relatively smaller controller effort. This could be seen in Figure 9b, which depicts the absolute value (not expressed in

550


Figure 9. (a) Closed-loop control of the continuous fermenter without adaptation. (b) Profile of the manipulated variable without adaptation.

Figure 10. (a) Closed-loop control of the continuous fermenter with adaptation. (b) Profile of the manipulated variable with adaptation. (c) Change in CN with set-point change.

the deviation form) of the manipulated variable around each set point. The change in the set point from 5.88 to 4.5 at the 50th instant and later to 7.2 at the 200th instant is seen to result in good control performance at the set-point value of 4.5, thus indicating that, even though adaptation is not active, the model is reasonably adequate. However, the set-point change to 7.2 is shown to result in poor tracking and substantial variation

around the set point, indicating the incentive for model updation. Furthermore, the application of an impulselike disturbance at the 275th instant is seen to result in runaway conditions in the reactor. It could seen in Figure 9a that the control scheme breaks down at this point as the controlled variable deviates from the target value and tends to move toward zero. At the same time there is a ramp growth in the manipulated variable; the


cause for this behavior when investigated was found to be the change in the biomass concentration gain directionality. Figure 10 shows the closed-loop performance with adaptation for the same set-point and disturbance changes as those for parts a and b of Figure 9. The adaptation is seen to yield faster set-point tracking and better disturbance rejection. The perturbations introduced due to the change in the operating region lead to the fall in CN. The maximum threshold CN for the reliable model estimation was found to be 560. Unlike the previous case, Figure 10b shows that the disturbance introduced at the 275th instant was rejected with a much smaller controller effort. Thus, the proposed adaptation and control methodology performed quite well under the closed-loop conditions, as shown by the above case studies. In a summary of the evaluation of the proposed strategy, it is important to note the following two points. First, closed-loop identification can greatly be enhanced if a dither signal is used assuming the resulting output variability is acceptable. However, the strategy proposed here presents a metric in terms of CN of the information matrix, which is helpful to estimate the minimum variance of the dither signal necessary for smoother adaptation. Because the set-point changes by themselves bring in innovations and enhance the CN of the information matrix, only a minimal level of excitation as determined by CN is necessary via the dither signal route. Second, in the method proposed here, CN can change even if there is a sustained disturbance. Obviously, under such circumstances, the occurrence of a disturbance needs to be detected and the adaptation should be turned off. Some recently proposed methods11,23 of distinguishing the occurrence of the disturbance from parametric changes could be used to turn off adaptation in the presence of a sustained disturbance.

If all of the latent vectors are estimated, then k ) r and E ) 0.

X ) TPT

Using the SVD, the X block can be decomposed as

X ) UΣVT

XTX ) PPT

X ) TkPkT + E; 1 e k < r where r ) min(m,n) ) n (44)

(47)

From eq 46, it follows that

XTX ) (UΣVT)T(UΣVT)

(48)

However, UTU ) I ∈ Rn×n. Hence, XTX ) V(ΣTΣ)VT. Let ΣTΣ ) Λ where Λ ) diag(λh1, λh2, ..., λhn) such that λh1 > λh2 > ... > λhn. It could be noted that λhi’s are the eigenvalues of XTX. Hence,

XTX ) VΛVT

(49)

The CN of X is given by

κ(X) )

σ1 ) σn

x

λh1 λhr

(50)

However, from eqs 47 and 49, it follows that

PPT ) VΛVT

(51)

where eq 51 is the spectral decomposition of the matrix PPT such that the matrices V and Λ are the eigenvectors and the eigenvalues of PPT, respectively. Hence, from eqs 50 and 51

κ(X) ) κ(PT)

(52)

Because r ) n, κ(P) ) κ(PT). Therefore,

κ(X) ) κ(P)

(53)

Appendix B Proof. If the number of latent vectors selected is equal to k (k < n), then eq 7 can be expanded as follows:

X ) t1p1T + t2p2T + ... + tkpkT + E; 1 e k < r (54) where the residual matrix E can be further decomposed as

Appendix A Proof. From eq 8, the full rank matrix X can be decomposed as

(46)

where U ∈ Rm×m and V ∈ Rr×r are orthonormal matrices. Though the actual dimensions of the left singular matrix, U, are (m, m), only the first r columns of the matrix are chosen. Hence, the dimension of the singular value matrix (Σ) are taken as (r, r) instead of (m, n) such that Σ ) diag(σ1, σ2, ..., σr) ∈ Rr×x. Using eq 45, XTX ) (TPT)T(TPT). However, because T T T ) I ∈ Rr×r,

7. Conclusion In this paper a new modified blockwise RPLSD technique, based on segregation of relatively old and new data, for online estimation and adaptation was presented. A metric, based on the CN of the loadings matrix, to assess the suitability of the existing data toward model updation was proposed and evaluated. This metric could also be used to determine the minimal level of dither signals necessary in the absence of sufficient information during transitions and is currently being investigated. Considering the hitherto heuristic nature of the exponential forgetting strategy, the method proposed in this paper provides a more quantitative basis for the choice of the forgetting factor. The strategy proposed was found to be efficient for both the online and offline adaptations. In particular, the proposed methodology was shown to be in line with the JIT modeling philosophy. Results obtained from simulations involving representative chemical process systems demonstrated the efficacy of the proposed methodology.

(45)

E ) tk+1pk+1T + tk+2pk+2T + ... + tnpnT

(55)

XTX ) PkPkT + 2PTTE + ETE

(56)

Hence,

Because the scores Tk are orthogonal to E

552


XTX ) PkPkT + ETE

EiTY ) (Ei-1 - tipiT)TY ) Ei-1TY - pitiTY (68)

(57)

but tiTY ) bi. Hence,

Upon substitution of eq 55 into eq 57,

XTX ) PkPkT + pk+1pk+1T + pk+2pk+2T + ... + pnpnT (58)

EiTY ) Ei-1TY - bipi

When eq 69 is rearranged, pi can be expressed as

The pipiT matrices are each of rank 1, and PkPkT is the matrix of rank k. Therefore, Pk+1Pk+1T given by eq 59 is a matrix of rank k + 1.

Pk+1Pk+1T ) PkPkT + pk+1pk+1T

(59)

eigen decomposition

pi )

wi )

Because λh1 > λh2 > ... > λhk+1, the CN of matrix Pk+1 is

(61)

Similarly, let λhi and δi be the eigenvalues of PkPkT and pk+1pk+1T such that λhk+1 ) λhk+2 ) ... ) λhr ) 0 and δ2 ) δ3 ) ... ) δr ) 0. λhi and λhi can be related as15

λhi - λhi ) fiδ1

such that (a) 0 e f e 1 and ) 1 and (b) fk+1δ1 < λhk. This leads to an important result as follows:

λh1 > λh1 and λhk > λhk+1

|Ei-1Ei-1 Y|

pi )

(

ciwi - ci-1wi+1 bici-1ci

)(

piTpj )

(64)

If i ) j, piTpj )

(66)

(71)

1 |EiEiTY|

ciwi - ci-1wi+1 bici-1ci

(63)

When the above analogy is extended to any number of latent vectors less than the rank r of the matrix X, it can be shown that

) ci-1Ei-1TY

(72)

Therefore, T

)

cjwj - cj-1wj+1 bjcj-1cj

ci2wiTwi + ci-12wi+1Twi+1 bi2ci-12ci2

If j ) i + 1, piTpj )

(65)

κ(X) > κ(Pk)

T

Hence,

Comparing eqs 61 and 64 and using the condition shown in eq 63, it can be seen that the CN of Pk+1 is always greater than the CN of Pk.

κ(Pk+1) > κ(Pk)

Ei-1TY

ci )

The CN of Pk is given as

κ(Pk) ) xλh1/λhk

(70)

where

(62)

k+1 ∑i)1 fi

Ei-1TY - EiTY bi

From Table 1, any ith weight vector can be estimated as

Pk+1Pk+1T 98 λhi, i ) 1, 2, ..., r; λhk+2 ) λhk+3 ) ... ) λhr ) 0 (60)

κ(Pk+1) ) xλh1/λhk+1

(69)

-wjTwj bibjci2

(73)

(74)

(75)

If j > i + 1, piTpj ) 0 (because the columns of the weighing matrix are always orthogonal, wiTwj ) 0 if i * j) (76) Equations 74-76 completely satisfy conditions a and b. Hence, PTP is a tridiagonal matrix.

Appendix C

[

Given P ) [p1, p2, ..., pn], then

p1Tp1 p1Tp2 ... p1Tpn p Tp p Tp ... p2Tpn PTP ) 2 1 2 2 · ·· l l l T T pn p1 pn p2 ... pnTpn

]

Appendix D

(67)

The method to estimate the raw regression coefficients for the mean-corrected data is shown here. Let the mean of the matrix X be denoted by the row vector Mx. The mean matrix X h is thus estimated as shown

1m1mT X h ) 1mMx ) T X 1m 1m

If PTP has to be tridiagonal, then (a) piTpj ) cij, ∀ j e i + 1 and 1 e i < n and (b) piTpj ) 0, ∀ i + 1 < j e n and 1 ei < n. As discussed in Table 1, the PLS residuals of X block are given by

where 1m ) [1, 1, ..., 1]T ∈ Rm. The mean-centered matrix Xm ) X - X h . Therefore,

Ei ) Ei-1 - tipiT

XmTXm ) XTX - 2XTX h +X h TX h

Hence

but

(77)

(78)


X h )

[1m1mT]X 1mT1m

)

[1m1mT]T[1m1mT] X [1mT1m]2

(79)

dard deviation of the columns of the X block. Similarly, σ is the standard deviation of the Y block. From Table 1, the first vector of the weighing matrix W is given by

Upon substitution of the second term into the third, T T

X h )

[1m1m ] X h [1mT1m]

w1 )

T

)

[1m1m ]X h [1mT1m]

(80)

X X h )

)

1m1mT XT T X h 1m 1m

w1,AS )

[ ]

XASTYAS T

|XASXAS YAS|

) CXASTYAS

(91)

where

1m1mT T X X h 1mT1m

h )X h TX

(90)

Hence, w1 corresponding to the autoscaled data blocks X and Y is given by

Therefore, T

XTY |XXTY|

C) (81)

1 |XASXASTYAS|

w1,AS can further be expanded as

Therefore,

w1,AS ) CΣ-1XmTYmσ-1

XmTXm ) XTX - 2X h TX h +X h TX h h TX h ) XTX - X

(82)

However,

XmTYm ) XTY - X h TY h (from eq 83)

Similarly, it can be proved that

XmTYm ) XTY - X h TY h

T

X X ) PP

T

and X Y ) PB

(84)

Substituting eq 84 into eqs 82 and 83 yields

(85)

h TY h XmTYm ) PmBm ) PB - X

(86)

where Pm and Bm are the matrices corresponding to the mean-centered data blocks. From eq 45, the regression coefficient for the raw data is given by

h ) (PP ) PB

where

|XXTY| ) (PTPB)T1n (see Appendix F for the derivation)

w1,AS ) CΣX-1[w1(PTPB)T1n - X h TY h ]σY-1

T -1

H ) (PmPm ) PmBm

Proof that |XXTY| ) [PTPB]1n. Let η ) XXTY, and η is a vector of length m × 1. The norm of η is given by

|η| ) xηTη ) x(XXTY)T(XXTY)

(95)

Upon expansion and rearrangement of the terms in the above expression

|η| ) xYTX(XTX)XTY

Substituting the above results, we obtain

(88)

(94)

Appendix F

(87)

Hence, for the mean-centered data blocks, it follows that

H ) (PPT - X h TX h )(PB - X h TY h)

(93)

Substituting the above results for w1,AS, we get

h TX h XmTXm ) PmPmT ) PPT - X

T -1

h TY h ) w1|XXTY| - X

(83)

The normalization procedure deployed in the RPLS scheme yields T

(92)

(96)

However, XTX ) PPT and XTY ) PTTY ) PB. Therefore

|η| ) xBTPTPPTPB

Appendix E

h )Σ-1 XAS ) (X - X

or

|η| ) x(PTPB)T(PTPB)

and

YAS ) (Y - Y h )σ-1

(89)

where Σ is the diagonal matrix representing the stan-

(97)

which means the norm of XTXY is equal to the norm of PTPB. However, as shown in Appendix C, PTP is atridiagonal matrix, such that

554

[


PTP )

c12w1Tw1 + c02w2Tw2 b12c02c12 w2Tw2 b1b2c12 0

-

w2Tw2

b1b2c12 c22w2Tw2 + c12w3Tw3 b22c12c22 0

l

l

0

0

and B ) [b1, b2, ..., bn]T. Thus, PTPB, when expanded and simplified, yields the vector of length n × 1 as follows:

T

P PB )

[

w1Tw1 b1c02

0 0 ... 0

]

T

(99)

Therefore,

η ) |PTPB| ) [PTPB]1n )

w1Tw1 b1c02

(100)

Hence,

|XXTY| ) [PTPB]1n Literature Cited (1) Ljung, L. System Identification: Theory for the User; Prentice Hall: New York, 1987. (2) Fortescue, T. R.; Kershenbaum, L. S.; Ydstie, B. E. Implementation of Self-Tuning Regulators with Variable Forgetting Factors. Automatica 1981, 17, 831-835. (3) Sripada, N. R.; Fisher, D. G. Improved Least Squares Identification. Int. J. Control 1987, 46 (6), 1889-1913. (4) Niu, S.; Fisher, D. G.; Xiao, D. An Augmented UD Identification Algorithm. Int. J. Control 1992, 56 (1), 193-211. (5) Geladi, P.; Kowalski, B. R. Partial Least-Squares Regression: A Tutorial. Anal. Chim. Acta 1986, 185, 1-17. (6) Ricker, N. L. The Use of Biased Least-Squares Estimators for Parameters in Discrete-Time Pulse-Response Models. Ind. Eng. Chem. Res. 1988, 27, 343-350. (7) Lakshminarayanan, S.; Sirish, S. L.; Nandakumar, K. Modeling and Control of Multivariable Processes: The Dynamic Projections to Latent Structures Approach. AIChE J. 1997, 43, 2307-2322.

...

0

...

0

...

0

··· -

wn-1Twn-1 bn-2bn-1cn-12

T

-

wn-1 wn-1 bn-2bn-1cn-12 wnTwn bn2cn-12

]

(98)

(8) Dayal, B. S.; MacGregor, J. F. Recursive Exponentially Weighted PLS and its Applications to Adaptive Control and Prediction. J. Process Control 1997, 7, 169-179. (9) Qin, S. J. A Recursive PLS algorithm for Adaptive Data Modeling. Comput. Chem. Eng. 1998, 48, 503-514. (10) Helland, K.; Berntsen, H. E.; Borgen, O. S.; Martens, H. Recursive Algorithm for Partial Least Squares Regression. Chemom. Intell. Lab. Syst. 1992, 14, 129-137. (11) Huang, B. On-line Closed-loop Model Validation and Detection of Abrupt Parameter Changes. J. Process Control 2001, 11, 699-715. (12) Basseville, M.; Nikiforov, I. V. Detection of Abrupt Changes; Prentice Hall: Englewood Cliffs, NJ, 1993. (13) Cybenko, G. Just-in-Time Learning and Estimation. In Identification, Adaptation, and Learning; Bittani, S., Picci, G., Eds.; NATO ASI Series; Springer: New York, 1996; pp 423-434. (14) Hoskuldsson, A. PLS Regression Methods. J. Chemom. 1988, 2, 211-228. (15) Wilkinson, J. H. The Algebraic Eigenvalue Problem; Oxford University Press: New York, 1992. (16) Li, W.; Yue, H. H.; Valle-Cervantes, S.; Qin, S. J. Recursive PCA for Adaptive Process Monitoring. J. Process Control 2000, 10, 471-486. (17) Stenman, A.; Gustafsson, F.; Ljung, L. Just in Time Models for Dynamical Systems. Proceedings of the 35th IEEE Conference on Decision and Control; IEEE Press: Piscataway, NJ, 1996; pp 1115-1120. (18) Dayal, B. S.; MacGregor, J. F. Identification of Finite Impulse Response Models: Methods and Robustness Issues. Ind. Eng. Chem. Res. 1996, 35, 4078-4090. (19) Box, G. E. P.; Jenkins, G. W. Time Series Analysis: Forecasting and Control; Holden Day: San Francisco, 1976. (20) Doyle, F. J.; Ogunnaike, B. A.; Pearson, R. K. Nonlinear Model-based Control Using Second-order Volterra Models. Automatica 1995, 31 (5), 697-714. (21) Ranganathan, S.; Raghunathan, R. Use of Inverse Repeat Sequence (IRS) for Identification in Chemical Process Systems. Ind. Eng. Chem. Res. 1999, 38, 3420-3432. (22) Henson, M. A.; Seborg, D. E. Nonlinear Control Strategies for Continuous Fermenters. Chem. Eng. Sci. 1992, 47, 821-835. (23) Huang, B. Multivariable Model Validation in the Presence of Time Varying Disturbance Dynamics. Chem. Eng. Sci. 2000, 55, 4583-4595.

Received for review January 16, 2002 Revised manuscript received June 19, 2002 Accepted October 31, 2002 IE020042R

Identification on Demand Using a Blockwise Recursive Partial Least

Recommend Documents