A systematic method to identify nonlinear dynamics of BWR by using the reactor noise

Progress in Nuclear Energy, Vol. 43, No. 1-4. pp. 209-216, 2003 Availableonline at www.sciencedirect.com ~} se,B.¢s o~.R¢'ro

Pergamon

© 2003 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0149-1970/03/$ - see front matter

www.elsevier.com/locate/pnucene

do|: 10.1016/S0149-1970(03)00030-1

A SYSTEMATIC M E T H O D TO IDENTIFY N O N L I N E A R DYNAMICS O F B W R B Y U S I N G T H E R E A C T O R

NOISE 1

LIU WENFENG, LUO ZHENGPEI, LI FU, WANG YAQI Institute of Nuclear Energy Technology, Tsinghua University, Beijing 100084, China

ABSTRACT For the identification of the dynamics of the Vermont Yankee BWR with the reactor noise, different parametric models have been tested. The widely used ARMA model is unable to identify the nonlinearity in the noise data. A systematic method by using the NARMA model, which takes advantage of both the ANN and ARMA, is developed. Comparisons are made between the identification results with ARMA and NARMA model. The advantages of identification with NARMA model over ARMA model are demonstrated. The linear-kernels of the identified NARMA models are extracted so that the natural frequency, damping ratio and time constants of the BWR are obtained. The values of those characteristics are well corresponded with the eigenvalues calculated by the differential equations of the Vermont Yankee BWR. The damping ratio with negative value is found to be a criterion for the existence of limit-cycle, which can be seen from the impulse response on the (X,, X,_~) plane, in stable nonlinear system. © 2003 Elsevier Science Ltd. All rights reserved. KEYWORDS ARMA, NARMA, ANN, limit-cycle, system identification 1. INTRODUCTION The BWR is susceptible to instability when operated under high power and low flow conditions (Gialdi et al., 1985). It has been found for decades that the observed power oscillation in BWR is mostly the limitcycle attractor in the nonlinear system (March-Leuba et al., 1986a, March-Leuba et al., 1986b). Deducing differential equations under appropriate assumptions is a general approach to study the affecting factors on the limit-cycle oscillation. Using different differential models and focusing research on various facets of the nuclear reactor system, abundant results have been obtained. For the instability of a heated l Supported by the Natural Science Foundation of China (Project number 19875028) 209

210

Liu Wen/eng et al.

channel caused by the two-phase flow, the effects of the pressure drop and the axial heat flux on the nonlinear dynamics of the system were derived by using a homogeneous equilibrium model (Narayanan et al., 1997). The effects of the void distribution and the drift velocity were obtained by using a drift flux model (Rizwan-uddin et al., 1986). Using the bifurcation theory, Van Bragt studied how the void distribution parameter and axial power profile affect the limit cycle in the heated channel as well as in the nuclear reactor (with the nuclear feedback) (Van Bragt et al., 2000), Masashi analyzed effects of various parameters such as the pressure loss coefficient, subcooling of inlet coolant, void reactivity feedback coefficient on the stability boundary (Masashi et al., 1993). The differential models can give a clear physical interpretation for the nonlinear dynamics of the system. However, for the problem to identify the nonlinear dynamics of the reactor in operation, they may not be eligible. The identification problem assumes that the reactor system is a black box and the various stochastic signals detected in the reactor are the responses for the inherent disturbance on the black box. Constructing models representing the dynamics of the system by using reactor noise data is an inverse task without prior knowledge of the parameters of the system. Traditionally, parametric models such as AR or ARMA are constructed to analyze the dynamics (Upadhyaya et al., 1981, Luo et al., 1996). Those approaches for identifying the dynamics of nuclear reactors are based on the assumption that the nuclear reactor is a linear system. It's true that those models with appropriate order can represent the dynamics of nuclear reactors in many cases. In effect, nuclear reactor systems are not linear, especially in the BWR when the power feedback is strong. To identify the nonlinear dynamics of an unknown system, this paper will introduce the NARMA model. It is a kind of model developed from a combination of neural networks and ARMA model. ANN (Artificial Neural Networks) as a good tool to solve nonlinear problem, has been used widely in diagnostics and surveillance problems in the area of nuclear engineering. ANN though has already been employed to identify the decay ratio by using the noise data (Van der Hagen et al., 1995), it can not reveal the nonlinear dynamics of the system because the ANN is only used to represent the nonlinear relationship between the decay ratio and the autocorrelation function. The ANN used in this paper can reflect the relationship between the input and output of a nonlinear system thus is able to represent the nonlinear dynamics of the identified system. Unlike the feedforward networks, it is in fact a kind of feedback networks because it will only take use of the reactor noise data series. Taking the superior ability of neural networks to approximate a nonlinear system, the NARMA model has been used successfully in nonlinear prediction (Bonnet, 1997). Here the NARMA model is simplified to do the identification. The ARMA model is also employed to testify the identification results with NARMA model. In this paper, the Vermont Yankee BWR model is chosen to generate the simulated noise data. Adjusting the void-temperature coefficient k, the noise data with different nonlinear dynamics can be obtained. When k=-0.00374 the bifurcation occurs in the system and a limit-cycle is shown. It is found that ARMA model will represent the dynamics of the system with an order higher than the intrinsic order of the system. As a result, many pseudo roots, which might conceal the important dynamics of the system, are derived. The NARMA model, however, could satisfy the model construction criteria in appropriate order. And the correspondent dynamics extracted from the linear kernel of NARMA model is consistent with the inherent dynamics of the Vermont Yankee Reactor. When the limit-cycle occurs in the noise data, the identified damping ratio turns out to be negative. It may be against the common sense because that it represents divergence in linear system, but is found to be a criterion for the existence of limit-cycle in a stable nonlinear system. In Section 2, the dynamics of the Vermont Yankee model is interpreted through the differential equations. In Section 3, the NARMA model will be introduced. The calculation results will be brought forth in Section 4. Section 5 is the conclusion part.

