Parameter Estimator Design with Application to a Chemical Reactor

The parameter estimator includes a left inverse of process model and at each time ... application and performance of the parameter estimation method a...
0 downloads 0 Views 130KB Size
Download PDF
Ind. Eng. Chem. Res. 1998, 37, 455-463

455

Parameter Estimator Design with Application to a Chemical Reactor Srinivas Tatiraju and Masoud Soroush* Chemical Engineering Department, Drexel University, Philadelphia, Pennsylvania 19104

A new method of on-line parameter estimation for a class of nonlinear systems is presented in this paper. The parameter estimator includes a left inverse of process model and at each time instant calculates least-squared-error estimates of parameters by using readily available online measurements. The estimator is computationally efficient and easy to implement. The application and performance of the parameter estimation method are illustrated and compared with those of parameter estimation via state estimation, by using a jacketed chemical reactor example for which the reaction-heat rate and overall heat-transfer coefficient are estimated. The simulation results show that the inversion-based estimator outperforms the state-estimationbased parameter estimator in the presence of measurement noise and plant-model mismatch. 1. Introduction In every physically meaningful process, sufficient data on important process parameters and state variables are essential for an efficient operation of the process. For example, accurate information on process parameters such as overall heat-transfer coefficient and reactionheat rate can ensure safe operation of exothermic chemical reactors. State variables of a process determine uniquely the state of the process and are either measured directly or estimated using a state estimator. On the other hand, process parameters provide a mathematical model with flexibility to fit process measurements, are often of great physical importance, and are usually not measured directly. Information on unknown process parameters can be obtained indirectly by means of a parameter estimator. In addition to its application in process monitoring, parameter estimation is an essential component of every adaptive control scheme. The general problem of parameter estimation is that of “fitting” a model to a set of measurements. It involves postulating a model with some unknown parameters and then calculating the values of the parameters that render the model-predicted values of process outputs “closest” to the measured values of the process outputs. In off-line parameter estimation (as in linear regression), a model is fitted “optimally” to the process measurements from one or several completed process runs (Gallant, 1987). However, in on-line parameter estimation, a model is fitted optimally to the past and present process measurements while the process is in operation (Narendra and Annaswamy, 1988; Bastin and Dochain, 1990; Slotine and Li, 1991; Yin, 1993). Existing methods of on-line parameter estimation include the following: (a) Prediction-error-based estimation methods such as the gradient estimator, the standard least-squares estimator, and least-squares estimator with exponential forgetting (Narendra and Annaswamy, 1988; Bastin and Dochain, 1990; Slotine and Li, 1991). * To whom correspondence should be addressed. E-mail: [email protected]. Telephone: (215) 895-1710. Fax: (215) 895-5837.

(b) Parameter estimation via state estimation. This method requires a dynamic model for each of the unknown parameters to be estimated. Simple parameter models such as “random walk” and “randon ramp” have usually been used in chemical and biochemical engineering (Stephanopoulos and San, 1984; Kosanovich et al., 1995). Once appropriate parameter models are chosen, a state estimator, such as an extended Kalman filter (de Valliere and Bonvin, 1989; Gudi et al., 1994; Elicabe et al., 1995; Semino et al., 1996; Sirohi and Choi, 1996) or a reduced-order observer (Soroush, 1997), is used to estimate the process parameters that appear as a subset of the state variables of the combined process and parameter models. This method has been used widely in chemical and biochemical engineering. (c) Parameter estimation via on-line optimization. The approach involves formulation of the parameter estimation problem as a minimization problem, wherein the parameter values are obtained by solving on-line a minimization problem such as the sum of the squared errors (Muske and Rawlings, 1995; Robertson et al., 1996). This method is difficult to implement and is computationally expensive. However, it allows one to include constraints in the estimation. (d) Calorimetric methods, in which process thermal and kinetic parameters are estimated on-line via simple mass and energy balances (Schuler and Schmidt, 1992; Carloff et al., 1994; Barandiaran et al., 1995; Regnier et al., 1996). These methods are quite popular in the process industry and now are available in the form of commercial software packages. In this paper, an inversion-based parameter estimation method for a class of nonlinear systems is presented. The application and performance of this estimation method are shown and compared with those of parameter estimation via state estimation, using a jacketed stirred-tank reactor example. The parameters representing the reaction-heat rate and overall heattransfer coefficient are estimated using available measurements of reactor and jacket temperatures. Parts of this paper were presented in the work by Tatiraju and Soroush (1997). The scope of this paper is first described, followed by a brief review of the method of parameter estimation via state estimation. The inversion-based parameter

S0888-5885(97)00536-8 CCC: $15.00 © 1998 American Chemical Society Published on Web 01/09/1998

456 Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998

estimation method is presented in section 3. The application and performance of the proposed parameter estimation method are illustrated by a chemical reactor example in section 4. 2. Scope and Preliminaries We consider the general class of nonlinear processes described by a linearly-parametrized mathematical model of the form

y˘ ) f(y,u) + g(y,u) p

(1)

