Sensitivity analysis of an Advanced Gas-cooled Reactor control rod model

A model has been made of the primary shutdown system of an Advanced Gas-cooled Reactor nuclear power station. The aim of this paper is to explore the use of sensitivity analysis techniques on this model. The two motivations for performing sensitivity analysis are to quantify how much individual uncertain parameters are responsible for the model output uncertainty, and to make predictions about what could happen if one or several parameters were to change. Global sensitivity analysis techniques were used based on Gaussian process emulation; the software package GEM-SA was used to calculate the main effects, the main effect index and the total sensitivity index for each parameter and these were compared to local sensitivity analysis results. The results suggest that the system performance is resistant to adverse changes in several parameters at once.


Introduction
The United Kingdom has seven Advanced Gas-cooled Reactor nuclear power stations (AGRs), which provide around 20% of its electricity. The AGR design was developed in the 1970s and is unique to Britain. The primary shutdown mechanism is provided by the control rods which absorb the neutrons needed to sustain a chain reaction of uranium fissions in the reactor core. The control rods are continually raised and lowered in order to maintain a critical reaction. A reactor typically has around 80 control rods, each with its own actuator. Should the reactor exceed its normal operating conditions, the control rods will be released by an electromagnetic clutch and insert into the core under gravity, shutting down the reactor. This system was designed experimentally and is regularly tested to ensure the rods will enter the core quickly enough to shutdown the reactor with a sufficient safety margin. The large amounts of collected data and modern modelling techniques, give an opportunity to understand and monitor the primary shut down system performance at a more detailed level; this is beneficial for managing the plant commercially and for giving early warning of any potential performance issues. The objective of this paper is to develop a mathematical model of the system and explore the use of probabilistic sensitivity analysis on this model. A future intended use for the model is as a condition monitoring tool, where system identification techniques are used to estimate parameter values which fit the model output to a set of insertion test data. The parameter values will then be used to attempt to assess the condition of the system components.
Sensitivity analysis is concerned with how a model's inputs affect its output. In the context of modelling control rods there are two main uses for sensitivity analysis. The first is to investigate how the uncertainty of individual model parameters is responsible for the uncertainty of the model output. This is useful when developing and refining the model, as effort can be focused on the most important parameters. The second use for sensitivity analysis is for making predictions about the effects of changing parameters on the performance of the system, i.e. what will be the effect on model output if one or several model parameters deviate from their original values?
It is relatively straightforward to assess the local sensitivity of a model to its parameters, by partially differentiating the model output with respect to different parameters. However this method is not particularly informative because it fails to take into account nonlinear responses. A slightly more informative technique involves running the model for a range of parameter values, keeping the others constant. Both of these methods fail to take into account the fact that sensitivity to a parameter could vary as other parameters change. For a model with many parameters assessing the effects of all possible parameter combinations is challenging (Saltelli et al., 2000). Global sensitivity analysis techniques investigate the entire range of the possible input space using statistical methods; these can be time consuming and still require careful interpretation.
Monte Carlo analysis can be used to sample from the probability distribution of model outputs, given a set of probability distributions of model inputs. These output distributions can be used to infer global sensitivity qualities (Doubilet et al., 1985;Helton, 1993). Whilst this is effective, it can be extremely computationally expensive, especially if there are many parameters of interest. A way of reducing the expense of global sensitivity analysis given a computationally expensive model is to use a surrogate model -a model of a model. This still requires many model runs but far fewer than Monte Carlo analysis. A technique based on Fourier amplitude sensitivity testing (FAST) provides an elegant way of estimating the contribution of input uncertainty to output uncertainty, however this method is limited to investigating the main effects of parameters and does not give information regarding interactions (McRae et al., 1982;Saltelli et al., 1999).
Choosing sensitivity analysis techniques requires a compromise between robustness, computational cost, ease of implementation and conceptual simplicity. Within the nuclear industry, robustness is generally preferred at the expense of computational cheapness (Cornell, 1986). In the current case, the purpose of the model and sensitivity measures is to assist in decision making and increase understanding of the system. It is desirable that the meaning of the measures used and the concepts behind them are sufficiently intuitive that someone with little knowledge of sensitivity analysis is able to confidently use them.
The technique chosen in this investigation is a Bayesian approach to surrogate modelling developed in Oakley and O'Hagan (2002), which will be introduced in more detail below. It was chosen as it is both computationally efficient and robust, and although the mathematics behind it is relatively demanding, the basic principles are not difficult to understand.
The paper is laid out as follows: Section two introduces the model of the control rod system. Section 3 describes the system identification techniques used to estimate unknown model parameters. Section 4 gives an outline of the sensitivity analysis techniques used. The results are presented in Section 5 and the paper finishes with some conclusions in Section 6.

