A Modified Extended Recursive Least-Squares Method for Closed

Ind. Eng. Chem. Res. 2009, 48, 6327–6338

6327

A Modified Extended Recursive Least-Squares Method for Closed-Loop Identification of FIR Models Christopher L. Betts, Srinivas Karra, M. Nazmul Karim, and James B. Riggs* Department of Chemical Engineering; Texas Tech UniVersity, Lubbock, Texas 79409

Performing a plant test under closed-loop conditions is desirable for model identification because production loss and safety problems may result when control loops are opened during plant testing. However, identification of models from closed-loop data is more difficult compared to identifying models from open-loop data because of the correlation between the colored noise and the process inputs created by the feedback. A novel method for identifying models using closed-loop data is proposed, which employs a time-varying bias term with a moving average dynamic component as the model structure. Then identification for this process model is performed using a modified extended recursive least-squares algorithm to eliminate the bias from the process parameter estimates. Evaluation of the proposed algorithm is performed using simulation case studies involving multivariable processes controlled by either diagonal PI controllers or a model predictive controller (DMCPlus). The simulation results showed that the proposed method is robust with respect to measurement noise and is able to identify high quality models from closed-loop data. 1. Introduction Model predictive control (MPC), which is a multivariable controller based on empirical models, has been widely accepted by the chemical processing industries (CPI) because of its ability to control near the economic optimum operating conditions.1 MPC has a number of advantages over conventional PID-type controllers because it combines optimization and control using empirical process models.2 MPC is effective for most processes including processes with numerous constraints, time-delay processes, inverse response processes, and highly coupled processes.2 Dynamic matrix control (DMC), which was proposed and developed by Charles R. Cutler in the 1970s, is the first known and most popular form of MPC in the CPI. DMC uses step response models (SRM) of a process to calculate control action. SRMs are linear models of a process that describe the relationship between the input and output data. Usually, SRMs are derived from other types of models (e.g., finite impulse response (FIR) models, output error (OE) models, autoregressive with exogenous inputs (ARX) models, and autoregressive moving average with exogenous inputs (ARMAX) models). FIR models are used in this research due to their simplicity, flexibility, and widespread industrial use. Typically, identification of FIR models used for control is performed by applying a plant test to the process under openloop conditions to prevent correlation in the training data. Conversely, if the plant test data is obtained under closed-loop operation, it exhibits correlation between the colored noise on the process measurements and process inputs because of the feedback from the controller, and therefore, the identified model parameters can be biased. However, identification of process models from closed-loop data is advantageous for the following reasons:3,4 (1) better models for control are possible, (2) closedloop identification can allow for closer operation to the economic optimum operating conditions during the plant testing, (3) testing can be conducted in less time, and (4) some plants may be unstable in open-loop operation. The fundamental problem with closed-loop identification is that the process noise is correlated with the process inputs because of the feedback of the closed* To whom correspondence should be addressed. Tel.: (806)7421765. E-mail: [email protected].

loop system.5 The correlated noise is usually neglected with the assumption that the closed-loop signals are noise free, or the input spectrum dominates the noise spectrum.3 If the correlation is not properly handled, it introduces bias in the model estimate.4 Therefore, research in the field of closed-loop identification is largely focused on solving this correlation issue. Over the years, three main approaches to closed-loop identification have been introduced, the direct, indirect, and joint input-output approaches. In the direct approach,6 the process output and the process input are used the same way as in the open-loop identification approach, ignoring any feedback. The control structure does not need to be known a priori. In the indirect approach,3 the closed-loop system is identified between a reference signal to the regulator and the process output. Then, the open-loop model of the plant is retrieved making use of the known control structure. In the joint input-output approach,5 the process output and the process input are considered outputs of a system driven by the reference input and noise. Knowledge of the controller and the plant are recovered. Forssell and Ljung5 revisit closed-loop identification by discussing common approaches that are classified within the direct, indirect, and joint input-output methods used in this field of system identification in the prediction error framework. Prediction error methods are aimed at improving parameter estimation with respect to consistency, minimum variance in the estimate, and estimating a disturbance model.4 Forssell and Ljung state that the prediction error methods can always be used in direct identification with any system, for example, stable or unstable systems, as long as the predictor is stable. For indirect methods, the control structure must be known and the set point must be measurable. For the joint input-output methods, the basic assumption is that the system input is generated using a regulator of some form. Therefore, using a closed-loop expression for the system output and input, estimates for the openloop system are computed. Most researchers in this field assume the controller is linear, unconstrained, and the process is single input single output (SISO). However, industrial applications are not linear because the plant inputs and outputs are constrained6 and the process is nonlinear and is not SISO due to the MIMO

10.1021/ie800922w CCC: $40.75  2009 American Chemical Society Published on Web 06/05/2009

6328

Ind. Eng. Chem. Res., Vol. 48, No. 13, 2009

nature of the system. According to Forssell and Ljung, the indirect and joint input-output approaches are usually applied when the feedback is linear. As a result, the direct closed-loop identification approach is the obvious choice when MPC controllers are used.6 Therefore, only direct methods are considered in this work. Karimi and Landau compare several common closed-loop identification methods in terms of the bias distribution.3 They conclude that the classical direct methods for closed-loop identification do not approach the desired white noise distribution in the presence of colored process noise by determining the bias distribution. However, the recently developed family of output-error algorithms for plant model identification in closed-loop has a bias distribution which is not influenced by noise. To handle the bias in the estimated parameters, some researchers have suggested using a prefilter to identify the process model by filtering the process input and output data before identifying the process model. MacGregor and Fogal discuss the role of the noise model and prefilters for closedloop identification.7 Prefiltering the process data will provide an unbiased estimate of the process parameters due to the reduction of the effect of colored noise on the closed-loop system. In practice, however, the noise structure may not be known a priori. In this case, Gopaluni suggests the identification method can be performed in two steps to determine the proper structure of the noise model.8 A bias-correction method for indirect identification of closedloop systems uses a bias-correction procedure to apply a leastsquares type method that is capable of handling the bias problem in indirect model identification.9 The proposed method requires fewer computations compared to other methods and, according to simulation case studies, is attractive with regard to the accuracy of the estimate. A similar approach can be applied to cases where the linear process systems are corrupted with noisy input and output observations.10 In this approach, the variance of the noise is estimated directly instead of the covariance vector. Once the noise variances are estimated, the parameter bias can be estimated and removed from the problem. The previous proposed method for indirect identification of closed-loop systems assumes that the order of the controller was not less than the order of the open-loop plant. However, most controllers are low-order. Therefore, Zheng11 modified his original indirect method to remove this assumption. Simulations conducted for testing the new modified BELS algorithm have yielded favorable results. A new version12 of the BELS algorithm is introduced which eliminates the design of the prefilter. This new method aims at providing an identification algorithm that is more numerically efficient and practical without using a prefilter, which uses dynamic error-in-variable (DEV) models. Two similar indirect identification methods presented in the literature are the two-stage13 and projection identification methods.14 The basis for these methods is to obtain process input data that is uncorrelated with the process noise indirectly using an identified sensitivity function for the regulator. Both methods yield favorable results. However, the projection method yields slightly more accurate results compared to the two-stage method according to Forssell.14 Other methods presented in the literature include a family of recursive algorithms for model identification in closed-loop.15 A recursive algorithm is usually used in adaptive control where the model is updated online.4 However, a recursive identification algorithm can be applied offline to consider time-varying parameters.16 Owing to this feature, many researchers are turning toward recursive parameter estimation for closed-loop identi-

