Sequential estimation for nonlinear differential and algebraic systems—theoretical development and application

Pergamon

Computerschem. EngngVol.21, No. 9, pp. 1051-1067,1997 © 1997ElsevierScienceLtd Printedin GreatBritain.All rightsreserved PII: S0098-1354(96)00335-3 0098-1354/97$17.00+0.00

Sequential estimation for nonlinear differential and algebraic systems-theoretical development and application Y. S. Cheng, T. Mongkhonsi* and L. S. Kershenbaum# Centre for Process Systems Engineering, and Department of Chemical Engineering and Chemical Technology, Imperial College of Science, Technology and Medicine, London, SW7 2BY, U.K.

(Received 9 March 1996)

Abstract Based on a minimum least squares criterion and without any assumption on process model error and measuring device noise, a sequential estimation algorithm has been developed for nonlinear differential and algebraic systems with the use of variational calculus. Estimators for both continuous and discrete time data are derived, and can be used without explicitly considering the indices of differential and algebraic systems. The estimation algorithm has been applied to a large pilot scale fixed-bed reactor with decaying catalysts for simultaneous estimation of catalyst activity, temperature and concentration profiles. Although a simplified lumped reaction kinetics network and an approximate catalyst deactivation model were utilised by the estimator, the catalyst activity profile was accurately inferred and the dynamic behaviour of the reactor was excellently tracked using only two or three real-time temperature measurements. The estimated variation of the catalyst activity enabled us to fully understand the change of the hot spot magnitude and the shift of the temperature profile. Computational experience indicated that the estimator can be implemented on-line for reactor control and optimisation. © 1997 Elsevier Science Ltd Introduction Differential and algebraic systems, sometimes also referred to as generalised state space or descriptor systems, are characterised by a set of mixed differential and algebraic equations. The differential equations represent the dynamics of the system and the algebraic equations usually describe the constraints among the variables. Differential and algebraic systems have been studied by a number of investigators (Luenberger, 1977; Luenberger, 1978; Verghese et al., 198 l), and are widely used for the dynamic modelling of chemical process (Pantelides et al., 1988; Pantelides and Barton, 1993). The state estimation problem for linear differential and algebraic systems has been studied by several researchers (Nikoukhah et al., 1992; Chisci and Zappa, 1992; Darouach and Zasadzinski, 1992). However, no progress seems to have been made in developing an estimation algorithm for nonlinear differential and algebraic systems. Although, in many cases, the algebraic equations can be eliminated by substitution, it is * Present address: Departmentof Chemical Engineering,ChulalongkornUniversity,Bangkok,Thailand f Author to whom correspondence should be addressed. Fax: (+44) 171 594 6606. E-mail:[email protected]

not always convenient to do so and, therefore, it would be useful to have an estimator for nonlinear differential and algebraic systems. In this paper we present a sequential estimation algorithm for general nonlinear differential and algebraic systems. There are many feasible criteria which could be used for estimation, such as minimal least squares, maximum likelihood and minimum maximum error. In this study the minimum least squares criterion, i.e., the minimisation of the integral of the weighted square residual errors in the process model and the measuring device, is employed. This criterion has the advantage of not requiring knowledge of the nature of the unknown inputs to the system and the measurement error of the output. The estimation problem is reformulated as a deterministic optimal control problem, a "formal" approach already used in ordinary differential systems (Detchmendy and Sridhar, 1966) and partial differential systems (Seinfeid, 1969). The calculus of variations for differential and algebraic systems (Jonckheere, 1988) is used to obtain a sequential estimation algorithm. The estimation algorithm is applicable to cases in which measurement data are available either continuously or at discrete sampling instants. It should be noted that the application of the nonlinear programming (NLP) approach based on a moving

1051

1052

Y. S. CHENGet al.

horizon for state and parameter estimation has been widely reported (Kim et al., 1990; Liebman et ai., 1992; Ramamurthi et al., 1993; Robertson and Lee, 1995; Cheng et al., 1996a). By comparison with the extended Kalman filter (EKF), the main advantages of the NLP approach are its robustness and its capability of handling inequality constraints including variable bounds. On the other hand, for on-line implementation, the computing time required by the NLP approach to obtain state and parameter estimates becomes prohibitive with growing dimension as more measurement data are collected (Cheng et al., 1996a). By contrast, the estimation algorithm developed here can be implemented recursively and solved in real-time for on-line applications; this will be demonstrated in the case study of a pilot scale fixed-bed reactor with catalyst deactivation. Fixed-bed reactors are highly nonlinear distributed parameter systems. The theory of state estimation for distributed parameter systems was developed during the 1970s (Tzafestas, 1978), but the progress in application of these ideas to real systems has been slow (Ray, 1978). Perhaps the main reasons are the implementation complexity of distributed parameter estimation algorithms, the computing time and memory limit for solving these in real-time, and the lack of robust numerical methods. Two methodologies, early lumping and late lumping, have been widely used for temperature and concentration estimation of fixed-bed reactors. In early lumping, the distributed parameter model is first converted to the lumped parameter one with the use of discretisation techniques and then the EKF is applied (Jutan et al., 1977; Sorenson, 1977; Clement et al., 1980; Wallman and Foss, 1981; Siliva et al., 1979). By contrast, the alternative late lumping applies distributed parameter estimation theory first and then solves the full distributed estimation equations by lumping approximation techniques (Windes et al., 1989). When deactivation of the catalyst occurs, the catalyst will generally have a variable activity which is a function of both space position and operating time. In almost all industrial reactors there is no means available to directly measure the catalyst activity. However, it may be possible to infer the catalyst activity from measurements of a small number of temperatures or/and concentrations within the reactor. Some researchers have estimated the catalyst activity in the reactors while assuming the temperature and concentration were at quasi-steady state and the catalyst activity was spatially uniform (Joffe and Sargent, 1972; Ajinka et al., 1974; Kurnoglu et al., 1981). In recent work, the authors have developed a new scheme for estimation of spatially varying dynamic catalyst activity and spatial distributions of quasi-steady state temperature and concentration (Cheng et al., 1993). In many cases, the transient behaviour of temperature can not be neglected; nor is the approximation of spatially uniform catalyst activity reasonable. In this paper we apply the sequential estimator for nonlinear differential and algebraic systems to simultaneously estimate the catalyst activity, temperature and concentra-

tion profiles in a fixed-bed reactor. The reactor studied here is a near-industrial scale (25 mm diameter x 3 m length), jacket-cooled, tubular reactor with twenty six thermocouples along the axial direction of the reactor for temperature measurements. The reactor is designed to carry out the partial oxidation of o-xylene to phthalic anhydride using a V2OsITiO2 catalyst. Because the catalyst can deactivate both reversibly and irreversibly (as described below), the estimation of catalyst activity would be extremely useful in deciding an optimum schedule either for catalyst reactivation on-line or for disposal and reloading with fresh catalyst. The details of the experimental system are described elsewhere (Kershenbaum and Lopez-Isunza, 1986; Lopez-lsunza and Kershenbaum, 1992). Problem

statement

We consider the class of systems governed by the nonlinear differential and algebraic equations

