HTGR Micro-Core Quadrant Power

Input

  • theta1: Angle of control drum in quadrant 1 (radians)

  • theta2: Angle of control drum in quadrant 1 (radians)

  • theta3: Angle of control drum in quadrant 2 (radians)

  • theta4: Angle of control drum in quadrant 2 (radians)

  • theta5: Angle of control drum in quadrant 3 (radians)

  • theta6: Angle of control drum in quadrant 3 (radians)

  • theta7: Angle of control drum in quadrant 4 (radians)

  • theta8: Angle of control drum in quadrant 4 (radians)

Output

  • fluxQ1 : Neutron flux in quadrant 1 (\(\frac{neutrons}{cm^{2} s}\))

  • fluxQ2 : Neutron flux in quadrant 2 (\(\frac{neutrons}{cm^{2} s}\))

  • fluxQ3 : Neutron flux in quadrant 3 (\(\frac{neutrons}{cm^{2} s}\))

  • fluxQ4 : Neutron flux in quadrant 4 (\(\frac{neutrons}{cm^{2} s}\))

The pyMAISE.datasets.load_HTGR data features 751 samples with eight inputs and five outputs. It is the raw data from [PRK22]. Using reactor symmetry, [PRK22] boosts the sample size to 3004.

The data set was built using MCNP simulations of the HOLOS-Quad reactor design, which is depicted below [PRK22]. This reactor implements modular construction where separate units can be transported independently and assembled at the plant. The HOLOS-Quad core is a 22 MWt high-temperature gas-cooled microreactor (HTGR) controlled by eight cylindrical control drums. It utilizes TRISO fuel particles contained in hexagonal graphite blocks as a moderator. These graphite blocks have channels where helium gas can pass through for cooling.

Using machine learning (ML), we aim to predict the neutron distribution given the angles of each of the control drums \(\in[0, 2\pi]\). The drums control reactivity by rotating to vary the proximity of \(B_4C\), indicated as the red strip on the drums shown below, to the fuel. Perturbations in the control drum angle can drastically affect the power distribution in the core. Therefore, predictions of control drum reactivity worth for arbitrary configurations make this problem nontrivial.

microreactor_imag.png

[1]:

# Importing Packages
import time
import numpy as np
import pandas as pd

# Set display option to show all rows and columns
pd.set_option('display.max_rows', None)
pd.set_option('display.max_columns', None)

# Set the width of the columns
pd.set_option('display.width', None)

# See the full content of each column
pd.set_option('display.max_colwidth', None)

import xarray as xr
import matplotlib.pyplot as plt
from scipy.stats import uniform, randint
from sklearn.preprocessing import Normalizer, MinMaxScaler

# pyMAISE specific imports
import pyMAISE as mai
from pyMAISE.datasets import load_HTGR
from pyMAISE.preprocessing import scale_data, train_test_split, correlation_matrix

# Plot settings
matplotlib_settings = {
    "font.size": 12,
    "legend.fontsize": 11,
    "figure.figsize": (8, 8)
}
plt.rcParams.update(**matplotlib_settings)

Training ML Models on the Raw Data (751 Samples)

Before applying symmetry to boost the samples, we evaluate the performance of several ML models on the original data set.

Preprocessing

We define a regression problem and load the data using pyMAISE.datasets.load_HTGR.

[2]:

# Initializing pyMAISE settings
global_settings = mai.init(
    problem_type=mai.ProblemType.REGRESSION,   # Define a regression problem
    cuda_visible_devices="-1"                  # Use CPU only
)

# Get data
data, inputs, outputs = load_HTGR()

The data consists of 751 samples with 8 inputs:

[3]:

inputs

[3]:

<xarray.DataArray (index: 751, variable: 8)>
array([[5.91952603, 2.36950271, 2.92365595, ..., 4.00890489, 4.97036758,
        2.98796636],
       [2.16238049, 0.273624302, 0.927740883, ..., 0.170166579,
        2.12404792, 4.98020926],
       [0.450100101, 0.006300644, 2.51221749, ..., 3.58225182,
        0.280763918, 4.88859544],
       ...,
       [5.96700387, 5.38006724, 5.20081594, ..., 5.65905406, 5.88623681,
        1.86231118],
       [1.34959071, 1.19391098, 4.65205108, ..., 3.19313598, 5.31477415,
        3.84849236],
       [3.60605462, 5.47085601, 1.75207802, ..., 5.07434654, 2.41128238,
        5.64650044]], dtype=object)
Coordinates:
  * index     (index) int64 0 1 2 3 4 5 6 7 ... 743 744 745 746 747 748 749 750
  * variable  (variable) object 'theta1' 'theta2' 'theta3' ... 'theta7' 'theta8'

and 4 outputs:

[4]:

outputs

[4]:

<xarray.DataArray (index: 751, variable: 4)>
array([[2.58e+19, 2.59e+19, 2.67e+19, 2.56e+19],
       [2.55e+19, 2.53e+19, 2.51e+19, 2.51e+19],
       [2.57e+19, 2.58e+19, 2.52e+19, 2.52e+19],
       ...,
       [2.57e+19, 2.53e+19, 2.52e+19, 2.57e+19],
       [2.64e+19, 2.52e+19, 2.59e+19, 2.61e+19],
       [2.51e+19, 2.55e+19, 2.61e+19, 2.55e+19]], dtype=object)
