- Email: [email protected]
Application of autoregressive moving average model in reactor noise analysis
Inn. Nucl. Energy, Vol. 20, No. 12, pp. 815~822, 1993
0306-4549/93 $6.00+0.00 Pergamon Press Ltd
)rinted in Great Britain
APPLICATION
OF AUTOREGRESSIVE MOVING AVERAGE IN REACTOR NOISE ANALYSIS
MODEL
TRAN DINH TRI CRIP Institute for Atomic Energy Research Applied Reactor Physics Department H-1525 Budapest 114, P.O.B. 49, Hungary
(Received 11 March 1993) A b s t r a c t - - In this paper the system of the generalized Yule-Walker equations is derived from the equation of an ARMA model, then a method for its solution is given. Numerical results show the applications of the method proposed.
1. INTRODUCTION The application of the AR model to estimating noise measurements has achieved many successes in reactor noise analysis in the last ten years. The physical processes that take place in the nuclear reactor, however, are described, as shown by Kishida (1982, 1984), by an ARMA model rather then by an AR model. Consequently more correct results could be received by applying the ARMA model insted of the AR model. An ARMA model with innovation can be written as n
p
y(t) = Z A(i)y(t - i) + 7(t) + Z B(i)7(t - i) , i=1
(1)
i=1
where y(t) is a q-dimensional observable vector, A(i) and B(i) are parameters of the ARMA model, n and p axe the order of the AR and MA parts respectively, 7(t) is the innovation which is in fact (Tran Dinh Tri, 1992a) the residual noise of an AR model: NAR
7(t) = y(t) - Z
M ( j ) y ( t - j) ,
(2)
j=l
where M ( j ) are parameters of an AR model with an order NAR NAIl
M(z)C(k) = [I- Z
M(J)z-J]C(k) = 0 ;
k = 1,2,...,NAR ,
(3)
j=l
I is the q-dimensional unit matrix, C(k) are correlation matrices: C(k) = < y(t)yT(t--k) >, C ( - k ) = c T ( k ) and z -i is the time shift operator: z - i C ( k ) = C(k - i). To estimate the parameters A(i) and B(i) we can use the least-squares method by taking the ensemble average < . - . . y T ( t -- k) > of both sides of Eq. (1). This brings about the following Yule-Walker equation (Tran Dinh Tri, 1992a): NAR
[A(z)- B(z)(I-
Z
M(J)z-J)]C(k) = 0 ;
j=l
815
k = 1,2,...,n ,
(4)
816
TRAN DINHTRI where
n
p
A(z) = I - ~ A(Oz-' ;
B(z) = X + ~ B(i)z-' .
i=1
(5)
i=1
If the operator function B(z) in Eq. (4) is equal to the unit matrix, that corresponds to the P
case of the AR model when the MA part, i.e. the term ~ B(i)7(t - i), is truncated from i=1
the Eq. (1), then taking into account Eq. (3) and assuming that n <_ NAn we can see that Eq. (4) becomes the conventional Yule-Walker equation of the AR model. Therefore, Eq. (4) is called the generalized Yule-Walker equation for the ARMA model. Note that the generalized Yule-Walker equation (4) corresponds to the least-squares condition only for A(i) but not for B(i). In section 2 it will be shown that the least-squares condition for B(i) results in another Yule-Walker equation that has the following form:
A(z)[M(z)C(-k)] T - B ( k ) r = 0 ;
k = 1 , 2 , . . . ,p ,
(6)
where P = M(z)C(O), B(k) is the k-th matrix of the operator function B ( z ) in (5). In a preceding paper (Tran Dinh Tri, 1992a) a recursive algorithm for A(i) and B(i) satisfying the generalized Yule-Walker Eq. (4) for k = 1, 2 , . . . , n has been proposed . In section 3 we will give a method for the solution .4,,,p(i) and Bn,n(i) of the system of Eqs. (4) and (6) based on A(i) and B(i) computed by the recursive algorithm mentioned above. In section 4 calculational results will be presented concerning the application of the ARMA model to evaluating noise measurements in various nuclear power reactors. 2. SYSTEM OF THE YULE-WALKER EQUATIONS FOR A(i) AND B(i) The least-squares method can be used for estimating parameters A(i) and B(i) of the ARMA model (1). For this, let us consider t h e / - t h residual time-series which follows from t h e / - t h component of y(t) in the ARMA model (1): n
q
p
q
;
7,(t)=yt(t)-~atj(i)yj(t-i)-~Btj(i)Tj(t-i) i=l j=l
l=l,2,...,q
,
(7)
i=1 j=l
where t is the discrete time t = 1, 2 , . . . , Nt and Nt is the total number of samples. The least-squares estimates of Atj (i) and Btj (i) will be obtained by minimizing the quantity u, = y ~ ~,~(t) = ~--~[y,(t) t=l
t=l
~
A,i(i)yj(*
- i)-
Btj(i)~,~(t
i=1 j=l
- i)l 2 ,
i = 1 ./=1
l= 1,2,...,q .
(S)
This requirement leads to the following conditions:
N, ~ Out -- --2 ~--~[yl(t)OAlm(k) t=l i=1
pq Atj(i)yj(t - i) - ~ ~
- -
j=l
I,m=l,2,...,q
Bij(i)Tj(t - i)]ym(t - k) = 0 ,
i=1 j=l
;
k=l,2,...,n
,
(9.a)
Use of the ARMA model in reactor noise analysis
817
and -
Out -
A,j(il~j(t - i) - ~ ~
- -2 ~[~,(t)t----1
i = 1 j----1
e,~(i)~(t
- ilbm(t
- k) = 0 ,
i=1 j = l
l,m=l,2,...,q
;
k=l,2,...,p
.
(9.b)
It is easy to see that Eqs. (9.a) and (9.b) can be rewritten in the matrix form as
A(i)y(t-i)-~-~B(i)~/(t-i)lyT(t-k)>=O
; k=l,2,...,n
(10.a)
; k=l,2,...,p
(10.b)
i=1
and p
A(i)y(t--i)--EB(i)7(t--i)lTT(t-k)>=O
i=1
Since < 7 ( t ) y T ( t - - k ) > = M ( z ) C ( k ) = 0 for k = 1 , 2 , . . . , N A n , obviously Eq. (10.a) can be written as < a ( z ) y ( t ) - B ( z ) 7 ( t ) l y T ( t -- k) > = 0 . The generalized Yule-Walker equation (4) follows from the above equation (Tran Dinh Tri, 1992a), namely, it is obtained from Eq. (10.a). Substituting 7T(t - - k ) a s determined by Eq. (2) into Eq. (10.b) and using the following expressions: NAR
< y(t - i)77"(t - k) > = C ( k - i) - E
C(k - i +j)MT(j)
= [M(z)C(-k + i)]T(ll)
j=l
< 7(t - i)TT(t - k) > = F6ik ,
(12)
yield the result NAR
n
C ( k + j ) M T ( j ) ] -- E
[C(k) - E j=l
i=1
NAR
A(i)[C(k - i ) - E
C(k - i + j ) M T ( j ) ] --
B(k)r =
0
j=l
k = 1,2,...,p.
(13)
Since NA R
C ( k - i) - E
C ( k - i + j ) M T ( j ) = [ M ( z ) C ( - k + i)] T = z - i [ M ( z ) C ( - k ) ] T ,
(14)
j=l
we can rewrite Eq. (13) as
[ M ( z ) C ( - k ) ] T - ~ - ~ A ( i ) z - i [ M ( z ) C ( - k ) ] T - B ( k ) r = 0; i=l
k = 1 , 2 , . . . ,p .
(15)
818
TRAN DINH TRI
Taking into account the definition (5) of the operator function A(z), from Eq. (15) we will get another Yule-Walker equation (6) that corresponds to the least-squares condition for the matrices B(i). Furthermore, using Eq.(3) we can now rewrite the generalized Yule-Walker equation (4) in the following form:
I B,,,p(1)M(z)C(O) + Bn,p(2)M(z)C(-1) + ... + B,,,p(p)M(z)C(1 - p) for k = l , B,,,p(2)M(z)C(O) + ... + B,,,p(p)M(z)C(2 - p) A,,,p(z)C(k) = for k = 2 , B,,,v(p)M(z)C(O) for k = p , ifk>p,
0
(16)
and the other Yule-Walker equation (6) as
[M(z)C(-1)] T - A,, v(1)F - B , p(1)r = 0 [M(z)C(-2)] T {A:,p(1)[M(z)C(-1)] T +
[M(z)C(-p)] T
for k = 1 , d.,p(2)r}
{A.,,p(1)[M(z)C(-p + 1)] T
+...
- B.,p(2)r
+ A.,p(p)r}
= 0 Ior - B.,p(p)r
k = 2 , = 0
for k = p , (17)
here the relation F T = F was exploited. In the next section we shall show a method for determining the solution of the system of the Yule-Walker equations (16) and (17) based on the solution of Eq. (4) for k = 1, 2 , . . . , n computed by the recursive algorithm described in the preceding paper (Trail Dinh Tri, 1992a). 3. SOLUTION OF THE SYSTEM OF THE YULE-WALKER EQUATIONS A solution satisfying Eq. (4) or (16) for k = 1, 2 . . . . ,n and also satisfying Eq. (6) or (17) for k = 1 , 2 , . . . ,p will be denoted by -4n,p(i) and/)n,p(i), and the corresponding operator functions by .4n,p(z) and/),,,v(z), while the notations A,,,v(i), B,,,p(i) and A,,,v(z), B,,,v(z) now will used to indicate the solution of Eq. (4) which satisfies this equation only for k = 1, 2 , . . . , n and is computed by the recursive algorithm mentioned above: The main differences between -4n,p(i),/)n,v(i) and An,v(i), B,,,v(i) are the followings: + First, Ji,,,v(i) and/),,,v(/) are the least-squares estimates of the ARMA model parameters, that are estimated not only from the least-squares condition (9.a) for the AR part but also from the least-squares condition (9.b) for the MA part, while An,v(/) and B,,,v(i) are calculated only from the least-squares condition (9.a) for the AR part. + Second, since An,p(i) and B,,,v(i) not satisfy the least-squares condition (9.b) for the MA part, they are not the least-squares estimates of the ARMA model parameters. However, the use of these A,,,v(i) and Bn,v(i) allows us to determine the least-squares estimates of the ARMA model parameters as shown below.
Use of the ARMA model in reactor noise analysis
819
Let us determine A.,p(z) and/},.,v(z) first for p = 1 in the following form: An,l(Z ) = aoAn_l,l(z ) + ffl-'An,l(Z)Z -n ~- a2an,l(Z )
*}.,1(z) = 0-oB._l,,(z) + 0-~-9.,,(z)z -~ + 0-~B.,~(z) ,
(18)
where An,1 (z) and Bn,l (z) are adjoint operator functions that satisfy the adjoint equation NAR [An'l(Z)
-
-Bn'l(Z)zn-l(I - Z
MJz-J)]C(k)
= O ,
for
k = -1,-2,...,-n
(19)
,
j=l
and 0-0, 0-1, 0"2 are constant matrices to be found. These matrices can be determined by substituting expressions (18) into Eq. (4) for k = n and into Eq. (6) for k = 1, then the following equations are obtained: 0-0Wn-l,l(n) + O'l~n,l(0)
= 0 w I 0-0w'_,,,(1) + 0 - ~ . , 1 ( 1 ) + 0-2 .,~(1) = 0 .
(20)
Here we use the following definitions: w n , v ( k ) = [an,v(z) - B , , v ( z ) M ( z ) ] C ( k ) W.,l(0) = [~.,1(z) - - 9 . , , ( z ) z " - I M ( z ) ] C ( O )
(21.a)
(21.b)
tbn,p(k) = [-4n,p(z) - 1 3 , , p ( z ) M ( z ) ] C ( k )
(21.c)
w',,v(k) = A n , v ( z ) [ M ( z ) C ( - k ) ]
(21.d)
~.,1(1)
T - B,,p(k)F
= --~.,l(n)[M(z)C(-1)]
~ - ].,~(n
ffo',v(k) = A , , p ( z ) [ M ( z ) C ( - k ) ]
T -/3,,p(k)F
- 1)r - r
.
(21.~)
(21.f)
It is seen from Eq. (20) that 0-i = aovi (i -= 1, 2) with V1 = - - W n _ l , l ( n ) w n l l ( 0 ) /]2 = --[Wtn--l,l(1) "~ /]lW-'~n,l(1)]Wtn.ll(I) •
(22)
Substituting these ui into the condition I = ao - 0 - , A n , l ( n ) + 0-2 = ao + a,B,,,a(1) + a2 ,
(23)
and taking into account that - A . , a ( n ) = Bn,l(1) we will get 0-0 as
oo = (I + . , ~ . , , ( 1 )
+ ~2) -~ .
(24)
Thus -4n,l(z) a n d / ~ n , l ( z ) in Eq. (18) can be computed with the aid of a0, ul and u2: .Z[n,l(Z ) = 0-o[dn_l,l(Z) + vl-An,l(Z)Z -n + v2An,l(Z)] Bn,l(Z) = 0-o[Bn-l,l( z)'~ b'lBn,l(Z) Z-1 + tl2Bn,l(Z)] •
(25)
820
TRAN DINH TRI
Let us go on to determine An,p(z) and/~n,p(z) for p = 2, 3,... in the following form: an,p(z ) = otoAn,l(Z ) + o Q A n , I ( Z ) + ... + otp-lan,p-l(Z) + apAn,p(z) JBn,p(Z) = otoBn,l(Z ) + o t l n n , l ( Z ) + ... + olp_lnn,p_l(Z ) + OlpBn,p(Z ) ,
(26)
here o~0, 31,..., a v are constant matrices that can be determined from the condition (27)
Or0 "31-O~1 "~ ... "~- O/p = I ,
and from the following equations received by putting Eq. (26) into Eq. (6) for k = 1, 2, ..., p: = -avwln,v(1) = -apw',,(2)
{ a0w~,l(1 )
a0w',x(2) + a l e " 1(2) a0w'.,l(3) +
+
=
Ot0Wtn,l (P) + OtlWtn,l(p) + Ot2Wln,2(p) + ... + Otp--l~3tn,p--l(p)
(28)
= --OtpWtn,p(p)
,
where w;,v(k ) mad Vo~,v(k ) are given in Eqs. (21.d) mad (21.f), respectively. To determine the matrices at, let at = avTt with l = 0, 1, 2,... , p - 1 then by substituting them into Eqs. (28) we obtain the following Gaussima linear algebraic equations: { 70w~,1(1 ) 7OWln,l(2) + 71Win,l(2) '70Win,l(3) + "71Wtn,l(3) + 72Wtn,2(3)
= --w~n,v(1) = --Wtn,p(2) = --w:,p(3)
70Wn,I(P) + "71Wn,l(P) + "72 n 2(P) + "'" + 7p-lWn,p-l(P)
(29)
= --Wn,p(P)
Obviously Eqs. (29) now could be solved directly from the first equation for 70, then to the second one for 71 with the given 70, etc. The final form of the solution of Eqs. (29) can be expressed as t 1-1 "70 = -w,,v(1)wn,1 (1) "71 = -- ['70Wtn,l( 2 ) + Wn,p(2)]Wn, t ~t--11 (2)
t -t-1 (3) 5"2 = --["7ow'0(3) + "71W',1(3) + w,,,p(3)]w.,2 I
~1
l
~1-1
"Tp--1 = --["70Wn,I(P) "{- 71Wln,I(P) + .-. "~- "Tp--2Wn,p--2(P ) "~ Wn,p(p)]Wn,p_l(p)
,
(30)
and the condition (27) gives a p = [70 31- 71 ~- 72 3t- "." "~- "7p-1 Ji- I ] - 1
.
(31)
Finaiy, the operator functions .4,,p(z) and Bn,p(z) in (26) will be computed by the expressions
A,,,p(z) = ap["7oA,,l(z) + 71/i,,a(z) + ... + 7,-lfi--,,,p-x(z) + Amp(z)] JBn,p(z) = Olp[7oSn,l(Z) + "71Bn,1(z) "gt--.. 3I- 7 p - l n n , p - l ( Z ) + Bn,p(Z)] .
(32)
Use of the ARMA model in reactor noise analysis 4. S O M E RESULTS O F T H E A P P L I C A T I O N W i t h the aim of illustration in this section we show some results of the evaluation of the noise measurements carried out at various NPPs. The numerical results were obtained by using the program package named A R M A A N A , which was built up by the author for evaluating the P W R and B W R reactor noise diagnostic measurement d a t a based on the recursive algorithm (Tran Dinh Tri, 1992a), the method developed in the previous section and the generalized S T P method (Tran Dinh Tri, 1992b). The measurement data mentioned above were processed into two d a t a files named " T I M . . . . . l A W ' and " T I M ...... 1BR". Some features of the two d a t a files are given in the following table: name of data file " T I M ...... 1AU . . . . T I M ...... 1BR number of signals 64 16 sampling interval A t , sec 0.01 0.01 fma~ = 2-~ ' H z 50 50 considered detector combinations: first signal pressure Q3 inlet pressure PB second signal vibration detector V3 neutron detector A7 third signal neutron detector C3 ionisation chamber I2 fourth signal ionisation chamber I1 thermocouple T3
A m o n g basic input parameters there are some important ones: - number of points per block
: NPOINT
- number of blocks number of frequency points
: N B L O C K = 10 : N F = 160
-
= 1024
To select the optimal order of the A R model, a certain information criterion, for example the Akaike's Information Criterion (AIC), can be applied (Fukunishi, 1977). For the given detector combinations of the data sets shown in the above table, the A I C ( n ) function for the A R model was computed for model order up to n -- 20. However, any m i n i m u m of the A I C ( n ) function up to n = 20 is not found in the case of the first combination of the " T I M . . . . . 1AU" data set, and in the case of the second combination of the " T I M ...... 1BR" data set, a minimum of the A I C ( n ) function is found at n = 19. For the both case of the given detector combinations the A I C ( n , p) funtion for the A R M A model however still were not considered. An examination of the optimal order of the A R M A for paractical noise data sets should be done in a separate paper. Therefore, in order to present some calculatinal results of the A R M A model abitraz7 indices n = 14, p = 14 and N A n = 20 are taken without any consideration of the optimal problem. So these indices are not any optimal model order. The ARMA(14,14) model parameters fi,,,p(i) and/~n,p(i) that are calculated by the m e t h o d given in section 3, are used to compute the spectra and coherences. T h e calculation results of autospectra, ordinary and partial coherences for the first d a t a set " T I M ...... •AU" are given in figs. 2, 3, 5 and 6, while the corresponding results for the second d a t a set
821
822
TRAN DINHTm " T I M . . . . . 1BR" are shown in figs. 8, 9, 11 and 12. The F F T autospectra and coherences of both of the given data sets also are presented in figs. 1, 4, 7 and 10 for comparision. Although the used indices (n = 14, p = 14 and WAR 20) are not any optimal model order, based on the received results of numerical calculations we can nevertheless conclude that in the case of the ARMA model the shape of the autospectral and coherence functions resembles more to the F F T spectrum and coherence than in the AR case, and this conclusion is valid even in the case when the order of AR part of the ARMA model is lower than that of the AR model. For example, the autospectral (APSD) and coherence (COH) functions in the ARMA(14,14) model of both data sets are better than that of the AR(20) model. For the partial quantities, similar conclusions can be drawn. =
ACKNOWLEDGMENTS The author would like to thank Dr I. Lux for his help. This work was partly supported by Hungarian Government Funds OTKA 1845 and OMFB 91-97-42-0344.
REFERENCES Fhkunishi K. (1977) Nucl. Sei. and Eng., 62, pp. 215-225. Kishida, K. (1982) Physical Langevin model and the time-series model in systems far from equilibrium. Physical Review A, Vol. 25, No. 1, pp. 496-507. Kishida, K. (1984) Equivalent random force and time-series model in system far from equilibrium. J. Math. Phys., Vol. 25, No. 5, pp. 1308-1313. Tran Dinh Tri (1990) Determination of coefficient matrices for ARMA model. KFKI Report, 1 9 9 0 - 4 5 / G , Budapest, Hungary. Tran Dinh Tri (1992a) A new algorithm for recursive estimation of ARMA parameters in reactor noise analysis. Annals of Nuclear Energy, Vol. 19, pp. 287-301. Tran Dinh Tri (1992b) Generalization of the signal transmission path method for ARMA model in reactor noise analysis. Annal, of Nuclear Energy, Vol. 19, pp. 341-345.
0306-4549/93 $6.00+0.00 Pergamon Press Ltd
)rinted in Great Britain
APPLICATION
OF AUTOREGRESSIVE MOVING AVERAGE IN REACTOR NOISE ANALYSIS
MODEL
TRAN DINH TRI CRIP Institute for Atomic Energy Research Applied Reactor Physics Department H-1525 Budapest 114, P.O.B. 49, Hungary
(Received 11 March 1993) A b s t r a c t - - In this paper the system of the generalized Yule-Walker equations is derived from the equation of an ARMA model, then a method for its solution is given. Numerical results show the applications of the method proposed.
1. INTRODUCTION The application of the AR model to estimating noise measurements has achieved many successes in reactor noise analysis in the last ten years. The physical processes that take place in the nuclear reactor, however, are described, as shown by Kishida (1982, 1984), by an ARMA model rather then by an AR model. Consequently more correct results could be received by applying the ARMA model insted of the AR model. An ARMA model with innovation can be written as n
p
y(t) = Z A(i)y(t - i) + 7(t) + Z B(i)7(t - i) , i=1
(1)
i=1
where y(t) is a q-dimensional observable vector, A(i) and B(i) are parameters of the ARMA model, n and p axe the order of the AR and MA parts respectively, 7(t) is the innovation which is in fact (Tran Dinh Tri, 1992a) the residual noise of an AR model: NAR
7(t) = y(t) - Z
M ( j ) y ( t - j) ,
(2)
j=l
where M ( j ) are parameters of an AR model with an order NAR NAIl
M(z)C(k) = [I- Z
M(J)z-J]C(k) = 0 ;
k = 1,2,...,NAR ,
(3)
j=l
I is the q-dimensional unit matrix, C(k) are correlation matrices: C(k) = < y(t)yT(t--k) >, C ( - k ) = c T ( k ) and z -i is the time shift operator: z - i C ( k ) = C(k - i). To estimate the parameters A(i) and B(i) we can use the least-squares method by taking the ensemble average < . - . . y T ( t -- k) > of both sides of Eq. (1). This brings about the following Yule-Walker equation (Tran Dinh Tri, 1992a): NAR
[A(z)- B(z)(I-
Z
M(J)z-J)]C(k) = 0 ;
j=l
815
k = 1,2,...,n ,
(4)
816
TRAN DINHTRI where
n
p
A(z) = I - ~ A(Oz-' ;
B(z) = X + ~ B(i)z-' .
i=1
(5)
i=1
If the operator function B(z) in Eq. (4) is equal to the unit matrix, that corresponds to the P
case of the AR model when the MA part, i.e. the term ~ B(i)7(t - i), is truncated from i=1
the Eq. (1), then taking into account Eq. (3) and assuming that n <_ NAn we can see that Eq. (4) becomes the conventional Yule-Walker equation of the AR model. Therefore, Eq. (4) is called the generalized Yule-Walker equation for the ARMA model. Note that the generalized Yule-Walker equation (4) corresponds to the least-squares condition only for A(i) but not for B(i). In section 2 it will be shown that the least-squares condition for B(i) results in another Yule-Walker equation that has the following form:
A(z)[M(z)C(-k)] T - B ( k ) r = 0 ;
k = 1 , 2 , . . . ,p ,
(6)
where P = M(z)C(O), B(k) is the k-th matrix of the operator function B ( z ) in (5). In a preceding paper (Tran Dinh Tri, 1992a) a recursive algorithm for A(i) and B(i) satisfying the generalized Yule-Walker Eq. (4) for k = 1, 2 , . . . , n has been proposed . In section 3 we will give a method for the solution .4,,,p(i) and Bn,n(i) of the system of Eqs. (4) and (6) based on A(i) and B(i) computed by the recursive algorithm mentioned above. In section 4 calculational results will be presented concerning the application of the ARMA model to evaluating noise measurements in various nuclear power reactors. 2. SYSTEM OF THE YULE-WALKER EQUATIONS FOR A(i) AND B(i) The least-squares method can be used for estimating parameters A(i) and B(i) of the ARMA model (1). For this, let us consider t h e / - t h residual time-series which follows from t h e / - t h component of y(t) in the ARMA model (1): n
q
p
q
;
7,(t)=yt(t)-~atj(i)yj(t-i)-~Btj(i)Tj(t-i) i=l j=l
l=l,2,...,q
,
(7)
i=1 j=l
where t is the discrete time t = 1, 2 , . . . , Nt and Nt is the total number of samples. The least-squares estimates of Atj (i) and Btj (i) will be obtained by minimizing the quantity u, = y ~ ~,~(t) = ~--~[y,(t) t=l
t=l
~
A,i(i)yj(*
- i)-
Btj(i)~,~(t
i=1 j=l
- i)l 2 ,
i = 1 ./=1
l= 1,2,...,q .
(S)
This requirement leads to the following conditions:
N, ~ Out -- --2 ~--~[yl(t)OAlm(k) t=l i=1
pq Atj(i)yj(t - i) - ~ ~
- -
j=l
I,m=l,2,...,q
Bij(i)Tj(t - i)]ym(t - k) = 0 ,
i=1 j=l
;
k=l,2,...,n
,
(9.a)
Use of the ARMA model in reactor noise analysis
817
and -
Out -
A,j(il~j(t - i) - ~ ~
- -2 ~[~,(t)t----1
i = 1 j----1
e,~(i)~(t
- ilbm(t
- k) = 0 ,
i=1 j = l
l,m=l,2,...,q
;
k=l,2,...,p
.
(9.b)
It is easy to see that Eqs. (9.a) and (9.b) can be rewritten in the matrix form as
A(i)y(t-i)-~-~B(i)~/(t-i)lyT(t-k)>=O
; k=l,2,...,n
(10.a)
; k=l,2,...,p
(10.b)
i=1
and p
A(i)y(t--i)--EB(i)7(t--i)lTT(t-k)>=O
i=1
Since < 7 ( t ) y T ( t - - k ) > = M ( z ) C ( k ) = 0 for k = 1 , 2 , . . . , N A n , obviously Eq. (10.a) can be written as < a ( z ) y ( t ) - B ( z ) 7 ( t ) l y T ( t -- k) > = 0 . The generalized Yule-Walker equation (4) follows from the above equation (Tran Dinh Tri, 1992a), namely, it is obtained from Eq. (10.a). Substituting 7T(t - - k ) a s determined by Eq. (2) into Eq. (10.b) and using the following expressions: NAR
< y(t - i)77"(t - k) > = C ( k - i) - E
C(k - i +j)MT(j)
= [M(z)C(-k + i)]T(ll)
j=l
< 7(t - i)TT(t - k) > = F6ik ,
(12)
yield the result NAR
n
C ( k + j ) M T ( j ) ] -- E
[C(k) - E j=l
i=1
NAR
A(i)[C(k - i ) - E
C(k - i + j ) M T ( j ) ] --
B(k)r =
0
j=l
k = 1,2,...,p.
(13)
Since NA R
C ( k - i) - E
C ( k - i + j ) M T ( j ) = [ M ( z ) C ( - k + i)] T = z - i [ M ( z ) C ( - k ) ] T ,
(14)
j=l
we can rewrite Eq. (13) as
[ M ( z ) C ( - k ) ] T - ~ - ~ A ( i ) z - i [ M ( z ) C ( - k ) ] T - B ( k ) r = 0; i=l
k = 1 , 2 , . . . ,p .
(15)
818
TRAN DINH TRI
Taking into account the definition (5) of the operator function A(z), from Eq. (15) we will get another Yule-Walker equation (6) that corresponds to the least-squares condition for the matrices B(i). Furthermore, using Eq.(3) we can now rewrite the generalized Yule-Walker equation (4) in the following form:
I B,,,p(1)M(z)C(O) + Bn,p(2)M(z)C(-1) + ... + B,,,p(p)M(z)C(1 - p) for k = l , B,,,p(2)M(z)C(O) + ... + B,,,p(p)M(z)C(2 - p) A,,,p(z)C(k) = for k = 2 , B,,,v(p)M(z)C(O) for k = p , ifk>p,
0
(16)
and the other Yule-Walker equation (6) as
[M(z)C(-1)] T - A,, v(1)F - B , p(1)r = 0 [M(z)C(-2)] T {A:,p(1)[M(z)C(-1)] T +
[M(z)C(-p)] T
for k = 1 , d.,p(2)r}
{A.,,p(1)[M(z)C(-p + 1)] T
+...
- B.,p(2)r
+ A.,p(p)r}
= 0 Ior - B.,p(p)r
k = 2 , = 0
for k = p , (17)
here the relation F T = F was exploited. In the next section we shall show a method for determining the solution of the system of the Yule-Walker equations (16) and (17) based on the solution of Eq. (4) for k = 1, 2 , . . . , n computed by the recursive algorithm described in the preceding paper (Trail Dinh Tri, 1992a). 3. SOLUTION OF THE SYSTEM OF THE YULE-WALKER EQUATIONS A solution satisfying Eq. (4) or (16) for k = 1, 2 . . . . ,n and also satisfying Eq. (6) or (17) for k = 1 , 2 , . . . ,p will be denoted by -4n,p(i) and/)n,p(i), and the corresponding operator functions by .4n,p(z) and/),,,v(z), while the notations A,,,v(i), B,,,p(i) and A,,,v(z), B,,,v(z) now will used to indicate the solution of Eq. (4) which satisfies this equation only for k = 1, 2 , . . . , n and is computed by the recursive algorithm mentioned above: The main differences between -4n,p(i),/)n,v(i) and An,v(i), B,,,v(i) are the followings: + First, Ji,,,v(i) and/),,,v(/) are the least-squares estimates of the ARMA model parameters, that are estimated not only from the least-squares condition (9.a) for the AR part but also from the least-squares condition (9.b) for the MA part, while An,v(/) and B,,,v(i) are calculated only from the least-squares condition (9.a) for the AR part. + Second, since An,p(i) and B,,,v(i) not satisfy the least-squares condition (9.b) for the MA part, they are not the least-squares estimates of the ARMA model parameters. However, the use of these A,,,v(i) and Bn,v(i) allows us to determine the least-squares estimates of the ARMA model parameters as shown below.
Use of the ARMA model in reactor noise analysis
819
Let us determine A.,p(z) and/},.,v(z) first for p = 1 in the following form: An,l(Z ) = aoAn_l,l(z ) + ffl-'An,l(Z)Z -n ~- a2an,l(Z )
*}.,1(z) = 0-oB._l,,(z) + 0-~-9.,,(z)z -~ + 0-~B.,~(z) ,
(18)
where An,1 (z) and Bn,l (z) are adjoint operator functions that satisfy the adjoint equation NAR [An'l(Z)
-
-Bn'l(Z)zn-l(I - Z
MJz-J)]C(k)
= O ,
for
k = -1,-2,...,-n
(19)
,
j=l
and 0-0, 0-1, 0"2 are constant matrices to be found. These matrices can be determined by substituting expressions (18) into Eq. (4) for k = n and into Eq. (6) for k = 1, then the following equations are obtained: 0-0Wn-l,l(n) + O'l~n,l(0)
= 0 w I 0-0w'_,,,(1) + 0 - ~ . , 1 ( 1 ) + 0-2 .,~(1) = 0 .
(20)
Here we use the following definitions: w n , v ( k ) = [an,v(z) - B , , v ( z ) M ( z ) ] C ( k ) W.,l(0) = [~.,1(z) - - 9 . , , ( z ) z " - I M ( z ) ] C ( O )
(21.a)
(21.b)
tbn,p(k) = [-4n,p(z) - 1 3 , , p ( z ) M ( z ) ] C ( k )
(21.c)
w',,v(k) = A n , v ( z ) [ M ( z ) C ( - k ) ]
(21.d)
~.,1(1)
T - B,,p(k)F
= --~.,l(n)[M(z)C(-1)]
~ - ].,~(n
ffo',v(k) = A , , p ( z ) [ M ( z ) C ( - k ) ]
T -/3,,p(k)F
- 1)r - r
.
(21.~)
(21.f)
It is seen from Eq. (20) that 0-i = aovi (i -= 1, 2) with V1 = - - W n _ l , l ( n ) w n l l ( 0 ) /]2 = --[Wtn--l,l(1) "~ /]lW-'~n,l(1)]Wtn.ll(I) •
(22)
Substituting these ui into the condition I = ao - 0 - , A n , l ( n ) + 0-2 = ao + a,B,,,a(1) + a2 ,
(23)
and taking into account that - A . , a ( n ) = Bn,l(1) we will get 0-0 as
oo = (I + . , ~ . , , ( 1 )
+ ~2) -~ .
(24)
Thus -4n,l(z) a n d / ~ n , l ( z ) in Eq. (18) can be computed with the aid of a0, ul and u2: .Z[n,l(Z ) = 0-o[dn_l,l(Z) + vl-An,l(Z)Z -n + v2An,l(Z)] Bn,l(Z) = 0-o[Bn-l,l( z)'~ b'lBn,l(Z) Z-1 + tl2Bn,l(Z)] •
(25)
820
TRAN DINH TRI
Let us go on to determine An,p(z) and/~n,p(z) for p = 2, 3,... in the following form: an,p(z ) = otoAn,l(Z ) + o Q A n , I ( Z ) + ... + otp-lan,p-l(Z) + apAn,p(z) JBn,p(Z) = otoBn,l(Z ) + o t l n n , l ( Z ) + ... + olp_lnn,p_l(Z ) + OlpBn,p(Z ) ,
(26)
here o~0, 31,..., a v are constant matrices that can be determined from the condition (27)
Or0 "31-O~1 "~ ... "~- O/p = I ,
and from the following equations received by putting Eq. (26) into Eq. (6) for k = 1, 2, ..., p: = -avwln,v(1) = -apw',,(2)
{ a0w~,l(1 )
a0w',x(2) + a l e " 1(2) a0w'.,l(3) +
+
=
Ot0Wtn,l (P) + OtlWtn,l(p) + Ot2Wln,2(p) + ... + Otp--l~3tn,p--l(p)
(28)
= --OtpWtn,p(p)
,
where w;,v(k ) mad Vo~,v(k ) are given in Eqs. (21.d) mad (21.f), respectively. To determine the matrices at, let at = avTt with l = 0, 1, 2,... , p - 1 then by substituting them into Eqs. (28) we obtain the following Gaussima linear algebraic equations: { 70w~,1(1 ) 7OWln,l(2) + 71Win,l(2) '70Win,l(3) + "71Wtn,l(3) + 72Wtn,2(3)
= --w~n,v(1) = --Wtn,p(2) = --w:,p(3)
70Wn,I(P) + "71Wn,l(P) + "72 n 2(P) + "'" + 7p-lWn,p-l(P)
(29)
= --Wn,p(P)
Obviously Eqs. (29) now could be solved directly from the first equation for 70, then to the second one for 71 with the given 70, etc. The final form of the solution of Eqs. (29) can be expressed as t 1-1 "70 = -w,,v(1)wn,1 (1) "71 = -- ['70Wtn,l( 2 ) + Wn,p(2)]Wn, t ~t--11 (2)
t -t-1 (3) 5"2 = --["7ow'0(3) + "71W',1(3) + w,,,p(3)]w.,2 I
~1
l
~1-1
"Tp--1 = --["70Wn,I(P) "{- 71Wln,I(P) + .-. "~- "Tp--2Wn,p--2(P ) "~ Wn,p(p)]Wn,p_l(p)
,
(30)
and the condition (27) gives a p = [70 31- 71 ~- 72 3t- "." "~- "7p-1 Ji- I ] - 1
.
(31)
Finaiy, the operator functions .4,,p(z) and Bn,p(z) in (26) will be computed by the expressions
A,,,p(z) = ap["7oA,,l(z) + 71/i,,a(z) + ... + 7,-lfi--,,,p-x(z) + Amp(z)] JBn,p(z) = Olp[7oSn,l(Z) + "71Bn,1(z) "gt--.. 3I- 7 p - l n n , p - l ( Z ) + Bn,p(Z)] .
(32)
Use of the ARMA model in reactor noise analysis 4. S O M E RESULTS O F T H E A P P L I C A T I O N W i t h the aim of illustration in this section we show some results of the evaluation of the noise measurements carried out at various NPPs. The numerical results were obtained by using the program package named A R M A A N A , which was built up by the author for evaluating the P W R and B W R reactor noise diagnostic measurement d a t a based on the recursive algorithm (Tran Dinh Tri, 1992a), the method developed in the previous section and the generalized S T P method (Tran Dinh Tri, 1992b). The measurement data mentioned above were processed into two d a t a files named " T I M . . . . . l A W ' and " T I M ...... 1BR". Some features of the two d a t a files are given in the following table: name of data file " T I M ...... 1AU . . . . T I M ...... 1BR number of signals 64 16 sampling interval A t , sec 0.01 0.01 fma~ = 2-~ ' H z 50 50 considered detector combinations: first signal pressure Q3 inlet pressure PB second signal vibration detector V3 neutron detector A7 third signal neutron detector C3 ionisation chamber I2 fourth signal ionisation chamber I1 thermocouple T3
A m o n g basic input parameters there are some important ones: - number of points per block
: NPOINT
- number of blocks number of frequency points
: N B L O C K = 10 : N F = 160
-
= 1024
To select the optimal order of the A R model, a certain information criterion, for example the Akaike's Information Criterion (AIC), can be applied (Fukunishi, 1977). For the given detector combinations of the data sets shown in the above table, the A I C ( n ) function for the A R model was computed for model order up to n -- 20. However, any m i n i m u m of the A I C ( n ) function up to n = 20 is not found in the case of the first combination of the " T I M . . . . . 1AU" data set, and in the case of the second combination of the " T I M ...... 1BR" data set, a minimum of the A I C ( n ) function is found at n = 19. For the both case of the given detector combinations the A I C ( n , p) funtion for the A R M A model however still were not considered. An examination of the optimal order of the A R M A for paractical noise data sets should be done in a separate paper. Therefore, in order to present some calculatinal results of the A R M A model abitraz7 indices n = 14, p = 14 and N A n = 20 are taken without any consideration of the optimal problem. So these indices are not any optimal model order. The ARMA(14,14) model parameters fi,,,p(i) and/~n,p(i) that are calculated by the m e t h o d given in section 3, are used to compute the spectra and coherences. T h e calculation results of autospectra, ordinary and partial coherences for the first d a t a set " T I M ...... •AU" are given in figs. 2, 3, 5 and 6, while the corresponding results for the second d a t a set
821
822
TRAN DINHTm " T I M . . . . . 1BR" are shown in figs. 8, 9, 11 and 12. The F F T autospectra and coherences of both of the given data sets also are presented in figs. 1, 4, 7 and 10 for comparision. Although the used indices (n = 14, p = 14 and WAR 20) are not any optimal model order, based on the received results of numerical calculations we can nevertheless conclude that in the case of the ARMA model the shape of the autospectral and coherence functions resembles more to the F F T spectrum and coherence than in the AR case, and this conclusion is valid even in the case when the order of AR part of the ARMA model is lower than that of the AR model. For example, the autospectral (APSD) and coherence (COH) functions in the ARMA(14,14) model of both data sets are better than that of the AR(20) model. For the partial quantities, similar conclusions can be drawn. =
ACKNOWLEDGMENTS The author would like to thank Dr I. Lux for his help. This work was partly supported by Hungarian Government Funds OTKA 1845 and OMFB 91-97-42-0344.
REFERENCES Fhkunishi K. (1977) Nucl. Sei. and Eng., 62, pp. 215-225. Kishida, K. (1982) Physical Langevin model and the time-series model in systems far from equilibrium. Physical Review A, Vol. 25, No. 1, pp. 496-507. Kishida, K. (1984) Equivalent random force and time-series model in system far from equilibrium. J. Math. Phys., Vol. 25, No. 5, pp. 1308-1313. Tran Dinh Tri (1990) Determination of coefficient matrices for ARMA model. KFKI Report, 1 9 9 0 - 4 5 / G , Budapest, Hungary. Tran Dinh Tri (1992a) A new algorithm for recursive estimation of ARMA parameters in reactor noise analysis. Annals of Nuclear Energy, Vol. 19, pp. 287-301. Tran Dinh Tri (1992b) Generalization of the signal transmission path method for ARMA model in reactor noise analysis. Annal, of Nuclear Energy, Vol. 19, pp. 341-345.