[o j=LAu,u).~u),t)I+L;~u)j ,,u)] rf,(x,(t),xz(t),t)]

r~,(t)]

(1,

with initial condition x, (0) =Xl.o+ 5, (0)

(2)

where the n-vector state x(t) can be split into an n:vector differential state x~(t), whose time derivatives exist explicitly, and an n2-vector algebraic state x2(t). ~(t) is a zero mean random process with unknown statistical characteristics. By reference to xdt) and x2(t), ((t) can be decomposed into an n:vector (t(t) and an n2-vector (2(0. Observations of the systems are available either continuously in time or at discrete sampling instants and consist of an m-vector y(t) or y(tk). In the former case, the measurements are related to the state by

y(t) =h(xj(t),x2(t),t)+ rl(t)

(3)

and in the latter case by

y(tk)=h(x,(tk),X2(tk),tk)+r/(tk); k= 1,2,...,M

(4)

where r/(t) or r/(tk) is an m-vector of unknown measurement noise. Based on the measurements y(t) or y(tk) in the time interval 0-t--
1= ~ [ x j ( 0 ) - xl.0]rp ~ I(0)[Xl(0) --

Xl.O]

+~1 fr { [~'(t)]

rfl(x'(t)'X2(t)t)'l]r r. I ' ] ~ R - (t) Lf,.(x,(t),x:(t),t) J J

{ [0,,,,]_

},,,

Lfz(x,(t)~c,.(t),t) j

1 (r [y(t) -- h(xl(t),x2(t),t)] rQ -i(t ) / Jo [y(t) -- h(x I(t),x2(t),t)]dt

+~

Sequential estimation for nonlinear differential and algebraic systems where the first term minimises the square error of the initial estimate of xt(0), the second term minimises the integral square system modelling error and third term minimises the integral square measurement error. The weighting factors P~ J(0), R - ~(t) and Q - I(t) are symmetric and positive definite matrices which can be chosen to reflect the errors in the initial estimate of differential state, the process model and the measurement device. It should be noted that these weighting matrices have no strict statistical significance for the nonlinear case, but are analogous to the corresponding linear case with Gaussian white noise in which Pu(0) is the covariance of the initial differential state errors, R(t) is the covariance of the process noise and Q(t) is the covariance of the measurement errors. In the discrete time data case Q - I(t) is replaced by Qd- I(t) M

M

0,7 '(t)= kE=mQ -~(tk)6(t- tk)Ak= ~l Qk-16(t- t~) (6) where 6(t-tk) is Dirac delta function, and Ak is the sampling interval at sample k. The estimation problem is to determine x~(t) and x2(t) such that the objective defined by (5) is minimised.

In this section, a sequential estimator will be derived for the case in which measurements are available continuously in time. First of all, the estimation problem is reformulated as an optimal control problem. This approach has the advantage of not requiring statistical assumptions on the system error and measurement noise, and was first used in ordinary differential systems (Detchmendy and Sridhar, 1966) and later extended to partial differential systems (Seinfeld, 1969). Then, the calculus of variations for differential and algebraic systems (Jonckheere, 1988) is employed to obtain the sequential estimator. We define

u(t)=[~,(t)]_rfl(x,(t),x2(t),t)l Lf2(x,(t),x2(t),t) j

(7)

r

H= 1_[y(t) - h(xl(t),x2(t),t)] rQ - i(t)[y(t) - h(xl(t),x2(t),t)] 2 1 r t(t)u(t)+FAl(t!] r +~ u(t) R[-h2(t)J [fl(xl(t)'x2(t)'t)l+u(t)} [f2(x,(t),x,_(t),t) j

(9)

The condition aH/au=0 yields

j

(lO)

and the adjoint equations are

Q -'(t)[y(t) - h(x,(t),x2(t),t)] of,

[0xafz 0x l r (t) 1 of~ j [.a~(t)j ax,

(11)

ax2

with two boundary conditions A~(t)=0

(12)

xl(0) =Xl.o- Pu(0)At (0)

(13)

Upon substitution of (10) into (7) we have

[o't']=r"x'""')'l Lf2(x,(t),x2(t),t) j

(14) L,~2(t/j

(l l) to (i 4) form a two-point boundary value problem which can he solved for xl(t ), x2(t), At(t) and As(t), and thus generate the estimate. To be more explicit, we denote ,~j(ttT), .~2(tIT)and fi(tlT)

and rewrite the least squares error functional (5) as I

systems (Jonckheere, 1988), we can obtain the Hamiltonian

afl

Estimator for continuous time data

1053

as the estimates and control at time t, with measurement data y(t) up to time T. Thus

i

I = ~ [xl(0) - x,.0] P/] (0)Ix,(0) - x,.0]

.~,(tlT) and ~2(ttT) +~

u(t) rR - I(t)u(t)dt

1 r r [y(t) - h(x,(t),x:(t),t)] rQ -'(t) + 2 Jo [y(t) - h(xj(t),x2(t),t)]dt

(8)

Now the estimation problem can be considered as a deterministic optimal control problem, i.e., select the control u(t) such that (8) is minimised subject to the constraint described by (7). Applying the first order necessary condition for a week optimum in nonlinear differential and algebraic

are the estimates obtained from solving the two-point boundary value problem of (l l) to (141. According to (13) and using a variant of the formal method suggested by Ray (1981), we make a transformation

Yq(tlT)=wt(t) - P~ i(t)Al(t)

(15)

and use the following relation for .~,_(tlT)

r.,,)]

~2(tW)=wa(t) - [P~,(t) PMt)I [,~:(t)]

(16)

where w(t) is an n-vector and P(t) is an n X n-matrix

Y. S. CH~NO

1054

which are to be determined. Based on &(t) and x2(t), w(t) and P(t) can be split into ..

P(t) =

r d

[P2,(t) P2~(t).i

(18)

et al.

r~(.~i(flT), .t2(tlT), t 1 L~2(2t(tlT), ~2(tll,'), t)

=[~,(wl(t),w2(t),t)]+[Ox, OXzoh ] L~:(w,(t), w2(t), t) 0h OXl Ox2

{[~,(,m]_[w,(t)]~

h(Ycl(tlT)'2"(tlT)'t)'~h(wl(3x---~l t)'woh]~2(tx=3h)'t)+[-

By combining (15) with (16), we have

[,,,,q

(22)

Lw~(t)J J

L~2(~r)J

Note that P(t) is a symmetric matrix, i.e., Pi2(t)=P21 r (t).

r-..<..1

<,(tlr)j=Lw~(t)jLP,,(t) G2(t) La~(t)J (19)

By differentiating (15) with respect to time t and rewriting, we obtain

[ '¢Col(tlT)]=[~l(t)]-[P2]F*"t' oU(t]-L~2(t) ) [:,,<,> [:,<,,]

{ [~,(t,~]_ [w,(t)]~, [,#r)j [w~(t)jj

(23)

By substituting (22) and (23) into (21) and using (19), then collecting terms on

Ai(t)

J

&(t) j'

(20)

we obtain

[o,(t)]=r~I(wi(t),w2(t),t)] L~2(w,(t),w2(t),t) [

Upon substituting (11) and (14) into (20), we get

[~2(Yq(tIT)'Yc2(tRrprAi IT)'t)I(t)l+i([~tlT)] _[:,,<,,:}r,,
+

Ox~l

P2z(t)] h(wl(t), w2(t), t)]

Q-'(t)[y(t)-

(.~,(,1/3,&(tit), t)J : "L&(t)J

0Z [ ~ ii(t)

J

0~,

: ] _ r 0xi 0 x : l [ , ( t > 0 - L 0.Jb2 372 ..i LP2i(t)

O&

-

]

P2z(t)

OX I

r........~.,. 1 [ ox +[0 P22(t)J o7=

(24)

OX2

0.2 l'r a,(t)l

o~,

o?,

072 _] La2(t)]

+Lo

Ox2

PMt)J Ox~

[y(t) - h(.~l(tlT), ,~2(tlT), t)}

!

(kl(tlT), &,(tiT), t)] (Yfi(tlT), 22(tiT), t)j

and h(;c,(tlT), 22(tiT), t) about the point w,(t) w,(t) J

Ox~ O&

[o,,<,,,,-
PMt)J L 7x, 0x~

(2 )

By applying Taylor series expansions truncated after the first order term, we linearise the vector function

[~

-

,'2:,,uoz ~,~

3x~i~ J LP2,(t) Pzz(t)

(25)

with the initial conditions WI(0)=XL0

(26)

Pu(0)=Pit.o

(27)

Note that the 2(TIT), the minimal least square estimate at time T conditional on all the measurement data up to time T, is always equal to w(T) (Ray, 1981). Thus the

Sequential estimation for nonlinear differential and algebraic systems estimates are determined from (24) where t is always the current time.

= L.f2(2,(tlt), 21(tl t),t)j

P22(t)J

Ox--~l~

each time step. However, this will increase significantly the computing time which is crucial for online applications.

Estimator for discrete time data

[P,,(t) P,2(t)l[ a~ta~tlr +

1055

In deriving the discrete data estimator, we employ the same procedure as above, except that we use the t~ function in (6) to incorporate data at discrete time points.

j

Q -I(t)[y(t) - h(2t(tlt), 22(tit), t)]

(28)

Prediction between measurements subject to initial condition 21(0)=Xl.o

(29)

The resulting estimator, which consists of (28) and (25) with initial conditions (29) and (27), is in the form of differential and algebraic equations and can be solved by standard differential and algebraic equations solvers such as DASOLV (Jarvis, 1994). Based on the derivation of the estimator, we observe the following: 1. (28) and (25) can be applied to any differential and algebraic system in the semi-explicit form described by (1) without considering the index of the system. However, for a system with high index (greater than one), the corresponding estimation equation is also a high index system. At the present time, generalpurpose algorithms for the numerical solution of differential and algebraic systems are restricted to index one systems; therefore, special techniques will be required in order to apply the estimation algorithm directly to higher index systems. 2. When applied to ordinary differential systems, (28) and (25) simplify to the standard EKF which can be expressed in terms of ~ and PH. 3. In the same way as in ordinary differential systems, the matrix P(t) in linear differential and algebraic systems corresponds to the covariance of the estimated error. For nonlinear differential and algebraic systems, although the matrix P(t) can be loosely associated to the uncertainty of the estimated state, it has no clear statistical meaning. 4. No general method is available for selecting the P.(0), Q(t) and R(t). The use of first order Taylor approximations for deriving the estimator would imply that these weighting matrices are linearised approximations to the covariance matrices and that, if these are available, they would provide a good working choice for these weighting matrices. The practical guideline of selecting these weighting matrices will be discussed later in the case study of a fixed-bed reactor. 5. The first order linearisation used in the derivation of the estimator is the main difference between the sequential estimator approach and the NLP approach. The robustness of the sequential estimator can be further improved by using repeated linearisations for

Between two successive data points, t,_ ~< t< tk, there are no measurement data. By noting that Q2~=0 and from (28) and (25) the following prediction equations can be obtained

(30)

o J - [?2(~,(t),~:(t),

0h 0h ] [P.(t) o LP2,(t)

Oxj

+

P22(t)]

P22(t)]

Ox2

0h

0~2

(31)

oqXl The initial conditions are the same as (29) and (27).

Updating at sampling points Because measurement data arrive at discrete points in time tk, .~l(t), 22(0, Pu(t), Pi2(t) and P22(t) must be updated at this sampling point. With the use of the procedure developed by Ajinkya et al. (1975), these updating equations can be obtained by carefully considering how the 8 function in (6) enters the continuous estimation equations derived in the last section. By integrating (28) from t, - to t~" across the sampling point, we get

oh oh ]"

ox--~,ox--~ Q['[y(t,)- h(~,(t,), ~2(t~), tN

(32)

In order to obtain a practical updating expression, by representing the second term on the right hand side of (32) with a first order Taylor series expansion about t, - , we get the updating equations for .t~ and Yc2

! 056

Y. S. CHENGet al.

']

2

I,

,

OX

r.,,.,

+[

P~2(t;)jt

o

o, ox, TAx,. Q;

Formulation of a sequential estimator for a fixedbed reactor

PMtf)J t ~2(ti)J

(33)

As in the derivation of (33), from (25) we can obtain the updating equations for PH, P~2 and P22:

[,,,:) :]__[,,:,00] -

[,.,,;:) ev_(t;)j ...2.:,1rob o.]T L~ ,ro

Q; L ox, ~

LP2,(t;)

= PA

~" COM

_ [Pll(otk-) P,2(t/)]

Ok- [ ~Xl ~x2

5

) PT

Fig. I. The complex kinetic scheme.

[y(tk) -- h(.~l(t k-), 22(t k-), tk)]

--Ox, ~x2 ]

) OT

o] PE2(tk-)

_[p,,;r,)P,4t[)] £(,:)j

In this section the reaction network, catalyst deactivation scheme and reactor model will be briefly outlined and a sequential estimator for the reactor is formulated with the use of the general algorithm for differential and algebraic systems.

Reaction network Although the production of phthalic anhydride by the partial oxidation of o-xylene has been studied for over three decades, there is no general agreement on the reaction network. Among many reaction networks proposed so far, one suggested by Caiderbank (Calderbank et al., 1977) has been widely used and is given in Fig. 1 where OX indicates o-xylene; OT o-tolualdehyde; PT phthalide; PA phthalic anhydride and COM combustion products (CO+ C02). It is well known that the computing time for solving estimation equations is crucial for the on-line application of the estimator. In order to achieve rapid solution, the reaction network shown in Fig. 1 is lumped into a simpler one shown in Fig. 2 by assuming infinitely fast rates in step 4 and step 5 and thus, involving only a single component (o-xylene). At 1 bar total pressure, the reaction rates (kmoll(kg.s)) can be described by

R~.'(O,c)=d~k.c; s= 1, 2, 3

0]

SP.(tD

8P.,i(tk)8P2z(tk)

(34)

(37)

where the subscript s represents the sth step reaction. The reaction rate constant k,,, and the Redox parameter 4) are expressed as

k.=k.,oexp( - E J 0 ) ; s= 1, 2, 3

(38)

~b=

(39)

where

~t:(tD j = L~2(t;) - . h ( t f ) J 8P,,(t,)

(35)

8e,,(tk)]

8P2,(tD

8P£(tDJ

[P,,(t;)-P,,(t;)

P,2(t~[)-P,:(t~')]

=[P2,(t~)-P2.(t[) P22(tD -

P~2(tf)J

(36)

(33) and (34) are quite complicated. However, in particular applications, these equations will become considerably simplified as will be illustrated in the case study of a fixed-bed reactor which follows. Simulation studies have also been carried out for a CSTR with discrete time measurement data (Cheng and Kershenbaum, 1995a).

k,,,p,,: k ..p o:+ ( ~lk,,, + v2k,,..+ v3k,.)c

where ~ is the number of oxygen molecules involved in step s, /9 indicates dimensionless temperature which is defined relative to the nominal feed temperature of 449 K and c denotes dimensionless o-xylene concentration which is defined relative to the nominal feed composition of 1 mole % o-xylene. Here, K.oPo.- can be regarded as approximately constant because of the large excess of 2

OX

~ PA

= COM Fig. 2. The simplifiedkineticscheme.

Sequential estimation for nonlinear differential and algebraic systems

1057

Table 1. The kinetic parameters values. vt = 3.0 v: =3.0 va=6.5

k,,= 0.0383 kmoll(kg.s) k,~,= 0.0130 lanoll(kg.s) k,, =0.0036 kmoll(kg.s) k,,,,Po:"--'7 22 × I 0 - ~'kmoll(kg.s)

Reactor dynamic model

E,, = 16.41 E,:= 14.59 E,, = 13.70

oxygen in the feed. The values of the kinetic parameters are given in Table 1. Although the simplification of the kinetics would introduce error in the reactor modelling, estimation algorithms, which are able to handle uncertainty, could effectively compensate for the inaccurate model.

Catalyst deactivation scheme The deactivation of the V2Ofl'iO2 catalyst can be irreversible and reversible. The main cause of the irreversible deactivation is the phase transformation of TiO., from anatase to rutile. However, some researchers (Haber et al., 1986; Nikolov et al., 1987) found that unless the catalyst was subjected to high temperature for a long period, this form of irreversible deactivation could be ignored. Two possible causes of the reversible deactivation are the reduction in the oxidation states of vanadium and the deposition of some strongly adsorbed species on the catalyst surface. In this study we only consider the reversible deactivation caused by some adsorbed compounds which are formed from the nonselective oxidation of o-xylene. These adsorbed compounds deposit more rapidly at high hydrocarbon concentration and can be partially removed by increasing oxygen partial pressure. This simplified deactivation-reactivation scheme has been justified by microreactor experiments (Mongkhonsi, 1993; Mongkhonsi et al., 1994). The following deactivation-reactivation model has been proposed to describe the variation of catalyst activity a with time.

Many models have been proposed to describe the dynamic and steady-state behaviour of fixed-bed catalytic reactors. Selecting a model must be made based on the requirements of a specific application. In this study our objective is to develop a sequential estimator which could be used for the purpose of control and on-line optimisation. An acceptable computational time for solving estimation equations is a very important factor. Therefore, we should choose as simple a model as possible. Extensive dynamic experiments (Mongkhonsi, 1993) demonstrated that a two-dimensional pseudohomogeneous model without axial dispersion could satisfactorily predict the dynamic behaviour of the jacket-cooled fixed-bed reactor used in this work. The reactor dynamic model consists of the following set of partial differential equations describing the mass and energy balances in dimensionless form.

O0

-:-+ Ot

80 Oz

)

1 //020

~

OC 0=-~Z+~.

/

1 00~

3

~r 2 + r Or +a E=,B~R,(O,c) (41)

1 / 02C

1 OC \

3

~ 0 r - - ~ + - r Or ) + a ~, = Rs(O,c) (42)

where Rs is the dimensionless rate of the sth reaction, z represents the dimensionless axial co-ordinate, r is the dimensionless radial co-ordinate, and t is the dimensionless time, t=-t'/r, where r is the relative thermal time constant. The catalyst deactivation dynamics described by (40) can also be written in terms of dimensionless time. The initial condition is 0(r,z,0) = Oo(r,z)

(43)

Boundary conditions are at r=O and 0 - z - I 00 0r

dt--7 = - kd,.oeXp -- -~- ca

--

0c

--

+kd:flo2exp( - -~ )(am-- a)

Or

(40)

=0

(44)

=0

(45)

at r= I and 0--
-'

--

Or

=0

(46) (47)

and at z=0 and 0-
(48) (49)

where subscripts in and w represent reactor inlet and wall respectively. The values for the model parameters are given in Table 3. The time scale for the concentration Table 3. T h e reactor parameters values.

R~uRI'I1.25 Ee, = 8.02

× 10 - 5 Ph,=8.98 × 10 - : P,r=7.81 × l0 - :

Ee:=66.82

Bi,=0.80

Table 2. T h e deactivation parameters values.

kd = 0.226 S - ~ kd":=2 4 × 10 *'~(kPa.s)

00

- - = - Biw(O- 0~) Or Oc

BI = 84.19

82=84.19 B3= 156.57 I'=-2.77 × l03 s

1058

Y.S.

CHENG

et al. k __ k + l a~c.2a,

transient is much smaller than those for catalyst deactivation and temperature transients; thus concentration can be assumed to be at quasi-steady state.

ok

The formulation of a sequential estimator for the reactor follows the strategy of early lumping, i.e., discretizing the distributed parameter system happens before the implementation of the lumped parameter estimation algorithm. Alternatively, late lumping can be used, i.e., applying distributed parameter estimation theory to obtain the estimation equations in the form of partial differential equations before numerical solution by a lumping approximation technique. The comparison of early lumping with late lumping can be found elsewhere (Cheng et al., 1996b). Because the temperature profile in the axial direction in the reactor is very steep, and the location of the hot spot may vary significantly with small disturbances in operating conditions, orthogonal collocation on finite elements is used for discretisation in the axial direction. It has been shown, however, that there is no computational advantage in using moving rather than fixed finite elements (Mongkhonsi, 1993). For the radial direction, due to the symmetry and smoothness of the profile, only one collocation point is used. The reactor model described by (40) to (42) is lumped into following differential and algebraic system.

- ~i / c'°ti

e,,.) (a,,.- af)

+k,,:.eo:~Xp - ~ aO~ I d~- = - ~

(50)

NC+2 8Bi.. O~ J~=' AO0~+ ehr(~-Bi..) (0..- ,) 3

B,R,(Oi,ci) 1

(51)

NC+2

3

k +ot ki .,=1 Z 0= - - -AZ k j =~,l A.~: u ~

k k Rs(Oi,Ci)

(52)

where superscript k indicates the kth finite element with the length of Azk; subscript i is the ith point of NC interior collocation points and A represents the approximation matrix for the axial spatial derivative. Initial conditions are

a~(0) a~o =

k

_

k

~c÷2- vl

(58)

k

k÷l

(59)

C NC+ 2 -_- C I

Sequential estimator for the reactor

d°t--~=-kd"°rexPdt

(57)

~k+ I

_

(53)

Oi(O)- 0j.o

(54)

0[(t) = 0~,,(t)

(55)

c l(t) = G,(t)

(56)

Inlet conditions are

In addition, the continuity of activity, temperature and concentration between the boundary of finite elements is required.

In defining the vector notation -- r

I

I

2

2

--

I

I

2

2

NE

NE

1T

a ' - - t Ot r " tXNc+ 20~ 2"- ot Nc+ 2""t1'2 .-.otNc+21 NE

NE

T

0-[02""0~c÷20".'"0uc+2""0,. '"0uc.2] _ r cl

C--

L

2""C

I 2 2 uE NE , r N c + 2C 2 " " C N c + 2 " " C 2 " " C N c + 2 j

(60)

(61) (62)

where NE is the number of finite elements, and considering modelling errors, (50) to (59) can be put in the vector form, dot

dt =f=(a,O,c)+ (.

(63)

~ =fn(a,O,c)+ (o

(64/

O=f~(a,O,c)+L

(65)

where nonlinear functions f,, fo and f~ correspond to catalyst deactivation, energy balance and material balance respectively. (,, (o and L represent corresponding model errors. The initial conditions are a(0) = ao

(66)

0(0) = 0o

(67)

In order to obtain the final estimation equations for the reactor, the following procedure was used. !. Identify the differential and algebraic variables. From (63) to (65) it is obvious that the differential variables are catalyst activity and temperature and that the algebraic variable is concentration. Hence, we use the following decomposition

x=[& ~x2]r=[a O'c] r

(68)

f=[f~ if,.] r =[f~fo~f.] r

(69)

2. Split the P matrix and R matrix based on differential and algebraic variables. With reference to (68) and (69) we can obtain the following partition

P..

P~o

~ Po~

P,.,

P,.0 ~

P,.,.

A similar partition has been done for R. 3. Choose a measurement structure. Because catalyst activity can not be directly measured and time delay exists in obtaining concentration measurements from a gas chromatograph, we only consider temperature measurements at discrete axial positions,

y(tk) = h( O(t,)) + rl(tk)

(71 )

where h is the function which interpolates temperature at any axial position from the values of temperature at the collocation points. "q is measurement noise. With the use of (68) to (71) and on the substitution of

1059

Sequential estimation for nonlinear differential and algebraic systems (63) to (65) into the discrete time data estimator developed earlier, we obtain the reactor estimator. Between two success sampling points, tk-,
(72)

d___~=~(&,~,~.)

(73)

dt

dt

~(t~) = &(t~-)

+Po~(t~) ~

an _n(~t~))l_P~o(t;)[.~]r -

0 =f~(&,~7,~)

(74)

dP*"= aJt"p,,.+~o[P~#lr+~c [P~lr dt Oa

Qi"[y(t,)

Qf ,On ~[a(t,) -- +

~(t[)]

(86)

~t~)= ~(t[) b1 T +p~,,)[F o~] Q: , L~(,,) +

oh ] , oh ^ + -~ ] Q; , ~[00~)

+R,,.

(75) -

~(t[)]

(87)

dP'° = Of~ Oa P~°+ ~o P°°+ ~ on 1 T

+R,o dt

Oct

~n

(76)

O0

- P.8(t[)] r P.o(tk)-P.~(tk )

Oc

+

(88)

_

ah 1 T +Roo 0=

T

oh

(77)

[°~] ~b T '°n

Of. p~ Oc

-P~o(t[)

O; ~o[Poe,(t~[)

- Pogtf)] PogtD=Po4tf ) +R~.

(78)

OLp,.~ 0--~c +te=l~[ oZ ~

oo j

[Oh] T ° n -Poo(tD 0-0 Q;' ~Poo(tf)

L oc

+R~

Q? O0 [Pdt~) (79)

(8O) with the initial conditions

~0)=&o ~(o) = ~,, P~a(0) = P~a.o

(81) (82)

Pa~(O)= Po..o

(83) (84)

Pod0) = Poo.o

(85)

At the sampling point tk, the updating equations are

(89)

-

Pso(t[)]

(90)

Implementation of the reactor estimator In carrying out the implementation of the sequential reactor estimator described in the previous section, we start by integrating (72) to (80) from time t0 with the initial conditions (81) to (85) to the first sampling instant t~ to obtain the prediction values, then calculate the updating values from (86) to (90) with the use of the measurement y(fi), and repeat these two steps within any succeeding time interval (t,- ~, t~). In this section, we will describe the numerical method and how to choose the values for the parameters in the estimator.

1060

Y. S. CHENG

Solution method The number of finite elements and interior collocation points in each element will be determined by a compromise between accuracy and computational load. If we use 5 finite elements and 2 interior collocation points within each element, we will have 16 catalyst deactivation equations; 15 temperature equations and 15 concentration equations (the inlet conditions for both temperature and concentration can be eliminated by simple algebraic manipulation). In addition, we also need to solve 256 equations of pa~; 240 equations of p,O; 240 equations of P~; 225 equations of pOo; 225 equations of P e~.and 225 equations of P". Because P~, p~o, p~,., poo and P " are nonlinear functions of the estimated state variables &,~ and ~, we must solve (72)-(80) simultaneously. Therefore, the number of equations which need to be solved simultaneously is 1141 which is quite large and the solution process can be very time-consuming. In order to reduce the computational burden for the real-time implementation of nonlinear estimation algorithms, some researchers (Windes et al., 1989) solved the P matrix equations in advance off-line using a set of nominal values of the state variables. In this study, because of significant changes in the catalyst activity and the large shift in the position of the temperature profile, it is not adequate to use one set of nominal P for the estimation of activity, temperature and concentration. Fortunately, the structure of (72)-(80) is quite sparse and can be solved in realtime by using a differential and algebraic system solver (Jarvis, 1994). In our work, five elements and two interior collocation points within each element were found to provide satisfactory accuracy in all situations. In more time-critical situations, it may be possible to update P less frequently than at every time interval.

Choice of estimator parameters When solving the reactor estimation problem, several parameters must be selected in advance. These are: 1. 2. 3. 4.

number and position of temperature measurements. initial values for P matrix. measurement error weighting matrix Qkmodel noise weighting matrix R.

To investigate the importance of selecting measurement set, we considered different numbers of temperature measurements at different axial locations. In a realistic industrial situation, this would be fixed by the physical location of a few temperature sensors. We were in the fortunate situation of having a large number of possible measurements available, so we could investigate effectively the importance of this parameter. The magnitude of the P matrix is an indicator of the estimation quality (Kumar and Seinfeld, 1978). Thus, for the given operating conditions, if desired, alternative measurements can be compared for a set of estimator parameters selected previously, and a lower overall P matrix value indicates a better measurement set. The matrix P initial value reflects the magnitude of the errors in initial conditions, and thus should be large

et al.

when little is known about the initial state and small when better initial conditions are available. Generally, the initial value of P is not critical. The Qk weights the measuring device errors. The value of Qk should be large, when measurement errors are large, and small in cases of small measurement errors. In this study, it was assumed that all temperature measurements were independent and relatively accurate; therefore, Qk was taken as a diagonal matrix with all elements having the value of ~ 10 -~0 (dimensionless). The selection of the R value reflects the degree of confidence in the model. Usually, the R value greatly influences the performance of the estimator, and thus we must choose with caution. In this paper, we selected values of matrix R as to be symmetric and decreasing with increasing distance between zi and zj R~(ij)=3.0 × 10-14(1 - ( z i - z j ) 2)

(91)

R~(ij)=6.7 × 10-'5(1 - ( z i - zj) 2)

(92)

Ro~(iJ) =4.0 × 10 -

13(1 - (zi - zj) 2)

(93)

R~.(ij) = 1.0 × 10 - ~4(1 - (zi - zj) 2)

(94)

--(Zi--Zj) 2)

(95)

Rcc(id) = 1.0 × 10 - 1o(1 - (z~- Zj) z)

(96)

R~(id')=2.2 × 10-t3(l

Estimator performance The performance of the reactor estimator developed in last two sections has been examined using real-time dynamic experimental data. Since at this stage of the work no signals were fed back to the reactor, the reactor estimator was tested through off-line implementation. Except for possible constraints of computing time, the results are identical for on-line estimation, using software which we have developed to collect and process the data in real-time. Only temperatures within the reactor are used as measurements. Twenty-six temperatures at different axial location of the reactor are measured at 12 s intervals. However, only two or three of those are utilised in the estimator; the remainder are used to represent the actual temperature profile in the reactor. The effectiveness of the estimator can be demonstrated by comparing the actual temperature profile to the estimated one generated by the estimator. In order to investigate the performance of the estimator, the following dynamic experiments have been carried out I. reactor start-up. 2. ramp changes in the salt bath (coolant) temperature

o.. 3. step changes in the o-xylene feed composition ci,,. Due to limitation of space, only the experimental and estimated results of a dynamic change of the coolant temperature will be given in detail. The other results are reported elsewhere (Cheng et ai., 1995b). This experiment involved beginning at steady-state operation and performing a ramp change in O,. Before the experiment the catalyst bed had been reactivated at

Sequential estimation for nonlinear differential and algebraic systems an air flow rate of 4 Nm3/hr and O.,=378°C for 19 h. The reactor was started with ci.=0.5, then. ci. was increased to 0.7 after 72 min and to 1.0 after i10 min (data not shown). During the start-up of the reactor, (9.. was kept constant at 383°C. After operating the reactor under the above conditions at steady-state for 30 min, a ramp increase in O.. from 383 to 403°C at a rate of 1°C/rain was made. The reactor was operated at Ow =403°C for 70 min, then (9., was decreased from 403 to 383°C at a rate of 0.22°C/min. The disturbance pattern for (9., is plotted in Fig. 3. Three temperature measurements located at 0.4 m, 0.8 m and 1.0 m respectively from the entrance of the reactor were used by the estimator. Other number of measurement points were considered; these represented the best compromise between accuracy and computational speed (Cheng et al., 1995b). The estimator assumed steadystate profiles as initial values. Due to the errors of both the reaction and deactivation models, the initial guess of the temperature profile was quite different from the experimental data. In addition to testing the dynamic response of the estimator to a ramp change in Ow this experiment also examines the efficiency of the estimator in the case of incorrect initial conditions. The estimated and experimental temperature profiles are shown in Fig. 4 for the period of the ramp increase (9.. and in Fig. 5 for the period of the ramp decrease in O~. Although the initial temperature profile was wrong (t=0), the estimated temperature profile converged to the true data very quickly (t=30 rain). The estimated steadystate temperature profile accurately predicted the presence of the shoulder in the region where the reactor temperature was almost equal to (9.. During the ramp increase in (9.. for the first several minutes, the hot spot gradually flattened and harmonised with the rest of the

1061

temperature. Then, the hot spot increased very rapidly, developed a peak again and moved upstream. All of these dynamic features were excellently tracked by the estimated temperature profiles shown in Fig. 4. For the ramp decrease in O., the hot spot moved very smoothly compared with that in the ramp increase in O.. The hot spot decreased and moved downstream. At the end of the experiment, the magnitude and location of the hot spot were different from those before the ramp change in O.. The estimator was able to characterise the shifting temperature profiles. The estimation of activity profiles during the ramp change in (9., is given in Figs 6 and 7. The difference between the initial guess of the activity profile (t=0) and the steady-state activity profile before the ramp increase in Ow (t=30 rain) explains why the initial guess of the temperature profile was different from the true steadystate temperature profile. It was noticed that, due to the deactivation mechanism there is an abrupt increase of catalyst activity in the region of the reactor entrance. This was reason that the temperature profile exhibited a shoulder. During the ramp increase in O., the activity of the catalyst in the front of the hot spot was increased due to the reactivating and the activity profile moved upstream. This caused the upstream movement of the hot spot. For the ramp decrease in O~, the activity profile moved downstream very slowly compared with that in the ramp increase of (9.. At the end of the experiment, a new steady-state activity profile, which is different from that before the ramp change in (9., was initiated. The behaviour could be attributed to the irreversible and reversible deactivation (Lopez-Isunza and Kershenbaum, 1992) and resulted in the different steady-state temperature profiles. The estimated concentration profiles are shown in Figs 8 and 9. The concentration profile

410

405

400

395 UJ n"

W n

I-

385

380

375

370

30

Fig. 3. Ramp changes in saltbath temperature.

i

i

I

I

i

90

120 TIME (MIN)

150

180

210

240

1062

Y.S. CHENG et al. 600

i

i

data estimate time(min) 1 0 o 2 30 x 3 40 + 4 50

550 500

*

0o

5

60

450

v

:~ 400

~- 300

200 I

15°o

,

I

,:s

2

2s

REACTOR LENGTH (M)

Fig. 4. Comparison of estimated temperature with experimental temperature during ramp increase of salt bath temperature. moved upstream or downstream with the ramp increase

or decrease in O.. Two different steady-state concentration profiles were obtained before and after the ramp change in (9.. These phenomena were consistent with the variation in temperature and activity profiles. The usefulness of the reactor estimator can be demonstrated by comparing experimental data with the results obtained by the estimator on the one hand and by open-loop model predictions on the other hand. Figure l0 shows the estimated temperature profiles and the 600

i

i

predicted/simulated temperature profiles obtained by using the identical reactor model both with and without catalyst deactivation. In both simulations, the initial activity and temperature profiles are assumed to be known. From Fig. 10 it can be seen that when the catalyst deactivation is included into the reactor model, the open-loop model prediction is substantially enhanced; however, it is still unable to track the dynamics of the reactor. By contrast, the estimator has successfully tracked its intricate dynamic behaviour. i

i

data estimate time(mln) o 1 120

550 500 1

A450 Ul

rr 400 uJ 350 UJ

I- 300 250 200 15(

0

i

I

I

L

0 5

I

1.5

2

2.5

REACTOR LENGTH (M)

Fig. 5. Comparison of estimated temperature with experimental temperature during ramp decrease of salt bath temperature.

Sequential estimation for nonlinear differential and algebraic systems

1063

estimate time(min) 1 0 2 30 3 40 4 50 5 60

0.9

0.8 0.7 0.6 ~0.5 O<0.4 0.3 0.2 0.1 t

0

015

1

I

115

2

2.5

REACTOR LENGTH(M) Fig. 6. Estimationof catalyst activityprofileduringramp increaseof salt bath temperature. Conclusions A sequential estimation algorithm has been derived for nonlinear differential and algebraic systems. By choosing the minimum least squares criterion, the estimation problem is converted to an optimal control problem without any assumption of the statistical nature of system modelling and measurement errors. With the use of the calculus of variations for differential and algebraic systems, we obtain a two-point boundary value problem. The linearisation approximation is made to obtain the sequential estimator. Both continuous and discrete time data estimators are derived. This nonlinear estimation algorithm can be used for general nonlinear differential and algebraic systems without consideration of the index. When applied to ordinary differential 1

systems the estimator simplifies to the standard form of EKF. The sequential estimator has been applied to a nearindustrial scale fixed-bed reactor with decaying catalyst to simultaneously estimate dynamic changes in catalyst activity, temperature and concentration profiles. Considering the computational limit for the on-line application of the estimator, the reaction network was lumped into a simpler one which required only the material balance of the main reactant (o-xylene) and a simple, experimentally verified catalyst deactivation scheme was used. Orthogonal collocation on finite elements was used to discretise the distributed parameter reactor model into a differential and algebraic system for which the sequential estimation algorithm could be implemented. The final estimator for the reactor was

•

estimate time(min)

1 2 3

0.9

0.8

120 150 180 215

4

0.7 O.E

~-o.G 0.4 0.3

0.1

%

t

I

t

2

REACTOR LENGTH (M)

Fig. 7. Estimationof catalyst activityprofile duringramp decrease of salt bath temperature.

2.5

i 064

Y.S. CHEN~ et al. 1[~-----~~ x 0.91~ ..~ I ~, ~o8L ~ z '1 \\ o I \\

estimate time(rain) 1 0 2 30 3 40 4 50 5 60

\\ \ \\ \ \ \ \ \ \ \ \ \ \

\\\

\\\

~ o.41-

\\

zo.21

4~5

\ \

\\

\ \ 1~,'

.~

~0.1 0 0

0.5

~

1 1.5 REACTOR LENGTH (M)

2

2.5

Fig. 8. Estimation of concentrationprofile daring ramp increaseof salt bath temperature.

quite sparse and could be solved in real-time by using a differential and algebraic system solver. Using two or three properly located temperature measurements, excellent estimation has been obtained for a range of situations, such as reactor start-up, ramp changes in coolant temperature and step changes in oxylene feed concentration. The estimator rapidly converged when given a poor initial guess. The complex dynamic behaviour of the reactor (including wrong-way behaviour and non-unique steady-state) has been successfully tracked. The catalyst activity profile was accurately inferred. The estimated variation of the catalyst activity enable us to fully understand the dynamic changes of the hot spot magnitude and the movement of the temperature profile. The sequential reactor estimator can be easily applied on-line in an analogous way. The computational feasi-

bility of applying the estimator in real-time is demonstrated by the fact that the required computing time is one-seventh of real operation time on our SUN SPARC 10 workstation. Nomenclature

A=approximationmatrix for axial spatial derivatives B,--dimensionlessheat of reaction for reaction s Biw=Biotnumber c=dimensionless concentration of o-xylene (relative to nominal feed) E,, =dimensionless activation energy of sth reaction step Ed,,Ed=dimensionless activation energy for deactivation/ reactivation f=nonlinear function H=Hamiltonian h =measurement function

I estimate time(min) 1 120 2 150 3 180 4 215

0.9 U-

0.8 O

0.7

~ 0.6 O

0.5

O

0.4 ~ (1.3 ~0.2 ~o.t 0

0.5

1 1.5 REACTOR LENGTH (M)

2

Fig. 9. Estimation of concentrationprofile during ramp decreaseof salt bath temperature.

2.5

Sequential estimation for nonlinear differential and algebraic systems

1065

600 A

to

*---500 I11 rr <~ 400 nr

/ ' A

120 min

Ilc

~ 300 I- 200

I- 200 0.5

1

1.5

2

2.5

0

REACTOR LENGTH (M)

0.5

1

1.5

2

2.5

REACTOR LENGTH (M)

600

60O A

*---500 40o I1: ~. 300

/ P/~ ~ d ~ ,

to

180 min

*---500

,, -/,~ "-2~'~ I

S ....

/ ./

~400

kU

215 min

/

i..- 2 ~

I.- 200 0

0:5 i 1:5 ;~ REACTOR LENGTH (MI

2.5

0

0.5

1

1.5

2

2.5

REACTOR LENGTH (M)

Fig. 10. Comparison of experimental data with estimation and simulation during ramp changes of salt bath temperature. * data, estimate, - - - - - - simulated without catalyst deactivation, simulated with catalyst deactivation.

- -

/=least squares functional

feed temperature) u,=number of oxygen molecules involved in sth reaction step r=thermal time constant, s 4~=redox parameter ,~=adjoint variable vector

k,,=sth step reaction rate constant, kmol/(kg, s) k,,~ =pre-exponential factor of sth step reaction rate constant, kmol/(kg • s) ka,=deactivation rate constant, s kd,, =pre-exponential factor of deactivation rate constant, S-I

k,/2=reactivation rate constant, (kPa-s) k,/:o=pre-exponential factor for reactivation rate constant, (kPa-s) k,., po =constant for adsorption of oxygen, kmol/(kg, s) M=sampling number NC=number of interior collocation points NE=number of finite elements P=state estimation error matrix Ph,=Peclet number for heat transfer P,,r=Peclet number for mass transfer Po:=partial pressure of oxygen, kPa Q=measurement error weighting matrix R=model error weighting matrix R',=reaction rate of step s, kmol/(kg.s) R,=dimensionless reaction rate of step s r=dimensionless radial co-ordinate T=time t=dimcnsionless time t' =time, s u =control variable vector w=transformation variable vector x=state variable vector y=measurement vector z=dimensionless axial co-ordinate

Greek letters a=catalyst activity am=maximum available activity A =sampling interval ,dzk=length of kth element 8=delta function g'=model error 'r/=measurement error O---dimensionless temperature (relative to nominal

Subscripts d=discrete time i=index of interior collocation point in =inlet condition j=index of interior collocation point k=sampling index s=reaction step index w=reactor wall 0=initial condition 1=differential variables or equations 2=algebraic variables or equations

Superscripts k=index of finite element +=immediate after sampling instant - =immediate after sampling instant ^ =estimate •=time derivative

References

Ajinka, M. B., Ray, W. H. and Froment, G. E (1974) Online estimation o f catalyst activity profiles in packed-bed reactors having catalyst decay. Ind. Engng Chem. Process Des. Dev. 13, 107-112. Ajinkya, M. B., Ray, W. H., Yu, T. K. and Seinfeld, J. H. (1975) The application o f an approximate nonlinear filter to systems governed by coupled ordinary and partial differential equations. Int. J. System Sci. 6, 313-332. Calderbank, P. H., Chandrasekharan, K. and Fumagalli,

1066

Y. S. CHENGet al.

C. (1977) The prediction of the performance of packed-bed catalytic reactors in the air-oxidation of o-xylene. Chem. Engng. Sci. 32, 1435-1443. Cheng, Y. S., Abi, C. E and Kershenbaum, L. S. (1996a) On-line estimation for a fixed-bed reactor with catalyst deactivation using nonlinear programming techniques. Computers chem. Engng. 20, $793$798. Cheng, Y. S. and Kershenbaum, L. S. (1995a) Sequential estimation for nonlinear differential and algebraic systems--I. Theoretical development, Technical Report B95.29, Centre for Process Systems Engineering, Imperial College. Cheng, Y. S., Mongkhonsi, T. and Kershenbaum, L. S. (1993) Estimation of catalyst activity profiles in fixed-bed reactor with decaying catalysts. Appl. Catal. 106, 193-199. Cheng, Y. S., Mongkhonsi, T. and Kershenbaum, L. S. (1995b) Sequential estimation for nonlinear differential and algebraic systems--II, Application to fixed-bed reactors with decaying catalysts. Technical Report B95.30, Centre for Process Systems Engineering, Imperial College. Cheng, Y. S., Mongkhonsi, T. and Kershenbaum, L. S. (1996) Nonlinear dynamic estimation for a fixedbed reactor with decaying catalysts. Chem. Engng Sci. 51, 1909-1918. Chisci, L. and Zappa, G. (1992) Square-root Kalman filtering of descriptor systems. Syst. Control Lett. 19, 325-334. Clement, K., Jorgensen, S. B. and Sorensen, J. E (1980) Fixed-bed reactor Kalman filtering-experimental investigation of discrete time case without stochastic disturbances. Chem. Engng Sci. 35, 123 l-1236. Darouach, M. and Zasadzinski, M. (1992) State estimation for a class of singular systems. Int. J. System Sci. 23, 517-530. Detchmendy, D. M. and Sridhar, R. (1966) Sequential estimation of state and parameters in noisy nonlinear dynamic systems. J. Basic Engng. 88, 362-368. Haber, J., Kozlowska, A. and Kozlowski, R. J. (1986) The structure and redox properties of vanadiumoxide surface-compounds. J. Catal. 102, 52-63. Jarvis, R. B. (1994) Robust Dynamic Simulation of Chemical Engineering Processes. PhD thesis, University of London. Joffe, B. L. and Sargent, R. W. H. (1972) The design of an on-line control scheme for a tubular catalytic reactor. Trans. Inst. Chem. Engrs 50, 270--282. Jonckheere, E. (1988) Variational calculus for descriptor problems. IEEE Trans. Auto. Control 33, 491-495. Jutan, A., Wright, J. D. and MagGregor, J. F. (1977) Multivariable computer control of a butane hydrogenolysis reactor: Part III. On-line linear quadratic stochastic control studies. AIChE J123, 751-758. Kershenbaum, L. S. and Lopez-Isunza, F. (1986) Measurement and estimation in catalytic reactors for partial oxidation. Trans. Inst. Msmt. Control 8, 137-143. Kim, I.-W., Liebman, M. J. and Edgar, T. E (1990) Robust error-in-variables methods for nonlinear dynamic systems. AIChE J136, 985-993. Kumar, S. and Seinfeld, J. H. (1978) Optimal location of measurements in tubular reactors. Chem. Engng Sci. 33, 1057-1516.

Kuruoglu, N., Ramirez, W. E and Clough, D. E. (1981) Distributed parameter estimation and identification for systems with fast and slow dynamics--a tubular fixed-bed catalytic reactor to form styrene. Chem. Engng Sci. 36, 1357-1364. Liebman, M. J., Edgar, T. E and Lasdon, L. S. (1992) Efficient data reconciliation and estimation for dynamic processes using nonlinear programming techniques. Computers chem. Engng 16, 963-986. Lopez-Isunza, E and Kershenbaum, L. S. (1992) The role of reversible changes in catalyst activity in the observed multiple steady states during partial oxidation dynamics. Chem. Engng Sci. 47, 2817-2822. Luenberger, D. G. (1977) Dynamic systems in descriptor form. IEEE Trans. Auto. Control 22, 312-32 I. Luenberger, D. G. (1978) Time-invariant descriptor systems. Automatica 14, 473-480. Mongkhonsi, T. (1993) Dynamic Behaviour of a FixedBed Catalytic Reactor with Catalyst Deactivation. PhD thesis, University of London. Mongkhonsi, T., Abi, C. F. and Kershenbaum, L. S. (1994) Dynamic behaviour of a fixed-bed catalytic reactor with catalyst deactivation. In AIChE Annual Meeting, Paper 88g, San Francisco. Nikolov, V. A., Klissurski, D. G. and Hadjiivanov, K. I. (1987) Deactivation of a V2Osfrio2 catalyst for the oxidation of o-xylene to phthalic anhydride. In Delmon, B. and G. E Froment, editors, Catalyst Deactivation, Elsevier, Amsterdam, 173-182. Nikoukhah, R., Willsky, A. S. and Levy, B. C. (1992) Kalman filtering and Riccati equations for descriptor systems. IEEE Trans. Auto. Control 37, 1325-1342. Pantelides, C. C. and Barton, E I. (1993) Equationoriented dynamic simulation current status and future perspectives. Computers chem. Engng 17, 263-285. Pantelides, C. C., Critsis, D., Moison, K. R. and Sargent, R. W. H. (1988) The mathematical-modelling of transient systems using differential algebraic equations. Computers chem. Engng 12, 449-454. Ramamurthi, Y., Sistu, E B. and Bequette, B. W. (1993) Control-relevant dynamic data reconciliation and parameter estimation. Computers chem. Engng 17, 41-59. Ray, H. W. (198 l) Advanced Process Control, McGrawHill Book Company, 2nd edition. Ray, W. H. (1978) Some recent application of distributed parameter systems theory--a survey. Automatica 14, 281-287. Robertson, D. G. and Lee, J. H. (1995) A least squares formulation for state estimation. J. Proc. Cont. 5, 291-299. Seinfeld, J. H. (1969) Nonlinear estimation for partial differential equations. Chem. Engng Sci. 24, 75-83. Siliva, J. M., Wallman, P. H. and Foss, A. S. (1979) Multi-bed catalytic reactor control systems: configuration development and experimental testing. Ind. Engng Chem. Fund. 18, 383-391. Sorenson, J. E (1977) Experimental investigation of the optimal control of a fixed-bed reactor. Chem. Engng Sci. 32, 763-774. Tzafestas, S. G. (1978) Distributed parameter state estimation. In Ray, W. H. and Laintiotis, D. G.,

Sequential estimation for nonlinear differential and algebraic systems editors. Distributed Parameter Systems, chapter 3, pages 135-207, Dekker. Verghese, G. C., Levy, B. C. and Kailath, T. (1981) A generalised state-space for singular systems. IEEE Trans. Auto. Control 26, 811-831. Wallman, P. and Foss, A. S. (1981) Experiences with

1067

dynamic estimators for fixed-bed reactors. Ind. Engng Chem. Fund. 20, 234-239. Windes, L. C., Cinar, A. and Ray, W. H. (1989) Dynamic estimation of temperature and concentration profiles in a packed bed reactor. Chem. Engng Sei. 44, 2087-2106.