Coordinates:
  * index     (index) int64 0 1 2 3 4 5 6 7 ... 743 744 745 746 747 748 749 750
  * variable  (variable) object 'fluxQ1' 'fluxQ2' 'fluxQ3' 'fluxQ4'

Below is the correlation for this data set.

[5]:

correlation_matrix(data)
plt.show()
../_images/benchmarks_HTGR_microreactor_9_0.png

We see there are no noticeable linear correlations between the inputs and outputs. There are some slightly positive correlations for theta6 with fluxQ3 and theta7 with fluxQ4.

We train/test split the data (70%/30%) and scale the inputs and outputs. The inputs are min-max scaled, and the outputs are normalized across each sample.

[6]:

# Train test split data
train_data, test_data = train_test_split(data=data, test_size=0.30)

theta_cols = [f"theta{i + 1}" for i in range(8)]
flux_cols = [f"fluxQ{i + 1}" for i in range(4)]

# Min-max scaling inputs and normalize and min-max scale outputs
xtrain, xtest, _ = scale_data(train_data.loc[:, theta_cols], test_data.loc[:, theta_cols], MinMaxScaler())
ytrain, ytest, _ = scale_data(train_data.loc[:, flux_cols], test_data.loc[:, flux_cols], Normalizer(norm="l1"))

Model Initialization

We initialize ten regression models:

  • Linear: Linear,

  • Lasso: Lasso,

  • Ridge: RD,

  • ElasticNet: EN,

  • Decision tree: DT,

  • Random forest: RF,

  • K-nearest neighbors: KN,

  • Stacking: Stacking,

  • Gradient Boosting: GB

  • Feedforward dense neural: FNN.

All the classical models are initialized as Scikit-learn defaults, except Stacking which requires you to provide at least one model for an estimator. We provided AdaBoosting and ElasticNet as our models, setting multi-output on. We also have to set multi-output on for Gradient Boosting since it is not supported natively. The FNN features a dense hidden network with three sections: input, hidden, and output. We hyperparameter tune the number of hidden layers and the number of nodes and sublayers for both the input and dense layers. We also tune the learning rate of the Adam optimizer and the batch size. We set epochs to 50 but will increase this after hyperparameter tuning.

[7]:

# Initializing classical models and define hyperparameter search space for FNN
model_settings = {
    "models": ["Linear", "Lasso", "Stacking", "DT", "RF", "KN", "GB", "RD", "EN", "FNN"],
    "FNN": {
        "structural_params": {
            "Dense_hidden": {
                "num_layers": mai.Int(min_value=0, max_value=2),
                "units": mai.Int(min_value=25, max_value=400),
                "activation": "relu",
                "kernel_initializer": "normal",
                "sublayer": mai.Choice(["Dropout_hidden", "None"]),
                "Dropout_hidden": {
                    "rate": mai.Float(min_value=0.2, max_value=0.6),
                },
            },
            "Dense_output": {
                "units": ytrain.shape[-1],
                "activation": "linear",
                "kernel_initializer": "normal",
            },
        },
        "optimizer": "Adam",
        "Adam": {
            "learning_rate": mai.Float(min_value=1e-5, max_value=0.001),
        },
        "compile_params": {
            "loss": "mean_absolute_error",
            "metrics": ["mean_absolute_error"],
        },
        "fitting_params": {
            "batch_size": mai.Choice([8, 16, 32]),
            "epochs": 50,
            "validation_split": 0.15,
        },
    },
    "GB":{
        "multi_output": True,
    },
    "Stacking":{
        "estimators": ["KN", "SVM"],
        "_multi_output": True,
    },
}
# Construct pyMAISE.Tuner object
tuner = mai.Tuner(xtrain, ytrain, model_settings=model_settings)

Hyperparameter Tuning

Below we define the search spaces for all the classical models except linear regression. For linear regression, we will use only the Scikit-learn default configuration. We utilize a random search with 200 iterations and 5-fold cross-validation, given their training speed. Given that the FNNs are more computationally expensive to train, we use a Bayesian search with 50 iterations, each with 5-fold cross-validation. We hope to converge on a performant model based on the Bayesian search.

[8]:

# Classical model search space
random_search_spaces = {
    "Lasso": {
        "alpha": uniform(loc=0.0001, scale=0.0099), # 0.0001 - 0.01
    },
    "DT": {
        "max_depth": randint(low=5, high=50), # 5 - 50
        "max_features": [None, "sqrt", "log2", 2, 4, 6],
        "min_samples_split": randint(low=2, high=20), # 2 - 20
        "min_samples_leaf": randint(low=1, high=20), # 1 - 20
    },
    "RF": {
        "n_estimators": randint(low=50, high=200), # 50 - 200
        "criterion": ["squared_error", "absolute_error", "poisson"],
        "min_samples_split": randint(low=2, high=20), # 2 - 20
        "min_samples_leaf": randint(low=1, high=20), # 1 - 20
        "max_features": [None, "sqrt", "log2", 2, 4, 6],
    },
    "KN": {
        "n_neighbors": randint(low=1, high=20), # 1 - 20
        "weights": ["uniform", "distance"],
        "leaf_size": randint(low=1, high=30), # 1 - 30
        "p": randint(low=1, high=10), # 1 - 10
    },
    "RD":{
        "alpha": uniform(loc=0.0001, scale=0.0099), # 0.0001 - 0.01
    },
    "EN":{
        "alpha": uniform(loc=0.0001, scale=0.0099), # 0.0001 - 0.01
    },
    "GB": {
        "estimator__learning_rate": uniform(loc=0.1, scale=1e-2)
    },
    "Stacking":{
        "estimator__KN__n_neighbors": randint(low=1, high=20), # 1 - 20
        "estimator__KN__weights": ["uniform", "distance"],
        "estimator__KN__leaf_size": randint(low=1, high=30), # 1 - 30
        "estimator__KN__p": randint(low=1, high=10), # 1 - 10
        "estimator__SVM__C": uniform(loc=0.01, scale=99.99), # 0.01 - 100

    },
}