fication. Landau and Karimi present a unified approach and evaluation using recursive algorithms for identification in closedloop.15 Recursive output error (OE) models are studied in the context of direct closed-loop identification. Simulation case studies were conducted to compare the different forms of the OE model. The results for all model forms were favorable. A new recursive algorithm for indirect ARMAX model identification in closed-loop is also proposed by Landau and Karimi.17 The purpose of developing the recursive algorithm for closedloop identification for ARMAX models is to ensure that global asymptotic stability is achieved, an asymptotic optimal predictor of the true closed-loop system is obtained, and an asymptotically unbiased estimate of the plant model parameters are obtained. Thus far, the developed methodologies for closed-loop identification found in the literature are devoted to solving the problem of the correlation between the noise and process inputs. However, most of these methodologies involving indirect or joint input/output approaches are not practical in the sense that industry standards use model predictive controllers (MPC) with FIR models. These approaches require a structured controller that is described by a known transfer function a priori. Unfortunately, MPC is not a structured controller. Therefore, it is well-known that only the direct or two-step closed-loop identification methods which do not require any information on the controller transfer function are useful for these cases. It was also reported in the literature that both the direct and twostep methods can be used to provide asymptotically unbiased results with parsimonious or nonparsimonious model structures. The underlying assumption is that the model structures contain the true process model and either the true disturbance model (direct method), or the true sensitivity function model (twostep method).29 However, with finite data sets and for disturbances that approach nonstationarity, both the latter techniques can produce biased results. In practical situations the noises (or disturbances) are nonstationary in nature with properties that are time-varying. In chemical process industries, external disturbances in the form of limit cycles caused by aggressive tuning of control loops, valve nonlinearities, interacting processes and/or hydrodynamic instabilities are prevalent, resulting in properties that change with time if there is a change in operating conditions of the process unit from which they are originating.30 Other sources of nonstationary disturbances are changes in ambient weather conditions, varying feed stocks or feed flow rates, changes in operating conditions or product grades, process variations such as heat exchanger fouling and catalyst deactivation in upstream process units, transients created by the events such as bypassing flows, switching pumps, changes in cooling or heating loads, sudden plugging and unplugging of pipelines, etc. Closed-loop systems employing linear/nonlinear control strategies with constraints may also generate nonstationary noise in input signal even if the measurement noise in the controlled variable is stationary in nature. These facts motivated us to devise a new system identification technique. The proposed methodology falls under the category of direct closed-loop identification techniques for identifying FIR models. The difficulty associated with identifying and checking for adequate parsimonious model structures was avoided by using nonparsimonious FIR models to represent process dynamics. To account for the nonstationary disturbances a nonparsimonious disturbance model is augmented with FIR model structure and combined parametrization is performed to obtain accurate process model. Therefore, the intended contribution of this work is to provide a means of identifying MIMO FIR models from


6329

closed-loop plant test data in the presence of colored noise on the process measurements. 2. The MX-FIR Method This section provides a derivation of the proposed modified extended finite impulse response (MX-FIR) algorithm for closed-loop identification as well as a theoretical analysis of it. Assume the true linear process is represented as an FIR-type model and is given by y(t) ) Gu(t) + e(t) ) φ(t)θ + e(t)

(1)

where G is the process model, u(t) is the process input, θ ) [g1 g2 · · · gn]T is the vector of process model parameters, φ(t) ) [u(t - 1) u(t - 2) · · · u(t - n)] is the regressor vector of process inputs, and e(t) is zero mean white noise. The objective is to estimate the parameters of the process model using closed-loop data. The least-squares solution to eq 1 in vector notation is θˆ ) [ΦTΦ]-1ΦTY where

[ ] [ ]

y(n + 1) y(n + 2) Y) l y(t) φ(n + 1) φ(n) Φ) l φ(t - 1) u(n) u(n + 1) ) l u(t - 1)

[

(2)

Figure 1. Process under closed loop. (a) Case 1 with a stationary zero mean white noise entering the system, es; (b) case 2 with a nonstationary external disturbance entering the system, ens.

Letting the bias term µ(t) ) ∆Gu(t) + β(t), eq 5 can be rewritten as, y(t) ) Gu(t) + µ(t) + e(t)

and

(6)

Remark 1. The linear predictor for the process shown in eq 6 is given by

u(n - 1) u(n) l u(t - 2)

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

u(1) u(2) l u(t - n)

]

ˆ u(t) + µˆ (t) yˆ(t) ) G

(3)

This solution is a good approximation for FIR models in openloop identification if the following assumptions are true. (A) The process noise and unmeasured disturbances are independent of the input data. (B) The process is open-loop stable. (C) Adequate signal-to-noise ratios are obtained during the plant test. (D) The prediction error, ε(t) ) y(t) - yˆ(t), is zero mean white noise. In other words, E{ε(t)} ) 0, where E is the expectation operator. (E) The magnitudes of unmeasured disturbances are low. However, under closed-loop feedback control, ε(t), which is the open-loop residual or prediction error, is not zero mean white noise, E{ε(t)} * 0. Even though the measurement noise assumes a white noise structure, because of (linear/nonlinear with or without constraints) feedback the noise properties that cycle through the control loop will be different than the properties specified for the measurement noise. If the model identification is based on the model form of eq 1, biased model parameters are obtained in prediction error methods. To remove this bias, an additional time varying term, β(t), is added to the original model form as shown in eq 4. β(t) models the nonstationary component of the noise in the system caused by feedback. This bias term shifts the closed-loop prediction error, εcl(t), to a zero mean. y(t) ) Gu(t) + β(t) + e(t)

(4)

(7)

which yields consistent predictions of the process in closed loop, ˆ is the process model where yˆ(t) is the process output estimate, G estimate, u(t) is the process input, and µˆ (t) is the bias estimate. The following development is provided in support of this statement. Consider two cases shown in Figure 1a,b. Case one (Figure 1a) with the process under closed loop with a zero mean white noise and case two (Figure 1b) with the process under closed loop with a zero mean white noise and nonstationary disturbance entering the system; u is the process input, y˜ is the true process output, es is the stationary zero mean white noise, ens is the nonstationary effect of the external disturbance entering the system on y, ysp is the controller set point, and ε is the error between the set point and the measured process output y. First, consider case one with zero mean white noise entering the system according to Figure 1a. The process output, y, at time t is represented as y(t) ) y˜(t) + es(t)

(8)

The actual process output is represented as y˜(t) ) φ(t)θ

(9)

assuming that the actual process dynamics are truly linear. θ is the true process parameters and φ(t) is the regressor vector of inputs as described before. The measured output y at time t is expressed as y(t) ) φ(t)θˆ + ε(t)