where y ) [y1, ..., yn]T, u ) [u1, ..., um]T, and p ) [p1, ..., pq]T represent the vectors of measurable outputs, measurable inputs, and process parameters, respectively, f:Rn × Rm f Rn is a smooth vector function, and g:Rn × Rm f Rn × Rq is a smooth matrix function. The parameters p1, ..., pq are unknown and can be timevarying. It is assumed that the n × q matrix g(y,u) has a generic rank of q on Rn × Rm; that is, rank[g(y,u)] ) q, on a subset of Rn × Rm. The rank assumption requires the number of measurable outputs to be higher than or equal to the number of process parameters. Remark 1: A general mathematical model of the form

y˘ ) F (y,u,p˜ )

(2)

can always be recast in the form of the linearlyparametrized model of eq 1. For this reason in this paper we only focus on the class of process models of the form of eq 1. It is worth mentioning that, by recasting a model of the form of eq 2 in the form of eq 1, the nonlinearities are hidden in the new parameters, but they should be dealt with when a set of nonlinear algebraic equations have to be solved to calculate the original model parameters p˜ 1, ..., p˜ q from p1, ..., pq. This requires the nonsingularity of the q × q matrix ∂p/∂p˜ . 2.1. Parameter Estimation via State Estimation. In addition to the model of eq 1 which describes the dependence of the measurable outputs on the process parameters, this parameter estimation approach requires a dynamic model to describe the variations of the process parameters p1, ..., pq over a very short time horizon into the future. Let us assume that the process parameters are governed by a general model of the form

p˘ ) F(y,u,t,p)

(3)

where t represents time and F(., ., ., .) is a smooth vector function. This parameter model together with the process model of eq 1 takes the form

{

p˘ ) F(y,u,t,p) y˘ ) f(y,u) + g(y,u) p

(4)

Simple parameter models such as p˘ i ) 0 and p˘ i ) ait, i ) 1, ..., q, where a1, ..., aq are constants, are usually used. If the preceding augmented system is observable from the measurable outputs, y1, ..., yn, then a state estimation method with locally asymptotically stable error dynamics and with an adjustable rate of convergence can be employed to estimate the values of these parameters (which are a subset of the state variables of the augmented system of eq 4). For example, if the nonlinear reduced-order state observer design method presented in the work of Soroush (1997) is used, we obtain a parameter estimator of the form:

{

z˘ ) F(y,u,t,z+Ly) - L[f(y,u) + g(y,u) (z + Ly)] pˆ ) z + Ly

(5)

where pˆ is the estimated value of p, and L is the estimator gain. The corresponding observer error dynamics are governed by

e˘ ) F(y,u,t,e+p) - F(y,u,t,p) - Lg(y,u) e

(6)

where e ) pˆ - p. The constant observer gain, L, should be chosen such that the observer error dynamics of eq 6 are asymptotically stable. 3. Parameter Estimation via Inversion Let yˆ represent the value of the vector of measurable outputs when the process model is driven by an estimated value of the vector of process parameters pˆ ; that is, yˆ is governed by

yˆ˘ ) f(yˆ ,u) + g(yˆ ,u) pˆ

(7)

Comparing the dynamic system of eq 1 with that of eq 7, we see that, if yˆ f y asymptotically, then pˆ f p asymptotically when the matrix g is full rank. Thus, if there is a mechanism for forcing yˆ to go to y asymptotically, then the estimated value of the vector of process parameters, pˆ , will converge to its actual value, p, asymptotically. This is an analogue of synthesizing a model-based controller wherein, by forcing the controlled outputs to go to their set-point values and using a process model inverse, the controller action is calculated (Valluri et al., 1997). We request each yˆ l to be governed by a dynamic system of the form

γlyˆ˘ l + yˆ l ) yl where γl is an adjustable positive scalar constant, which ensures that each yˆ l goes to yl asymptotically at an adjustable rate set by γl. Thus, by the process model of eq 7

{}

diag

1 [y - yˆ ] ) f(yˆ ,u) + g(yˆ ,u) pˆ γl

Since rank(g) ) q, the q × q matrix gTg is invertible, and using the preceding algebraic equation, we can solve for pˆ :

( {}

pˆ ) [g(yˆ ,u)Tg(yˆ ,u)]-1g(yˆ ,u)T diag

)

1 [y - yˆ ] - f(yˆ ,u) γl

The following theorem presents the resulting parameter estimator and describes precisely the theoretical properties of the estimator. Theorem 1: For a process of the form of eq 1 with rank[g(y,u)] ) q, ∀(y, u) ∈ Rn × Rm, the parameter estimator

{

1 yˆ˘ 1 ) (y1 - yˆ 1) γ1 l 1 yˆ˘ n ) (yn - yˆ n) γn

(8)

pˆ ) [g(yˆ ,u)Tg(yˆ ,u)]-1g(yˆ ,u)T × 1 diag [y - yˆ ] - f(yˆ ,u) γl

( {}

)

has the following properties: (a) pˆ f p asymptotically, in the limit that the values of the parameters γ1, ..., γr all go to zero.

Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998 457

of the parameters that render yˆ as close as possible to y. This appealing property is a consequence of the fact that eq 8 is an nth-order state-space realization of the left inverse of the dynamic system of eq 1 and an nthorder filter.

a

b