start = time.time()

# Run classical model random search
random_search_configs = tuner.random_search(
    param_spaces=random_search_spaces,
    n_iter=200,
    n_jobs=6,
    cv=5,
)

# Run FNN Bayesian search
bayesian_search_configs = tuner.nn_bayesian_search(
    objective="r2_score",
    max_trials=50,
    cv=5,
)
print("Hyperparameter tuning took " + str((time.time() - start) / 60) + " minutes to process.")

Hyperparameter tuning took 9.143632805347442 minutes to process.

We can visualize the hyperparameter tuning with Bayesian search of the FNN from the convergence plot below.

[9]:

fig, ax = plt.subplots()
ax.set_ylim([-10, 1])
ax = tuner.convergence_plot(model_types="FNN")
plt.show()
../_images/benchmarks_HTGR_microreactor_17_0.png

The Bayesian search finds several configurations close to 0.4 validation \(R^2\).

Model Postprocessing

With the top 5 configurations saved, we can pass these to the pyMAISE.PostProcessor for model comparison and analysis. For the FNNs, we increase the epochs to 300 for better performance.

[10]:

new_model_settings = {"FNN": {"fitting_params": {"epochs": 300}}}

postprocessor = mai.PostProcessor(
    data=(xtrain, xtest, ytrain, ytest),
    model_configs=[random_search_configs, bayesian_search_configs],
    new_model_settings=new_model_settings,
)

Below we computed the performance metrics based on the performance of each model on all quadrants.

[11]:

postprocessor.metrics().drop("Parameter Configurations", axis=1)

