Parameter identification for spatial xenon transient analysis and control

Annals of Nuclear Energy, Vol. 6, pp. 369 to 374 Pergamon Press Ltd 1979. Printed in Great Britain

PARAMETER I D E N T I F I C A T I O N FOR SPATIAL X E N O N TRANSIENT ANALYSIS A N D C O N T R O L R. J. ONEGA Virginia Polytechnic Institute and State University, Blacksburg, VA 24061, U.S.A. and R. A. KISNER

ORNL, Oak Ridge, TN 37830, U.S.A. Abstract--The real-time, on-line control of xenon-induced spatial flux oscillations in a large PWR requires the use of fairly simple models of the process for minicomputer implementation. Model simplification requires that the analyst reduce the physical phenomena to mathematical terms with coarser approximations and greater inaccuracies. In order to make a model's behavior fit that of a reactor as closely as possible, the estimation of parameters in the model requires special care. Three parameters were chosen for identification: the diffusion coefficient, D; the power reactivity coefficient, Cry; and the microscopic taSXe cross section, a x. The method of maximum likelihood was used to obtain estimates of these parameters by comparing calculated results with actual plant data. Numerical values of these parameters are presented.

1. INTRODUCTION A mathematical description of the dynamics of a nuclear power plant (PWR) subject to xenon-induced flux oscillations is necessary for the on-line control of the plant. Simple models can have a measure of realism if the parameters are adjusted so that the model generated reactor response is matched as closely as possible to the experimentally determined response. The parameter identification method has various i m p o r t a n t practical consequences, since uneven burnup, thermal cycling, and departure from nucleate boiling are concerns of the plant owner. A simple model with properly adjusted parameters could be used with an on-line computer (or operator) in a feedforward loop (anticipatory) or as a state observer to estimate the levels of the iodine and xenon in the reactor. In the feedforward mode, the model can be used to predict the response to perturbations in the state of the reactor such as an enexpected power transient. For example, a perturbation in the inlet coolant temperature could induce an axial fluctuation in the neutron flux which may go unnoticed for several hours due to the long time constants involved. The model described could easily predict the xenon oscillations induced as a result of this condition. 2. THE NEED FOR PARAMETER

IDENTIFICATION The need for parameter estimation is apparent. For a few group calculations, one may think that the few Asf~.6 7.8 A

group parameters are accurate representations of the cross sections, etc. involved in the model formulation. However, simple models "truncate the physical p h e n o m e n a " after a point. A complete model to describe spatial xenon oscillations would include the neutronics and thermal hydraulics. Most models omit the short neutronic time constants (delayed neutrons) and long time constants (fuel burnup) and lump the thermal effects in the p r o m p t power feedback parameter, ~r (Onega, 1978; Bell, 1970). The neutronics of the xenon spatial oscillation p h e n o m e n o n are complicated even with the lumping of thermal effects. The model generally does not include shielding of the xenon by the fuel. Average neutron cross sections are spectrum dependent. The neutron spectrum hardens as the oscillations continue, due to thermal neutron absorption. Also, even if spectral changes due to the presence of xenon itself are ignored, there is a need to evaluate the cross sections for the proper neutron spectrum which changes with temperature, i.e. as the power density changes, the spectrum changes locally. Admittedly, some of the effects will not be highly significant but it is good to at least check the model parameters to determine the adequacy of the parameters. One question which immediately confronts one when wishing to estimate the parameters is "which parameters should be estimated?" Physical insight a n d experience with a particular model provide the answer to this question. M a t h e m a t i c a l analysis may be helpful but it generally will not provide the insight needed. 369

370

R.J. ONEGA and R. A. KISNER

The model developed by Onega and Kisner (1978) was used for describing the axial spatial xenon oscillations in a reactor. Data were obtained for the neutron flux following an induced xenon oscillation. The three parameters that were estimated were the thermal diffusion coefficient D, the xenon absorption microscopic cross section ax, and the power coefficient of reactivity ~F. The reason for choosing these three parameters are as follows. First, the estimation of as few parameters as is practical is desired and, second, the parameters are to have physical significance. The diffusion coefficient has a marked effect on the period of oscillation and is not easily estimated from first principles for a power reactor. The xenon cross section G is neutron spectrum dependent and is also modified due to the shielding of tile xenon atoms by the fuel. Since the xenon density is low, xenon self-shielding is probably not important. The power coefficient of reactivity is again very difficult to estimate but it does take (to some extent) the thermal effects into account. Therefore, these were the three parameters that were estimated. Previously we obtained these parameters using one group averages over the reactor core.