The parameter estimator of eq 8 cannot be used for a process of the form of eq 1 that crosses a singular manifold on which rank[g(y,u)] < q, although the process has a g(y,u) matrix with a generic rank of q. To address this issue, the parameter estimation problem is formulated in the form of the following optimization problem that has an analytical solution:

c

Figure 1. (a) Right process inverse. (b) Left process inverse. (c) Inversion-based parameter estimator.

(b) pˆ f p asymptotically, if y˘ 1, ..., y˘ r f 0 asymptotically. Proof: Part (a): The vector of the predicted measurable outputs, yˆ , is governed by

{}

1 yˆ˘ ) diag (y - yˆ ) γl

(9)

and thus in the limit that γ1, ..., γr f 0, after an infinitesimal period of time, yˆ ) y. The quantities y and yˆ also satisfy (by eqs 1 and 7)

g(yˆ ,u)Tg(yˆ ,u) pˆ ) g(yˆ ,u)T[yˆ˘ - f(yˆ ,u)] g(y,u)Tg(y,u) p ) g(y,u)T[yˆ˘ - f(y,u)] The q × q matrix gTg is nonsingular ∀(y, u) ∈ Rn × Rm and after an infinitesimal period of time y ) yˆ , and therefore pˆ ) p after an infinitesimal period of time. Part (b): In the case that y˘ 1, ..., y˘ r f 0 asymptotically, since the linear system of eq 9 is asymptotically stable, yˆ f y asymptotically and thus pˆ f p asymptotically. QED The parameter estimator of eq 8 inherently includes a low-pass filter for each of the measurable outputs: each yˆ l is indeed the output of a first-order low-pass filter which has a time constant of γl and whose input is the measurable output yl. A more noisy measurement of yl requires a higher value for γl, and a higher value of γl attenuates more the noise in the yl measurement. The inversion-based estimator calculates the parameter estimates on the basis of filtered measurable outputs. The inversion-based parameter estimator of eq 8 includes a left inverse of the process. The notions of right and left inverses are depicted in parts a and b of Figure 1, respectively. As in the case of a linear static process described by matrices, a left inverse and/or a right inverse of a dynamic process may not exist. In the case of nonlinear processes of the form of eq 1, a right inverse may not exist. In model-based control, a controller includes an inverse of a process model to calculate the desired values of manipulated inputs on the basis of controlled output measurements. Note that a model-based controller usually incorporates a right inverse of a process model. A block diagram of the inversion-based parameter estimator of eq 8 is depicted in Figure 1c. At each time instant the inversion-based parameter estimator calculates least-squared-error estimates of the parameters; at each time instant it provides the values

min (pˆ -p j)

{|

{} |

f(y,u) + g(y,u)pˆ - diag

1

γl

(y - yˆ )

2

+

}

q

βl(pˆ l - p j l)2 ∑ l)1

where p j is the nominal value of the vector of process parameters, |‚| denotes the Euclidean norm, and β1, ..., βq are positive parameters. Note that, if the deviation of pˆ from p is not penalized in the preceding optimization problem (if βl ) 0, l ) 1, ..., q) and the matrix g is full rank everywhere, then the solution to the preceding optimization problem will be given by eq 8. This approach to addressing the singular points where the matrix g is not full rank is an analogue of the approach to synthesizing an input-output linearizing controller for a process with a locally singular characteristic matrix, described in the work of Valluri et al. (1997). Solving the preceding optimization problem results in the parameter estimator given in the following remark. Remark 2: Consider a process of the form of eq 1 whose matrix g(y,u) has a generic rank of q but loses a rank on a manifold and whose nominal steady-state pair, denoted by (yj, u j ), is not on the singular manifold. For such a process, one can use the parameter estimator:

{}

1 yˆ˘ ) diag (y - yˆ ) γl pˆ ) [g(yˆ ,u)Tg(yˆ ,u) + diag{βl}]-1g(yˆ ,u)T × 1 diag [y - yˆ ] - f(yˆ ,u) - g(yˆ ,u) p j +p j (10) γl

( {}

)

where the values of the scalar tunable parameters β1, ..., βq should be chosen such that the matrix (gTg + diag{βl}] is nonsingular for every (y, u) ∈ Rn × Rm. The values of these tunable parameters are recommended to be set according to βl ) gll(yj,u j )2/10000, l ) 1, ..., q. 3.1. Implementation of the Parameter Estimator. In order to implement the parameter estimator on a digital computer, a discrete-time version of the parameter estimator is needed. The dynamic component of the estimator is linear and thus it can be timediscretized exactly. The exact time-discretized versions of the parameter estimators of eqs 8 and 10 respectively have the following forms:

458 Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998

yˆ (k+1) ) diag{λl}yˆ (k) + [I - diag{λl}]y(k)

(11)

pˆ (k) ) [g(yˆ (k),u(k))Tg(yˆ (k),u(k))]-1g(yˆ (k),u(k))T × 1 diag (y(k) - yˆ (k)) - f(yˆ (k),u(k)) γl

[ {}

]

yˆ (k+1) ) diag{λl}yˆ (k) + [I - diag{λl}]y(k)

Schuler and Schmidt, 1992; Barandiaran et al., 1995; Elicabe et al., 1995; Kosanovich et al., 1995; Regnier et al., 1996; Semino et al., 1996). 4.1. The Reactor. We consider the continuous stirred-tank reactor depicted in Figure 2, in which the following irreversible parallel reactions take place: k1