[11]:
Model Types Train R2 Train MAE Train MAPE Train RMSE Train RMSPE Test R2 Test MAE Test MAPE Test RMSE Test RMSPE
29 GB 0.983467 0.000413 0.164976 0.000540 0.214139 0.911719 0.000986 0.393361 0.001248 0.495324
26 GB 0.983750 0.000410 0.163706 0.000535 0.211759 0.910683 0.000999 0.398571 0.001259 0.500008
30 GB 0.983528 0.000412 0.164426 0.000538 0.213496 0.910662 0.000992 0.395955 0.001259 0.499872
28 GB 0.983264 0.000416 0.166319 0.000543 0.215309 0.909914 0.000996 0.397611 0.001263 0.501404
27 GB 0.983881 0.000409 0.163136 0.000532 0.211078 0.909896 0.001000 0.399252 0.001263 0.501905
42 FNN 0.885324 0.001157 0.462625 0.001423 0.560177 0.841728 0.001316 0.526366 0.001661 0.658946
44 FNN 0.876756 0.001222 0.488402 0.001493 0.570846 0.812901 0.001425 0.568708 0.001784 0.692570
45 FNN 0.788750 0.001604 0.642872 0.001949 0.755871 0.731481 0.001703 0.682048 0.002129 0.830303
41 FNN 0.793487 0.001543 0.617249 0.001899 0.722030 0.690805 0.001886 0.755055 0.002367 0.922341
43 FNN 0.777917 0.001699 0.678816 0.001994 0.771722 0.685179 0.001908 0.761591 0.002326 0.913173
17 RF 0.891503 0.001099 0.438799 0.001378 0.547905 0.625088 0.002050 0.818008 0.002584 1.023890
16 RF 0.936638 0.000844 0.337117 0.001052 0.418109 0.624696 0.002053 0.819313 0.002584 1.024092
19 RF 0.876376 0.001169 0.466503 0.001471 0.584289 0.614802 0.002086 0.832203 0.002617 1.037431
20 RF 0.876483 0.001175 0.468892 0.001472 0.585309 0.609549 0.002105 0.839952 0.002637 1.045834
18 RF 0.910339 0.000994 0.397080 0.001251 0.497497 0.609524 0.002095 0.836107 0.002635 1.044466
24 KN 1.000000 0.000000 0.000000 0.000000 0.000000 0.605093 0.002088 0.833138 0.002652 1.052266
22 KN 1.000000 0.000000 0.000000 0.000000 0.000000 0.599978 0.002099 0.837563 0.002668 1.058594
23 KN 1.000000 0.000000 0.000000 0.000000 0.000000 0.599978 0.002099 0.837563 0.002668 1.058594
21 KN 1.000000 0.000000 0.000000 0.000000 0.000000 0.597157 0.002115 0.843946 0.002677 1.062344
25 KN 1.000000 0.000000 0.000000 0.000000 0.000000 0.594733 0.002127 0.848653 0.002686 1.065676
15 DT 0.557221 0.002198 0.877738 0.002785 1.107820 0.337739 0.002738 1.093277 0.003423 1.361120
14 DT 0.567213 0.002192 0.875144 0.002754 1.094594 0.310598 0.002755 1.100043 0.003467 1.378080
13 DT 0.623282 0.002014 0.804163 0.002568 1.020332 0.290526 0.002809 1.121057 0.003551 1.407312
12 DT 0.547270 0.002218 0.885596 0.002817 1.120458 0.237935 0.002905 1.159861 0.003679 1.460174
11 DT 0.660487 0.001948 0.777785 0.002436 0.968025 0.228219 0.002924 1.167367 0.003727 1.465348
31 RD 0.259911 0.002925 1.167418 0.003613 1.436461 0.225834 0.002967 1.184048 0.003699 1.469195
32 RD 0.259911 0.002925 1.167418 0.003613 1.436461 0.225834 0.002967 1.184048 0.003699 1.469195
33 RD 0.259911 0.002925 1.167418 0.003613 1.436461 0.225834 0.002967 1.184048 0.003699 1.469195
34 RD 0.259911 0.002925 1.167418 0.003613 1.436461 0.225834 0.002967 1.184048 0.003699 1.469195
35 RD 0.259911 0.002925 1.167418 0.003613 1.436461 0.225833 0.002967 1.184048 0.003699 1.469195
0 Linear 0.259911 0.002925 1.167418 0.003613 1.436462 0.225821 0.002967 1.184056 0.003699 1.469206
36 EN 0.217982 0.003000 1.197600 0.003713 1.476293 0.200210 0.003012 1.201854 0.003765 1.494888
1 Lasso 0.211628 0.003011 1.201968 0.003728 1.482299 0.195644 0.003021 1.205523 0.003776 1.499379
37 EN 0.205563 0.003022 1.206181 0.003742 1.488004 0.191287 0.003029 1.208842 0.003787 1.503648
38 EN 0.190648 0.003048 1.216575 0.003777 1.501873 0.180971 0.003048 1.216431 0.003813 1.513724
2 Lasso 0.187193 0.003054 1.219003 0.003785 1.505076 0.178543 0.003053 1.218215 0.003819 1.516071
3 Lasso 0.167250 0.003087 1.232355 0.003831 1.523457 0.162762 0.003082 1.229799 0.003856 1.530998
4 Lasso 0.157135 0.003105 1.239255 0.003855 1.532682 0.154534 0.003097 1.235927 0.003876 1.538642
39 EN 0.111916 0.003180 1.269452 0.003956 1.573006 0.115404 0.003172 1.265967 0.003965 1.574074
5 Lasso 0.102890 0.003195 1.275457 0.003976 1.580861 0.106665 0.003188 1.272305 0.003985 1.581847
40 EN 0.099469 0.003201 1.277730 0.003983 1.583814 0.103271 0.003195 1.274838 0.003993 1.584854
6 Stacking 0.011404 0.003350 1.337361 0.004174 1.659549 0.006224 0.003369 1.344298 0.004205 1.668890
7 Stacking 0.011404 0.003350 1.337361 0.004174 1.659549 0.006224 0.003369 1.344298 0.004205 1.668890
8 Stacking 0.010531 0.003352 1.337941 0.004176 1.660274 0.005244 0.003370 1.344865 0.004207 1.669714
10 Stacking 0.008668 0.003355 1.339182 0.004180 1.661836 0.005000 0.003371 1.345281 0.004208 1.669944
9 Stacking 0.008668 0.003355 1.339182 0.004180 1.661836 0.005000 0.003371 1.345281 0.004208 1.669944

While the most performant models are the FNNs and the GB’s, all models performed poorly on this data set. Most models were overfit, with the FNNs offering the best performance. Below are the hyperparameter configurations.

[12]:

for model in ["Linear", "Lasso", "DT", "RF", "KN", "GB", "RD", "EN", "Stacking", "FNN"]:
    postprocessor.print_model(model_type=model)
    print()

Model Type: Linear
  copy_X: True
  fit_intercept: True
  n_jobs: None
  positive: False

Model Type: Lasso
  alpha: 0.00010412251455299124

Model Type: DT
  max_depth: 26
  max_features: 6
  min_samples_leaf: 11
  min_samples_split: 2

Model Type: RF
  criterion: squared_error
  max_features: log2
  min_samples_leaf: 2
  min_samples_split: 4
  n_estimators: 130

Model Type: KN
  leaf_size: 13
  n_neighbors: 7
  p: 1
  weights: distance

Model Type: GB
  estimator__learning_rate: 0.10881748050180497

Model Type: RD
  alpha: 0.009989838843071154

Model Type: EN
  alpha: 0.00019133574600236952

Model Type: Stacking
  estimator__KN__leaf_size: 23
  estimator__KN__n_neighbors: 1
  estimator__KN__p: 1
  estimator__KN__weights: distance
  estimator__SVM__C: 76.43781578179083

Model Type: FNN
  Structural Hyperparameters
    Layer: Dense_hidden_0
      units: 326
      sublayer: None
    Layer: Dense_output_0
  Compile/Fitting Hyperparameters
    Adam_learning_rate: 0.0008197275012268372
    batch_size: 8
Model: "FNN"
_________________________________________________________________
 Layer (type)                Output Shape              Param #
=================================================================
 Dense_hidden_0 (Dense)      (None, 326)               2934

 Dense_output_0 (Dense)      (None, 4)                 1308