The model
A schematic of the system of interest is shown in Fig. 1 and a more detailed sketch of the brake system is given in Fig. 2. The governor shaft is connected to the motor by an electromagnetic clutch (not shown). If power to the clutch is lost then the governor shaft will be released and the rods will insert into the core. A two-stage braking mechanism is attached to the governor shaft. The primary brake is driven by flyweights and is dependent on the rod velocity.
The secondary brake is driven by a lead screw, which is connected to the bevel shaft by gears (not shown). Key assumptions made during the derivation of the model structure are: The effects of the chain friction, bearing friction, gear/sprocket efficiency and the friction between the side walls of the core and the rods are lumped together in a single, scaled friction parameter F f which acts as a constant force resisting the rod motion. The coefficient of friction between the governor and the brake is a constant. The drag forces from the coolant gas are directly proportional to the velocity of the rods and do not depend on the displacement. All components are assumed to be fully rigid, except the springs in the brake mechanism.
The action of the brake is sufficiently complicated that there are 9 distinct stages of rod motion. Each stage is described by 2 simultaneous differential equations, one for the motion of the rods and one for the motion of the flyweights. The points at which the model switches between stages are dictated by the position of the flyweights and the rods. The equations describing the rod's motion are shown below and a description of how they were derived is given in the appendix.
x refers to the position of the rod and h f is the flyweight angle. The rod acceleration, € x, is given by, the acceleration of the flyweights, € h f is given by, the values c and d are dependent on h f .
when h f reaches h spring the rod begins stage 2 of its motion. Stage 2.
The rods acceleration, € x, is still given by Eq. (1). During this phase the flyweights are stationary at h f ¼ h spring . when, the rod will begin stage 3 of its motion. Stage 3. The rods acceleration, € x, is still given by Eq.
(1). The flyweight acceleration is given by, when h f ¼ h max the rod will begin stage 4. Stage 4. The rod acceleration, € x, is given by, where, during this phase the flyweights are stationary at h f ¼ h fmax . When the rod reaches the final 1.5 m of its travel the rod begins stage 5 of its travel.

Stage 5.
The rod acceleration, € x, is given by, where, when Eq. (20) becomes negative the rod will begin stage 6 of its travel. Stage 6. The rod acceleration, € x, is given by, where, the acceleration of the flyweights is given by, where, when the flyweight angle reaches h spring the rod will begin stage 7. Stage 7.
The rod acceleration is given by, where, when the reaction force Frs is equal to zero the rod will begin stage 8 of its motion Stage 8. The rod acceleration is given by, where, the acceleration of the flyweights is given by, when h f ¼ h fbrake2 the rod will begin the last stage of its travel. Stage 9.
The rod acceleration is given by, where, the flyweights are stationary at h fbrake2 .
The model relies on 28 parameters which represent physical attributes of the system e.g. spring stiffnesses, masses, coefficients of friction etc. These parameters are listed in Table 1. Some of the parameters (the gear ratios and the mass of the rods) are accurately known quantities. Many of the parameters were estimated using 3D models of the system which are accurate to a tolerance of 1mm. There were some parameters which are not possible to measure directly or estimate analytically with any accuracy; these were the friction terms, spring stiffnesses, and the viscous drag coefficient.
The possible values of the unknown parameters were inferred using a Bayesian framework, which combines prior knowledge of the parameters with measured data from the system to give probability distributions of parameter values (the techniques used are described in Section 3). The values of the unknown parameters could vary while a control rod mechanism is in service, because it is possible that they are sensitive to temperature, pressure, wear or contaminants; for this reason, these parameters are estimated on an insertion-by-insertion basis.
The sensitivity analysis methods used in this investigation require the model to give a univariate output. The value chosen here is the distance the rod has inserted 4.5 s after it has been released, which was chosen as it is used as a key measure of how well the primary shutdown system is performing. The design specification is that the rods must have inserted at least 6.8 m after 4.5 s to shutdown the reactor with a sufficient safety margin. It is desirable that the rods enter the core as quickly as possible, while not traveling fast enough to cause any damage.