(10)

In the case of slight additive nonlinearity (represented by ∆G) present in the process, the true process model is shown in eq 5.

where θˆ is the parameter estimate, ε(t) is the prediction error. Transforming eq 10 into vector notation becomes

y(t) ) Gu(t) + ∆Gu(t) + β(t) + e(t)

Y ) Φθ + E

(5)

(11)

6330


The measurement vector Y, observation matrix Φ, are defined in eq 3. The prediction error vector is given by E ) [ε(n + 1) ε(n + 2) · · · ε(t)]T. The least-squares solution to eq 11 is given by θˆ ) [ΦTΦ]-1ΦTY

(12)

If the model is consistent with the process, the least-squares solution yields a zero mean white prediction error sequence. Substituting eq 11 into the above equation yields θˆ ) [ΦTΦ]-1ΦT(Y˜ + Es)

) [ΦTΦ]-1ΦT(Φθ + Es)

(13)

) [ΦTΦ]-1ΦTΦθ + [ΦTΦ]-1ΦTEs

Since the noise vector is zero mean white noise, Es is uncorrelated with the process inputs. Hence, [ΦTΦ]-1ΦTEs is a null vector and [ΦTΦ]-1ΦTΦ is an identity matrix which yields θˆ ) θ. Consider case two with a nonstationary disturbance entering the system according to Figure 1b. Nonstationary disturbance is defined as additive and time-varying component of the disturbance. Assuming the same process described by eq 9, the measured process output y is defined for this case as Y ) Φθ + Ens + Es

(14)

where Ens is the nonstationary external disturbance entering the process and all other variables are defined as before. The leastsquares estimator for linear models assumes the prediction error is zero mean white noise which is a good assumption for openloop disturbance-free training data. However, if any colored or nonstationary component is present in the measured process output, and model structure is not inclusive, the process parameters are modeled as if the colored or nonstationary component reflects the actual process dynamics and this introduces a bias into the process parameter estimates.18 Therefore, the measured process output is defined as Y ) Yˆ + E ) Φθˆ + E

(15)

where θˆ is now the biased estimate of process parameters in this case. Using eqs 14 and 15, the following expression is derived: Φθ + Ens + Es ) Φθˆ + E

(16)

Simplifying,

n

ens(t) )

∑ce

i ns(t

- i) + µ(t)

Thus, the new model form is defined as follows, y(t) ) φ(t)θ + C(q-1) e(t) + µ(t)

(17)

The true process output is rewritten as Y˜ ) ) )

Φ(θˆ - [ΦTΦ]-1ΦTEns) Φθˆ - Φ[ΦTΦ]-1ΦTEns Yˆ - Ens

θˆ (0) ) 0,

P(0) ) FI

φ(t) ) [u(t - 1) u(t - 2) · · · u(t - n) e(t - 1) e(t - 2) · · · e(t - n) ] (23)

ii. Estimation of the Gain Vector. The gain vector, K(t), is estimated according to the following equation. This vector can be viewed as a weighting vector for parameter estimation. K(t) ) P(t - 1) φ(t)T[1 + φ(t) P(t - 1) φ(t)T]-1

y˜(t) ) φ(t)θˆ + ens(t)

(19)

Here, the nonstationary component of the above equation is modeled as a moving average process with mean µ(t) as follows,

(24)

iii. Estimation of the Residual Vector. The residual vector, ε(t), is then estimated by the following equation: (25)

iv. Estimation of the Process Parameters. The process parameters, θˆ (t), are updated by use of the following equation: (26)

v. Estimation of the Covariance Matrix. Lastly, the covariance matrix, P(t), is estimated using the following equation. P(t) ) P(t - 1) - P(t - 1) φ(t)T × φ(t) P(t - 1)[1 + φ(t) P(t - 1) φ(t)T]-1

Finally, the true process output is written as

(22)

where F is a large number, P is the covariance matrix, and I is the identity matrix. F is chosen as a large number so that the algorithm does not trust the initial parameter estimates. The regressor vector, φ(t), is defined as

θˆ (t) ) θˆ (t - 1) + K(t) ε(t)

(18)

(21)

where φ(t) is the regressor vector of inputs, y(t) is the measured process output, θ is the vector of process parameters, C(q-1) is the transfer operator defined as C(q-1) ) 1 + ∑kn ) 1ckq-k, and µ(t) is defined as the bias parameter as a function of time as before. Remark 2. The above model form is nonparsimonious in terms of modeling process dynamics and disturbance dynamics. The φ(t)θ component is equivalent to FIR dynamics, C(q-1) e(t) captures the zero-mean quasi-stationary disturbance, and µ(t) captures the disturbance that causes a slow drift in output as bias. 2.1. Identification Algorithm. The solution of eq 21 for θˆ cannot be solved by the normal least-squares method. Therefore, a recursive algorithm18 is used to search for the solution of this equation. Recursive algorithms are used mainly for identifying models online or, in other words, adaptive control. However, recursive algorithms are also used for identification offline if certain model parameters vary with time.18 The version of the pseudolinear recursive least-squares algorithm used to solve eq 21 in this work is shown below. i. Initiation. The initial values, at time, t, equal to zero for the recursive algorithm are taken as

ε(t) ) y(t) - φ(t) θˆ (t - 1) θ ) θˆ - [ΦTΦ]-1ΦTEns

(20)

i)1

(27)

Steps ii to v are repeated until t is equal to N, the size of the data set. Once closed-loop data is generated, the above algorithm is used to identify the process model. Persistent excitation in the


6331

external signal is needed for good convergence in the parameter estimates. The parameter estimation can become unstable because the term [1 + φ(t) P(t - 1) φ(t)T] may approach zero if there is no significant movement in the process input. If the above term does go to zero, the gain vector, K(t), will go to infinity which will cause unstable parameter estimation. 2.2. Consistency of the Predictor. Assume the following process

disturbance, C describes the transfer function for the noise model, and e is zero mean white noise. First, an expression for the input, u, to the process must be derived in terms of all excitation variables. The excitation variables for this system are the reference signal or the set point, r, and the noise, e. The following expression was derived using block algebra.

y(t) ) Gu(t) + Ce(t) + µ

where S is a sensitivity function and Tr is a transfer function of the form

(28)

where y ∈ Rm and u ∈ Rn. And also it is assumed that the system has persistent excitation. G, C, and µ are the true process model (m × n transfer function matrix), stationary disturbance model (m × 1 transfer function vector) and nonstationary additive disturbances (m × 1 vector), respectively. In SISO case, m ) n ) 1. The one step ahead predictor, yˆ(t/θ), ˆ u(t) + µˆ ] + [1 - Cˆ-1]y(t) yˆ(t/θ) ) Cˆ-1[G

(30)

The loss function to be minimized is N

∑ ε(t/θ)

2

t)1

)

1 N

N

∑ (Cˆ

-1

ˆ u(t) - µˆ ])2 [y(t) - G

t)1

(31) Substituting y(t) into ε(t/θ) yields ˆ -1