=================================================================
Total params: 4242 (16.57 KB)
Trainable params: 4242 (16.57 KB)
Non-trainable params: 0 (0.00 Byte)
_________________________________________________________________

Below is the network plot for the FNN.

[13]:

postprocessor.nn_network_plot(
    to_file="./supporting/HTGR_microreactor_network_0.png",
    show_shapes=True,
    show_layer_names=True,
    expand_nested=True,
    show_layer_activations=True,
)

[13]:
../_images/benchmarks_HTGR_microreactor_25_0.png

To visualize the performance of these models we can use diagonal validation and validation plots as shown below.

[14]:

models = np.array([["Linear", "Lasso"], ["DT", "KN"], ["RF", "FNN"], ["GB", "RD"], ["Stacking", "EN"]])
fig, axarr = plt.subplots(models.shape[0], models.shape[1], sharex=True, sharey=True, figsize=(15,20))
for i in range(models.shape[0]):
    for j in range(models.shape[1]):
        plt.sca(axarr[i, j])
        axarr[i, j] = postprocessor.diagonal_validation_plot(model_type=models[i, j])
        axarr[i, j].set_title(models[i, j])

plt.show()
../_images/benchmarks_HTGR_microreactor_27_0.png 
[15]:

models = np.array([["Linear", "Lasso"], ["DT", "KN"], ["RF", "FNN"], ["GB", "RD"], ["Stacking", "EN"]])
fig, axarr = plt.subplots(models.shape[0], models.shape[1], sharex=True, sharey=True, figsize=(15,20))
for i in range(models.shape[0]):
    for j in range(models.shape[1]):
        plt.sca(axarr[i, j])
        axarr[i, j] = postprocessor.validation_plot(model_type=models[i, j])
        axarr[i, j].set_title(models[i, j])

plt.show()
../_images/benchmarks_HTGR_microreactor_28_0.png

The diagonal validation and validation plots agree with the conclusions gathered from the performance metrics. Each model overfit the training data, and none offers good performance. Below is the learning curve for the best-performing FNN, according to test \(R^2\).

[16]:

fig, ax = plt.subplots()
ax = postprocessor.nn_learning_plot()

plt.show()
../_images/benchmarks_HTGR_microreactor_30_0.png

The top-performing FNN does not overfit the training data, as shown by the overlapping of the training and validation loss curves.

Increasing Samples Using Reactor Symmetry (3004 Samples)

We need more data to achieve better results. As shown in [PRK22], we can apply symmetry rules to multiply the sample size.

Preprocessing

To multiply our samples, we use the mult_sym and g21 functions from https://github.com/deanrp2/MicroControl/blob/main/pmdata/utils.py. These functions apply reflections over quadrant symmetries to multiply our samples by \(4\times\). These functions have been reimplemented and simplified in mult_samples using the Xarray data structures pyMAISE uses. Additionally, rather than moving the drum range to \([-\pi, \pi]\), we leave it at \([0, 2\pi]\) to keep consistent with the previous section.

[17]:

# Credit to mult_sym and g21 from https://github.com/deanrp2/MicroControl/blob/main/pmdata/utils.py#L51
def mult_samples(data):
    # Create empty arrays
    ht = xr.DataArray(
        np.zeros(data.shape),
        coords={
            "index": [f"{idx}_h" for idx in data.coords["index"].values],
            "variable": data.coords["variable"],
        },
    )
    vt = xr.DataArray(
        np.zeros(data.shape),
        coords={
            "index": [f"{idx}_v" for idx in data.coords["index"].values],
            "variable": data.coords["variable"],
        },
    )
    rt = xr.DataArray(
        np.zeros(data.shape),
        coords={
            "index": [f"{idx}_r" for idx in data.coords["index"].values],
            "variable": data.coords["variable"],
        },
    )

    # Swap drum positions
    hkey = [f"theta{i}" for i in np.array([3, 2, 1, 0, 7, 6, 5, 4], dtype=int) + 1]
    vkey = [f"theta{i}" for i in np.array([7, 6, 5, 4, 3, 2, 1, 0], dtype=int) + 1]
    rkey = [f"theta{i}" for i in np.array([4, 5, 6, 7, 0, 1, 2, 3], dtype=int) + 1]

    ht.loc[:, hkey] = data.loc[:, theta_cols].values
    vt.loc[:, vkey] = data.loc[:, theta_cols].values
    rt.loc[:, rkey] = data.loc[:, theta_cols].values

    # Adjust angles
    ht.loc[:, hkey] = (3 * np.pi - ht.loc[:, hkey].loc[:, hkey]) % (2 * np.pi)
    vt.loc[:, vkey] = (2 * np.pi - vt.loc[:, hkey].loc[:, vkey]) % (2 * np.pi)
    rt.loc[:, rkey] = (np.pi + rt.loc[:, hkey].loc[:, rkey]) % (2 * np.pi)

    # Fill quadrant tallies
    hkey = [2, 1, 4, 3]
    vkey = [4, 3, 2, 1]
    rkey = [3, 4, 1, 2]

    ht.loc[:, [f"fluxQ{i}" for i in hkey]] = data.loc[:, flux_cols].values
    vt.loc[:, [f"fluxQ{i}" for i in vkey]] = data.loc[:, flux_cols].values
    rt.loc[:, [f"fluxQ{i}" for i in rkey]] = data.loc[:, flux_cols].values

    sym_data = xr.concat([data, ht, vt, rt], dim="index").sortby("index")

    # Normalize fluxes
    sym_data.loc[:, flux_cols].values = Normalizer().transform(sym_data.loc[:, flux_cols].values)

    # Convert global coordinate system to local
    loc_offsets = np.array(
        [3.6820187359906447, 4.067668586955522, 2.2155167202240653 - np.pi, 2.6011665711889425 - np.pi,
         0.5404260824008517, 0.9260759333657285, 5.3571093738138575 - np.pi, 5.742759224778734 - np.pi]
    )

    # Apply correct 0 point
    sym_data.loc[:, theta_cols] = sym_data.loc[:, theta_cols] - loc_offsets + 2 * np.pi

    # Reverse necessary angles
    sym_data.loc[:, [f"theta{i}" for i in [3,4,5,6]]] *= -1

    # Scale all to [0, 2 * np.pi]
    sym_data.loc[:, theta_cols] = sym_data.loc[:, theta_cols] % (2 * np.pi)

    return sym_data