(12)

A 98 U1

pˆ (k) ) [g(yˆ (k),u(k))Tg(yˆ (k),u(k)) + 1 diag{βl}]-1g(yˆ (k),u(k))T diag (y(k) - yˆ (k)) γl

A 98 U2

[ {}

]

k2

kP

A 98 P

f(yˆ (k),u(k)) - g(yˆ (k),u(k)) p j +p j where ∆t is the sampling period of the measurements, y(k) and u(k) are respectively the measured values of y and u at a time instant k, yˆ (k) is the value of yˆ at the time instant k, and λl ) exp(-∆t/γl), l ) 1, ..., n.

where P is the desired product, and U1 and U2 are the undesirable side products. The reactor model has the form

C˙ A ) -k1CA3 - k2CA0.5 - kPCA +

4. Application to a Chemical Reactor For tight temperature control and reliable monitoring of extremely exothermic chemical reactors, such as bulk polymerization reactors, accurate information on the reaction-heat-production rate and overall heat-transfer coefficient of the reactor/cooling equipment is essential (Schnelle and Richards, 1986). There are no sensors that can measure these important parameters directly. The rate of heat of reaction strongly depends on the type and purity of reactants, and the overall heat-transfer coefficient usually decreases with time because of fouling and/or an increase in the viscosity of the reacting mixture. Features such as these have motivated a great number of studies on on-line estimation of these parameters (Wu, 1985; de Valliere and Bonvin, 1989;

C˙ P ) kPCA T˙ )

T˙ j )

CAi - CA τ

CP τ

H1k1CA3

+

(13) H2k2CA0.5 Frcr

+ HPkPCA

+

Ti - T SU (T - T) + τ FrcrVr j

SU Q (T - Tj) + mjcj mjcj

where the reaction rate constants k1 ) z1 exp(-E1/RT), k2 ) z2 exp(-E2/RT), and kP ) zP exp(-EP/RT). In this reactor, the reactor-jacket overall heat-transfer coefficient is assumed to decrease linearly with the product concentration according to

U ) U0(1 - 0.15CP)

Figure 2. Schematic diagram of the simulated reactor and estimator system. Table 1. Parameter Values and Operating Conditions of the Chemical Reactor Example cr 4.2 kJ‚kg-1‚K-1 cj 4.2 kJ‚kg-1‚K-1 mj 3.0 kg Vr 1.0 × 10-2 m3 U0 5.0 × 10-1 kJ‚m-2‚K-1‚s-1 S 2.0 × 10-2 m2 S∞ 2.0 × 10-2 m2 U∞ 2.1 × 10-1 kJ‚m-2‚K-1‚s-1 T∞ 3.00 × 102 K T(0) 2.90 × 102 K Tj(0) 3.00 × 102 K CA(0) 1.0 × 10-1 kmol‚m-3 CP(0) 0.0 × 100 kmol‚m-3 Qss -1.013 kJ‚s-1

where U0 is the value of the reactor-jacket overall heattransfer coefficient when the concentration of the product in the reactor is zero. This decrease in the overall heat-transfer coefficient is quite common in many reactors wherein generation of viscous products causes a decrease in the overall heat-transfer coefficient (Chylla and Haase, 1993; Soroush and Kravaris, 1992b). The parameter values of the reactor model are the same as those in the work of Soroush and Kravaris (1992a), unless given in Table 1. Steady state analysis of the reactor shows that, at a steady-state temperature of 400 K, the steady-state concentration of the desirable product is maximum (Soroush and Kravaris, 1992a). We use the following proportional-integral (PI) controller to maintain the reactor temperature at 400 K by manipulating the heat input to the jacket (Q):

η˘ ) 400 - T Q ) Kc(400 - T) +

Kc η + Qss τI

(14)

with Kc ) 0.1 kJ‚K-1‚s-1 and τI ) 1000 s. The dynamic system of eq 13 along with the controller of eq 14 is used to represent the actual process. An integration step size

Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998 459

of 5 s is used in the numerical integrations of the differential equations of the model and controller. 4.2. Estimation Problem. The objective is to estimate on-line the rate of heat release by the reactions and the overall heat-transfer coefficient. In particular, it is desirable to estimate the following two parameters:

H1k1CA3 + H2k2CA0.5 + HPkPCA p1 ) , Frcr SU0(1 - 0.15CP) p2 ) (15) FrcrVr The uncertain parameters p1 and p2 include the rate of heat production by the reactions and the overall heattransfer coefficient, respectively. The model of eq 13 can be recast in the following form:

{

T˙ ) p1 + p2(Tj - T) + R1 T˙ j ) R2p2(T - Tj) + R3

Ti - T , τ

R2 )

FrcrVr , mjcj

R3 )

Q mjcj

Thus, given the temperature measurements T, Ti, and Tj and the rate of heat input Q (i.e., known values of Ri, i ) 1, 2, 3), we would like to estimate the process parameters (p1 and p2). Equation 16 is of the form of eq 1 with

y ) [T

Tj]T,

[ ]

R f(y,u) ) R1 , 3

