CC BY
24
In the case of coastal sites and sites on tidal waters, the base assumption is a storm-tide water level with this probability value. A procedure for determining such a seldom flood level is presented in Section 6.3.
On the basis of the flood runoff or of the storm tide water level, the corresponding site specific water level in the vicinity of the plant components to be protected and the protective structures of the nuclear power plant shall be determined, e.g. by hydraulics calculations.
6.2. Determination of water runoffs for a flood with a probability value of 10-4/a for inland water sites
To determine the decisive water runoff of floods for inland water sites, a statistical extrapolation based on the convention introduced in [13] covering the simultaneous occurrence of unfavourable influences shall normally be employed. In this case the following standardized distribution function shall be employed in its expanded form:
HQ(10-i) = MHQ + Shq • k(10-4)5
where
HQ(10-4): peak-level water runoff of a flood with a probability value of 10-4, in m3/sec,
MHQ: average peak-level water runoff of a flood over an extended measurement period, in m3/sec, sHQ: standard deviation of peak-level water runoff of a flood over an extended measurement period, in m3/sec,
k -4): frequency factor for an event with the probability value 10-4/a.
In this procedure the peak-level water runoff of a flood event with a probability value of 10-4/a is extrapolated from the peak-level water runoff of a flood event with a probability value of 10-2/a. Hereby, it is assumed that the peak-level water runoff of a flood event with a probability value of 10-2/a is determined using standard statistical procedures [5]. The extended extrapolation is then performed using the Pearson-III probability distribution. This is the basis on which the necessary frequency factors are determined. The convention introduced in [13] calls for a maximization of the skewness coefficient, c, to the value of c = 4.
The statistical parameters MHQ and sHQ and the actual skewness coefficient, c, shall be calculated from the observed data of a representative flood level.
The frequency factor, k _4), shall be calculated as the product of the frequency factor, k 2), and a
quotient, f, as follows:
k (10-4) = k (10-2) f.
The frequency factor, k 2), for a flood with the probability value of 10-2/a shall be interpolated from
Table 2 based on the actual skewness coefficient, c, of the observed data. The frequency factor may, alternatively, be calculated with sufficient accuracy from the following equation
k(10-2) = 2.3183 + 0.7725 x c - 0.0650 x c2.
The quotient, f, shall be calculated for a maximized skewness coefficient, c = 4, from the frequency factor, k (10-4)max, and from the frequency factor, k (10-2)max, as follows
f = k (10-4) / k (10-2) = 12.36/4.37 = 2.8.
(10 ) max (10 ) max
Both frequency factors are independent of site-specific data.
6.3. Derivation of water levels for a storm tide with a probability value of 10-4/a for coastal sites and sites on tidal waters
The storm tide water levels for nuclear power plants on coastal sites and sites on tidal waters shall normally be derived employing the following statistical extrapolation procedure. The water level for a storm tide with a probability factor of 10-4/a, SFWH(10-4), shall be determined as the sum of a base value,
BHWH(10-2), and an extrapolation difference, ED, as follows:
SFWH(10-4} = BHWH(10-2} + ED ,
where
BHWH(10-2): base value of the water level for a storm tide with a probability value of 10-2/a at the site,
ED: extrapolation difference representing the water level difference between the water level of a
storm tide with a probability value of 10-4/a and the base value.
The base value, BHWH^), shall be determined on the basis of a quantitative statistical extreme-value
analysis (in accordance, e.g., with [9] and [10]) taking relevant parameters [5] into consideration. The quality of the data shall also be taken into consideration.
The base value can be determined employing suitable statistical procedures, because
- the spread of the base values, bhwh (10-2), is relatively small due to the usually extensive and high
quality water-level time series available for coasts and tidal waters,
- the bhwh 2) water level as a function of the observation duration of the individual time series still
is partly in the interpolation region or in the near extrapolation region,
- the bhwh 2)water level is assured by extensive investigations and is verifiable by physical as
well as numerical models. The water-level data shall be homogenized considering that the storm-tide water levels are dependent on the development of the water level at the coast - especially the secular rise of the sea level - as well as on the anthropogenic changes to the tidal waters.
The extrapolation difference for coasts or for the mouths of tidal rivers shall be determined, e.g., in accordance with [9] and [10].
The local tide-related excessive wave amplitude is not included in the extrapolation difference.
Table 2. Frequency factors, k, for an event with a probability factor of 10-2/a and the actual skewness coefficient, c, of the observed data
c 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
k 2.326 2.399 2.472 2.544 2.615 2.685 2.755 2.823 2.891 2.957
c 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9
k 3.022 3.086 3.149 3.211 3.271 3.330 3.388 3.444 3.499 3.552
c 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9
k 3.605 3.656 3.705 3.753 3.800 3.845 3.889 3.931 3.973 4.012
c 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0
k 4.051 4.088 4.124 4.159 4.192 4.224 4.255 4.285 4.314 4.341 4.367
7. Results of a sensitivity study for a flood event with extreme waves in the German north sea
PSA regulations consider extreme events of recurrence intervals of 10000 years. Beside the frequently occurring extreme storm flood events, it has been investigated to which extent other events have to be considered. One example is the possible impact of an extreme wave triggered by an offshore landslide. Geotechnical records give evidence for three tsunamis in the North Sea between 8000 and 1500 years ago [3]. One well-explored source region is the Storegga slide, which was released approximately 8100 cal years bp [14].
In the framework of a dedicated study on behalf of BfS, a numerical model was applied by the Centre of Marine Environmental Sciences (MARUM) of University of Bremen to simulate the propagation and development of extreme waves in the North Sea towards the German Bight.
Based on an implicit finite differences modelling system, a hydrodynamic numerical model of the European continental shelf sea has been set-up in order to provide high-resolution data on the hydrodynamics of the North Sea. The rectilinear spherical grid covers the region between W13/N48 and E13/N62 with a resolution of 2.5nm (1/24°) in the latitudinal and 3.75nm (1/16°) in the longitudinal direction (Figure. 1). The model bathymetry was interpolated from sea floor topography derived by satellite altimetry and digitised sea-charts [15]. For this study the propagation of an extreme wave event (tsunami) initiated by a hypothetical slide at the continental margin off the Norwegian continental margin has been simulated. Soliton waves were prescribed as water level boundary conditions at the northern open sea boundary of the model.
Figure 1. Domain of interest: Model grid nodes are indicated as blue dots. The red line denotes the position of open model boundaries. Tidal gauge stations are chosen for further analysis
As the real height of a possible wave cannot be defined, a range of different wave heights was tested. Simulations show the propagation of the wave across the model domain, considering uniform mean sea level as initial surface elevation condition: After entering the North Sea through the northern boundary, the wave is partly deflected towards the West, because of Coriolis force effects, and partly moves in southern direction through the Norwegian deep. The deflected wave then approaches the British East coast and partly reflects back into the North Sea. Here the primary wave and the reflected wave super-impose into complex patterns. It takes about 8.5 hours for the first wave to reach the German Bight.
Figure 2. Extreme waves as calculated at the coastal stations featuring the first direct wave and reflected wave
The heights and characteristics of the waves at the three coastal stations are similar, all featuring the first direct wave, and about four hours later the reflected wave, which then reaches higher maximum water levels (Figure 2). Generally a significant reduction in wave height from the boundary to the German Bight due to bottom friction can be observed.
The characteristics of the wave triggered by the ancient Storegga event were simulated in [12]. Considering their calculated wave height of 3 meters at the Northern boundary of the model, results in maximum deviations of about 0.5 to 0.7m at the tidal gauges in the German Bight.
In contrast to the simulations described above, the natural hydrodynamics of the North Sea are driven by tidal and meteorological forcing. Thus the super-position of the extreme wave with the astronomical tidal conditions of the North Sea has been simulated (Figure 3). Although non-linear effects are present, generally a linear superposition of tidal elevation and extreme wave dimensions based on uniform mean sea level seem to be possible. It is noted that in the German Bight the transformed extreme wave is of much smaller height than the astronomical tidal signal: The effect of an extreme wave at the gauges Helgoland and Cuxhaven results in less than 10% of the tidal range and only one fifth of the expected surface elevation of a light storm flood, as defined by German hydrographic agencies. Similarly at gauge "Alte Weser", the extreme wave is damped to 0.55m, which is about 17 percent of the tidal range and less than one third of a light storm flood.
Considering the natural hydrodynamic conditions as tides and storm surges of the German Bight, the modelled impact of an extreme event that could be triggered by mass slide events at the northern continental margin, seems negligible.
8. Concluding remarks
The approach for a probabilistic assessment of external hazards to be applied within comprehensive safety reviews of NPP in Germany starts with a screening process, which should not be too conservative so that the number of scenarios and buildings remains manageable for the detailed quantitative analysis. However, it has to be ensured that all relevant areas are investigated within the quantitative analysis. These screening procedures are specific according to the different types of hazards.
However, for those areas which have not been screened out or where a coarse meshed analysis is not sufficient it is compulsory to perform a quantitative analysis as a second step. Finally, the frequency of initiating events induced by the respective hazard, the main contributors and the calculated core damage frequency are determined.
On international level, as already mentioned earlier, there exist some standards and guidelines [1], [8], but they are on a very general level and do not allow to perform a PSA of external flooding in a comparable manner for different plants. Moreover a full scope PSA for external flooding of a nuclear power plant is not available to date.
In Germany the graded process defines only one NPP for which no analysis will be necessary because of its high-grounded level compared with the surroundings. For the other plants probabilistic considerations will be necessary with a different extent of detail.
Compared with other external events (e.g. unintended airplane crash and external pressure wave), which can have frequencies as low as 10-7/a the occurrence frequency of external flooding can be expected substantially higher.
In the case of tidal-river NPPs the value will be higher than the risk of a seismic event due to the seismic situation (Intensity < 6) of the respective sites. The results of a simulation study have shown that an extreme wave in the North Sea towards to the German Bight triggered by an offshore landslide did not indicate significant impacts on the flooding risk of coastal sites. It is not expected that these conditions will be different compared to tidal-river NPP sites, this has, however, to be answered by flood hazard analyses for these sites.
For NPPs in the Southwest of Germany, the contribution of the seismic hazard to the total core damage frequency is expected to be higher compared with external flooding, but the overall core damage frequency remains dominated by internal events and internal hazards.
It should be underlined that the probabilistic assessment of external hazards, although an important part of PSA, has not yet achieved the same level of methodological maturity as being typical for other
disciplines of PSA. Therefore, it is intended to conduct a kind of pilot study to get feedback from these analyses for an improvement of the German guidance documents.
However independently from NPPs and other industrial facilities floods from rivers, estuaries and the sea threaten many millions of people in Europe. Flooding is the most widely distributed of all natural hazards across Europe, causing distress and damage wherever it happens.
Previous research has improved understanding of individual factors but many complex interactions need to be addressed for flood mitigation in practice. Thus the first round of the Sixth Framework Programme of the European Commission (2002-2006) included an "Integrated Project" on flood risk management, called FLOODsite.
To achieve the goal of integrated flood risk management, the FLOODsite project has brought together managers, researchers and practitioners from a range of governmental, commercial and research organisations, all devoted to various, but complementary, aspects of flood risk management.
The FLOODsite project covers the physical, environmental, ecological and socio-economic aspects of floods from rivers, estuaries and the sea. The project is arranged into seven themes covering:
- Risk analysis - hazard sources, pathways and vulnerability of receptors.
- Risk management - pre-flood measures and flood emergency management.
- Technological integration - decision support and uncertainty.
- Pilot applications - for river, estuary and coastal sites.
- Training and knowledge uptake - guidance for professionals, public information and educational material.
- Networking, review and assessment.
- Co-ordination and management.
Within these themes there are over 30 project tasks including the pilot applications in Belgium, the Czech Republic, France, Germany (in particular for flood event measures and pilot application sites), Hungary, Italy, the Netherlands, Spain, and the UK. Published results are expected in 2007.
References
[1] ANS, American Nuclear Society (2003). Draft of the External Events PRA Methodology Standard. BSR/ANS 58.21-200X.
[2] Berg, H.-P. & Görtz, R. (2006). Probabilistic Safety Assessment of External Flooding of Nuclear Power Plants. Proc. of the European Safety and Reliability Conference ESREL'06, Safety and Reliability for Managing Risk, Estoril, Portugal, Vol. 2, Taylor & Francis, London, 1341 - 1346.
[3] Bondevik, S., Lovholt, F., Harbitz, C.B., Mangerud, J., Dawson, A.G. & Svendsen, J.I., (2005). The Storegga Slide Tsunami - Comparing Field Observations with Numerical Simulations. Marine and Petroleum Geology 22, 195-208.
[4] Bundesministerium für Umwelt, Naturschutz und Reaktorsicherheit (BMU) (2005). Safety Review for NPP According to § 19a of the Atomic Energy Act - Probabilistic Safety Assessment Guide (Sicherheitsüberprüfung für Kernkraftwerke gemäß §19a des Atomgesetzes - Leitfaden Probabilistische Sicherheitsanalyse, 31. Januar 2005, Bekanntmachung vom 30. August 2005, Bundesanzeiger Nr. 207a vom 03.
[5] Deutscher Verband für Wasserwirtschaft und Kulturbau (1999). Statistical Analysis of Peak Level Water Runoffs (Statistische Analyse von Hochwasserabflüssen), DVWK-Merkblatt 251.
[6] Facharbeitskreis Probabilistische Sicherheits-analyse für Kernkraftwerke (2005). Methods for PSA for NPPs (Methoden zur probabilistischen Sicherheitsanalyse für Kernkraftwerke, Stand: August 2005), BfS-SCHR - 37/05, Salzgitter.
[7] Facharbeitskreis Probabilistische Sicherheits-analyse für Kernkraftwerke (2005). Data for PSA for NPPs (Daten zur probabilistischen Sicherheits-analyse für Kernkraftwerke, Stand: August 2005), BfS-SCHR -38/05, Salzgitter.
[8] International Atomic Energy Agency (2003). Flood Hazard for Nuclear Power Plants on Coastal and River Sites. Safety Guide, No. NS-G-3.5, Vienna.
[9] Jensen, J. (2000). Probability of Occurrence of Storm Floods - Statistical View, HANSA, Vol 137, Nr. 12, 60 -66.
[10] Jensen, J & Frank, T. (2003). On the Determination of Water Levels from Storm-Floods with a very Small Probability Value, Die Küste, Spezialausgabe, Nr. 67.
[11] Jensen, J. et. al. (2003). New Procedures for the Assessment of Rare Water Levels from Storm Floods (Neue Verfahren zur Abschätzung von seltenen Sturmflutwasserständen), HANSA, Vol. 140, Nr. 11,6879.
[12] Kerntechnischer Ausschuss (2004). Flood Protection of NPP (Schutz von Kernkraftwerken gegen Hochwasser), KTA 2207.
[13] Kleeberg, H.-B. & Schumann, A. H. (2001). Derivation of Water Runoffs with Small Exceedance Frequencies (Ableitung von Bemessungsabflüssen kleiner Überschreitungswahrscheinlichkeiten), Wasserwirtschaft, Vol 21, Nr. 2, 2001, 90 - 95.
[14] Lovholt, F., Harbitz, C.B. & Haugen, K.B. (2005). A Parametric Study of Tsunamis Generated by Submarine Slides in the Ormen Lange/Storegga Area off Western Norway. Marine and Petroleum Geology 22, 219-231.
[15] Smith, W. & Sandwell, D. (1997). Global Sea Floor Topography from Satellite Altimetry and Ship Depth Soundings. Science 277, 1956-1962
RECURRENT NEURAL NETWORKS FOR DYNAMIC RELIABILITY ANALYSIS
Cadini Francesco, Zio Enrico, Pedroni Nicola
Department of Nuclear Engineering, Polytechnic of Milan, Milan, Italy
Keywords
dynamic reliability analysis, infinite impulse response-locally recurrent neural network, long-term non-linear dynamics, system state memory, simplified nuclear reactor
Abstract
A dynamic approach to the reliability analysis of realistic systems is likely to increase the computational burden, due to the need of integrating the dynamics with the system stochastic evolution. Hence, fast-running models of process evolution are sought. In this respect, empirical modelling is becoming a popular approach to system dynamics simulation since it allows identifying the underlying dynamic model by fitting system operational data through a procedure often referred to as 'learning'. In this paper, a Locally Recurrent Neural Network (LRNN) trained according to a Recursive Back-Propagation (RBP) algorithm is investigated as an efficient tool for fast dynamic simulation. An application is performed with respect to the simulation of the non-linear dynamics of a nuclear reactor, as described by a simplified model of literature.
1. Introduction
Dynamic reliability aims at broadening the classical event tree/ fault tree methodology so as to account for the mutual interactions between the hardware components of a plant and the physical evolution of its process variables. The dynamical aspects concern the ordering and timing of events in the accident propagation, the dependence of transition rates and failure criteria on the process variable values, the human operator and control actions. Obviously, a dynamic approach to reliability analysis would not bear any significant added value to the analysis of systems undergoing slow accidental transients for which the control variables do not vary in such a way to affect the component transition rates and/or to demand the intervention of the control.
Dynamic reliability methods are based on a powerful mathematical framework capable of integrating the interactions between the components and the environment in which they function. These methods perform a more realistic modelling of the system and hence improve the quality and accuracy of risk assessment studies. A formal approach to incorporating the dynamic behaviour of systems in risk analysis was formulated under the name Probabilistic Dynamics [10]. Several methods for tackling the solution to the dynamic reliability problem have been formulated over the past ten years [1], [9], [13], [15], [16], [20]. Among these, Monte Carlo methods have demonstrated to be particularly efficient in taking up the numerical burden of such analysis, while allowing for flexibility in the assumptions and for a thorough uncertainty and sensitivity analysis [14], [16].
For realistic systems, a dynamic approach to reliability analysis is likely to require a significant increase in the computational efforts, due to the need of integrating the dynamic evolution, with its characteristic times, with the system stochastic evolution characterized by very different time constants. The fast increase in computing power has rendered, and will continue to render, more and more feasible the incorporation of dynamics in the safety and reliability models of complex engineering systems. In particular, as mentioned above, the Monte Carlo simulation framework offers a natural environment for estimating the reliability of systems with dynamic features. However, the high reliability of systems and components favours the adoption of forced transition schemes and leads, correspondingly, to an increment of the integration of physical models in each trial. Thus, the time-description of the dynamic processes may render the Monte Carlo simulation quite burdensome and it becomes mandatory to resort to fast-running models of
process evolution. In these cases, one may resort to either simplified, reduced analytical models, such as those based on lumped effective parameters [2], [7], [8], or empirical models. In both cases, the model parameters have to be estimated so as to best fit to the available plant data.
In the field of empirical modelling, considerable interest is devoted to Artificial Neural Networks (ANNs) because of their capability of modelling non-linear dynamics and of automatically calibrating their parameters from representative input/output data [16]. Whereas feedforward neural networks can model static input/output mappings but do not have the capability of reproducing the behaviour of dynamic systems, dynamic Recurrent Neural Networks (RNNs) are recently attracting significant attention, because of their potentials in temporal processing. Indeed, recurrent neural networks have been proven to constitute universal approximates of non-linear dynamic systems [19].
Two main methods exist for providing a neural network with dynamic behaviour: the insertion of a buffer somewhere in the network to provide an explicit memory of the past inputs, or the implementation of feedbacks.
As for the first method, it builds on the structure of feedforward networks where all input signals flow in one direction, from input to output. Then, because a feedforward network does not have a dynamic memory, tapped-delay-lines (temporal buffers) of the inputs are used. The buffer can be applied at the network inputs only, keeping the network internally static as in the buffered multilayer perceptron (MLP) [11], or at the input of each neuron as in the MLP with Finite Impulse Response (FIR) filter synapses (FIR-MLP) [4]. The main disadvantage of the buffer approach is the limited past-history horizon, which needs to be used in order to keep the size of the network computationally manageable, thereby preventing modelling of arbitrary long time dependencies between inputs and outputs [12]. It is also difficult to set the length of the buffer given a certain application.
Regarding the second method, the most general example of implementation of feedbacks in a neural network is the fully recurrent neural network constituted by a single layer of neurons fully interconnected with each other or by several such layers [18]. Because of the required large structural complexity of this network, in recent years growing efforts have been propounded in developing methods for implementing temporal dynamic feedback connections into the widely used multi-layered feedforward neural networks. Recurrent connections can be added by using two main types of recurrence or feedback: external or internal. External recurrence is obtained for example by feeding back the outputs to the input of the network as in NARX networks [5], [17]; internal recurrence is obtained by feeding back the outputs of neurons of a given layer in inputs to neurons of the same layer, giving rise to the so called Locally Recurrent Neural Networks (LRNNs) [6].
The major advantages of LRNNs with respect to the buffered, tapped-delayed feedforward networks and to the fully recurrent networks are [6]: 1) the hierarchic multilayer topology which they are based on is well known and efficient; 2) the use of dynamic neurons allows to limit the number of neurons required for modelling a given dynamic system, contrary to the tapped-delayed networks; 3) the training procedures for properly adjusting the network weights are significantly simpler and faster than those for the fully recurrent networks.
In this paper, an Infinite Impulse Response-Locally Recurrent Neural Network (IIR-LRNN) is adopted together with the Recursive Back-Propagation (RBP) algorithm for its batch training [6]. In the IIR-LRNN the synapses are implemented as Infinite Impulse Response digital filters, which provide the network with system state memory.
The proposed neural approach is applied to a highly non-linear dynamic system of literature, the continuous time Chernick model of a simplified nuclear reactor [8]: the IIR-LRNN is devised to estimate the neutron flux temporal evolution only knowing the reactivity forcing function. The IIR-LRNN ability of dealing with both the short-term dynamics governed by the instantaneous variations of the reactivity and the long-term dynamics governed by Xe oscillations is verified by extensive simulations on training, validation and test transients.
The paper is organized as follows: in Section 2, the IIR-LRNN architecture is presented in detail together with the RBP training algorithm; in Section 3, the adopted neural approach is applied to simulate the reactor neutron flux dynamics. Finally, some conclusions are proposed in the last Section.
2. Locally Recurrent Neural Networks
2.1. The IIR-LRNN architecture and forward calculation
A LRNN is a time-discrete network consisting of a global feed-forward structure of nodes interconnected by synapses which link the nodes of the k-th layer to those of the successive (k + 1)-th layer, k = 0, 1, ..., M, layer 0 being the input and M the output. Differently from the classical static feed-forward networks, in an LRNN each synapse carries taps and feedback connections. In particular, each synapse of an IIR-LRNN contains an IIR linear filter whose characteristic transfer function can be expressed as ratio of two polynomials with poles and zeros representing the AR and MA part of the model, respectively.
For simplicity of illustration, and with no loss of generality, we start by considering a network constituted by only one hidden layer, i.e. M = 2, like the one in Figure 1. At the generic time t, the input to
N 0
the LRNN consists of a pattern x(t) e ^ , whose components feed the nodes of the input layer 0 which simply transmit in output the input received, i.e. x0m(t) = xm(t), m = 1, 2, ..., N0. A bias node is also typically inserted, with the index m = 0, such that x00(t) = 1 for all values of t. The output variable of the m-th input node at time t is tapped a number of delays L1nm - 1 (except for the bias node output which is not tapped, i.e. L1n0 - 1 = 0) so that from each input node m ^ 0 actually L1nm values, x0m(t), x0m(t - 1), x0m(t - 2), ..., x0m(t -L1nm + 1) are processed forward through the synapses connecting input node m to the generic hidden node n = 1, 2, ... N1. The L1nm values sent from the input node m to the hidden node n are first multiplied by the respective synaptic weights w1nm(p), p = 0, 1, ., L1nm - 1 being the index of the tap delay (the synaptic weight w1n0(p) connecting the bias input node m = 0 is the bias value itself) and then processed by a summation operator to give the MA part of the model with transfer function
wL(0) + wU)B + wL(2) b 2 +... + w'mLnm BLnm (1)
B being the usual delay operator of unitary step. The finite set of weights w1nm(p) which appear in the MA model form the so called impulse response function and represent the components of the MA part of the synaptic filter connecting input node m to hidden node n. The weighed sum thereby obtained, y1nm, is fed back, for a given number of delays I1nm (I1n0 = 0 for the bias node) and weighed by the coefficient v1nm(p) (the AR part of the synaptic filter connecting input node m to hidden node n, with the set of weights v1nm(p) being the so-called AR filter's impulse response function), to the summation operator itself to give the output quantity of the synapse ARMA model:
LL-1
L (t) = Z WUp) X0(t - P) + Z VUp) yL(t - P)
(2)
p=0
P=1
This value represents the output at time t of the IIR-filter relative to the nm-synapse, which connects the m-th input neuron to the n-th hidden neuron. The first sum in (2) is the MA part of the synaptic filter and the second is the AR part. As mentioned above, the index m = 0 usually represents the bias input node, such that x00(t) is equal to one for all values of t, L1n0 - 1 = I'n0 = 0 and thus, y'nC(t) = w1n0(0). The quantities y1nm(t), m = 0, 1, ..., N0, are summed to obtain the net input s1n(t) to the non-linear activation function f (), typically a sigmoid, Fermi function, of the n-th hidden node, n = 1, 2, ...N1:
N0
s1(t) = X ylm (t). (3)
m=0
The output of the activation function gives the state of the n-th hidden neuron, x1n(t):
x1(t) = f1 k(t)]. (4)
The output values of the nodes of the hidden layer 1, x1n(t), n = 1, 2, ..., N1, are then processed forward along the AR and MA synaptic connections linking the hidden and output nodes, in a manner which is
absolutely analogous to the processing between the input and hidden layers. A bias node with index n = 0 is also typically inserted in the hidden layer, such that x10(t) = 1 for all values of t.
The output variable of the n-th hidden node at time t is tapped a number of delays LMrn - 1 ( = 0 for the bias node n = 0) so that from each hidden node n actually LMrn values, x1n(t), x1n(t - 1), x1n(t - 2), ..., x1n(t -LMrn + 1), are processed forward through the MA-synapses connecting the hidden node n to the output node r = 1, 2, ..., NM. The LMrn values sent from the hidden node n to the output node r are first multiplied by the respective synaptic weights w1Mrn(P), p = 0, 1, ..., LMrn - 1 being the index of the tap delay (the synaptic weight
w r0 connecting the bias hidden node n = 0 is the bias value itself) and then processed by a summation operator to give the MA part of the model with transfer function
, ,.,M D , „,M R2 M T>Ln-1
w ,„-, + w ,,,D + w +... + w ,rM , B . (5)
rn (0) rn(1) rn (2) rn(IM-1) v '
The sum of these values, yMrn, is fed back, for a given number of delays t^m (]Mr0 = 0 for the bias node) and weighed by the coefficient VMrn(p) (the AR part of the synaptic filter connecting hidden node n to output node r, with the set of weights VMrn(p) being the corresponding impulse response function), to the summation operator itself to give the output quantity of the synapse ARMA model:
(t) = y w rn( p ) n (t - p) + y vrn(p)ym (t - p).
(6)
p=0
p=1
As mentioned before, the index n = 0 represents the bias hidden node, such that x10(t) is equal to one for all values of t, LMr0 - 1 = tIr0 = 0 and thus, yir0(t) = wMr0(0).
The quantities yMrn(t), n = 0, 1, ..., N1, are summed to obtain the net input SMr(t) to the non-linear activation function fM(-), also typically a sigmoid, Fermi function, of the r-th output node r = 1, 2, ..., NM:
N
sM (t) = y yM (t). (7)
n=0
The output of the activation function gives the state of the r-th output neuron, xMr(t):
xM (t) = fM [sM (t)]. (8)
The extension of the above calculations to the case of multiple hidden layers (M > 2) is straightforward. The time evolution of the generic neuron j belonging to the generic layer k = 1, 2, ..., M is described by the following equations:
xk (t) = fk [sk (t)] (= 1 for the bias node, j = 0), (9)
Nk-1
s) (t) = Z yj (t), (10)
I=0
L"t-1 'k yj(t) = Iwki(p)xk-1(t-p) + Xvki(p)yj(t-P) . (U)
p=0 p=1
Note that if all the synapses contain only the MA part (i.e., f^ = 0 for all j, k, l), the architecture reduces to a FIR-MLP and if all the synaptic filters contain no memory (i.e., Lj - 1 = 0 and fji = 0 for all j, k, I), the classical multilayered feed-forward static neural network is obtained.
2.2. The Recursive Back-Propagation (RBP) algorithm for batch training
The Recursive Back-Propagation (RBP) training algorithm [6] is a gradient - based minimization algorithm which makes use of a particular chain rule expansion rule expansion for the computation of the necessary derivatives. A thorough description of the RBP training algorithm is given in the Appendix at the end of the paper.
INPUT (k = 0) HIDDEN (k = 1) OUTPUT (k = 2 = M)
N0 = 1 N1 = 2 NM = 1
L1ii - 1 = 1 L12i - 1 = 1 LMii - 1 = 2 LMi2 - 1 = 1
1 1 i i I I I I i I M = 1 IMi2 = 0
Figure 1. Scheme of an IIR-LRNN with one hidden layer
3. Simulating reactor neutron flux dynamics by LRNN
In general, the training of an ANN to simulate the behaviour of a dynamic system can be quite a difficult task, mainly due to the fact that the values of the system output vector y(t) at time t depend on both the forcing functions vector x(-) and the output y() itself, at previous steps:
y(t) = F(x(t), x(t - 1),..., y(t - 1),..., 0), (12)
where 0 is a set of adjustable parameters and F( ) the non-linear mapping function describing the system dynamics.
In this Section, a locally recurrent neural network is trained to simulate the dynamic evolution of the neutron flux in a nuclear reactor.
3.1. Problem formulation
The reference dynamics is described by a simple model based on a one group, point kinetics equation with non-linear power reactivity feedback, combined with Xenon and Iodine balance equations [8]:
a d®
A-
dt
(p -YO
cL f
O
dXe = YxeLf O + A'' - AXeXe - cr^Xe® (13)
^ = 7l L f O-Ajl dt
where O, Xe and l are the values of flux, Xenon and Iodine concentrations, respectively.
The reactor evolution is assumed to start from an equilibrium state at a nominal flux level 00= 4.661012 n/cm2s. The initial reactivity needed to keep the steady state is p0 = 0.071 and the Xenon and Iodine concentrations are Xe0 = 5.73 1015 nuclei/cm3 and l0 = 5.811015 nuclei/cm3, respectively. In the following, the values of flux, Xenon and Iodine concentrations are normalized with respect to these steady state values.
The objective is to design and train a LRNN to reproduce the neutron flux dynamics described by the system of differential equations (13), i.e. to estimate the evolution of the normalized neutron flux O(t), knowing the forcing function p(t).
Notice that the estimation is based only on the current values of reactivity. These are fed in input to the locally recurrent model at each time step t: thanks to the MA and AR parts of the synaptic filters, an estimate of the neutron flux O(t) at time t is produced which recurrently accounts for past values of both the network's inputs and the estimated outputs, viz.
O(t) = F(p(t),p(t -1),..., O(t -1),..., 0) (14)
where 0 is the set of adjustable parameters of the network model, i.e. the synaptic weights.
On the contrary, the other non-measurable system state variables, Xe(t) and l(t), are not fed in input to the LRNN: the associated information remains distributed in the hidden layers and connections. This renders the LRNN modelling task quite difficult.
3.2. Design and training of the LRNN
The LRNN used in this work is characterized by three layers: the input, with two nodes (bias included); the hidden, with six nodes (bias included); the output with one node. A sigmoid activation function has been adopted for the hidden and output nodes.
The training set has been constructed with Nt = 250 transients, each one lasting T = 2000 minutes and sampled with a time step At of 40 minutes, thus generating np = 50 patterns. Notice that a temporal length of 2000 minutes allows the development of the long-term dynamics, which are affected by the long-term Xe oscillations. All data have been normalized in the range [0.2, 0.8].
Each transient has been created varying the reactivity from its steady state value according to the following step function:
P(t) = iPo t"T (15)
U+ap t > ts
where Ts is a random steady-state time interval and Ap is random reactivity variation amplitude. In order to build the 250 different transients for the training, these two parameters have been randomly chosen within the ranges [0, 2000] minutes and [-5 10-4, +5 10-4], respectively.
The training procedure has been carried out on the available data for nepoch = 200 learning epochs (iterations). During each epoch, every transient is repeatedly presented to the LRNN for nrep = 10 consecutive