Bayesian system identification
Bayesian inference provides a powerful method of estimating probability distributions for parameters using prior knowledge about the system and measured data. It uses Bayes theorem: where w is a vector of model parameters and D is a set of measured data. PðwÞ is the set of prior probability distributions representing belief in the parameter values before the data has been seen.
PðDjwÞ is the likelihood: the probability that the measured response D will have been witnessed given that one believes the system is represented by the model with parameters w. The likelihood depends on the level of noise the data is corrupted with. If it is assumed that the noise is Gaussian then its variance can be included as an unknown parameter in w. PðDÞ is a normalising constant, which ensures the posterior distribution integrates to one.
PðwjDÞ is the posterior probability distribution for the parameters: the probability that the parameters take a certain value given that the system, correctly represented by the model, produced the data D. To evaluate PðwjDÞ analytically is problematic due to the difficulty of finding PðDÞ and the nonlinear nature of the model. It is possible to draw samples from the posterior allowing a picture of it to be built up, but this is often difficult due to the complex geometry of the posterior and the small size of the posterior relative to the prior. Markov Chain Monte Carlo (MCMC) methods provide a popular way of attempting to overcome these difficulties.
The Metropolis algorithm (Metropolis et al., 1953) is one of the most widely used MCMC methods, and is the technique employed in the current work. It works on the principle that it is possible to evaluate a function f ðwÞ which is proportional to the posterior, it is therefore possible to calculate the relative probability of two parameter vectors. The Markov chain in the metropolis algorithm is generated as follows: 1. An initial parameter vector, w current is generated. 2. A variation on w current is proposed: w new by taking a random sample from the proposal density, a Gaussian distribution with a mean of w current . 3. The ratio R m ¼ f ðw new Þ=f ðw current Þ is calculated, if R m > 1 then w new is accepted as the new w current . If R m < 1 then there is a probability of R m that the new vector will be accepted. 4. Repeat from step two.
It can be shown that the resulting Markov chain is ergodic and has a stationary distribution equal to PðwjDÞ. Consequently, once converged, the Markov chain can generate samples from the posterior.
MCMC techniques provide a very powerful tool; however, they have several drawbacks: Many iterations of the Markov chain are required and each iteration requires a run of the model, so MCMC methods are infeasible if the model is time consuming to run. It can take many iterations for the Markov chain to converge on its stationary distribution. It is possible for the chain to become stuck in local traps: locally optimal but globally sub-optimal areas of the parameter space. Choosing the correct proposal density width for each parameter is a crucial factor with regard to the successful application of MCMC algorithms.
In an attempt to address some of these problems a novel variant of Simulated Annealing which was proposed in Green (2015) has been employed. The algorithm utilises a ''temperature" variable b which controls the likelihood's effect on the posterior. If b ¼ 0 then the posterior is equal to the prior and when b ¼ 1, the posterior is as defined in Eq. (33). When the Markov chain is first initiated b is close to zero, after every N iterations its value is increased until it reaches 1. This gradual transition from prior to posterior increases the chances of the Markov chain escaping from local traps.
Choosing the rate at which b is varied (the annealing schedule) is important. Annealing too fast will make the chain less likely to escape local traps, but annealing too slowly will take an unnecessary amount of time. The optimum rate of annealing is not a constant, it often needs to be slower towards the start of the process and faster towards the end. It is hypothesised in Green (2015) that the optimum schedule is that where the change in information content of the target distribution is constant. This is achieved by selecting the next b value as that which increases the absolute value of the Shannon entropy by a pre-defined constant, DS, which is negative. Choosing an appropriate value for DS requires some trial and error as it will be different for different problems. In the current case, the best value of DS was found to be À0.2.
The algorithm used here also automatically tunes the width of the proposal densities. For each value of b the variance of the posterior is used to give the variance of the proposal density for the next stage. At the start of the process the proposal density is equal to the prior, so the Markov chain searches the entire possible parameter space. As b increases the posterior narrows and the Markov chain focuses on a smaller area.
Using the mean values of the Markov chain from the final stage of the annealing process as its initial values, and the tuned proposal densities; the Metropolis algorithm has a much better chance of quickly converging on a stationary distribution in the correct region of the parameter space.
Parameter distributions were estimated using data from eleven different test insertions. In order to estimate a parameter distribution, the simulated annealing algorithm was run, followed by 200,000 iterations of the Metropolis algorithm. An example histogram showing the posterior distribution of the main spring stiffness for one of the test insertions is given in Fig. 3. Using the mean values of the posterior distributions, the model was able to give a good fit to each of the sets of test data, a plot of the modelled and measured rod positions during an insertion is shown in Fig. 4. The maximum and minimum values from all of the distributions were used to inform the parameter ranges used for sensitivity analysis.