(32) This shows that, as N f ∞, ε(t/θ) ) e(t) + (a term independent of e(t)). Then, ˆ u(t) - µˆ ])2 VN(θ) f Eε(t/θ)2 ) E(Cˆ-1[y(t) - G V∞(θ) ) Eε(t/θ) g Ee(t) ) R 2

2

The one-step ahead predictor for the given model form is yˆ(t/θ) ) C-1[Gu(t) + µ] + [1 - C-1]y(t)

ˆ - G) u(t) + (µˆ - µ)] + yˆ(t/θ) ) C-1[(G (Gu(t) + Ce(t) + µ) - e(t)

(39)

(40)

The expression for the prediction error is as follows, ˆ - G) u(t) + (µˆ - µ)] + εCL(t/θ) ) y(t) - yˆ(t/θ) ) C-1[(G e(t) (41)

ˆ - G)(Tr r(t) + Sµu µ(t) + Seu Ce(t)) + εCL(t/θ) ) C-1[(G (µˆ - µ)] + e(t) (42) Using Parseval’s relation, the bias distribution is expressed as follows.

(33) θ* ) (34)

where R is the variance of the error, e(t).

||

ˆ | 2Eu(t)2 + |Cˆ-1 | 2 |µ - µˆ | 2 + C 2R Eε(t/θ)2 ) |Cˆ-1 | 2 |G - G Cˆ (35) For an infinite number of samples the parameter estimates are obtained by the minimization of the expectation of the squared prediction error given by eq 35, θ ) arg min Eε(t/θ)2 θ

ˆ | ) 0, In such case, the global minimum to be obtained, |G - G ˆ ) G, µˆ ) µ, |µ - µˆ | ) 0, and |C/Cˆ| ) 1. This implies that G and Cˆ ) C. Therefore, the model is consistent. 2.3. Bias Distribution. The bias distribution is important in explaining how the bias is distributed over the frequency domain for a given model by showing how all of the excitation signals affect the estimation of the model parameters. Consider the following process under closed loop, y(t) ) GPu(t) + Ce(t) + µ(t)

(38)

Substituting the expression for u(t) into the above equation results in

ˆ ]u(t) + C [µ - µˆ ] + C e(t) ε(t/θ) ) C [G - G Cˆ ˆ -1

GC GC and Tr ) 1 + GCGP 1 + GCGP

(37)

Substituting y(t) into the one-step ahead predictor yields the following result.

ˆ u(t) - µˆ ] ε(t/θ) ) Cˆ-1[y(t) - G

1 N

Sµu ) Seu )

(29)

and the prediction error, ε(t/θ) ) y(t) - yˆ(t/θ), will be

VN(θ) )

u(t) ) Trr(t) + Sµuµ(t) + SeuCe(t)

(36)

Let, GC and GP be the transfer functions for the controller and process, respectively, µ describes the constant or drifting

∫

π

-π

ˆ (eiω)| 2Tr2 φr(ω) + |C(eiω)-1 | 2[|G(eiω) - G

ˆ (eiω))Sµu + 1]2µ(eiω)2 + µˆ (eiω)2 + [(G(eiω) - G S2euC(eiω)2 φe(ω)] dω (43) where Γ(eiω) is the frequency domain representation of discrete ˆ , µ, or µˆ . According to the transfer function Γ. Γ can be C, G, G above bias distribution, the bias of the process model, G(eiω) ˆ (eiω), defined in eq 21, is uncorrelated with the spectrum of G the noise, φe(ω), because both terms are separated within the bias distribution which is desirable for closed-loop identification. If the excitations in the reference signals, φr(ω), are designed in a way that is uncorrelated with the excitation due to noise, Seu2 C(eiω)2 φe(ω), the bias of the process model will be uncorrelated with the noise spectrum. The terms [(G(eiω) ˆ (eiω))Sµu + 1]2 µ(eiω)2 and µˆ (eiω)2 are also uncorrelated with G the noise spectrum, the former term defined as the bias caused by any unmodeled nonlinearities and the latter term defined as ˆ (eiω) ) 0, the bias estimate would the bias estimate. If G(eiω) - G equal the actual bias, µ(eiω)2. Furthermore, from this analysis, when performing a closed-loop plant test, the spectrum of the reference signal must dominate the spectrum of the noise. In other words, the signal-to-noise ratio level must be large enough so that the spectrum of the noise does not affect the bias distribution (e.g., a 5 to 1 signal-to-noise ratio3). 2.4. Stability Analysis. A stability analysis of the recursive algorithm is performed to prove the stability of the algorithm

6332


in estimating the model parameters offline. This proof was adopted from the approach used by Landau.17 Theorem. Provide the following assumptions: a. The open-loop system is in the model set and the closedloop system is stable. b. The controller is known and asymptotically stable. c. The external excitation, r(t), is bounded and |e(t + 1) ≡ 0| is defined. The recursive algorithm assures for all initial conditions lim ε(t + 1) ) 0 tf∞

φ(t) < c,

| |

0 < c < ∞ for all t

Proof. Rewrite the prediction error given in eq 41, as E{e(t)} ) 0, ˆ - G) u(t) + (µˆ - µ)] εCL(t/θ) ) C-1[(G

(44)

In vector notation εCL(t + 1) ) [θ - θˆ (t + 1)]T φ(t)

(45)

To prove this theorem, one has to show that φ(t) is bounded. The components of φ(t) are u(t - i), ε(t - i), and 1. Taking account of eq 37 and assumptions B and C, u(t) is bounded. By assumption B, ε(t) is bounded.; 1 is also bounded because it is a constant value. Therefore, φ(t) is bounded which results in a stable recursive algorithm. 3. Illustrative Example 1: 2 by 2 Transfer-Function Model In this section, the Wood-Berry column19 is used to evaluate the MX-FIR method using closed-loop data generated by a centralized PI controller. 3.1. Process Description. The Wood-Berry distillation column separates methanol from water, where methanol is the main product of the system. This column is a two product column with an (L, V) control configuration,1,19 that is, reflux and reboiler duty as the MVs and the overhead and bottoms compositions as the CVs. A methanol/water stream enters the column and a methanol rich stream leaves the column as distillate. The water-rich waste stream leaves the column as a bottoms product. The process model for the Wood-Berry column is as follows:19 The compositions of the top and bottom products (xD and xB) expressed in weight % of methanol are the controlled variables. The reflux and the reboiler steam flow rates (R and S) are the manipulated inputs expressed in lb/min.

[

12.8e-s xD(s) 16.7s + 1 ) xB(s) 6.6e-7s 10.9s + 1

[ ]

-18.9e-3s 21.0s + 1 -19.4e-3s 14.4s + 1

]

[ ] R(s) S(s)

Table 1. The Calculated Residual Statistics for the Wood-Berry Column with White Noise and with Colored Noise white noise case

colored noise case

transfer function

SSE

maximum percent error

SSE


xD(t)/R(t) xD(t)/S(t) xB(t)/R(t) xB(t)/S(t)

17.85 35.5 7.6 13.3

9.0% 8.3% 9.1% 3.7%

25.5 48.9 9.86 17.7