u ) [R1 g(y,u) )

[

1 0

R2

{

T˙ ) p′1 + p′2 + R1 T˙ j ) -R2p′2 + R3

R3]T

(Tj - T) R2(T - Tj)

]

4.3. Parameter Estimation via Inversion. In this example, the matrix g(y,u) is generically nonsingular (rank(g) ) 2 ) q), but it is singular where T ) Tj (rank(g) ) 1 * q when T ) Tj). Because the reactor will cross the singular manifold during the temperature control, we implement the discrete-time version of the parameter estimator of Remark 2, given by eq 12, with β1 ) 0. The discrete-time estimator has the form

T ˆ (k+1) ) λ1T ˆ (k) + [1 - λ1]T(k),

T ˆ (0) ) T(0)

ˆ j(k) + [1 - λ2]Tj(k), T ˆ j(k+1) ) λ2T

T ˆ j(0) ) Tj(0)

pˆ 2(k) ) Tj(k) - T ˆ j(k) - γ2[R3(k) + R2p j 2(T ˆ (k) - T ˆ j(k))] × γ2 ˆ (k) - T ˆ j(k)) R2(T +p j2 2 β2 + R2 (T ˆ (k) - T ˆ j(k))2 T(k) -T ˆ (k) - γ1[(T ˆ j(k) - T ˆ (k))pˆ 2(k) + R1(k)] pˆ 1(k) ) γ1 (17) where λ1 ) exp(-∆t/γ1) and λ2 ) exp(-∆t/γ2). We assume that the measurements are available every 5 s. The tunable parameters γ1 and γ2 of the estimator are chosen according to the guidelines described in section 3. However, in the presence of noise in the measurable outputs these values have to be

(18)

where p′1 ) p1 and p′2 ) p2(Tj - T). Assuming that the parameter p′1 and p′2 are constant within a small interval of time, we can write the following models for the process parameters:

p˘ ′2 ) 0

p˘ ′1 ) 0,

(16)

where

R1 )

increased to reduce the effect of the noise on the parameter estimates. Note that the first two equations of eq 17 represent two first-order low-pass filters; the higher the values of γ1 and γ2, the more the attenuation of the noise in the temperature measurements. The other tunable parameter, β2, is set to 0.1 K2 in all the simulation studies. 4.4. Parameter Estimation via State Estimation. In order to apply the state estimation method, we first parametrize the model into the form

(19)

The augmented system consisting of the systems of eqs 18 and 19 is observable from the measurable outputs T and Tj, and thus we use the reduced-order observer described in the work of Soroush (1997) to estimate the two state variables p′1 and p′2 only. The observer is of the form

[] [

z˘ 1 -L11 z˘ 2 ) -L21

][ ]

-L11 + R2L12 z1 -L21 + R2L22 z2 + [L11 + L21]T + [L12 + L22]Tj + R1 -R2[L21T + L22Tj] + R3

[

]

pˆ ′1 ) z1 + L11T + L12Tj pˆ ′2 ) z2 + L21T + L22Tj

(20)

The corresponding estimator error dynamics are given by

e˘ 1 ) -L11e1 - e2[L11 - L12R2] e˘ 2 ) -L21e1 - e2[L21 - L22R2]

(21)

where e1 ) pˆ 1′ - p′1 and e2 ) pˆ 2′ - p′2. The gain matrix [Lij] should be chosen such that the error dynamics are asymptotically stable; that is,

L22R2 < L11 + L21,

L11L22 < L21L12

In all the simulation studies the observer gain components are set according to L11 ) 0.10, L21 ) 0.00, L12 ) 0.33, and L22 ) -0.10. The continuous-time system of eq 20 is linear in z, and thus it can be time-discretized exactly:

z(k+1) ) exp(∆tA)z(k) + [exp(∆tA) - I]A-1b(k) pˆ ′1(k) ) z1(k) + L11T(k) + L12Tj(k)

(22)

pˆ ′2(k) ) z2(k) + L21T(k) + L22Tj(k) where

[

-L A ) -L11 21

]

-L11 + R2L12 -L21 + R2L22 , [L11 + L21]T(k) + [L12 + L22]Tj(k) + R1 b(k) ) -R2[L21T(k) + L22Tj(k)] + R3(k)

[

]

460 Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998

The preceding time-discretized version of the parameter estimator of eq 20 is used in the simulation studies. In this approach, due to the reduced-order form of the observer, the parameter estimates are directly dependent on the temperature measurements (the last two equations of the system of eq 22 have the T and Tj measurements as input); any noise in the temperature measurements will directly affect the estimates of the parameters. Therefore, noisy measurements should be filtered before they are fed to the estimator. In the simulation case studies (presented in the next subsection) the noise in each measurement is attenuated by a first-order filter with a time constant of 0.9 s, and then the filtered measurements are fed to the estimators. 4.5. Results and Discussion. The performances of the two estimators (eqs 17 and 22) are studied in the following four cases: Case I: No measurement noise and no plant-model mismatch (nominal case). Case II: White noise (with a standard deviation of 0.5 K) in the measurable outputs (T and Tj). Case III: Model-plant mismatch. Case IIIA: A term representing heat loss from the reactor to the surrounding is added to the system of eq 13 without changing the parameter estimator. Case IIIB: The efficiency of the jacket heater/cooler is set to 90%; which is not accounted for in the parameter estimator. Case IV: Both the white noise of case II and the model-plant mismatch of case III. In all these case studies, the process is represented by a combination of eqs 13 and 14, the reactor and jacket temperatures are assumed to be sampled every 5 s, and in the state-estimation approach, the initial estimate of each of the parameters is set to 50% of the actual value of the corresponding parameter. Note that no such initialization is needed for the inversion-based estimation method. In all the simulation studies the performances of the two parameter estimators are presented and compared on the basis of the relative error between the estimated and actual parameter values. The relative error, denoted by Rei, is defined as