Bayesian sensitivity analysis
This section introduces the sensitivity measures which are used in this analysis and gives a brief description of the theory behind the GEM-SA package which was used to infer them (O'Hagan and Kennedy, 2004). For a more detailed description of these measures and how they are inferred, (see Oakley andO'Hagan, 2002 andWorden andBecker, 2012).

The emulator
The sensitivity analysis technique used here involves the use of an emulator -a model of the model. The model is treated as an unknown function, with the possible ranges of the input parameters specified by probability distributions. A selection of input vectors are sampled from these distributions, using a Maximin Latin Hypercube design to ensure complete, even coverage of the input space. The model is then run using these vectors as inputs to provide the training data for creating the emulator.
If it is assumed that the model is a smooth function of its inputs, then a response surface can be fitted to the training data using a least squares regression and the output can be estimated for any set of inputs. Early use of emulators in sensitivity analysis involved using the response surface to perform Monte Carlo analysis at a reduced cost (Helton, 1993). In the current case the response surface is used to provide the mean of the multivariate Gaussian probability distribution which represents the prior belief in the value of the model output. The prior distribution is then conditioned on the training data to give a posterior distribution over functions, which can be used to infer many global sensitivity values, the ones of interest will be described below. This technique is described in detail in Oakley and O'Hagan (2002).

Main effects and interactions
The model output, y can be decomposed into main effects and interactions of its input parameters, x (x denotes the vector of n input parameters fx 1 ; . . . ; x n g) where, z i;j ðx i;j Þ ¼ EðYjx i ; jÞ À z i ðx i Þ À z j ðx j Þ À EðYÞ ð 36Þ z i;j;k ðx i;j;k Þ ¼ EðYjx i ; j; kÞ À z i;j ðx i;j Þ À z i;k ðx i;k Þ À z j;k ðx j;k Þ À z i ðx i Þ À z j ðx j Þ À z k ðx k Þ ð 37Þ z i ðx i Þ is the main effect of x i ; z i;j ðx i;j Þ is the first order interaction between x i and x j ; z i;j;k ðx i;j;k Þ is the second order interaction etc. Y is the random variable corresponding to the function output, EðYÞ is the expected value of the output considering all possible combinations of inputs. The main effect of a parameter is the output of the model with the parameter held constant, averaged over all of the other parameters' possible values. This can be plotted over the parameter's possible range and gives a good visual representation of the model's sensitivity to that parameter. A plot of the interactions shows the effect of varying two or more parameters simultaneously (in addition to their main effects) averaged over the rest of the parameter space.

Variance based measures
The variance of the main effect is known as the main effect index (MEI) and it can be written as, This is the expected amount that the uncertainty of the model output would be reduced if the true value of x i was known.
The total sensitivity index (TSI) is the variance caused by a parameter and any interaction involving that parameter, It can also be thought of as the remaining variance if the true values of all of the parameters except x i are known (Ài refers to the complement of the subset i).

The first run of sensitivity analysis
When performing global sensitivity analysis, choosing the shape and width of the parameter distributions is important. Choosing an incorrectly wide distribution means that a parameter's importance could be overestimated and an incorrectly narrow one would underestimate importance. All of the parameters are assumed to have uniform distributions.
The purpose of the first run of sensitivity analysis performed here is to investigate how the uncertainty in individual model parameters is responsible for the uncertainty of the model output. This can be used to decide which parameters are most important when developing the model, as well as giving insight into how the system behaves. There are 25 parameters which were investigated here out of 28 in total (the mass of the rods, the gear ratios and the gravitational constant are all known accurately and are not going to change). 18 of the parameters were measured from 3D drawings of the system which were made to a tolerance of 1mm, this tolerance was used to estimate the parameter ranges. Probability distributions for the other 7 parameters were estimated using Bayesian system identification, however these distributions were estimated assuming that the other parameter values were accurate, so the ranges used have been doubled.
At this point it is worth noting the difference between subjective and objective uncertainty in parameter values. Subjective uncertainty results from a lack of accurate knowledge of the system, e.g. a dimension which is not accurately known. Objective uncertainty results from the fact that some elements in the system behave in a stochastic way, for instance, brake pad friction coefficients have been shown to vary unpredictably (Ostermeyer, 2003).
The MEIs and TSIs of the parameters which were responsible for more than 1% of the output variance are shown in Table 2. It can be seen that more than half of the output variance arises from the brake friction term l. It is also clear from the table that the vast majority of the variance arises from the main effects of parameters, since the values of the MEIs are close to the values of the TSIs. The sum of the MEIs is 95%, so interactions account for only around 5% of the output variance.
Plots of the main effects for selected parameters are shown in Fig. 5-8. Alongside these are plots of the model output with all of the parameters held at their expected values, except the parameter of interest which is varied across a range of values. It should be noted that the y-axis limits are not the same on the main effects plots and the corresponding one-at-a-time (1AAT) plots. A common theme across all of the parameters is that the expected model outputs from the main effects plots are higher (the rod has inserted further) than the corresponding model output from the 1AAT plots. This is because, on average, the model is more sensitive to the change in a parameter when it increases the distance the rod inserts than when it decreases it. This can be seen in Fig. 7 which shows that the output is more sensitive to a decrease in brake friction than an increase. The fact that the system is generally less sensitive to parameters when they slow down the rod's insertion suggests that the system is more likely to remain safe, but it is not the case for all parameters. Fig. 8 shows that plotting the main effects can obscure a local nonlinearity in the response to a change in a parameter. While in the current case the model is fairly insensitive to the parameter, it does highlight the fact that when looking at the global behaviour it is possible to miss details in the local behaviour.

