A Heuristic Extended Kalman Filter Based Estimator for Fault

May 31, 2003 - In the present study, an extended Kalman filter (EKF) based estimator for a complex ... filtering theory to nonlinear systems by linear...
1 downloads 0 Views 202KB Size
Download PDF
Ind. Eng. Chem. Res. 2003, 42, 3361-3371

3361

A Heuristic Extended Kalman Filter Based Estimator for Fault Identification in a Fluid Catalytic Cracking Unit Yuanjie Huang, G. V. Reklaitis, and Venkat Venkatasubramanian* School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47907-1283

In the present study, an extended Kalman filter (EKF) based estimator for a complex chemical process, namely, the Amoco model IV fluid catalytic cracking unit (FCCU), is investigated. This model is multivariable, strongly interacting, and highly nonlinear and is represented by an index one differential and algebraic equation system. The EKF has been modified to be able to handle algebraic state variables. A heuristic using pseudomeasurements is presented to reduce the model linearization errors in the EKF implementation. The performance of the estimator when applied to the Amoco model IV FCCU case study is presented in terms of early fault detection, speed, and estimation accuracy. The results show that, by introduction of pseudomeasurements, the accuracy and robustness of the estimator are improved significantly. 1. Introduction Fault diagnosis is an important aspect of process engineering. It is important not only from a safety viewpoint but also for the maintenance of yield and quality of a process. Consequently, the area has received significant attention from industry and academia alike because of the safety and economic impact involved. One review estimates that the loss to petrochemical industries in the U.S. alone is $20 billion/year.1 The problem of fault diagnosis and subsequent control is made much more difficult by the scale and complexity of modern plants. Typically, the fault diagnosis procedure has three tasks: (1) fault detection indicates that something is going wrong in the monitored system; (2) fault isolation determines the location of the fault; and (3) fault identification determines the magnitude of the fault. After a fault is detected and diagnosed, then in some applications, it is required that the fault is selfcorrected, usually through control reconfiguration. This is referred to as fault accommodation. The existing techniques for fault detection, isolation, and identification (FDI) can be broadly divided into process history based and process model based methods. Each of these can further be classified into qualitative and quantitative methods. There exist a number of techniques aimed at accurate and timely FDI utilizing state and parameter estimation techniques. Extensive reviews on these methods can be found in refs 2-5. State and parameter estimation methods for nonlinear systems are briefly discussed here. The estimation methods for nonlinear systems can be roughly classified into the following three categories: 1. Derived through direct extension of the linear filtering theory to nonlinear systems by linearization. Examples are the extended Kalman filter (EKF) and the extended Luenberger observer.6,7 2. Derived through some nonlinear state transformations using Lie algebra. These include the geometric observer (GO)8-11 and the high-gain (HG) observer.12-14 3. Formulated as a parameter estimation problem using optimization techniques.15-17 * To whom correspondence should be addressed. Tel.: (765)4940734. Fax: (765)4940805. E-mail: [email protected].

The GO is obtained by transforming a nonlinear state space model into one that is linear in the state variables but nonlinear in the input and output models under certain conditions. For this linearized state variables model, a linear observer is then designed. The GO guarantees convergence with linear output error dynamics but applies to an extremely restrictive class of nonlinear processes.18 The HG observer has been applied to biological and free-radical polymerization reactors.19,20 Under certain conditions, the HG observer has guaranteed convergence, but it involves a complex tuning procedure. By formulation of state and parameter estimation as a minimization problem, state estimates can be obtained by solving an online optimization problem, such as minimizing the sum of weighted squared errors.15,16 This method is statistically equivalent to the EKF,17 but it allows one to include constraints in the estimator. Because a dynamic optimization problem needs to be solved in every sampling period, this method may be difficult to implement and is computationally prohibitive for online use at the present time. Among the existing quantitative model-based methods, the EKF has found widespread use and is one of the most popular methods because of its simplicity in implementation and its ability to handle a reasonable degree of nonlinearities.21-23 There are, however, a number of practical challenges in designing such systems as a result of several factors such as the complexity of process dynamics, lack of adequate models, incomplete and uncertain data, diverse sources of knowledge, and amount of effort and expertise required to develop and maintain them. Hence, in the chemical engineering community, application of the EKF has only been reported for small (in terms of complexity and scale) case studies, most often for the continuously stirred tank reactor and the bioreactor.4,24-30 For large deviations from the reference state trajectory, the EKF performs poorly or becomes unstable. In the present study, first, various issues encountered in the application of the EKF to a complex process, the Amoco model IV FCCU,31 are investigated. This case study is chosen to evaluate the practical feasibility of the approach in terms of speed, accuracy, and computational complexity not only because it is a highly

10.1021/ie010659t CCC: $25.00 © 2003 American Chemical Society Published on Web 05/31/2003

3362

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