Method to identify nonlinear dynamics

211

2. THE NONLINEAR DYNAMICS IN BWR

2.1 Vermont Yankee BWR model To generate the simulated nonlinear noise data, a phenomenological model for Vermont Yankee BWR proposed by March-Leuba is used. The model is represented in the following equations (March-Leuba et al., 1986a): dn(t) _ p(t) - t3 n(t) + 2c(t) ÷ p ( t ) dt A A

(l)

de(t) _ 13 n ( t ) - 2c(t) dt A dT(t) -

dt

(2)

aln(t ) - a2T(t )

(13)

d~p~(t) dp~(t) dt 2 I-a 3 dt + a 4 P ~ ( t ) = k T ( t )

(4)

p ( t ) = p~ (t) + D T ( t )

(5)

Where n(t) = excess neutron density normalized to the steady-state neutron density c(t) = excess delayed neutron precursors concentration, also normalized to the steady-state neutron density T(t) = excess average fuel temperature 13 = delayed neutron fraction 2 = decay constant o f delayed neutron precursors A = neutron generation time p(t) = excess reactivity p~ (t) = excess void reactivity D = Doppler reactivity coefficient k, a~, a2, a 3

and

a 4 =

empirically determined constant parameters

The parameters take the following values: TABLE

Parameter Value

aI 25.04 Ks l

az 0.23 sl

1 -

a3 2.25 s"1

Parameters for BWR Model a4 -2.52x10 5 s 2

D 6.82 K l

2 0.08 sl

13 0.0056

A 4.00x10Ss

2.2 Analysis and simulation o f neutron noise It has been proved that k 0 = -0.00374 is the H o p f bifurcation point o f the system (Suzudo, 1993). When k < k0 , the limit-cycle occurs in the evolution o f the dynamics o f the system. To analyze its dynamics, we study the linearization part o f the equations.

212

Liu Wen/en£, et al. --j3/A

A

p/a A =

0 1/A

D/A

0

0

0

0

0

a,

0

-a 2

0

0

k

0

0

0

- a3 - a4 1

0

According to the concept in the linear system, some dynamics can be calculated when k takes different values as follows: TABLE 2 - Dynamics for BWR Model k

Time constant r I (sec.)

Time constant r 2 (see.)

Time constant r 3 (sec.)

Natural frequency (Hz)

Damping ratio

-0.0032

0.0071446

0.42752

13.357

2.6248

0.034219

-0.00374

0.0071446

0.39521

13.235

2.7058

-0.002241

Though the dynamic characteristics are calculated through the orthogonal decomposition of A thus they are independent, the original physical parameters can still be shown in them. The smallest time constant rlis mainly affected by the f l / A = 140/sec. The largest one is approximating to the half-time of the precursor of the delayed neutrons 12.5 sec. r2 is approximating to the temperature feedback time constant in Eq.(3). The oscillation dynamics is mainly determined by the void reactivity changes in Eq.(4). The neutron density fluctuation is simulated with temperature disturbance to the system as March-Leuba had done (March-Leuba et.al., 1986a). Taking the variance with the value of 0.01, external temperature disturbance is added as white noise to the system. When the sampling interval is taken as 0.25s and the whole sampling period 250s, the output neutron density is obtained for different k. To identify those characteristics such as time constant, natural frequency and damping ratio, NARMA model is presented in the following section.

3. NARMA MODEL Traditionally, AR or ARMA models are used to identify the characteristic of system with neutron noise. But they are based on the theory for linear system, and may not be qualified to identify the nonlinear noise here. Therefore, we will introduce the NARMA model, a kind of nonlinear model, to perform the identification. For the time series {X, }NI, The NARMA model takes as the following form: X , = f ( X , _ , , X,_2,...X,_p, c,q,,~,_2,...C,_q) + ~,