Rei )

pˆ i - pi , pi

i ) 1, 2

4.5.1. Case I: The Nominal Case. Figure 3 depicts the performances of the inversion-based parameter estimator and the state estimation approach for this case. The values of the tunable parameters (γ1 and γ2) of the inversion-based estimator are both set to 0.8 s. It can be clearly seen that the relative errors, Re1 and Re2, both decay to zero, indicating the rapid convergence of estimates of the parameters to their actual values. As Figure 3 shows, in the neighborhood of the point where T ) Tj, at t ≈ 680 s, Re2 deviates slightly from zero (for both estimators). This is due to the fact that when T ) Tj, according to the model representing temperature dynamics, the temperature measurements do not include any information on the unknown process parameter p2 (the term containing p2 vanishes in the system of eq 16). In such a situation, for the inversionbased estimator, it is suggested to use the previous estimate of the parameter p2 when |Tj - T| < 1.0; that is, pˆ 2(k) ) pˆ 2(k-1) when |Tj - T| < 1.0. Once this manifold is crossed, again estimates of the parameters

Figure 3. Actual values of the process parameters and the relative errors in the estimated values for the nominal case (case I).

are found to correspond to the actual values. The error in the estimate of the parameter p1 is not as much as that in the estimate of p2 (deviation in Re2 is much more pronounced than that in Re1). Similarly in the case of the state-estimation-based parameter estimator, it can be seen that the estimate of p2 deviates from its actual value at t ≈ 680 s (when |Tj - T| < 1.0). In this case, when T ) Tj, the parameter (state variable) p2 is not observable from the temperature measurements. However, under such conditions p1 is observable, and hence there is no deviation in its estimate. Figure 3 also shows that the parameter representing the overall heattransfer coefficient (p2) reduces by 60% of its initial value during the course of the reaction. 4.5.2. Case II: Noise in the Measurable Outputs. The noisy reactor temperature measurement gives rise to a noisy controller output (Q). Thus, in this case both measurable outputs and one measurable input are corrupted with noise. The values of tunable parameters (γ1 and γ2) of the inversion-based estimator are both increased to 0.95 s; larger values of these two tunable parameters attenuate more the noises in the temperature measurements, resulting in less noisy estimates of the parameters. For this case, the performances of the two parameter estimators are depicted in Figure 4. The relative errors in the parameter estimates rapidly converge to zero, thereby indicating that the estimates of the parameters correspond to their actual values. In this case, also the estimates of the parameter p2 calculated by the inversion-based and state-estimation-based estimator deviate slightly from their actual values in the neighborhood of the point where T ) Tj. This deviation was explained sufficiently in the discussion of case I. A procedure similar to that explained earlier is followed in estimat-

Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998 461

Figure 4. Actual values of the process parameters and the relative errors in the estimated values in the presence of the measurement noise (case II).

Figure 6. Actual values of the process parameters and the relative errors in the estimated values in the presence of 90% jacket-heater efficiency (case IIIB).

Figure 5. Actual values of the process parameters and the relative errors in the estimated values in the presence of heat loss to the ambient from the reactor (case IIIA).

Figure 7. Actual values of the process parameters and the relative errors in the estimated values in the presence of heat loss from the reactor and 90% jacket-heater efficiency (case III).

ing the parameters in the neighborhood where T ) Tj. The state-estimation-based approach exhibits a higher

sensitivity to the measurement noise, since the noisy temperature measurements (though filtered) are di-

462 Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998

rectly used in calculating the estimates. However, the inversion-based estimator filters out the noise and gives less noisy T ˆ and T ˆ j which are then used to calculate the parameters. 4.5.3. Case III: The Model-Plant Mismatch. In order to study the performances of the two estimation methods in the presence of the model-plant mismatch, the following changes are made in the process without modifying the parameter estimators of eqs 17 and 20: (A) A heat loss term accounting for heat dissipation from the reactor to the ambient atmosphere is added to the system of eq 13: the third equation of the system of eq 13 is replaced with

T˙ )

H1k1CA3 + H2k2CA0.5 + HPkPCA Ti - T + + Frcr τ S∞U∞ SU0(1 - 0.15CP) (Tj - T) + (T - T) (23) FrcrVr FrcrVr ∞

Figure 5 depicts the performances of the two estimation methods in the presence of this plant-model mismatch. In this case, the estimated value of p1 is underestimated asymptotically, because the heat loss term is always negative during the operation. However, the plantmodel mismatch does not cause a considerable deviation in the estimate of the parameter p2. (B) The heat input to the process is 90% of the controller output (Q): the fourth equation of the system of eq 13 is replaced with