10% 8.9% 9.5% 4.0%

used to excite the process and generate the closed-loop data. Unless otherwise specified, the type of test signal used for the plant test in all the case studies presented in this paper, is a step-change test signal where the intervals for the step changes are 1/4Tss, 1/2Tss, 3/4Tss, Tss, 5/4Tss, where Tss is the time to steady state.1 In this case, the amplitudes of the steps in xD and xB set points are (0.7 and (0.5, respectively. The control configuration consists of the reflux and steam flow rates as the MVs and the methanol concentration in the distillate and bottoms product as the CVs. The identified SRMs obtained agreed well with the process model defined above. Both the process gains and time constants were calculated for this case from the SRM and compared to the actual process parameters by using heuristics.1 The average and maximum percent errors for the identified process gain, Kp, compared to the actual Kp are 0.30% and 0.44%, respectively. The average and maximum percent errors for the identified process time constant, τp, compared to the actual τp are 0.20% and 0.37%, respectively. Table 1 contains the calculated residual statistics, for example, sum of the square errors (SSE) and the maximum percent errors, for this case. Next, a simulation case study with the Wood-Berry column is performed with colored noise added to the simulator as filtered white noise. The colored noise model used in this case study is shown below, where V(s) is the colored noise and e(s) is zero mean white noise. From observing the autocorrelation of this noise model, V(s) is highly correlated. 1 V(s) ) 2 e(s) 2s + 10s + 1

(47)

The level of noise applied to the CVs is (0.10% (σe ) 0.01) on CV1 and (0.30% (σe ) 0.03) on CV2. Figures 2 and 3 show the identified SRM and the calculated residual, respectively, for this case. The residual is defined as the value of the difference between the identified model predictions of the process output and the true process values for a unit step change in input. The average and maximum

(46)

Distillation columns experience significant coupling and can have a wide range of dead-time to time constant ratios. For this case, the dead-time to time constant ratio ranges between 0.06 and 0.64. The above process is controlled by a diagonal PI controller; xD is controlled by manipulating R (Kc ) 0.5, ti ) 0.01), xB is controlled by manipulating S (Kc ) -0.03, ti ) -0.005). 3.2. Results and Discussion. First, a simulation case study with the Wood-Berry column is performed with no noise added to the simulator. The set points to diagonal PI controllers were

Figure 2. Identified SRM for the Wood-Berry column case with colored noise added: (s) identified SRM, ( · · · ) true SRM.


6333

the main product for this system. The column has an overhead product, a side stream product, and a bottoms product. The overhead product has a high concentration of ethanol and the bottoms product has a low concentration of ethanol. The side stream product is a product with an intermediate methanol concentration. Close control of the overhead product is important because it is the feed to other high purity columns. The bottoms product is also important because it is the feed material to a final ethanol clean up column. The temperature of tray number 19 is used as an inferential measurement of the concentration of ethanol in the bottoms product because the concentration of ethanol is strongly correlated with this tray temperature. The process model for the ethanol/water column is as follows:20 Figure 3. The calculated residual for the Wood-Berry column case with colored noise added.

percent errors for the identified process gain, Kp, compared to the actual Kp are 0.90% and 1.2%, respectively. The average and maximum percent errors for the identified process time constant, τp, compared to the actual τp are 0.56% and 0.78%, respectively. According to Figure 3, the calculated residual is noisy for this case compared to the previous case because of the addition of colored noise. The residual in the previous case study is relatively smooth and damps out over some time. However, from Table 1, the SSE and the maximum percent error of the calculated residuals are relatively low, indicating a relatively accurate model. The closed-loop test that was conducted for this case used a signal-to-noise ratio of 7 for CV1 and 5 for CV2. To evaluate the magnitude of the variation in the colored noise, the level of noise is increased to (0.20% (σe ) 0.02) on CV1 and (0.60% (σe ) 0.06) on CV2. The average and maximum percent errors for the process gain, Kp, are 1.1% and 1.6%, respectively. The average and maximum percent errors for the identified process time constant, τp, compared to the actual τp are 0.61% and 0.84%, respectively. Finally, the noise is increased again to (0.40% (σe ) 0.04) and (1.2% (σe ) 0.12) on CV1 and CV2, respectively. The average and maximum percent errors for the process gain, Kp, are 1.3% and 1.9%, respectively. The average and maximum percent errors for the identified process time constant, τp, compared to the actual τp are 0.66% and 0.91%, respectively. For each of these cases, the signal-to-noise ratio was constant and equal to that used for the original case (i.e., Figures 2 and 3). Therefore, from these results, the magnitude of the colored noise does not significantly affect the accuracy of the models even though the magnitude of the errors does increase with an increase in the magnitude of the colored noise. 4. Illustrative Example 2: 3 by 3 Transfer-Function Model In this section, an empirical model of a pilot scale distillation column, which is used for the separation of ethanol and water,20 is used to evaluate the performance of the MX-FIR method for identifying FIR-type models from closed-loop data. A simulation was built using SimuLink to generate closed-loop data that is used in the model identification algorithm. The goal of this case study is to evaluate the MX-FIR method by using a higher dimensioned system. The ethanol/water distillation column used in this case study is a 3 by 3 system (i.e., 3 MVs and 3CVs). 4.1. Process Description. The ethanol/water distillation column has 19 trays, numbered from the top to the bottom of the column, where it separates ethanol from water. Ethanol is

[ ]

xD(s) xS(s) ) T19(s) -0.61e-3.5s -0.0049e-s 0.66e-2.6s 6.7s + 1 8.64s + 1 9.06s + 1 -2.36e-3s -0.012e-1.2s 1.11e-6.5s × 3.25s + 1 5s + 1 7.09s + 1 0.87(11.61s + 1)e-s -33.68e-9.2s 46.2e-9.4s 8.15s + 1 10.9s + 1 (3.89s + 1)(18.8s + 1) FR(s) FS(s) (48) S(s)

[

][ ]