dard deviation ai, so that

1

fi =

__ e x p ( - ½q/~/a2).

The probability for having the particular set {r/i } is the product of the probability density functions, so N

N

fill, O) = ~I f, = I] i: 1

i-. 1

PARAMETER

ESTIMATION

N

0,.] r,

i=l

(5)

i=l

Then, since tI is a function of 0 and we assume the ai are independent of O, we have (the gradient of L in parameter space 0) OL

~ (

1

ao - ao I - b-? [¢" - ¢(t, 0)]TIC" - 0 % 0)]

~ a,}.

,6,

Equation (6) indicates that maximizing L(O) is the same as minimizing 22(0) where we define )~2(0) = [¢"(t) - Oc(t, 0)] r [ 0 " ( t ) - •(t, 0)].

(1)

where Oj are the physical parameters in the system model. Of the various analytical methods available, the method of maximum likelihood is used in this analysis to obtain 0 from a comparison of the experimental flux at time t~ with the generated flux at the same time. The measured relative flux ~7' at a time t~ is corrupted by noise, so that (for N components or N experimental data points) (2)

where ~0~ is the calculated relative flux using the model and r/is the error or residual vector with components r/i. Both O~(t, 0) and tt(t, 0) are dependent upon the parameter reactor 0. The q~ are assumed to satisfy a Gaussian probability density function with zero mean and a stan-

(7)

So, differentiating Z2(0) with respect to 0, we have for c o m p o n e n t j, (setting al = 1 for all i) n

c'Z. ~ - t "0"/ = - v 2 [ ~ # ullj

~//m (t) = 0 c (t) -J- ~,

(4)

N

L = l n f ( 0 , r/) = - ½ ~ qi2/a~ - In \ / 2 n ~ ai.

~

02 . . . . .

__ e-'/2 ~ q{/a{. ffix//2n i =l

METHOD

The determination of the system parameters is an example of an inverse problem in engineering (Bekey, 1970). Observations of the relative neutron flux ~7' at time ti for tiE[0, t] after a transient is initiated are used to determine the parameter vector

0 = [01,

1

The basic assumption of the maximum likelihood f mlJ method is that the observed set of measurements ,0i is more likely to have come from the model with parameter vector 0 which maximizes f(O). Since the logarithm is a monotonic function, we define (Hsia, 1977)

-ln\/2~i=, 3. T H E

(3)

ai\/2n

i:

x

0 % , 0)3

1

Ol~(ti, O) - O, . j = 1 , 2 , . . . , M . ~0~

(8) The derivative c~O"(ti, O)/OOj must be evaluated for the M different parameters to be estimated. 4. T H E S Y S T E M

MODEL

The xenon oscillation model (Onega, 1978) is described in detail in the literature. However, notation must be specified. A P W R is assumed to be represented by two points. The flux, xenon and iodine concentrations are assumed to be expressible as n 2n ~(z, t) = c o s ~ z + A(t) sin ~ z, (9) 2n

x(z, t) = cos ~ z + B(t) sin --~ z H

(10)

Parameter identification for spatial xenon transient analysis and control and n

2n y(z, t) = cos ~ z + C(t) sin ~ z,

(11)

respectively. T h e average flux, x e n o n a n d iodine over the lower half of the core is o b t a i n e d by i n t e g r a t i o n of (9)-(11). T h e average relative flux ~k(t) (which c a n be c o m p a r e d with the experimentally m e a s u r e d flux in the lower half of the core) is 0

O(t) = f]