nonlinear, strongly coupled, multivariable but also because it has a significant economic impact. Second, a heuristic-based EKF (HEKF) using pseudomeasurements is presented, applied to the same process, and shown to improve the performance of the estimator significantly. The paper is structured as follows. First, a brief review of the existing model-based quantitative techniques for FDI is provided. Next, the EKF and HEKF estimators are discussed. Then, the model IV FCCU is briefly described, followed by discussions on the simulations of some typical cases including both single- and multiple-fault scenarios. The paper ends with conclusions. 2. EKF-Based Estimator The EKF, a heuristic filter based on the linearized dynamics of a system, has become the standard for state estimation of nonlinear systems. Many successful applications to chemical processes have been reported.24,25,28,29 The EKF assumes that (1) deviations from the reference state trajectory are small, (2) the mathematical description of the system dynamics and observations is accurate, and (3) the conditional probability density function of the state is Gaussian. For large deviations from the reference state trajectory, the EKF performs poorly or becomes unstable. Some studies have shown that the EKF is inadequate for certain processes.18,29,32 There has been extensive work on alleviating this problem. One of the approaches, called the iterated EKF, iteratively relinearizes the model equations about the current estimate.33-35 The idea is that, by locally making each estimate converge, the performance of the EKF will eventually be improved because the linearization error will be somewhat reduced. Consider the following nonlinear model:

xk+1 ) f(xk,uk) + k

(1)

yk ) h(xk) + νk

(2)

where xk, uk, and yk are state, input, and output vectors, respectively, and k and νk are zero-mean white noise. Assume that the zero-mean noise and the initial state estimation error are described by

{[ ]

} [ ]

 Q 0 δ E νk [jT νjT ] ) 0 R kj k

(3)

E{e0e0T} ) P0

(4)

The EKF iteration constitutes the following equations:

Lk ) FkPkHkT[HkPkHkT + JkRJkT]-1

(5)

Pk+1 ) FkPkFkT + GkQGkT - Lk[HkPkHkT + JkRJkT]LkT (6) xk+1|k ) f(xk,uk)

(7)

xk+1 ) xk+1|k + Lk[yk - h(xk+1|k)]

(8)

where Fk ) ∂f(xk)/∂xk, Gk ) ∂xk+1/∂k, Hk ) ∂h(xk)/∂xk, and Jk ) ∂yk/∂νk. More information on the EKF for the above model can be found in books.6,36

We consider the following general nonlinear process model:

x3 ) f(x,z,u,d)

(9)

0 ) g(x,z,u,d)

(10)

y ) h(x,z)

(11)

where x ∈ Rnx is the differential state vector, z ∈ Rnz is the algebraic state vector, u ∈ Rnu is the input vector, d ∈ Rnd is the unmeasured disturbance and parameter vector, and y ∈ Rny is the measured output vector. For ease of implementation, u and d are assumed to be piecewise constants. The dynamic process model is composed of sets of differential and algebraic equations, where the differential equations describe the behavior of the system and the algebraic equations represent the thermodynamic and kinetic relationships of the system. A large family of real chemical processes can be described by eqs 9-11. The discrete version of eqs 9-11 can be written as follows:

xk+1 ) ˜f(xk,zk,uk,dk)

(12)

0 ) g(xk,zk,uk,dk)

(13)

yk ) h(xk,zk)

(14)

where function ˜f(xk,zk,uk,dk) denotes the state at time tk+1 obtained by integrating eqs 9 and 10 from tk to tk+1 with the initial conditions of xk and zk for given uk and dk. There are two strategies in implementing the EKF for the system (9)-(11). One strategy is as follows: the differential algebraic equation (DAE) system is first transformed into an ordinary differential equation (ODE) system by differentiating eq 10. Then the state vector is augmented as x˜ ) [xT, zT]T. The discrete version of the resulting model has the form of eqs 1 and 2. The standard EKF iterations (5)-(8) can now be directly implemented. Using this approach, the differential and algebraic state variables in the original system are both updated simultaneously in each iteration. However, by the combination of x and z into one state vector, the dimension of the state vector will be increased dramatically for most of the process models in the chemical industry. This is not favorable from the numerical point of view because, in each EKF iteration step, a highdimensional matrix needs to be inverted. The alternative approach is that both eqs 12 and 13 are used directly in the implementation and the EKF only updates the differential state vector x, while z’s are updated recursively using the updated x’s. In this approach, the estimator is implemented for the state space with dimension nx. Matrix inversion only takes place in this reduced space of the original DAE system. This reduced space update approach is adopted in this work. The underlying dynamic system from which eqs 12 and 13 are derived is assumed to be an index one DAE system. This implies that the Jacobian matrix ∂g/∂zk is nonsingular. To update zk, the nonlinear algebraic equation (13) needs to be solved at each iteration. In addition, some of the algebraic variables are measured, and therefore these measurements should ideally also be used for state and parameter estimation. It is obvious