The second run of sensitivity analysis
The purpose of the second run was to investigate what could happen if parameters were to change. Only four parameters were investigated; the brake friction coefficient l, the general friction term F f , the main spring stiffness and the reaction spring stiffness. These parameters were chosen because they could feasibly change during the lifetime of the reactor, and because the model was not shown to be totally insensitive to them during the first run of sensitivity analysis. The authors currently do not have quantitative information regarding how these parameters could change, so the ranges chosen are speculative. Each parameter was varied, from its expected value in the direction that would have a detri-mental effect on the systems performance, i.e. which would slow the speed of the rod's insertion. Table 3 shows the MEIs and TSIs from the second run of sensitivity analysis where fewer parameters were considered with an extended range. Again it can be seen that there is a relatively small contribution from the interactions. The results suggest that the most influential parameters are the brake friction coefficient and the main spring stiffness. However, it is not known how likely these parameters are to change, and by how much; it is not possible to truly state which are the most influential parameters without this information.
Figs. 9-12 show the main effects and 1AAT plots. The main effects show the rod inserting less far than the 1AAT plots. This      is to be expected since the main effects are the outputs averaged over the other parameters ranges, and the other parameters are varied from their expected values in the direction which slows the rod's insertion. The main effects plots also show much lower sensitivity to each of the parameters compared to the 1AAT plots. This shows that if parameters change causing the rod insertion to become slower, then the system will become less sensitive to other parameters changing. This suggests that unless a component drastically fails, the system is likely to stay safe.

Conclusions
The aim of this paper was to develop a model of the AGR primary shutdown system and explore the use of probabilistic sensitivity analysis techniques on this model, with the objective of shedding light on the contribution of parameters to model uncertainty and the system's performance.
The results show that the majority of model uncertainty comes from the brake friction, this is because the value of the friction coefficient is uncertain and the model is sensitive to changes in brake friction. However, the brake friction is an objectively uncertain parameter, it can vary due to changes in reactor conditions or wear of the brake disks. It is intended to further develop the model as a diagnostic tool, which interprets test data in order to monitor the values of the objectively uncertain parameters. In order to refine the model for this purpose, only the accuracy of the subjectively uncertain parameters can be improved. In this case the mass and dimensions of the flyweights are the most important parameters and their values should be verified, whereas parameters such as the return spring stiffness and the flyweight rotational inertia have very little effect on the model output and require no further investigation.
The results from the second round of sensitivity analysis imply that the system ought to be resistant to changes in several parameters at once; a component's properties would have to change dramatically before the system becomes unsafe. Extension of operating lives is a strategy in use to optimise the AGR fleet as part of the UK energy mix, and some initial periods of life extension have already been implemented. These results demonstrate the robustness of the primary shutdown system, suggesting that control rod system is unlikely to be a limitation for any further life extension.
It is made clear in these results that while these global sensitivity analysis techniques use information from the entire range of the possible parameter space, they do not fully describe the model's sensitivity to these parameters, because details (such as nonlinearities) can be lost when averaging over the other possible parameters. This study has shown that it can be useful to show local sensitivity analysis results alongside the global ones as it provides context for comparison.