xD is the overhead ethanol concentration (mole %), xS is the ethanol concentration in the side stream (mole %), T19 is the temperature on tray no. 19 (°F), FR is the overhead reflux flow rate (lb/min), FS is the side stream draw-off rate (lb/min) and S is the reboiler steam pressure (psia). The above process is controlled by a diagonal PI controller; xD is controlled by manipulating FR (Kc ) 1.2, ti ) 0.19), xS is controlled by manipulating FS (Kc ) -0.001, ti ) -0.01), and T19 is controlled by manipulating S (Kc ) 0.005, ti ) 0.05). 4.2. Results and Discussions. First, a simulation case study with the ethanol/water column is performed with no colored noise added to the simulator. The set points to diagonal PI controllers were used to excite the process and generate the closed-loop data. In this case, the amplitudes of the steps in xD, xS, and T19 set points are (0.7, (0.5 and (20, respectively. The control configuration consists of the reflux flow rate (FR), side streamflow rate (FS), and reboiler steam pressure (S) as the MVs and the overhead concentration (xD), side stream concentration (xS), and tray number 19 temperature (T19) as the CVs. The identified SRM agreed well with the process model defined above. The process gains and time constants were calculated and compared to the actual process parameters from the identified SRM. The average and maximum percent errors for the identified process gain, Kp, compared to the actual Kp are 0.26% and 0.31%, respectively. The average and maximum percent errors for the identified process time constant, τp, compared to the actual τp are 0.18% and 0.28%, respectively. The SSE and maximum percent error are also low according to Table 2. The maximum percent error corresponds to the calculated maximum percent of the maximum residual obtained in this case study. Next, a simulation with the ethanol/water column is performed with colored noise added to the simulator. The level of noise applied to the CVs is (0.3% (σe ) 0.003) on xD, (0.5% (σe ) 0.005) on xS, and (0.5 OF (σe ) 0.05) on T19. The colored noise model used in case 1 is also used in this case. Figures 4 and 5

6334


Table 2. The Calculated Residual Statistics for the Ethanol/Water Column with No Noise case with white noise transfer function xD(s)/FR(s) xD(s)/FS(s) xD(s)/S(s) xS(s)/FR(s) xS(s)/FS(s) xS(s)/S(s) T19(s)/FR(s) T19(s)/FS(s) T19(s)/S(s)

SSE -4

7.6 × 10 2.5 × 10-3 7.7 × 10-7 4.6 × 10-3 0.15 2.56 × 10-6 2.16 39.2 0.15

case with colored noise


SSE


2.1% 0.5% 1.2% 2.5% 2.2% 0.8% 3.3% 5.6% 1.9%

7.9 × 10-4 5.2 × 10-3 8.2 × 10-7 1.3 × 10-2 0.31 4.9 × 10-6 4.1 52.6 0.32

2.8% 1.1% 1.6% 2.9% 2.6% 2.2% 2.8% 10.2% 3.2%

show the identified SRM and calculated residual, respectively, with colored noise added to the simulator. The process gains and time constants were calculated from the identified SRM and compared to the actual process parameters. The average and maximum percent errors for the identified Kp compared to the actual Kp are 0.88% and 1.7%, respectively. The average and maximum percent errors for the identified τp compared to the actual τp are 1.0% and 1.9%, respectively. The calculated residual shown in Figure 5 is noisy compared to the simulation with no colored noise, which is relatively smooth and damps out over some time. However, the SSE and maximum percent error (Table 2) are also relatively low. The closed-loop test that was conducted for this case used a signal-to-noise ratio of 6 for xD, 4 for xS, and 8 for T19. To show the effect of the magnitude of the variation of the colored noise, the noise is increased to (1.2% to xD, (2.0% to xS, and (2.0 °F to T19. The average and maximum percent errors for the process gain are 1.3% and 2.2%, respectively. The average and maximum percent errors for the identified process time constant compared to the actual process time constant are 1.4% and 2.0%, respectively. For each of these cases, the signalto-noise ratio was constant and equal to that used for the original case (i.e., Figures 4 and 5). Once again the effect of a relatively large increase in the colored noise resulted in a relatively small increase in the errors for the calculated process gain and time constant. 5. Illustrative Example 3: Nonlinear FCC Simulator A detailed simulator of a fluidized catalytic cracking (FCC) process is used as a case study to further evaluate the proposed method. 5.1. FCC Process Description. A typical side-by-side FCC unit is considered in this work as shown in Figure 6. First, hydrocarbon feedstock (e.g., gas oil) is heated in the feed preheater (G) prior to entering the riser reactor (E) where it is mixed and vaporized with the hot catalyst from the regenerator (RG). In the reactor riser, endothermic catalytic reactions occur as the hydrocarbon vapor rises along the length of the reactor. The products from the reactor include lighter hydrocarbons and coke which deposits onto the catalyst surface which causes catalyst deactivation. The main products, which are light gases, gasoline, and heavier cuts, are then sent to the main-fractionator (MF). The MF separates the feed into a variety of product streams ranging from light gases to residual oil. The deactivated catalyst is separated by the cyclone (F) and falls into the stripping section (D) where steam strips any residual products off the deactivated catalyst. Then the spent catalyst is transported to the regenerator (RG) where the coke is burnt off of the catalyst. The reactivated catalyst is sent back to the reactor riser

to vaporize and react with the feedstock through a transport line. Thus, the regenerated catalyst continuously enters the riser and the spent catalyst is sent to the regenerator yielding a continuous cyclic reaction process. In industry, auxiliary units are included along with the reactor, regenerator, and main fractionators, such as the following. Wet gas compressor (A): The wet gas from the mainfractionator is compressed by the wet gas compressor prior to being sent to the gas unit for separation. Feed preheater (G): Prior to the feedstock entering the reactor, it must be raised to a desired temperature, therefore, it is heated by a furnace. The waste light gases from the main-fractionator are used as the fuel for the furnace. Catalyst cooler (H): A catalyst cooler is used to remove any excess heat energy from the catalyst using a heat exchanger. More than one cooler can be used. Consequently, steam is produced by cooling the catalyst. Stack gas expander (I): High-pressure stack gas is employed to generate electricity or to drive the air blower using the expander. Air blower (J): The air blower delivers a large amount of air to the regenerator in order to feed the combustion reactions occurring in the regenerator. The regenerator not only removes coke from the surface of the catalyst but it also generates enough heat to drive the endothermic reactions in the riser reactor. CO boiler (L): Because of excess amounts of carbon monoxide in the stack gas, a CO boiler is used to reduce the amount of this pollutant by converting it into carbon dioxide. The process model of a side-by-side FCC process was adapted from Han.21 In this research, this detailed FCC model was used to test the proposed closed-loop identification methodology. The FCC simulator was constructed using a modular structure which was developed by Han and Chung22 implemented as a Fortran code. The FCC model has the following features: (1) a 10-lump model22 is used to model the conversion in the riser reactor; (2) the simulator employs momentum balance equations to predict molar expansion and catalyst slip in the riser; (3) the simulator uses the two-region, two-phase theory for the regenerator model;2 (4) the models for most auxiliary units (the feed preheat system, stack gas expander, catalyst cooler, blower, wet gas compressor, and the CO boiler) were added to capture their dynamics; and (5) the model employs the momentum balance to predict the molar expansion and catalyst slip in the riser. The entire dynamic model in the simulator comprises a mixed system of differential and algebraic equations.22 The HYBRD1 routine of MINPACK23 was used in the solution of the nonlinear algebraic equations. It uses the Powell’s hybrid algorithm, which is based on a Newton’s method.22 The purpose of the HYBRD1 routine is to find the zeros of n nonlinear functions with n unknown variables. The method of lines was used to solve the partial differential equations at each time step. In addition, the LSODE24 routine, which is based on the Gears algorithm, is used to solve the ordinary differential equations resulting from the application of the method of lines as well as other model equations. These numerical methods are used in the FCC simulator to calculate the state variables at each time step. The simulator communicates with an IP-21 database, which communicates with the DMCPlus control software.25,26 For each control cycle, the output of the DMC controller is written back to the simulator input file and the outputs are then calculated by the simulator.


6335

Figure 4. Identified SRM for the ethanol/water column case with colored noise added: (-) identified SRM, ( · · · ) true SRM.

Figure 5. The calculated residual for the ethanol/water column case with colored noise added.

5.2. Regulatory Control Profile. Four regulatory control loops, which are PID-type controllers, are modeled in the simulator. Level Control (LC Loop). This control loop adjusts the level of the disengaging-stripping section of the reactor by manipulating the slide valve position on the spent catalyst transport line. The LC loop has a direct effect on the riser temperature since the temperature depends on the flow rate of the catalyst. Pressure Difference Control (PDC Loop). This control loop adjusts the pressure difference between the regenerator and the disengaging-stripping section by manipulating the valve position on the stack gas line. Feed Temperature Control (FTC Loop). This control loop adjusts the temperature of the feed prior to entering the riser by manipulating the flow of fuel to the furnace. Main-Fractionator Pressure Control (MFPC Loop). This control loop adjusts the pressure of the overhead drum of the fractionator by manipulating the valve position of the wet gas transport line to the wet gas compressor.

All four of the regulatory controllers are highly interactive creating significant coupling in the FCC simulator resulting in an additional layer of colored noise. 5.3. DMC Controller Design. The control objective for the FCC process is to maintain the operation of the process at the specified set points in full combustion mode.2 The following controlled variables (CVs) and manipulated variables (MVs) are chosen for this work. CV1. The riser outlet temperature is controlled at a set point, which is selected on the basis of yield compromises (gasoline versus light gas yields) using an economic analysis. CV3. The stack gas oxygen concentration is controlled within a range or at set point to maintain the system in full combustion mode. MV4. Air blower suction valve. MV5. Slide valve position on the regenerator to reactor catalyst transport line. 5.4. The “True” FCC Model. The “true” model is defined as the best available model for this process. After some

6336


Figure 6. Diagram of a typical FCC process: (RX) reactor; (RG) regenerator; (MF) main-fractionator; (A) wet gas compressor; (B) overhead drum; (C) disengager; (D) stripper; (E) reactor riser; (F) cyclones; (G) feed preheater; (H) catalyst cooler; (I) stack gas expander; (J) air blower; (K) motor; (L) CO boiler [Han, 2004].

preliminary tests, the time to steady state was determined to be 60 min. No noise or disturbances were used during the plant test. Also, piecewise linearization27 on the process output data was performed to reduce the effects of nonlinearity. Therefore, the model obtained from this open-loop data set is the best model available. This model is called the “true” model of the FCC process in this work even though, in reality, no perfect model exists. This model is used to compare models identified from the proposed methodology to assess identification performance. 5.5. Results and Discussion. First, the effect of colored noise on the model identification is evaluated by adding colored noise to the CV measurements and collecting closed-loop data. Finally, a comparison between the proposed methodology and the twostage method is performed. 5.5.1. Effect of Colored Noise. Colored sensor noise is added to the CV measurements in this case. A FORTRAN subroutine developed by Fox28 and modified for this work is used to generate the colored noise on the FCC simulator. The noise model is based on Gaussian colored noise which is generated from Gaussian white noise with an exponential correlation function. The level of noise (variance) added to the CVs for this case is as follows; (2 × 10-3 to CV3 and (1.0 K to CV1. Figures 7 and 8 show the identified SRMs from closed-loop data using the DMCPlus controller for this case using the normal least-squares approach and the MX-FIR method, respectively. The dynamics and the steady-state gains of the FCC process are modeled accurately from the MX-FIR method. The average percent and maximum percent errors in the process gain are 1.1% and 2.3%, respectively. The expectations of the prediction errors, defined by eq 3, in this case, center on a zero mean which indicates an unbiased estimate. The calculated expectations are 9.69 × 10-6 and 5.80 × 10-3 for CV3 and CV1, respectively. The closed-loop test that was conducted for this case used a signal-to-noise ratio of 4 for CV1 and 8 for CV3.

Figure 7. Identified SRM from the closed-loop data using the traditional least-squares solution for noise levels of (2 × 10-3 for CV3 and (1.0 K for CV1: (s) SRM from closed-loop data, (- · -) “true” SRM.

The autocorrelation is calculated for the estimated model bias defined by eq 6. This bias term is designed to model the nonstationary components of the colored noise which eliminates any bias present in the model parameters due to external disturbances and nonlinearities. According to the results, the modeled bias term is nonstationary during the identification procedure. The autocorrelation ranges linearly from 1.0 to 0.8 for CV3 and 1.0 to 0.9 over 400 lags, where a lag is defined as a time shift. Therefore, the correlation of the training data has minimal affect on the process parameter estimates because the added bias model captures the dynamics of the additive nonstationary components. For the MX-FIR method, the calculated residual is assumed white. To test this assumption, the autocorrelation is calculated for the prediction errors, also known as the residuals. For a set


Figure 8. Identified SRM from closed-loop data using the MX-FIR method with colored noise added ((2 × 10-3 to CV3 and (1.0 K to CV1): (s) SRM from closed-loop data, (- · -) “true” SRM.

6337

Therefore, the prediction error at time equal t has the characteristics of white noise. Furthermore, the whiteness tests show that the bias associated with the correlated data is removed from the residual and process parameter estimates. 5.5.2. Comparison between the MX-FIR Method and the Two-Stage Method. In this section, the MX-FIR method is compared to a two-stage method13 which is available in the literature. First, an overview of the two-stage method is discussed and, then, the results are compared. The basis for the two-stage method is to obtain process input data that is uncorrelated with the process noise indirectly using an identified sensitivity function of the regulator. This method was first proposed by Van den Hof.13 An FIR-type model is used to identify a sensitivity function between the reference signal to the regulator and the process input since the process noise is uncorrelated with the reference signal. Predictions for the process input are calculated by projecting the reference signal onto the process input using the identified regulator model. Since this is a projection, the process input predictions are uncorrelated. The uncorrelated predictions are used in identifying the process model. The two-stage method is shown in the following steps. (1) Identify the sensitivity function that satisfies u(t) ) S(q-1) r(t) + e(t) where S(q-1) is the sensitivity function, r(t) is the reference signal, and u(t) is the process input. Then simulate u(t) ) Sˆ(q-1) r(t). (2) Identify the process model using a model of the form y(t) ) G(q-1, θ) uˆ(t) + e(t)

Figure 9. Sample autocorrelation of the calculated residual for both CVs to evaluate whiteness for noise levels of (2 × 10-3 for CV3 and (1.0 K for CV1.

where G(q-1,θ) is the process model estimate and uˆ(t) is the uncorrelated input. The above methodology is used to identify an SRM using training data generated by the FCC simulator with the DMCPlus controller and compared to an SRM identified using the MXFIR method. Figure 10 shows the identified SRMs using the MX-FIR method and the two-stage method. The two-stage method does improve the quality of the model compared to Figure 7. However, the SRM identified using the two-stage method still does not compare well to the “true” model. Therefore, the MXFIR method is shown to outperform the two-stage method for identification of process models using closed-loop data for this case study. 6. Conclusions

Figure 10. Comparison of the SRMs identified using the MX-FIR method (s) and the two-stage method ( · · · ) to the “true” model (- · -).

of data to be considered white, the autocorrelation of the first lag should be equal to one and at or close to zero elsewhere. Figure 9 shows the whiteness test for both CVs. According to this figure, the prediction errors have white characteristics since the moving average process, C(q-1) e(t), shown in eq 19, models the colored noise and whitens the prediction error due to the characteristics of the C polynomial. The first term in the C polynomial, C(0), is equal to one which means that C(0) is not estimated during the identification procedure. The other C(t 1), C(t - 2),.., C(t - n) parameters are uncorrelated with C(0).

A new methodology for closed-loop identification has been proposed. A noise model is introduced into the assumed model form to remove the effects of the parameter bias where a recursive algorithm is used to find the identification solution. This proposed method whitens and shifts the prediction error to a zero mean. From a consistency proof, the new model is consistent as N goes to infinity. Also, the proposed identification algorithm is bounded according to a stability analysis. The derived bias distribution for the new model shows that the process model bias is uncorrelated with the noise spectrum which is desirable for closed-loop identification. The proposed methodology is tested using two test models found in the literature. From the results, colored noise is the primary cause of parameter bias in closed-loop identification. However, accurate models are produced for all simulation case studies using the MX-FIR algorithm. The effect of relatively

6338


large increases in the colored noise resulted in relatively small increases in the errors in the identified models. Finally, a detailed simulator of an industrial FCC process is used to evaluate the performance of the proposed methodology. The FCC simulator uses colored noise on the process measurements as well as regulatory PID loops, which also contribute to the level of colored noise on the process measurements. The bias term models the nonstationary components of the colored noise which is seen by analysis of the autocorrelation results. Also, the results show the prediction error is whitened and shifted to a zero mean for all simulation case studies. Finally, the MX-FIR method is shown to outperform an existing technique for closed-loop identification of FIR models. Acknowledgment This work was funded by the Texas Tech University Process Control and Optimization Consortium. The authors would like to thank Aspen Technology and Cutler Technology for permission to use their software products. Literature Cited (1) Riggs, J. B.; Karim, M. N. Chemical and Bio-Process Control, 3rd ed.; Ferret Publishing: Lubbock, TX, 2006. (2) Han, I. S.; Riggs, J. B.; Chung, C. B. Multivariable control of a fluidized catalytic cracking process under full and partial combustion modes. J. Chem. Eng. Jpn. 2002, 35, 830–9. (3) Karimi, A.; Landau, I. D. Comparison of the closed-loop identification methods in terms of the bias distribution. Syst. Control Lett. 1998, 34, 159–67. (4) Zhu, Y. MultiVariable System Identification; Pergamon: Amsterdam, 2001. (5) Forssell, U.; Ljung, L. Closed-loop identification revisited. Automatica 1999, 35, 1215–41. (6) Klerk, E. de; Craig, I. K. Closed-loop system identification of a MIMO plant controlled by a MPC controller. IEEE AFRICON Conf. 2002, 1, 85–90. (7) MacGregor, J. F.; Fogal, D. T. Closed-loop identification: the role of the noise model and prefilters. J. Process Control 1995, 5, 163–71. (8) Gopaluni, R. B.; Patwardhan, R. S.; Shah, S. L. The nature of data pre-filters in MPC relevant identification-open- and closed-loop issues. Automatica 2003, 39, 1617–26. (9) Zheng, W. X. A Bias-correction method for indirect identification of closed-loop systems. Automatica 1995, 31, 1019–24. (10) Zheng, W. X. Modified least-squares identification of linear systems with noisy input and output observations. Proc. 35th Conf. Decision Control 1996, 1067-8. (11) Zheng, W. X. Identification of closed-loop systems with low-order controllers. Automatica 1996, 32, 1753–57.

(12) Zheng, W. X. A bias correction method for identification of linear dynamic errors-in-variables models. IEEE Trans. Autom. Control 2002, 47, 1142–7. (13) Van den Hof, P. M. J.; Schrama, R. An indirect method for transfer function estimation from closed-loop data. Automatica 1993, 46, 1523–7. (14) Forssell, U.; Ljung, L. A projection method for closed-loop identification. IEEE Trans. Autom. Control 2000, 45, 2101–06. (15) Landau, I. D.; Karimi, A. Recursive Algorithms for Identification in Closed Loop: A Unified Approach and Evaluation. Automatica 1997, 33, 1499–1523. (16) Isermann, R.; Lachmann, K.-H.; Matko, D. AdaptiVe Control Systems; Prentice Hall: New York, 1992. (17) Landau, I. D.; Karimi, A. A recursive algorithm for ARMAX model identification in closed loop. IEEE Trans. Autom. Control 1999, 44, 840– 3. (18) Ljung, L. System Identification Theory for the User; Prentice Hall: NJ, 1987. (19) Wood, R. K.; Berry, M. W. Terminal composition control of a binary distillation column. Chem. Eng. Sci. 1973, 28, 1707–17. (20) Ogunnaike, B. A.; Ray, W. H. Process Dynamics, Modeling, and Control; Oxford University Press: New York, 1994. (21) Han, I.; Chung, C. B.; Riggs, J. B. Modeling and optimization of a fluidized catalytic cracking process under full and partial combustion modes. Chem. Eng. Process. 2004. (22) Han, I.; Chung, C. B. Dynamic modeling and simulation of a fluidized catalytic cracking process. Part I: Process modeling. Chem. Eng. Sci. 2001, 56, 1951–71. (23) Garbow, B. S.; Hillstrom, K. E.; More, J. J. MINIPACK Project ; Argonne National Laboratory: Argonne, IL, 1980. (24) Riggs, J. B. An Introduction to Numerical Methods for Chemical Engineers, 2nd ed.; Texas Tech University Press: Lubbock, TX, 1994. (25) Aspen Technology, Inc. DMCPlus TM MultiVariable Control Software; Aspen Technology, Inc.: Houston, TX, 1999. (26) Aspen Technology, Inc. Cim-IO User’s Manual; Aspen Technology, Inc.: Houston, TX, 1998. (27) Dong, H.; Riggs, J. B. Detection of unmeasured disturbances in model predictive control (MPC) plant test data. Ind. Eng. Chem. Res. 2008, 47, 332–43. (28) Fox, R. F.; Gatland, I. R.; Roy, R.; Vemuri, G. Fast, accurate algorithm for numerical simulation of exponentially correlated colored noise. Phys. ReV. A 1988, 38, 5938–40. (29) Esmaili, A.; MacGregor, J. F.; Taylor, P. A. Direct and two-step methods for closed-loop identification: a comparison of asymptotic and finite data set performance. J. Process Control 2000, 10, 525–37. (30) Bauer, M.; Thornhill, N. F. A practical method for indentifying the propagation path of plant wide disturbances. J. Process Control 2008, 18 (7-8), 707–19.

ReceiVed for reView June 11, 2008 ReVised manuscript receiVed April 26, 2009 Accepted April 28, 2009 IE800922W

A Modified Extended Recursive Least-Squares Method for Closed

Recommend Documents