Ind. Eng. Chem. Res., Vol. 42, No. 14, 2003 3363

that the errors in the estimates and measurements are used as feedback to the estimator with respect to eq 12; the fact that these errors are also used as feedback to the estimator with respect to eq 13 is not that obvious. The reason goes as follows: because of various errors in the estimates and measurements, they will not satisfy eq 13 in general. However, by our approach which solves z from this equation, all estimation errors and measurement errors are preserved in the obtained z values, which, in turn, are used in the estimator update steps through eq 12. Measurement function h(xk,zk) is assumed to be a linear function in xk and zk. If h(xk,zk) is nonlinear, then by introduction of additional algebraic variables, h(‚) can be easily transformed into a linear function. 2.1. Gradient Matrix Evaluation. The EKF utilizes the following gradient matrices. The gradient of xk+1 to xk is

Fk )

∂xk+1 ∂f˜ ∂f˜ ∂zk ) + ∂xk ∂xk ∂zk ∂xk

[]

[

]

I ∂f˜ ∂f˜ ∂z k ) ∂x ∂z k k ∂xk

(15)

dz/dl ) φ(z,l,p), for 0 e l e L

∂yk ∂h ∂h ∂zk ) + ∂xk ∂xk ∂zk ∂xk

[

[]

]

I ∂h ∂h ) ∂x ∂z ∂zk k k ∂xk

( )

-1

∂g ∂xk

(19)

∂ψ[z(l)L)]/∂p From the chain rule,

(16)

where ∂h/∂xk and ∂h/∂zk are constant matrices because function h(‚) is linear. Matrix ∂zk/∂xk can be computed using the following formula:

∂zk ∂g )∂xk ∂zk

(18)

with initial condition z(l)0) ) z0, where p is a parameter vector. For a given function ψ(z), the goal is to evaluate the sensitivity:

The gradient of yk to xk is

Hk )

availability of cheap and easy to use parallel computing systems, a direct perturbation- and integration-based method is becoming feasible even for large-scale problems. In addition, at each iteration, the algebraic variable zk needs to be determined from eq 13. Usually the model structure can be exploited such that most of the algebraic variables can be computed in an a priori specified order. For those that must be solved simultaneously, a Newton-Raphson type of second-order nonlinear algorithm can be used because here a good estimate can be obtained as the initial guess for the solver. If the process model contains a distributed subsystem, then a few of the variables in eq 13 are implicitly governed by a set of ODE/partial differential equation (PDE) or DAE in spatial dimension(s). Such a situation is often seen in chemical process models. For such cases, ∂zk/∂xk cannot be computed using eq 17. Methods such as the finite difference approximation discussed next can be used. 2.2. Calculation of ∂zk/∂xk in a Distributed System. This procedure is explained using the following ODEs:

(17)

Matrix Fk measures the sensitivity of the DAE system (12) and (13) with respect to the initial condition xk. The calculation of Fk is worth further elaboration because the accuracy of Fk is directly related to the accuracy of the first-order Taylor expansion of the nonlinear model. In general, there is no closed form for function ˜f(‚), which means that the Jacobian matrices ∂f˜/∂xk and ∂f˜/∂zk in eq 15 can only be computed numerically. The methods used to solve the combined DAE system and these sensitivity matrices can be classified into two categories: (i) the staggered direct scheme in which the linear systems for the sensitivity corrector steps are solved directly after the convergence of the nonlinear corrector step;37-39 (ii) the simultaneous corrector method which solves the DAE system and its sensitivity equations simultaneously, obtaining the solution from a corrector iteration on the combined nonlinear system40 or two corrector iterations, one for the DAE system and the other for the sensitivities.41 From the computational point of view, simultaneous methods are better than the staggered scheme because the number of matrix factorizations that need to be performed along the solution trajectory is minimized. With the increasing power of modern computers and the

∂ψ[z(l)L)] ∂ψ ∂z ) ∂p ∂z ∂p

|

(20)

l)L

Thus, we only need to compute ∂z(l)L)/∂p. Suppose z ) φ(l,p) is the solution to eq 18; then ∂z(l)L)/∂p can be approximated by

∂z ∂pi

|

l)L

≈

φ(L,p+∆pei) - φ(L,p) ∆p

(21)

i ) 1, 2, ..., N where N is the number of elements in vector p, ei is the unit vector with all elements equal to 0 except the ith element, which is 1. Using this method to compute ∂zk/ ∂xk, eq 18 must be integrated N + 1 times. In the case of large N, all integrations should be performed simultaneously to achieve efficiency. 2.3. Pseudomeasurements Strategy. In the EKF, the state distribution is approximated by a Gaussian random variable, which is then propagated analytically through “first-order” linearization of the nonlinear system. As such, the EKF can be viewed as providing first-order approximation for the optimal terms in the Kalman filter of linear systems. These approximations, however, can introduce large errors in the true posterior mean and covariance of the transformed random variable, which may lead to a suboptimal estimate and sometimes even divergence of the filter. A closer look at the linearization process indicates that the linearization error consists of the error due to the deviation of the current estimate from the reference state and the error caused by the nonlinearity of the system. Because of the iterative nature of the EKF, the linearization

