- Email: [email protected]
Real time application of simultaneous estimation theory for parameters and state to one-dimensional heat conduction
Nuclear Engineering and Design 81 (1984) 437-449 North-Holland, Amsterdam
437
REAL TIME APPLICATION OF SIMULTANEOUS ESTIMATION THEORY FOR PARAMETERS AND STATE TO ONE-DIMENSIONAL HEAT CONDUCTION Ryoichi T A K A H A S H I a n d A k i o A R A K A W A * Research Laboratory for Nuclear Reactors, Tokyo Institute of Technology, O-okayama, Meguro- ku, Tokyo 152, Japan Received 18 June 1983
This paper applies filtering theory of the parameters and states of a distributed system to one-dimensional heat conduction in order to obtain information on the real time application of nonlinear estimation theory. This involves a trial to study a simple and fast algorithm for the state estimation of power reactors limited to a small number of sensors by noisy observations. The authors estimated the heat transfer coefficients and the transient temperature profile in a copper rod simultaneously in real time by using a small computer, NOVA 3/12, where the parameter and state of the heat conduction-type P.D.E. could serve as a simplified alternative to the nuclear constant and power profile. The authors also examined the feasibility of a nonlinear filter provided for simultaneous estimation, considering unknown parameters as alternative state variables. The filter tested converged during the sampling period and generated reasonable estimates of the parameter and state from noisy observations with a small number of sensors and at a discrete number of times.
I. Introduction In J a p a n , light water reactors have served as suppliers to the base-load. These capacities are being used to generate a b o u t thirty percent of all the electricity, especially in the d a i l y o f f - p e a k d u r a t i o n . Since this t r e n d m a y be growing in the n e t w o r k of suppliers, l o a d - f o l l o w i n g o p e r a t i o n of L W R s shall be necessary for reasons of e c o n o m y in the n e a r future. A 1 2 - 3 - 6 - 3 cycle is representative of the average daily l o a d cycle for J a p a n e s e utility systems. This designation m e a n s that the system o p e r a t e s for 12 hours at 100% load. This is followed b y a decrease r a m p d o w n to 50% l o a d in 3 hours. T h e l o a d r e m a i n s at 50% for 6 hours a n d is then followed b y a l o a d increase r a m p b a c k to 100% in 3 hours. These are designed to use the c o n t r o l r o d s for r a p i d p o w e r change. T h e linear p o w e r d e n s i t y of the fuel r o d deviates when the p o w e r generation increases a n d decreases. It will be i m p o r t a n t to e x a m i n e the p r o b a b i l i t y of fuel fialure in the case where the r a p i d m o v e m e n t of the c o n t r o l r o d induces a change in p o w e r d e n s i t y so that the p e l l e t - c l a d interaction c a n occur in the n e i g h b o u r i n g a r e a of the c o n t r o l rod. T h u s the l o a d - f o l l o w i n g has at once the a d v a n t a g e of o p t i m a l l y e m p l o y i n g the reactor c a p a c i t y for off-peak d u r a t i o n a n d the d i s a d v a n t a g e of reinforcing the reliability of the i n s t r u m e n t a l system a n d avoiding fuel failure. W h e n using o n - p o w e r diagnostics as an aid to r e a c t o r control, it is useful to develop a simple a n d fast technique for e s t i m a t i n g nuclear constants a n d p o w e r density before the l o a d - f o l l o w i n g is driven b y necessity. By a n a l o g y with the case where the t i m e - d e p e n d e n t d i s t r i b u t i o n of p o w e r d e n s i t y is d e s c r i b e d a p p r o x i m a t e l y b y the heat c o n d u c t i o n equation, t e m p e r a t u r e profiles can be i n t r o d u c e d as a r e a s o n a b l e alternative to p o w e r density. T h e p r e s e n t authors, A r a k a w a a n d T a k a h a s h i [1], w h o e s t i m a t e d t e m p e r a t u r e profiles of a c o p p e r rod in heat conduction, showed that g o o d estimates were o b t a i n e d with a few sensors b y noisy observations. Here, e s t i m a t i o n of p a r a m e t e r s in the e q u a t i o n shall be c o n s i d e r e d b y experiments, for example, in e s t i m a t i n g t h e r m a l diffusivity a n d heat transfer coefficient of a r o d surface or in estimating the a b s o r p t i o n cross section of the fuel a s s e m b l y with uncertainty. * At present, Nippon Atomic Industry Group Co., Ltd., Ukishima-cho, Kawasaki-ku, Kawasaki 210, Japan. 0 0 2 9 - 5 4 9 3 / 8 4 / $ 0 3 . 0 0 © Elsevier Science Publishers B.V. ( N o r t h - H o l l a n d Physics Publishing Division)
R. TakahashL A. Arakawa / Real time estimation of heat conduction
438
Since estimating the parameters and state is a problem in adopting filtering theory for use with a distributed system, many works [2] have been widely published. However, most applications of theory have been restricted to computer simulations with data which were generated computationally by a physical model and the off-line filtering of real process data. For instance, Chen and Seinfeld [3] indicate from computer simulation that the steepest descent method provides for good estimation of an unknown parameter and state and also show that this technique requires us to solve a two-point boundary value problem, which usually results in a large and time-consuming calculation. The major disadvantage of these filters comes from their being nonlinear by nature. The filters should be devised to reduce the sensitivity to the noise level and to convert the algorithm into an efficient computation, so that these become available for actual processes. Their feasibility, however, will not always be clarified until applications are implemented in real time. This paper formulated a nonlinear filter for the parameter and state in a closed form from criteria of the weighted least squares estimation over space and time, where the state variables are recomposed by considering unknown parameters in the equation as the alternative state variables. The filtering algorithm was not iterative, but recursive to facilitate real time implementation. A solution of the present filter was obtained by expanding it in terms of eigenfunctions and polynomials. We traced an analogy between the heat conduction of a copper rod and time-dependent linear power density in the load-following process in order to study several questions related to the present filter simultaneous estimation of the parameter and state: (1) the effect of an initial guess of a parameter, (2) the effect of the level of the observation noise, and (3) the effect of the choice of the measurement scheme.
2. Statement of the problem
The experimental apparatus, a one-dimensional heat conduction system, is shown in Appendix A, which was recomposed from ref. [1]. In the present section, the parameters and state in a partial differential equation will be estimated simultaneously from noisy observations. The heat conduction system tested is assumed to be governed by the equation of the form,
OT(z't)= 3z ~-~{a(z) OT(z't) 3z } + yl(z)ul( t) + Y2u2(t) for
O<~z~z I,
=0,
for
z~
y2(z)=0,
for
0~
for
z.~
~fi(z)=Vo,
=Y0,
+
w(z, t),
(1)
where T(z, t) is the temperature profile, t is the time variable, and z is the spatial variable defined on the normalized domain [0,1]. Parameters a and fl were aggregated thermal diffusivity and heat transfer coefficient, respectively. The variable y~u1 + Y2u2 was the heater input. Here the system error w(z, t) was considered as resulting from modelling the physical process by a partial differential equation. The boundary conditions are expressed in the form, 3T O---ff= 0,
( z = 0, 1).
(2)
R. Takahashi, A. A r a k a w a / R e a l time estimation o f heat conduction
439
In the present analysis, observations yj (t~) may be carried out only at a discrete number of sensor locations z*, ( j = 1, 2 . . . . . M ) and at sampling instants k~ (i = 1, 2 .... ), by the form
(3)
yj(ti) = T(zT(t,), t,) +rlj(ti),
where T(z 7, ti) and ~lj(t~) represent the true temperature and the measurement noise appeared at z7 and t~, respectively. The problem is to determine the state T(z, t), the parameters a,/3 and the manipulates Ul(t ), u2(t ) from some criterion of least squares estimation. First, considering the unknown parameters as alternative state variables to recompose the system equation, the expanded state vector X(z, t) is defined by the form,
X(z, t)= col[V(z,
(4)
Hereafter, unknown parameters shall be limited to a and fl for reasons of simplicity of description. Furthermore, assuming Oa/Bt and OB/3t, the new state vector may be predicted by the equations, ~X a-7 = F ( T , ,~, /3, ul, u 2 ) + w,
z(O zl
F=
w ,
(5)
W ~
This emphasizes that we can never provide a model for the parameters of the partial differential equations. A scheme of estimation of the parameters and state will be derived from those statements of the problem. These filtered estimates will be determined so that a functional dependent on deviations of the estimates from observations is minimized. This criterion function J was chosen for a distributed system of continuous time by the quadratic form, M M
[y,(t)- T(z*(t), t)lo,j(t)[ yj(t ) -
J= 0
i
,)]dt
j
+ - - -I"lXllw(z,0R(z, z; t)w(z', t)dz'dzdt, Jo ,,o Jo
(6)
where Qij(t) represented the non-negative weight of observations. The weight R(z, z', t) was defined in such a manner that the inverse of R, R+(z, z', t) was written by f01 R + (z, z t , t ) R ( z t , s, t ) d z ' = 8 ( z - s ) .
(7)
The first term of eq. (6) minimizes the measurement error and provides the correlation at the point z ~ [0,1]. The second term does the same for the model error. The model error w(z, t), which represents the deviation of the modelling process from real heat conduction, could be introduced by fluctuations of coolant temperature, heat generation, property of a rod, and heat transfer characteristics. The boundary errors are neglected. The estimate of the parameter and state will be obtained by minimizing the functional eq. (6) subject to the constraint, eqs. (2) and (5). Applying variational calculus to the present problem and simultaneously reformulating the two-point boundary value problem into the initial value problem by using the sweep
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
440
method [1], [4], the filter equations which will generate an estimate
M
=
the forms,
,)ld,,
• (z, , ; + , ) : , ( z , 2(z, t;+,)
)((z, t~+ 1) take
M
5:(z, t ~ + a ) + E E e ( ~ , z*(tk+l),t£+,)hQkiy +if[ y j ( l k + l ) - i
~"(ZT, t k + l ) ] ,
(9)
i
e(~, ~', ,;+~)= e(z, ~', ,2)+ +R(z, z',
(8)
ti
A~e(~, ~', , ) + e ( z , ~', ,)A~,
t)]dt,
(10)
where derivation of these equations shall be given in Appendix B. Here, eqs. (8)-(10) will provide a predictor and corrector for obtaining the filtered estimates, since eq. (8) will predict the value of the process at the instant k£+ 1 based on the estimate X(z, t~£) by means of the model of process and then eq. (9) will correct the prediction by multiplying the innovations [y(t~.+~)-T(z, kk+l) ] by the time-dependent feedback gain P(z, z', t).
3. Lumping an filter for a distributed system
In the present paper, the space-dependent properties of the problem were postponed until numerical implementation was necessary for eqs. (8)-(10), while ordinary differential equations a n d / o r finite difference equations could be assumed for a distributed system at the beginning of the formulation. Thus the present filter had the advantage of considering the space-dependent problem in nature and introducing an approximation at the end. Here it was a problem of lumping this filter into a minicomputer-oriented scheme, which involved a trial to obtain a real-time estimation of the parameter and state in a copper rod by a small number of sensors. The solutions could be approximated by expanding the space-dependent variables which appeared in eqs. (8)-(10) in terms of eigenfunctions a n d / o r polynomials. The solutions are assumed with truncation of the series at N as follows: N
7r(z, t)=Y'.T,(t)q,i(z ),
(11)
i N
N
eTa(Z, z', t) = E E ePj(z ) p~e~,( z') . i
(12)
j
The trial functions used are the eigenfunctions of the homogeneous part of eq. (1) with the boundary conditions (2). These take the forms of
~i(z)
= 1,
for i = 1,
= v~ c o s ( / - 1)~rz,
for i = 2,3 . . . . . N.
(13)
The polynomials ( z ~- 1 } are assumed for representing parameters with truncation of the series at N~ and N~ as follows: Na
NO
& ( z ) = E a , z '-1,
~(z)= y'jS, z ~-a.
i
i
(14)
R. Takahashi, A. A r a k a w a / R e a l time estimation o f heat conduction
441
The remaining cross solutions take the forms of
NNo PT.( z, z', t) = P.T( Z, z', t) = E ~-,+,(z)P~"(t) z j=', i
j
NNo Pr~(Z, z', t) = Par(z, z', t) = E ~-.¢',(z)P~a(t) zj-~, i j N. No P.~(z,z t ,t)=~_.~_.z i - 1Pijecoe(t) zj-l, i
(15)
j N
N~
P~a(z, z', t)= Pp.(z, z', t)= E Ezi-lpSI3( t) z j-l, i j NaN a
Paa( z, z', t) = E E z ' - i p ~ ( i
t ) z+- '
i
We can solve for the coefficients of expansion, and hence the solution at any spatial point will be obtained by the inverse procedure. For simplicity of description, the several vectors and matrices were recomposed by the unknown coefficients in the forms, T(t) = col.[ .... ~ ( t ) .... ],
~ = c o L [ .... + .... ],
~ = c o l . [ .... ~j .... ],
,, = col.[ u ~ ( t / , u : ( t ) ] ,
p(t)= poT(t)
(16)
p~o(t) pro(t)] p~(t)
p~t~(t) ,
(17)
,,,~.o(t) ,,,~a(t)
t,a~(t)
T
J'oT=J'To = [ P ~ ( ' ) ] , Paa
oo( t )] ,
=
p~r=pra= [ p y ( ,)],
Pij
(18)
For the sake of simplicity, those vectors and matrices were defined by the expanding bases in the forms,
+=col.[+,(z)..... +N(z)], ~a = co1.[1 . . . . . z j - ' . . . . . z~=], ok#= col.[1 . . . . . z / - ' ..... zN~],
~=
(19)
f° o1 +T
0
.
g
After we substituted eqs. (11) and (14) into eq. (8) and integrated the equation on the domain, z ~ [0,1], after multiplying of the equation by g,j(z) at the left, we obtained the coupled ordinary differential
442
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
equations in the vector forms, , ( t ; + l ) = ~r(t;) + ft;+,{ [ a ( & ) - B ( / ~ ) ] i"(t) = C ( u ) } d t ,
a(,~) = < , , { , > , ~ }z>,
(20)
B ( k ) =
= (~,vT.),
•It --__col.[ yl(z),y2(z)] , where the symbol (x, y ) represented the spatial integration of the inner product of x and y on the domain z ~ [0,1]. After we substituted eqs. (15) and (19) into eq. (10) and integrated the equations on the domains z ~ (0,1] and z' ~ [0,1], after multiplying the equation by I/' at the left'and q,T at the right, we obtained the coupled ordinary differential equations in the matrix forms, p ( t ; + l ) = p (t~-)
+ f';+'[Gp(t)
G = [ G1
G3]'
0G2
~ = - (q,, f T ( t ) , t , e ) ,
+p(t)G T + R(t)]dt,
GI=A(&)-B(#)' ~ =v,
A similar procedure enabled us to obtain the solution for the coefficient .~ = col.[T, & fl], and hence the solution at any point was determined by ,,~'(Z, t~-+l ) = ~T(z)~(t~-+I), Jr(tk'+l) = 3C(tk+ 1) + P ( t [ + I ) F - I [
(22) Yk+l - H(tk+l)-fc(k;+l)],
p(k[+x)=p(k;+l)-p(k;+l)HV(tk+l)F 1H(tk+l)p(kk+l), F-' = I-l(tk+l)P(tk+l)HT(tk+l)
+ Q(tk+l) ,
(23)
(24) (25)
with the definitions of matrices by the forms,
u = [¢(,~),o],
[
<(z~'(,~)),
~,*(t~) =
(Pl( z ~ ( tk)),
(26)
... ~,N(z~(tk))]. (26)
.... *,,(zZ,(tk))
4. Real time estimation of parameter and state
Several weights were chosen for the present filter to estimate the parameter and state for heat conduction in a copper rod. The weighting function P(z, z', O) at the time zero was assumed to be
R. TakahashL A. Arakawa / Real time estimation of heat conduction
443
diagonal and space-independent, as given by the form,
PTT( Z,
Z', O) =
P2o3( z - z'),
PBt~(O)= P~,,
(27)
where the value of 10(° C) was used for P20. The weights R and Q in the criterion function were assumed to be space-independent, as given by the form, RTT(Z,
Z', t) =
- z')als,
Qu ( t ) = Qo23ij,
(28)
where the following parameters were used: R 0 = 0.002(°C/s), Q0 = 2.0 ( ° C ) , and At s = 5(s). The measurement error of the thermocouples was found to be small (--- _+ 1 ° C) in those experiments. The value of the error enabled us to define the performance of the filter at the instant t k in the form,
~(tk)=max[yj(tk)--
J'(Zj, tk)l,
(29)
J
where yj(tk) was the observation from the thermocouple at the location z7. Since nonlinearity of the filter cannot provide a mathematical tool for evaluating the performance by nature, a look at the time-dependent curves of c(tk) may be an alternative-solution. This procedure shall be used to study the efficiency of the filter throughout the present paper. Implementation of the simultaneous estimation of/3 and T(z, t) was discussed here as an example. At first the problem was to examine the effect of an initial guess of a parameter fl(0) on the estimated value/~. This may be directed toward simulating the following technical aspect. Linear power density may be evaluated with uncertainty by nuclear parameters of the cold clean state as, for instance, in the calculation of temperature profiles with some uncertainty based on the heat transfer coefficient cited in handbooks. The parameter and state in actual processes are not always obtained until operational data are accumulated through technical experiences. A simplified experiment demonstrates that estimates converge to the actual parameter and state when we started with a rough initial guess. Initially, the copper rod was at room temperature in the wind tunnel. The system was excited by a step change in the heater input and was regulated so that the noisy observation came at the single location z~" at every five-second sampling. The calculated coefficient by Hilpert's formula [5], 1.01 × 10 -2 ( l / s ) , was assumed for the initial guess of the heat transfer coefficient fl0. The axially uniform profile was postulated for an initial guess of the rod temperature. A value of 30( o C) was assumed for the initial guess J~(z, 0) throughout the present paper. The filter generated simultaneous estimates in fig. 1 as a function of/3o. There is an interesting curve in the case of the initial guess/3o- At an earlier transient, the filter was confused into thinking that an estimate of the coefficient was not close to the converged value. This illusion may be attributed to the condition where an earlier estimate of the temperature profile can generate a large apparent heat flux so that the large value is filtered for the heat transfer coefficient. This initial condition says alternatively that J'(z, 0) was higher than room temperature. It was shown in fig. 2 that the temperature profile was becoming corrected at the vicinity of z~' by the innovations and that the estimation error was being removed from left to right in the case of fl(0) = 2f10 and the single measurement location z~'. The second problem was to examine the effect of the noise level of the measurements on the estimates. The nonlinear filters are not always stabilized, since their nonlinearity is very sensitive to noisy observations. The stability condition for the nonlinear filters cannot always be devised until the filters are implemented in generating reasonable estimates. Since the temperature error for thermocouples was found to be very small (--- _+1°C) in those experiments, the computer-generated artificial noise was added to examine the effect of the noise level o on the estimated values of/~ and 7"(z, t). The pseudorandom signal, M-sequences with the zero-mean and the o-standard deviation was provided for the artificial noise. Only
444
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
the single location z~ was specified. The three cases o = 0; a = 2.5 and o = 5.0 ( ° C ) were studied as shown in fig. 3, where only the estimated value of/7 is illustrated for simplicity. The smoothed curve corresponds to the artificial error-free observations. There is no significant difference between the two noisy cases and then the estimates converge to the value labelled error free. Fig. 3 demonstrates that the filter generates reasonable estimates even in the presence of high-level noise. It was necessary here to study the effect of the choice of measurement strategies on the filter performance. A difficulty encountered in real time estimation is that it is time-consuming to recieve observations and solve coupled nonlinear equations, since the filter requires that all of the processes, prediction-innovation-correction, have been completed before the next sampling instant comes. The predictor inevitably depends on a large n u m b e r of computations which can remove the numerical instability involved in the nonlinearity of equations. Best estimates are essential for obtaining the observations continuously over space and time. This practice resulted in high computer costs. Here we will propose an alternate measurement scheme. The reciprocal traverse of measurement location may be introduced by providing that the filter receives the single observation from the new single location when the new sampling instant comes. This strategy shall be c o m p a r e d with the most expensive measurement which was performed at eight fixed locations at 5-second samplings. U n d e r the condition where strategies were made to equal the real time calculational burden on the computer, the following data were used in generating those results, as shown in fig. 4: a single location z~' at 5-second samplings, eight locations z~-z~ at 40-second samplings, a single location traversed from z~-z~, z~ to z~', z*-z* 6 3, z~' at 10-second samplings. Fig. 4 indicates that the reciprocal traverse of the sensing location produced the converged value of the parameter and state even at the earlier transient. Numerical divergence appeared at estimates when the large sampling period was used, as shown by the dotted curves in fig. 4. This, however, was eliminated as the parameter approached the value o f / 3 o . Reduction of error m a y be attributed to the condition that I
I
I
I[
I
I
i
Estimat:dat u I
I
I
1
/I
I
I
13(0)=13o
6o\
10.
[.4
5.
I
O.
0
I
I
100 I
I
I
200 I
I
I
TIME(sec) I
2. ~
"~ (X10-2 sec-1)
1.
0,
o
100
200
TIMEI sec)
Fig. 1. Effect of the initial guess of rio.
o
O.
I
I
I
i[i .5
I
]
1
1. Axial Direction
@ t=20 0 t=320 [] t=lO [] t=80 a, measured z~ t=40 at t=O ........ initial guess Fig. 2. Estimated time-dependent temperature profile 7"(~, t).
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
445
the coefficient of the partial differential equation converged to the true value, and consequently the converged parameter encouraged the mathematical model to identify with true physics.
5. Conclusion
By analogy with the practice which estimates the macroscopic nuclear constants and power distribution of the reactor core during the load-following transient, this paper applied filtering theory to one-dimensional heat conduction to obtain information on estimating the parameter and state of the partial differential equation. Utility of these filters will not be clarified because of their nonlinearity until the estimates are generated precisely during the time available. This paper was directed toward implementing the simultaneous estimation of the heat transfer coefficient fl and temperature profile T(z, t) of the copper rod. Experimental results demonstrated that the nonlinear filter tested converged in real time and then produced accurate estimates,/} and ~r(z, t), by noisy observations even at a single sensor. The sweeps of sensing locations provided a useful suggestion, so that the filter was devised for minimizing a heavy burden on the small computer. While this filter was implemented by the small computer NOVA 3/12, it will be important to optimize the capacity of the process computers for reasons of economy and safety in the case where estimation of the fast transient state associated with large changes in parameters is required.
2.1
i
~
j
I
L
i
I
c
i, ~J,l.~
b
.. 13 ( X l O -2 sec-,) -
I
I
I
I
10.
I
It
¢ (°C)
~'~v/~-a
O.
I
I
100
I
200
I
O.
I
0
TIME(sec) 2.
15.
~_1 I , b \, c
~,. C
Y
-
a
c (*C)
_
I
TIME(sec)
J
J
J
~ (X10_2 s e c _ l ) -
.,c
1.
o.
'~' .~ /4. ~:',l..~ ,,/,:-',~
5.
200
b
I/" 10.
100
~1 e i~ 0
100
I1
?1 200
?1 ? 19 TIME(sec)
Sampling instant in "c" 0.,
f
0
I
100
I
I
200
I
I
TIME(sec)
Fig. 3. Effect of the noise level; (a) o = 0, (b) o = 2.5, (c) o = 5.0.
Fig. 4. Effect of the measurement scheme; (a) scanning sensing, (b) single fixed measurement point ( M = 1, zl* , At, = 5 s), (c) 8 fixed measurement points ( M = 8, z{~-z~, At s = 40 s), (d) 8 fixed measurement points ( M = 8, z~'-z~, At s = 5 s).
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
446
Simultaneous estimation of the parameters and state, however, should give an applicable first step for further studies of optimal manipulation of reactor power.
Appendix A. Experimental apparatus A test section was a copper rod 35 cm long and 0.6 cm in diameter, as shown in fig. 5. The rod was electrically heated at both ends and thermally insulated at both faces. This was placed normal to the flow in a wind tunnel to assure uniform axial heat transfer because of the fact that the heater generated a relatively low heat flux. The wind tunnel was 170 cm long and 45 x 45 cm 2 in cross section and employed two straightening vanes at the entrance and in front of a fan. There were eight thermocouples soldered on the surface in an axial direction, avoiding at the flow separation point. Finally the Biot number was found from calculation to be 8.7 × 10 -4 ( << 1). This low Biot number enabled us to regard the present system as one-dimensional. The system was interfaced with a small computer, N O V A 3 / 1 2 through the input single multiplexer and amplifier leading to the A / D converter, while the heater, controlled in the interface, received a normalized voltage between 0 and 10 V from the D / A converter. The system is sketched in fig. 6, and the parameters are shown in table 1. Table 1 System parameters heat transfer coefficient a H = 46.5 k c a l / m 2h ° C coolant temperature T. = 22.0 o C
a = 0.816×10 -3 I / s /3 = 0.967 X i 0 -2 1/s "Yo= 2.68 × I 0 - 1 m3o C / W s u I,
uz=4.91×106W/m 3
4200
1500
K,
\1
TEST SECTION
I
-.,.':
~
II
I
--
f
l~
.:..! I
,7
I
I •P ~
'
I ~'q'
l
)t
T
I
]-z2:5.0
I~]-I:Lz~J/LI~ 1 -
I
,,
-,,, ':"
"
I
I|
1.1
l Jw I J /" I
.'f
I, ! L]_:, '--~I -~
',
ICE BOx
'
rT-----~ ~,r > u ° - - - f f - ~
i
INSULATOR
L_A_C, j_.~_v_] "XSIRAIGHTENIN~
i
;
VANES
[
(CM)
Fig. 5. A sketch of t h e test section and the wind tunnel.
t
HEATER CONTROLLER
I
o--@ FLOPPY DISK
I ,
J
NOVA 3
Fig. 6. System block diagram.
Appendix B. Derivation of the practical expression for filtering The practical expressions given by eqs. (8)-(10) will be derived to present the filter for a distributed system. This involved a trial to show procedures of prediction, innovation, and correction of the filter in the following forms of equations. The estimate of the parameter and state will be obtained by a straightforward
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
447
application of variational calculus, so that the quadratic functional, eq. (6), is minimized with the constraint, eq. (5). Thus the Kalman-type filter for a distributed system with continuous time and discrete measurement points takes the form, M
aX(z,t) at
M
F[2(z,t)]+Y'~_,P(z,z*(t) i j
'
t ) h Q , j ( t ) [ y j ( t ) - 7 " ( z T ( t ) t)] ' '
(30)
where J'(z, t) is filtered estimate of the parameter and stae X(z, t). The feedback gain, the so-called covariance value, P(z, z', t), could be found by the symmetric solution of the Riccati-type equation,
ae(z, ~', t) at
- A z e ( z , z', t)+e(z, z', t)AT,+ R(z, z', t) M
M
- ~_, ~_,P(z, z*(t), t)hQij(t)hTP(z~(t), i
z', t),
(31)
j
with the definitions of variables,
e(z, z', t)=
-PTT(Z, z',t )
eTo(z, z', t)
eTe(z, z', k)
coT(z, z', t)
P ~ ( z , z', t)
P~e(z, z', t)
PeT(z, z', t)
Pe~(z, z', t)
Pee(z, z', t)
IRTT(Z, Z', t) o
. ( z z,
(32)
EIj
and the definition of the operator,
ul(t)(')+uz(t)(') °
0 0
0 0
(33)
Boundary conditions concerned with ~" and P are given by
of" 0--~=0,
(z=O, 1),
OPrr az
OP~r aPr~ az az
aPer az
aPTT
aPar
aPra
OPTe
az'
az'
az'
az'
aPr~ az
0.
(z=0,1),
( z = 0 , 1).
(34)
Supposing that noisy measurements are obtained at the sampling instants, t = tk+ 1 (k = 0, 1, 2 .... ), the weight of observation Qij(t) could be expressed by the form, K
Qij(t) = ~-, ~ij°k+13(t - tk+l).
(35)
k
Considering eq. (35) and the property of the 8-function, the time integration of eq. (30) in the domain
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
448
(t~- ~< t ~
2(z,
t k++ , )
(36)
M
= 2(z, t;+l)+EEe(z, t 1"
z*(t~+,),
t+k + , , ~htOk+l[Yj(tk ~ij + 1)--J'(zi(tk+l), tk+l)]-
(37)
A look at fig. 7 may be helpful. The eq. (36) provides a prediction for the parameter and state. The noisy observations yj(tk+l) ( j = 1 . . . . . M) will be available for generating innovations [ y j ( t k + l ) J'(zj(t~-+ 1), tk-+1)]- The predicted value ~'(z, tk+ 1) can be corrected by combining the innovations and the feedback gain, P(z, z*(tk+l), tk+ + 1), as shown in eq. (37). Eq. (31) is easily integrated in the domain (t~ ~ t ~< t~+l) by considering the assumption, eq. (35). The resulting equation is in the form, ,;+,):
)+
f'k+'[ A z e ( z, z ,' t) + e ( z, z ,' t)A~,+ n( z, z ,' t)]dt,
(38)
where P(z, z', tk+l) represents prediction. Several formulations, however, are required to obtain P(z, z', t k++1) at the instance t k++a explicitly. Here the inverse of P, P+(z, z', t), is defined by the forms,
fo~P+(s, z, t)P(z, s', t)dz=lS(s-s'),
foP
(39)
(s, z, t)P+(z, s', t)dz=lS(s-s').
By multiplying eq. (31) at both the right and left by P+(z, z', t) and integrating the resulting equation triply with the regions z ~ [0,1], z' ~ [0,1] and the time domain, the equation can be formulated as M
M
P+(z, z', t;+l) = P( z, z', /k+l)q-E E I ~ ( S -- zt( t ) )hQk+ lhT~( S t - zt( k ) )"
(40)
i j The lemma of matrix inversion [6] leads us to find the forms, M
M
e( z, z', t ; + l ) = P ( z , z', tk+ l) - E E P ( z, z*( tk+l), tk+x)hFijl( tk+l)hTe( z/( lk+l), Z, tk+l) , l
/
(41)
['i]( tk+,) = PTT( Z*( t;+I), zt( tk+l), tk-+l) + Qij( tk+,)Finally, the filtered estimate of the parameter and state, X(z, t), will be produced by combining eqs. (8)-(10) in section 2.
observation(t~, I) .......Q 0
O...~estimate(t~',t) prediction(t~a)
..
I
I
t¢
t¢,+l
I time
Fig. 7. Prediction,observation,and estimate.
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
References [1] [2] [3] [4] [5] [6]
A. Arakawa and R. Takahashi, Nucl. Engrg. Des. 70 (1982) 215-221. W.H. Ray and D.G. Lainiotis, Distributed Parameter Systems (Marcel Dekker, Inc., New York, 1978). W.H. Chen and J.H. Seinfeld; Internat. J. Control 15 (1972) 487-495. J.S. Meditch, Automatica 7 (1971) 315-322. Dennetsu Kohgaku Shiryo (Heat Transfer Handbook), JSME 3rd-Edition (1977) (in Japanese). S.G. Tzafestas, Intemat. J. Control 15 (1972) 273-295.
449
437
REAL TIME APPLICATION OF SIMULTANEOUS ESTIMATION THEORY FOR PARAMETERS AND STATE TO ONE-DIMENSIONAL HEAT CONDUCTION Ryoichi T A K A H A S H I a n d A k i o A R A K A W A * Research Laboratory for Nuclear Reactors, Tokyo Institute of Technology, O-okayama, Meguro- ku, Tokyo 152, Japan Received 18 June 1983
This paper applies filtering theory of the parameters and states of a distributed system to one-dimensional heat conduction in order to obtain information on the real time application of nonlinear estimation theory. This involves a trial to study a simple and fast algorithm for the state estimation of power reactors limited to a small number of sensors by noisy observations. The authors estimated the heat transfer coefficients and the transient temperature profile in a copper rod simultaneously in real time by using a small computer, NOVA 3/12, where the parameter and state of the heat conduction-type P.D.E. could serve as a simplified alternative to the nuclear constant and power profile. The authors also examined the feasibility of a nonlinear filter provided for simultaneous estimation, considering unknown parameters as alternative state variables. The filter tested converged during the sampling period and generated reasonable estimates of the parameter and state from noisy observations with a small number of sensors and at a discrete number of times.
I. Introduction In J a p a n , light water reactors have served as suppliers to the base-load. These capacities are being used to generate a b o u t thirty percent of all the electricity, especially in the d a i l y o f f - p e a k d u r a t i o n . Since this t r e n d m a y be growing in the n e t w o r k of suppliers, l o a d - f o l l o w i n g o p e r a t i o n of L W R s shall be necessary for reasons of e c o n o m y in the n e a r future. A 1 2 - 3 - 6 - 3 cycle is representative of the average daily l o a d cycle for J a p a n e s e utility systems. This designation m e a n s that the system o p e r a t e s for 12 hours at 100% load. This is followed b y a decrease r a m p d o w n to 50% l o a d in 3 hours. T h e l o a d r e m a i n s at 50% for 6 hours a n d is then followed b y a l o a d increase r a m p b a c k to 100% in 3 hours. These are designed to use the c o n t r o l r o d s for r a p i d p o w e r change. T h e linear p o w e r d e n s i t y of the fuel r o d deviates when the p o w e r generation increases a n d decreases. It will be i m p o r t a n t to e x a m i n e the p r o b a b i l i t y of fuel fialure in the case where the r a p i d m o v e m e n t of the c o n t r o l r o d induces a change in p o w e r d e n s i t y so that the p e l l e t - c l a d interaction c a n occur in the n e i g h b o u r i n g a r e a of the c o n t r o l rod. T h u s the l o a d - f o l l o w i n g has at once the a d v a n t a g e of o p t i m a l l y e m p l o y i n g the reactor c a p a c i t y for off-peak d u r a t i o n a n d the d i s a d v a n t a g e of reinforcing the reliability of the i n s t r u m e n t a l system a n d avoiding fuel failure. W h e n using o n - p o w e r diagnostics as an aid to r e a c t o r control, it is useful to develop a simple a n d fast technique for e s t i m a t i n g nuclear constants a n d p o w e r density before the l o a d - f o l l o w i n g is driven b y necessity. By a n a l o g y with the case where the t i m e - d e p e n d e n t d i s t r i b u t i o n of p o w e r d e n s i t y is d e s c r i b e d a p p r o x i m a t e l y b y the heat c o n d u c t i o n equation, t e m p e r a t u r e profiles can be i n t r o d u c e d as a r e a s o n a b l e alternative to p o w e r density. T h e p r e s e n t authors, A r a k a w a a n d T a k a h a s h i [1], w h o e s t i m a t e d t e m p e r a t u r e profiles of a c o p p e r rod in heat conduction, showed that g o o d estimates were o b t a i n e d with a few sensors b y noisy observations. Here, e s t i m a t i o n of p a r a m e t e r s in the e q u a t i o n shall be c o n s i d e r e d b y experiments, for example, in e s t i m a t i n g t h e r m a l diffusivity a n d heat transfer coefficient of a r o d surface or in estimating the a b s o r p t i o n cross section of the fuel a s s e m b l y with uncertainty. * At present, Nippon Atomic Industry Group Co., Ltd., Ukishima-cho, Kawasaki-ku, Kawasaki 210, Japan. 0 0 2 9 - 5 4 9 3 / 8 4 / $ 0 3 . 0 0 © Elsevier Science Publishers B.V. ( N o r t h - H o l l a n d Physics Publishing Division)
R. TakahashL A. Arakawa / Real time estimation of heat conduction
438
Since estimating the parameters and state is a problem in adopting filtering theory for use with a distributed system, many works [2] have been widely published. However, most applications of theory have been restricted to computer simulations with data which were generated computationally by a physical model and the off-line filtering of real process data. For instance, Chen and Seinfeld [3] indicate from computer simulation that the steepest descent method provides for good estimation of an unknown parameter and state and also show that this technique requires us to solve a two-point boundary value problem, which usually results in a large and time-consuming calculation. The major disadvantage of these filters comes from their being nonlinear by nature. The filters should be devised to reduce the sensitivity to the noise level and to convert the algorithm into an efficient computation, so that these become available for actual processes. Their feasibility, however, will not always be clarified until applications are implemented in real time. This paper formulated a nonlinear filter for the parameter and state in a closed form from criteria of the weighted least squares estimation over space and time, where the state variables are recomposed by considering unknown parameters in the equation as the alternative state variables. The filtering algorithm was not iterative, but recursive to facilitate real time implementation. A solution of the present filter was obtained by expanding it in terms of eigenfunctions and polynomials. We traced an analogy between the heat conduction of a copper rod and time-dependent linear power density in the load-following process in order to study several questions related to the present filter simultaneous estimation of the parameter and state: (1) the effect of an initial guess of a parameter, (2) the effect of the level of the observation noise, and (3) the effect of the choice of the measurement scheme.
2. Statement of the problem
The experimental apparatus, a one-dimensional heat conduction system, is shown in Appendix A, which was recomposed from ref. [1]. In the present section, the parameters and state in a partial differential equation will be estimated simultaneously from noisy observations. The heat conduction system tested is assumed to be governed by the equation of the form,
OT(z't)= 3z ~-~{a(z) OT(z't) 3z } + yl(z)ul( t) + Y2u2(t) for
O<~z~z I,
=0,
for
z~
y2(z)=0,
for
0~
for
z.~
~fi(z)=Vo,
=Y0,
+
w(z, t),
(1)
where T(z, t) is the temperature profile, t is the time variable, and z is the spatial variable defined on the normalized domain [0,1]. Parameters a and fl were aggregated thermal diffusivity and heat transfer coefficient, respectively. The variable y~u1 + Y2u2 was the heater input. Here the system error w(z, t) was considered as resulting from modelling the physical process by a partial differential equation. The boundary conditions are expressed in the form, 3T O---ff= 0,
( z = 0, 1).
(2)
R. Takahashi, A. A r a k a w a / R e a l time estimation o f heat conduction
439
In the present analysis, observations yj (t~) may be carried out only at a discrete number of sensor locations z*, ( j = 1, 2 . . . . . M ) and at sampling instants k~ (i = 1, 2 .... ), by the form
(3)
yj(ti) = T(zT(t,), t,) +rlj(ti),
where T(z 7, ti) and ~lj(t~) represent the true temperature and the measurement noise appeared at z7 and t~, respectively. The problem is to determine the state T(z, t), the parameters a,/3 and the manipulates Ul(t ), u2(t ) from some criterion of least squares estimation. First, considering the unknown parameters as alternative state variables to recompose the system equation, the expanded state vector X(z, t) is defined by the form,
X(z, t)= col[V(z,
(4)
Hereafter, unknown parameters shall be limited to a and fl for reasons of simplicity of description. Furthermore, assuming Oa/Bt and OB/3t, the new state vector may be predicted by the equations, ~X a-7 = F ( T , ,~, /3, ul, u 2 ) + w,
z(O zl
F=
w ,
(5)
W ~
This emphasizes that we can never provide a model for the parameters of the partial differential equations. A scheme of estimation of the parameters and state will be derived from those statements of the problem. These filtered estimates will be determined so that a functional dependent on deviations of the estimates from observations is minimized. This criterion function J was chosen for a distributed system of continuous time by the quadratic form, M M
[y,(t)- T(z*(t), t)lo,j(t)[ yj(t ) -
J= 0
i
,)]dt
j
+ - - -I"lXllw(z,0R(z, z; t)w(z', t)dz'dzdt, Jo ,,o Jo
(6)
where Qij(t) represented the non-negative weight of observations. The weight R(z, z', t) was defined in such a manner that the inverse of R, R+(z, z', t) was written by f01 R + (z, z t , t ) R ( z t , s, t ) d z ' = 8 ( z - s ) .
(7)
The first term of eq. (6) minimizes the measurement error and provides the correlation at the point z ~ [0,1]. The second term does the same for the model error. The model error w(z, t), which represents the deviation of the modelling process from real heat conduction, could be introduced by fluctuations of coolant temperature, heat generation, property of a rod, and heat transfer characteristics. The boundary errors are neglected. The estimate of the parameter and state will be obtained by minimizing the functional eq. (6) subject to the constraint, eqs. (2) and (5). Applying variational calculus to the present problem and simultaneously reformulating the two-point boundary value problem into the initial value problem by using the sweep
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
440
method [1], [4], the filter equations which will generate an estimate
M
=
the forms,
,)ld,,
• (z, , ; + , ) : , ( z , 2(z, t;+,)
)((z, t~+ 1) take
M
5:(z, t ~ + a ) + E E e ( ~ , z*(tk+l),t£+,)hQkiy +if[ y j ( l k + l ) - i
~"(ZT, t k + l ) ] ,
(9)
i
e(~, ~', ,;+~)= e(z, ~', ,2)+ +R(z, z',
(8)
ti
A~e(~, ~', , ) + e ( z , ~', ,)A~,
t)]dt,
(10)
where derivation of these equations shall be given in Appendix B. Here, eqs. (8)-(10) will provide a predictor and corrector for obtaining the filtered estimates, since eq. (8) will predict the value of the process at the instant k£+ 1 based on the estimate X(z, t~£) by means of the model of process and then eq. (9) will correct the prediction by multiplying the innovations [y(t~.+~)-T(z, kk+l) ] by the time-dependent feedback gain P(z, z', t).
3. Lumping an filter for a distributed system
In the present paper, the space-dependent properties of the problem were postponed until numerical implementation was necessary for eqs. (8)-(10), while ordinary differential equations a n d / o r finite difference equations could be assumed for a distributed system at the beginning of the formulation. Thus the present filter had the advantage of considering the space-dependent problem in nature and introducing an approximation at the end. Here it was a problem of lumping this filter into a minicomputer-oriented scheme, which involved a trial to obtain a real-time estimation of the parameter and state in a copper rod by a small number of sensors. The solutions could be approximated by expanding the space-dependent variables which appeared in eqs. (8)-(10) in terms of eigenfunctions a n d / o r polynomials. The solutions are assumed with truncation of the series at N as follows: N
7r(z, t)=Y'.T,(t)q,i(z ),
(11)
i N
N
eTa(Z, z', t) = E E ePj(z ) p~e~,( z') . i
(12)
j
The trial functions used are the eigenfunctions of the homogeneous part of eq. (1) with the boundary conditions (2). These take the forms of
~i(z)
= 1,
for i = 1,
= v~ c o s ( / - 1)~rz,
for i = 2,3 . . . . . N.
(13)
The polynomials ( z ~- 1 } are assumed for representing parameters with truncation of the series at N~ and N~ as follows: Na
NO
& ( z ) = E a , z '-1,
~(z)= y'jS, z ~-a.
i
i
(14)
R. Takahashi, A. A r a k a w a / R e a l time estimation o f heat conduction
441
The remaining cross solutions take the forms of
NNo PT.( z, z', t) = P.T( Z, z', t) = E ~-,+,(z)P~"(t) z j=', i
j
NNo Pr~(Z, z', t) = Par(z, z', t) = E ~-.¢',(z)P~a(t) zj-~, i j N. No P.~(z,z t ,t)=~_.~_.z i - 1Pijecoe(t) zj-l, i
(15)
j N
N~
P~a(z, z', t)= Pp.(z, z', t)= E Ezi-lpSI3( t) z j-l, i j NaN a
Paa( z, z', t) = E E z ' - i p ~ ( i
t ) z+- '
i
We can solve for the coefficients of expansion, and hence the solution at any spatial point will be obtained by the inverse procedure. For simplicity of description, the several vectors and matrices were recomposed by the unknown coefficients in the forms, T(t) = col.[ .... ~ ( t ) .... ],
~ = c o L [ .... + .... ],
~ = c o l . [ .... ~j .... ],
,, = col.[ u ~ ( t / , u : ( t ) ] ,
p(t)= poT(t)
(16)
p~o(t) pro(t)] p~(t)
p~t~(t) ,
(17)
,,,~.o(t) ,,,~a(t)
t,a~(t)
T
J'oT=J'To = [ P ~ ( ' ) ] , Paa
oo( t )] ,
=
p~r=pra= [ p y ( ,)],
Pij
(18)
For the sake of simplicity, those vectors and matrices were defined by the expanding bases in the forms,
+=col.[+,(z)..... +N(z)], ~a = co1.[1 . . . . . z j - ' . . . . . z~=], ok#= col.[1 . . . . . z / - ' ..... zN~],
~=
(19)
f° o1 +T
0
.
g
After we substituted eqs. (11) and (14) into eq. (8) and integrated the equation on the domain, z ~ [0,1], after multiplying of the equation by g,j(z) at the left, we obtained the coupled ordinary differential
442
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
equations in the vector forms, , ( t ; + l ) = ~r(t;) + ft;+,{ [ a ( & ) - B ( / ~ ) ] i"(t) = C ( u ) } d t ,
a(,~) = < , , { , > , ~ }z>,
(20)
B ( k ) =
= (~,vT.),
•It --__col.[ yl(z),y2(z)] , where the symbol (x, y ) represented the spatial integration of the inner product of x and y on the domain z ~ [0,1]. After we substituted eqs. (15) and (19) into eq. (10) and integrated the equations on the domains z ~ (0,1] and z' ~ [0,1], after multiplying the equation by I/' at the left'and q,T at the right, we obtained the coupled ordinary differential equations in the matrix forms, p ( t ; + l ) = p (t~-)
+ f';+'[Gp(t)
G = [ G1
G3]'
0G2
~ = - (q,, f T ( t ) , t , e ) ,
+p(t)G T + R(t)]dt,
GI=A(&)-B(#)' ~ =v,
A similar procedure enabled us to obtain the solution for the coefficient .~ = col.[T, & fl], and hence the solution at any point was determined by ,,~'(Z, t~-+l ) = ~T(z)~(t~-+I), Jr(tk'+l) = 3C(tk+ 1) + P ( t [ + I ) F - I [
(22) Yk+l - H(tk+l)-fc(k;+l)],
p(k[+x)=p(k;+l)-p(k;+l)HV(tk+l)F 1H(tk+l)p(kk+l), F-' = I-l(tk+l)P(tk+l)HT(tk+l)
+ Q(tk+l) ,
(23)
(24) (25)
with the definitions of matrices by the forms,
u = [¢(,~),o],
[
<(z~'(,~)),
~,*(t~) =
(Pl( z ~ ( tk)),
(26)
... ~,N(z~(tk))]. (26)
.... *,,(zZ,(tk))
4. Real time estimation of parameter and state
Several weights were chosen for the present filter to estimate the parameter and state for heat conduction in a copper rod. The weighting function P(z, z', O) at the time zero was assumed to be
R. TakahashL A. Arakawa / Real time estimation of heat conduction
443
diagonal and space-independent, as given by the form,
PTT( Z,
Z', O) =
P2o3( z - z'),
PBt~(O)= P~,,
(27)
where the value of 10(° C) was used for P20. The weights R and Q in the criterion function were assumed to be space-independent, as given by the form, RTT(Z,
Z', t) =
- z')als,
Qu ( t ) = Qo23ij,
(28)
where the following parameters were used: R 0 = 0.002(°C/s), Q0 = 2.0 ( ° C ) , and At s = 5(s). The measurement error of the thermocouples was found to be small (--- _+ 1 ° C) in those experiments. The value of the error enabled us to define the performance of the filter at the instant t k in the form,
~(tk)=max[yj(tk)--
J'(Zj, tk)l,
(29)
J
where yj(tk) was the observation from the thermocouple at the location z7. Since nonlinearity of the filter cannot provide a mathematical tool for evaluating the performance by nature, a look at the time-dependent curves of c(tk) may be an alternative-solution. This procedure shall be used to study the efficiency of the filter throughout the present paper. Implementation of the simultaneous estimation of/3 and T(z, t) was discussed here as an example. At first the problem was to examine the effect of an initial guess of a parameter fl(0) on the estimated value/~. This may be directed toward simulating the following technical aspect. Linear power density may be evaluated with uncertainty by nuclear parameters of the cold clean state as, for instance, in the calculation of temperature profiles with some uncertainty based on the heat transfer coefficient cited in handbooks. The parameter and state in actual processes are not always obtained until operational data are accumulated through technical experiences. A simplified experiment demonstrates that estimates converge to the actual parameter and state when we started with a rough initial guess. Initially, the copper rod was at room temperature in the wind tunnel. The system was excited by a step change in the heater input and was regulated so that the noisy observation came at the single location z~" at every five-second sampling. The calculated coefficient by Hilpert's formula [5], 1.01 × 10 -2 ( l / s ) , was assumed for the initial guess of the heat transfer coefficient fl0. The axially uniform profile was postulated for an initial guess of the rod temperature. A value of 30( o C) was assumed for the initial guess J~(z, 0) throughout the present paper. The filter generated simultaneous estimates in fig. 1 as a function of/3o. There is an interesting curve in the case of the initial guess/3o- At an earlier transient, the filter was confused into thinking that an estimate of the coefficient was not close to the converged value. This illusion may be attributed to the condition where an earlier estimate of the temperature profile can generate a large apparent heat flux so that the large value is filtered for the heat transfer coefficient. This initial condition says alternatively that J'(z, 0) was higher than room temperature. It was shown in fig. 2 that the temperature profile was becoming corrected at the vicinity of z~' by the innovations and that the estimation error was being removed from left to right in the case of fl(0) = 2f10 and the single measurement location z~'. The second problem was to examine the effect of the noise level of the measurements on the estimates. The nonlinear filters are not always stabilized, since their nonlinearity is very sensitive to noisy observations. The stability condition for the nonlinear filters cannot always be devised until the filters are implemented in generating reasonable estimates. Since the temperature error for thermocouples was found to be very small (--- _+1°C) in those experiments, the computer-generated artificial noise was added to examine the effect of the noise level o on the estimated values of/~ and 7"(z, t). The pseudorandom signal, M-sequences with the zero-mean and the o-standard deviation was provided for the artificial noise. Only
444
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
the single location z~ was specified. The three cases o = 0; a = 2.5 and o = 5.0 ( ° C ) were studied as shown in fig. 3, where only the estimated value of/7 is illustrated for simplicity. The smoothed curve corresponds to the artificial error-free observations. There is no significant difference between the two noisy cases and then the estimates converge to the value labelled error free. Fig. 3 demonstrates that the filter generates reasonable estimates even in the presence of high-level noise. It was necessary here to study the effect of the choice of measurement strategies on the filter performance. A difficulty encountered in real time estimation is that it is time-consuming to recieve observations and solve coupled nonlinear equations, since the filter requires that all of the processes, prediction-innovation-correction, have been completed before the next sampling instant comes. The predictor inevitably depends on a large n u m b e r of computations which can remove the numerical instability involved in the nonlinearity of equations. Best estimates are essential for obtaining the observations continuously over space and time. This practice resulted in high computer costs. Here we will propose an alternate measurement scheme. The reciprocal traverse of measurement location may be introduced by providing that the filter receives the single observation from the new single location when the new sampling instant comes. This strategy shall be c o m p a r e d with the most expensive measurement which was performed at eight fixed locations at 5-second samplings. U n d e r the condition where strategies were made to equal the real time calculational burden on the computer, the following data were used in generating those results, as shown in fig. 4: a single location z~' at 5-second samplings, eight locations z~-z~ at 40-second samplings, a single location traversed from z~-z~, z~ to z~', z*-z* 6 3, z~' at 10-second samplings. Fig. 4 indicates that the reciprocal traverse of the sensing location produced the converged value of the parameter and state even at the earlier transient. Numerical divergence appeared at estimates when the large sampling period was used, as shown by the dotted curves in fig. 4. This, however, was eliminated as the parameter approached the value o f / 3 o . Reduction of error m a y be attributed to the condition that I
I
I
I[
I
I
i
Estimat:dat u I
I
I
1
/I
I
I
13(0)=13o
6o\
10.
[.4
5.
I
O.
0
I
I
100 I
I
I
200 I
I
I
TIME(sec) I
2. ~
"~ (X10-2 sec-1)
1.
0,
o
100
200
TIMEI sec)
Fig. 1. Effect of the initial guess of rio.
o
O.
I
I
I
i[i .5
I
]
1
1. Axial Direction
@ t=20 0 t=320 [] t=lO [] t=80 a, measured z~ t=40 at t=O ........ initial guess Fig. 2. Estimated time-dependent temperature profile 7"(~, t).
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
445
the coefficient of the partial differential equation converged to the true value, and consequently the converged parameter encouraged the mathematical model to identify with true physics.
5. Conclusion
By analogy with the practice which estimates the macroscopic nuclear constants and power distribution of the reactor core during the load-following transient, this paper applied filtering theory to one-dimensional heat conduction to obtain information on estimating the parameter and state of the partial differential equation. Utility of these filters will not be clarified because of their nonlinearity until the estimates are generated precisely during the time available. This paper was directed toward implementing the simultaneous estimation of the heat transfer coefficient fl and temperature profile T(z, t) of the copper rod. Experimental results demonstrated that the nonlinear filter tested converged in real time and then produced accurate estimates,/} and ~r(z, t), by noisy observations even at a single sensor. The sweeps of sensing locations provided a useful suggestion, so that the filter was devised for minimizing a heavy burden on the small computer. While this filter was implemented by the small computer NOVA 3/12, it will be important to optimize the capacity of the process computers for reasons of economy and safety in the case where estimation of the fast transient state associated with large changes in parameters is required.
2.1
i
~
j
I
L
i
I
c
i, ~J,l.~
b
.. 13 ( X l O -2 sec-,) -
I
I
I
I
10.
I
It
¢ (°C)
~'~v/~-a
O.
I
I
100
I
200
I
O.
I
0
TIME(sec) 2.
15.
~_1 I , b \, c
~,. C
Y
-
a
c (*C)
_
I
TIME(sec)
J
J
J
~ (X10_2 s e c _ l ) -
.,c
1.
o.
'~' .~ /4. ~:',l..~ ,,/,:-',~
5.
200
b
I/" 10.
100
~1 e i~ 0
100
I1
?1 200
?1 ? 19 TIME(sec)
Sampling instant in "c" 0.,
f
0
I
100
I
I
200
I
I
TIME(sec)
Fig. 3. Effect of the noise level; (a) o = 0, (b) o = 2.5, (c) o = 5.0.
Fig. 4. Effect of the measurement scheme; (a) scanning sensing, (b) single fixed measurement point ( M = 1, zl* , At, = 5 s), (c) 8 fixed measurement points ( M = 8, z{~-z~, At s = 40 s), (d) 8 fixed measurement points ( M = 8, z~'-z~, At s = 5 s).
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
446
Simultaneous estimation of the parameters and state, however, should give an applicable first step for further studies of optimal manipulation of reactor power.
Appendix A. Experimental apparatus A test section was a copper rod 35 cm long and 0.6 cm in diameter, as shown in fig. 5. The rod was electrically heated at both ends and thermally insulated at both faces. This was placed normal to the flow in a wind tunnel to assure uniform axial heat transfer because of the fact that the heater generated a relatively low heat flux. The wind tunnel was 170 cm long and 45 x 45 cm 2 in cross section and employed two straightening vanes at the entrance and in front of a fan. There were eight thermocouples soldered on the surface in an axial direction, avoiding at the flow separation point. Finally the Biot number was found from calculation to be 8.7 × 10 -4 ( << 1). This low Biot number enabled us to regard the present system as one-dimensional. The system was interfaced with a small computer, N O V A 3 / 1 2 through the input single multiplexer and amplifier leading to the A / D converter, while the heater, controlled in the interface, received a normalized voltage between 0 and 10 V from the D / A converter. The system is sketched in fig. 6, and the parameters are shown in table 1. Table 1 System parameters heat transfer coefficient a H = 46.5 k c a l / m 2h ° C coolant temperature T. = 22.0 o C
a = 0.816×10 -3 I / s /3 = 0.967 X i 0 -2 1/s "Yo= 2.68 × I 0 - 1 m3o C / W s u I,
uz=4.91×106W/m 3
4200
1500
K,
\1
TEST SECTION
I
-.,.':
~
II
I
--
f
l~
.:..! I
,7
I
I •P ~
'
I ~'q'
l
)t
T
I
]-z2:5.0
I~]-I:Lz~J/LI~ 1 -
I
,,
-,,, ':"
"
I
I|
1.1
l Jw I J /" I
.'f
I, ! L]_:, '--~I -~
',
ICE BOx
'
rT-----~ ~,r > u ° - - - f f - ~
i
INSULATOR
L_A_C, j_.~_v_] "XSIRAIGHTENIN~
i
;
VANES
[
(CM)
Fig. 5. A sketch of t h e test section and the wind tunnel.
t
HEATER CONTROLLER
I
o--@ FLOPPY DISK
I ,
J
NOVA 3
Fig. 6. System block diagram.
Appendix B. Derivation of the practical expression for filtering The practical expressions given by eqs. (8)-(10) will be derived to present the filter for a distributed system. This involved a trial to show procedures of prediction, innovation, and correction of the filter in the following forms of equations. The estimate of the parameter and state will be obtained by a straightforward
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
447
application of variational calculus, so that the quadratic functional, eq. (6), is minimized with the constraint, eq. (5). Thus the Kalman-type filter for a distributed system with continuous time and discrete measurement points takes the form, M
aX(z,t) at
M
F[2(z,t)]+Y'~_,P(z,z*(t) i j
'
t ) h Q , j ( t ) [ y j ( t ) - 7 " ( z T ( t ) t)] ' '
(30)
where J'(z, t) is filtered estimate of the parameter and stae X(z, t). The feedback gain, the so-called covariance value, P(z, z', t), could be found by the symmetric solution of the Riccati-type equation,
ae(z, ~', t) at
- A z e ( z , z', t)+e(z, z', t)AT,+ R(z, z', t) M
M
- ~_, ~_,P(z, z*(t), t)hQij(t)hTP(z~(t), i
z', t),
(31)
j
with the definitions of variables,
e(z, z', t)=
-PTT(Z, z',t )
eTo(z, z', t)
eTe(z, z', k)
coT(z, z', t)
P ~ ( z , z', t)
P~e(z, z', t)
PeT(z, z', t)
Pe~(z, z', t)
Pee(z, z', t)
IRTT(Z, Z', t) o
. ( z z,
(32)
EIj
and the definition of the operator,
ul(t)(')+uz(t)(') °
0 0
0 0
(33)
Boundary conditions concerned with ~" and P are given by
of" 0--~=0,
(z=O, 1),
OPrr az
OP~r aPr~ az az
aPer az
aPTT
aPar
aPra
OPTe
az'
az'
az'
az'
aPr~ az
0.
(z=0,1),
( z = 0 , 1).
(34)
Supposing that noisy measurements are obtained at the sampling instants, t = tk+ 1 (k = 0, 1, 2 .... ), the weight of observation Qij(t) could be expressed by the form, K
Qij(t) = ~-, ~ij°k+13(t - tk+l).
(35)
k
Considering eq. (35) and the property of the 8-function, the time integration of eq. (30) in the domain
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
448
(t~- ~< t ~
2(z,
t k++ , )
(36)
M
= 2(z, t;+l)+EEe(z, t 1"
z*(t~+,),
t+k + , , ~htOk+l[Yj(tk ~ij + 1)--J'(zi(tk+l), tk+l)]-
(37)
A look at fig. 7 may be helpful. The eq. (36) provides a prediction for the parameter and state. The noisy observations yj(tk+l) ( j = 1 . . . . . M) will be available for generating innovations [ y j ( t k + l ) J'(zj(t~-+ 1), tk-+1)]- The predicted value ~'(z, tk+ 1) can be corrected by combining the innovations and the feedback gain, P(z, z*(tk+l), tk+ + 1), as shown in eq. (37). Eq. (31) is easily integrated in the domain (t~ ~ t ~< t~+l) by considering the assumption, eq. (35). The resulting equation is in the form, ,;+,):
)+
f'k+'[ A z e ( z, z ,' t) + e ( z, z ,' t)A~,+ n( z, z ,' t)]dt,
(38)
where P(z, z', tk+l) represents prediction. Several formulations, however, are required to obtain P(z, z', t k++1) at the instance t k++a explicitly. Here the inverse of P, P+(z, z', t), is defined by the forms,
fo~P+(s, z, t)P(z, s', t)dz=lS(s-s'),
foP
(39)
(s, z, t)P+(z, s', t)dz=lS(s-s').
By multiplying eq. (31) at both the right and left by P+(z, z', t) and integrating the resulting equation triply with the regions z ~ [0,1], z' ~ [0,1] and the time domain, the equation can be formulated as M
M
P+(z, z', t;+l) = P( z, z', /k+l)q-E E I ~ ( S -- zt( t ) )hQk+ lhT~( S t - zt( k ) )"
(40)
i j The lemma of matrix inversion [6] leads us to find the forms, M
M
e( z, z', t ; + l ) = P ( z , z', tk+ l) - E E P ( z, z*( tk+l), tk+x)hFijl( tk+l)hTe( z/( lk+l), Z, tk+l) , l
/
(41)
['i]( tk+,) = PTT( Z*( t;+I), zt( tk+l), tk-+l) + Qij( tk+,)Finally, the filtered estimate of the parameter and state, X(z, t), will be produced by combining eqs. (8)-(10) in section 2.
observation(t~, I) .......Q 0
O...~estimate(t~',t) prediction(t~a)
..
I
I
t¢
t¢,+l
I time
Fig. 7. Prediction,observation,and estimate.
R. Takahashi, A. Arakawa / Real time estimation of heat conduction
References [1] [2] [3] [4] [5] [6]
A. Arakawa and R. Takahashi, Nucl. Engrg. Des. 70 (1982) 215-221. W.H. Ray and D.G. Lainiotis, Distributed Parameter Systems (Marcel Dekker, Inc., New York, 1978). W.H. Chen and J.H. Seinfeld; Internat. J. Control 15 (1972) 487-495. J.S. Meditch, Automatica 7 (1971) 315-322. Dennetsu Kohgaku Shiryo (Heat Transfer Handbook), JSME 3rd-Edition (1977) (in Japanese). S.G. Tzafestas, Intemat. J. Control 15 (1972) 273-295.
449