T˙ j )

SU0(1 - 0.15CP) 0.9Q (T - Tj) + mjcj mjcj

(24)

without accounting for this change in the parameter estimators. The performances of the two estimation methods in the presence of this plant-model mismatch are shown in Figure 6. It is evident that the estimates of the parameter p1 converge to the actual value since Re1 rapidly converges to zero. However, the estimates of p2, representing the overall heat-transfer coefficient, deviate slightly from the actual value. Figure 7 shows the estimator performances when both the aforementioned model errors are present. In this case, the estimates of both Re1 and Re2 converge to values slightly different from zero, since the parameters calculated by both estimators differ from their actual values. Again the inversion-based estimator outperforms the other estimator; it provides least-squarederror estimates of the parameters. In all three model error case studies, γ1 and γ2 of the inversion-based estimator are both set to 0.8 s. 4.5.4. Case IV: Noise in Measurable Outputs and Plant-Model Mismatch. Figure 8 depicts the performance of the two estimation methods when the conditions of cases II and III are present together. In this case, the tunable parameters γ1 and γ2 of the inversion-based estimator are both set to 0.95 s. In this case also, the inversion-based estimator performs better in the presence of the measurement noise and plantmodel mismatch, because at each time instant the inversion-based estimator provides least-squared-error estimates of the parameters. In the presence of the process-model mismatch, parameter estimates with zero steady-state errors may be obtained, if we add adaptive capabilities to the estimators.

Figure 8. Actual values of the process parameters and the relative errors in the estimated values in the presence of the measurement noise and the process-model mismatch (case IV).

5. Conclusions A simple model-inversion-based estimation method for on-line parameter estimation in a class of nonlinear systems was presented. At each time instant the estimator provides least-squared-error estimates of parameters. The parameter estimator is easy to implement and computationally efficient. The implementation and performance of the estimator were demonstrated and compared with those of parameter estimation via state estimation by a chemical reactor example. Overall, the inversion-based parameter estimator was found to outperform the state estimation method in the presence of measurement noise and model errors. Acknowledgment Partial financial support from the National Science Foundation through Grant CTS-970-3278 is gratefully acknowledged. Acknowledgment is also made to the donors of the Petroleum Research Fund, administered by the ACS, for partial support of this research. Notation CA ) outlet concentration of reactant A, kmol‚m-3 CAi ) inlet concentration of reactant A, kmol‚m-3 CP ) outlet concentration of product P, kmol‚m-3 cr ) heat capacity of the reacting mixture, kJ‚kg-1‚K-1 cj ) heat capacity of the jacket fluid, kJ‚kg-1‚K-1 e ) observer error El ) activation energy for the undesired reactions, l ) 1, 2, kJ‚kmol-1 EP ) activation energy for the desired reactions, kJ‚kmol-1 Hl ) heat of reaction for the undesired reaction l, l ) 1, 2, kJ‚kmol-1

Ind. Eng. Chem. Res., Vol. 37, No. 2, 1998 463 HP ) heat of reaction for the desired reaction, kJ‚kmol-1 k1 ) rate constant for the undesired reaction 1, m6‚kmol-2‚s-1 k2 ) rate constant for the undesired reaction 2, m-1.5‚ kmol0.5‚s-1 kP ) rate constant for the desired reaction, s-1 L ) observer gain matrix mj ) mass of coolant in the reactor jacket, kg pi ) unknown process parameter, i ) 1, ..., m pˆ i ) estimated value of a parameter pi Q ) heat input to the reactor jacket, kJ‚s-1 Qss ) steady-state value of heat input to the reactor jacket, kJ‚s-1 R ) universal gas constant, kJ‚kmol-1‚K-1 Re1 ) relative error in the estimated value of p1, Re1 ) (pˆ 1p1)/p1 Re2 ) relative error in the estimated value of p2, Re2 ) (pˆ 2p2)p2 S ) jacket-reactor heat-transfer surface area, m2 S∞ ) jacket-surrounding heat-transfer surface area, m2 T ) reactor temperature, K Ti ) temperature of the inlet stream to the reactor, K Tj ) jacket temperature, K T∞ ) ambient temperature, K t ) time, s u ) vector of measurable inputs U0 ) overall jacket-reactor heat-transfer coefficient when CP ) 0, kJ‚m-2‚s-1‚K-1 U∞ ) overall ambient-reactor heat-transfer coefficient, kJ‚m-2‚s-1‚K-1 Vr ) reactor volume, m3 Vj ) reactor jacket volume, m3 y ) vector of measurable outputs yˆ ) predicted value of y z1 ) Arrhenius factor for the undersired reaction 1, m6‚kmol-2‚s-1 z2 ) Arrhenius factor for the undesired reaction 2, m-1.5‚kmol0.5‚s-1 zP ) Arrhenius factor for the desired reaction, s-1 Greek Letters βj ) estimator tunable parameter (j ) 1, ..., q) ∆t ) sampling period, s γj ) estimator tunable parameter (j ) 1, ..., q), s Fr ) density of the reaction mixture, kg‚m-3 τ ) reactor residence time, s Math Symbols diag{βl} ) q × q diagonal matrix with diagonal elements, β1, ..., βq diag{1/γl} ) n × n diagonal matrix with diagonal elements 1/γ1, ..., 1/γn diag{λl} ) n × n diagonal matrix with diagonal elements λ1, ..., λn gT ) transpose of a matrix g I ) n × n identity matrix