3364

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

Figure 1. Pseudomeasurements and linearization errors.

error is mostly caused by the error in the initial state estimate and the nonlinear structure of the model. The error introduced in the process of linearization is one of the key reasons that makes EKF inadequate for some processes. Any strategy that results in a smaller linearization error will potentially improve the performance of the EKF significantly. Ideally, if the sampling interval can be arbitrarily decreased, then the linearization error from model structures can be made arbitrarily small if the current estimate is accurate enough. On the basis of these observations, a pseudomeasurementsbased estimator is proposed. In each iteration, once a new measurement is available, this new measurement and the previous measurement are interpolated to obtain nM > 0 pseudomeasurements during the past sampling interval. These pseudomeasurements are treated as true measurements and used in the measurement update step. The ideas behind this are shown in Figure 1. Two measurements are taken at times t1 and t5, and pseudomeasurements are introduced at times t2, t3, and t4. Generally, the linearization errors E, e1, e2, e3, and e4 will have the following relationship:

e1 + e2 + e3 + e4 < E

(22)

There are a number of ways for obtaining the pseudomeasurements. The simplest method is using linear interpolation. Once such an interpolation strategy is chosen, the next decision is choosing the number of pseudomeasurements in each sampling period. Usually this number can be determined through simulation. Another interpolation method is to use a spline function, which generally gives a better approximation of process dynamics. For the case when the samples are not taken very often and where the process trend may be nonmonotonic in one sampling interval, a spline function instead of a linear function should be used for the interpolation. 3. FCCU Case Study In this section, an overview is given for the Amoco model IV FCCU. A detailed description of the FCCU is given in the original reference by McFarlane et al.31 The schematic of the FCCU is shown in Figure 2. The FCCU is essentially a catalytic reactor which cracks high boiling gas oil feed into lower boiling more valuable hydrocarbons in the presence of a catalyst. Most of this endothermic reaction takes place in the reactor riser.

The reactor riser temperature is very close to the metallurgical limits for optimum production. As a result of cracking, carbonaceous products (called coke) get deposited on the catalyst, which decreases the effectiveness and the lifetime of the expensive catalyst. This is continuously regenerated in the regenerator by blowing in air. Coke is combusted to CO, CO2, and H2O. The amount of CO vented out through the stack gas is very crucial from an environment point of view, and one of the challenges for the model IV FCCU is not to violate the environmental threshold and to achieve an optimum performance in the face of disturbances. There is a constant flow of regenerated and spent catalyst between the reactor and the regenerator, as shown in Figure 2. This flow is partly driven by the pressure differential between the reactor and regenerator, and the remaining momentum is supplied by the lift air blower. The fractionator section separates the product hydrocarbons for further processing and is not modeled in detail. A feed system consisting of low-level flow controllers and a preheating furnace preprocess the feed for cracking. A simulator of the model IV FCCU, called CATSIM,42 which is based on the model equations given in work by McFarlane et al.,31 is used in this work. The nomenclature and nominal operating conditions used in this work are identical to those discussed by McFarlane et al. The following three low-level proportionalintegral (PI) controllers (not shown in Figure 2) are essential for stable operation of the FCCU: (1) antisurge controller, which maintains a constant total air flow (combustion air + lift air + spill air) to the regenerator; (2) reactor pressure controller, which maintains a constant reactor pressure; (3) reactor/regenerator differential pressure controller, which maintains a constant pressure differential. CATSIM can simulate 87 different faults. These include 9 disturbances, 18 parameter faults, 8 flow regulatory controller faults, 34 sensor biases, and 18 regulatory controller faults. The severity range of each of the faults was in the range of 5-15% in both the positive and negative direction. There are 20 measured variables in the process, and their response time varied from instantaneous in some cases to 30 min in others. Gaussian noise with zero-mean and known standard deviations are added to the measured values. The nominal value and standard deviation of each measured variable are shown in Table 1. The standard deviations for the measured variables have been chosen to give a reasonable S/N ratio and based on the typical noise levels encountered in the instruments used to measure them. In CATSIM, the three PI controllers adjust their outputs every second to achieve stable and smooth operation, while the measurements are taken every 30 s. The FCCU dynamic model can be expressed as a DAE system (9)-(11). The model used for state and parameter estimation consists of 17 differential state variables, 20 outputs, and 85 algebraic variables. Among these 85 algebraic variables, the four in the regenerator, xCO,sg, xCO2,sg, xO2,sg, and Tcyc, are governed by a set of DAEs in spatial dimension. The integration of this distributed subsystem is the most time-consuming step in the estimation computation because of its large scale in spatial dimension. The EKF uses only real measurements, and the HEKF uses both real measurements and pseudomeasurements. Both approaches are implemented based on the discrete model (12)-(14). During each sampling period, 29 pseudomeasurements were

