Critical Heat Flux (CHF) Prediction

Inputs

D (m): Diameter of the test section (\(0.002 - 0.016~m\))
L (m): Heated length (\(0.07 - 15.0~m\))
P (kPa): Pressure (\(100 - 20,000~kPa\))
G (kg m-2s-1): Mass flux (\(17.7 - 7712.0~\frac{kg}{m^2\cdot s}\))
Tin (C): Inlet temperature length (\(9.0 - 353.62 ^\circ C\))
Xe (-): Outlet equilibrium quality (\(-0.445 - 0.986\))

Output

CHF (kW m-2): Critical heat flux (\(130.0 - 13345.0~\frac{kW}{m^2}\))

The CHF dataset was taken from the Nuclear Regulatory Commission CHF database. This dataset was generated from vertical water-cooled uniformly heated tubes, producing 24,579 samples from 59 different sources of measurements. The dataset was collected over experimental measurements spanning 60 years of CHF data collection methods such as visual identification, physical burnout, changes in the test section resistances, and the usage of thermocouples. The parameters collected consist of pressure (\(P\)), test section diameter (\(D\)), heated length (\(L\)), mass flux (\(G\)), inlet temperature (\(T_{in}\)), outlet equilibrium quality (\(X\)), and \(CHF\).

Negative quality \(X\) could represent a subcooled fluid. The database was limited in terms of diameter (2 \(<\) D \(<\) 25 mm), \(L/D\) ratio (\(L/D\) \(>\) 50 for X \(>\) 0, \(L/D\) \(>\) 25 for X \(<\) 0), pressure (100 \(\leq\) P \(\leq\) 21,000 kPa) and mass flux (0 \(\leq\) G \(<\) 8,000 kg/m2/s). The measured data was also not equally distributed over the whole span, and no data beyond \(D\) = 16 mm was found in the database. Our goal is to identify optimal ML models with pyMIASE that can predict CHF accurately given other geometrical parameters and boundary conditions (:math:`D, L, P, G, T_{in}, X`).

Important Note: This dataset is confidential, it is not open-sourced with pyMAISE and can be requested directly from the OECD Nuclear Energy Agency (NEA). In this notebook, we provide 2500 synthetic samples that mimic the real ones (split into 2000 train and 500 test samples). The synthetic samples were generated by adding random noise to the real experiments. The user only needs to replace the synthetic csv files below with the ones from NEA. More information about the benchmark can be found at this link: https://www.oecd-nea.org/jcms/pl_89619/benchmark-on-artificial-intelligence-and-machine-learning-for-scientific-computing-in-nuclear-engineering-phase-1-critical-heat-flux-exercise-specifications?details=true

[1]:

# Import other packages
import time
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from scipy.stats import uniform, randint
from sklearn.preprocessing import MinMaxScaler

# pyMAISE specific imports
from pyMAISE.datasets import load_chf
from pyMAISE.preprocessing import correlation_matrix, scale_data, read_csv
import pyMAISE as mai

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

Model Initialization

We will examine the performance of 7 regression models:

linear: Linear,
lasso: lasso,
support vector machine SVM,
decision tree DT,
random forest RF,
k-nearest neighbors KN,
dense feedforward neural networks FNN.

For hyperparameter tuning, we initialize all classical models as scikit-learn defaults. For the FNN we define input and output layers with possible dropout layers. These layers include hyperparameter tuning of their number of nodes, use of sublayers, and the rate of dropout. The dense hidden layers include tuning of their depth.

[7]:

model_settings = {
    "models": ["Linear", "Lasso", "SVM", "DT", "RF", "KN", "FNN"],
    "FNN": {
        "structural_params": {
            "Dense_hidden": {
                "num_layers": mai.Int(min_value=0, max_value=4),
                "units": mai.Int(min_value=25, max_value=300),
                "activation": "relu",
                "kernel_initializer": "normal",
                "sublayer": mai.Choice(["Dropout_hidden", "None"]),
                "Dropout_hidden": {
                "rate": mai.Float(min_value=0.4, 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([16, 32, 64]),
            "epochs": 30,
            "validation_split": 0.15,
        },
    },
}
tuner = mai.Tuner(xtrain, ytrain, model_settings=model_settings)

Hyperparameter Tuning

The hyperparameter tuning spaces are defined below. We use random search for the classical models as their training is quick, and random search can cover ample parameter space. Bayesian search is utilized for the FNNs as their training is more computationally expensive. Both search methods use 5-fold cross-validation to mitigate bias from the training data. We train 200 classical models from random search and 50 iterations of Bayesian search for neural networks.

[8]:

random_search_spaces = {
    "Lasso": {
        "alpha": uniform(loc=0.0001, scale=0.0099), # 0.0001 - 0.01
    },
    "SVM": {
        "kernel": ["linear", "poly", "rbf", "sigmoid"],
        "degree": randint(low=1, high=5),
        "gamma": ["scale", "auto"],
    },
    "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
    },
}

start = time.time()
random_search_configs = tuner.random_search(
    param_spaces=random_search_spaces,
    n_iter=200,
    n_jobs=6,
    cv=5,
)

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 23.15457481543223 minutes to process.

After the training, we can see training results for each iteration using the pyMAISE.Tuner.convergence_plot function. For the neural networks, we have the following:

[9]:

ax = tuner.convergence_plot(model_types="FNN")
ax.set_ylim([0, 1])
plt.show()

The Bayesian search explores the search space but finds several suitable hyperparameter configurations.

Model Postprocessing

With the top pyMAISE.Settings.num_configs_saved, we can pass these parameter configurations to the pyMAISE.PostProcessor for model comparison and analysis. For the FNNs, we define the "epochs" parameter as 200 for better performance.

[10]:

postprocessor = mai.PostProcessor(
    data=(xtrain, xtest, ytrain, ytest),
    model_configs=[random_search_configs, bayesian_search_configs],
    new_model_settings={
        "FNN": {"fitting_params": {"epochs": 200}},
    },
    yscaler=yscaler,
)

Below we computed the performance metrics based on the performance of each model on the CHF dataset.

[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
26	FNN	0.990508	94.516480	6.649731	141.619037	10.694737	0.984016	112.056938	7.155676	168.414268	9.765656
29	FNN	0.981713	130.973074	11.190134	196.569625	19.033209	0.980012	133.205079	10.444191	188.331572	16.682770
28	FNN	0.981788	123.758863	9.232552	196.167621	15.062968	0.978672	129.940655	8.839379	194.544377	12.666816
30	FNN	0.980160	137.554895	10.049621	204.747434	17.117261	0.976495	146.234950	9.795921	204.231296	13.868115
27	FNN	0.983183	112.064923	7.608649	188.505708	12.652984	0.976474	125.963180	7.825601	204.320630	11.469703
16	RF	0.993368	65.585293	4.510185	118.381224	7.187412	0.972953	153.625124	10.616529	219.079658	14.802359
20	RF	0.987236	81.077412	5.414823	164.227892	8.852843	0.969859	159.234441	10.801584	231.269415	14.995095
17	RF	0.986148	91.061165	6.335240	171.081482	10.739503	0.968772	161.034753	11.308920	235.404268	15.945121
19	RF	0.984688	94.123667	6.631856	179.871080	11.169824	0.961154	174.426685	12.140031	262.551185	17.195610
18	RF	0.984847	96.270814	6.930750	178.937088	12.111184	0.959509	177.228360	12.392103	268.052022	17.778772
21	KN	1.000000	0.000000	0.000000	0.000000	0.000000	0.933615	215.336825	13.840165	343.222384	20.584267
22	KN	1.000000	0.000000	0.000000	0.000000	0.000000	0.933615	215.336825	13.840165	343.222384	20.584267
24	KN	0.967781	154.904106	9.857600	260.917217	14.309083	0.931029	223.432334	14.515832	349.843983	22.164982
25	KN	0.967781	154.904106	9.857600	260.917217	14.309083	0.931029	223.432334	14.515832	349.843983	22.164982
14	DT	0.989358	82.580507	5.293493	149.951876	7.964006	0.923636	254.440601	17.471253	368.117166	25.283433
13	DT	0.973614	130.161606	8.147538	236.122261	11.623378	0.922344	261.998800	17.658795	371.217626	24.205510
11	DT	0.974340	117.200649	7.370938	232.849637	10.797160	0.922214	261.855286	17.658778	371.529430	24.188682
23	KN	1.000000	0.000000	0.000000	0.000000	0.000000	0.921545	227.629803	14.753816	373.122665	21.699285
12	DT	0.960467	156.372291	9.907215	289.020105	14.289249	0.918241	267.786020	18.265848	380.899293	24.892493
15	DT	0.966189	145.913486	9.179704	267.284736	13.291821	0.917124	269.810126	18.446146	383.491269	25.067672
10	SVM	0.806954	542.860286	60.626177	638.673390	103.960390	0.821172	473.896169	45.830185	563.324304	72.300791
9	SVM	0.806954	542.860286	60.626177	638.673390	103.960390	0.821172	473.896169	45.830185	563.324304	72.300791
8	SVM	0.806954	542.860286	60.626177	638.673390	103.960390	0.821172	473.896169	45.830185	563.324304	72.300791
7	SVM	0.806954	542.860286	60.626177	638.673390	103.960390	0.821172	473.896169	45.830185	563.324304	72.300791
6	SVM	0.806954	542.860286	60.626177	638.673390	103.960390	0.821172	473.896169	45.830185	563.324304	72.300791
0	Linear	0.700590	567.152874	57.689799	795.392588	122.844808	0.702795	512.110104	44.459503	726.221750	74.579545
1	Lasso	0.700212	564.042244	56.414526	795.893993	118.768912	0.702055	508.729649	43.538637	727.125794	72.654205
2	Lasso	0.700203	563.995461	56.395515	795.905770	118.712196	0.702041	508.686986	43.525811	727.142169	72.628259
3	Lasso	0.699871	562.646548	55.829363	796.347214	117.001227	0.701573	507.432404	43.145331	727.713187	71.855605
4	Lasso	0.699842	562.554779	55.789611	796.384942	116.880006	0.701535	507.347537	43.119137	727.759704	71.801809
5	Lasso	0.699836	562.536258	55.781555	796.392753	116.855255	0.701527	507.330195	43.113784	727.769300	71.790839

This data set is not linear, with a poor performance from linear and lasso regression. The remaining four models in order are FNN > RF > KN > DT > SVM with all five models providing reliable performance. KN overfits with a perfect training \(R^2\).

Using the pyMAISE.PostProcessor.print_model function, we can see all the top performing (based on test \(R^2\)) models’ hyperparameters.

[12]:

for model in ["Lasso", "SVM", "DT", "RF", "KN", "FNN"]:
    postprocessor.print_model(model_type=model)
    print()

Model Type: Lasso
  alpha: 0.00017553613590585732

Model Type: SVM
  degree: 4
  gamma: scale
  kernel: rbf

Model Type: DT
  max_depth: 26
  max_features: 6
  min_samples_leaf: 1
  min_samples_split: 7

Model Type: RF
  criterion: squared_error
  max_features: None
  min_samples_leaf: 1
  min_samples_split: 2
  n_estimators: 174

Model Type: KN
  leaf_size: 24
  n_neighbors: 5
  p: 1
  weights: distance

Model Type: FNN
  Structural Hyperparameters
    Layer: Dense_hidden_0
      units: 231
      sublayer: None
    Layer: Dense_hidden_1
      units: 138
      sublayer: None
    Layer: Dense_hidden_2
      units: 267
      sublayer: Dropout_hidden
    Layer: Dense_hidden_2_sublayer_Dropout_hidden_0
      rate: 0.4995897609454529
    Layer: Dense_output_0
  Compile/Fitting Hyperparameters
    Adam_learning_rate: 0.0009311391232267503
    batch_size: 64
Model: "FNN"
_________________________________________________________________
 Layer (type)                Output Shape              Param #
=================================================================
 Dense_hidden_0 (Dense)      (None, 231)               1617

 Dense_hidden_1 (Dense)      (None, 138)               32016

 Dense_hidden_2 (Dense)      (None, 267)               37113

 Dense_hidden_2_sublayer_Dr  (None, 267)               0
 opout_hidden_0 (Dropout)

 Dense_output_0 (Dense)      (None, 1)                 268

=================================================================
Total params: 71014 (277.40 KB)
Trainable params: 71014 (277.40 KB)
Non-trainable params: 0 (0.00 Byte)
_________________________________________________________________

Below is the network plot for the FNN.

[13]:

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

[13]:

We can use the pyMAISE.PostProcessor.diagonal_validation_plot and pyMAISE.PostProcessor.validation_plot functions to visualize the performance of these models. Diagonal validation plots show the actual versus predicted outcome. A well-performing model follows \(y=x\) on this plot. You can choose what outputs to display on the diagonal validation plot; however, since we have a single output (CHF), a single plot will be generated. Both functions use the top-performing model on test \(R^2\) by default.

[14]:

models = np.array([["Linear", "Lasso"], ["DT", "KN"], ["RF", "SVM"]])
fig = plt.figure(figsize=(15, 25))
gs = fig.add_gridspec(4, 2)
for i in range(models.shape[0]):
    for j in range(models.shape[1]):
        ax = fig.add_subplot(gs[i, j])
        ax = postprocessor.diagonal_validation_plot(model_type=models[i, j])
        ax.set_title(models[i, j])
ax = fig.add_subplot(gs[3, :])
ax = postprocessor.diagonal_validation_plot(model_type="FNN")
ax.set_title("FNN")
plt.show()

All models display a close spread to \(y=x\); however, given the more extensive spread, linear models’ lack of performance is apparent.

Validation plots show the absolute error of the predicted output relative to the actual outputs of the testing data set.

[15]:

fig = plt.figure(figsize=(15, 25))
gs = fig.add_gridspec(4, 2)
for i in range(models.shape[0]):
    for j in range(models.shape[1]):
        ax = fig.add_subplot(gs[i, j])
        ax = postprocessor.validation_plot(model_type=models[i, j])
        ax.set_title(models[i, j])
ax = fig.add_subplot(gs[3, :])
ax = postprocessor.validation_plot(model_type="FNN")
ax.set_title("FNN")
plt.show()

Finally, we can see if the FNN is overfit based on its learning curve.

[16]:

postprocessor.nn_learning_plot()
plt.show()

The FNN is not overfitting, as the validation curve closely follows the training curve.

Finally, using the best FNN model, we can generate new CHF predictions for new experimental setups that do not exist in our data set.

Critical Heat Flux (CHF) Prediction

Initialization and Preprocessing

Model Initialization

Hyperparameter Tuning

Model Postprocessing