- Email: [email protected]
Estimation of states and parameters by invariant imbedding technique
"JChemicalEnoineerino Science, Vol. 44, No. 8, pp. 1665 1674, 1989.
0009 2509/89 $3.00+0.00 © 1989 Pergamon Press plc
Printed in Great Britain.
ESTIMATION OF STATES A N D PARAMETERS BY INVARIANT IMBEDDING TECHNIQUE WU YUAN Department of Chemical Engineering, Wuhan Institute of Chemical Engineering, Wuhan, China (Received 24 October 1987; accepted for publication 28 November 1988) Abstract--A set of sequential estimator equations has been developed using the invariant imbedding
concept. The equations can be used for both the filteringestimation of states of non-linear dynamic systems and the estimation of parameters in kinetics research and other areas, and are of potential in simplifying experimental procedures and reducing the number of necessary experiments for the latter purpose. Results calculated for severalcases show that the convergenceof the estimated values for states, the stability of those for parameters, and the accuracy are all very good, suggestingthat the estimator equations may be expected to be a useful tool in chemical engineering and even in many other fields.
INTRODUCTION Kinetics research and process control are important in many areas. To solve them requires the estimation of states and parameters. Establishing a model is one of the main jobs in kinetics research and is an important basis for process analysis, design, scale-up, operation and control. However, the experimental investigation for this purpose is usually complex, difficult, expensive and timeconsuming. Therefore, it is useful to develop new methods for data interpretation which require less data and may allow a simplified experimental procedure. Since the 1970s, great attention has been devoted to "sequential experiment planning" [e.g. Rose (1981) and Diamond (1981)], for which the basic idea is fully using the information provided by the last experiments to appropriately plan the next ones in order to enable each experiment to give a high responsibility and thus a high efficiency. However; the experimental procedure required by this method is also considerably complex. The correct estimation of states is the prerequisite for effective process control. Therefore, the filtering estimation of states of dynamic systems is continuously one of the main research projects in process control. It is general for non-linear systems to "linearize" them, as the linear differential equations may easily be solved with the Laplace transforms (Stephanopoulos, 1984). However, the linearization usually introduces certain errors into the system, and, in some cases, the errors introduced may exceed the acceptable range. For the estimation of the states of non-linear systems themselves, a number of methods have been developed (Omatu, 1986), in which the real physical systems are usually described with stochastic models, that is, each of the equations describing dynamic behavior of the states includes a white noise excitation. Such a treatment may yield a more accurate analysis, but it also creates difficulties in computation and concept as well (Bullin and Dukler, 1975). On the other hand, a definite (non-stochastic)
model may be accurate enough for a number of processes in chemical engineering. Therefore, it is also of interest to study the filtering estimation of the states of non-linear dynamic systems with non-stochastic models. The invariant imbedding technique is one of the effective methods for solving boundary value problems of non-linear differential equations, for which certain applications in chemical engineering have been developed (Lee, 1968a; Bellman and Lee, 1977). Lee (1968b, c) used the invariant imbedding concept for non-linear filtering and parameter estimation. Unfortunately, because of omissions in the mathematical treatment, the estimator equations Lee developed are very complex, with no definite initial values, and thus have little value for practical applications. In fact, no successful application or progress in improving Lee's equations has been reported. In the present paper, a set of new sequential filtering estimator equations will be developed, and two applications to kinetics research will be given. ONE-DIMENSIONALESTIMATOREQUATION Consider a system with a single state variable whose kinetic behavior can be described by the non-linear equation dx
- - = f ( x , t). dt
(1)
From the initial time, t =0, to the current instant, t = ty, the state x is measured, but random errors must exist. The value measured for state x at time t, z(t), can be represented by z(t) = x(t) + er(t)
(2)
where x is the true state, and er is the measuring error. The requirement is to estimate the true state x(t) using the measurement within the period from 0 to tf. Provided that the state x at a certain instant is exactly estimated, then the definite condition of eq. (1) is
1665
1666
W U YUAN
determined, and thus the value for x at any time can be easily calculated. If one first estimates the current state x(tl), then the problem becomes the filtering estimation. From the classical least-squares criterion, the estimated value for x(tl) should minimize the integral
J=f
[x(t)-z(t)]2dt.
(3)
Substituting c by e and using eq. (12), eq. (15) becomes
f(e, a)rcc(e, a) + rca(e, a) = 2 [ e - z(a)] rca(e, a) ----=f(e, rcc(e, a)
a) -~
2[z(a)- e]
(17)
roe(e, a)
Combining eq. (17) with eq. (14):
de 2[z(a)-e] ~a=f(e,a)+ rcc(e,a~
Defining y(t) = I t [ x ( t ) - z(t)] 2 dt 30
(4)
gives dy --
dt
~- I x ( t ) - - z(t)'] 2
y(ty)=J.
(5)
(6)
In order to use the invariant imbedding concept for estimation, consider a varying interval of 0 ~
fo
Ix(t)-- z(t)]2 dt.
x(a)=c
(7)
(8)
(9)
because various values of c define various integral curves, i.e. various x(t). There is now a set.of equations consisting of eqs (1) and (6) with the final conditions eqs (8) and (9). From eqs (7)--(9), the invariant imbedding equation can be obtained as
dr da
dr(c, a) dc Or(c, a) Oc da F =[c-~(a)] 2
(10)
or
f(c,a)rc(c,a)+ro(c,a)=[c--z(a)] 2.
(11)
The job of filtering estimation is, among innumerable possible values, to find a value of c which minimizes y(a). If such a value of c is denoted by e(a): r~(c, a)lc=e=r~(e, a)=0.
r(c, a)=
dr~(e, a)--- rc~(e, a) de+r,o(e, a) d a = 0
(13)
de da -
r,o(e, a) r~(e, a)"
(14)
On the other hand, the derivative of eq. (11) with respect to c leads to
f (c, a)r¢~(c,a) +f~(c, a)r~(c, a) + rca(c, a) = 2[c - z(a)].
(15)
(t9)
where x(t) can be represented by
x(t)=c--
f(x, t) dt
(20)
t
in principle. The integral of the second term on the right side of eq. (20) is independent of c as c is arbitrarily given. Substituting eq. (20) into eq. (19) gives
r(c,a)--
f:[
;
c-z(t)-
f(x, t) dt
t
]2
dt.
(21)
For a given interval, the partial derivative of r with respect to c equals its total derivative. Thus
r~(c,a)lo=--=-- I dc dc Jo
c-z(t)-
dt
=fi2[c-z(t)-flf(x,t)dt]dt.
(22)
Similarly, the second-order partial derivative can be obtained as r c a
or
[x(t)-z(t)] 2 dt o
(12)
The total differential of eq. (12) is given by
(18)
which describes the change rate of the optimal-estimated value for x, e(a), with the terminal of the varying interval, a, where both functions f and z are known; and the function rc~ remains to be determined. Equation (18) is what Lee developed. However, he had difficulties in determining the function r~c. The following treatment is totally different from Lee's. In order to determine Function r~, it is essential to distinguish the variables. As mentioned above, r is a function of both a and c. However, c and a are quite different in nature, a is the terminal of the varying interval O<~t<~a, and is a one-dimensional variable, while c is an arbitrarily assumed value for x at a, and is a free variable in a two-dimensional space. For a given interval, a is constant, while c is a variable. From the definition:
then y is a function not only of a but also of c, which can be represented by
y(a) = r(c, a)
(16)
or
dr~
d
a
o
=2a.
(23)
Substituting eq. (23) into eq. (18) gives the onedimensi~onal filtering estimator equation as de
1
da =f(e, a) + a [ z ( a ) - el.
(24)
1667
Estimation of states and parameters by invariant imbedding technique INITIAL
Define
VALUE
If an initial time of ao > 0 is taken for an estimation, a point of x = ~ in the interval 0 <~t <~ao can necessarily be found so that
y(t)=
[xi(t)-zi(t)] 2 dr.
(33)
[xi(t)-- Zi(t)] 2
(34)
01=1
Thus
y(ao)= f
ao [x(t)--z(t)] 0
dy =
2 dt=ao[x(~)-z(¢)] 2. (25)
dt
For a small ao
y(ao) = ao[x(~ ) - z(~)] z _ ao[X(ao)- Z(ao)]z
n
i=1
y(t:)=J.
(35)
For the terminal of an interval, a, let x(a)=e
(36)
y(a) = r(c, a).
(37)
(26) then
from the continuity of the function. Since ao > 0 and [x(ao)-z(ao)]2>~O, the condition minimizing y(ao) must be
e(ao)= z(ao)
(27)
where ao should have a small positive value. Theoretically, ao could be zero, because the second term on the right side of eq. (24) has also a finite value at ao~O as [z(ao)-e(ao)]2~O, too, at the same time. However, a zero value can not be selected for ao from the numerical calculation point of view. The problem of one-dimensional estimation thus becomes the definite one consisting of eqs (24) and (27).
Now the requirement is to find a vector c which minimizes the integral
[x(ti)-z(Q] z dr.
y(a)=
(38)
0 i=1
If such a vector is represented by e(a):
[Or(c,a)l ~ey(a)jc=e=L~Jc=e=rc,(e, a)=0 (i= 1, 2 . . . . . n). (39) For a certain i, the total differential of eq. (39) can be obtained as
N-DIMENSIONAL
ESTIMATOR
EQUATIONS
The above treatment can be easily expanded to the N-dimensional problem. Consider a complex system with the following nonlinear kinetics vector equation: dx
--=f(x, t) dt
d r c -_ ~ '
~ T-rc,(e,a)dej+~arc,(e,a)da=O
(40)
j= 1 6 e j
and, for all i, the total differential can be represented by
roe(e, a) de + rat(e, a) da = 0
(28)
(41)
where --I
where x is the column vector consisting of N state variables of the system: x r = [xl x 2 . . . x,]
(29)
rcc =
and, naturally, f is also a N-dimensional vector: f r = I f , f 2 . . ,f,].
(30)
Assume that all the states can be measured, or that some of the states can not be measured directly, but the indirectly-measured values for them can be determined from the results directly measured for the others with certain relationships such as stoichiometry, kinetics and others. Thus, the measured values of all the states, z(t), can be represented as an N-dimensional vector equation: z(t) = x(t) + er(t)
(31)
where x is the true state vector of the system, and er is the measuring error vector. To estimate states with the measured values, the estimated values should minimize the integral
[x,(t)-z,(t)]2dt.
(32)
rcl£1
rc1¢2
" ' "
reich 1
r....
r....
...
r.... ]
rcnct
rc c2
. ••
re c
(42)
de r = [de 1 de e . . . de,]
(43)
r r = [rat1 rac2.., r,c,-I.
(44)
From eq. (41): de
- - = - [roe(e, a)] - ~rac(e, a) da
(45)
where [roe(e, a)]- 1 is the reverse of matrix rode, a). The invariant imbedding equation for the missing final condition r can be obtained as
ra(C, a ) + ~ fi(e, a)rc,(e, a ) = ~ [c,-z,(a)] 2. (46) i=1
i=1
The derivation ofeq. (46) with respect to vector e leads to rae(C, a) + [fc(e, a)] rre(c, a) + ree(e, a)f(e, a) = 2 [ e - z(a)] (47)
1668
W U YUAN
where
and the second-order derivative can be obtained as
I Lc, Lo2 ...
Lc.] rc~(C, a ) =
LS.c,
s.,
...
and r c is the column vector of the partial derivative of function r with respect to vector e: r [ = [re, re2--, re,].
(49)
Substituting c by e and using the extreme value conditions, rc(e, a ) = 0, eq. (47) becomes
I
de
(50)
(51)
Now the function matrix [rcc(e, a)] - 1 remains to be determined. From the definition:
[x~(t)-zi(t)] 2 dt
(52)
0i
where x~ can be represented by
t)dt 0 = 1 , 2 . . . . n)
(53)
in principle. Substituting eq. (53) into eq. (52) leads to
f~(x, t)
cl-zi(t)Oi=l
dt dt.(54)
ro,(e,a)=f i 2[ c,- z,(t)- f lf~(x,t) dt ] dt.
(55)
The derivative of function r with respect to vector c can be written as
re(e,a)=
0 2 c2-z2(t)-
0
...
2a
(57)
1
(58)
(59)
F o r estimations, no essential difference exists between a state and a parameter. The N-dimensional estimator equations can be used for the estimation of both states and parameters. The estimator equations developed in this paper are simple. To estimate N states and/or parameters, the total number of differential equations to be solved is only N; and each of the equations is simple, of definite initial value, and convenient for calculation. Since no approximations were made in the development, the estimator equations can give very accurate results• By contrast with Lee's (1968b) method of estimating N states and/or parameters, the number of equations to be simultaneously solved is equal to (N + N2).
t
As discussed for the one-dimensional problem, care must be taken to distinguish between the variables c~ and a. For a given interval, a is constant, while c~ are free variables. Since each c~ is arbitrarily given, the integrals of the second term on the fight side of eq. (53) are independent of any c~, and all e~are independent of each other. Thus, for the terminal of a given interval, da = 0, and the partial derivative of r with respect to a certain c~ is equal to its total derivative. From eq. (54):
o 2 cl-zl(t)-
0
1
e(ao)= Z(ao).
de - - = f(e, a) + 2[rcc(e, a)] - ~[z(a)-e]. da
;[ ;[ ;[
00
and its initial condition can be easily obtained as
Combining it with eq. (45) leads to
r(¢, a)=
.. .. ..
~ = f ( e , a) + -a [z(a) - e]
+2[rcc(e, a)]-l[z(a)--e].
xi(t)=ci+
02a
Substituting the reverse of the matrix of eq. (57) into eq. (51) and rearranging the expression gives the Ndimensional filtering estimator equation as
-[rcc(e, a)]-lrac(e, a)= f(e, a)
r(e, a)=
02a
; ; ;
1 ] }
fl(x,t) dt dt f2(x,t)dt
A P P L I C A T I O N S IN KINETICS RESEARCH
In order to illustrate the method of estimating parameters with the N-dimensional filtering estimator equations in kinetics research, two cases are given below.
Case 1 Hellin and Jungers [see Levenspiel (1972, p. 88)] studied the kinetics of the reaction of sulphuric acid with diethyl sulphuric acid: (C2Hs)2SO, + H2SO 4 = 2C2HsSO,H.
(A)
(B)
(60)
(R)
The initial concentrations of the components in the mixture were CA0=CB0=5.5 tool/l,
Cao=0.
The reaction occurred at 22.9°C and the measured data are listed in Table 1. If the system has secondorder reversible kinetics:
dt
dC~ = kx CACn-- k2C~"
o2 c.--z.(t)-- ,f"(x't) dt dt
dt
(56)
(61)
Determine the parameters with the estimator equations.
1669
Estimation of states and parameters by invariant imbedding technique
e2(ao)= Zz(ao)
Table 1. Measured data for the concentration of the product, CR
t (rain)
CR (tool/l)
t (rain)
CR (mol/l)
0 41 48 55 75 96 127 146 162
0 1.18 1.38 1.63 2.24 2.75 3.31 3.76 3.81
180 194 212 267 318 368 379 410 (oo)
4.11 4.31 4.45 4.86 5.15 5.32 5.35 5.42 (5.80)
= zl(a°+At)-zl(a°) [zl(a°)]2 ~ - ~ 4.9763 J
ao=41 min.
el(ao) = e1(41) = 1.18 mol/1 e2(a°)=e2(41)
5"82 = 4.9763. (5.5-5.0/2) 2
Substituting eq. (62) into eq. (61) and using the value for K, eq. (61) can be rewritten as 4.97631 C2 ]
(63)
1.38--1.18 48 --41
=0.
(64)
Thus, the estimator equations for the system can be written as follows: ~
1
= e2[ (5-5 - e,/2)2 -4.9763
1
e~ +~[zl(a)-el] (65)
de2
1
da
a
--=-[z2(a)-e2]
5.5--
2
4.9763 = 0.001199 1/(mol min).
(73)
There is now a definite problem consisting of eqs (65), (66), (68), (71), (72) and (73). If a step size of 1 min is chosen, then
ai+ 1 =ai.-b 1 and the "measured" value z~ at each ai at which no actual measurement has been made, zt(ai), can be determined simply by a linear interpolation procedure from the actually measured values. The calculated results are shown in Fig. 1. It can be seen that, after t = 250 min, the estimated value for CR, el, converges well to that measured; and both the convergence and the stability of the value estimated for the rate constant k ~, e2, are very good. The value finally calculated for kl is
(66)
where z~ is the measured value for CR, for which data are listed in Table 1. There is no value directly measured for the rate constant k~, while using the kinetics relation, eq. (63), the indirectly-measured value for it, z2, can be represented in term of z~ as
E( 1,8)2
(72)
1.182 1-1
and the kinetics of parameter k~ can be described by
dt
(71)
This is reasonable since 41 min is short enough in comparison with the whole measuring time, 410 min. Thus, from the data listed in Table 1:
F r o m the measured data at 0% the reaction equilibrium constant can be determined as
dkt
(70)
As mentioned above, the initial time ao should be a small positive value. However, it is more convenient for this case to take the time of 41 min, at which the first measurement was made, as the initial one, that is
Ca = Ca = (Cao + CBo- CR)/2 = 5.5 - CR/2. (62)
d(C5'"=5k- lc[' / 2 ) 2 - d t
2
At
Solution: a material balance yields
K = [C*]2 [C*]2 CA*Ca* (5.5-C*/2) 2
{E 5 . 5 -
k 1 =0.001333 1/(mol min). The rate constant of the reverse reaction, k2, can be found according to the relationship between k 1, k 2 and K, as
k 2 = kl/K = 0.001333/4.9763 = 2.799 × 10 -4) 1/(tool min).
4.9763 z
.
(67)
However, for a numerical calculation, using the optimal-estimated value e~ to represent z2, that is, using the relation z2(a)=
(5.5-et/2) 2-
1 e~ 4.9763
l-'
o~o.o~o
_
•
(68)
• --
= estimated
el(ao)=zl(ao)
(69)
0
k~
of
×
~"
is more convenient and may improve the convergence. The initial values for el and e2 can be determined from eqs (59) and (68), as
1::
E
/
-
1O0
200 a =
300
400
t (min)
Fig. 1. Estimating and calculated results (case 1).
e
1670
Wu YUAN
Thus, the kinetics of the system studied can be expressed by dCR
dt
Table 2. Measured data for the concentration of A, Ca
= 1.333 x 10-3C~CB-2.799 x 10-4C~. (74)
The simulation calculated results by eq. (74) are shown in Fig. 2. It can be seen that the equation fits the measured data well. The standard error of the calculated results to the experimental data is found as s=0.1080mol/l, showing that the estimated parameters are appropriate. It is interesting to note that, in the estimation above, no particular treatment of the experimental data was done before putting them into the computer and the "measured" values at the points where no actual measurements were made were determined simply by a linear interpolation. This is very convenient. For the same case, Levenspiel (1975) used the integral expression of eq. (61):
In[ x* -(2x*-l x*-xa )xal=_]2k,(-~a-l)Caot
dC--~g=6.742 x IO-4CaCB - 1.355 x 10-4C~. (76) dt The standard error of the results calculated by eq. (76) to the measured data is sL= 1.164 mol/1. Obviously, estimating the parameters with the N-dimensional estimator equations on a computer is more convenient and can result in higher accuracy. Case 2 Sucrose was catalyzed with the enzyme sucrase to hydrolyze at room temperature (Levenspiel, 1972): (E)
sucrose
, product.
CA (mmol~)
t (h)
CA (mmol~)
0 1 2 3 4 5
1.0 0.84 0.68 0.53 0.38 0.27
6 7 8 9 10 11
0.16 0.09 0.04 0.018 0.006 0.0025
the Michaelis-Menten equation dCA
(77)
The initial concentration of sucrose is Cao = 1.0 mmol/l, and that of sucrase is Ceo =0.01 mmol/1. The data listed in Table 2 were obtained with a batch reactor. The kinetics of the system can be described by
k3CEoC.t
dt
(78)
CA+M
Determine the parameters k 3 and M. Solution: for convenience, eq. (78) can be rewritten as
(75)
to determine the parameters k~ and k 2 by a graphic method, and obtained the following rate expression:
sucrose (A)
t (h)
dCA
-- k4Ca
dt
1 + ks CA
(79)
where k4=
k3CEo M
and
1 ks=~.
Since there is one set of data measured for state C a and one restrictive relationship (a total of two only), it is impossible to estimate state CA and parameters k4 and ks (a total of three) directly and simultaneously, and a trial-and-error procedure is required. There are two possible ways: (1) trying k4, and (2) trying ks. Which one should be selected could be determined from the sensitivity of the kinetics model to it. F o r the case discussed here, the model is more sensitive to k4 than to ks, or, in other words, the effect of the change in k4 on the reaction rate is larger than that of the change in k s . Therefore, it may be expected to get more accurate results if trying k4, that is, estimating Ca and k5 with various assumed values for k4. Let CA = x l and k 5 = x 2, then the kinetics equation can be rewritten as the following set: dx 1
k4x 1
.
dt
(81)
1 + x2x 1
dx 2
--=0
6
and the corresponding estimator equations are j 0 *0" -
4 0
EVa
cY
2 1 0
(82)
dt
5
(80)
/
/°
/°
de1
da
CR calculated by eq. (74)
k4et
1 I-- [ z l ( a ) - el] l+e2e 1 a
de 2
1 -
da 1 100
I
I
21)11
300
I 4OO
(84)
[z2(a)-e2].
The indirectly-measured values for ks, zz(a), can be represented by
t (mini
Fig. 2. Results simulated by eq. (74).
a
(83)
z2(a)=
-k,(del~-I \ da ]
1 el
(85)
Estimation of states and parameters by invariant imbedding technique = k4/0.16 - - 1/0.84.
from the kinetics equation, eq. (81). If the time of 1 h after the reaction begins is taken as the initial time, that is, if ao= 1 h
then the initial values for e 1 and e z can be obtained from the data listed in Table 2 as (87)
. [zl(ao+At)--za(ao)-] -x 1 e2(a")=zz(a°)= --k4" L ~ - z,(ao)
J
(a)
(88)
The results calculated with a step size of 0.1 h and ming various values of k4 are shown in Fig. 3. It seems that the convergence of the estimated values for CA, el, with various values of k4 assumed in a considerably large range are all very good; while the stabilities of those for k 5 are different. There is a visual difference even between the two results obtained with assuming k 5 =0.9994 and 0.9997. This suggests that the way of trying k4 is sensitive enough and consequently can result in good accuracy. The best results
(86)
e~(ao) = 0.84 mmol/l
k 4 = 0.999
k s = 5.054377 10
\\
0.I
¢o
0.6
8 e2
/
6
vE
E v
o,4
4
Zl
~
81
Z
0.2
n ] I 6
- - I l l L L I
I
I
L I 4
J n J I 8
J I
0
I I 10
a=t(h)
(b)
k 4 = 0.9994
k 5 = 5.061156 -10
1m
-8
0.8
0
E E ---= L~
"\
0.6
/
e2
E 0.4
zl/°~
°
je
II
1
0..~
I I I J l
--I
0
4
2
6
8
10
a = t (h)
t
0.8
O.E E E == (.J
(C)
k 5 = 5.071345
k4 = 1
-10
J
-\ "\ -
82
/
-
\
O.Z
~
el
0.2
'e ~o~e,....," 0
"6 E E II
zl--°~o
JlLll 2
1671
4
6 a =t(h)
Fig. 3. (a)-(c).
8
10
1672
Wu (d)
YUAN
k 4 == 1.1
ks =
7.124979
1 o.~
3 e2
E E
0.6
E ¢
It
0.2
L
I
J
I
,
I
,
I
2
J
I
,
I
4
a
I
J
6
I
t
I
~
8
n"T
10
a=t(h) (e)
1(4 =
0.9997
ks =
5.066253
10
1
"\\
0.8
-6 E E
O.E
=
r
8
/ e2
''0~
0.~
"\
Zl
6
-6 E vE
4 el
0.; 0 ~- I ) ' [ ' 0
I
,
2
l 4
J
I
~
1
6
,
l
~
-
I
,
[
i
~
£
I
10
a = t (h)
(f)
k4 =
k5 =
0.9
3.626837
-d10
1
0.8
\
0.6
E
8
"\
?-
"\
/
e2
,j= 0.4--
-6 E E
4
~II
2
0.2 ~
o- I l l l l o 2
-- 6
~0 4
6
10
8
a=t(h)
Fig. 3. Estimated and calculated results with various assumed values for ka.
are those with k4=0.9994, and the corresponding value estimated for ks is e2=5.0612 [see Fig. 3(b)]. More accurate results might be possible by further tries, but that is no longer necessary for a kinetics model. Therefore, the parameters are determined as
k+M =,0.9994,,0.19758,z/t~ 19.746h-1. k3 = CEo 0.01 -The kinetics of the system can thus be expressed as dCa
19.746CroCA -
ka=0.9994h -1,
k5=5.0612 l/mmol.
Correspondingly, from eq. (80): M = 1/k 5 = 1/5.0612=0.19758 mmol/l
dt
(89)
0.19758 + C A
The results simulated by eq. (89) are given in Fig. 4, and the standard error of the measured data is s = 0.005011 mmol/l.
Estimation of states and parameters by invariant imbedding technique
0.8
•
3E o.e E ¢~
e~e'=
Calculatedbyeq.(89)
0.4
0.;
ieasu red~e~e~e~e~e,..
0.0-
j
I 2
I
I 4
I
I 6 t(h)
I
~ I 8
"~'~'~-10
11
distinguishing models, the only requirement is one set of data on changes in states such as composition, temperature, pressure etc., with time or an equivalent independent variable. The experiments may be carried out simply with a batch reactor, provided that the measuring methods used are fast and accurate enough. Therefore, the estimator equations are of potential in simplifying experimental procedures and reducing the number of necessary experiments. In some cases, not all the states and the parameters can be directly and simultaneously estimated. In general, the criterion below should be satisfied:
Fig. 4. Results simulated by eq. (89)
Ne<,Nz+N r
The parameters k4 and k 5 can also be determined by trying k5, and the estimated values in such a way are k4=0.9686 h - 1 and k5 = 5.035 1/mmol, respectively. The standard error of the results calculated by the kinetics equation with these values of parameters is found to be s'-0.01362mmol/l, showing that the accuracy is lower than that by trying k4. By a graphical procedure, Levenspiel (1975) obtained the following kinetics equation for the same case: dCA d-T-=
16.4CEoCa 0.164 + C a
(90)
A simulation of eq. (90) leads to a standard error of sL = 0.05365 mmol/l. The accuracy is significantly lower than those obtained with the estimator equations for both ways of trying k4 and trying k5 and, in fact, a graphical procedure is usually troublesome. DISCUSSION In all of the cases discussed above, the kinetics equations selected are of the proper forms. The results calculated for these cases show that the N-dimensional estimator equations developed with the invariant imbedding concept can be used for the determination of parameters in kinetics research and, for a kinetics model of proper form, can result in a good accuracy, although the equations themselves are formally simple. To determine whether a kinetic model is proper or not, a procedure for distinguishing models is generally required. The estimator equations can also be used for such a procedure, although that has not been discussed in this paper. For each of the models to be distinguished, the principle and the method of the estimation calculation are the same as for the others. For complex systems, the estimator equations may present advantages over other methods currently used in kinetics research. Since all of the estimation calculation can be carried out sequentially on a computer, the estimator equations are more convenient for practical applications, and may eliminate some difficulties in experiments and interpreting data. The estimator equations require fewer and simpler data. For example, when including the procedure of
1673
(91)
where Ne is the total number of states and parameters to be estimated, N~ the number of states measured, and N, the number of independent restrictive relationships. If N e > ( N z + Nr), a trial-and-error procedure is needed. Which parameters should be selected for trying can be determined according to the sensitivities of the kinetics model to them in order to obtain more accurate results. It appears that a trial-and-error procedure requires a larger amount of calculation. However, the estimation calculation on a computer is rather convenient, even if the trial-and-error method has to be used. The estimator equations themselves do not require that the independent variable must be time and do not specify what variables the x~ must be. It can be expected that the N-dimensional estimator equations may find more and more applications in chemical engineering and other fields by appropriately setting the variables.
a C, C e, e
J r t X, X
Y Z, Z
NOTATION final value of independent variable t final value of state x or state vector x optimal-estimated value for state x or state vector x function to be minimized = y(a), missing final condition independent variable dependent variable, state or state vector of the system integral function defined by eq. (7) or (33) value measured for state x or state vector x
REFERENCES
Bellman, R. and Lee, E. S., 1977, Invariant imbedding and solid-fluid reaction--I. The unreacted core model• Chem. Engn9 Sci. 32, 1051-1054. Bullin, J. A. and Dukler, A. E., 1975, Hybrid computer modelling of a stochastic nonlinear dynamic system. Chem. Engn 9 Sci. 30, 631-35. Diamond, W. J., 1981, Practical Experimental Desion, pp. 76-83. Wadsworth, Belmont, CA. Lee, E. S., 1968a, Invariant imbedding, a versatile computation concept. Ind. Engng Chem. 60(9), 55-64. Lee, E. S., 1968b, Quasilinearization and the estimation of parameters in differential equations. Ind. Enon9 Chem. Fundam. 7, 152-158. Lee, E. S., 1968c, Invariant imbedding nonlinear filtering, and
1674
WU YUAN
parameter estimation. Ind. Engng Chem. Fundam. 7, 164-171. Levenspiel, O., 1972, Chemical Reaction Engineering, 2nd edition, pp. 88-89. John Wiley, New York. Levenspiel, O., 1975, Solutions to Problems in Chemical Reaction Engineering, 2nd edition. Oregon State University.
Omatu, J. A., 1986, Optimal estimation theory for distributed parameter systems, Control Dynamic Systems 23, 195-241. Rose, L. M., 1981, Chemical Reactor Design in Practice, pp. 127-137. Elsevier, Oxford. Stephanopoulos, G., 1984, Chemical Process Control, an Introduction to Theory and Practice, pp. 113-115. PrenticeHall, Englewood Cliffs, NJ.
0009 2509/89 $3.00+0.00 © 1989 Pergamon Press plc
Printed in Great Britain.
ESTIMATION OF STATES A N D PARAMETERS BY INVARIANT IMBEDDING TECHNIQUE WU YUAN Department of Chemical Engineering, Wuhan Institute of Chemical Engineering, Wuhan, China (Received 24 October 1987; accepted for publication 28 November 1988) Abstract--A set of sequential estimator equations has been developed using the invariant imbedding
concept. The equations can be used for both the filteringestimation of states of non-linear dynamic systems and the estimation of parameters in kinetics research and other areas, and are of potential in simplifying experimental procedures and reducing the number of necessary experiments for the latter purpose. Results calculated for severalcases show that the convergenceof the estimated values for states, the stability of those for parameters, and the accuracy are all very good, suggestingthat the estimator equations may be expected to be a useful tool in chemical engineering and even in many other fields.
INTRODUCTION Kinetics research and process control are important in many areas. To solve them requires the estimation of states and parameters. Establishing a model is one of the main jobs in kinetics research and is an important basis for process analysis, design, scale-up, operation and control. However, the experimental investigation for this purpose is usually complex, difficult, expensive and timeconsuming. Therefore, it is useful to develop new methods for data interpretation which require less data and may allow a simplified experimental procedure. Since the 1970s, great attention has been devoted to "sequential experiment planning" [e.g. Rose (1981) and Diamond (1981)], for which the basic idea is fully using the information provided by the last experiments to appropriately plan the next ones in order to enable each experiment to give a high responsibility and thus a high efficiency. However; the experimental procedure required by this method is also considerably complex. The correct estimation of states is the prerequisite for effective process control. Therefore, the filtering estimation of states of dynamic systems is continuously one of the main research projects in process control. It is general for non-linear systems to "linearize" them, as the linear differential equations may easily be solved with the Laplace transforms (Stephanopoulos, 1984). However, the linearization usually introduces certain errors into the system, and, in some cases, the errors introduced may exceed the acceptable range. For the estimation of the states of non-linear systems themselves, a number of methods have been developed (Omatu, 1986), in which the real physical systems are usually described with stochastic models, that is, each of the equations describing dynamic behavior of the states includes a white noise excitation. Such a treatment may yield a more accurate analysis, but it also creates difficulties in computation and concept as well (Bullin and Dukler, 1975). On the other hand, a definite (non-stochastic)
model may be accurate enough for a number of processes in chemical engineering. Therefore, it is also of interest to study the filtering estimation of the states of non-linear dynamic systems with non-stochastic models. The invariant imbedding technique is one of the effective methods for solving boundary value problems of non-linear differential equations, for which certain applications in chemical engineering have been developed (Lee, 1968a; Bellman and Lee, 1977). Lee (1968b, c) used the invariant imbedding concept for non-linear filtering and parameter estimation. Unfortunately, because of omissions in the mathematical treatment, the estimator equations Lee developed are very complex, with no definite initial values, and thus have little value for practical applications. In fact, no successful application or progress in improving Lee's equations has been reported. In the present paper, a set of new sequential filtering estimator equations will be developed, and two applications to kinetics research will be given. ONE-DIMENSIONALESTIMATOREQUATION Consider a system with a single state variable whose kinetic behavior can be described by the non-linear equation dx
- - = f ( x , t). dt
(1)
From the initial time, t =0, to the current instant, t = ty, the state x is measured, but random errors must exist. The value measured for state x at time t, z(t), can be represented by z(t) = x(t) + er(t)
(2)
where x is the true state, and er is the measuring error. The requirement is to estimate the true state x(t) using the measurement within the period from 0 to tf. Provided that the state x at a certain instant is exactly estimated, then the definite condition of eq. (1) is
1665
1666
W U YUAN
determined, and thus the value for x at any time can be easily calculated. If one first estimates the current state x(tl), then the problem becomes the filtering estimation. From the classical least-squares criterion, the estimated value for x(tl) should minimize the integral
J=f
[x(t)-z(t)]2dt.
(3)
Substituting c by e and using eq. (12), eq. (15) becomes
f(e, a)rcc(e, a) + rca(e, a) = 2 [ e - z(a)] rca(e, a) ----=f(e, rcc(e, a)
a) -~
2[z(a)- e]
(17)
roe(e, a)
Combining eq. (17) with eq. (14):
de 2[z(a)-e] ~a=f(e,a)+ rcc(e,a~
Defining y(t) = I t [ x ( t ) - z(t)] 2 dt 30
(4)
gives dy --
dt
~- I x ( t ) - - z(t)'] 2
y(ty)=J.
(5)
(6)
In order to use the invariant imbedding concept for estimation, consider a varying interval of 0 ~
fo
Ix(t)-- z(t)]2 dt.
x(a)=c
(7)
(8)
(9)
because various values of c define various integral curves, i.e. various x(t). There is now a set.of equations consisting of eqs (1) and (6) with the final conditions eqs (8) and (9). From eqs (7)--(9), the invariant imbedding equation can be obtained as
dr da
dr(c, a) dc Or(c, a) Oc da F =[c-~(a)] 2
(10)
or
f(c,a)rc(c,a)+ro(c,a)=[c--z(a)] 2.
(11)
The job of filtering estimation is, among innumerable possible values, to find a value of c which minimizes y(a). If such a value of c is denoted by e(a): r~(c, a)lc=e=r~(e, a)=0.
r(c, a)=
dr~(e, a)--- rc~(e, a) de+r,o(e, a) d a = 0
(13)
de da -
r,o(e, a) r~(e, a)"
(14)
On the other hand, the derivative of eq. (11) with respect to c leads to
f (c, a)r¢~(c,a) +f~(c, a)r~(c, a) + rca(c, a) = 2[c - z(a)].
(15)
(t9)
where x(t) can be represented by
x(t)=c--
f(x, t) dt
(20)
t
in principle. The integral of the second term on the right side of eq. (20) is independent of c as c is arbitrarily given. Substituting eq. (20) into eq. (19) gives
r(c,a)--
f:[
;
c-z(t)-
f(x, t) dt
t
]2
dt.
(21)
For a given interval, the partial derivative of r with respect to c equals its total derivative. Thus
r~(c,a)lo=--=-- I dc dc Jo
c-z(t)-
dt
=fi2[c-z(t)-flf(x,t)dt]dt.
(22)
Similarly, the second-order partial derivative can be obtained as r c a
or
[x(t)-z(t)] 2 dt o
(12)
The total differential of eq. (12) is given by
(18)
which describes the change rate of the optimal-estimated value for x, e(a), with the terminal of the varying interval, a, where both functions f and z are known; and the function rc~ remains to be determined. Equation (18) is what Lee developed. However, he had difficulties in determining the function r~c. The following treatment is totally different from Lee's. In order to determine Function r~, it is essential to distinguish the variables. As mentioned above, r is a function of both a and c. However, c and a are quite different in nature, a is the terminal of the varying interval O<~t<~a, and is a one-dimensional variable, while c is an arbitrarily assumed value for x at a, and is a free variable in a two-dimensional space. For a given interval, a is constant, while c is a variable. From the definition:
then y is a function not only of a but also of c, which can be represented by
y(a) = r(c, a)
(16)
or
dr~
d
a
o
=2a.
(23)
Substituting eq. (23) into eq. (18) gives the onedimensi~onal filtering estimator equation as de
1
da =f(e, a) + a [ z ( a ) - el.
(24)
1667
Estimation of states and parameters by invariant imbedding technique INITIAL
Define
VALUE
If an initial time of ao > 0 is taken for an estimation, a point of x = ~ in the interval 0 <~t <~ao can necessarily be found so that
y(t)=
[xi(t)-zi(t)] 2 dr.
(33)
[xi(t)-- Zi(t)] 2
(34)
01=1
Thus
y(ao)= f
ao [x(t)--z(t)] 0
dy =
2 dt=ao[x(~)-z(¢)] 2. (25)
dt
For a small ao
y(ao) = ao[x(~ ) - z(~)] z _ ao[X(ao)- Z(ao)]z
n
i=1
y(t:)=J.
(35)
For the terminal of an interval, a, let x(a)=e
(36)
y(a) = r(c, a).
(37)
(26) then
from the continuity of the function. Since ao > 0 and [x(ao)-z(ao)]2>~O, the condition minimizing y(ao) must be
e(ao)= z(ao)
(27)
where ao should have a small positive value. Theoretically, ao could be zero, because the second term on the right side of eq. (24) has also a finite value at ao~O as [z(ao)-e(ao)]2~O, too, at the same time. However, a zero value can not be selected for ao from the numerical calculation point of view. The problem of one-dimensional estimation thus becomes the definite one consisting of eqs (24) and (27).
Now the requirement is to find a vector c which minimizes the integral
[x(ti)-z(Q] z dr.
y(a)=
(38)
0 i=1
If such a vector is represented by e(a):
[Or(c,a)l ~ey(a)jc=e=L~Jc=e=rc,(e, a)=0 (i= 1, 2 . . . . . n). (39) For a certain i, the total differential of eq. (39) can be obtained as
N-DIMENSIONAL
ESTIMATOR
EQUATIONS
The above treatment can be easily expanded to the N-dimensional problem. Consider a complex system with the following nonlinear kinetics vector equation: dx
--=f(x, t) dt
d r c -_ ~ '
~ T-rc,(e,a)dej+~arc,(e,a)da=O
(40)
j= 1 6 e j
and, for all i, the total differential can be represented by
roe(e, a) de + rat(e, a) da = 0
(28)
(41)
where --I
where x is the column vector consisting of N state variables of the system: x r = [xl x 2 . . . x,]
(29)
rcc =
and, naturally, f is also a N-dimensional vector: f r = I f , f 2 . . ,f,].
(30)
Assume that all the states can be measured, or that some of the states can not be measured directly, but the indirectly-measured values for them can be determined from the results directly measured for the others with certain relationships such as stoichiometry, kinetics and others. Thus, the measured values of all the states, z(t), can be represented as an N-dimensional vector equation: z(t) = x(t) + er(t)
(31)
where x is the true state vector of the system, and er is the measuring error vector. To estimate states with the measured values, the estimated values should minimize the integral
[x,(t)-z,(t)]2dt.
(32)
rcl£1
rc1¢2
" ' "
reich 1
r....
r....
...
r.... ]
rcnct
rc c2
. ••
re c
(42)
de r = [de 1 de e . . . de,]
(43)
r r = [rat1 rac2.., r,c,-I.
(44)
From eq. (41): de
- - = - [roe(e, a)] - ~rac(e, a) da
(45)
where [roe(e, a)]- 1 is the reverse of matrix rode, a). The invariant imbedding equation for the missing final condition r can be obtained as
ra(C, a ) + ~ fi(e, a)rc,(e, a ) = ~ [c,-z,(a)] 2. (46) i=1
i=1
The derivation ofeq. (46) with respect to vector e leads to rae(C, a) + [fc(e, a)] rre(c, a) + ree(e, a)f(e, a) = 2 [ e - z(a)] (47)
1668
W U YUAN
where
and the second-order derivative can be obtained as
I Lc, Lo2 ...
Lc.] rc~(C, a ) =
LS.c,
s.,
...
and r c is the column vector of the partial derivative of function r with respect to vector e: r [ = [re, re2--, re,].
(49)
Substituting c by e and using the extreme value conditions, rc(e, a ) = 0, eq. (47) becomes
I
de
(50)
(51)
Now the function matrix [rcc(e, a)] - 1 remains to be determined. From the definition:
[x~(t)-zi(t)] 2 dt
(52)
0i
where x~ can be represented by
t)dt 0 = 1 , 2 . . . . n)
(53)
in principle. Substituting eq. (53) into eq. (52) leads to
f~(x, t)
cl-zi(t)Oi=l
dt dt.(54)
ro,(e,a)=f i 2[ c,- z,(t)- f lf~(x,t) dt ] dt.
(55)
The derivative of function r with respect to vector c can be written as
re(e,a)=
0 2 c2-z2(t)-
0
...
2a
(57)
1
(58)
(59)
F o r estimations, no essential difference exists between a state and a parameter. The N-dimensional estimator equations can be used for the estimation of both states and parameters. The estimator equations developed in this paper are simple. To estimate N states and/or parameters, the total number of differential equations to be solved is only N; and each of the equations is simple, of definite initial value, and convenient for calculation. Since no approximations were made in the development, the estimator equations can give very accurate results• By contrast with Lee's (1968b) method of estimating N states and/or parameters, the number of equations to be simultaneously solved is equal to (N + N2).
t
As discussed for the one-dimensional problem, care must be taken to distinguish between the variables c~ and a. For a given interval, a is constant, while c~ are free variables. Since each c~ is arbitrarily given, the integrals of the second term on the fight side of eq. (53) are independent of any c~, and all e~are independent of each other. Thus, for the terminal of a given interval, da = 0, and the partial derivative of r with respect to a certain c~ is equal to its total derivative. From eq. (54):
o 2 cl-zl(t)-
0
1
e(ao)= Z(ao).
de - - = f(e, a) + 2[rcc(e, a)] - ~[z(a)-e]. da
;[ ;[ ;[
00
and its initial condition can be easily obtained as
Combining it with eq. (45) leads to
r(¢, a)=
.. .. ..
~ = f ( e , a) + -a [z(a) - e]
+2[rcc(e, a)]-l[z(a)--e].
xi(t)=ci+
02a
Substituting the reverse of the matrix of eq. (57) into eq. (51) and rearranging the expression gives the Ndimensional filtering estimator equation as
-[rcc(e, a)]-lrac(e, a)= f(e, a)
r(e, a)=
02a
; ; ;
1 ] }
fl(x,t) dt dt f2(x,t)dt
A P P L I C A T I O N S IN KINETICS RESEARCH
In order to illustrate the method of estimating parameters with the N-dimensional filtering estimator equations in kinetics research, two cases are given below.
Case 1 Hellin and Jungers [see Levenspiel (1972, p. 88)] studied the kinetics of the reaction of sulphuric acid with diethyl sulphuric acid: (C2Hs)2SO, + H2SO 4 = 2C2HsSO,H.
(A)
(B)
(60)
(R)
The initial concentrations of the components in the mixture were CA0=CB0=5.5 tool/l,
Cao=0.
The reaction occurred at 22.9°C and the measured data are listed in Table 1. If the system has secondorder reversible kinetics:
dt
dC~ = kx CACn-- k2C~"
o2 c.--z.(t)-- ,f"(x't) dt dt
dt
(56)
(61)
Determine the parameters with the estimator equations.
1669
Estimation of states and parameters by invariant imbedding technique
e2(ao)= Zz(ao)
Table 1. Measured data for the concentration of the product, CR
t (rain)
CR (tool/l)
t (rain)
CR (mol/l)
0 41 48 55 75 96 127 146 162
0 1.18 1.38 1.63 2.24 2.75 3.31 3.76 3.81
180 194 212 267 318 368 379 410 (oo)
4.11 4.31 4.45 4.86 5.15 5.32 5.35 5.42 (5.80)
= zl(a°+At)-zl(a°) [zl(a°)]2 ~ - ~ 4.9763 J
ao=41 min.
el(ao) = e1(41) = 1.18 mol/1 e2(a°)=e2(41)
5"82 = 4.9763. (5.5-5.0/2) 2
Substituting eq. (62) into eq. (61) and using the value for K, eq. (61) can be rewritten as 4.97631 C2 ]
(63)
1.38--1.18 48 --41
=0.
(64)
Thus, the estimator equations for the system can be written as follows: ~
1
= e2[ (5-5 - e,/2)2 -4.9763
1
e~ +~[zl(a)-el] (65)
de2
1
da
a
--=-[z2(a)-e2]
5.5--
2
4.9763 = 0.001199 1/(mol min).
(73)
There is now a definite problem consisting of eqs (65), (66), (68), (71), (72) and (73). If a step size of 1 min is chosen, then
ai+ 1 =ai.-b 1 and the "measured" value z~ at each ai at which no actual measurement has been made, zt(ai), can be determined simply by a linear interpolation procedure from the actually measured values. The calculated results are shown in Fig. 1. It can be seen that, after t = 250 min, the estimated value for CR, el, converges well to that measured; and both the convergence and the stability of the value estimated for the rate constant k ~, e2, are very good. The value finally calculated for kl is
(66)
where z~ is the measured value for CR, for which data are listed in Table 1. There is no value directly measured for the rate constant k~, while using the kinetics relation, eq. (63), the indirectly-measured value for it, z2, can be represented in term of z~ as
E( 1,8)2
(72)
1.182 1-1
and the kinetics of parameter k~ can be described by
dt
(71)
This is reasonable since 41 min is short enough in comparison with the whole measuring time, 410 min. Thus, from the data listed in Table 1:
F r o m the measured data at 0% the reaction equilibrium constant can be determined as
dkt
(70)
As mentioned above, the initial time ao should be a small positive value. However, it is more convenient for this case to take the time of 41 min, at which the first measurement was made, as the initial one, that is
Ca = Ca = (Cao + CBo- CR)/2 = 5.5 - CR/2. (62)
d(C5'"=5k- lc[' / 2 ) 2 - d t
2
At
Solution: a material balance yields
K = [C*]2 [C*]2 CA*Ca* (5.5-C*/2) 2
{E 5 . 5 -
k 1 =0.001333 1/(mol min). The rate constant of the reverse reaction, k2, can be found according to the relationship between k 1, k 2 and K, as
k 2 = kl/K = 0.001333/4.9763 = 2.799 × 10 -4) 1/(tool min).
4.9763 z
.
(67)
However, for a numerical calculation, using the optimal-estimated value e~ to represent z2, that is, using the relation z2(a)=
(5.5-et/2) 2-
1 e~ 4.9763
l-'
o~o.o~o
_
•
(68)
• --
= estimated
el(ao)=zl(ao)
(69)
0
k~
of
×
~"
is more convenient and may improve the convergence. The initial values for el and e2 can be determined from eqs (59) and (68), as
1::
E
/
-
1O0
200 a =
300
400
t (min)
Fig. 1. Estimating and calculated results (case 1).
e
1670
Wu YUAN
Thus, the kinetics of the system studied can be expressed by dCR
dt
Table 2. Measured data for the concentration of A, Ca
= 1.333 x 10-3C~CB-2.799 x 10-4C~. (74)
The simulation calculated results by eq. (74) are shown in Fig. 2. It can be seen that the equation fits the measured data well. The standard error of the calculated results to the experimental data is found as s=0.1080mol/l, showing that the estimated parameters are appropriate. It is interesting to note that, in the estimation above, no particular treatment of the experimental data was done before putting them into the computer and the "measured" values at the points where no actual measurements were made were determined simply by a linear interpolation. This is very convenient. For the same case, Levenspiel (1975) used the integral expression of eq. (61):
In[ x* -(2x*-l x*-xa )xal=_]2k,(-~a-l)Caot
dC--~g=6.742 x IO-4CaCB - 1.355 x 10-4C~. (76) dt The standard error of the results calculated by eq. (76) to the measured data is sL= 1.164 mol/1. Obviously, estimating the parameters with the N-dimensional estimator equations on a computer is more convenient and can result in higher accuracy. Case 2 Sucrose was catalyzed with the enzyme sucrase to hydrolyze at room temperature (Levenspiel, 1972): (E)
sucrose
, product.
CA (mmol~)
t (h)
CA (mmol~)
0 1 2 3 4 5
1.0 0.84 0.68 0.53 0.38 0.27
6 7 8 9 10 11
0.16 0.09 0.04 0.018 0.006 0.0025
the Michaelis-Menten equation dCA
(77)
The initial concentration of sucrose is Cao = 1.0 mmol/l, and that of sucrase is Ceo =0.01 mmol/1. The data listed in Table 2 were obtained with a batch reactor. The kinetics of the system can be described by
k3CEoC.t
dt
(78)
CA+M
Determine the parameters k 3 and M. Solution: for convenience, eq. (78) can be rewritten as
(75)
to determine the parameters k~ and k 2 by a graphic method, and obtained the following rate expression:
sucrose (A)
t (h)
dCA
-- k4Ca
dt
1 + ks CA
(79)
where k4=
k3CEo M
and
1 ks=~.
Since there is one set of data measured for state C a and one restrictive relationship (a total of two only), it is impossible to estimate state CA and parameters k4 and ks (a total of three) directly and simultaneously, and a trial-and-error procedure is required. There are two possible ways: (1) trying k4, and (2) trying ks. Which one should be selected could be determined from the sensitivity of the kinetics model to it. F o r the case discussed here, the model is more sensitive to k4 than to ks, or, in other words, the effect of the change in k4 on the reaction rate is larger than that of the change in k s . Therefore, it may be expected to get more accurate results if trying k4, that is, estimating Ca and k5 with various assumed values for k4. Let CA = x l and k 5 = x 2, then the kinetics equation can be rewritten as the following set: dx 1
k4x 1
.
dt
(81)
1 + x2x 1
dx 2
--=0
6
and the corresponding estimator equations are j 0 *0" -
4 0
EVa
cY
2 1 0
(82)
dt
5
(80)
/
/°
/°
de1
da
CR calculated by eq. (74)
k4et
1 I-- [ z l ( a ) - el] l+e2e 1 a
de 2
1 -
da 1 100
I
I
21)11
300
I 4OO
(84)
[z2(a)-e2].
The indirectly-measured values for ks, zz(a), can be represented by
t (mini
Fig. 2. Results simulated by eq. (74).
a
(83)
z2(a)=
-k,(del~-I \ da ]
1 el
(85)
Estimation of states and parameters by invariant imbedding technique = k4/0.16 - - 1/0.84.
from the kinetics equation, eq. (81). If the time of 1 h after the reaction begins is taken as the initial time, that is, if ao= 1 h
then the initial values for e 1 and e z can be obtained from the data listed in Table 2 as (87)
. [zl(ao+At)--za(ao)-] -x 1 e2(a")=zz(a°)= --k4" L ~ - z,(ao)
J
(a)
(88)
The results calculated with a step size of 0.1 h and ming various values of k4 are shown in Fig. 3. It seems that the convergence of the estimated values for CA, el, with various values of k4 assumed in a considerably large range are all very good; while the stabilities of those for k 5 are different. There is a visual difference even between the two results obtained with assuming k 5 =0.9994 and 0.9997. This suggests that the way of trying k4 is sensitive enough and consequently can result in good accuracy. The best results
(86)
e~(ao) = 0.84 mmol/l
k 4 = 0.999
k s = 5.054377 10
\\
0.I
¢o
0.6
8 e2
/
6
vE
E v
o,4
4
Zl
~
81
Z
0.2
n ] I 6
- - I l l L L I
I
I
L I 4
J n J I 8
J I
0
I I 10
a=t(h)
(b)
k 4 = 0.9994
k 5 = 5.061156 -10
1m
-8
0.8
0
E E ---= L~
"\
0.6
/
e2
E 0.4
zl/°~
°
je
II
1
0..~
I I I J l
--I
0
4
2
6
8
10
a = t (h)
t
0.8
O.E E E == (.J
(C)
k 5 = 5.071345
k4 = 1
-10
J
-\ "\ -
82
/
-
\
O.Z
~
el
0.2
'e ~o~e,....," 0
"6 E E II
zl--°~o
JlLll 2
1671
4
6 a =t(h)
Fig. 3. (a)-(c).
8
10
1672
Wu (d)
YUAN
k 4 == 1.1
ks =
7.124979
1 o.~
3 e2
E E
0.6
E ¢
It
0.2
L
I
J
I
,
I
,
I
2
J
I
,
I
4
a
I
J
6
I
t
I
~
8
n"T
10
a=t(h) (e)
1(4 =
0.9997
ks =
5.066253
10
1
"\\
0.8
-6 E E
O.E
=
r
8
/ e2
''0~
0.~
"\
Zl
6
-6 E vE
4 el
0.; 0 ~- I ) ' [ ' 0
I
,
2
l 4
J
I
~
1
6
,
l
~
-
I
,
[
i
~
£
I
10
a = t (h)
(f)
k4 =
k5 =
0.9
3.626837
-d10
1
0.8
\
0.6
E
8
"\
?-
"\
/
e2
,j= 0.4--
-6 E E
4
~II
2
0.2 ~
o- I l l l l o 2
-- 6
~0 4
6
10
8
a=t(h)
Fig. 3. Estimated and calculated results with various assumed values for ka.
are those with k4=0.9994, and the corresponding value estimated for ks is e2=5.0612 [see Fig. 3(b)]. More accurate results might be possible by further tries, but that is no longer necessary for a kinetics model. Therefore, the parameters are determined as
k+M =,0.9994,,0.19758,z/t~ 19.746h-1. k3 = CEo 0.01 -The kinetics of the system can thus be expressed as dCa
19.746CroCA -
ka=0.9994h -1,
k5=5.0612 l/mmol.
Correspondingly, from eq. (80): M = 1/k 5 = 1/5.0612=0.19758 mmol/l
dt
(89)
0.19758 + C A
The results simulated by eq. (89) are given in Fig. 4, and the standard error of the measured data is s = 0.005011 mmol/l.
Estimation of states and parameters by invariant imbedding technique
0.8
•
3E o.e E ¢~
e~e'=
Calculatedbyeq.(89)
0.4
0.;
ieasu red~e~e~e~e~e,..
0.0-
j
I 2
I
I 4
I
I 6 t(h)
I
~ I 8
"~'~'~-10
11
distinguishing models, the only requirement is one set of data on changes in states such as composition, temperature, pressure etc., with time or an equivalent independent variable. The experiments may be carried out simply with a batch reactor, provided that the measuring methods used are fast and accurate enough. Therefore, the estimator equations are of potential in simplifying experimental procedures and reducing the number of necessary experiments. In some cases, not all the states and the parameters can be directly and simultaneously estimated. In general, the criterion below should be satisfied:
Fig. 4. Results simulated by eq. (89)
Ne<,Nz+N r
The parameters k4 and k 5 can also be determined by trying k5, and the estimated values in such a way are k4=0.9686 h - 1 and k5 = 5.035 1/mmol, respectively. The standard error of the results calculated by the kinetics equation with these values of parameters is found to be s'-0.01362mmol/l, showing that the accuracy is lower than that by trying k4. By a graphical procedure, Levenspiel (1975) obtained the following kinetics equation for the same case: dCA d-T-=
16.4CEoCa 0.164 + C a
(90)
A simulation of eq. (90) leads to a standard error of sL = 0.05365 mmol/l. The accuracy is significantly lower than those obtained with the estimator equations for both ways of trying k4 and trying k5 and, in fact, a graphical procedure is usually troublesome. DISCUSSION In all of the cases discussed above, the kinetics equations selected are of the proper forms. The results calculated for these cases show that the N-dimensional estimator equations developed with the invariant imbedding concept can be used for the determination of parameters in kinetics research and, for a kinetics model of proper form, can result in a good accuracy, although the equations themselves are formally simple. To determine whether a kinetic model is proper or not, a procedure for distinguishing models is generally required. The estimator equations can also be used for such a procedure, although that has not been discussed in this paper. For each of the models to be distinguished, the principle and the method of the estimation calculation are the same as for the others. For complex systems, the estimator equations may present advantages over other methods currently used in kinetics research. Since all of the estimation calculation can be carried out sequentially on a computer, the estimator equations are more convenient for practical applications, and may eliminate some difficulties in experiments and interpreting data. The estimator equations require fewer and simpler data. For example, when including the procedure of
1673
(91)
where Ne is the total number of states and parameters to be estimated, N~ the number of states measured, and N, the number of independent restrictive relationships. If N e > ( N z + Nr), a trial-and-error procedure is needed. Which parameters should be selected for trying can be determined according to the sensitivities of the kinetics model to them in order to obtain more accurate results. It appears that a trial-and-error procedure requires a larger amount of calculation. However, the estimation calculation on a computer is rather convenient, even if the trial-and-error method has to be used. The estimator equations themselves do not require that the independent variable must be time and do not specify what variables the x~ must be. It can be expected that the N-dimensional estimator equations may find more and more applications in chemical engineering and other fields by appropriately setting the variables.
a C, C e, e
J r t X, X
Y Z, Z
NOTATION final value of independent variable t final value of state x or state vector x optimal-estimated value for state x or state vector x function to be minimized = y(a), missing final condition independent variable dependent variable, state or state vector of the system integral function defined by eq. (7) or (33) value measured for state x or state vector x
REFERENCES
Bellman, R. and Lee, E. S., 1977, Invariant imbedding and solid-fluid reaction--I. The unreacted core model• Chem. Engn9 Sci. 32, 1051-1054. Bullin, J. A. and Dukler, A. E., 1975, Hybrid computer modelling of a stochastic nonlinear dynamic system. Chem. Engn 9 Sci. 30, 631-35. Diamond, W. J., 1981, Practical Experimental Desion, pp. 76-83. Wadsworth, Belmont, CA. Lee, E. S., 1968a, Invariant imbedding, a versatile computation concept. Ind. Engng Chem. 60(9), 55-64. Lee, E. S., 1968b, Quasilinearization and the estimation of parameters in differential equations. Ind. Enon9 Chem. Fundam. 7, 152-158. Lee, E. S., 1968c, Invariant imbedding nonlinear filtering, and
1674
WU YUAN
parameter estimation. Ind. Engng Chem. Fundam. 7, 164-171. Levenspiel, O., 1972, Chemical Reaction Engineering, 2nd edition, pp. 88-89. John Wiley, New York. Levenspiel, O., 1975, Solutions to Problems in Chemical Reaction Engineering, 2nd edition. Oregon State University.
Omatu, J. A., 1986, Optimal estimation theory for distributed parameter systems, Control Dynamic Systems 23, 195-241. Rose, L. M., 1981, Chemical Reactor Design in Practice, pp. 127-137. Elsevier, Oxford. Stephanopoulos, G., 1984, Chemical Process Control, an Introduction to Theory and Practice, pp. 113-115. PrenticeHall, Englewood Cliffs, NJ.