Ind. Eng. Chem. Res., Vol. 42, No. 14, 2003 3365

Figure 2. Schematic of Amoco model IV FCCU. Table 1. Nominal Value and Standard Deviation of Measured Variables in Model IV FCCU sensor index

measured variable

nominal value

std dev

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

combustion air flow to regen (F7) lift air flow to regen (F9) spill air flow to regen (F10) reactor pressure (P4) main frac pressure (P5) regen pressure (P6) regen stack gas temp at cyclone (Tcyc) regen temp (Treg) reactor riser temp (Tr) comb air blower inlet flow (Fsucncomb) lift air blower inlet flow (Fsucnlift) lift air blower surge flow (Fsurgelift) O2 concentration in stack gas (CO2,sg) CO concentration in stack gas (CCO,sg) actual speed of lift air blower (Sa) level of catalyst in standpipe (Lsp) wet gas compressor inlet flow (Fsucnwg) temp of fresh feed (T2) furnace firebox temp (T3) air flow into regen (Ft)

61.15 (lb/s) 14.31 (lb/s) 0 (lb/s) 33.02 (psia) 23.52 (psia) 29.64 (psia) 1283.42 (°F) 1271.79 (°F) 994.99 (°F) 49536.23 (ICFM) 11547.59 (ICFM) 9534.78 (ICFM) 1.52 (mol %) 67.45 (ppm) 5442.42 (rpm) 11.29 (ft) 19067.87 (ICFM) 667.26 (°F) 1607.55 (°F) 75.46 (lb/s)

0.01 0.03 0.05 0.02 0.03 0.05 0.06 0.06 0.06 4.5 4.5 0.8 0.01 0.08 0.6 0.08 4.5 0.08 0.08 0.04

introduced. In total, more than 20 faults are simulated and correctly identified. Several typical case runs are discussed next; these include both single- and multiplefault scenarios. Unless explicitly pointed out, all faults are introduced at the 10th sample time because the process is almost at steady state at that time. In all of the case runs, no random noise was added to the state equations, and the noise was only added to the measurements. Thus, the covariance matrix Q was set to 0. The R matrix was estimated using a wavelet threshold shrinkage method.43 Case 1: Coke Factor. A 5% increase in the coke factor (ψf) of the feed was simulated. This disturbance affects the FCCU severely and leads to violation of almost all operating constraints. More coke in the feed quickly leads to an increase in the weight fraction of coke on spent catalyst (Csc) in the reactor. More coke enters the regenerator, and because the total air flow is held constant, insufficient O2 is available for combustion, resulting in a higher stack CO concentration. A

higher reactor temperature (due to higher Treg) increases the wet gas production. The reactor pressure (P4) is held constant by opening of valve V11, which saturates quickly, and then P4 increases. The differential pressure controller increases the regenerator pressure (P6), which reduces the throughput of the combustion air blower. There is a large time delay before the disturbance affects the blower system. Quick detection is of utmost importance in this fault. The EKF both using and without using the pseudomeasurements are applied to detect the fault, the coke factor (psif). Both methods are able to estimate the correct magnitude of the fault. However, the HEKF tracks down the fault in less than 10 sampling periods. This is much less than 40 sampling periods taken by the EKF. The EKF fails to converge for the state variable Wreg, while the HEKF successfully captures the true value of Wreg. The simulation results are shown in Figures 3 and 4. Case 2: Heat-Transfer Coefficient. This corresponds to a 10% decrease in the furnace heat-transfer

3366

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

Figure 3. 5% increase in ψf, state & fault (EKF).

Figure 4. 5% increase in ψf, state & fault (HEKF).

coefficient (UAf). A decrease in the overall heat-transfer coefficient results in a decrease in the fresh feed temperature and an increase in the flue gas temperature leaving the furnace. Lower feed temperature decreases the endothermic reaction and leads to a decrease in the riser temperature. Asymptotically, both

approaches give a correct fault parameter (UAf) estimate. However, the estimate using the HEKF is much more accurate than the estimate obtained using the EKF. In terms of efficiency, the HEKF only needs about 6 sampling intervals to converge to the real value, while the EKF needs about 60 sampling intervals. The

Ind. Eng. Chem. Res., Vol. 42, No. 14, 2003 3367

Figure 5. 10% decrease in UAf, state & fault (EKF).

Figure 6. 10% decrease in UAf, state & fault (HEKF).

improvement in timely identification is about 1 order of magnitude in this case. Estimation results obtained from EKF and HEKF are shown in Figures 5 and 6. Case 3: Catalyst Particle Density. This corresponds to a 5% increase in the catalyst particle density (Fpart). An increase in Fpart increases the material density in the reactor (Fris), through a decrease in the riser velocity. An increase in Fris increases the reactor bottom