W(z,t)dzf] dz=-2[l-A(t,]

T h e one g r o u p n e u t r o n diffusion e q u a t i o n was satisfied as closely as possible by the functional forms (9)-(11) by m e a n s of a v a r i a t i o n a l procedure. The e x p a n s i o n coefficient for the flux A(t) can be d e t e r m i n e d from the diffusion equation, a n d is given by A(t)

HI2

(12,

=

v -

~/'v 2 + 1

(17)

and (~A

v -

1 -

~3v

0

H/2

371

(18)

x/v 2 + 1

where

It

- Xocrx + 7~bo~rZo

[8XoaxB(t)].

a n d the average x e n o n a n d iodine are (19)

x(t) = -211 - B(t)]

(13)

It

and 2 y(t) = - [1 - C(t)], It

(14)

respectively. T h e differential e q u a t i o n s describing the time b e h a v i o u r of the x e n o n a n d iodine can be written in terms of B(t) a n d c(t) as

T h e p a r a m e t e r s in the expression are presented in T a b l e 1. Initially, they were estimated using one g r o u p diffusion theory. F o r example, the diffusion coefficient was estimated by the relation D = 1/3Xr, a n d E r , was a v o l u m e weighted average over the core. O n e of the objectives of the study was to determine h o w accurate these o n e g r o u p estimates were. T h e m a x i m u m likelihood (7) can be rewritten as c9X2 00j

c9Z2 ~h~ t~A ~v ~B t95~ ~A c~v c~B t30j'

(20)

or (7) b e c o m e s - (~xq~o)

(A + / 3 ) + ~ AB

(15)

c~g2 _

2 ~ 2[0~,_ _ ¢ c ( t , , 0 ) ] [

~Oj

n ~ : l ai

1

=v

]

\IV 2 + 1:

and

dt

71Ey

A(t) - 21C(t).

(16)

x

v dB(t) B(ti) t?O#

0

(21)

Table 1. Definition of symbols and parameter values obtained from one-group diffusion theory for O C O N E E II plant.

Symbol

Unit

Parameter value obtained from one-group diffusion theory

)'1 ),~ 21 2x ax D H Ef red cte Eo (~- E~o) Ax* P

--sec -1 sec -I cm 2 cm cm cm- 1 cm- 1 cm 2 sec cm- 1 % MW

0.061 0.003 2.87 x 10 -s 2.09 × 10 5 2.72 × 10 -18 0.395 365.8 0.64 1.56 3.64 x 10 16 1.53 o.18 2500

G

MWsec fission

8.456 x 10 -21

-

Z ooJ

Z~o-

l

x 100.

Meaning Iodine fission yield Xenon fission yield 135I decay constant 135Xe decay constant 135Xe absorption cross section Diffusion coefficient Reactor core height Macroscopic fission cross section Neutrons emitted in fission/cm Prompt-power coefficient Average neutron absorption cross section Change in fixed absorber from equilibrium Eao Reactor operating power (thermal) Conversion factor

D~tred flux level ~o

'@

I

values

I

I

Xo '

4)0

-j

v

Calculate

I

X0

aV

Calculate

_•

Calculo~d ftuxes

Y(t), C[t)

v (t)

Catc~ta~e ~A

Calculate a___~A

Calculate ~A

B (t)

Fig. 1. Block diagram of the calculational procedure used in determining parameter space norm of the gradient Z2.

J

tJon

equilibrium

l

AZa--J .......

Solve variational estimate far flux

YO

~ez

Calculate ~X2

~102

Calculate ~X 2

catculate ~x 2

Obl"ain norm of VX 2

Parameter identification for spatial xenon transient analysis and control when the index j goes from 1 to M, and in our case we choose M to be 3. All quantities in (21) are determined except for ~B/~Oj. In order to obtain this quantity we must investigate the system sensitivity equations (Brogan, 1974). The system sensitivity equations tell how sensitive the dependent variables A(t), B(t) and C(t) are to a change in parameter 0i . The appropriate equations for the flux are

,?A _ ~ A

~v.

~30~

~v ~0~

~A

dA

c~02

?,v

~A 1 [8B__~) ] ~v 0, + v ,

(22)

oo/..0..,.,,,

(23)

and c~A

~03

-

c3A

c~v

(24)

[7dpoZ,~/8Xoa~OtB(t)].

373

where the superscript (°) indicates the original estimate of the parameter. The determination of 0 was done by using a computer program. The flow of the information was as indicated in Fig. 1. The desired step change in absorption cross section was calculated to simulate the insertion or removal of rods in either half of the core. The expansion coefficients A(t), B(t) and C(t) were then calculated and from these the sensitivities dA/~O~, dB/~Oj and ~C/OOj were obtained. Equation (21) was then used, since all information including the experimental data @7' was then available. The calculation was carried out for a given set of parameters for the 53.5 h transient from which the @7' were obtained. The best estimates for 0i were those for which the norm of the gradient was a minimum; i.e. +

ilvz ~il=,,~0,.~,, i

\ ~o,)

+

\ Jo3)

"

(28)

The xenon sensitivity equations are 5. RESULTS

d,L

¢J

+

+., +_;.,,j o

- °.

F2 / 8A

o°,btoo

+

n( dB B OA~7 +%,Ao¢+ ~)j,i=l,2

..... M, (25)