Rather than using the careful_split function in https://github.com/deanrp2/MicroControl/blob/main/pmdata/utils.py, we apply mult_samples to the training and testing data split, train_data and test_data, in the previous section. This ensures that reflected samples are in the same split as their original sample and keeps the split consistent across sections.

[18]:

sym_train_data = mult_samples(train_data)
sym_test_data = mult_samples(test_data)
print(f"Multiplied training shape: {sym_train_data.shape}, Multiplied testing shape: {sym_test_data.shape}")

Multiplied training shape: (2100, 12), Multiplied testing shape: (904, 12)

We see our total samples increased from 751 to 3004.

Consistent with [PRK22] and the previous section, we normalize the fluxes for each sample. Additionally, we min-max the inputs.

[19]:

# Min-Max scaling data
xtrain, xtest, _ = scale_data(sym_train_data.loc[:, theta_cols], sym_test_data.loc[:, theta_cols], MinMaxScaler())
ytrain, ytest, _ = scale_data(sym_train_data.loc[:, flux_cols], sym_test_data.loc[:, flux_cols], Normalizer(norm="l1"))

Model Initialization and Hyperparameter Tuning

Here, we tune the models with the same hyperparameter search spaces defined in the previous section. We also use the same search methods with consistent iterations and cross-validation.

[20]:

# Initialize pyMAISE.Tuner
tuner = mai.Tuner(xtrain, ytrain, model_settings=model_settings)

start = time.time()

# Run classical model random search
random_search_configs = tuner.random_search(
    param_spaces=random_search_spaces,
    n_iter=200,
    n_jobs=6,
    cv=5,
)

# Run FNN Bayesian search
bayesian_search_configs = tuner.nn_bayesian_search(
    objective="r2_score",
    max_trials=50,
    cv=5,
)
print("Hyperparameter tuning took " + str((time.time() - start) / 60) + " minutes to process.")

Hyperparameter tuning took 179.93817623853684 minutes to process.

Once again we can show the convergence of Bayesian search.

[21]:

fig, ax = plt.subplots()
ax.set_ylim([0, 1])
ax = tuner.convergence_plot(model_types="FNN")
plt.show()
../_images/benchmarks_HTGR_microreactor_40_0.png

The Bayesian search provides FNN configurations that perform noticeable better than the FNNs tuned with 751 samples.

Model Postprocessing

With the models tuned and the top 5 saved, we can now pass these models to the pyMAISE.PostProcessor for model comparison and analysis. Consistent with before we change the epochs of the FNNs to 300.

[22]:

postprocessor = mai.PostProcessor(
    data=(xtrain, xtest, ytrain, ytest),
    model_configs=[random_search_configs, bayesian_search_configs],
    new_model_settings=new_model_settings,
)

The average performance metrics are shown below.

[23]:

postprocessor.metrics().drop("Parameter Configurations", axis=1)