pressure (Prb). This leads to the sudden drop in the reactor riser temperature (Tr). Similarly, a reduction in Fris reduces the dense phase density (Fc,dense), which causes an increase in the pressure at the bottom of the lift pipe. This, in turn, leads to a sudden drop in the total air flow rate. This fault is only identified by the HEKF method. The EKF does not converge for both the state variable (Wreg) and the fault parameter (Fpart)

3368

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

Figure 7. 5% increase in Fpart, state & fault (EKF).

Figure 8. 5% increase in Fpart, state & fault (HEKF).

because of the sudden dramatic change in the state vector caused by the fault parameter change. The HEKF satisfactorily detects the magnitude of the occurring fault in about 4 sampling intervals, or 2 min. The estimation results are shown in Figures 7 and 8.

Case 4: Ambient Temperature and Coke Factor. In this run a multiple fault occurs: a 10% increase in the ambient temperature (Tatm) and a 5% increase in the coke factor (ψf). An increase in Tatm decreases F6, which decreases the combustion air blower discharge

Ind. Eng. Chem. Res., Vol. 42, No. 14, 2003 3369

Figure 9. 10% increase in Tatm and 5% increase in ψf, state & fault (HEKF).

Figure 10. 10% decrease in F45000 and 10% increase in furgc, state & fault (HEKF).

pressure (P2). The lift air blower discharge pressure also decreased. However, this only lasts for a very short duration because the FT controller quickly increases the lift air flow rate F9. This also causes Tr to drop suddenly and then quickly recover and begin to increase slowly.

Because in this case the coke factor fault is also taking place, the increase in Tr and Treg speeds up, which causes valve V11 to saturate sooner. Again the HEKF works well for this run, it converges in less than 10 sampling periods with a relative error to the change

3370

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

Table 2. Case Runs Summary EKF

HEKF

index

fault description

stability

timeliness (min)

stability

timeliness (min)

5 8 10 1 and 5 2 and 9

coke factor catalyst particle density furnace heat-transfer coefficient air blower capacity and coke factor ambient temperature and friction factor

converge diverge converge diverge diverge

20

converge converge converge converge converge

4 2 3 3 5

of the parameter at about 0.5%, while the EKF fails to converge to the true state. The result obtained using the HEKF is plotted in Figure 9. Case 5: Combustion Air Blower Capacity and Catalyst Friction Factor. In this final case, a 10% decrease in the combustion air blower capacity (F45000) and a 10% increase in the regenerated catalyst friction factor (furgc) are introduced. The sudden drop and quick recovery in P5 and P7 are captured by the HEKF. In about 10 sampling periods, or 5 min, the estimator converges to the true values of the fault parameters. The EKF diverges in this multiple-fault case. Figure 10 shows the estimation result computed using the HEKF. The simulations show that, by adding pseudomeasurements, the estimator becomes much more stable. The applicable class of process models is enlarged significantly. As for the two multiple-fault cases, the EKF diverges in both cases while the HEKF converges to the true fault values in both cases. One of the reasons that the EKF fails in both multiple-fault cases is because when multiple faults are introduced in the process, the sudden changes in the state variables are much more severe than those in single-fault cases. Because of the linearization nature of the EKF, the situation is even worse. It cannot recover from these dramatic variations. However, for the HEKF, the linearization effect is reduced significantly; this reduces the effect caused by the sudden variations in the state variables. In terms of accuracy and convergence rate, the performance of the HEKF is much better. For those two cases where both HEKF and EKF converge, the HEKF converges to the real value more quickly and with smaller errors. In general, the HEKF improves the convergence rate for about 1 order of magnitude. Statistics of these case runs are summarized in Table 2. 4. Conclusions In this paper, an EKF-based estimator was described in some detail. For a general nonlinear DAE system, the standard EKF is first extended to handle the algebraic variables in the system by partitioning the state variables into dependent and independent sets. Then, the EKF update steps are performed in the independent state (reduced) space, which greatly reduces the iteration computation time because a matrix inversion now only needs to be carried out in a much smaller dimensional space. This makes the approach suitable even for large-scale models. One of the main factors that leads to the failure of the EKF for certain nonlinear systems is the linearization process. To reduce the linearization errors introduced in the EKF procedure, a HEKF using pseudomeasurements was presented. The approach is based on the hypothesis that more measurements generally lead to better estimates. By introduction of more measurements, nonlinear models reduce to linear models locally around each measurement with higher accuracy. During

30

