Applying Machine Learning to Nuclear Reactor Kinetics
A primer on nuclear physics, reactor operation, and predictive modeling using Python
According to the Nuclear Regulatory Commission, there are 31 nuclear research reactors in the United States. I happen to have a license to operate one of them, and in this article I’ll be demonstrating how I applied machine learning and general data analytics techniques to predict pulse power levels and increase the repeatability of our experiments.
Background
A fission nuclear reactor works by harnessing the power of fissile atoms. When uranium-235 absorbs a neutron it has a chance to fission and split, releasing fission products, neutrons, and kinetic energy. This energy heats a coolant medium which is usually piped to a heat exchanger and then a steam turbine, which generates electricity. My facility is home to a TRIGA reactor that doesn’t produce any electricity — it’s used purely for research and experimentation.
Fun fact: one kilogram of U-235 has about 3 million times more energy than one kilogram of coal. Oh, and the fission reaction produces no carbon emissions. (I’m not biased at all.)
A TRIGA reactor is unique in many ways when compared to a commercial reactor, one of which being its capability to perform a “pulse”. The highly negative reactivity coefficient of the fuel means that as temperature increases, the reactivity — and thus the rate of fission chain reactions — decreases. This means that the reactor is self-limiting in terms of power level and, due to the design of the fuel, will physically shut itself down after a pulse with no input from the operator. (That’s why they let people like me operate it!)
A pulse works by pneumatically (using pressurized air) ejecting one of the control rods from the core of the reactor, which causes a rapid increase in power level. The following happens in the span of a few milliseconds:
- The control rod is shot vertically upwards out of the core.
- The power level increases from around 50 W to up to 2000 MW.
- At this high power level, the prompt negative feedback effect of the fuel provides negative reactivity to the core, which shuts itself down.
At this point the control rod may still be on its way out of the core but will subsequently fall back down due to gravity. You end up with a power response function P(t) that looks something like Figure 2 below — power increases as the control rod is ejected, then rapidly decreases with the addition of negative reactivity.
While pulsing has numerous research applications, it works especially well to simulate, on a very small scale, the radiation emitted from a nuclear explosion. Specifically, we’re observing how neutrons and gamma rays interact with electronics on an atomic level. Say you were designing a new electronic component to be used to control some portion of a nuclear weapon system. How certain are you that your electronics will survive a nearby nuclear detonation? What damage will the integrated circuitry sustain by sitting next to a radioactive warhead in storage for 10 years? How will that damage affect the functionality of your component? Or what if, instead, you’re designing a CPU chip to be used in a new fighter jet. Is the flight control system that’s enabled by your CPU going to fail when exposed to a certain amount of radiation? You can see how these questions are of serious importance. Combine that with the fact that a true nuclear detonation is a bit of an ordeal to set up and has a number issues (not to mention banned by the Comprehensive Nuclear-Test-Ban Treaty) and you’ll soon agree that reactor pulsing is a key capability.
NOTE: It’s worth elucidating that this pulse is only generating radiation. In this type of experiment, nothing is being exploded. We’re simply exposing a sample to high levels of radiation in an extremely controlled environment that’s designed for such a task.
The maximum power level reached by a pulse can be chosen by the operator based on the reactivity “worth” of the control rods. For many nuclear reactors this worth is measured in units of dollars, and since it’s outside of the scope of this article to explain why we use this seemingly weird unit, you can read all about it here if you wish.
Depending on the experiment at hand, a specific pulse power level might be desired. The positions of the control rods in the core determine the power level at any given time. The dataset that I’ve compiled from my own work at this reactor consists of these values and a handful of others, described below:
- Date. The specific datetime formatted date that the pulse occurred on.
- Estimated Reactivity. The estimated reactivity insertion amount, in dollars, of the pulse. This estimate is found by consulting the integral worth of a certain control rod. Example: an experimenter will request a pulse worth $2.00. The operator will find the position to place the control rod at based on the total integral worth of said control rod.
- Rod Positions. There are four control rods in the core — transient, shim 1, shim 2, and the regulating rod. Labeled “Trans, S1, S2, Reg” in my dataset, these values are related to the physical location of the control rods in the core. They range from 0 to 960, where 0 is fully inserted into the core and 960 is fully removed.
- Peak Power. Measured in MW.
- Total Energy. Measured in MW seconds.
- Peak Temperature. Measured in degrees Celsius, this is the peak temperature reached in an instrumented fuel element (IFE). A few fuel rods have thermocouples embedded inside of them to monitor the temperature of the core.
- Calculated Reactivity. This is the “true” reactivity value, which generally differs from the estimated reactivity by a certain amount. It’s calculated automatically by the reactor console.
In the rest of this post, I’ll analyze the dataset to see what insights I can glean. In addition, I’ll apply a linear regression machine learning model to predict calculated reactivity based on estimated reactivity and rod positions. This predictive model will help us to be more accurate with our pulse powers (to a certain extent — more thoughts on that later) and increase the efficacy and veracity of our experiments. My end goal is to make our pulses more accurate and repeatable in terms of the amount of radiation exposure given to a specific sample.
Exploratory Data Analysis (EDA)
I went ahead and cleaned the data beforehand, so I’ll spare you the details. In most cases, the cleaning process consists of removing incorrect or incomplete entries from the dataset. There were a lot of manual entry errors (typos) that I had to identify and cleanse, as well as some missing data points that I either completely removed or replaced with average values.
Once I have my clean data, the first thing I’ll do is import the relevant libraries and load up my DataFrame:
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as npfrom sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn import metrics
from sklearn.metrics import r2_scoredf = pd.read_excel('Pulse_Data_2021_NO_NULL.xlsx', )
After that, I always like to use three main functions to learn the general characteristics of my data — .head(), .info(), and .describe().
.head() just shows the top few rows of the DataFrame, and we can see how the data is generally structured:
.info() outputs the number of entries in the DataFrame, the columns names, the number of null entries in each column, and the data type for each column:
And .describe() provides some summary analysis of the data itself — mean, standard deviation, quartiles, etc. of each column. I like to transpose it for readability:
The next step of my EDA is to just start graphing whatever might make sense. For instance, I thought it would be interesting to see the relationship between calculated reactivity and peak power, so I made a simple scatterplot:
There’s a clear exponential relationship between the two variables which is consistent with the principles of reactor kinetics that drive this operation. The graph is also useful in that it identifies potential outlier points — either weird behaviors of the reactor or, more likely, incorrect entries. From this graph, I found additional fat-fingered entries that needed to be cleaned or removed.
A correlation heatmap is another useful tool that can be applied to most numerical datasets. Figure 7 below depicts which columns are closely correlated to each other, with values ranging from -1 to 1.
To interpret this heatmap, look for values that are high (close to one). The diagonal line of ones is an artifact of the layout of the heatmap, and just tells me that each feature is directly correlated to itself. Makes sense. S1 and S2 have a correlation of 0.99, which also makes sense because these two control rods are almost always in the exact same vertical position when pulsing. The feature I’m most interested in, calculated reactivity, is highly correlated with peak temperature, total energy, and peak power. This also seems logical, because as reactivity increases, I expect each of those features to increase as well.
Since I’m trying to improve the predictive accuracy of the estimated reactivity value, my next step is to compare a distribution of it versus the calculated reactivity. To accomplish this, I overlay two histograms and do some formatting to make it look nice:
There are discrete pulse reactivity values that we generally use when deciding how “large” of a pulse we want. Figure 8 reflects this clearly, showing that $1.50, $2.00, $2.50, and $3.00 are common estimated values. One would expect a somewhat normal distribution of calculated reactivity values around each of the estimated reactivity values, which is shown (albeit loosely) by the blue graph.
Again looking at Figure 8, it appears that the estimated reactivity is skewed slightly higher than the calculated counterpart. What this means is that generally speaking, if you request a $2.00 pulse, you’ll actually get a value a bit less than that. I can quantify this by simply subtracting the two columns of data and finding the mean:
Which demonstrates that, on average, the “true” calculated reactivity is $0.16 lower than the estimated reactivity.
Lastly, I generated a graph that serves as an interesting proxy for reactor utilization over time. This reactor first started up (AKA “went critical”) in 1992, and a graph of the number of pulses per year since then offers some deep rabbit holes to go down:
Questions that could arise from Figure 10 include:
- Why were there so few pulses from 1994–1996, and from 2013–2014?
- Was there an administrative change that affected the type of experiments run at this facility?
- Were there any new national or university research developments that required additional pulses?
- What experiments were happening in 2000 and 2020–20201 (now) that required so many pulses, and why were there no such experiments in between?
It’s amazing how much insight can be garnered from a few small graphs.
Predictive Modelling
The task of predicting the calculated reactivity is fit (pun intended) for a linear regression model. This is considered a supervised machine learning model because the data is already labeled (x and y values are provided to train the model). It’s technically a multiple regression model (the definition of which encompasses linear regression) because it uses multiple independent variables (Est_Reactivity, Trans, S1, S2, Reg) to predict the outcome of the dependent variable (Calc_Reactivity).
NOTE: this linear regression model assumes that the relationship between the dependent variable and the independent variables are linear. Estimated reactivity is clearly linearly related to calculated reactivity, but the same cannot be said about our other independent variables (Trans, S1, S2, Reg). There is also clear multicollinearity among the independent variables, as rod placements (S1, S2, Reg) are highly correlated. For these reasons, a linear regression model may not be the best option. Further exploration of the independent variables is needed to know for sure, but polynomial regression or ridge regression might be a good fit. The bottom line is that it’s important to understand the assumptions and limitations of your model and how they vary across different techniques.
I first split the data into a training set and testing set to ensure an unbiased evaluation of the model. I then instantiate the model and fit the data:
Once the model has been fit, I use the testing data to compare the model output against the expected values. A completely straight line would indicate a perfect model:
With a few notable outliers, the model does a great job of accurate prediction. Depending on the use case of the model, it may be worth investigating the data even further to determine where these outliers are coming from and how to mitigate them to improve model accuracy. In my case, this model is sufficient. I can even quantify the “goodness of fit” by calculating the R-squared value for my model:
R-squared is the proportion of variance in the dependent variable that’s predictable from the independent variable. Expressed as a percentage, 91% of the variance in the estimated reactivity of any given pulse can be explained with our input values. The takeaway is that this is a great model for our purposes. That being said, there’s one more check I can do to learn about the model. Figure 14 to the left is a histogram of residuals, or the distance between any given data point and the regression line of best fit. The graph indicates that there is a normal distribution of the random error around zero, which is good. If this were not the case, then we might have an issue with our model and/or dataset. Again there are notable outliers, but nothing to be concerned about when backed up with our impressive R-squared score.
Conclusions (and limitations)
With this model, I can now accurately predict a pulse reactivity value based on reactor parameters such as control rod positions and rod worth. Repeatability is important for all types of experimentation, and this model will help tighten our pulse values and ensure similar radiation exposure for each component being irradiated.
NOTE: The model does have a number of limitations that might be obvious to those with reactor kinetics and operations knowledge. Namely, the model has no knowledge of previous operations and the extremely important buildup of fission product poisons. Xenon-135 is created (and burned off) over time as the reactor operates and significantly affects neutron absorption and reactor behavior.
EDIT: After testing the model with brand new data before a pulse, the predicted reactivity was less than 3% off from the “true” reactivity calculated by the reactor console. Success!
If you enjoyed this post, feel free to check out some of my others and let me know what you think in the comments. Thanks for reading!