[23]:
Model Types Train R2 Train MAE Train MAPE Train RMSE Train RMSPE Test R2 Test MAE Test MAPE Test RMSE Test RMSPE
42 FNN 9.869901e-01 0.000372 0.148485 0.000480 0.190612 9.783814e-01 0.000489 0.195313 0.000621 0.247027
41 FNN 9.839465e-01 0.000426 0.170394 0.000533 0.212061 9.746946e-01 0.000533 0.213109 0.000672 0.267268
43 FNN 9.790173e-01 0.000476 0.190050 0.000609 0.240936 9.726791e-01 0.000543 0.216569 0.000698 0.276566
45 FNN 9.610011e-01 0.000667 0.266749 0.000831 0.325246 9.588864e-01 0.000691 0.276460 0.000856 0.336885
44 FNN 9.600366e-01 0.000657 0.262652 0.000841 0.303945 9.552796e-01 0.000707 0.282787 0.000893 0.330861
26 GB 9.662669e-01 0.000582 0.231956 0.000773 0.306755 9.337358e-01 0.000839 0.334627 0.001087 0.431100
28 GB 9.659397e-01 0.000587 0.233863 0.000776 0.308170 9.327820e-01 0.000845 0.336798 0.001095 0.434185
27 GB 9.655059e-01 0.000590 0.235065 0.000781 0.309942 9.314079e-01 0.000853 0.339980 0.001106 0.438374
30 GB 9.648083e-01 0.000596 0.237737 0.000789 0.313182 9.302537e-01 0.000860 0.342765 0.001115 0.442052
29 GB 9.648018e-01 0.000596 0.237758 0.000789 0.313210 9.302495e-01 0.000860 0.342768 0.001115 0.442066
16 RF 9.617788e-01 0.000652 0.260076 0.000822 0.326965 7.181760e-01 0.001776 0.708767 0.002242 0.890274
17 RF 9.356989e-01 0.000843 0.336655 0.001067 0.424249 7.056837e-01 0.001825 0.728162 0.002291 0.909529
18 RF 8.918179e-01 0.001094 0.437154 0.001384 0.550749 6.990946e-01 0.001837 0.733207 0.002316 0.920578
22 KN 1.000000e+00 0.000000 0.000000 0.000000 0.000000 6.864486e-01 0.001850 0.737823 0.002364 0.937956
21 KN 1.000000e+00 0.000000 0.000000 0.000000 0.000000 6.864486e-01 0.001850 0.737823 0.002364 0.937956
19 RF 8.996727e-01 0.001046 0.417502 0.001332 0.529705 6.847126e-01 0.001868 0.745232 0.002371 0.940717
24 KN 1.000000e+00 0.000000 0.000000 0.000000 0.000000 6.829574e-01 0.001860 0.741633 0.002378 0.942836
20 RF 8.821701e-01 0.001145 0.456940 0.001444 0.574139 6.792046e-01 0.001890 0.753917 0.002392 0.948955
23 KN 1.000000e+00 0.000000 0.000000 0.000000 0.000000 6.776922e-01 0.001877 0.748579 0.002397 0.950374
25 KN 1.000000e+00 0.000000 0.000000 0.000000 0.000000 6.730756e-01 0.001892 0.754388 0.002414 0.957432
11 DT 6.080479e-01 0.002087 0.833224 0.002633 1.047142 3.954722e-01 0.002627 1.048927 0.003283 1.304359
12 DT 6.080479e-01 0.002087 0.833224 0.002633 1.047142 3.954722e-01 0.002627 1.048927 0.003283 1.304359
13 DT 6.080479e-01 0.002087 0.833224 0.002633 1.047142 3.954722e-01 0.002627 1.048927 0.003283 1.304359
14 DT 5.592260e-01 0.002212 0.883038 0.002793 1.110065 3.722579e-01 0.002682 1.070737 0.003345 1.328524
15 DT 5.503775e-01 0.002235 0.892255 0.002821 1.120969 3.658446e-01 0.002696 1.076524 0.003363 1.335077
6 Stacking 4.677475e-02 0.003296 1.315596 0.004107 1.633356 3.044742e-02 0.003329 1.328519 0.004158 1.652109
8 Stacking 4.639922e-02 0.003297 1.315855 0.004108 1.633677 2.998532e-02 0.003330 1.329002 0.004159 1.652483
7 Stacking 4.639922e-02 0.003297 1.315855 0.004108 1.633677 2.998532e-02 0.003330 1.329002 0.004159 1.652483
9 Stacking 4.639922e-02 0.003297 1.315855 0.004108 1.633677 2.998532e-02 0.003330 1.329002 0.004159 1.652483
10 Stacking 4.634154e-02 0.003297 1.315895 0.004108 1.633727 2.912378e-02 0.003331 1.329389 0.004161 1.653213
36 EN 5.469788e-03 0.003366 1.343430 0.004195 1.668423 1.134387e-03 0.003379 1.348633 0.004220 1.677128
37 EN 5.302387e-03 0.003366 1.343568 0.004195 1.668563 1.133508e-03 0.003379 1.348649 0.004220 1.677123
38 EN 3.303298e-03 0.003370 1.345302 0.004199 1.670220 9.474134e-04 0.003380 1.349052 0.004220 1.677208
39 EN 1.418452e-03 0.003374 1.346637 0.004203 1.671774 4.972900e-04 0.003381 1.349406 0.004221 1.677508
5 Lasso -7.216450e-16 0.003376 1.347487 0.004206 1.672950 1.054712e-15 0.003382 1.349631 0.004222 1.677881
1 Lasso -7.216450e-16 0.003376 1.347487 0.004206 1.672950 1.054712e-15 0.003382 1.349631 0.004222 1.677881
2 Lasso -7.216450e-16 0.003376 1.347487 0.004206 1.672950 1.054712e-15 0.003382 1.349631 0.004222 1.677881
3 Lasso -7.216450e-16 0.003376 1.347487 0.004206 1.672950 1.054712e-15 0.003382 1.349631 0.004222 1.677881
4 Lasso -7.216450e-16 0.003376 1.347487 0.004206 1.672950 1.054712e-15 0.003382 1.349631 0.004222 1.677881
40 EN -7.216450e-16 0.003376 1.347487 0.004206 1.672950 1.054712e-15 0.003382 1.349631 0.004222 1.677881
31 RD 1.053331e-02 0.003357 1.339893 0.004184 1.664425 -1.754346e-03 0.003387 1.351598 0.004226 1.679776
32 RD 1.053331e-02 0.003357 1.339893 0.004184 1.664425 -1.754346e-03 0.003387 1.351598 0.004226 1.679776
33 RD 1.053331e-02 0.003357 1.339893 0.004184 1.664425 -1.754346e-03 0.003387 1.351598 0.004226 1.679776
34 RD 1.053331e-02 0.003357 1.339893 0.004184 1.664425 -1.754348e-03 0.003387 1.351598 0.004226 1.679776
35 RD 1.053331e-02 0.003357 1.339893 0.004184 1.664425 -1.754351e-03 0.003387 1.351598 0.004226 1.679776
0 Linear 1.053331e-02 0.003357 1.339893 0.004184 1.664425 -1.754943e-03 0.003387 1.351598 0.004226 1.679776