each sampling period, pseudomeasurements are inserted using simple function interpolation methods. These pseudomeasurements are then used by the estimator just like they are real measurements from sensors. An estimation step is performed on each measurement (pseudo and real). These approaches were applied to a highly nonlinear, tightly coupled system, namely, the Amoco model IV FCCU process. The results from the applications of the EKF and HEKF to both single- and multiple-fault scenarios showed that the EKF could detect only some of the single-fault scenarios and fails to converge for all of the multiple-fault cases. By contrast, the HEKF could detect the faults in all scenarios. The HEKF not only obtains the true estimate more quickly but also converges to more accurate state values. These case runs showed that the HEKF is a good candidate for largescale nonlinear model-based state and parameter estimation, especially for use in fault detection, isolation, and identification of complex systems. Nomenclature Csc ) weight fraction of coke on spent catalyst (lb of coke/ lb of catalyst) FT ) air flow rate into the regenerator (lb/s) furgc ) regenerated catalyst friction factor (lb/ft2) F6 ) combustion air blower throughput (lb/s) F9 ) lift air flow to the regenerator (lb/s) F45000 ) constant in the combustion air blower inlet flow equation (ICFM) Prb ) pressure at the bottom of the reactor riser (psia) psif ) effective coke factor for gas oil feed P2 ) combustion air blower discharge pressure (psia) P4 ) reactor pressure (psia) P5 ) reactor fractionator pressure (psia) P6 ) regenerator pressure (psia) P7 ) wet gas compressor suction pressure (psia) Tatm ) atmospheric temperature (°F) Tcyc ) regenerator stack gas temperature at cyclone (°F) Tr ) temperature of the reactor riser (°F) Treg ) temperature of the regenerator bed (°F) UAf ) furnace overall heat-transfer coefficient (Btu/s) V11 ) wet gas compressor suction valve position Wreg ) inventory of the catalyst in the regenerator (lb) xCO,sg ) molar ratio of CO to air in stack gas (mol of CO/ mol of air) xCO2,sg ) molar ratio of CO2 to air in stack gas (mol of CO2/ mol of air) xO2,sg ) molar ratio of O2 to air in stack gas (mol of O2/mol of air) Fc,dense ) density of the catalyst in the dense bed (lb/ft3) Fpart ) settled density of the catalyst (lb/ft3) Fris ) average density of the material in the reactor riser (lb/ft3) ψf ) effective coke factor for gas oil feed

Literature Cited (1) Nimmo, I. Adequately address abnormal situation operations. Chem. Eng. Prog. 1995, 91 (9), 36.

Ind. Eng. Chem. Res., Vol. 42, No. 14, 2003 3371 (2) Frank, P. M. Fault diagnosis in dynamic systems using analytical and knowledge-based redundancysA survey and some new results. Automatica 1990, 26 (3), 459. (3) Gertler, J. Analytical redundancy in fault detection and isolation - survey and synthesis. IFAC Symposium on Online Fault Detection and Supervision in the Chemical Process Industries (SAFEPROCESS ’91), Baden-Baden, Germany, Sept 10-13, 1991; IFAC: 1992; pp 9-21. (4) Soroush, M. State and parameter estimations and their applications in process control. Comput. Chem. Eng. 1998, 23, 229. (5) Gertler, J. J. Fault detection and diagnosis in Engineering systems; Marcel Dekker: New York, 1998. (6) Jazwinski, A. H. Stochastic Processes and Filtering Theory; Academic Press: New York, 1970. (7) Zeitz, M. The extended luenberger observer for nonlinear systems. Syst. Control Lett. 1987, 9, 149. (8) Nijmeijer, H. Observability of autonomous discrete-time nonlinear systems: A geometric approach. Int. J. Control 1982, 36, 867. (9) Bestle, D.; Zeitz, M. Canonical form observer design for nonlinear time-variable systems. Int. J. Control 1983, 38 (2), 419. (10) Krener, A. J.; Isidori A. Linearization by output injection and nonlinear observers. Syst. Control Lett. 1983, 3 (1), 47. (11) Xia, X.; Gao, W. Nonlinear observer design by observer error linearization. SIAM J. Control Optim. 1989, 27, 199. (12) Tsinias, J. Observer design for nonlinear systems. Syst. Control Lett. 1989, 13 (2), 135. (13) Gauthier, J. P.; Hammouri, H.; Othman, S. A simple observer for nonlinear systems applications to bioreactors. IEEE Trans. Autom. Control 1992, 37 (6), 875. (14) Deza, F.; Busvelle, E.; Gauthier, J. P.; Rakotopara, D. High gain estimation for nonlinear systems. Syst. Control Lett. 1992, 18 (4), 295. (15) Liebman, M. J.; Edgar, T. J.; Lasdon, L. S. Efficient date reconciliation and estimation for dynamic processes using nonlinear programming techniques. Comput. Chem. Eng. 1992, 16, 963. (16) Michalska, H.; Mayne, D. Q. Moving horizon observers and observer-based control. IEEE Trans. Autom. Control 1995, 40 (6), 995. (17) Robertson, D. G.; Lee, J. H.; Rawlings, J. B. A moving horizon-based approach for least-squares estimation. AIChE J. 1996, 42 (8), 2209. (18) Kantor, J. C. Finite dimensional nonlinear observer for an exothermic stirred-tank reactor. Chem. Eng. Sci. 1989, 44 (7), 1503. (19) van Dootingh, M.; Viel, F.; Rakotopara, D.; Gauthier, J. P.; Hobbes, P. Nonlinear deterministic observer for state estimation: application to a continuous free radical polymerization reactor. Comput. Chem. Eng. 1992, 16 (8), 777. (20) Gibon-Fargeot, A. M.; Hammouri, H.; Celle, F. Nonlinear observers for chemical reactors. Chem. Eng. Sci. 1994, 49 (14), 2287. (21) Watanabe, K.; Himmelblau, D. M. Fault diagnosis in nonlinear chemical processspart I. Theory. AIChE J. 1983, 29 (2), 243. (22) Watanabe, K.; Himmelblau, D. M. Fault diagnosis in nonlinear chemical processspart II. Application to a chemical reactor. AIChE J. 1983, 29 (2), 250. (23) Chang C. T.; Hwang, J. I. Simplification techniques for EKF computations in fault diagnosis: Model decomposition. AIChE J. 1998, 44 (6), 1392.