(6)

N can be viewed as the nonlinear response Where f is a nonlinear function. {c,} N, is white noise. {X , },=.

with the input of{c, }~l. To generalize the nonlinear approximating ability of NARMA model, f usually takes the form of ANN. It has been proved that a nonlinear multivariant function can be well approximated by a multi-layer networks (Poggio and Gerose, 1990). A simple form with three layers: one input layer, one hidden layer and one output layer can be used as follows: M

~..~.~,j

q

x, = Z w / , ( i=1

(7) j=l

k=l

where o-(.) is the Sigmoid function, M is the number of neurons in the hidden layer.

Method to identify nonlinear dynamics

2:13

Suppose the white noise input is separated from the nonlinear dynamics o f the system, therefore, {X, }i~_~ is the nonlinear response for any linear combination of white noise L 0 j ~ , _ j + c , . The model can be t=l

simplified in the following form: M

p

q

i-1

j~l

k-I

Given M, p and q, the system identification with the above model is a process of parameter estimation with series {X, }iv1 . The Nonlinear Least Square algorithm can be applied to the estimation of the parameters. To identify the dynamics of the BWR noise with NARMA model, a criterion for choosing model order should be setup. Because usually the smaller value the covariance is, the more appropriate model it is, here we simply use the Covariance Comparison Criterion for choosing order. From Eq.(8) the order of the NARMA model is determined by (p, q, M). q is fixed as p - 1, thus a 2-dimensional search for (p, M) should be made to select the model order. The order-searching process is as follows: (i) For any p, calculate the covariance o f the NARMA(p, q, M) model denoted by cov(p, q, M) from M=2p-1 to any number when]l-cov(p, q, M)/cov(p, q, M-1)]
4. IDENTIFICATION RESULTS

4.1 Failure o f ARMA identification The ARMA(p, p - l ) model is applied to identify the simulated noise data for k = -0.0032 and k -- -0.00374 respectively. For the linear system identification, the necessary condition is that the residual should be white. The Chi-square test with a significant level at 0.05 is used to check the whiteness for the residual. Some identification results are obtained as follows: TABLE 3 - Identification Performance with ARMA Model 3

4

5

6

7

8

9

Whiteness check N N Y Y Y Y Y (k = -0.0032) Whiteness check N N N N N N Y (k = -0.00374) Residual covariance 0.008675 0.008548 0.008024 0.008006 0.007976 0.007906 0.007900 (k = -0.0032) Residual covariance 0.1175 0.082415 0.078521 0.076922 0.063764 0.062751 0.054658 (k = -0.00374) * N represents the residual fails to pass the whiteness check while Y represents a pass. From the TABLE 2, the order of this system should be 5, but because the sampling interval is 0.25s, the dynamics r~ and r 3 can not be detected in the simulated noise data, the real order of the system should be 3. But to pass the whiteness check for the residual, for k = -0.0032, the order o f ARMA should at least be 5 and 9 for k = - 0 . 0 0 3 7 4 . Thus many other pseudo characteristics are obtained. When oscillation occurs in the reactor system because of the limit cycle, the ARMA model failed to identify the dynamics

Liu Wenfeng et al.

214 4.2 Identification with NARMA model

The identification results with NARMA model and ARMA model are listed comparatively in the following table: TABLE 4 - Identification Results with ARMA and NARMA Model k : -0.0032

k : -0.00374

NARMA ARMA NARMA ARMA Model order (3, 2, 6) (5, 4) (3, 2, 7) (9, 8) Whiteness check Y Y Y Y Residual covariance 0.002650 0.008024 0.000965 0.054658 The identification results with NARMA are much better than ARMA in the sense of the covariance and whiteness check of the residuals, especially when the limit-cycle tends to occur in the system i.e. k = -0.00374. To analyze the dynamics of the system identified by NARMA model, the autoregressive part in the model is firstly expanded as a Taylor series on the bias of different neurons. Then the expanded series were combined in one series as follows: M

S

,--o ,:1

t-j) n

n!

(9)

The first order part in Eq.(9) called linear kemel is used to generate the characteristics polynomial. Thus p is still the determinant factor of the order of the system and the pseudo dynamics in the ARMA model can be avoided in the NARMA model. From the characteristics polynomial, the eigenvalues can be calculated and corresponding dynamics such as time constant, natural frequency and damping ratio can be solved (Pandit and Wu, 1984). In the following table, the dynamics for different k is listed. TABLE 5 - Identified Dynamics k = -0.0032 NARMA

Vermont Yankee model 0.42752

k = -0.00374 NARMA

Vermont Yankee model 0.39521