Literature Cited Barandiaran, M. J.; De Arbina, L. L.; De La Cal, J. C.; Gugliotta, L. M.; Ausa, J. M. Parameter Estimation in Emulsion Copolymerization using Reaction Calorimeter Data. J. Appl. Polym. Sci. 1995, 55, 1231. Bastin, G.; Dochain, D. Online Estimation and Adaptive Control of Bioreactors: Elsevier: New York, 1990.

Carloff, R.; Probss, A.; Reichert, K. H. Temperature Oscillation Calorimetry in Stirred Tank Reactors with Variable Heat Transfer. Chem. Eng. Technol. 1994, 17, 406. Chylla, R. W.; Haase, D. R. Temperature Control of Semibatch Polymerization Reactors. Comput. Chem. Eng. 1993, 17, 257. de Valliere, Ph.; Bonvin, D. Application of Estimation Techniques to Batch ReactorssI. Modeling Thermal Effects. Comput. Chem. Eng. 1989, 13 (1/2), 1. Elicabe, G. E.; Ozdeger, E.; Georgakis, C.; Cordeiro, C. On-line Estimation of Reaction Rates in Semicontinuous Reactors. Ind. Eng. Chem. Res. 1995, 34, 1219. Gallant, A. R. Nonlinear Statistical Models; John Wiley & Sons: New York, 1987. Gudi, R. D.; Shah, S. L.; Gray, M. R. Multirate State and Parameter Estimation in an Antibiotic Fermentation with Delayed Measurements. Biotechnol. Bioeng. 1994, 44, 1271. Kosanovich, K. A.; Piovoso, M. J.; Rokhlenko, V.; Guez, A. Nonlinear Adaptive Control with Parameter Estimation of a CSTR. J. Process Control. 1995, 5, 137. Muske, K. R.; Rawlings, J. B. Nonlinear Moving Horizon State Estimation. In Methods of Model Based Process Control; Berber, R., Ed.; Kluwer Academic Publisher: Amsterdam, The Netherlands, 1995. Narendra, K. S.; Annaswamy, A. M. Stable Adaptive Systems; Prentice-Hall: Englewood Cliffs, NJ, 1988. Regnier, N.; Defaye, G.; Caralp, L.; Vidal, C. Software Sensor Based Control of Exothermic Batch Reactors. Chem. Eng. Sci. 1996, 51, 5125. Robertson, D. G.; Lee, J. H.; Rawlings, J. A. A Moving HorizonBased Approach for Least-Squares State Estimation. AIChE J. 1996, 22 (8), 2209. Schnelle, P. D., Jr.; Richards, J. R. A Review of Industrial Reactor Control: Difficult Problems and Workable Solutions. In Chemical Process Control III; Morari, M., McAvoy, T. J., Eds.; Elsevier: New York, 1986; p 749. Schuler, H.; Schmidt, Ch.-U. Calorimetric State Estimators for Chemical Reactor Diagnosis and Control: Review of Methods and Applications. Chem. Eng. Sci. 1992, 47, 899. Semino, D.; Morretta, M.; Scali, C. Parameter Estimation in Extended Kalman Filters for Quality Control in Polymerization Reactor. Comput. Chem. Eng. 1996, 20, 913. Sirohi, A.; Choi, K. Y. On-line Parameter Estimation in a Continuous Polymerization Process. Ind. Eng. Chem. Res. 1996, 35, 1332. Slotine, J. J.; Li, W. Applied Nonlinear Control; Prentice-Hall: Englewood Cliffs, NJ, 1991. Soroush, M. Nonlinear State-Observer Design with Application to Reactors. Chem. Eng. Sci. 1997, 52, 387. Soroush, M.; Kravaris, C. Discrete-Time Nonlinear Controller Synthesis by Input/Output Linearization AIChE J. 1992a, 38, 1923. Soroush, M.; Kravaris, C. Nonlinear Control of a Batch Polymerization Reactor: an Experimental Study. AIChE J. 1992b, 38, 1429. Stephanopoulos, G.; San, K. Y. Studies on On-Line Bioreactor Identification. Biotechnol. Bioeng. 1984, 26, 1176. Tatiraju, S.; Soroush, M. Parameter Estimation via Inversion with Application to a Chemical Reactor. Proc. ACC 1997, 2429. Valluri, S.; Soroush, M.; Nikravesh, M. Shortest-PredictionHorizon Nonlinear Model Predictive Control. Chem. Eng. Sci. 1997, in press. Wu, R. S. H. Dynamic Thermal Analyzer for Monitoring Batch Processes. Chem. Eng. Prog. 1985, Sept, 57. Yin, K. W. Parameter Estimation and Related Problems for External Input-Output Models. Chem. Eng. Sci. 1993, 48, 3780.

Received for review July 31, 1997 Revised manuscript received November 12, 1997 Accepted November 13, 1997 IE9705362