(24) Jo, J. H.; Bankoff, S. G. Digital monitoring and estimation of polymerization reactors. AIChE J. 1976, 22 (2), 361. (25) Schuler, H.; Zhang, S. Z. Real-time estimation of the chain length distribution in a polymerization reactor. Chem. Eng. Sci. 1985, 40 (10), 1891. (26) Ramirez, W. F. Optimal state and parameter identification: an application to batch fermentation. Chem. Eng. Sci. 1987, 42, 2749. (27) Choi, K. Y.; Khan, A. A. Optimal state estimation in the transesterification stage of a continuous poly(ethylene terephthalate) condensation polymerization process. Chem. Eng. Sci. 1988, 43, 749. (28) Kim, K. J.; Choi, K. Y. On-line estimation and control of a continuous stirred tank polymerization reactor. J. Process Control 1991, 1, 96. (29) Kozub, D. J.; MacGregor, J. F. State estimation for semibatch polymerization reactors. Chem. Eng. Sci. 1992, 47 (5), 1047. (30) Sriniwas, G. R.; Arkun, Y.; Schork, F. J. Estimation and control of an alpha-olefin polymerization reactor. J. Process Control 1995, 5, 303. (31) McFarlane, R. C.; Reineman, R. C.; Bartee, J. F.; Georgakis, C. Dynamic simulator for a model IV fluid catalytic cracking unit. Comput. Chem. Eng. 1993, 17 (3), 275. (32) Gudi, R. D.; Shah, S. L.; Gray, M. R. Adaptive multirate state and parameter estimation strategies with application to a bioreactor. AIChE J. 1995, 41 (11), 2451. (33) Ljung, L. Asymptotic behavior of the extended kalman filter using filtered state estimates. IEEE Trans. Autom. Control 1979, 24 (1), 36. (34) Galkowski, P. J.; Islam, M. A. The iterated kalman filter update as a gauss-newton method. IEEE Trans. Autom. Control 1991, 36 (11), 1323. (35) Bell, B. M. The iterated kalman smoother as a gaussnewton method. SIAM J. Optim. 1994, 4 (8), 626. (36) Goodwin, G. C.; Sin, K. S. Adaptive Filtering Prediction and Control; Prentice Hall: Englewood Cliffs, NJ, 1984. (37) Caracotsios, M.; Stewart, W. E. Sensitivity analysis of initial value problems with mixed ODEs and algebraic equations. Comput. Chem. Eng. 1985, 9 (4), 359. (38) Haug, E. J.; Ehle, P. E. Second-order design sensitivity analysis of mechanical system dynamics. Int. J. Numer. Methods Eng. 1982, 18 (1), 1699. (39) Rabitz, H.; Kramer, M.; Dacol, D. Sensitivity analysis in chemical-kinetics. Annu. Rev. Phys. Chem. 1983, 34, 419. (40) Maly, T.; Petzold, L. R. Numerical method and software for sensitivity analysis of differential-algebraic systems. Appl. Numer. Math. 1996, 20, 57. (41) Feehery, W. F.; Tolsma, J. E.; Barton, P. I. Efficient sensitivity analysis of large-scale differential-algebraic systems. Appl. Numer. Math. 1997, 25 (1), 41. (42) Mylaraswamy, D. Reference Manual for CATSIM, 2.1 ed.; LIPS, Purdue University: West Lafayette, IN, 1995. (43) Huang, Y. J.; Reklaitis, G. V.; Venkatasubramanian, V. Wavelet shrinkage based measurement error covariance matrix estimation from process measurements. Proceedings of ESCAPE11, Kolding, Denmark, 2001.

Resubmitted for review June 21, 2002 Revised manuscript received September 25, 2002 Accepted April 18, 2003 IE010659T