0.35854 0.6115 Time constant (sec.) Natural frequency (Hz) 3.0816 2.6248 3.3207 2.7058 0.022516 0.034219 -0.031599 -0.002241 Damping ratio It is interesting to find that the negative damping ratio can be calculated by the NARMA model for k = 0.00374, which is similar to the linearization part of the Eq.(1)-Eq.(5). It's reasonable to guess that in such a BWR system, the negative damping ratio is the criterion for the existence of limit-cycle. For a p-Dimensional discrete system, the impulse response on the sub phase space (Xt, Xt_ ~) can reflect some characteristics. Input an impulse whose amplitude is the square root of the residual covariance, the response of NARMA model is plotted on the (Xt, X,_~)plane as shown in Fig. 1. From Fig. I.(B), the limit-cycle can clearly be observed, which can be distinguished from Fig. I.(A) where the response of the system tends to zero with the increase of time. Because the negative damping ratio can't be obtained for any linear system, the advantages of the NARMA model over ARMA model are obvious.

Method to identify nonlinear dynamics

215

O.06r

"

0 04.

1

,,s.

I

1i

002 i

.

/

oi

f*-

ost -002

i * I

"i

/

°i

-0 06.

!

'

I

*]

-004

"k

t

/ J

*

,

-0.08 -0 1

-012

12

471

-008

-006

<)0,4

4J02

0

002

004

-1

006

4]5

0

O5

1

(B). Impulse response for k = -0.00374

(A). Impulse response for k = -0.0032 FIGURE 1

5. CONCLUSION This paper presents a kind of NARMA model in the dynamics identification with nuclear reactor noise data. Comparing with the traditional ARMA model, the NARMA model has the following advantages:(i) Avoid generating the pseudo dynamics by identifying the system with a low order model. (ii) Give the criterion for the existence of limit-cycle qualitatively. The current analysis is still applied to the first order part of the Taylor expansion of the NARMA model. The analysis for its high order might reveal more information about the system. And the further work will focus on the higher-order analysis and studying its application area.

REFERENCES Upadhyaya B. R. and Kitamura M. (1981), Stability Monitoring of Boiling Water Reactors by Time Series Analysis of Neutron Noise. Nucl. Sci. & Eng. 4, 480. Bonnet D., Labouisse V. and Grumbach A. (1997), 6-NARMA Neural NetWorks: A New Approach to Signal Prediction. 1EEE Transactions on Signal Processing 45, 28 78. Van Bragt D.D.B., Rizwan-uddin. and Van der Hagen T. H. J. J. (2000), Effect of Void Distribution Parameter and Axial Power Profile on Boiling Water Reactor Bifurcation Characteristics. Nucl. Sci. & Eng 134, 227. Gialdi E., Grifoni S., Parmegginai C. and Tricoli C. (1985), Core Stability in Operating BWR: Operational Experience. Prog. in Nucl. Energy 15, 447. March-Leuba J., Cacuci D. G. and Perez R. B. (1986a), Nonlinear Dynamics and Stability of Boiling Water Reactors: Partl -Qualitative Analysis. Nucl. Sci & Eng 93, 111. March-Leuba J., Cacuci D. G. and Perez R. B. (1986b), Nonlinear Dynamics and Stability of Boiling Water Reactors: Part 2-Quantititive Analysis. Nucl. Sci & Eng 93, 124 Rizwan-uddin and Doming J. J. (1986), Some Nonlinear Dynamics of a Heated Channel. Nuclear Engineering and Design 93, 1. Pandit S. M. and Wu S. M. (1984), Time Series and System Analysis with Application, John Wiley & Sons Narayanan S., Srinivas B., Pushpavanarn S. and Murty Bhallamudi S. (1997), Non-linear Dynamics of a Two Phase Flow System in an Evaporator: the Effects of (i) a Time Varying Pressure Drop (ii) an Axially Varying Heat Flux. Nuclear Engineering and Design 178, 294.

216

Liu Wenfeng et al.

Van der Hagen T. H. J. J. (1995), Artificial Neural Networks versus Conventional Methods for Boiling Water Reactor Stability Monitoring. Nucl. Technology 109, 286. Suzudo T. (1993), Reactor Noise Analysis Based on Nonlinear Dynamics Theory-Application to Power Oscillation. Nucl. Sci. & Eng 113, 145. Poggio T., Gerose F. (1990), Network for Approximation and Learning. Proc. IEEE 78, 1481. Masashi T., Katsuhisa N. and Masakuni N. (1993), Stability Analysis of BWRs Using Bifurcation Theory. Journal of Nucl. Sci. and Tech 30, 1107. Luo Z. P., Li F. and Li B. Y. (1996), Some New Viewpoints in Reactor Noise Analysis. Nucl. Sci. & Tech 7, 202.