There is a significant performance increase in the FNNs, random forest models, and k-nearest neighbor models on the testing data set. All FNN models produce training and testing \(R^2\) above 0.95, with the best above 0.973. The remaining models did not significantly improve performance from the increased sample size. All but the FNNs overfit to the training data set.

Below are the optimal hyperparameter configurations based on test \(R^2\).

[24]:

for model in ["Linear", "Lasso", "DT", "RF", "KN", "GB", "RD", "EN", "Stacking", "FNN"]:
    postprocessor.print_model(model_type=model)
    print()

Model Type: Linear
  copy_X: True
  fit_intercept: True
  n_jobs: None
  positive: False

Model Type: Lasso
  alpha: 0.0001138097715342958

Model Type: DT
  max_depth: 22
  max_features: None
  min_samples_leaf: 11
  min_samples_split: 7

Model Type: RF
  criterion: squared_error
  max_features: 6
  min_samples_leaf: 1
  min_samples_split: 2
  n_estimators: 145

Model Type: KN
  leaf_size: 6
  n_neighbors: 8
  p: 1
  weights: distance

Model Type: GB
  estimator__learning_rate: 0.10997304567272648

Model Type: RD
  alpha: 0.009906793976170773

Model Type: EN
  alpha: 0.0001074768498204962

Model Type: Stacking
  estimator__KN__leaf_size: 22
  estimator__KN__n_neighbors: 1
  estimator__KN__p: 2
  estimator__KN__weights: distance
  estimator__SVM__C: 69.42516231541238

Model Type: FNN
  Structural Hyperparameters
    Layer: Dense_hidden_0
      units: 199
      sublayer: None
    Layer: Dense_hidden_1
      units: 400
      sublayer: Dropout_hidden
    Layer: Dense_hidden_1_sublayer_Dropout_hidden_0
      rate: 0.3225718287912892
    Layer: Dense_output_0
  Compile/Fitting Hyperparameters
    Adam_learning_rate: 0.00011376283985074373
    batch_size: 8
Model: "FNN"
_________________________________________________________________
 Layer (type)                Output Shape              Param #
=================================================================
 Dense_hidden_0 (Dense)      (None, 199)               1791

 Dense_hidden_1 (Dense)      (None, 400)               80000

 Dense_hidden_1_sublayer_Dr  (None, 400)               0
 opout_hidden_0 (Dropout)

 Dense_output_0 (Dense)      (None, 4)                 1604

=================================================================
Total params: 83395 (325.76 KB)
Trainable params: 83395 (325.76 KB)
Non-trainable params: 0 (0.00 Byte)
_________________________________________________________________

Below is the network plot for the FNN.

[25]:

postprocessor.nn_network_plot(
    to_file="./supporting/HTGR_microreactor_network_1.png",
    show_shapes=True,
    show_layer_names=True,
    expand_nested=True,
    show_layer_activations=True,
)

[25]:
../_images/benchmarks_HTGR_microreactor_48_0.png

Once again we can visualize the performance of the optimal models using diagonal validation and validation plots.

[26]:

models = np.array([["Linear", "Lasso"], ["DT", "KN"], ["RF", "FNN"], ["GB", "RD"], ["Stacking", "EN"]])
fig, axarr = plt.subplots(models.shape[0], models.shape[1], sharex=True, sharey=True, figsize=(15,20))
for i in range(models.shape[0]):
    for j in range(models.shape[1]):
        plt.sca(axarr[i, j])
        axarr[i, j] = postprocessor.diagonal_validation_plot(model_type=models[i, j])
        axarr[i, j].set_title(models[i, j])

plt.show()
../_images/benchmarks_HTGR_microreactor_50_0.png 
[27]:

models = np.array([["Linear", "Lasso"], ["DT", "KN"], ["RF", "FNN"], ["GB", "RD"], ["Stacking", "EN"]])
fig, axarr = plt.subplots(models.shape[0], models.shape[1], sharex=True, sharey=True, figsize=(15,20))
for i in range(models.shape[0]):
    for j in range(models.shape[1]):
        plt.sca(axarr[i, j])
        axarr[i, j] = postprocessor.validation_plot(model_type=models[i, j])
        axarr[i, j].set_title(models[i, j])

plt.show()
../_images/benchmarks_HTGR_microreactor_51_0.png

Both plots reaffirm the conclusions from the performance metrics. The FNN model performs significantly better than the other models, with most test data predictions off by less than 1%. Below is the learning curve of the optimal FNN model.

[28]:

fig, ax = plt.subplots()
ax = postprocessor.nn_learning_plot()
plt.show()
../_images/benchmarks_HTGR_microreactor_53_0.png

The validation loss curve is consistently below the training curve, indicating the FNN is not overfitting.

pyMAISElogo.png