where 3~t is the Kronecker delta. Finally, the iodine sensitivity is given by

d r~C7

~tL~J

:[ °°IAo

, OC .

1,2 . . . . . M.

(26)

We choose M = 3 as indicated by (22-24). The parameter vector is given by

0 =

D/D° /

(27)

The results of the analysis are summarized in Table 2. The parameters are 'best' only in a local sense. To determine if the minimum in }rVZ2 It is global or local is very difficult as there may be other starting points for which smaller values of IIVz211 exist. The univariate search method (Cooper, 1970) was used for obtaining the optimal 0i . A few general remarks concerning the parameters can be made. The diffusion coefficient has the least effect on IlVzZlp. For example, when ctr = 3.64 × 10 -i6 and ax = 2.6 x 10 -18 , and if D is 0.355, then ]]Vz2ii is 0.068724, whereas if D is 0.425, then IIV•2]] is 0.086026. If other values of ~r and ax are used, the minimum ]IVz2]] changes for a given change in D, but not drastically. The curves ]WX2[] vs D for given values of ax and ~r do not have very steep slopes, although IIVz/]] is not a monotonic function of D in any reasonable range of D values. For example, using the above values of ~v and ax, the minimum ppvzZJJ for a D of 0.385 is 0.06619. Changes in the parameters ax and c~r change the minimum ]]Vz2[] drastically after certain values are reached. As an example of this, if D is 0.375 and a~

Table 2. Comparison of parameter estimates. Component* of 0

01 = ax/a~ 02 = D/D ° 03 = ~xe/o~~

Quantity

One group diffusion theory estimate

Optimized estimate

ax(cm 2) D(cm) o~i:(cm 2 sec) IIVzz II

2.72 × 10 18 0.395 3.64 × 10 -16 0.11367

2.60 x 10 -18 0.375 3.60 × 10 -16 0.024991

* The superscripted value indicates the one group diffusion theory estimates, e.g. D° = 0.395 cm.

374

R.J. ONEGA and R. A. KISNER

1.20 ~I~

~

~lO

/

oE I.O0

Experimental data ~.\ /

][

~r

O90

Model generated flux

~'~.

Modelgenerated flux II~xZll =0.11567

E zo

0.80 TrQnsieil- induced 0.70 O0

I IO0

200

I 300 Time, h

I 400

500

I 600

Fig. 2. Comparison of the time response of the calculated and experimentally obtained fluxes for the lower haft of the reactor core.

is 2.60 x 10 -18, where ~tr changes from 3.6 x 10 -~8 to 3.4 x 10 -18, then the minimum IIVz211 varies from 0.027 to 0.410. Figure 2 shows the experimental flux, the flux determined from the one-group diffusion theory estimates, and the flux determined from the optimized parameters. The previous article (Onega, 1978) used preliminary experimental data from a different run than the data used in this analysis. These data indicate that a simple, one-group one-dimensional model cannot duplicate the experimental results. Even though the model does not accurately duplicate the experimental data in every detail, the improvement of the optimized parameter results is apparent. 6. S U M M A R Y

The estimation of parameters is only one aspect of the identification of dynamic systems (Frogner, 1977). The structure of the model, whether to retain the nonlinearities involved, etc., are all important considerations in the mathematical description of the system. It is necessary to choose the correct number and the correct parameter vector in the estimation of parameters. Physical as well as mathematical consider-

ations dictate this choice. Several important aspects of the analysis are: (i) the parameters can be identified operationally under noisy conditions for this model; (ii) the parameters of interest may be optimized for particular purposes in order to provide input for taylored responses. For example, the proposed method for control of the xenon oscillation can be used to initiate an oscillation from which the parameters are to be optimized; and (iii) the parameters may be updated as fuel burnup occurs by running another xenon oscillation experiment at 'some convenient time'. Acknowledoement--The authors express their gratitude to J. Teachman for his considerable help in programing. REFERENCES

Bekey G. A. (1970) Simulation 15, 151. Bell G. I. and Glasstone S. (1970) Nuclear Reactor Theory. Von Nostrand Reinhold, New York, p. 556. Brogan W. L. Modern Control Theory. Quantum Publishers, New York, p. 340. Cooper L. and Steinberg, D. (1970) Introduction to Methods of Optimization. W. B. Saunders, Philadelphia, p. 158. Hsia T. C. (1974) System Identification. D. C. Heath, Lexington, p. 33. Onega R. J. and Kisner R. A. (1978) Ann. nuclear Energy 5, 13.